[WIP] fix for writing of virtuals in Desmond

intermol/desmond/desmond_parser.py

@@ -1,3 +1,4 @@
+from collections import OrderedDict, defaultdict
 import logging
 from warnings import warn
 import math
@@ -111,6 +112,24 @@ class DesmondParser(object):
     lookup_desmond_bonds = create_lookup(desmond_bonds)  # not unique
     desmond_bond_types = create_type(desmond_bonds)
 
+    # desmond virtuals are not quite the same as gromacs virtuals.  
+    # lc2, lc3, out3 are the same. 
+    # fdat3 is generalization of gromacs 3fd and 3fda. Eventually this conversion
+    # should be moved to a canonicalization file.
+    # 4fdn is not defined in desmond
+
+    desmond_virtuals = {
+        'lc2': TwoVirtual,
+        'lc3': ThreeLinearVirtual,
+        'fdat3': ThreeFdVirtual,
+        'fdat3': ThreeFadVirtual,
+        'out3': ThreeOutVirtual,
+        None: FourFdnVirtual
+    }
+
+    lookup_desmond_virtuals = create_lookup(desmond_virtuals)  # not unique
+    desmond_virtual_types = create_type(desmond_virtuals)
+
     def canonical_bond(self, bond, params, direction='into', name=None):
 
         if direction == 'into':
@@ -1162,15 +1181,21 @@ def write_vdwtypes_and_sites(self, molecule):
         combrule = self.system.combination_rule
         for atom in molecule.atoms:
             i+=1
+            if atom.ptype == 'A':
+                atom_type = 'atom'
+            elif atom.ptype == 'D':
+                atom_type = 'pseudo'
+            else:
+                raise DesmondError('The physical type of an atom is {0}, which is not allowed'.format(atom.ptype))
             if atom.residue_index:
