Simulation with 2 `espaloma` molecules explodes

[JSLJ23]

I am trying to run OpenMM simulations with espaloma for the SAMPL8 Cyclodextrin host guest complexes, specifically TEMOA + G1. I am trying to use espaloma as it is fast enough to parameterise the large host molecule, whereas other methods like GAFF and SMIRNOFF pretty much take forever. Adding the host and guest molecules into the same system (docked pose) somehow causes the simulation to exploded even after running energy minimisation. Any idea why this might be the case?

import sys
from pathlib import Path

from openff.toolkit.topology import Molecule
from openmm import LangevinMiddleIntegrator
from openmm.app import Modeller, Simulation, StateDataReporter
from openmm.app.forcefield import ForceField
from openmm.openmm import Platform
from openmm.unit import femtosecond, kelvin, kilojoule_per_mole, nanometers, picosecond, molar
from openmmforcefields.generators import EspalomaTemplateGenerator

host_sdf_path = Path("./TEMOA.sdf")
guest_sdf_path = Path("./TEMOA_G1_a_docked.sdf")

system_name = "TEMOA_G1_a"

host_off_mol = Molecule.from_file(host_sdf_path)
host_mol_topology = host_off_mol.to_topology().to_openmm()
host_mol_positions = host_off_mol.to_topology().get_positions().to_openmm()

guest_off_mol = Molecule.from_file(guest_sdf_path)
guest_mol_topology = guest_off_mol.to_topology().to_openmm()
guest_mol_positions = guest_off_mol.to_topology().get_positions().to_openmm()

modeller = Modeller(host_mol_topology, host_mol_positions)
modeller.add(guest_mol_topology, guest_mol_positions)

force_field = ForceField("amber14-all.xml", "amber14/tip3p.xml")
template_generator = EspalomaTemplateGenerator(molecules=[host_off_mol, guest_off_mol], forcefield="espaloma-0.3.2")
force_field.registerTemplateGenerator(template_generator.generator)

modeller.addSolvent(force_field, padding=1.0 * nanometers, ionicStrength=0.15 * molar)

system = force_field.createSystem(modeller.topology)

total_time_steps = 1_000_000
n_fs = 1
time_in_ns = total_time_steps * n_fs / 1_000_000
temperature_k = 298.15

temperature = temperature_k * kelvin
friction_coeff = 1 / picosecond
time_step = n_fs * femtosecond
integrator = LangevinMiddleIntegrator(temperature, friction_coeff, time_step)
integrator.setRandomNumberSeed(42)

platform = Platform.getPlatformByName("CUDA")
properties = {"DeviceIndex": "0"}
simulation = Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)
simulation.context.setVelocitiesToTemperature(temperature)

state = simulation.context.getState(getEnergy=True, enforcePeriodicBox=False)
energy = state.getPotentialEnergy().value_in_unit(kilojoule_per_mole)
print(f"Initial energy: {energy} kj/mol")

reporter_console = StateDataReporter(
    file=sys.stdout,
    reportInterval=100_000,
    step=True,
    time=True,
    potentialEnergy=True,
    temperature=True,
)

simulation.reporters.append(reporter_console)

simulation.minimizeEnergy()
simulation.step(total_time_steps)

Initial energy: -10295.931816101074 kj/mol

ValueError                                Traceback (most recent call last)
Cell In[11], line 4
      1 simulation.minimizeEnergy()
      3 start = perf_counter()
----> 4 simulation.step(total_time_steps)
      5 end = perf_counter()
      6 print(f"Time taken for MD: {end - start} seconds.")

File ~/mambaforge/envs/ab_initio_gradients_py310/lib/python3.10/site-packages/openmm/app/simulation.py:147, in Simulation.step(self, steps)
    145 def step(self, steps):
    146     """Advance the simulation by integrating a specified number of time steps."""
--> 147     self._simulate(endStep=self.currentStep+steps)

File ~/mambaforge/envs/ab_initio_gradients_py310/lib/python3.10/site-packages/openmm/app/simulation.py:249, in Simulation._simulate(self, endStep, endTime)
    247     self._generate_reports(wrapped, True)
    248 if len(unwrapped) > 0:
--> 249     self._generate_reports(unwrapped, False)

File ~/mambaforge/envs/ab_initio_gradients_py310/lib/python3.10/site-packages/openmm/app/simulation.py:269, in Simulation._generate_reports(self, reports, periodic)
    265 state = self.context.getState(getPositions=getPositions, getVelocities=getVelocities, getForces=getForces,
    266                               getEnergy=getEnergy, getParameters=True, enforcePeriodicBox=periodic,
    267                               groups=self.context.getIntegrator().getIntegrationForceGroups())
    268 for reporter, next in reports:
