Merge branch 'nbfix' of https://github.com/ChayaSt/openmm into lj

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

@@ -2157,6 +2157,146 @@ class NonbondedGenerator(object):
 
 parsers["NonbondedForce"] = NonbondedGenerator.parseElement
 
+class LennardJonesGenerator(object):
+    """A NBFix generator to construct the L-J force with NBFIX implemented as a lookup table"""
+
+    def __init__(self, forcefield, lj14scale):
+        self.ff = forcefield
+        self.nbfix_types = {}
+        self.lj14scale = lj14scale
+        self.lj_types = ForceField._AtomTypeParameters(forcefield, 'LennardJonesForce', 'Atom', ('sigma', 'epsilon'))
+
+    def registerNBFIX(self, parameters):
+        types = self.ff._findAtomTypes(parameters, 2)
+        if None not in types:
+            type1 = types[0][0]
+            type2 = types[1][0]
+            epsilon = _convertParameterToNumber(parameters['epsilon'])
+            sigma = _convertParameterToNumber(parameters['sigma'])
+            self.nbfix_types[(type1, type2)] = [sigma, epsilon]
+            self.nbfix_types[(type2, type1)] = [sigma, epsilon]
+
+    def registerLennardJones(self, parameters):
+        self.lj_types.registerAtom(parameters)
+
+
+    @staticmethod
+    def parseElement(element, ff):
+        existing = [f for f in ff._forces if isinstance(f, LennardJonesGenerator)]
+        if len(existing) == 0:
+            generator = LennardJonesGenerator(ff, float(element.attrib['lj14scale']))
+            ff.registerGenerator(generator)
+        else:
+            # Multiple <LennardJonesForce> tags were found, probably in different files
+            generator = existing[0]
+
+        for LJ in element.findall('Atom'):
+            generator.registerLennardJones(LJ.attrib)
+        for Nbfix in element.findall('NBFixPair'):
+            generator.registerNBFIX(Nbfix.attrib)
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+
+        # Ewald and PME will be interpreted as CutoffPeriodic for the CustomNonbondedForce
+
+        # We need a CustomNonbondedForce to implement the NBFIX functionality.
+        nbfix_type_set = set().union(*self.nbfix_types)
+        # First derive the lookup tables
+        lj_indx_list = [None for atom in data.atoms]
+        num_lj_types = 0
+        lj_type_list = []
+        type_map = {}
+        for i, atom in enumerate(data.atoms):
+            atype = data.atomType[atom]
+            values = (self.lj_types.paramsForType[atype]['sigma'], self.lj_types.paramsForType[atype]['epsilon'])
+            if lj_indx_list[i] is not None: continue # ljtype has been assigned an index
+            if values in type_map and atype not in nbfix_type_set:
+                # only non NBFIX types can be compressed
+                lj_indx_list[i] = type_map[values]
+            else:
+                type_map[values] = num_lj_types
+                lj_indx_list[i] = num_lj_types
+                num_lj_types += 1
+                lj_type_list.append(atype)
+        reverse_map = [0 for type_values in range(len(type_map))]
+        for type_value in type_map:
+            reverse_map[type_map[type_value]] = type_value
+
+        # Now everything is assigned. Create the A- and B-coefficient arrays
+        acoef = [0 for i in range(num_lj_types*num_lj_types)]
+        bcoef = acoef[:]
+        for m in range(num_lj_types):
+            for n in range(num_lj_types):
+                pair = (lj_type_list[m], lj_type_list[n])
+                if pair in self.nbfix_types:
+                    wdij = self.nbfix_types[pair][1]
+                    rij = self.nbfix_types[pair][0] * 2**(1.0/6) / 2
+                    rij6 = rij**6
+                    acoef[m+num_lj_types*n] = wdij * rij6**2
+                    bcoef[m+num_lj_types*n] = 2 * wdij * rij6
+                    continue
+                else:
+                    rij = (reverse_map[m][0] + reverse_map[n][0]) * 2**(1.0/6) / 2
+                    rij6 = rij**6
+                    wdij = math.sqrt(reverse_map[m][-1] * reverse_map[n][-1])
+                    acoef[m+num_lj_types*n] = wdij * rij6**2
+                    bcoef[m+num_lj_types*n] = 2 * wdij * rij6
+
+        self.force = mm.CustomNonbondedForce('a/r^12-b/r^6;'
+                                'a=acoef(type1, type2);'
+                                'b=bcoef(type1, type2)')
+        self.force.addTabulatedFunction('acoef',
+                mm.Discrete2DFunction(num_lj_types, num_lj_types, acoef))
+        self.force.addTabulatedFunction('bcoef',
+                mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef))
+        self.force.addPerParticleParameter('type')
+        if (nonbondedMethod is PME or nonbondedMethod is Ewald or
+                nonbondedMethod is CutoffPeriodic):
+            self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
+        elif nonbondedMethod is NoCutoff:
+            self.force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff)
+        elif nonbondedMethod is CutoffNonPeriodic:
+            self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
+        else:
+            raise AssertionError('Unrecognized nonbonded method [%s]' %
+                                 nonbondedMethod)
+        # Add the particles
+        for i in lj_indx_list:
+            self.force.addParticle((i,))
+        self.force.setUseLongRangeCorrection(True)
+        self.force.setCutoffDistance(nonbondedCutoff)
+
+        sys.addForce(self.force)
+
+    def postprocessSystem(self, sys, data, args):
+
+        bondIndices = []
+        for bond in data.bonds:
+            bondIndices.append((bond.atom1, bond.atom2))
+
+        # If a virtual site does *not* share exclusions with another atom, add a bond between it and its first parent
+        # atom.
+        for i in range(sys.getNumParticles()):
+            if sys.isVirtualSite(i):
+                (site, atoms, excludeWith) = data.virtualSites[data.atoms[i]]
+                if excludeWith is None:
+                    bondIndices.append((i, site.getParticle(0)))
+
+        # Certain particles, such as lone pairs and Drude particles, share exclusions with a parent atom.
+        # If the parent atom does not interact with an atom, the child particle does not either.
+
+        for atom1, atom2 in bondIndices:
+            for child1 in data.excludeAtomWith[atom1]:
+                bondIndices.append((child1, atom2))
+                for child2 in data.excludeAtomWith[atom2]:
+                    bondIndices.append((child1, child2))
+            for child2 in data.excludeAtomWith[atom2]:
+                bondIndices.append((atom1, child2))
+
+        # Create the exceptions.
+        self.force.createExclusionsFromBonds(bondIndices, 3)
+
+parsers["LennardJonesForce"] = LennardJonesGenerator.parseElement
 
 ## @private
 class GBSAOBCGenerator(object):

