openforcefield
Make HMR simulations easier
What's the status/priority level on this? We wanted to try looking at the difference in membranes with and without HMR.
We are happy to test it out. What issues do you potentially see?
No interest yet, though it might be useful as we better integrate with OpenFE. (One of their settings objects has HMR turned on by default, I think with 3 amu mass transfer.)
A few issues that immediately come to mind, prior to having written any code:
- I have no idea how non-OpenMM codes implement this. Are you looking to run in OpenMM or another engine?
- The toolkit doesn't allow atomic masses to diverge from what the element says the mass should be, so I'd first try to implement HMR as late as possible in the pipeline, i.e. as a final step in
to_openmm_systembut not before. So the.topologyobject would always have 1 amu masses, but the relevant OpenMM objects wouldn't. Fortunately this seems to be pretty much what OpenMM already does.
In fact, it looks like HMR could be applied totally as a post-processing step? This script shifts hydrogen masses around on the OpenMM objects only. I didn't test it past looking at the masses in the system (it doesn't seem to me like the topology needs its masses changed). Could you & your students adapt this into a workflow and see if it produces the sort of results you'd expect?
Not here I'm skipping waters for simplicity (in one-off uses you could just skip all waters, for something hitting production we'd want to check if it's water and also if it's fully constrained); the OpenMM code I linked above checks if the water is rigid, which would also be done in a final implementation here.
import openmm
import openmm.app
from openff.toolkit import ForceField, Molecule
def apply_hmr(
system: openmm.System,
topology: openmm.app.Topology,
hydrogen_mass: float,
):
"""
Apply HMR via transferring mass from heavy atoms to bonded hydrogens.
Follows OpenMM's own implementation in ForceField.createSystem:
https://github.com/openmm/openmm/blob/71f4b3fcb0b07d0fe41302df03f5373c8424cefa/wrappers/python/openmm/app/forcefield.py#L1257-L1288
"""
import openmm.unit
from openmm.app.element import hydrogen
hydrogen_mass *= openmm.unit.dalton
for atom1, atom2 in topology.bonds():
if atom1.element is hydrogen:
(atom1, atom2) = (atom2, atom1)
if atom2.element is hydrogen and atom1.element not in (hydrogen, None):
transfer_mass = hydrogen_mass - system.getParticleMass(atom2.index)
system.setParticleMass(
atom2.index,
hydrogen_mass,
)
system.setParticleMass(
atom1.index,
system.getParticleMass(atom1.index) - transfer_mass,
)
molecule = Molecule.from_smiles("CCO")
molecule.generate_conformers(n_conformers=1)
interchange = ForceField("openff-2.1.0.offxml").create_interchange(
molecule.to_topology(),
)
system = interchange.to_openmm_system()
topology = interchange.to_openmm_topology()
hydrogen_mass = 3.0
apply_hmr(system, topology, hydrogen_mass)
for particle_index, atom in enumerate(
topology.atoms(),
):
print(
atom.element.symbol,
system.getParticleMass(particle_index),
)
"""
$ python quick_hmr.py
C 6.034621 Da
C 8.026674 Da
O 14.007377 Da
H 3.0 Da
H 3.0 Da
H 3.0 Da
H 3.0 Da
H 3.0 Da
H 3.0 Da
"""
It looks like @chapincavender is using a similar solution in his work (except not skipping waters out of laziness like I did): https://github.com/openforcefield/proteinbenchmark/blob/70e55c8311c0b862cc6fde3634d3f49f85d5bb26/proteinbenchmark/system_setup.py#L626-L648
I have no idea how non-OpenMM codes implement this. Are you looking to run in OpenMM or another engine?
GROMACS, for various reasons. You just rewrite the masses in the .top files to get it to work. So one really doesn't have a programmatic way of doing it other than writing python scripts to rewrite the .top files - there are no API's to handle that.
So, ideally it would be nice to do this upstream, so that one doesn't have to have a separate script for each output format. All MD simulations that I am aware of allow one to rewrite the masses essentially arbitrarily - they don't enforce any type of particle having a constrained atomic mass.
The toolkit doesn't allow atomic masses to diverge from what the element says the mass should be,
Interesting. One could imagine that one would want to do this like, have different isotopes. Is there a strict reason for this? To catch errors later? One could imagine the masses being set to atomic masses by default on creation, but allowing them to be overwritten.
It looks like @chapincavender is using a similar solution
I guess it looks like what is happening there is chapin creates an OpenMM system from interchange, does the HMR, and then uses interchange to convert from OpenMM to GROMACS. Which is a little indirect, but I guess it works!
For the GROMACS engine in the proteinbenchmark repo, @ajfriedman22 developed a workflow that exports to OpenMM and then uses the parmed python API to apply HMR and export to GROMACS. My understanding is that the limitation here is the fixed mass expectation that @mattwthompson pointed out. The OpenFF toolkit / Interchange cannot read the OpenMM system after HMR is applied, so we can't use Interchange to do the OpenMM -> GROMACS step.
Agree with above; briefly my impression is that getting it in GROMACS would involve some sort of "thing done at the last step" similar to how it's done in OpenMM, but not directly using the OpenMM object.
@chapincavender thanks for letting me know what is going on in my group, lol :)
We'll figure out something shorter term (probably editing the files directly in GROMACS).
GIVEN that HRM is a "thing the cool kids do" now, it SEEMS like it would be nice to do this in Interchange eventually, and not have a bunch of different scripts post-conversion, given that all of the MD simulation codes we have been considering exporting to are cool with alternate masses. Especially if one has a simulation in OpenMM that happens to be using HMR, and wants to ingest it to interchange to do something to it, it doesn't seems like "No I don't like your masses" is a useful blocking error?
Thoughts?
The main branch now has partial support for HMR in both OpenMM and GROMACS. If anybody has time to take a test those in their workflows that would be much appreciated, otherwise they'll land in the next release by sometime early next week.
The general pattern is to add a hydrogen_mass=3.0 argument to the main export method. Following the convention in OpenMM, this argument specifies how much mass the hydrogens should end up having, not necessarily how much should be transferred away. There are some nuances that happen under the hood, but so far I haven't come across anything surprising (hence my asking here).
0.3.20 is released and includes the above changes
I have not looked into how Amber does HMR and I'm unlikely to look into how it could be done in LAMMPS until somebody shows interested in it (which isn't to say it would necessarily be hard to do.
I'm comfortable assessing that a good stab at this is complete - IIUC, people have been using Interchange to enable HMR simulations in OpenMM for a few months now, possibly also with GROMACS.
As described above, doing the same thing in Amber or LAMMPS has not been looked into - if these are things people want, please get in touch i.e. open a new issue.