-                sites.append(' %3d %5s %9.8f %9.8f %2s %1d %4s\n' % (
-                        i, 'atom',
+                sites.append(' %3d %6s %9.8f %9.8f %2s %1d %4s\n' % (
+                        i, atom_type,
                         atom._charge[0].value_in_unit(units.elementary_charge),
                         atom._mass[0].value_in_unit(units.atomic_mass_unit),
                         atom.atomtype[0], atom.residue_index, atom.residue_name))
             else:
                 sites.append(' %3d %5s %9.8f %9.8f %2s\n' % (
-                        i, 'atom',
+                        i, atom_type,
                         atom._charge[0].value_in_unit(units.elementary_charge),
                         atom._mass[0].value_in_unit(units.atomic_mass_unit),
                         atom.atomtype[0]))
@@ -1649,6 +1674,114 @@ def write_constraints(self, moleculetype):
         hlines.extend(dlines)
         return hlines
 
+    def write_virtuals(self, moleculetype):
+
+        #ADDING VIRTUALS
+        logger.debug("   -Writing virtuals...")
+
+        virtuals = defaultdict(list)
+        virtuallist = sorted(list(moleculetype.virtual_forces),key=lambda x: x.atom1)
+
+        dlines = list()
+        hlines = list()
+
+        for virtual in virtuallist:
+            if hasattr(virtual, 'atom5'):
+                virtuals[4].append(virtual)
+            elif hasattr(virtual, 'atom4'):
+                virtuals[3].append(virtual)
+            else:
+                virtuals[2].append(virtual)
+
+        virtuals[2] = sorted(virtuals[2], key=lambda x: (x.atom1, x.atom2, x.atom3))
+        virtuals[3] = sorted(virtuals[3], key=lambda x: (x.atom1, x.atom2, x.atom3, x.atom4))
+        virtuals[4] = sorted(virtuals[4], key=lambda x: (x.atom1, x.atom2, x.atom3, x.atom4, x.atom5))
+
+        alen_max = 2 # figure out the maximum number of atoms in all of the virtuals.
+        for j in range(2,5):
+            if virtuals[j] > 0:
+                alen_max = j
+
+        i = 0
+
+        clen_max = 3 # maximum number of parameters involved. We won't try to automatically find. 
+        for n_body_type, vsites in virtuals.items():
+            for vsite in vsites:
+                i += 1
+                cline = '      %d ' % i
+                for n in range(1, n_body_type + 2):
+                    atom = getattr(vsite, 'atom{}'.format(n))
+                    cline += '{0:7d} '.format(atom)
+                for j in range(n_body_type + 1, alen_max):
+                    cline += ' <>'
+                cline += ' {:4s} '.format(self.lookup_desmond_virtuals[vsite.__class__])
+
+                vsite_params = self.get_parameter_list_from_force(vsite)
+                param_units = self.unitvars[vsite.__class__.__name__]
+                for param, unit in zip(vsite_params, param_units):
+                    cline += '{0:18.8e}'.format(param.value_in_unit(unit))
+                for j in range(len(vsite_params),clen_max):
+                    cline += ' <>'
+                cline += '\n'
+                dlines.append(cline)
+
+        hlines.append("    ffio_virtuals[%d] {\n" % (i))
+        if i == 0:
+            hlines.append("      :::\n")
+        else:
+        #decide on alen and clen based on the type of virtual
+            letters = ['i','j','k','l','m']
+            hlines.append('      i_ffio_index\n')
+            for j in range(alen_max-1):
+                hlines.append('      i_ffio_a%s\n' % letters[j])
+            hlines.append('      s_ffio_funct\n')
+            for j in range(clen_max):
+                hlines.append('      r_ffio_c%d\n' %(j+1))
+            hlines.append("      :::\n")
+
+        dlines.append("      :::\n")
+        dlines.append("    }\n")
+        hlines.extend(dlines)
+        return hlines
+
+    def write_pseudos(self, moleculetype):
+
+        dlines = list()
+        hlines = list()
+
+        i = 0
+        for molecule in moleculetype.molecules:
+            for atom in molecule.atoms:
+                if atom.ptype == 'D':
+                    i += 1
+                    line = '      %d ' % i
+                    for j in range(3):
+                        line += " %10.8f" % (float(atom._position[j].value_in_unit(units.angstroms)))
+                    line +=   "     %2d %4s  %2s" % (
+                        atom.residue_index,
+                        '"%s"'%atom.residue_name,
+                        '"%s"'%atom.name)
+                    dlines.append(line + '\n')
+
+        hlines.append("    ffio_pseudo[%d] {\n" % (i))
+        if i == 0:
+            hlines.append("      :::\n")
+        else:
+            hlines.append('      r_ffio_x_coord\n')
+            hlines.append('      r_ffio_y_coord\n')
+            hlines.append('      r_ffio_z_coord\n')
+            hlines.append('      i_ffio_residue_number\n')
+            hlines.append('      s_ffio_pdb_residue_name\n')
+            hlines.append('      s_ffio_atom_name\n')
+            hlines.append("      :::\n")
+
+        dlines.append("      :::\n")
+        dlines.append("    }\n")
+
+        hlines.extend(dlines)
+
+        return hlines
+
     def write(self):
 
 #        Write this topology to file
@@ -1703,6 +1836,8 @@ def write(self):
         for moleculetype in self.system._molecule_types.values():
             for molecule in moleculetype.molecules:
                 for atom in molecule.atoms:
+                    if atom.ptype == 'D': # dummy (virtual) atoms do not go here.
+                        continue
                     i += 1
                     line = '    %d        %d' % (i,1) #HAVE TO PUT THE 1 HERE OR ELSE DESMOND DIES, EVEN THOUGH IT DOESN'T USE IT
                     for j in range(3):
@@ -1783,7 +1918,8 @@ def write(self):
 
         for molecule_name, moleculetype in self.system.molecule_types.items():
             logger.debug('Writing molecule block %s...' % (molecule_name))
-            #BEGINNING BLOCK
+
+            #BEGINNING BLOCK - should be separate function.
 
             logger.debug("  Writing f_m_ct...")
             lines.append('f_m_ct {\n')
@@ -1798,7 +1934,8 @@ def write(self):
                 for bj in range(3):
                     lines.append('%22s\n' % float(bv[bi][bj].value_in_unit(units.angstroms)))
             lines.append('  solute\n')
-            #M_ATOMS
+
+            #M_ATOMS - should be separate function
 
             logger.debug("  Writing m_atoms...")
             apos = len(lines) #pos of where m_atom will be; will need to overwite later based on the number of atoms
@@ -1812,6 +1949,8 @@ def write(self):
             i = 0
             for molecule in moleculetype.molecules:
                 for atom in molecule.atoms:
+                    if atom.ptype == 'D': # dummy (virtual) atoms do not go here.
+                        continue
                     i += 1
                     #NOT SURE WHAT TO PUT FOR MMOD TYPE; 1 is currently used.
                     #This can't be determined currently from the information provided,
@@ -1836,7 +1975,8 @@ def write(self):
             lines.append('    :::\n')
             lines.append('  }\n')
 
-            #M_BONDS
+
+            #M_BONDS - should be separate function
             logger.debug("  Writing m_bonds...")
 
             hlines = list()
@@ -1918,6 +2058,8 @@ def write(self):
             lines += self.write_torsion_torsion(moleculetype)
             lines += self.write_exclusions(moleculetype)
             lines += self.write_pairs(moleculetype)
+            lines += self.write_virtuals(moleculetype)
+            lines += self.write_pseudos(moleculetype)
             lines += self.write_constraints(moleculetype)
 
             #STILL NEED TO ADD RESTRAINTS

intermol/gromacs/gromacs_parser.py

@@ -398,6 +398,7 @@ def read(self):
         self.pairtypes = dict()
         self.cmaptypes = dict()
         self.nonbondedtypes = dict()
+        self.virtualtypes = dict()
 
         # Parse the top_filename into a set of plain text, intermediate
         # TopMoleculeType objects.
@@ -797,6 +798,7 @@ def create_atom(self, temp_atom):
                                 "{0}. Visually inspect before using.".format(atom))
             atom.sigma = (state, intermol_atomtype.sigma)
             atom.epsilon = (state, intermol_atomtype.epsilon)
+            atom.ptype = intermol_atomtype.ptype
 
         self.current_molecule.add_atom(atom)
         self.n_atoms_added += 1
@@ -1557,6 +1559,7 @@ def process_nonbond_params(self, line):
         self.system.nonbonded_types[tuple(nonbonded_vars)] = nonbonded_type
 
     def process_virtual_sites(self, line, v_site_type):
+
         """Process a line in a [ virtual_sites? ] category."""
         if v_site_type == 'n':
             raise UnimplementedSetting('Parsing of [ virtual_sitesn ] directives'