Parlsey 1.2.0 fails on HMR simulations with propyne-containing molecules

[j-wags]

Progress

• Make reproducing case for reported issue
• Ensure case reproduces for openff-1.2.0, but not openff-1.0.0 or openff-1.1.0
• Make openff-1.2.1 candidate (both constrained and unconstrained), with changes described below
• Ensure candidate FF does not reproduce error
• Follow Standard Release Process

Problem Description

On Aug 11, 2020:

Christopher Bayly 11:09
We have just uncovered a really bad symptom with Parsley 1.2 when running MD on a molecule containing a a propyne substituent on a phenyl, as in c1ccccc1C#CC. The symptom is that the MD blows up and returns a NaN, but it only happens when Hydrogen Mass Repartitionlng (HMR) is turned on. We looked into it and the 4 fs timestep with HMR is the cause: the methyl bond with the triple-bond carbon hits a resonant frequency and ultimately explodes causing the NaN. This should not happen and I think I can see why: I think the force constants for the C#C triple bond and the attached methyl C-CH3 do not make sense. For the latter (single bond), bond type b24 has a super-high force constant 1.2977e+03 that you would expect for a triple bond, not a single bond. Conversely the triple bond with bond type b27 has a weak force constant 8.23e+02 that you would expect for a single bond not a triple bond. You can check this against the other "single bond" and "double/triple bond" force constants in the same section of the offxml. The smoking gun is that when we use GAFF2 which has more-or-less expected k's (high for the triple bond, regular low for the single), there is no issue with MD using HMR. (edited)

(gif from Gaetano Calabro)

Proposed resolution

Per decision in FF-release meeting on September 3, 2020

Decisions:

We will modify 1.2.0 by hand to assign b24’s k to 800, b27’s k to 1000 → Rename to b24-MANand b27-MAN → test against propyne containing molecule (jax 15 beta-secretase case from DH) → Release 1.2.1

[jchodera]

Do we know if this happens in both OpenMM and gromacs, or just OpenMM? Also, which integrator was used?

[j-wags]

OpenMM presumably. Unsure about integrator. @cbayly13 indicated that the error happens reliably and quickly, but is unable to share code since it's internal. I'm working on a reproducing case and will upload that once I have it.

[j-wags]

@jchodera Do you have a code snippet using OpenMM+HMR? DIgging through docs now but having an example would save me some time.

[jchodera]

Sure!

# Create an openforcefield Molecule object from SMILES                                                                                                                                                                                         
from openforcefield.topology import Molecule
molecule = Molecule.from_smiles('c1ccccc1C#CC')
# Define the keyword arguments to feed to ForceField; use heavy hydrogens and constrained X-H bonds                                                                                                                                            
# (Note that this differs a bit from allowing SMIRNOFF to control masses and cosntraints)                                                                                                                                                      
from simtk import openmm, unit
from simtk.openmm import app
forcefield_kwargs = { 'constraints' : app.HBonds, 'rigidWater' : True, 'removeCMMotion' : False, 'hydrogenMass' : 4*unit.amu }
# Initialize a SystemGenerator using GAFF                                                                                                                                                                                                      
from openmmforcefields.generators import SystemGenerator
system_generator = SystemGenerator(small_molecule_forcefield='openff-1.2.0', forcefield_kwargs=forcefield_kwargs, molecules=molecule)
# Create an OpenMM System from an Open Force Field toolkit Topology object                                                                                                                                                                     
system = system_generator.create_system(molecule.to_topology().to_openmm())
# Run a simulation                                                                                                                                                                                                                             
temperature = 300 * unit.kelvin
collision_rate = 1.0 / unit.picoseconds
timestep = 4.0 * unit.femtoseconds
# Use BAOAB                                                                                                                                                                                                                                    
integrator = openmm.LangevinMiddleIntegrator(temperature, collision_rate, timestep)
context = openmm.Context(system, integrator)
molecule.generate_conformers()
context.setPositions(molecule.conformers[0])
# Integrate                                                                                                                                                                                                                                    
niterations = 1000
nsteps_per_iteration = 250
for iteration in range(niterations):
    integrator.step(nsteps_per_iteration)
    state = context.getState(getEnergy=True)
    print(f'{state.getTime()/unit.picoseconds:12.3f} ps :  potential {state.getPotentialEnergy()/unit.kilocalories_per_mole:16.3f} kcal/mol')

