search

openforcefield / openff-toolkit

The Open Forcefield Toolkit provides implementations of the SMIRNOFF format, parameterization engine, and other tools. Documentation available at http://open-forcefield-toolkit.readthedocs.io
http://openforcefield.org
MIT License
311 stars 91 forks source link

Add integration tests for hbond constraints #391

Open j-wags opened 5 years ago

j-wags commented 5 years ago

Now that hbond constraints are becoming an official part of the smirnoff99Frosst line of FFs (starting with 1.1.0), I put together a quick local integration test to ensure that they're being applied correctly as written. This integration test should be put into our CI, so we know if changes to OpenMM break our constraints.

This goes a step further than our existing constraint tests and actually runs a brief simulation to ensure that OpenMM interprets the constraints the way we want. This could easily be expanded with ParmEd round trips to/from different formats to ensure that constraints are transferred to other packages as well.

I think this test belongs in the toolkit repo (rather than the smirnoff99Frosst repo) since the toolkit is responsible for interpreting and applying the constraints.

Reproducing code here:

from openforcefield.topology import Molecule
from openforcefield.typing.engines.smirnoff import ForceField
from simtk import openmm, unit
import copy

# Create methane
mol = Molecule.from_smiles('C')
mol.generate_conformers(n_conformers=1)
off_top = mol.to_topology()
omm_top = off_top.to_openmm()

ff109 = ForceField('smirnoff99Frosst-1.0.9.offxml')
ff110 = ForceField('smirnoff99Frosst-1.1.0.offxml')

sys109 = ff109.create_openmm_system(off_top)
sys110 = ff110.create_openmm_system(off_top)

integrator109 = openmm.VerletIntegrator(1.0 * unit.femtoseconds)
integrator110 = copy.deepcopy(integrator109)

# Propagate the System with Langevin dynamics.
time_step = 2*unit.femtoseconds  # simulation timestep
temperature = 300*unit.kelvin  # simulation temperature
friction = 1/unit.picosecond  # collision rate
integrator = openmm.LangevinIntegrator(temperature, friction, time_step)

# Set up an OpenMM simulation.
sim109 = openmm.app.Simulation(omm_top, sys109, integrator109)
sim110 = openmm.app.Simulation(omm_top, sys110, integrator110)

# Set the initial positions.
positions = mol.conformers[0] 
sim109.context.setPositions(positions)
sim110.context.setPositions(positions)

# Randomize the velocities from a Boltzmann distribution at a given temperature.
sim109.context.setVelocitiesToTemperature(temperature)
sim110.context.setVelocitiesToTemperature(temperature)

#sim109.step(100)
#sim110.step(100)

sim109.step(1)
sim110.step(1)

pos109 = sim109.context.getState(getPositions=True).getPositions()
pos110 = sim110.context.getState(getPositions=True).getPositions()
ch_len109_0 = sum([(d/unit.angstrom)**2 for d in pos109[0] - pos109[1]]) ** 0.5
ch_len110_0 = sum([(d/unit.angstrom)**2 for d in pos110[0] - pos110[1]]) ** 0.5

sim109.step(1)
sim110.step(1)

pos109 = sim109.context.getState(getPositions=True).getPositions()
pos110 = sim110.context.getState(getPositions=True).getPositions()
ch_len109_1 = sum([(d/unit.angstrom)**2 for d in pos109[0] - pos109[1]]) ** 0.5
ch_len110_1 = sum([(d/unit.angstrom)**2 for d in pos110[0] - pos110[1]]) ** 0.5

print(ch_len109_0 - ch_len109_1)
print(ch_len110_0 - ch_len110_1)
assert abs(ch_len109_0 - ch_len109_1) > 1e-5
assert abs(ch_len110_0 - ch_len110_1) < 1e-5
mattwthompson commented 2 years ago

I don't think it makes sense anymore to put this in the toolkit, but I could yoink this into Interchange if you think it's important?