--> 269     reporter.report(self, state)

File ~/mambaforge/envs/ab_initio_gradients_py310/lib/python3.10/site-packages/openmm/app/dcdreporter.py:107, in DCDReporter.report(self, simulation, state)
    102 if self._dcd is None:
    103     self._dcd = DCDFile(
    104         self._out, simulation.topology, simulation.integrator.getStepSize(),
    105         simulation.currentStep, self._reportInterval, self._append
    106     )
--> 107 self._dcd.writeModel(state.getPositions(), periodicBoxVectors=state.getPeriodicBoxVectors())

File ~/mambaforge/envs/ab_initio_gradients_py310/lib/python3.10/site-packages/openmm/app/dcdfile.py:125, in DCDFile.writeModel(self, positions, unitCellDimensions, periodicBoxVectors)
    123     positions = positions.value_in_unit(nanometers)
    124 if any(math.isnan(norm(pos)) for pos in positions):
--> 125     raise ValueError('Particle position is NaN.  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan')
    126 if any(math.isinf(norm(pos)) for pos in positions):
    127     raise ValueError('Particle position is infinite.  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan')

ValueError: Particle position is NaN.  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan

TEMOA_G1.tar.gz

[JSLJ23]

I've tried this in the presence of solvent and also in a vacuum and they both result in particle positions going out of the physical range. However, when it is just the host molecule TEMOA alone, the simulation runs perfectly fine in explicit solvent and in a vacuum.

[ijpulidos]

@JSLJ23 Thanks for showing interest in our tools. I noticed that if you just minimize before equilibrating, the system behaves as expected. As far as I can tell minimizing before equilibrating is the standard practice. The code is just as follows, I just moved the call to minimizeEnergy to be done before setVelocitiesToTemperature. I hope this helps.

import sys
from pathlib import Path

from openff.toolkit.topology import Molecule
from openmm import LangevinMiddleIntegrator
from openmm.app import Modeller, Simulation, StateDataReporter, PDBReporter
from openmm.app.forcefield import ForceField
from openmm.openmm import Platform
from openmm.unit import femtosecond, kelvin, kilojoule_per_mole, nanometers, picosecond, molar
from openmmforcefields.generators import EspalomaTemplateGenerator

host_sdf_path = Path("./TEMOA.sdf")
guest_sdf_path = Path("./TEMOA_G1_a_docked.sdf")

system_name = "TEMOA_G1_a"

host_off_mol = Molecule.from_file(host_sdf_path)
host_mol_topology = host_off_mol.to_topology().to_openmm()
host_mol_positions = host_off_mol.to_topology().get_positions().to_openmm()

guest_off_mol = Molecule.from_file(guest_sdf_path)
guest_mol_topology = guest_off_mol.to_topology().to_openmm()
guest_mol_positions = guest_off_mol.to_topology().get_positions().to_openmm()

modeller = Modeller(host_mol_topology, host_mol_positions)
modeller.add(guest_mol_topology, guest_mol_positions)

force_field = ForceField("amber14-all.xml", "amber14/tip3p.xml")
template_generator = EspalomaTemplateGenerator(molecules=[host_off_mol, guest_off_mol], forcefield="espaloma-0.3.2")
force_field.registerTemplateGenerator(template_generator.generator)

modeller.addSolvent(force_field, padding=1.0 * nanometers, ionicStrength=0.15 * molar)

system = force_field.createSystem(modeller.topology)

total_time_steps = 1_000_000
n_fs = 1
time_in_ns = total_time_steps * n_fs / 1_000_000
temperature_k = 298.15

temperature = temperature_k * kelvin
friction_coeff = 1 / picosecond
time_step = n_fs * femtosecond
integrator = LangevinMiddleIntegrator(temperature, friction_coeff, time_step)
integrator.setRandomNumberSeed(42)

platform = Platform.getPlatformByName("CUDA")
properties = {"DeviceIndex": "0"}
simulation = Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)
simulation.minimizeEnergy()  # Minimize first
simulation.context.setVelocitiesToTemperature(temperature)

state = simulation.context.getState(getEnergy=True, enforcePeriodicBox=False)
energy = state.getPotentialEnergy().value_in_unit(kilojoule_per_mole)
print(f"Initial energy: {energy} kj/mol")

reporter_console = StateDataReporter(
    file=sys.stdout,
    reportInterval=100_000,
    step=True,
    time=True,
    potentialEnergy=True,
    temperature=True,
)

simulation.reporters.append(reporter_console)

simulation.reporters.append(PDBReporter('output.pdb', 1000))

simulation.step(total_time_steps)

[JSLJ23]

Thank you so much for your help, yes this works!