ParmEd / ParmEd

Parameter/topology editor and molecular simulator
https://parmed.github.io/ParmEd
390 stars 148 forks source link

OpenMM -> ParmEd, dealing with hydrogen bond constraints #1068

Open ahy3nz opened 4 years ago

ahy3nz commented 4 years ago

Hi everyone,

I've been messing around with smirnoff99frosst-1.1.0.offxml. It's making OpenMM systems great, but I noticed there's some issue with the conversion to parmed. Because bonds-containing-hydrogen are constrained, the openmm.HarmonicBondForce does not contain parameters for hydrogen bonds. During the openmm->parmed conversion, many parmed.bond objects will not have an associated parmed.bondtype because of the lack of parameters in the HarmonicBondForce.

When trying to write the resultant parmed.Structure to gromacs topology files for gromacs simulation, some gromacs topology bond lines have no bond parameters (because the parmed.Structure didn't have any). A gromacs error is then raised because no default bond types are detected for the gromacs bonds that have no bond parameters.

What is the proper way to deal with this? Modify the gromacs files prior to simulation? Add "empty" bond types to the parmed.Structure after conversion from openmm with hydrogen bond constraints?

swails commented 4 years ago

I would personally create a new bond type that you just assign to all of the constrained bonds.

Of course there's the understanding that if you don't constrain those bonds in GROMACS, you'll get garbage.

Another option, by the way, is to create the system with the flexibleConstraints keyword if that's an option. Many OpenMM createSystem calls take a flexibleConstraints keyword that will assign parameters to the constrained DoF.

ahy3nz commented 4 years ago

https://github.com/openforcefield/openforcefield/blob/master/openforcefield/typing/engines/smirnoff/forcefield.py#L1131

openforcefield.forcefield.create_openmm_system does not yet handle flexibleConstraints as a kwarg, so the only choice is to create the openmmsystem with the hydrogen bond constraints. Upon conversion to ParmEd, then a user would have to "post-process" to specify bond parameters for the hydrogen bonds, with the understanding that they should be constrained.

Is this something worth submitting a PR, or something best left to each user to handle?