Unable to Run in Microcanonical Ensemble

[fjclark]

Hello,

I am trying to run a simulation in the NVE ensemble to check that energy is conserved when I apply strong (200 kcal mol-1) restraints along with HMR and a 4 fs timestep. Initially, I tried without any restraints: my input is here. I tried two different inputs:

1. Specifying andersen = False in the input.

Looking at the debug output and checking the energies showed that the Andersen thermostat was in use. This was not unexpected as @jmichel80 mentioned that this was broken. I attempted to get round this with:

2. andersen = true with andersen frequency = 0

However, energies changed as before.

Could you please suggest how I could run in NVE with SOMD?

Thanks very much.

[lohedges]

Hi there,

If I recall correctly andersen is a completely redundant option (it triggers no code) and you should use thermostat for the option. Try setting thermostat = False and see what happens.

The following is from the parameter section of OpenMMD.py.

andersen = Parameter("thermostat", True,
    """Whether or not to use the Andersen thermostat (needed for NVT or NPT simulation).""")

barostat = Parameter("barostat", True, """Whether or not to use a barostat (needed for NPT simulation).""")

andersen_frequency = Parameter("andersen frequency", 10.0, """Collision frequency in units of (1/ps)""")

barostat_frequency = Parameter("barostat frequency", 25,
    """Number of steps before attempting box changes if using the barostat.""")

(The variable is named andersen, but the config parameter is thermostat.)

[lohedges]

Annoyingly, SOMD has no real error checking on invalid config file entries so things like this are currently silently ignored.

[jmichel80]

Yes @lohedges is correct. Also you may have to edit the crd rst7 input to remove box dimensions and run with a non periodic cutoff to avoid complaints about lack of thermostat.

[fjclark]

Thanks very much @lohedges @jmichel80.

With the original input I obtained the following:

I repeated with thermostat = False, and made sure that the Andersen Thermostat set line disappeared from the debug output. However, I obtained a very similar plot:

I also reran as above but removed box dimensions from the rst7 input and changed the cutoff type to cutoff type = cutoffperiodic. This ran for one cycle:

Before crashing:

Cycle = 2

Starting /home/finlayclark/sire_cartesian_euler.app/bin/somd-freenrg: number of threads equals 20 Traceback (most recent call last): File "/home/finlayclark/sire_cartesian_euler.app/share/Sire/scripts/somd-freenrg.py", line 146, in OpenMMMD.runFreeNrg(params) File "/home/finlayclark/sire_cartesian_euler.app/lib/python3.7/site-packages/Sire/Tools/init.py", line 176, in inner retval = func() File "/home/finlayclark/sire_cartesian_euler.app/lib/python3.7/site-packages/Sire/Tools/OpenMMMD.py", line 2067, in runFreeNrg system = moves.move(system, nmoves.val, True) UserWarning: Exception 'SireError::programbug' thrown by the thread 'master:main'. The running total energy (1.41121e+15 kcal mol-1) disagrees with the total energy as calculated from scratch (1.41121e+15 kcal mol-1). Energy components which show disagreements greater than 0.1 kcal mol-1 are listed below:E{total}: 1.41121e+15 versus 1.41121e+15 kcal mol-1. Difference equals 0.25 kcal mol-1

Thrown from FILE: /home/finlayclark/software/devel/Sire/corelib/src/libs/SireMove/moves.cpp, LINE: 582, FUNCTION: void SireMove::Moves::postCheck(SireSystem::System&) const Backtrace ( 0) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libSireError.so.2022 ([0x7f1c2ed08e6e] ++0x2e) -- SireError::getBackTrace()

( 1) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libSireError.so.2022 ([0x7f1c2ed06f53] ++0x83) -- SireError::exception::exception(QString, QString)

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libSireMove.so.2022(+0x10f38e) [0x7f1c32bcb38e] /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libSireMove.so.2022(+0xa9741) [0x7f1c32b65741] ( 4) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libSireMove.so.2022 ([0x7f1c32c63bba] ++0x33a) -- SireMove::WeightedMoves::move(SireSystem::System const&, int, bool)

/home/finlayclark/sire_cartesian_euler.app/lib/python3.7/site-packages/Sire/Move/_Move.so(+0x27e212) [0x7f1c1979b212] ( 6) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libboost_python37.so.1.74.0 ([0x7f1c32d8a4c5] ++0x295) -- boost::python::objects::function::call(_object, _object) const

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libboost_python37.so.1.74.0(+0x1f669) [0x7f1c32d8a669] ( 8) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libboost_python37.so.1.74.0 ([0x7f1c32d8fc63] ++0x63) -- boost::python::detail::exception_handler::operator()(boost::function0 const&) const

/home/finlayclark/sire_cartesian_euler.app/lib/python3.7/site-packages/Sire/Error/_Error.so(+0x2a51) [0x7f1c2e6e2a51] ( 10) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libboost_python37.so.1.74.0 ([0x7f1c32d8fc37] ++0x37) -- boost::python::detail::exception_handler::operator()(boost::function0 const&) const

/home/finlayclark/sire_cartesian_euler.app/lib/python3.7/site-packages/Sire/Error/_Error.so(+0x2a93) [0x7f1c2e6e2a93] ( 12) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libboost_python37.so.1.74.0 ([0x7f1c32d8fb7f] ++0x3f) -- boost::python::handle_exception_impl(boost::function0)

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libboost_python37.so.1.74.0(+0x1d3c3) [0x7f1c32d883c3] ( 14) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32f97dd9] ++0xd9) -- _PyObject_FastCallKeywords

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0(+0x643c3) [0x7f1c32e0e3c3] ( 16) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32e14c77] ++0x64d7) -- _PyEval_EvalFrameDefault

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0(+0x6ddea) [0x7f1c32e17dea] /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0(+0x64416) [0x7f1c32e0e416] ( 19) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32e12d92] ++0x45f2) -- _PyEval_EvalFrameDefault

