TestGromacsTopFile.py 14.3 KB
Newer Older
1
2
import unittest
from validateConstraints import *
3
4
5
6
7
from openmm.app import *
from openmm import *
from openmm.unit import *
from openmm.app.gromacstopfile import _defaultGromacsIncludeDir
import openmm.app.element as elem
8
from numpy.testing import assert_allclose
9

Robert McGibbon's avatar
Robert McGibbon committed
10
11
12
GROMACS_INCLUDE = _defaultGromacsIncludeDir()

@unittest.skipIf(not os.path.exists(GROMACS_INCLUDE), 'GROMACS is not installed')
13
14
15
class TestGromacsTopFile(unittest.TestCase):

    """Test the GromacsTopFile.createSystem() method."""
Robert McGibbon's avatar
Robert McGibbon committed
16

17
18
19
20
21
22
23
    def setUp(self):
        """Set up the tests by loading the input files."""

        # alanine dipeptide with explicit water
        self.top1 = GromacsTopFile('systems/explicit.top', unitCellDimensions=Vec3(6.223, 6.223, 6.223)*nanometers)

    def test_NonbondedMethod(self):
Peter Eastman's avatar
Peter Eastman committed
24
        """Test all six options for the nonbondedMethod parameter."""
25

Robert McGibbon's avatar
Robert McGibbon committed
26
27
28
        methodMap = {NoCutoff:NonbondedForce.NoCutoff,
                     CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic,
                     CutoffPeriodic:NonbondedForce.CutoffPeriodic,
Peter Eastman's avatar
Peter Eastman committed
29
30
31
                     Ewald:NonbondedForce.Ewald,
                     PME:NonbondedForce.PME,
                     LJPME:NonbondedForce.LJPME}
32
33
34
        for method in methodMap:
            system = self.top1.createSystem(nonbondedMethod=method)
            forces = system.getForces()
Robert McGibbon's avatar
Robert McGibbon committed
35
36
            self.assertTrue(any(isinstance(f, NonbondedForce) and
                                f.getNonbondedMethod()==methodMap[method]
37
38
                                for f in forces))

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
    def test_LJPME_mixing(self):
        """Test that LJPME is not allowed with non-standard LJ mixing or NBFIX."""

        tests = {
           'apgr.nbfix.pairs.comb1.top': False,
           'apgr.nbfix.pairs.comb2.top': False,
           'apgr.nbfix.pairs.comb3.top': False,
           'apgr.nbfix.nopairs.comb1.top': False,
           'apgr.nbfix.nopairs.comb2.top': False,
           'apgr.nbfix.nopairs.comb3.top': False,
           'apgr.nonbfix.pairs.comb1.top': False,
           'apgr.nonbfix.pairs.comb2.top': True,
           'apgr.nonbfix.pairs.comb3.top': False,
           'apgr.nonbfix.nopairs.comb1.top': False,
           'apgr.nonbfix.nopairs.comb2.top': True,
           'apgr.nonbfix.nopairs.comb3.top': False,
        }

        for topfile, ljpme_allowed in tests.items():
            top = GromacsTopFile(f'systems/{topfile}', unitCellDimensions=Vec3(30, 30, 30)*angstroms)
            if ljpme_allowed:
                system = top.createSystem(nonbondedMethod=LJPME)
                forces = system.getForces()
                self.assertTrue(any(isinstance(f, NonbondedForce) and
                                    f.getNonbondedMethod()==NonbondedForce.LJPME
                                    for f in forces))
            else:
                with self.assertRaisesRegex(ValueError, 'LJPME is not supported'):
                    top.createSystem(nonbondedMethod=LJPME)

69
70
71
72
73
74
75
76
77
    def test_ff99SBILDN(self):
        """ Test Gromacs topology #define replacement as used in ff99SB-ILDN """
        top = GromacsTopFile('systems/aidilnaaaaa.top')
        gro = GromacsGroFile('systems/aidilnaaaaa.gro')
        system = top.createSystem()
        for force in system.getForces():
            if isinstance(force, PeriodicTorsionForce):
                force.setForceGroup(1)
        context = Context(system, VerletIntegrator(1*femtosecond),
Peter Eastman's avatar
Peter Eastman committed
78
                          Platform.getPlatform('Reference'))
