[BUG] ForceFieldInfo sets incorrect cutoff from a TriclinicBox space

[Lactoria-cornuta]

When I load the system with truncated octahedral box using mols = sr.load("react.inpcrd","react.prmtop") and create OpenMM context for futher EMLE simulation with mols, engine = sr.qm.emle(mols, mols[0,1], calculator, "10A", 20) omm = sr.convert.to(mols, "openmm", {"platform":"CPU"}) I get OpenMMException: NonbondedForce: The cutoff distance cannot be greater than half the periodic box size. For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#boxsize

However the OpenMM context with the same cutoff is created normally when I create it using just OpenMM (no EMLE, just normal MD simulation). Also, the context is created normally with EMLE and sire when I use rectangular box: mols = sr.load("react_rect.inpcrd", "react_rect.prmtop")

OS: Ubuntu 20.04 Python: 3.10.14 Sire: 2024.1.0.dev from openbiosim/label/emle

Files to reproduce: react.zip code to reproduce: import sire as sr from emle.calculator import EMLECalculator calculator = EMLECalculator(device="cpu") mols = sr.load("tleap/react.prmtop", "tleap/react.inpcrd") mols, engine = sr.qm.emle(mols, mols[0,1], calculator, "10A", 20) omm = sr.convert.to(mols, "openmm", {"platform":"CPU"})

[lohedges]

What cut-off are you using when doing things via OpenMM directly, i.e. the MM non-bonded cutoff? The cutoff you are setting is just the one for the EMLE force, so is completely independent of the one used by OpenMM directly. If you pas a sensible value via the map, then it works:

omm = sr.convert.to(mols, "openmm", {"platform":"CPU", "cutoff":"10A"})

This works up to around 17A:

omm = sr.convert.to(mols, "openmm", {"platform":"CPU", "cutoff":"17A"})

Is okay, but:

omm = sr.convert.to(mols, "openmm", {"platform":"CPU", "cutoff":"17.5A"})

gives:

ValueError: There was a problem creating the OpenMM context. Perhaps the platform was not supported for this system, options or on this
computer? The error message is: NonbondedForce: The cutoff distance cannot be greater than half the periodic box size.  For more information,
see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#boxsize

I thought the cutoff was set to the default OpenMM value, which is presumably what you are using with your vanilla OpenMM code. I'll check what it's being set to and report back.

Note that converting a system to an OpenMM context via sr.convert will not allow you to perform ML/MM simulations with EMLE. This is only possible via mols.dynamics, passing in the appropriate qm_engine kwarg. If you check the forces in the OpenMM system, then you will see that the EMLE force is missing:

In [1]: import sire as sr

