[gmx-users] different result for entropy with normal mode analysis and schlitter-approximation

[oliver.kuhn at uni-duisburg-essen.de]

Dear Gromacs Users,
I'm trying to calculate entropies from a md trajectory using g_anaeig.
There are two ways to go (question at bottom ;-):

1. NMA and quasi-harmonic approximation: Use a bunch of snapshots (maybe
5-20), minimize each of them to very low maximum forces, calculate the
hessian matrix, diagonalize and use g_anaeig to calculate the entropy from
the resulting eigenvector-matrix assuring that there are no negative
eigenvalues in the eigenvectors 7 to N (first six eigenvectors will not be
part of the calculation). - as follows:

# Energy Minimization
grompp_d -f em_nma.mdp -t md.fitted.trr -time $t -c md.gro -p protein.top
-o $t.em.tpr
mdrun_d -v -deffnm $t.em -table table6-12_4r_doublePrecision.xvg -tablep
table6-12_4r_doublePrecision.xvg

# Hessian Matrix
grompp_d -f nma.mdp -t $t.em.trr -c md.gro -p protein.top -o $t.hessian.tpr
mdrun_d -v -deffnm $t.hessian -table table6-12_4r_doublePrecision.xvg
-tablep table6-12_4r_doublePrecision.xvg

# Diagonalization of the Hessian
g_nmeig_d -f $t.hessian.mtx -s $t.hessian.tpr -first 1 -last 10000 -v
$t.eigenvec.trr

# Entropy calculation - vibrational (without first 6 modes)
g_anaeig_d -v $t.eigenvec.trr -f $t.em.trr -s $t.hessian.tpr -temp 298.15
-nevskip 6 -entropy 2>&1 | tee out.anaeig.Svib.$t

grep 'The Entropy due to the Quasi Harmonic approximation is'
out.anaeig.Svib.$t | awk '{print $10}' >> result/Svib.nma

I use distance-dependent dielectric e=4r, but that doesn't make much
difference.

2. Schlitter approximation based on covariance: Use all snapshots of the
md trajectory, calculate the covariance matrix (g_covar), - diagonalized
matrix will be returned -, and subsequently calculate the entropy with
g_anaeig. - as follows:

# covariance matrix as time average over configurations
g_covar_d -f md$i.trr -s md.gro -v md$i.eigenvec.trr -mwa -av
average.$i.pdb -ascii covar.$i -xpm covar.$i -xpma covara.$i -l covar.$i
-o md$i.eigenval.xvg <<- EOF
0
0
EOF

# Analysis of the principal components (and entropy calculation)
g_anaeig -v md$i.eigenvec.trr -f md$i.trr -s md.gro -first 1 -last -1
-entropy > out.anaeig.schlitter.$i

grep 'The Entropy due to the Schlitter formula is' out.anaeig.schlitter.$i
| awk '{print $9}' >> result/Svib.schlitter

Somebody before mentioned, he would like to have the undiagonalized
covariance matrix as input for the entropy calculation, I think, that
doesn't make a difference, am I right?

So, practically, I tried to reproduce entropy from Schlitter 1993. A
simulation of a deca-alanine-helix in vacuo in the old gmx force-field
with vdw-cut-off etc. and I could reproduce the value of ca. 700
kJoule/mol K with the Schlitter approximation.
And now the question, why don't I get the same range of values when doing
normal-mode analysis (as described above)?

values of the Schlitter-approximation (for different simulation lengths):
667.365
685.594
681.259
680.269
values given by the quasi-harmonic approximation when calculating from
covariance:
582.731
596.97
590.71
589.07

values from NMA and quasi-harmonic approximation (for 3 snapshots):
21662.9
21674.9
21662.9

So, there is a factor of round 30 between hessian- and covariance-based
entropy!?

I'm totally stuck with this.
If anybody has experience with this phenomenon, any help is appreciated.
Thanx in advance

Oliver Kuhn