[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}
Reply to:
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
More information about the gromacs.org_gmx-developers
mailing list