[gmx-users] entropy calculation
tarak karmakar
tarak20489 at gmail.com
Wed Apr 15 16:25:18 CEST 2015
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
More information about the gromacs.org_gmx-users
mailing list