Support for OpenMM

[InnocentBug]

I am considering to use Foyer as tool for a tool to generate OPLSAA forcefields from smiles string definitions of molecules and run them OpenMM. Currently, I can only find examples that write in the GROMACS format. Are there examples, that I am missing that describe different file formats? I am most interested in OpenMM, and I thought that since you are using there file format to describe forcefields, that should be an easy step. Or am I mistaken and Foyer does not support this?

[justinGilmer]

Hey @InnocentBug , thank you for the issues raised today and welcome! At the moment, this should be totally possible with our workflow, using ParmEd as the backend data structure (currently the default option).

Here is a link to an example in the ParmEd docs that loads in a GRO and TOP file as the starting point. Also, it should be possible to take the example there, and once you have parametrized a system with foyer to go straight to an openMM system. See the snippet below.

# OpenMM Imports
import sys
import simtk.openmm as mm
import simtk.openmm.app as app

# ParmEd Imports
from parmed import load_file
from parmed.openmm.reporters import NetCDFReporter
from parmed import unit as u

# mosdef imports
import mbuild as mb
import foyer

eth = mb.load("CC", smiles=True)
eth_box = mb.packing.fill_box(eth, n_compounds=100, box=mb.Box(lengths=[1,1,1]))

oplsaa = foyer.forcefields.load_OPLSAA()

param_struct = oplsaa.apply(eth_box)

# Create the OpenMM system
print('Creating OpenMM System')
system = param_struct.createSystem(nonbondedMethod=app.NoCutoff,
                          constraints=app.HBonds, implicitSolvent=app.GBn2,
                          implicitSolventSaltConc=0.1*u.moles/u.liter,
)

# Create the integrator to do Langevin dynamics
integrator = mm.LangevinIntegrator(
                        300*u.kelvin,       # Temperature of heat bath
                        1.0/u.picoseconds,  # Friction coefficient
                        2.0*u.femtoseconds, # Time step
)

# Define the platform to use; CUDA, OpenCL, CPU, or Reference. Or do not specify
# the platform to use the default (fastest) platform
platform = mm.Platform.getPlatformByName('CPU')

# Create the Simulation object
sim = app.Simulation(param_struct.topology, system, integrator, platform,)

# Set the particle positions
sim.context.setPositions(param_struct.positions)

# Minimize the energy
print('Minimizing energy')
sim.minimizeEnergy(maxIterations=500)

# Set up the reporters to report energies and coordinates every 100 steps
sim.reporters.append(
        app.StateDataReporter(sys.stdout, 100, step=True, potentialEnergy=True,
                              kineticEnergy=True, temperature=True)
)
sim.reporters.append(
        NetCDFReporter('dhfr_gb.nc', 100, crds=True)
)

# Run dynamics
print('Running dynamics')
sim.step(10000)

[justinGilmer]

The major changes I added were: building system using mbuild, parametrizing using foyer's oplsaa xml, using the parametrized parmed structure from foyer to then set the rest.

[InnocentBug]

Thanks for the explanation. I think I know where the confusion cam from. I wanted to write files first that describe the topology and forcefields first. (Preferably in the native XML format of OpenMM.) And have those files on record for reproducibility. But your approach should work too. I can always write a state with OpenMM and read that one.