[gmx-developers] Re: printing to edr in source code

Inon Sharony InonShar at TAU.ac.IL
Sun Sep 13 00:56:47 CEST 2009

```Hi Berk (and anyone else following along)!

My system is a couple of atoms with weight m=32.06 g/mole (each)
connected by a harmonic bond with equilibrium length b0=0.238 nm and
force constant k=680,000.0 kJ/mole per nm. There are no other
interactions.I'm using g_energy to calculate the potential energy
(Vgmx), and also using g_traj to calculate it according to
Vcalc=0.5*k*(|xL-xR|-b0)^2 where xL and xR are the atoms' positions.
Initially the atoms are located such that their x and y coordinates
are equal, so the problem is essentially one-dimensional.

It turns out that the energetics (potential, kinetic and total) I got
for this oscillator are different whether I use g_energy with the
original trajectory, the rerun, or if I use g_traj to get the
positions and velocities, and try to calculate the energetics by myself.

#g_energy from md
#	time	potential		kinetic		total
0.000000000000   21.760000000001    0.010738129589   21.770738129590
0.000152532426   21.717068677863    0.053648265968   21.770716943831
0.000305064852   21.631375477836    0.139299178097   21.770674655933
0.000457597278   21.503258619639    0.267352813159   21.770611432798

#g_energy from rerun
#	time	potential		kinetic		total
0.000000000000   21.760000000001    0.021476259178   21.781476259179
0.000152532426   21.717068677863    0.085820272759   21.802888950622
0.000305064852   21.631375477836    0.192778083435   21.824153561272
0.000457597278   21.503258619639    0.341927542883   21.845186162522

#energies calculated by me from g_traj positions and velocities
#time	potential	kinetic	total
0.0e+00 21.760000 0.000000 21.760000
1.5e-04 21.760000 0.021476 21.781476
3.1e-04 21.651336 0.085820 21.737156
4.6e-04 21.542944 0.192778 21.735722

#positions and velocities from g_traj
#d	displacement from equilibrium bond length
#v	velocity of one of the two atoms (the other is just -1.0 times this)
#F	force (Newton's third law checks out)
#d_F	displacement from force (F=-k*d_F)
#d_E	displacement from energy (E=0.5*k*d_E^2)
#time	d	v	F	 d_F	d_E
0.0e+00 -0.008000 0.000000 5440.000000 -0.008000 0.008000
1.5e-04 -0.008000 -0.025882 5434.630000 -0.007992 0.007992
3.1e-04 -0.007980 -0.051738 5423.900000 -0.007976 0.007976
4.6e-04 -0.007960 -0.077544 5407.810000 -0.007953 0.007953

*****PROCEDURE SCRIPT*****

grompp_mpi -f md -c b4eq -n -maxwarn 0 -quiet

mdrun_mpi -e md -c md -g md -o md -pd -quiet -v

# extract Pos. Res. and Tot. Energies from ener.edr and output energy.xvg
echo | echo 4 5 6 | g_energy_mpi -f md -dp -quiet

# extract phase data (velocity) for the index groups named "1" and "0"
# note that the *.trr and not the *.xtc file needs to be read in order
to get velocities
echo 2 | g_traj_mpi -f md.trr -n -ov v_1 -quiet    # select index
group 2 in interactive menu of g_traj

echo 1 | g_traj_mpi -f md.trr -n -ov v_0 -quiet    # select index
group 1 in interactive menu of g_traj

echo 2 | g_traj_mpi -f md.trr -n -ox x_1 -quiet    # select index
group 2 in interactive menu of g_traj

echo 1 | g_traj_mpi -f md.trr -n -ox x_0 -quiet    # select index
group 1 in interactive menu of g_traj

# rerun last trajectory, and throw out all energetic information of
interactions other than relevant pair (Au-S)
grompp_mpi -f md_rerun -c md -n -p topol_rerun.top -t md -maxwarn 1
-quiet  # preprocess using existing traj.trr file from

previous md run, and "rerun" options in md_rerun.mdp & topol_rerun.top

mdrun_mpi  -rerun md.trr -g rerun -e rerun -o rerun -pd -quiet -v
# rerun using traj.trr and output force data into rerun.trr

#energetics from rerun trajectory
echo | echo 5 6 7 | g_energy_mpi -f rerun -o rerun -dp -quiet

# extract force data for the index group named "1"
echo 2 | g_traj_mpi -f rerun.trr -n -of f_1 -quiet # select index
group 2 in interactive menu of g_traj

# extract force data for the index group named "0"
echo 1 | g_traj_mpi -f rerun.trr -n -of f_0 -quiet # select index
group 1 in interactive menu of g_traj

*****MD.MDP FILE*****

integrator      =       md              ; stochastic dynamics (solve
Langevin equation)
dt              =       0.00015253242594417709812523865069445	; this
is half the value in the next line
;dt             =       0.00030506485188835419625047730138889	; this
is one hundredth of a period of the bond oscillation
nsteps          =       3
;
;------------------------------------------------------------------------------------------------
;
pbc             =       no              ; no periodic boundary
conditions (only one molecule)
;comm_mode      =       Angular         ; remove center of mass
translatory AND rotational motion
nstcomm         =       1               ; number of steps between
removal of center of mass motion
;
;------------------------------------------------------------------------------------------------
; number of steps between writing out of different phase data and energetics
nstxout         =       1
nstvout         =       1
nstfout         =       1
nstenergy       =       1
; groups to write energetics of
energygrps      =       AU      S
;
;-------------------------------------------------------------------------------------------------
; neighbor list search
nstype          =       simple          ; particle (as opposed to
domain) decomposition
;rlist          =       1.1             ; radius of neighbor list search
;
coulombtype     =       cut-off
;rcoulomb       =       1.2
;coulombtype    =       PME-switch
vdwtype         =       Shift
;
;-------------------------------------------------------------------------------------------------
;
; Langevin dynamics temperature coupling
;
tc-grps         =       AU      S
tau_t           =       100000  100000  ; mass/gamma
ref_t           =       10      10      ; reference (bath) temperature

;-------------------------------------------------------------------------------------------------
;
; Velocity generation
gen_vel         =       no
gen_temp        =       0               ; (300) [K]

*****MD_RERUN.MDP FILE*****

integrator      =       md              ; solve Newtonian equations of motion
dt              =       0.00015253242594417709812523865069445
nsteps          =       3
;
;------------------------------------------------------------------------------------------------
;
nstxout         =       1
nstvout         =       1
nstfout         =       1
nstenergy       =       1
;
energygrps      =       AU      S
;
;-------------------------------------------------------------------------------------------------
;
nstype          =       simple
rlist           =       1.1
;
coulombtype     =       PME-switch
vdwtype         =       Shift
;
;-------------------------------------------------------------------------------------------------
;
;Energy group exclusions for rerun
energygrp_excl  =       AU AU   S S

*****TARJECTORY FILE*****

md-11535530CFD49-1.trr frame 0:
natoms=         2  step=         0  time=0.0000000e+00  lambda=         0
box (3x3):
box[    0]={ 6.37511e+00,  0.00000e+00,  0.00000e+00}
box[    1]={ 0.00000e+00,  6.37511e+00,  0.00000e+00}
box[    2]={ 0.00000e+00,  0.00000e+00,  6.37511e+00}
x (2x3):
x[    0]={ 3.15000e+00,  3.15000e+00,  3.26500e+00}
x[    1]={ 3.15000e+00,  3.15000e+00,  3.03500e+00}
v (2x3):
v[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
v[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
f (2x3):
f[    0]={ 0.00000e+00,  0.00000e+00,  5.44000e+03}
f[    1]={ 0.00000e+00,  0.00000e+00, -5.44000e+03}
md-11535530CFD49-1.trr frame 1:
natoms=         2  step=         1  time=1.5253243e-04  lambda=         0
box (3x3):
box[    0]={ 6.37511e+00,  0.00000e+00,  0.00000e+00}
box[    1]={ 0.00000e+00,  6.37511e+00,  0.00000e+00}
box[    2]={ 0.00000e+00,  0.00000e+00,  6.37511e+00}
x (2x3):
x[    0]={ 3.15000e+00,  3.15000e+00,  3.26500e+00}
x[    1]={ 3.15000e+00,  3.15000e+00,  3.03500e+00}
v (2x3):
v[    0]={ 0.00000e+00,  0.00000e+00,  2.58820e-02}
v[    1]={ 0.00000e+00,  0.00000e+00, -2.58820e-02}
f (2x3):
f[    0]={ 0.00000e+00,  0.00000e+00,  5.43463e+03}
f[    1]={ 0.00000e+00,  0.00000e+00, -5.43463e+03}
md-11535530CFD49-1.trr frame 2:
natoms=         2  step=         2  time=3.0506485e-04  lambda=         0
box (3x3):
box[    0]={ 6.37511e+00,  0.00000e+00,  0.00000e+00}
box[    1]={ 0.00000e+00,  6.37511e+00,  0.00000e+00}
box[    2]={ 0.00000e+00,  0.00000e+00,  6.37511e+00}
x (2x3):
x[    0]={ 3.15000e+00,  3.15000e+00,  3.26501e+00}
x[    1]={ 3.15000e+00,  3.15000e+00,  3.03499e+00}
v (2x3):
v[    0]={ 0.00000e+00,  0.00000e+00,  5.17384e-02}
v[    1]={ 0.00000e+00,  0.00000e+00, -5.17384e-02}
f (2x3):
f[    0]={ 0.00000e+00,  0.00000e+00,  5.42390e+03}
f[    1]={ 0.00000e+00,  0.00000e+00, -5.42390e+03}
md-11535530CFD49-1.trr frame 3:
natoms=         2  step=         3  time=4.5759728e-04  lambda=         0
box (3x3):
box[    0]={ 6.37511e+00,  0.00000e+00,  0.00000e+00}
box[    1]={ 0.00000e+00,  6.37511e+00,  0.00000e+00}
box[    2]={ 0.00000e+00,  0.00000e+00,  6.37511e+00}
x (2x3):
x[    0]={ 3.15000e+00,  3.15000e+00,  3.26502e+00}
x[    1]={ 3.15000e+00,  3.15000e+00,  3.03498e+00}
v (2x3):
v[    0]={ 0.00000e+00,  0.00000e+00,  7.75438e-02}
v[    1]={ 0.00000e+00,  0.00000e+00, -7.75438e-02}
f (2x3):
f[    0]={ 0.00000e+00,  0.00000e+00,  5.40781e+03}
f[    1]={ 0.00000e+00,  0.00000e+00, -5.40781e+03}

You have not told system details nor how you calculate the energies from
the output, so I have no clue what could be the problem. In any MD
integrator the energies and forces are simply calculated directly from
the coordinates. In Gromacs these coordinates are written to trr. If
you do mdrun -rerun with a trr file with x(t) and v(t-dt/2) (which is
what mdrun writes) you get the exact energies back, including kinetic,
virial and pressure, up to numerical rounding errors.Note that the
total energy does not say anything about the forces.Berk

```