Zero interaction energy from mdrun -rerun in GMX 2020 vs non zero energies in GMX 2018

GROMACS version: 2020.6 / 2018.8
GROMACS modification: No

Dear Gromacs Users,

My task:
I am (re-)calculating the interaction energy (LJ SR + Coul SR) of two residues in a trajectory using mdrun -rerun with energy groups GRPX and GRPY containing indices of resid X and resid Y in separate index groups.

My problem:
In GMX version 2020 these interactions are all 0.0, while using the same procedure with GMX version 2018 results in reasonable (non zero) interaction energies. Interestingly, the self interactions residueX-residueX are non zero but huge (>80000 kJ/mol in 2020 vs 270 kJ/mol in 2018) in version 2020.

Solution?
I’d guess that it is somehow related to some changes in the way -rerun works in GMX 2020, though I could not find anything in the release notes. Thanks in advance for any information!

Best
Fabian

Sounds like your 2020 calculation was done on a GPU, which does not support energy groups. Make sure you’re running only on CPU for mdrun -rerun.

That’s a good point, as GMX actually notes that in the logs, but both versions are compiled without gpu support and use thread_mpi.

From the log file (shortened):

GROMACS: gmx mdrun, version 2020.6
Executable: /home/f/f_kell07/software/gromacs_analysis-2020.6/bin/gmx
Data prefix: /home/f/f_kell07/software/gromacs_analysis-2020.6
Working dir: /scratch/tmp/f_kell07/analysis_mclattice/dppc_330/testruns
Process ID: 119139
Command line:
gmx -nocopyright mdrun -ntomp 8 -ntmpi 8 -v -rerun /scratch/tmp/f_kell07/mdfiles_mclattice/md_trj/dppc_330_whole.xtc -s mdrerun_resid1_leaflet.tpr

GROMACS version: 2020.6
Verified release checksum is 2f568d8884e039acbc6b68722432516e0628be00c847969b7c905c8b53ef826f
Precision: single
Memory model: 64 bit
MPI library: thread_mpi
OpenMP support: enabled (GMX_OPENMP_MAX_THREADS = 64)
GPU support: disabled
SIMD instructions: AVX_512
FFT library: fftw-3.3.8-sse2-avx-avx2-avx2_128-avx512
RDTSCP usage: enabled
TNG support: enabled
Hwloc support: hwloc-2.2.0
Tracing support: disabled
C compiler: /Applic.HPC/Easybuild/skylake/2020a/software/GCCcore/9.3.0/bin/cc GNU 9.3.0
C compiler flags: -mavx512f -mfma -fexcess-precision=fast -funroll-all-loops -O3 -DNDEBUG
C++ compiler: /Applic.HPC/Easybuild/skylake/2020a/software/GCCcore/9.3.0/bin/c++ GNU 9.3.0
C++ compiler flags: -mavx512f -mfma -fexcess-precision=fast -funroll-all-loops -fopenmp -O3 -DNDEBUG

Running on 1 node with total 36 cores, 72 logical cores
Hardware detected:
CPU info:
Vendor: Intel
Brand: Intel(R) Xeon(R) Gold 6140 CPU @ 2.30GHz
Family: 6 Model: 85 Stepping: 4
Features: aes apic avx avx2 avx512f avx512cd avx512bw avx512vl clfsh cmov cx8 cx16 f16c fma hle htt intel lahf mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdrnd rdtscp rtm sse2 sse3 sse4.1 sse4.2 ssse3 tdt x2apic
Number of AVX-512 FMA units: 2

Multiple energy groups is not implemented for GPUs, falling back to the CPU. For better performance, run on the GPU without energy groups and then do gmx mdrun -rerun option on the trajectory with an energy group .tpr file.
Input Parameters:
integrator = md
tinit = 0
dt = 0.002
nsteps = 0
init-step = 0
simulation-part = 1
nstcalcenergy = 100
nstenergy = 100
cutoff-scheme = Verlet
nstlist = 20
pbc = xyz
periodic-molecules = false
verlet-buffer-tolerance = 0.005
rlist = 1.2
coulombtype = PME
coulomb-modifier = Potential-shift
rcoulomb-switch = 0
rcoulomb = 1.2
epsilon-r = 1
epsilon-rf = inf
vdw-type = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1
rvdw = 1.2
DispCorr = No
table-extension = 1
tcoupl = Nose-Hoover
nsttcouple = 20
nh-chain-length = 1
print-nose-hoover-chain-variables = false
pcoupl = Parrinello-Rahman
pcoupltype = Semiisotropic
nstpcouple = 20
tau-p = 5

And here is the output from gmx energy…
…from GMX 2018

