Merge pull request #2008 from sunhwan/master

handles virtual_sites2 and nbfix in gromacs topology

wrappers/python/simtk/openmm/app/charmmpsffile.py

@@ -342,7 +342,7 @@ class CharmmPsfFile(object):
             # We have a CHARMM PSF file; now do NUMLP/NUMLPH sections
             numlp, numlph = psfsections['NUMLP NUMLPH'][0]
             if numlp != 0 or numlph != 0:
-                raise NotImplemented('Cannot currently handle PSFs with lone '
+                raise NotImplementedError('Cannot currently handle PSFs with lone '
                                           'pairs defined in the NUMLP/NUMLPH '
                                           'section.')
         # Now do the CMAPs

wrappers/python/simtk/openmm/app/gromacstopfile.py

@@ -43,7 +43,8 @@ import math
 import os
 import re
 import distutils.spawn
-from collections import OrderedDict
+from collections import OrderedDict, defaultdict
+from itertools import combinations, combinations_with_replacement
 
 HBonds = ff.HBonds
 AllBonds = ff.AllBonds
@@ -115,7 +116,9 @@ class GromacsTopFile(object):
             self.exclusions = []
             self.pairs = []
             self.cmaps = []
+            self.vsites2 = []
             self.has_virtual_sites = False
+            self.has_nbfix_terms = False
 
     def _processFile(self, file):
         append = ''
@@ -249,6 +252,10 @@ class GromacsTopFile(object):
                 self._processPairType(line)
             elif self._currentCategory == 'cmaptypes':
                 self._processCmapType(line)
+            elif self._currentCategory == 'nonbond_params':
+                self._processNonbondType(line)
+            elif self._currentCategory == 'virtual_sites2':
+                self._processVirtualSites2(line)
             elif self._currentCategory.startswith('virtual_sites'):
                 if self._currentMoleculeType is None:
                     raise ValueError('Found %s before [ moleculetype ]' %
@@ -435,6 +442,22 @@ class GromacsTopFile(object):
             raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line)
         self._cmapTypes[tuple(fields[:5])] = fields
 
+    def _processNonbondType(self, line):
+        """Process a line in the [ nonbond_params ] category."""
+        fields = line.split()
+        if len(fields) < 5:
+            raise ValueError('Too few fields in [ nonbond_params ] line: '+line)
+        if fields[2] != '1':
+            raise ValueError('Unsupported function type in [ nonbond_params ] line: '+line)
+        self._nonbondTypes[tuple(sorted(fields[:2]))] = fields
+
+    def _processVirtualSites2(self, line):
+        """Process a line in the [ virtual_sites2 ] category."""
+        fields = line.split()
+        if len(fields) < 5:
+            raise ValueError('Too few fields in [ virtual_sites2 ] line: ' + line)
+        self._currentMoleculeType.vsites2.append(fields[:5])
+
     def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):
         """Load a top file.
 
@@ -484,6 +507,7 @@ class GromacsTopFile(object):
         self._implicitTypes = {}
         self._pairTypes = {}
         self._cmapTypes = {}
+        self._nonbondTypes = {}
         self._processFile(file)
 
         # Create the Topology from it.
@@ -643,6 +667,24 @@ class GromacsTopFile(object):
                 for types in dihedralTypeTable.values():
                     types.append(key)
 
+        # NBFIX
+
+        atom_types = []
+        for moleculeName, moleculeCount in self._molecules:
+            moleculeType = self._moleculeTypes[moleculeName]
+            for _ in range(moleculeCount):
+                for atom in moleculeType.atoms:
+                    atom_types.append(atom[1])
+        has_nbfix_terms = any([pair in self._nonbondTypes for pair in combinations_with_replacement(sorted(set(atom_types)), 2)])
+
+        if has_nbfix_terms:
+            # Build a lookup table and angle/dihedral indices list to
+            # let us handle exclusion manually.
+            angleIndices = []
+            torsionIndices = []
+            atom_partners = defaultdict(lambda : defaultdict(set))
+            atom_charges = []
+
         # Loop over molecules and create the specified number of each type.
 
         for moleculeName, moleculeCount in self._molecules:
@@ -846,12 +888,19 @@ class GromacsTopFile(object):
                     else:
                         q = float(params[4])
 
+                    if has_nbfix_terms:
+                        if self._defaults[1] != '2':
+                            raise NotImplementedError('NBFIX terms with LB combination rule is not yet supported')
+                        nb.addParticle(q, 1.0, 0.0)
+                        atom_charges.append(q)
+                    else:
                         if self._defaults[1] == '1':
                             nb.addParticle(q, 1.0, 0.0)
                             lj.addParticle([float(params[6]), float(params[7])])
                         elif self._defaults[1] == '2':
                             nb.addParticle(q, float(params[6]), float(params[7]))
 
+                for fields in moleculeType.atoms:
                     if implicitSolvent is OBC2:
                         if fields[1] not in self._implicitTypes:
                             raise ValueError('No implicit solvent parameters specified for atom type: '+fields[1])
@@ -860,6 +909,23 @@ class GromacsTopFile(object):
                 for fields in moleculeType.bonds:
                     atoms = [int(x)-1 for x in fields[:2]]
                     bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))
+                if has_nbfix_terms:
+                    for fields in moleculeType.bonds:
+                        atoms = [int(x)-1 for x in fields[:2]]
+                        atom_partners[baseAtomIndex+atoms[0]]['bond'].add(baseAtomIndex+atoms[1])
+                        atom_partners[baseAtomIndex+atoms[1]]['bond'].add(baseAtomIndex+atoms[0])
+                    for fields in moleculeType.angles:
+                        atoms = [int(x)-1 for x in fields[:3]]
+                        angleIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2]))
+                        for pair in combinations(atoms, 2):
+                            atom_partners[baseAtomIndex+pair[0]]['angle'].add(baseAtomIndex+pair[1])
+                            atom_partners[baseAtomIndex+pair[1]]['angle'].add(baseAtomIndex+pair[0])
+                    for fields in moleculeType.dihedrals:
+                        atoms = [int(x)-1 for x in fields[:4]]
+                        torsionIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3]))
+                        for pair in combinations(atoms, 2):
+                            atom_partners[baseAtomIndex+pair[0]]['torsion'].add(baseAtomIndex+pair[1])
+                            atom_partners[baseAtomIndex+pair[1]]['torsion'].add(baseAtomIndex+pair[0])
 
                 # Record nonbonded exceptions.
 
@@ -886,12 +952,66 @@ class GromacsTopFile(object):
                         if atom > atoms[0]:
                             exclusions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom))
 
+                # Record virtual sites
+
+                for fields in moleculeType.vsites2:
+                    atoms = [int(x)-1 for x in fields[:3]]
+                    c1 = float(fields[4])
+                    vsite = mm.TwoParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], (1-c1), c1)
+                    sys.setVirtualSite(baseAtomIndex+atoms[0], vsite)
 
         # Create nonbonded exceptions.
 
+        if not has_nbfix_terms:
             nb.createExceptionsFromBonds(bondIndices, fudgeQQ, fudgeLJ)
+        else:
+            excluded_atom_pairs = set() # save these pairs so we don't zero them out
+            for tor in torsionIndices:
+                # First check to see if atoms 1 and 4 are already excluded because
+                # they are 1-2 or 1-3 pairs (would happen in 6-member rings or
+                # fewer). Then check that they're not already added as exclusions
+                if 'bond' in atom_partners[tor[3]] and tor[0] in atom_partners[tor[3]]['bond']: continue
+                if 'angle' in atom_partners[tor[3]] and tor[0] in atom_partners[tor[3]]['angle']: continue
+                key = min((tor[0], tor[3]),
+                          (tor[3], tor[0]))
+                if key in excluded_atom_pairs: continue # multiterm...
+                params1 = self._atomTypes[atom_types[tor[0]]]
+                params4 = self._atomTypes[atom_types[tor[3]]]
+                q1 = atom_charges[tor[0]]
+                rmin1 = float(params1[6])
+                eps1 = float(params1[7])
+                q4 = atom_charges[tor[3]]
+                rmin4 = float(params4[6])
+                eps4 = float(params4[7])
+
+                charge_prod = (q1*q4)
+                epsilon = (math.sqrt(abs(eps1 * eps4)))
+                rmin14 = (rmin1 + rmin4) / 2
+                nb.addException(tor[0], tor[3], charge_prod, rmin14, epsilon)
+                excluded_atom_pairs.add(key)
+
+            # Add excluded atoms
+            for atom_idx, atom in atom_partners.items():
+                # Exclude all bonds and angles
+                for atom2 in atom['bond']:
+                    if atom2 > atom_idx:
+                        nb.addException(atom_idx, atom2, 0.0, 1.0, 0.0)
+                        excluded_atom_pairs.add((atom_idx, atom2))
+                for atom2 in atom['angle']:
+                    if ((atom_idx, atom2) in excluded_atom_pairs):
+                        continue
+                    if atom2 > atom_idx:
+                        nb.addException(atom_idx, atom2, 0.0, 1.0, 0.0)
+                        excluded_atom_pairs.add((atom_idx, atom2))
+                for atom2 in atom['dihedral']:
+                    if atom2 <= atom_idx: continue
+                    if ((atom_idx, atom2) in excluded_atom_pairs):
+                        continue
+                    nb.addException(atom_idx, atom2, 0.0, 1.0, 0.0)
+
         for exclusion in exclusions:
             nb.addException(exclusion[0], exclusion[1], 0.0, 1.0, 0.0, True)
+
         if self._defaults[1] == '1':
             # We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce
             # to handle the exceptions.
@@ -932,6 +1052,88 @@ class GromacsTopFile(object):
             lj.setNonbondedMethod(methodMap[nonbondedMethod])
             lj.setCutoffDistance(nonbondedCutoff)
 
+        if has_nbfix_terms:
+            if self._defaults[1] != '2':
+                raise NotImplementedError('NBFIX terms with LB combination rule is not yet supported')
+
+            atom_nbfix_types = set([])
+            for pair in self._nonbondTypes:
+                atom_nbfix_types.add(pair[0])
+                atom_nbfix_types.add(pair[1])
+
+            lj_idx_list = [0 for _ in atom_types]
+            lj_radii, lj_depths = [], []
+            num_lj_types = 0
+            lj_type_list = []
+            for i,atom_type in enumerate(atom_types):
+                atom = self._atomTypes[atom_type]
+                if lj_idx_list[i]: continue # already assigned
+                atom_rmin = atom[6]
+                atom_epsilon = atom[7]
+                num_lj_types += 1
+                lj_idx_list[i] = num_lj_types
+                ljtype = (atom_rmin, atom_epsilon)
+                lj_type_list.append(atom)
+                lj_radii.append(float(atom_rmin))
+                lj_depths.append(float(atom_epsilon))
+                for j in range(i+1, len(atom_types)):
+                    atom_type2 = atom_types[j]
+                    if lj_idx_list[j] > 0: continue # already assigned
+                    atom2 = self._atomTypes[atom_type2]
+                    atom2_rmin = atom2[6]
+                    atom2_epsilon = atom2[7]
+                    if atom2 is atom:
+                        lj_idx_list[j] = num_lj_types
+                    elif atom_type not in atom_nbfix_types:
+                        # Only non-NBFIXed atom types can be compressed
+                        ljtype2 = (atom2_rmin, atom2_epsilon)
+                        if ljtype == ljtype2:
+                            lj_idx_list[j] = num_lj_types
+
+            # Now everything is assigned. Create the A-coefficient and
+            # B-coefficient arrays
+            acoef = [0 for i in range(num_lj_types*num_lj_types)]
+            bcoef = acoef[:]
+            for i in range(num_lj_types):
+                namei = lj_type_list[i][0]
+                for j in range(num_lj_types):
+                    namej = lj_type_list[j][0]
+                    try:
+                        params = self._nonbondTypes[tuple(sorted((namei, namej)))]
+                        rij = float(params[3])
+                        wdij = float(params[4])
+                    except KeyError:
+                        rij = (lj_radii[i] + lj_radii[j]) / 2
+                        wdij = math.sqrt(lj_depths[i] * lj_depths[j])
+                    acoef[i+num_lj_types*j] = 2 * math.sqrt(wdij) * rij**6
+                    bcoef[i+num_lj_types*j] = 4 * wdij * rij**6
+            cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
+                                             'a=acoef(type1, type2);'
+                                             'b=bcoef(type1, type2)')
+            cforce.addTabulatedFunction('acoef',
+                    mm.Discrete2DFunction(num_lj_types, num_lj_types, acoef))
+            cforce.addTabulatedFunction('bcoef',
+                    mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef))
+            cforce.addPerParticleParameter('type')
+            if (nonbondedMethod in (ff.PME, ff.LJPME, ff.Ewald, ff.CutoffPeriodic)):
+                cforce.setNonbondedMethod(cforce.CutoffPeriodic)
+                cforce.setCutoffDistance(nonbondedCutoff)
+                cforce.setUseLongRangeCorrection(True)
+            elif nonbondedMethod is ff.NoCutoff:
+                cforce.setNonbondedMethod(cforce.NoCutoff)
+            elif nonbondedMethod is ff.CutoffNonPeriodic:
+                cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
+                cforce.setCutoffDistance(nonbondedCutoff)
+            else:
+                raise ValueError('Unrecognized nonbonded method')
+            for i in lj_idx_list:
+                cforce.addParticle((i - 1,)) # adjust for indexing from 0
+
+            for i in range(nb.getNumExceptions()):
+                ii, jj, q, eps, sig = nb.getExceptionParameters(i)
+                cforce.addExclusion(ii, jj)
+            sys.addForce(cforce)
+
         # Adjust masses.
 
         if hydrogenMass is not None:

wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py

@@ -531,7 +531,7 @@ class ResidueList(list):
         return atom
 
     def append(self, thing):
-        raise NotImplemented('Use "add_atom" to build a residue list')
+        raise NotImplementedError('Use "add_atom" to build a residue list')
 
     extend = append
 

wrappers/python/tests/TestGromacsTopFile.py

@@ -151,6 +151,24 @@ class TestGromacsTopFile(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
 
+    def test_VirtualParticle(self):
+        """Test virtual particle works correctly."""
+
+        top = GromacsTopFile('systems/bnz.top')
+        gro = GromacsGroFile('systems/bnz.gro')
+        system = top.createSystem()
+
+        self.assertEqual(26, system.getNumParticles())
+        self.assertEqual(1, len(top._moleculeTypes['BENX'].vsites2))
+
+        context = Context(system, VerletIntegrator(1*femtosecond),
+                          Platform.getPlatformByName('Reference'))
+        context.setPositions(gro.positions)
+        context.computeVirtualSites()
+        ene = context.getState(getEnergy=True).getPotentialEnergy()
+        # the energy output is from gromacs and it only prints out 6 sig digits.
+        self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), 1.88855e+02, places=3)
+
 if __name__ == '__main__':
     unittest.main()
 

wrappers/python/tests/systems/bnz.gro

+Gnomes, ROck Monsters And Chili Sauce
+   26
+    1BEN     CG    1   0.052  -0.016   0.455
+    1BEN     HG    2   0.064  -0.020   0.563
+    1BEN    CD1    3   0.004  -0.129   0.386
+    1BEN    HD1    4  -0.020  -0.219   0.441
+    1BEN    CD2    5   0.084   0.100   0.385
+    1BEN    HD2    6   0.120   0.187   0.438
+    1BEN    CE1    7  -0.012  -0.125   0.247
+    1BEN    HE1    8  -0.048  -0.211   0.194
+    1BEN    CE2    9   0.069   0.104   0.245
+    1BEN    HE2   10   0.094   0.194   0.191
+    1BEN     CZ   11   0.022  -0.008   0.177
+    1BEN     HZ   12   0.011  -0.005   0.069
+    1BEN    LPA   13   0.037  -0.012   0.316
+    1BEN     CG   14   0.036   0.038  -0.433
+    1BEN     HG   15   0.092   0.057  -0.523
+    1BEN    CD1   16  -0.054  -0.070  -0.428
+    1BEN    HD1   17  -0.068  -0.133  -0.515
+    1BEN    CD2   18   0.052   0.120  -0.321
+    1BEN    HD2   19   0.120   0.205  -0.326
+    1BEN    CE1   20  -0.127  -0.095  -0.311
+    1BEN    HE1   21  -0.196  -0.178  -0.307
+    1BEN    CE2   22  -0.019   0.094  -0.203
+    1BEN    HE2   23  -0.004   0.157  -0.116
+    1BEN     CZ   24  -0.108  -0.014  -0.198
+    1BEN     HZ   25  -0.163  -0.034  -0.108
+    1BEN    LPA   26  -0.036   0.012  -0.316
+   0.00000   0.00000   0.00000

wrappers/python/tests/systems/bnz.top

+;
+;	File '6lyz_gmx.top' was generated
+;	By user: sunhwan (40816)
+;	On host: kuta
+;	At date: Fri Nov 10 13:02:26 2017
+;
+;	This is a standalone topology file
+;
+;	Created by:
+;	           :-) GROMACS - gmx pdb2gmx, 2017-dev-20170822-762c6f0 (-:
+;	
+;	Executable:   /home/sunhwan/local/gromacs/src/gromacs/build/bin/gmx
+;	Data prefix:  /home/sunhwan/local/gromacs/src/gromacs (source tree)
+;	Working dir:  /data/sunhwan/work/lysozyme/1_setup
+;	Command line:
+;	  gmx pdb2gmx -f 6lyz_4pdb2gmx.pdb -o 6lyz_gmx.pdb -p 6lyz_gmx.top -ff charmm36 -water tip3p -ter -merge all
+;	Force field was read from current directory or a relative path - path added.
+;
+
+[ defaults ]
+; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
+1	2	yes	1.0	1.0
+
+[ atomtypes ]
+;type atnum         mass   charge  ?           sigma  epsilon
+CG2R61    6    12.011000    0.000  A  0.355005321205  0.29288
+HGR61     1     1.008000    0.000  A  0.242003727796  0.12552
+LP        0     0.000000    0.00   A  0.0             0.0    ; pram
+
+; NBFIX for LP, rmin=1.2/2^(1/6), eps=4.184*.01
+[ nonbond_params ]
+LP     LP  1  1.069078461768  0.041840000000                 ; pram
+
+[ bondtypes ]
+;      i        j  func           b0           kb
+  CG2R61   CG2R61     1   0.13750000    255224.00
+  CG2R61    HGR61     1   0.10800000    284512.00
+
+[ angletypes ]
+;      i        j        k  func       theta0       ktheta          ub0          kub
+  CG2R61   CG2R61   CG2R61     5   120.000000   334.720000   0.24162000     29288.00
+  CG2R61   CG2R61    HGR61     5   120.000000   251.040000   0.21525000     18409.60
+
+[ dihedraltypes ]
+;      i        j        k        l  func         phi0         kphi  mult
+  CG2R61   CG2R61   CG2R61   CG2R61     9   180.000000    12.970400     2
+  CG2R61   CG2R61   CG2R61    HGR61     9   180.000000    17.572800     2
+   HGR61   CG2R61   CG2R61    HGR61     9   180.000000    10.041600     2
+
+[ dihedraltypes ]
+; 'improper' dihedrals 
+;      i        j        k        l  func         phi0         kphi
+
+[ moleculetype ]
+; Name            nrexcl
+BENX              3
+
+[ atoms ]
+;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
+; residue   1 BENX rtp BENX q  0.0
+     1     CG2R61      1   BENX     CG      1     -0.115     12.011   ; qtot -0.115
+     2      HGR61      1   BENX     HG      1      0.115      1.008   ; qtot 0
+     3     CG2R61      1   BENX    CD1      1     -0.115     12.011   ; qtot -0.115
+     4      HGR61      1   BENX    HD1      1      0.115      1.008   ; qtot 0
+     5     CG2R61      1   BENX    CD2      1     -0.115     12.011   ; qtot -0.115
+     6      HGR61      1   BENX    HD2      1      0.115      1.008   ; qtot 0
+     7     CG2R61      1   BENX    CE1      1     -0.115     12.011   ; qtot -0.115
+     8      HGR61      1   BENX    HE1      1      0.115      1.008   ; qtot 0
+     9     CG2R61      1   BENX    CE2      1     -0.115     12.011   ; qtot -0.115
+    10      HGR61      1   BENX    HE2      1      0.115      1.008   ; qtot 0
+    11     CG2R61      1   BENX     CZ      1     -0.115     12.011   ; qtot -0.115
+    12      HGR61      1   BENX     HZ      1      0.115      1.008   ; qtot 0
+    13         LP      1   BENX    LPA      1      0.000      0.000   ; qtot 0
+
+[ bonds ]
+;  ai    aj funct            c0            c1            c2            c3
+    1     2     1 
+    1     3     1 
+    1     5     1 
+    3     4     1 
+    3     7     1 
+    5     6     1 
+    5     9     1 
+    7     8     1 
+    7    11     1 
+    9    10     1 
+    9    11     1 
+   11    12     1 
+
+[ pairs ]
+;  ai    aj funct            c0            c1            c2            c3
+    1     8     1 
+    1    10     1 
+    1    11     1 
+    2     4     1 
+    2     6     1 
+    2     7     1 
+    2     9     1 
+    3     6     1 
+    3     9     1 
+    3    12     1 
+    4     5     1 
+    4     8     1 
+    4    11     1 
+    5     7     1 
+    5    12     1 
+    6    10     1 
+    6    11     1 
+    7    10     1 
+    8     9     1 
+    8    12     1 
+   10    12     1 
+
+[ angles ]
+;  ai    aj    ak funct            c0            c1            c2            c3
+    2     1     3     5 
+    2     1     5     5 
+    3     1     5     5 
+    1     3     4     5 
+    1     3     7     5 
+    4     3     7     5 
+    1     5     6     5 
+    1     5     9     5 
+    6     5     9     5 
+    3     7     8     5 
+    3     7    11     5 
+    8     7    11     5 
+    5     9    10     5 
+    5     9    11     5 
+   10     9    11     5 
+    7    11     9     5 
+    7    11    12     5 
+    9    11    12     5 
+
+[ dihedrals ]
+;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
+    2     1     3     4     9 
+    2     1     3     7     9 
+    5     1     3     4     9 
+    5     1     3     7     9 
+    2     1     5     6     9 
+    2     1     5     9     9 
+    3     1     5     6     9 
+    3     1     5     9     9 
+    1     3     7     8     9 
+    1     3     7    11     9 
+    4     3     7     8     9 
+    4     3     7    11     9 
+    1     5     9    10     9 
+    1     5     9    11     9 
+    6     5     9    10     9 
+    6     5     9    11     9 
+    3     7    11     9     9 
+    3     7    11    12     9 
+    8     7    11     9     9 
+    8     7    11    12     9 
+    5     9    11     7     9 
+    5     9    11    12     9 
+   10     9    11     7     9 
+   10     9    11    12     9 
+
+
+[ virtual_sites2 ]
+; Vsite from                    funct   a 
+13       1       11       1      0.500
+
+[ system ]
+; Name
+two benzene
+
+[ molecules ]
+; Compound        #mols
+BENX     2