[BUG] Incorrect application of friction coefficient in SOMD

[lohedges]

SOMD via OpenMMD.py specifies and inverse friction coefficient with units of picoseconds, e.g. here. The default value can be found via:

somd-freenrg -H | grep -A2 inverse 
inverse friction = 0.1 ps
Inverse friction time for the Langevin thermostat.

This is set in the OpenMMFrEnergyST object by using setFriction(inverse_friction.val). This is then applied to the OpenMM integrator as the friction coefficient directly, which should be in units of inverse picoseconds:

const double converted_friction = convertTo(friction.value(), picosecond);
integrator_openmm = new OpenMM::LangevinMiddleIntegrator(converted_Temperature, converted_friction, dt);

Here's the OpenMM docs showing that it should be in inverse picoseconds.

The unit for the inverse friction coefficient is correct (picoseconds) but the value is not being inverted when creating the integrator. To work around this in BioSimSpace I have applied the thermostat time constant in an inverse fashion to the other engines, but have incorrectly stripped the unit which is causing SOMD to crash on load.

[lohedges]

In fact, as with most of the default SOMD configuration output, the specification for inverse friction is wrong. This must be specified as 0.1 picosecond, otherwise it will fail, i.e.:

inverse friction = 0.1 ps

Gives:

!!!FATAL!!!
Could not understand the following parameters!

inverse friction = 0.1 ps

CANNOT CONTINUE - PROGRAM EXITING

inverse friction = 0.1

Gives:

Traceback (most recent call last):
  File "/home/lester/.conda/envs/openbiosim/share/Sire/scripts/somd-freenrg.py", line 161, in <module>
    OpenMMMD.runFreeNrg(params)
  File "/home/lester/.conda/envs/openbiosim/lib/python3.10/site-packages/sire/legacy/Tools/__init__.py", line 175, in inner
    retval = func()
  File "/home/lester/.conda/envs/openbiosim/lib/python3.10/site-packages/sire/legacy/Tools/OpenMMMD.py", line 2558, in runFreeNrg
    moves = setupMovesFreeEnergy(
  File "/home/lester/.conda/envs/openbiosim/lib/python3.10/site-packages/sire/legacy/Tools/OpenMMMD.py", line 1825, in setupMovesFreeEnergy
    Integrator_OpenMM.setFriction(
Boost.Python.ArgumentError: Python argument types in
    OpenMMFrEnergyST.setFriction(OpenMMFrEnergyST, float)
did not match C++ signature:
    setFriction(SireMove::OpenMMFrEnergyST {lvalue}, SireUnits::Dimension::PhysUnit<0, 0, 1, 0, 0, 0, 0> arg0)

Using inverse friction = 0.1 picosecond runs.

[lohedges]

@jmichel80: Could you clarify the units for inverse friction in SOMD so that I can correct the configuration generation code in BioSimSpace. Obviously I need to write "picosecond" in the configuration file for it to work, but should the value refer to the result of 1 / tau-t (GROMACS) or gamma_ln (AMBER). (Basically I will ignore the unit in the file, but I need to get the value right.)

[jmichel80]

I think this was overlooked as historically simulations were done with a Verlet Integrator. This thread discusses recommended settings for the LangevinMiddle integrator. openmm/openmm#2520
Note that the discussion may be relevant to issue #47 opened by @fjclark

[lohedges]

Thanks, that agrees with my default setting. Should we fix SOMD so that it inverts the value before passing to OpenMM, or should I just do this in BioSimSpace, i.e. treat the configuration parameter as the actual friction coefficient? I just wonder if changing it in Sire might have unintended consequences for anyone using old SOMD configuration files.

[jmichel80]

Yes this may not be explicitly set by the users so could surprise someone, I would do this in BioSimSpace.

[lohedges]

I've reported this in BioSimSpace here and fixed here. Feel free to close as required. It's still obviously incorrect in the actual SOMD code, but most users will interface with this via BioSImSpace and it probably is something we can address with the re-write. Since the problem is now documented, things should be okay.

Cheers.

[chryswoods]

I will close this as it has been fixed in BioSimSpace. We will refer back to this when somd2 development ramps up, to make sure that it is correct in the new code.