[gmx-users] Bonds energies calculation

[Alexander]

Hi gromacs developers,
I have simply topology for Oxygen molecule with bond description:
 [ bonds ]
; atom0   atom1      optimal_length    energy_constant
    O      O00           0.10860          51730.6
I have a frame from trajectory(in pdb format):
MODEL        2
ATOM      1  O   LIG     1      30.001  29.984  30.000  1.00  0.00
ATOM      2  O00 LIG     1      29.999  30.016  31.080  1.00  0.00
TER
ENDMDL
I have try to calculate the bond energy in current frame:
length=sqrt((x[O]-x[O00])^2+(y[O]-y[O00])^2+(z[O]-z[O00])^2)
length=sqrt(0.000004+0.001024+1.1664)=1.08048 A -> 0.108048nm
the bond deviation:
deviation=(length-optimal_length)=0.0005524nm
and the bond energy is:
BOND_ENERGY=energy_constant*deviation^2
BOND_ENERGY=51730.6*(0.0005524)^2=0.0157854 kJ/mol
but, the same energy extracted by g_energy is:
energy.xvg
    0.006000    0.007892
It seems that round(0.0157854*0.5)=0.007892, but where 0.5 coefficient
taken from?

Thank you,
Best regards,

Alexander Yakovenko
Institute of Molecular Biology & Genetic of NAS of Ukraine
03143
acad.Zabolotnogo str. 150
Kiev
Ukraine
______________________________________________
E-mail: yakovenko_a at ukr.net