Closed j-wags closed 3 years ago
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.
parmed.Structure
to LAMMPS data files on diskThe 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 UserWarning
s while writing the data file (This data file assumes you have this line in your input script:
) and that's obviously not ideal.
The shortest term is going through InterMol, since InterMol can already convert to LAMMPS.
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.
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
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
Thanks, Matt!
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?
Yes, that sounds great! I suspect it is related to https://github.com/openforcefield/openff-interchange/issues/159 but I have yet to isolate it - please feel free to add to that issue or open a new one.
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.