output with openff-1.2.0:

       1.000 ps :  potential           39.453 kcal/mol
       2.000 ps :  potential           43.479 kcal/mol
       3.000 ps :  potential           48.913 kcal/mol
       4.000 ps :  potential           47.405 kcal/mol
       5.000 ps :  potential           43.038 kcal/mol
       6.000 ps :  potential           53.795 kcal/mol
       7.000 ps :  potential              nan kcal/mol
       8.000 ps :  potential              nan kcal/mol
       9.000 ps :  potential              nan kcal/mol
      10.000 ps :  potential              nan kcal/mol
      11.000 ps :  potential              nan kcal/mol

with openff-1.0.0:

       1.000 ps :  potential           73.354 kcal/mol
       2.000 ps :  potential           71.585 kcal/mol
       3.000 ps :  potential           74.965 kcal/mol
       4.000 ps :  potential           74.602 kcal/mol
       5.000 ps :  potential           78.396 kcal/mol
       6.000 ps :  potential           77.423 kcal/mol
       7.000 ps :  potential           79.104 kcal/mol
       8.000 ps :  potential           75.350 kcal/mol
       9.000 ps :  potential           75.802 kcal/mol
      10.000 ps :  potential           77.676 kcal/mol
      11.000 ps :  potential           76.635 kcal/mol
      12.000 ps :  potential           74.958 kcal/mol
      13.000 ps :  potential           74.385 kcal/mol
      14.000 ps :  potential           77.336 kcal/mol
      15.000 ps :  potential           76.727 kcal/mol
      16.000 ps :  potential           75.620 kcal/mol
      17.000 ps :  potential           71.890 kcal/mol
      18.000 ps :  potential           76.717 kcal/mol
      19.000 ps :  potential           78.543 kcal/mol
      20.000 ps :  potential           75.369 kcal/mol
      21.000 ps :  potential           78.076 kcal/mol
      22.000 ps :  potential           76.400 kcal/mol
      23.000 ps :  potential           74.812 kcal/mol
      24.000 ps :  potential           76.822 kcal/mol
      25.000 ps :  potential           77.451 kcal/mol
      26.000 ps :  potential           76.050 kcal/mol
      27.000 ps :  potential           74.618 kcal/mol
      28.000 ps :  potential           74.391 kcal/mol
      29.000 ps :  potential           79.433 kcal/mol
      30.000 ps :  potential           74.723 kcal/mol
      31.000 ps :  potential           83.331 kcal/mol
      32.000 ps :  potential           75.466 kcal/mol
      33.000 ps :  potential           75.762 kcal/mol
      34.000 ps :  potential           74.493 kcal/mol
      35.000 ps :  potential           74.044 kcal/mol
      36.000 ps :  potential           77.265 kcal/mol
...

[jchodera]

That said, I'm not certain this isn't an OpenMM specific numerical instability with angles near pi, or torsions when any one angle nears pi. I'm pretty sure we are still not using the numerically stable expressions derived by Bill Swope for these terms, which may be needed.

[j-wags]

Thanks @jchodera!

Here's the code to make both normal and unconstrained openff-1.2.1, reproduce the error, and show that it affects the expected force fields.

Making the new FF files