79
80
81
82
        context.setPositions(gro.positions)
        ene = context.getState(getEnergy=True, groups=2**1).getPotentialEnergy()
        self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), 341.6905133582857)

Stephen Constable's avatar
Stephen Constable committed
83
    def test_SMOG(self):
Stephen Constable's avatar
Stephen Constable committed
84
        """ Test to ensure that SMOG models can be run without problems """
Stephen Constable's avatar
Stephen Constable committed
85
86
87
88
89
        top = GromacsTopFile('systems/2ci2.pdb.top')
        gro = GromacsGroFile('systems/2ci2.pdb.gro')
        system = top.createSystem()

        context = Context(system, VerletIntegrator(1*femtosecond),
Peter Eastman's avatar
Peter Eastman committed
90
                          Platform.getPlatform('Reference'))
Stephen Constable's avatar
Stephen Constable committed
91
92
93
94
        context.setPositions(gro.positions)
        ene = context.getState(getEnergy=True).getPotentialEnergy()
        self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), -346.940915296)

95
96
97
98
99
100
101
102
103
104
    def test_ionic(self):
        """Test simulating an ionic liquid"""
        gro = GromacsGroFile('systems/ionic.gro')
        top = GromacsTopFile('systems/ionic.top', periodicBoxVectors=gro.getPeriodicBoxVectors())
        system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1.2)
        for f in system.getForces():
            if isinstance(f, CustomNonbondedForce):
                f.setUseLongRangeCorrection(True)

        context = Context(system, VerletIntegrator(1*femtosecond),
Peter Eastman's avatar
Peter Eastman committed
105
                          Platform.getPlatform('Reference'))
106
107
108
109
110
        context.setPositions(gro.positions)
        energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
        self.assertAlmostEqual(energy, 3135.33, delta=energy*0.005)
        self.assertEqual(1400, system.getNumConstraints())

111
112
113
    def test_Cutoff(self):
        """Test to make sure the nonbondedCutoff parameter is passed correctly."""

Peter Eastman's avatar
Peter Eastman committed
114
        for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, LJPME]:
Robert McGibbon's avatar
Robert McGibbon committed
115
116
            system = self.top1.createSystem(nonbondedMethod=method,
                                               nonbondedCutoff=2*nanometer,
117
118
119
120
121
122
123
124
                                               constraints=HBonds)
            cutoff_distance = 0.0*nanometer
            cutoff_check = 2.0*nanometer
            for force in system.getForces():
                if isinstance(force, NonbondedForce):
                    cutoff_distance = force.getCutoffDistance()
            self.assertEqual(cutoff_distance, cutoff_check)

125
126
    def test_SwitchingFunction(self):
        """Test using a switching function."""
127
128
129
130
131
132
133
134
135
136
        top = GromacsTopFile('systems/ionic.top')
        for distance in (None, 0.8*nanometers):
            system = top.createSystem(nonbondedMethod=CutoffNonPeriodic, switchDistance=distance)
            for f in system.getForces():
                if isinstance(f, NonbondedForce) or isinstance(f, CustomNonbondedForce):
                    if distance is None:
                        self.assertFalse(f.getUseSwitchingFunction())
                    else:
                        self.assertTrue(f.getUseSwitchingFunction())
                        self.assertEqual(distance, f.getSwitchingDistance())
137

