Determine steps to export to LAMMPS format

[j-wags]

Is your feature request related to a problem? Please describe.
In order to integrate SMIRNOFF-format FFs into OpenKim's benchmarking framework, we will need to specify a simulation engine to run OpenFF calculations. OpenKim is willing to use OpenMM if needed, but would prefer to use LAMMPS, since that's widely used across their stack.

Describe the solution you'd like
A working code example, starting from a molecule file (for example, 1_cyclohexane_1_ethanol.pdb), and producing the input files needed to run a LAMMPS simulation. This will almost certainly use ParmEd or InterMol.

If possible, we should identify what aspects of the System are preserved vs. lost, for example periodic box dimensions, electrostatic and Lennard Jones cutoffs, and nonstandard 1-4 scaling.

Describe alternatives you've considered
If this ends up being unreliable or it is difficult to preserve the relevant information, then OpenKim is OK with us using OpenMM.

Tagging @mattwthompson who has some previous experience with LAMMPS and parameter conversion, and @yafshar from OpenKim.

[mattwthompson]

Quick thoughts:

• Long-term, it's clear that this will be handled by the OpenFF system object

• The short-term strategy should be to find the cleanest path from OpenMM into LAMMPS.

  • Going from OpenMM to ParmEd, as in the linked example, is probably where I'd start. ParmEd does not natively know a thing about LAMMPS but mBuild has a LAMMPS writer that has been used in production for materials research for a couple of years, although I would not trust it 100%. It's named poorly; it doesn't really do any mBuild things - doesn't even take an mBuild object - it's just a writer from parmed.Structure to LAMMPS data files on disk
  • Something similar to above but built internally could provide some benefits, but it would be a lot of re-inventing the wheel
  • InterMol exists in similar space and, although I'm less familiar with it, it seems equipped to do conversions to LAMMPS from AMBER or GROMACS files on disk.

• The scope of what we need to cover here is crucial. LAMMPS supports maybe one or two orders of magnitude more and crazier physics than biophysics codes. If the scope is just what is supported by the OpenFF toolkit right now, and we are comfortable prohibiting things outside of that, that should be workable. Considering things in the SMIRNOFF spec that the toolkit does not yet support would take some investigating.

• Off the top of my head, I can't think of any information that is guaranteed to be lost in this conversion, but some details can't be forced into the data file. LAMMPS, at least the way I and my peers have used it in the past, has needs two files, a data file and an input script (I use data.foo and in.bar, respectively). The data file includes atomic positions and charges, force field parameters, valences terms, and box vectors. Non-bonded cutoffs, electrostatic cutoffs and methods, 1-4 scaling factors, and most anything more complicated are conventionally declared in the input script. In the past I have written out instructions to the user via raising UserWarnings while writing the data file (This data file assumes you have this line in your input script:) and that's obviously not ideal.

[mrshirts]

The shortest term is going through InterMol, since InterMol can already convert to LAMMPS.

[mrshirts]

If the scope is just what is supported by the OpenFF toolkit right now, and we are comfortable prohibiting things outside of that, that should be workable

Mostly, yes. Virtual sites will be supported soon. Other things can be extended when we reach them. We just need to map what OpenFF is doing onto LAMMPS.

[mattwthompson]

Here is a primitive prototype that demonstrates the above idea. It's not elegant but demonstrates that the basic idea is possible without too much hacking. (Note that it's not yet clear if this is generally feasible, but one step in that direction). Notable pain points here are how to keep track of positions before leaving OpenFF/OpenMM world and the strange dropping of a residue name somewhere along the way.

This pipeline should later on be re-imagined to take in a clearly-defined data structure/chemical topology; this just builds up a single molecule from scratch.

Rename the file prototype.ipynb: prototype.txt

[mattwthompson]

This effort has moved to the Interchange project, which provides an export to LAMMPS. It is not yet battle-tested and is under active development, but it does well enough on the bare minimum case presented a few comments above.

In:

from openff.interchange.components.interchange import Interchange
from openff.interchange.drivers import get_lammps_energies, get_openmm_energies
from openmm import app

from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField

ethanol = Molecule.from_smiles("CCO")
cyclohexane = Molecule.from_smiles("C1CCCCC1")

pdbfile = app.PDBFile(
    "examples/using_smirnoff_in_amber_or_gromacs/1_cyclohexane_1_ethanol.pdb"
)
openff_topology = Topology.from_openmm(
    pdbfile.topology, unique_molecules=[ethanol, cyclohexane]
)

sage = ForceField("openff_unconstrained-2.0.0.offxml")

interchange = Interchange.from_smirnoff(force_field=sage, topology=openff_topology)
interchange.positions = pdbfile.positions

print("OpenMM energies:\n")
print(get_openmm_energies(interchange))
print("LAMMPS energies:\n")
print(get_lammps_energies(interchange))
print("Energy comparison:\n")
get_lammps_energies(interchange).compare(get_openmm_energies(interchange))

Out:

OpenMM energies:

Energies:

Bond:          		0.44039154052734375 kJ / mol
Angle:         		155.7012939453125 kJ / mol
Torsion:       		19.563228607177734 kJ / mol
Nonbonded:     		-5.5570980420777545 kJ / mol
vdW:           		None
Electrostatics:		None

LAMMPS energies:

Energies:

Bond:          		0.10526246 kcal / mol
Angle:         		37.213514 kcal / mol
Torsion:       		4.675725 kcal / mol
Nonbonded:     		None
vdW:           		2.3018639467999997 kcal / mol
Electrostatics:		-3.6336973 kcal / mol

Energy comparison:

/Users/mwt/miniconda3/envs/openff-dev/lib/python3.8/site-packages/pandas/core/dtypes/cast.py:1990: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  result[:] = values
Traceback (most recent call last):
  File "run_lammps.py", line 28, in <module>
    get_lammps_energies(interchange).compare(get_openmm_energies(interchange))
  File "/Users/mwt/miniconda3/envs/openff-dev/lib/python3.8/site-packages/openff/interchange/drivers/report.py", line 142, in compare
    raise EnergyError(
openff.interchange.exceptions.EnergyError:
Some energy difference(s) exceed tolerances!
All values are reported in kJ/mol:
      key      diff   tol     ener1     ener2
Nonbonded -0.015293 0.001 -5.572391 -5.557098
Nonbonded -0.015293 0.001 -5.572391 -5.557098

0.015 kJ/mol energy differences do not pass the default tolerances I've assigned in the comparison method, but I'd argue it's a workable error given the feature provided and that the valence terms are within tolerance here. I'll still be working to improve this and other test cases.

Please direct future issues relating to LAMMPS interoperability with OpenFF force fields here: https://github.com/openforcefield/openff-interchange/issues

[j-wags]

Thanks, Matt!

[mrshirts]

There's a variety of reasons the nonbonded energies could be different because of how the bonded interactions - we should discuss a bit. Should I create an issue on the interchange project on resolving this?

[mattwthompson]

Yes, that sounds great! I suspect it is related to openforcefield/openff-interchange#159 but I have yet to isolate it - please feel free to add to that issue or open a new one.