from openforcefield.typing.engines.smirnoff import ForceField
ff = ForceField('openff-1.2.0.offxml')
bh = ff.get_parameter_handler('Bonds')
for param in bh.parameters:
    if param.id == "b24":
        param.k = 800. * unit.kilocalorie_per_mole / (unit.angstrom ** 2)
        print(param)
    if param.id == "b27":
        param.k = 1000. * unit.kilocalorie_per_mole / (unit.angstrom ** 2)
        print(param)
ff.to_file('openff-1.2.1.offxml')

<BondType with smirks: [#6X2:1]-[#6X4:2] length: 1.465575905254 A k: 800.0 kcal/(A**2 mol) id: b24 >
<BondType with smirks: [#6X2:1]#[#6X2:2] length: 1.218049770442 A k: 1000.0 kcal/(A**2 mol) id: b27 >

from openforcefield.typing.engines.smirnoff import ForceField
ff = ForceField('openff_unconstrained-1.2.0.offxml')
bh = ff.get_parameter_handler('Bonds')
for param in bh.parameters:
    if param.id == "b24":
        param.k = 800. * unit.kilocalorie_per_mole / (unit.angstrom ** 2)
        print(param)
    if param.id == "b27":
        param.k = 1000. * unit.kilocalorie_per_mole / (unit.angstrom ** 2)
        print(param)
ff.to_file('openff_unconstrained-1.2.1.offxml')

<BondType with smirks: [#6X2:1]-[#6X4:2] length: 1.465575905254 A k: 800.0 kcal/(A**2 mol) id: b24 >
<BondType with smirks: [#6X2:1]#[#6X2:2] length: 1.218049770442 A k: 1000.0 kcal/(A**2 mol) id: b27 >

Reproducing the error for 1.2.0, and showing it doesn't affect 1.0.0, 1.1.0, and 1.2.1

# Create an openforcefield Molecule object from SMILES                                                                                                                                                                                         
from openforcefield.topology import Molecule
molecule = Molecule.from_smiles('c1ccccc1C#CC')
# Define the keyword arguments to feed to ForceField; use heavy hydrogens and constrained X-H bonds                                                                                                                                            
# (Note that this differs a bit from allowing SMIRNOFF to control masses and cosntraints)                                                                                                                                                      
from simtk import openmm, unit
from simtk.openmm import app
forcefield_kwargs = { 'constraints' : app.HBonds, 'rigidWater' : True, 'removeCMMotion' : False, 'hydrogenMass' : 4*unit.amu }
# Initialize a SystemGenerator using GAFF                                                                                                                                                                                                      
from openmmforcefields.generators import SystemGenerator
for ff_name in ('openff-1.2.0', 'openff-1.1.0', 'openff-1.0.0', 'openff-1.2.1'):
    print()
    print()
    print(ff_name)
    system_generator = SystemGenerator(small_molecule_forcefield=ff_name, forcefield_kwargs=forcefield_kwargs, molecules=molecule)
    # Create an OpenMM System from an Open Force Field toolkit Topology object                                                                                                                                                                     
    system = system_generator.create_system(molecule.to_topology().to_openmm())
    # Run a simulation                                                                                                                                                                                                                             
    temperature = 300 * unit.kelvin
    collision_rate = 1.0 / unit.picoseconds
    timestep = 4.0 * unit.femtoseconds
    # Use BAOAB                                                                                                                                                                                                                                    
    integrator = openmm.LangevinIntegrator(temperature, collision_rate, timestep)
    context = openmm.Context(system, integrator)
    molecule.generate_conformers()
    context.setPositions(molecule.conformers[0])
    # Integrate                                                                                                                                                                                                                                    
    niterations = 10
    nsteps_per_iteration = 250
    for iteration in range(niterations):
        integrator.step(nsteps_per_iteration)
        state = context.getState(getEnergy=True)
        print(f'{state.getTime()/unit.picoseconds:12.3f} ps :  potential {state.getPotentialEnergy()/unit.kilocalories_per_mole:16.3f} kcal/mol')

openff-1.2.0
1.000 ps : potential 43.568 kcal/mol
2.000 ps : potential nan kcal/mol
3.000 ps : potential nan kcal/mol
4.000 ps : potential nan kcal/mol
5.000 ps : potential nan kcal/mol
6.000 ps : potential nan kcal/mol
7.000 ps : potential nan kcal/mol
8.000 ps : potential nan kcal/mol
9.000 ps : potential nan kcal/mol
10.000 ps : potential nan kcal/mol

openff-1.1.0
1.000 ps : potential 72.692 kcal/mol
2.000 ps : potential 74.271 kcal/mol
3.000 ps : potential 70.596 kcal/mol
4.000 ps : potential 73.193 kcal/mol
5.000 ps : potential 78.326 kcal/mol
6.000 ps : potential 81.247 kcal/mol
7.000 ps : potential 79.374 kcal/mol
8.000 ps : potential 77.016 kcal/mol
9.000 ps : potential 70.389 kcal/mol
10.000 ps : potential 72.727 kcal/mol

openff-1.0.0
1.000 ps : potential 75.120 kcal/mol
2.000 ps : potential 84.182 kcal/mol
3.000 ps : potential 78.065 kcal/mol
4.000 ps : potential 76.811 kcal/mol
5.000 ps : potential 81.001 kcal/mol
6.000 ps : potential 81.146 kcal/mol
7.000 ps : potential 77.968 kcal/mol
8.000 ps : potential 78.745 kcal/mol
9.000 ps : potential 78.205 kcal/mol
10.000 ps : potential 76.285 kcal/mol

openff-1.2.1
1.000 ps : potential 39.100 kcal/mol
2.000 ps : potential 47.156 kcal/mol
3.000 ps : potential 42.099 kcal/mol
4.000 ps : potential 43.874 kcal/mol
5.000 ps : potential 45.213 kcal/mol
6.000 ps : potential 40.509 kcal/mol
7.000 ps : potential 46.466 kcal/mol
8.000 ps : potential 48.286 kcal/mol
9.000 ps : potential 51.831 kcal/mol
10.000 ps : potential 48.733 kcal/mol

[j-wags]

We are deviating from the meeting notes above: They say to rename the modified parameters to b24-MAN, but @trevorgokey points out (and I agree) that this is confusing, and even b24-manual is inconvenient, because it makes this force field incomparable (by parameter labeling studies) to others in the Parsley line.

We will leave a manual comment in the OFFXML files explaining the change.

[j-wags]

Ah, testing this with the UNCONSTRAINED candidate FF still shows the explosion. Is this acceptable, or should we look at other solutions?

openff_unconstrained-1.2.1
1.000 ps : potential 42.302 kcal/mol
2.000 ps : potential 44.210 kcal/mol
3.000 ps : potential 41.396 kcal/mol
4.000 ps : potential 41.892 kcal/mol
5.000 ps : potential 50.392 kcal/mol
6.000 ps : potential 46.798 kcal/mol
7.000 ps : potential 50.715 kcal/mol
8.000 ps : potential 45.248 kcal/mol
9.000 ps : potential nan kcal/mol
10.000 ps : potential nan kcal/mol

[j-wags]

Verified with @davidlmobley -- This is expected, and only the constrained version of the new FF should be used with HMR.

Jeffrey Wagner  15:38
@John Chodera @dmobley @cbayly I've made the proposed changes to create a 1.2.1 release candidate. The normal version of this FF openff-1.2.1 successfully simulates the propyne system. However, openff_unconstrained-1.2.1 still sometimes explodes.
I'm unfamiliar with HMR. Is HMR always run with hbond constraints? If not, the changes in the new unconstrained FF are not adequate to fix the propyne issue.

David Mobley  15:39
That should be what we’re after! The unconstrained FF is not for MD.
15:40
When running MD without constraints on HBonds, one needs a smaller timestep. HMR does nothing to fix this. It simply makes it so that you can get by with an EVEN LARGER timestep than normal when running with HBond constraints.