mirror of
https://github.com/openmm/openmm
synced 2026-06-03 06:39:48 +09:00
e62bdf6adb
* Can use getPlatform() instead of getPlatformByName() * More concise arguments for getState()
222 lines
10 KiB
Python
222 lines
10 KiB
Python
import unittest
|
|
import warnings
|
|
import tempfile
|
|
from datetime import datetime, timedelta
|
|
from openmm import *
|
|
from openmm.app import *
|
|
from openmm.unit import *
|
|
import math, random
|
|
|
|
class TestIntegrators(unittest.TestCase):
|
|
"""Test Python Integrator classes"""
|
|
|
|
def testMTSIntegratorExplicit(self):
|
|
"""Test the MTS integrator on an explicit solvent system"""
|
|
# Create a periodic solvated system with PME
|
|
pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
|
|
ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
|
|
system = ff.createSystem(pdb.topology, nonbondedMethod=PME)
|
|
|
|
# Split forces into groups
|
|
for force in system.getForces():
|
|
if force.__class__.__name__ == 'NonbondedForce':
|
|
force.setForceGroup(1)
|
|
force.setReciprocalSpaceForceGroup(2)
|
|
else:
|
|
force.setForceGroup(0)
|
|
|
|
# Create an integrator
|
|
integrator = MTSIntegrator(4*femtoseconds, [(2,1), (1,2), (0,8)])
|
|
|
|
# Run a few steps of dynamics
|
|
context = Context(system, integrator)
|
|
context.setPositions(pdb.positions)
|
|
integrator.step(10)
|
|
|
|
# Ensure energy is well-behaved.
|
|
state = context.getState(getEnergy=True)
|
|
if not (state.getPotentialEnergy() / kilojoules_per_mole < 0.0):
|
|
raise Exception('Potential energy of alanine dipeptide system with MTS integrator is blowing up: %s' % str(state.getPotentialEnergy()))
|
|
|
|
def testMTSIntegratorConstraints(self):
|
|
"""Test the MTS integrator energy conservation on a system of constrained particles with no inner force (just constraints)"""
|
|
|
|
# Create a constrained test system
|
|
numParticles = 8
|
|
numConstraints = 5
|
|
system = System()
|
|
force = NonbondedForce()
|
|
for i in range(numParticles):
|
|
system.addParticle(5.0 if i%2==0 else 10.0)
|
|
force.addParticle((0.2 if i%2==0 else -0.2), 0.5, 5.0);
|
|
system.addConstraint(0, 1, 1.0);
|
|
system.addConstraint(1, 2, 1.0);
|
|
system.addConstraint(2, 3, 1.0);
|
|
system.addConstraint(4, 5, 1.0);
|
|
system.addConstraint(6, 7, 1.0);
|
|
system.addForce(force)
|
|
|
|
# Create integrator where inner timestep just evaluates constraints
|
|
integrator = MTSIntegrator(1*femtoseconds, [(1,1), (0,4)])
|
|
integrator.setConstraintTolerance(1e-5);
|
|
|
|
positions = [ (i/2., (i+1)/2., 0.) for i in range(numParticles) ]
|
|
velocities = [ (random.random()-0.5, random.random()-0.5, random.random()-0.5) for i in range(numParticles) ]
|
|
|
|
# Create Context
|
|
platform = Platform.getPlatform('Reference')
|
|
context = Context(system, integrator, platform)
|
|
context.setPositions(positions)
|
|
context.setVelocities(velocities)
|
|
context.applyConstraints(1e-5)
|
|
|
|
# Simulate it and see whether the constraints remain satisfied.
|
|
CONSTRAINT_RELATIVE_TOLERANCE = 1.e-4 # relative constraint violation tolerance
|
|
ENERGY_RELATIVE_TOLERANCE = 1.e-2 # relative energy violation tolerance
|
|
for i in range(1000):
|
|
state = context.getState(getPositions=True, getEnergy=True)
|
|
positions = state.getPositions()
|
|
for j in range(numConstraints):
|
|
[particle1, particle2, constraint_distance] = system.getConstraintParameters(j)
|
|
current_distance = 0.0 * nanometers**2
|
|
for k in range(3):
|
|
current_distance += (positions[particle1][k] - positions[particle2][k])**2
|
|
current_distance = sqrt(current_distance)
|
|
# Fail test if outside of relative tolerance
|
|
relative_violation = (current_distance - constraint_distance) / constraint_distance
|
|
if (relative_violation > CONSTRAINT_RELATIVE_TOLERANCE):
|
|
raise Exception('Constrained distance is violated by relative tolerance of %f (constraint %s actual %s)' % (relative_violation, str(constraint_distance), str(current_distance)))
|
|
# Check total energy
|
|
total_energy = state.getPotentialEnergy() + state.getKineticEnergy()
|
|
if (i == 1):
|
|
initial_energy = total_energy
|
|
elif (i > 1):
|
|
relative_violation = abs((total_energy - initial_energy) / initial_energy)
|
|
if (relative_violation > ENERGY_RELATIVE_TOLERANCE):
|
|
raise Exception('Total energy is violated by relative tolerance of %f on step %d (initial %s final %s)' % (relative_violation, i, str(initial_energy), str(total_energy)))
|
|
# Take a step
|
|
integrator.step(1)
|
|
|
|
def testBadGroups(self):
|
|
"""Test the MTS integrator with bad force group substeps."""
|
|
# Create a periodic solvated system with PME
|
|
pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
|
|
ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
|
|
system = ff.createSystem(pdb.topology, nonbondedMethod=PME)
|
|
|
|
# Split forces into groups
|
|
for force in system.getForces():
|
|
if force.__class__.__name__ == 'NonbondedForce':
|
|
force.setForceGroup(1)
|
|
force.setReciprocalSpaceForceGroup(2)
|
|
else:
|
|
force.setForceGroup(0)
|
|
|
|
with self.assertRaises(ValueError):
|
|
# Create an integrator
|
|
integrator = MTSIntegrator(4*femtoseconds, [(2,1), (1,3), (0,8)])
|
|
|
|
# Run a few steps of dynamics
|
|
context = Context(system, integrator)
|
|
context.setPositions(pdb.positions)
|
|
integrator.step(10)
|
|
|
|
def testMTSLangevinIntegrator(self):
|
|
"""Test the MTSLangevinIntegrator on an explicit solvent system"""
|
|
# Create a periodic solvated system with PME
|
|
pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
|
|
ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
|
|
system = ff.createSystem(pdb.topology, nonbondedMethod=PME)
|
|
|
|
# Split forces into groups
|
|
for force in system.getForces():
|
|
if force.__class__.__name__ == 'NonbondedForce':
|
|
force.setForceGroup(1)
|
|
force.setReciprocalSpaceForceGroup(2)
|
|
else:
|
|
force.setForceGroup(0)
|
|
|
|
# Create an integrator
|
|
integrator = MTSLangevinIntegrator(300*kelvin, 5/picosecond, 4*femtoseconds, [(2,1), (1,2), (0,4)])
|
|
|
|
# Run some equilibration.
|
|
context = Context(system, integrator)
|
|
context.setPositions(pdb.positions)
|
|
context.setVelocitiesToTemperature(300*kelvin)
|
|
integrator.step(500)
|
|
|
|
# See if the temperature is correct.
|
|
totalEnergy = 0*kilojoules_per_mole
|
|
steps = 50
|
|
for i in range(steps):
|
|
integrator.step(10)
|
|
totalEnergy += context.getState(getEnergy=True).getKineticEnergy()
|
|
averageEnergy = totalEnergy/steps
|
|
dofs = 3*system.getNumParticles() - system.getNumConstraints() - 3
|
|
temperature = averageEnergy*2/(dofs*MOLAR_GAS_CONSTANT_R)
|
|
self.assertTrue(290*kelvin < temperature < 310*kelvin)
|
|
|
|
def testMTSLangevinIntegratorFriction(self):
|
|
"""Test the MTSLangevinIntegrator on a force-free particle to ensure friction is properly accounted for (issue #3790)"""
|
|
# Create a System with a single particle and no forces
|
|
system = System()
|
|
system.addParticle(12.0*amu)
|
|
platform = Platform.getPlatform('Reference')
|
|
initial_positions = [Vec3(0,0,0)]
|
|
initial_velocities = [Vec3(1,0,0)]
|
|
nsteps = 125 # number of steps to take
|
|
collision_rate = 1/picosecond
|
|
timestep = 4*femtoseconds
|
|
|
|
def get_final_velocities(nsubsteps):
|
|
"""Get the final velocity vector after a fixed number of steps for the specified number of substeps"""
|
|
integrator = MTSLangevinIntegrator(0*kelvin, collision_rate, timestep, [(0,nsubsteps)])
|
|
context = Context(system, integrator, platform)
|
|
context.setPositions(initial_positions)
|
|
context.setVelocities(initial_velocities)
|
|
integrator.step(nsteps)
|
|
final_velocities = context.getState(getVelocities=True).getVelocities()
|
|
del context, integrator
|
|
return final_velocities
|
|
|
|
# Compare sub-stepped MTS with single-step MTS
|
|
for nsubsteps in range(2,6):
|
|
mts_velocities = get_final_velocities(nsubsteps)
|
|
self.assertAlmostEqual(math.exp(-timestep*nsteps*collision_rate), mts_velocities[0].x)
|
|
self.assertAlmostEqual(0, mts_velocities[0].y)
|
|
self.assertAlmostEqual(0, mts_velocities[0].z)
|
|
|
|
def testNoseHooverIntegrator(self):
|
|
"""Test partial thermostating in the NoseHooverIntegrator (only API)"""
|
|
pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
|
|
ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
|
|
system = ff.createSystem(pdb.topology, nonbondedMethod=PME)
|
|
|
|
integrator = NoseHooverIntegrator(1.0*femtosecond)
|
|
integrator.addSubsystemThermostat(list(range(5)), [], 200*kelvin, 1/picosecond, 200*kelvin, 1/picosecond, 3,3,3)
|
|
con = Context(system, integrator)
|
|
con.setPositions(pdb.positions)
|
|
|
|
integrator.step(5)
|
|
self.assertNotEqual(integrator.computeHeatBathEnergy(), 0.0*kilojoule_per_mole)
|
|
|
|
def testDrudeNoseHooverIntegrator(self):
|
|
"""Test the DrudeNoseHooverIntegrator"""
|
|
warnings.filterwarnings('ignore', category=CharmmPSFWarning)
|
|
psf = CharmmPsfFile('systems/ala3_solv_drude.psf')
|
|
crd = CharmmCrdFile('systems/ala3_solv_drude.crd')
|
|
params = CharmmParameterSet('systems/toppar_drude_master_protein_2013e.str')
|
|
# Box dimensions (cubic box)
|
|
psf.setBox(33.2*angstroms, 33.2*angstroms, 33.2*angstroms)
|
|
|
|
system = psf.createSystem(params, nonbondedMethod=PME, ewaldErrorTolerance=0.0005)
|
|
integrator = DrudeNoseHooverIntegrator(300*kelvin, 1.0/picosecond, 1*kelvin, 10/picosecond, 0.001*picoseconds)
|
|
con = Context(system, integrator)
|
|
con.setPositions(crd.positions)
|
|
|
|
integrator.step(5)
|
|
self.assertNotEqual(integrator.computeHeatBathEnergy(), 0.0*kilojoule_per_mole)
|
|
|
|
if __name__ == '__main__':
|
|
unittest.main()
|