( 20) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32ecbe80] ++0x920) -- _PyEval_EvalCodeWithName

( 21) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32f970ea] ++0x8a) -- _PyFunction_FastCallKeywords

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0(+0x64416) [0x7f1c32e0e416] ( 23) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32e14c77] ++0x64d7) -- _PyEval_EvalFrameDefault

( 24) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32ecbe80] ++0x920) -- _PyEval_EvalCodeWithName

( 25) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32ecb29f] ++0x3f) -- PyEval_EvalCodeEx

( 26) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32ecc93c] ++0x1c) -- PyEval_EvalCode

( 27) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32e99446] ++0xb6) -- PyRun_FileExFlags

( 28) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32e9b5cf] ++0xff) -- PyRun_SimpleFileExFlags

/home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0(+0xd0d63) [0x7f1c32e7ad63] ( 30) /home/finlayclark/sire_cartesian_euler.app/bin/../lib/libpython3.7m.so.1.0 ([0x7f1c32e7afed] ++0x3d) -- Py_Main

/home/finlayclark/sire_cartesian_euler.app/bin/somd-freenrg(+0x4ee9) [0x5632be7dbee9] ( 32) /lib/x86_64-linux-gnu/libc.so.6 ([0x7f1c2e3a90b3] ++0xf3) -- __libc_start_main

/home/finlayclark/sire_cartesian_euler.app/bin/somd-freenrg(+0x5514) [0x5632be7dc514] EndTrace Exception 'SireError::programbug' thrown by the thread 'master:main'. The running total energy (1.41121e+15 kcal mol-1) disagrees with the total energy as calculated from scratch (1.41121e+15 kcal mol-1). Energy components which show disagreements greater than 0.1 kcal mol-1 are listed below:E{total}: 1.41121e+15 versus 1.41121e+15 kcal mol-1. Difference equals 0.25 kcal mol-1

Thrown from FILE: /home/finlayclark/software/devel/Sire/corelib/src/libs/SireMove/moves.cpp, LINE: 582, FUNCTION: void SireMove::Moves::postCheck(SireSystem::System&) const

Thanks

[chryswoods]

The key part of the error here is

The running total energy (1.41121e+15 kcal mol-1) disagrees with the total energy as calculated from scratch (1.41121e+15 kcal mol-1). Energy components which show disagreements greater than 0.1 kcal mol-1 are listed below:E_{total}: 1.41121e+15 versus 1.41121e+15 kcal mol-1. Difference equals 0.25 kcal mol-1

Sire runs checks during a simulation to make sure that the energy is being accumulated correctly. It does this by periodically recalculating the total energy from scratch, and then comparing this to the accumulated energy (normally a running total from a Monte Carlo calculation, but can be a single-point energy in Sire compared to what came from OpenMM).

The tolerance for the check is less than 0.25 kcal mol-1, to account for numerical imprecision.

In your case, the difference is 0.25 kcal mol-1 because this is about the precision possible for a double precision number that is 10^15 in magnitude.

The trouble is that your total energy is 10^15 in magnitude. Something, somewhere, has allowed two atoms to overlap. That both the running total (likely from OpenMM?) and Sire both independently calculate 10^15 indicates that this is a real energy, and thus the overlap is real.

If possible, take a look at the coordinates after the MD and see if there is anything overlapping. My guess is that it may be something that overlaps after the periodic boundary conditions are applied?

[fjclark]

Thanks very much - I will check for overlapping coordinates. However, my main issue is that I am still failing to conserve energy despite running with thermostat and barostat == False. Are there any other config file options I should be considering?

The above simulation with a non periodic cut-off did not crash when I reran it, and showed similar behaviour to previous runs:

Thanks

[lohedges]

Just to check, you are showing the total system energy, right (potential plus kinetic)?

[fjclark]

Sorry, stupid mistake - I was in fact plotting the potential energies only.

I've made a couple of small changes to try and print out total energy instead. If I haven't missed something in the modified code, then switching to total energy removes most of the sharp drop at the start of the runs, but the drift is still present. This is improved by decreasing the timestep. Plotting runs until they complete or crash (same error as previously):

If my understanding is correct, plotting the reduced perturbed energies of the current value of lambda shows the same contributions as the total potential, minus the bond bending, torsional, and stretching terms in force group 1. Plotting these shows stable energies (potential only), as would be expected if neglecting the high-frequency torsions and bends (all bonds constrained):

Have I missed anything when modifying the code to print out the total energy, or is it more likely that my system is unstable with a 4 fs timestep, even with HMR and all bonds constrained? Thanks

[jmichel80]

I wouldn't be surprised if protocols optimised for speed (high integration time step) show energy drifts in NVE. Thermostating hides this behaviour, but could lead to issues over long MD simulations. I suppose you have read the relevant literature e.g. https://pubs.acs.org/doi/10.1021/ct5010406

There could also be issues with cutoffs and the reaction-field implementation. The reduced energies should not include the protein-protein, protein-solvent, solvent-solvent contributions which are large and dominate the total energy.

You could try increasing the cutoffs to see if that helps as well as lowering the timestep. You could also run in vacuum although that may be problematic for a protein-ligand simulation (the ligand alone should be ok).

[fjclark]

Thanks @jmichel80 - maybe this is not the best way to check for errors caused by the introduction of strong restraints with the 4 fs timestep.