Reproducibility of energy terms

GROMACS version: 2019.6
GROMACS modification: No

Hello, I am running some tests with my system and I have run a 1 ns NVT equilibration of a simple peptide in vacuum, using the Amber14SB forcefield. The system is neutralized with sodium ions. I have run this simulation twice, simply by submitting the same job twice because I wanted to see if the individual energy terms would vary from one simulation to another. As you can see below, there is some variability, particularly in the dihedral terms, LJ (SR), Coulomb (SR) (only by 0.1 here), Coulomb reciprocal. Why is there such variability when the starting point of the simulation is the same? Shouldn’t they be the same? How critical are these variations? I have done the exact same thing with a simpler molecule (butanol) and I obtain almost exactly the same value for the individual terms except for the Coulomb reciprocal, then why does it change with the peptide? Thanks in advance for your help!

RUN 1
Energies (kJ/mol)
Bond Angle Proper Dih. Improper Dih. LJ-14
2.18151e+02 4.15799e+02 5.99572e+02 3.42114e+01 1.49068e+02
Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
2.62039e+03 9.53013e+01 -5.37452e-01 -7.64046e+03 2.87425e+01
Potential Kinetic En. Total Energy Conserved En. Temperature
-3.47976e+03 6.36320e+02 -2.84344e+03 4.00436e+04 3.00125e+02
Pres. DC (bar) Pressure (bar)
-5.79191e-02 4.97935e-01

RUN 2

Energies (kJ/mol)
Bond Angle Proper Dih. Improper Dih. LJ-14
2.16605e+02 4.18275e+02 6.15346e+02 4.17750e+01 1.50802e+02
Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
2.60034e+03 8.33783e+01 -5.37452e-01 -7.52557e+03 3.73552e+01
Potential Kinetic En. Total Energy Conserved En. Temperature
-3.36223e+03 6.36027e+02 -2.72620e+03 4.41374e+04 2.99986e+02
Pres. DC (bar) Pressure (bar)

MD simulations diverge rapidly. See Managing long simulations — GROMACS 2021.4 documentation

Thank you, I will try to revise those parameters and see if I can get numbers that are close together. Is there a way of setting the seed to one specific number?

In terms of the system itself, is there any reason why two runs (same inputs, hardwarem etc) generate almost exactly the same numbers for a simpler molecule (in this case butanol), in comparison to the peptide (it is only 11 amino acids long in my case)?

To illustrate, these are two runs for butanol in vacuum (the results for the peptide are in the original message). As mentioned, the only difference is in the Coulomb reciprocal term.

RUN 1
Energies (kJ/mol)
Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14
9.48867e+00 1.80425e+01 1.83977e+01 6.30814e+00 -6.48597e+01
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
1.04040e+00 -1.02735e-05 3.76011e+01 3.56998e+00 2.95888e+01
Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
5.21985e+01 8.17873e+01 9.32121e+04 2.98954e+02

RUN 2

Energies (kJ/mol)
Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14
9.38603e+00 1.82814e+01 1.84250e+01 6.35797e+00 -6.48728e+01
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
1.02060e+00 -1.02732e-05 3.70608e+01 4.02204e+00 2.96810e+01
Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
5.24935e+01 8.21746e+01 9.13529e+04 3.00644e+02

Thanks for your time!

Yes, that’s what gen-seed is for.

I don’t follow. All the energy terms are different. You can’t expect any two MD simulations to produce identical results. The only way to expect equivalence between different systems is to do single-point evaluations on the same coordinate file and you should get the same result, assuming the same topology and .mdp settings.

Hi Justin,

Thank you for your reply. It’s not between different systems. It’s two runs on the same system (same input files, topology, etc). I have just done that with a peptide and with a butanol molecule. The butanol molecule produces the same energies in both runs (with differences after the second decimal point), but there are larger discrepancies between the two runs for the peptide. That is what confused me.

It makes a lot sense to do the single point though. Thank you for your explanation. It did clear some things to me.

Thanks for your help and your time! =)

The number of accessible configurations for butanol is very small compared to a polypeptide. No matter what your starting point is for butanol, it’s simple to converge it to the same result with enough simulation time. For polypeptides, you’ll never get two simulations to give you the exact same ensemble and therefore the same energies. There are too many configurations that are possible.

That makes a lot of sense. Thanks for clarifying that and for your patience =)