wrappers/python/tests/TestForceField.py

@@ -11,6 +11,7 @@ try:
 except ImportError:
     from io import StringIO
 import os
+import warnings
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
@@ -712,5 +713,66 @@ class AmoebaTestForceField(unittest.TestCase):
             diff = norm(f1-f2)
             self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
 
+    def test_LennardJones_generator(self):
+        """ Test the LennardJones generator"""
+        warnings.filterwarnings('ignore', category=CharmmPSFWarning)
+        psf = CharmmPsfFile('systems/ions.psf')
+        pdb = PDBFile('systems/ions.pdb')
+        params = CharmmParameterSet('systems/toppar_water_ions.str'
+                                    )
+
+        # Box dimensions (found from bounding box)
+        psf.setBox(12.009*angstroms,   12.338*angstroms,   11.510*angstroms)
+
+        # Turn off charges so we only test the Lennard-Jones energies
+        for a in psf.atom_list:
+            a.charge = 0.0
+
+        # Now compute the full energy
+        plat = Platform.getPlatformByName('Reference')
+        system = psf.createSystem(params, nonbondedMethod=PME,
+                                  nonbondedCutoff=5*angstroms)
+
+        con = Context(system, VerletIntegrator(2*femtoseconds), plat)
+        con.setPositions(pdb.positions)
+
+        # Now set up stystem from ffxml. Setting chareges to 0 so we only test the Lennard-Jones energies
+        xml = """
+<ForceField>
+ <AtomTypes>
+  <Type name="SOD" class="SOD" element="Na" mass="22.98977"/>
+  <Type name="CLA" class="CLA" element="Cl" mass="35.45"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="CLA">
+   <Atom name="CLA" type="CLA" charge="0.0"/>
+  </Residue>
+  <Residue name="SOD">
+   <Atom name="SOD" type="SOD" charge="0.0"/>
+  </Residue>
+ </Residues>
+ <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
+  <UseAttributeFromResidue name="charge"/>
+  <Atom type="SOD" sigma="1.0" epsilon="0.0"/>
+  <Atom type="CLA" sigma="1.0" epsilon="0.0"/>
+ </NonbondedForce>
+ <LennardJonesForce lj14scale="1.0">
+  <Atom type="CLA" sigma="0.404468018036" epsilon="0.6276"/>
+  <Atom type="SOD" sigma="0.251367073323" epsilon="0.1962296"/>
+  <NBFixPair type1="CLA" type2="SOD" sigma="0.664788623476" epsilon="0.350933"/>
+ </LennardJonesForce>
+</ForceField> """
+        ff = ForceField(StringIO(xml))
+        system2 = ff.createSystem(pdb.topology, nonbondedMethod=PME,
+                                  nonbondedCutoff=5*angstroms)
+        con2 = Context(system2, VerletIntegrator(2*femtoseconds), plat)
+        con2.setPositions(pdb.positions)
+
+        state = con.getState(getEnergy=True, enforcePeriodicBox=True)
+        ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        state2 = con2.getState(getEnergy=True, enforcePeriodicBox=True)
+        ene2 = state2.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        self.assertAlmostEqual(ene, ene2)
+
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/systems/ions.pdb

