Why is it different every time the interaction energy is calculated?

User discussions
,
1

GROMACS version:2021
GROMACS modification: Yes/No
I want to calculate the interaction energy between the two species. But every time I calculate the same simulation, the interaction energy varies between 250 and 600. I use the NVT ensemble. There is a vacuum inside the simulation box, and I use 148 core CPU.
forcefield=charmm36

gmx mdrun -v -deffnm nvt -nb cpu -pme cpu -pmefft cpu -bonded cpu -update cpu -pin auto -tunepme yes -reseed 123456789 -dlb yes

gmx grompp -f Interaction_energy.mdp -c nvt.gro -p topol.top -o Interaction.tpr
gmx mdrun -s Interaction.tpr -rerun nvt.xtc -nt 1

gmx energy -f ener.edr -o Interactionenergy.xvg > Interaction.log

nvt.mdp
integrator = md-vv
nsteps = 10000000
dt = 0.001
nstlist = 10
rlist = 1.4
cutoff-scheme = Verlet
vdw-type = PME
rvdw = 1.4
coulombtype = P3M-AD
DispCorr = Ener
lj-pme-comb-rule = Lorentz-Berthelot
rcoulomb = 1.4
tcoupl = v-rescale
tc-grps = system
tau-t = 0.1
ref-t = 353
pbc = xyz
gen-vel = yes
gen-temp = 353
gen-seed = 123456789
ld-seed = 123456789
constraints = none
freezegrps = ****
freezedim = Y Y Y Y Y Y
nstxtcout = 10000
nstxout = 10000
nstenergy = 10000

Interaction_energy.mdp

integrator = md-vv
nsteps = 10000000
dt = 0.001
nstlist = 10
rlist = 1.4
cutoff-scheme = Verlet
vdw-type = PME
rvdw = 1.4
coulombtype = P3M-AD
DispCorr = Ener
lj-pme-comb-rule = Lorentz-Berthelot
rcoulomb = 1.4
tcoupl = v-rescale
tc-grps = system
tau-t = 0.1
ref-t = 353
pbc = xyz
gen-vel = yes
gen-seed = 123456789
ld-seed = 123456789
gen-temp = 353
constraints = none
energygrps = ** **
nstxtcout = 10000
nstxout = 10000
nstenergy = 10000

2

Why nobody answers me? The question is not complicated.

3

@jalemkul
@alevilla
@hess

4

Before I answer I would like to say that your behavior is very unpolite. GROMACS is not a commercial software that you can demand support for, but free open-source software and many of the people answering questions on this forum are doing that to be helpful and are not begin paid for answering.

The answer to your question is that molecular dynamics is chaotic and mdrun is not binary reproducible (unless you use the -reprod flag and even then runs sometimes are not reproducible). So trajectories, also starting from the same tpr file, diverge exponentially in time.

5

Hi Reza
I have a question: Are you re-running both nvt and interaction.tpr each time? If you repeat just the second part (interaction.tpr) on a single nvt run do you get consistent result? In other words, at which stage is the inconsistency creeping in?

Best Wishes
Guy Vigers

6

Thank you very much, Hess. I did not mention you intentionally. I just saw that you answered a few questions a few hours ago, and I thought you could help solve this problem. Because I posted this question a long time ago, and no one answered. I apologize for my behavior that upset you.
I had a question about the -reprod flag. Does the interaction energy change decrease when I use this flag? The interaction energy I calculated is a variable range of 250-600. Does this change range decrease?

7

Dear Guy Vigers,
the command I use every time is below:

first simulation

gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt1.tpr
gmx mdrun -v -deffnm nvt1 -nb cpu -pme cpu -pmefft cpu -bonded cpu -update cpu -pin auto -tunepme yes -reseed 123456789 -dlb yes
gmx grompp -f Interaction_energy.mdp -c nvt1.gro -p topol.top -o Interaction.tpr
gmx mdrun -s Interaction.tpr -rerun nvt1.xtc -nt 1
gmx energy -f ener1.edr -o Interactionenergy1.xvg > Interaction1.log

Second one

gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt2.tpr
gmx mdrun -v -deffnm nvt2 -nb cpu -pme cpu -pmefft cpu -bonded cpu -update cpu -pin auto -tunepme yes -reseed 123456789 -dlb yes
gmx grompp -f Interaction_energy.mdp -c nvt2.gro -p topol.top -o Interaction.tpr
gmx mdrun -s Interaction.tpr -rerun nvt2.xtc -nt 1
gmx energy -f ener2.edr -o Interactionenergy2.xvg > Interaction2.log

8

So if you just rerun the second calculation you get the same results, right?
That tells you that your nvt runs are not identical. Have you looked at your coordinates at the ends of the nvt runs? Are they identical? I would guess not.
This could be because of numerical instabilities or because of different random seeds. See hess’s answer above for a more precise answer.
Guy

9

The interaction1 is not equal interaction2. The final coordinate about the same.