138
139
140
    def test_EwaldErrorTolerance(self):
        """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""

Peter Eastman's avatar
Peter Eastman committed
141
        for method in [Ewald, PME, LJPME]:
Robert McGibbon's avatar
Robert McGibbon committed
142
143
            system = self.top1.createSystem(nonbondedMethod=method,
                                               ewaldErrorTolerance=1e-6,
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
                                               constraints=HBonds)
            tolerance = 0
            tolerance_check = 1e-6
            for force in system.getForces():
                if isinstance(force, NonbondedForce):
                    tolerance = force.getEwaldErrorTolerance()
            self.assertEqual(tolerance, tolerance_check)

    def test_RemoveCMMotion(self):
        """Test both options (True and False) for the removeCMMotion parameter."""

        for b in [True, False]:
            system = self.top1.createSystem(removeCMMotion=b)
            self.assertEqual(any(isinstance(f, CMMotionRemover) for f in system.getForces()), b)

    def test_RigidWaterAndConstraints(self):
        """Test all eight options for the constraints and rigidWater parameters."""

        topology = self.top1.topology
        for constraints_value in [None, HBonds, AllBonds, HAngles]:
            for rigidWater_value in [True, False]:
Robert McGibbon's avatar
Robert McGibbon committed
165
                system = self.top1.createSystem(constraints=constraints_value,
166
                                                   rigidWater=rigidWater_value)
Robert McGibbon's avatar
Robert McGibbon committed
167
                validateConstraints(self, topology, system,
168
169
                                    constraints_value, rigidWater_value)

170
171
    def test_HydrogenMass(self):
        """Test that altering the mass of hydrogens works correctly."""
Robert McGibbon's avatar
Robert McGibbon committed
172

173
174
175
176
177
178
179
        topology = self.top1.topology
        hydrogenMass = 4*amu
        system1 = self.top1.createSystem()
        system2 = self.top1.createSystem(hydrogenMass=hydrogenMass)
        for atom in topology.atoms():
            if atom.element == elem.hydrogen:
                self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
180
181
182
183
                if atom.residue.name == 'HOH':
                    self.assertEqual(system1.getParticleMass(atom.index), system2.getParticleMass(atom.index))
                else:
                    self.assertEqual(hydrogenMass, system2.getParticleMass(atom.index))
184
185
186
187
        totalMass1 = sum([system1.getParticleMass(i) for i in range(system1.getNumParticles())]).value_in_unit(amu)
        totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
        self.assertAlmostEqual(totalMass1, totalMass2)

188
189
190
191
192
    def test_VirtualParticle(self):
        """Test virtual particle works correctly."""

        top = GromacsTopFile('systems/bnz.top')
        gro = GromacsGroFile('systems/bnz.gro')
193
194
195
196
197
198
199
        for atom in top.topology.atoms():
            if atom.name.startswith('C'):
                self.assertEqual(elem.carbon, atom.element)
            elif atom.name.startswith('H'):
                self.assertEqual(elem.hydrogen, atom.element)
            else:
                self.assertIsNone(atom.element)
200
201
202
203
204
205
        system = top.createSystem()

        self.assertEqual(26, system.getNumParticles())
        self.assertEqual(1, len(top._moleculeTypes['BENX'].vsites2))

        context = Context(system, VerletIntegrator(1*femtosecond),
Peter Eastman's avatar
Peter Eastman committed
206
                          Platform.getPlatform('Reference'))
207
        context.setPositions(gro.positions)
Sunhwan Jo's avatar
Sunhwan Jo committed
208
        context.computeVirtualSites()
209
210
        ene = context.getState(getEnergy=True).getPotentialEnergy()
        # the energy output is from gromacs and it only prints out 6 sig digits.
Sunhwan Jo's avatar
Sunhwan Jo committed
211
        self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), 1.88855e+02, places=3)
212

213
    def test_Vsite3Func1(self):
214
215
216
        """Test a three particle virtual site."""
        top = GromacsTopFile('systems/tip4pew.top')
        system = top.createSystem()
217
        self.assertEqual(3, system.getNumConstraints())
218
219
220
221
222
223
224
225
226
227
        self.assertTrue(system.isVirtualSite(3))
        vs = system.getVirtualSite(3)
        self.assertIsInstance(vs, ThreeParticleAverageSite)
        self.assertEqual(0, vs.getParticle(0))
        self.assertEqual(1, vs.getParticle(1))
        self.assertEqual(2, vs.getParticle(2))
        self.assertAlmostEqual(0.786646558, vs.getWeight(0))
        self.assertAlmostEqual(0.106676721, vs.getWeight(1))
        self.assertAlmostEqual(0.106676721, vs.getWeight(2))

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    def test_Vsite3Func4(self):
        """Test a three particle virtual site."""
        top = GromacsTopFile('systems/tip5p.top')
        system = top.createSystem()
        self.assertEqual(3, system.getNumConstraints())
        for i in (3, 4):
            self.assertTrue(system.isVirtualSite(i))
            vs = system.getVirtualSite(i)
            self.assertIsInstance(vs, OutOfPlaneSite)
            self.assertEqual(0, vs.getParticle(0))
            self.assertEqual(1, vs.getParticle(1))
            self.assertEqual(2, vs.getParticle(2))
            self.assertAlmostEqual(-0.344908, vs.getWeight12())
            self.assertAlmostEqual(-0.344908, vs.getWeight13())
            wc = -6.4437903493
            if i == 4:
                wc = -wc
            self.assertAlmostEqual(wc, vs.getWeightCross())
246
    
247
248
249
250
251
252
253
254
    def test_GROMOS(self):
        """Test a system using the GROMOS 54a7 force field."""

        top = GromacsTopFile('systems/1ppt.top')
        gro = GromacsGroFile('systems/1ppt.gro')
        system = top.createSystem()
        for i, f in enumerate(system.getForces()):
            f.setForceGroup(i)
Peter Eastman's avatar
Peter Eastman committed
255
        context = Context(system, VerletIntegrator(1*femtosecond), Platform.getPlatform('Reference'))
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
        context.setPositions(gro.positions)
        energy = {}
        for i, f in enumerate(system.getForces()):
            energy[f.getName()] = context.getState(getEnergy=True, groups={i}).getPotentialEnergy().value_in_unit(kilojoules_per_mole)

        # Compare to energies computed with GROMACS.

        assert_allclose(1.12797e+03, energy['GROMOSBondForce'], rtol=1e-4)
        assert_allclose(5.59066e+02, energy['GROMOSAngleForce'], rtol=1e-4)
        assert_allclose(3.80152e+02, energy['PeriodicTorsionForce'], rtol=1e-4)
        assert_allclose(9.59178e+01, energy['HarmonicTorsionForce'], rtol=1e-4)
        assert_allclose(2.75307e+02, energy['LennardJonesExceptions'], rtol=1e-4)
        assert_allclose(-7.53704e+02, energy['LennardJonesForce'], rtol=1e-4)
        assert_allclose(-6.23055e+03+4.36880e+03, energy['NonbondedForce'], rtol=1e-4)
        total = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
        assert_allclose(-1.77020e+02, total, rtol=1e-3)

273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
    def test_NBFIX(self):
        """Test systems using NBFIX with and without pairtypes and different combination rules."""

        gro = GromacsGroFile('systems/apgr.gro')

        energies = {
           'apgr.nbfix.pairs.comb1.top': -986.713,
           'apgr.nbfix.pairs.comb2.top': -986.485,
           'apgr.nbfix.pairs.comb3.top': -986.713,
           'apgr.nbfix.nopairs.comb1.top': -912.181,
           'apgr.nbfix.nopairs.comb2.top': -903.437,
           'apgr.nbfix.nopairs.comb3.top': -912.181,
           'apgr.nonbfix.pairs.comb1.top': -993.104,
           'apgr.nonbfix.pairs.comb2.top': -992.685,
           'apgr.nonbfix.pairs.comb3.top': -993.104,
           'apgr.nonbfix.nopairs.comb1.top': -918.572,
           'apgr.nonbfix.nopairs.comb2.top': -909.637,
290
           'apgr.nonbfix.nopairs.comb3.top': -918.572,
291
292
           'apgr.nbfix.nopairs.comb2.fudge.top': -856.721,
           'apgr.nbfix14.nopairs.comb2.top': -909.755
293
294
295
296
297
298
299
300
301
302
        }
        
        for topfile, expected in energies.items():   
           top = GromacsTopFile(f'systems/{topfile}')
           system = top.createSystem(nonbondedMethod=NoCutoff,switchDistance=None,constraints=None,useDispersionCorrection=False)
           context = Context(system, VerletIntegrator(1*femtosecond), Platform.getPlatform('Reference'))
           context.setPositions(gro.positions)
           energy=context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
           assert_allclose(energy, expected, rtol=1E-6)

303
304
305
if __name__ == '__main__':
    unittest.main()