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

[Ran Friedman]

Hi Oliver,

I think the eigenvalues in NMA are not the same (there used to be a
factor of 2PI and the mass weighting). Maybe you can try my script from
the user contributions and see if you get something more reasonable (use
to flag -n to indicate that your eigenvalues are from NMA).

Ran.

oliver.kuhn at uni-duisburg-essen.de wrote:
> 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
>
>
>
>
>
>
>
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>   


-- 
------------------------------------------------------
Ran Friedman
Postdoctoral Fellow
Computational Structural Biology Group (A. Caflisch)
Department of Biochemistry
University of Zurich
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
Tel. +41-44-6355593
Email: r.friedman at bioc.unizh.ch
Skype: ran.friedman
------------------------------------------------------