@ title “GROMACS Energies”
@ xaxis label “Time (ps)”
@ yaxis label “(kJ/mol)”
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend “Potential”
@ s1 legend “Coul-SR:resid_341-resid_341”
@ s2 legend “LJ-SR:resid_341-resid_341”
@ s3 legend “Coul-14:resid_341-resid_341”
@ s4 legend “LJ-14:resid_341-resid_341”
@ s5 legend “Coul-SR:resid_341-leaflet”
@ s6 legend “LJ-SR:resid_341-leaflet”
@ s7 legend “Coul-14:resid_341-leaflet”
@ s8 legend “LJ-14:resid_341-leaflet”
@ s9 legend “Coul-SR:resid_341-interleaflet”
@ s10 legend “LJ-SR:resid_341-interleaflet”
@ s11 legend “Coul-14:resid_341-interleaflet”
0.000000 -569471.812500 347.489532 -58.237930 -574.714905 18.380203 -174.961395 -257.103912 0.000000 0.000000 1.774202 -22.761156 0.000000
1000.000000 -572068.000000 228.242523 -47.072948 -449.503845 34.441422 -134.704391 -294.663940 0.000000 0.000000 1.343589 -15.537995 0.000000
2000.000000 -566054.812500 270.467377 -41.391750 -478.206879 27.504395 -120.776375 -233.636124 0.000000 0.000000 2.779344 -60.008224 0.000000
3000.000000 -568577.000000 259.353149 -47.965710 -460.656738 40.298637 -169.074280 -266.141602 0.000000 0.000000 4.071910 -74.605553 0.000000
4000.000000 -568047.437500 306.153503 -55.144253 -511.083923 34.962753 -97.457947 -227.578644 0.000000 0.000000 2.914246 -53.402748 0.000000

From version 2020:

@ title “GROMACS Energies”
@ xaxis label “Time (ps)”
@ yaxis label “(kJ/mol)”
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend “Potential”
@ s1 legend “Coul-SR:resid_341-resid_341”
@ s2 legend “LJ-SR:resid_341-resid_341”
@ s3 legend “Coul-14:resid_341-resid_341”
@ s4 legend “LJ-14:resid_341-resid_341”
@ s5 legend “Coul-SR:resid_341-leaflet”
@ s6 legend “LJ-SR:resid_341-leaflet”
@ s7 legend “Coul-14:resid_341-leaflet”
@ s8 legend “LJ-14:resid_341-leaflet”
@ s9 legend “Coul-SR:resid_341-interleaflet”
0.000000 -569472.875000 0.000000 0.000000 -574.715027 18.380203 0.000000 0.000000 0.000000 0.000000 0.000000
1000.000000 -572068.375000 0.000000 0.000000 -449.503754 34.441410 0.000000 0.000000 0.000000 0.000000 0.000000
2000.000000 -566056.687500 0.000000 0.000000 -478.206848 27.504377 0.000000 0.000000 0.000000 0.000000 0.000000
3000.000000 -568577.687500 0.000000 0.000000 -460.656616 40.298637 0.000000 0.000000 0.000000 0.000000 0.000000
4000.000000 -568048.500000 0.000000 0.000000 -511.084076 34.962761 0.000000 0.000000 0.000000 0.000000 0.000000
5000.000000 -567961.000000 0.000000 0.000000 -534.145630 22.905722 0.000000 0.000000 0.000000 0.000000 0.000000
6000.000000 -568862.187500 0.000000 0.000000 -478.079651 22.513742 0.000000 0.000000 0.000000 0.000000 0.000000
7000.000000 -569597.937500 0.000000 0.000000 -466.317657 31.928671 0.000000 0.000000 0.000000 0.000000 0.000000

Note that exactly the same input was used for both reruns. For GMX 2020 I tried using tpr files generated with 2020 and, additionally, generated with 2018.
I also tried manually disabling all GPU calculations (with mdrun flags like -bonded cpu -nb cpu …).

Is this maybe a bug?

Was there any resolution of this problem? I’m seeing the same thing when I try and calculate interaction energies for a subset of water molecules. gmx hbonds works fine on the subset, but gmx energy -rerun -nb cpu gives all zeros. The same problem occurs with gromacs 2021.2 and 2022.

Hey,
I haven’t received feedback on this issue on any platform… In the end, I went back calculating the energies using Gromacs version 2018.
Best
Fabian

Hi Fabian
Thank you for your reply. I have applied to join the developers mailing list to see if I can help with this problem. Gromacs is a spectacular set of programs, so I’d love to squash what appears to be a real bug :) I’ll let you know if I make any headway.

Best Wishes
Guy Vigers

Hi Fabian

The developers have found and fixed my problem. Many thanks to @berkhess and @erik.lindahl for the great work!

The fix should be out in 2022.1 Hopefully it will work for your case too.

Best Wishes

Guy Vigers

Great! Thanks for your efforts to find a working solution and also thanks to the developers! :-)