+CRYST1   12.009   12.338   11.510  90.00  90.00  90.00 P 1
+ATOM      1  SOD SOD I   1      -5.844   1.432   3.239  1.00  0.00      ION NA
+ATOM      2  CLA CLA I   2      -5.605  -3.153   1.213  1.00  0.00      ION CL
+END
\ No newline at end of file

wrappers/python/tests/systems/ions.psf

+PSF
+
+       4 !NTITLE
+ REMARKS original generated structure x-plor psf file
+ REMARKS topology ../../../../../Charmm36_FF/toppar/top_all36_cgenff.rtf 
+ REMARKS topology ../../../../../Charmm36_FF/toppar/toppar_water_ions.str 
+ REMARKS segment ION { first NONE; last NONE; auto angles dihedrals }
+
+       2 !NATOM
+       1 ION  1    SOD  SOD  SOD    1.000000       22.9898           0
+       2 ION  2    CLA  CLA  CLA   -1.000000       35.4500           0
+
+       0 !NBOND: bonds
+
+
+       0 !NTHETA: angles
+
+
+       0 !NPHI: dihedrals
+
+
+       0 !NIMPHI: impropers
+
+
+       0 !NDON: donors
+
+
+       0 !NACC: acceptors
+
+
+       0 !NNB
+
+       0       0
+
+       1       0 !NGRP
+       0       0       0
+

wrappers/python/tests/systems/methanol_ions.pdb

+CRYST1   12.009   12.338   11.510  90.00  90.00  90.00 P 1           1
+ATOM      1  CB  MEOHM   1      -0.748  -0.015   0.024  1.00  0.00      METH C
+ATOM      2  OG  MEOHM   1       0.558   0.420  -0.278  1.00  0.00      METH O
+ATOM      3  HG1 MEOHM   1       0.716   1.404   0.137  1.00  0.00      METH H
+ATOM      4  HB1 MEOHM   1      -1.293  -0.202  -0.901  1.00  0.00      METH H
+ATOM      5  HB2 MEOHM   1      -1.263   0.754   0.600  1.00  0.00      METH H
+ATOM      6  HB3 MEOHM   1      -0.699  -0.934   0.609  1.00  0.00      METH H
+ATOM      7  SOD SOD I   2      -5.844   1.432   3.239  1.00  0.00      ION NA
+ATOM      8  CLA CLA I   3      -5.605  -3.153   1.213  1.00  0.00      ION CL
+CONECT 1    2    4     5     6
+CONECT 2    3
+END
+END
\ No newline at end of file

wrappers/python/tests/systems/methanol_ions.psf

+PSF EXT
+
+       4 !NTITLE
+ REMARKS original generated structure x-plor psf file
+ REMARKS topology ../../../../../Charmm36_FF/toppar/top_all36_cgenff.rtf 
+ REMARKS topology ../../../../../Charmm36_FF/toppar/toppar_water_ions.str 
+ REMARKS segment ION { first NONE; last NONE; auto angles dihedrals }
+
+       8 !NATOM
+         1 METH     1        MEOH     CB       CG331   -0.040000       12.0110           0
+         2 METH     1        MEOH     OG       OG311   -0.650000       15.9994           0
+         3 METH     1        MEOH     HG1      HGP1     0.420000        1.0080           0
+         4 METH     1        MEOH     HB1      HGA3     0.090000        1.0080           0
+         5 METH     1        MEOH     HB2      HGA3     0.090000        1.0080           0
+         6 METH     1        MEOH     HB3      HGA3     0.090000        1.0080           0
+         7 ION      2        SOD      SOD      SOD      1.000000       22.9898           0
+         8 ION      3        CLA      CLA      CLA     -1.000000       35.4500           0
+
+       5 !NBOND: bonds
+         1         2         1         4         1         5         1         6
+         2         3
+
+       7 !NTHETA: angles
+         1         2         3         2         1         6         2         1         5
+         2         1         4         4         1         6         4         1         5
+         5         1         6
+
+       3 !NPHI: dihedrals
+         3         2         1         4         3         2         1         5
+         3         2         1         6
+
+       0 !NIMPHI: impropers
+
+
+       0 !NDON: donors
+
+
+       0 !NACC: acceptors
+
+
+       0 !NNB
+
+         0         0         0         0         0         0         0         0
+
+
+       1       0 !NGRP
+         0         0         0
+