[gmx-users] entropy calculation

[tarak karmakar]

Hi  Prof. David van der Spoel and Dr. Lemkul,

We had quite an interesting discussion related to the entropy calculation
in gromacs-4.5 version.
Lately, I resumed the configurational entropy calculations in the newer
(not latest) version of gromacs.

Presently, I am looking into the results of a paper [
http://ac.els-cdn.com/S0006349503747699/1-s2.0-S0006349503747699-main.pdf?_tid=2c9f47ec-e2be-11e4-ba68-00000aab0f01&acdnat=1429026894_1ff268149a3375671afa08320948af66
] and trying to reproduce Figure 15 (entropy vs simulation time). For this
I have simulated the closed ligand-free protein (as stated in the article)
in water. Conditions of the simulation (simulation time, temp, force field
etc.) are almost similar except few settings in gromacs-4.6.7.

After running the simulation for 1.5 ns, I tried to calculate the
configurational entropy for this system.
As usual I used g_covar and g_anaeig tools to do so. Few lines from the
script I use to calculate entropy are given below,

------------------------------------------------------------------------------------------
for (( N = 200; N <= 1500; N+=100 ))
   do
    echo "time =" $N "ps"
      g_covar_mpi_d -f ../protein_npt_prod -s ../protein_npt_prod -o
eigenval -v eigenvec.trr -av average.pdb -e $N << EOF
      3
      3
EOF
      g_anaeig_mpi_d -v eigenvec -entropy -temp 300
   done
grep 'Quasi Harmonic approximation' s_traj.out > QH.dat
grep 'Schlitter' s_traj.out > SC.dat
------------------------------------------------------------------------------------------

Surprisingly, I have got some doubtful plot, data are given below.

Frame Time (ps) S(J/K/mol)
400 200   2807.63
600 300   2807.63
800 400   2807.63
1000 500   2487.93
1200 600   2546.4
1400 700   2594.98
1600 800   2629.82
1800 900   2667.44
2000 1000 2697.97
2200 1100 2725.4
2400 1200 2748.65
2600 1300 2773.31
2800 1400 2792.37
3000 1500 2807.63
There are 331 C-alpha atoms that I chose to generate the covariance matrix.
Interestingly, only above 800 (I suppose above (331x3 = 993)) gmx_anaeig.c
code is behaving properly.
So, for 3N_atoms > N_frames, eigenvalues with 'zero' magnitude are creating
problem!

Does this indicate that the old bug https://gerrit.gromacs.org/#/c/3447/ is
still there in version-4.6.7?
This time instead of 'nan' some value is getting printed.

Regards,
Tarak