In [2]:
Do you really want to exit ([y]/n)? y
(sire-emle) ~/C/o/s/react ❯❯❯ ipython                                                                                          feature_emle ◼
Python 3.10.13 | packaged by conda-forge | (main, Dec 23 2023, 15:36:39) [GCC 12.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.20.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import sire as sr
   ...: from emle.calculator import EMLECalculator
   ...: calculator = EMLECalculator(device="cpu")
   ...: mols = sr.load("react.prmtop", "react.inpcrd")
   ...: mols, engine = sr.qm.emle(mols, mols[0,1], calculator, "10A", 20)
   ...: omm = sr.convert.to(mols, "openmm", {"platform":"CPU", "cutoff": "10A"})

In [2]: d = mols.dynamics(
   ...:     timestep="1fs",
   ...:     constraint="none",
   ...:     cutoff_type="pme",
   ...:     platform="cpu",
   ...:     qm_engine=engine,
   ...: )

In [3]: d._d._omm_mols.getSystem().getForces()
Out[3]:
[<openmm.openmm.Force; proxy of <Swig Object of type 'OpenMM::Force *' at 0x7df9cc9bf7b0> >,
 <openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7df9cc9bf960> >,
 <openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7df9cc9bfa80> >,
 <openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7dfb5c48e1f0> >,
 <openmm.openmm.PeriodicTorsionForce; proxy of <Swig Object of type 'OpenMM::PeriodicTorsionForce *' at 0x7dfb5c48d560> >,
 <openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7dfb5c48d350> >,
 <openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7dfb5c48daa0> >,
 <openmm.openmm.CustomBondForce; proxy of <Swig Object of type 'OpenMM::CustomBondForce *' at 0x7dfb5c48cea0> >]

In [4]: omm.getSystem().getForces()
Out[4]:
[<openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7dfb5c423d50> >,
 <openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7dfb5c4205a0> >,
 <openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7dfb5c421920> >,
 <openmm.openmm.PeriodicTorsionForce; proxy of <Swig Object of type 'OpenMM::PeriodicTorsionForce *' at 0x7dfb5c422070> >,
 <openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7dfb5c423f90> >,
 <openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7dfb5c508660> >,
 <openmm.openmm.CustomBondForce; proxy of <Swig Object of type 'OpenMM::CustomBondForce *' at 0x7df9cc9dd800> >]

Note that there is one extra force in the first set of forces:

<openmm.openmm.Force; proxy of <Swig Object of type 'OpenMM::Force *' at 0x7df9cc9bf7b0> >

This is the EMLEForce.

[lohedges]

For some reason, the FFInfo object for your molecule has a crazily large cutoff. (Printed at the point that the cutoff is set in the C++ layer.)

const auto cutoff = ffinfo.cutoff().to(SireUnits::nanometers);
qDebug() << ffinfo.toString();
qDebug() << cutoff;
cljff.setCutoffDistance(cutoff);

gives:

"ForceFieldInfo(\n  space=TriclinicBox( ( 41.8721, 0, 0 ), ( -13.9574, 39.4774, 0 ), ( -13.9574, -19.7387, 34.1884 ) ),\n  cutoff_type=CUTOFF,\n  cutoff=1.79769e+308 Å,\n  detail=MM ForceField{ amber::ff,\n               combining_rules = arithmetic,\n               1-4 scaling = 0.833333, 0.5,\n               nonbonded = coulomb, lj,\n               bond = harmonic, angle = harmonic,\n               dihedral = cosine }\n)"
1.79769e+307

Will try to figure out why it's so large.

[lohedges]

Okay, it seems to be a bug setting an ForceFieldInfo object from a system containing a triclinic space:

In [1]: import sire as sr

In [2]: mols = sr.load_test_files("ala.crd", "ala.top")

In [3]: ffinfo = sr.legacy.System.ForceFieldInfo(mols._system)

In [4]: ffinfo
Out[4]:
ForceFieldInfo(
  space=PeriodicBox( ( 26.6843, 28.9808, 24.8784 ) ),
  cutoff_type=CUTOFF,
  cutoff=11.4392 Å,
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)

In [5]: mols = sr.load("react.prmtop", "react.inpcrd")

In [6]: ffinfo = sr.legacy.System.ForceFieldInfo(mols._system)

In [7]: ffinfo
Out[7]:
ForceFieldInfo(
  space=TriclinicBox( ( 41.8721, 0, 0 ), ( -13.9574, 39.4774, 0 ), ( -13.9574, -19.7387, 34.1884 ) ),
  cutoff_type=CUTOFF,
  cutoff=1.79769e+308 Å,
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)

I'll double check a few for a few more systems then update the title accordingly.

Just to note that I am in the process of updating the sire-emle documentation to give more examples and to show how to get the force directly and use with OpenMM in a more traditional way.

[lohedges]

Yes, that's the issue:

In [1]: import BioSimSpace as BSS

In [2]: from sire.legacy.System import ForceFieldInfo

In [3]: system = BSS.IO.readMolecules(["amber/ala/ala.crd", "amber/ala/ala.top"])

In [4]: ForceFieldInfo(system._sire_object)
Out[4]:
ForceFieldInfo(
  space=PeriodicBox( ( 31.3979, 34.1, 29.273 ) ),
  cutoff_type=CUTOFF,
  cutoff=13.6365 Å,
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)

In [5]: box, angles = BSS.Box.truncatedOctahedron(30*BSS.Units.Length.angstrom)

In [6]: system.setBox(box, angles)

In [7]: ForceFieldInfo(system._sire_object)
Out[7]:
ForceFieldInfo(
  space=TriclinicBox( ( 30, 0, 0 ), ( -10, 28.2843, 0 ), ( -10, -14.1421, 24.4949 ) ),
  cutoff_type=CUTOFF,
  cutoff=1.79769e+308 Å,
  detail=MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }
)

[lohedges]

This is an easy fix. It's because there is no TriclinicBox::maximumCutoff method.

[Lactoria-cornuta]

Thanks for a fast reply! I didn't pay attention that the EMLE force dissappears after sr.convert.to(). Thanks again for your clarification! Off topic: how long will it take you to update the sire-emle docs? I'm looking forward to using EMLE in my research.

[lohedges]

I'm hoping to get this done by the end of the week.

[Lactoria-cornuta]

Great news! Using EMLE in OpenMM is an important feature for running metadynamics simulations as OpenMM is integrated with PLUMED (widely used metadynamics plugin for MD simulations) while sire isn't.