To add the previous post. I took a small fraction of the protein (a loop), follow the same protocol, and got the results.
The Entropy due to the Quasi Harmonic approximation is 2814 J/mol K The Entropy due to the Schlitter formula is 3000.03 J/mol K Regarding the g_anaeig.c code, do we need to increase the precision of all the floating point numbers, say from 'double' to 'long double'? Tarak On Thu, May 15, 2014 at 3:10 PM, tarak karmakar <tarak20...@gmail.com>wrote: > Dear Sir, > Thank you again for guiding me towards this. Got some clue. > The g_anaeig code is > > ......................................................................................................................................... > static void calc_entropy_qh(FILE *fp,int n,real eigval[],real temp,int > nskip) > { > int i; > double hwkT,w,dS,S=0; > double hbar,lambda; > > hbar = PLANCK1/(2*M_PI); > for(i=0; (i<n-nskip); i++) { > if (eigval[i] > 0) { > lambda = eigval[i]*AMU; > w = sqrt(BOLTZMANN*temp/lambda)/NANO; > hwkT = (hbar*w)/(BOLTZMANN*temp); > dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT))); > S += dS; > if (debug) > fprintf(debug,"i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n", > i,w,lambda,hwkT,dS); > } > else { > fprintf(stderr,"eigval[%d] = %g\n",i,eigval[i]); > w = 0; > } > } > fprintf(fp,"The Entropy due to the Quasi Harmonic approximation is %g > J/mol K\n", > S*RGAS); > } > > static void calc_entropy_schlitter(FILE *fp,int n,int nskip, > real *eigval,real temp) > { > double dd,deter; > int *indx; > int i,j,k,m; > char buf[256]; > double hbar,kt,kteh,S; > > hbar = PLANCK1/(2*M_PI); > kt = BOLTZMANN*temp; > kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO); > if (debug) > fprintf(debug,"n = %d, nskip = %d kteh = %g\n",n,nskip,kteh); > > deter = 0; > for(i=0; (i<n-nskip); i++) { > dd = 1+kteh*eigval[i]; > deter += log(dd); > } > S = 0.5*RGAS*deter; > > fprintf(fp,"The Entropy due to the Schlitter formula is %g J/mol K\n",S); > } > > .......................................................................................................................... > What I got from the 'g_anaeig_d.debug' file is given below > > i = 0 w = 1.06303e+12 lam = 3.66537e-27 hwkT = 0.0270659 dS = > 4.60951 > i = 1 w = 1.50699e+12 lam = 1.82384e-27 hwkT = 0.0383696 dS = > 4.26055 > i = 2 w = 1.62965e+12 lam = 1.55961e-27 hwkT = 0.0414928 dS = > 4.18231 > i = 3 w = 2.02202e+12 lam = 1.01306e-27 hwkT = 0.0514829 dS = > 3.96662 > i = 4 w = 2.80273e+12 lam = 5.27285e-28 hwkT = 0.0713605 dS = > 3.64022 > i = 5 w = 3.1712e+12 lam = 4.11871e-28 hwkT = 0.0807422 dS = > 3.51677 > i = 6 w = 3.45782e+12 lam = 3.46419e-28 hwkT = 0.0880401 dS = > 3.43029 > i = 7 w = 3.56568e+12 lam = 3.25778e-28 hwkT = 0.0907862 dS = > 3.39959 > i = 8 w = 4.00962e+12 lam = 2.57633e-28 hwkT = 0.102089 dS = > 3.28234 > i = 9 w = 4.4611e+12 lam = 2.08125e-28 hwkT = 0.113584 dS = > 3.17575 > i = 10 w = 4.6037e+12 lam = 1.95431e-28 hwkT = 0.117215 dS = > 3.14431 > i = 11 w = 5.01062e+12 lam = 1.64977e-28 hwkT = 0.127576 dS = > 3.05972 > i = 12 w = 5.31327e+12 lam = 1.46718e-28 hwkT = 0.135282 dS = > 3.00116 > i = 13 w = 5.51169e+12 lam = 1.36345e-28 hwkT = 0.140334 dS = > 2.96455 > i = 14 w = 5.70613e+12 lam = 1.27211e-28 hwkT = 0.145284 dS = > 2.92994 > i = 15 w = 5.8637e+12 lam = 1.20466e-28 hwkT = 0.149296 dS = > 2.90275 > i = 16 w = 6.01592e+12 lam = 1.14447e-28 hwkT = 0.153172 dS = > 2.87717 > ---------------------- > ....................... > i = 996 w = 9.2246e+14 lam = 4.86757e-33 hwkT = 23.4869 dS = > 1.54425e-09 > i = 997 w = 9.51359e+14 lam = 4.57635e-33 hwkT = 24.2226 dS = > 7.62131e-10 > i = 998 w = 9.61243e+14 lam = 4.48271e-33 hwkT = 24.4743 dS = > 5.98465e-10 > i = 999 w = 9.86259e+14 lam = 4.25819e-33 hwkT = 25.1113 dS = > 3.2445e-10 > i = 1000 w = inf lam = 0 hwkT = inf dS = -nan > i = 1001 w = inf lam = 0 hwkT = inf dS = -nan > i = 1002 w = inf lam = 0 hwkT = inf dS = -nan > ...................... > .................. > i = 1370 w = inf lam = 0 hwkT = inf dS = -nan > i = 1371 w = inf lam = 0 hwkT = inf dS = -nan > i = 1372 w = inf lam = 0 hwkT = inf dS = -nan > i = 1373 w = inf lam = 0 hwkT = inf dS = -nan > i = 1374 w = inf lam = 0 hwkT = inf dS = -nan > i = 1375 w = inf lam = 0 hwkT = inf dS = -nan > i = 1376 w = inf lam = 0 hwkT = inf dS = -nan > i = 1377 w = inf lam = 0 hwkT = inf dS = -nan > i = 1378 w = inf lam = 0 hwkT = inf dS = -nan > i = 1379 w = inf lam = 0 hwkT = inf dS = -nan > n = 1386, nskip = 6 kteh = 4569.58 > Opening library file /usr/share/gromacs/top/gurgle.dat > > > It is very interesting that the code printed 'nan' after 999 th term. What > does this mean? > Is there any upper limit of number of atoms (c-alpha here)? > > Tarak > > > > On Thu, May 15, 2014 at 2:54 PM, David van der Spoel <sp...@xray.bmc.uu.se > > wrote: > >> On 2014-05-15 10:51, tarak karmakar wrote: >> >>> | Is your system far out of equilibrium? >>> >>> It should not be as I have simulated the system for 100ns. Few properties >>> to check to get the feelings of equilibrium are looking good. >>> >>> |Do the eigenvalues look reasonable? >>> Yes. File is attached. >>> >>> | Or is your system very large? >>> My system is a protein dimer. >>> >>> It creates the following covariance matrix (c-alpha atoms) >>> >>> Constructing covariance matrix (1386x1386) >>> >> That is big. >> >> Try running g_anaeig -debug 1 and check the output file g_anaeig.debug. >> Then check the source code (gmx_anaeig.c) to see whether the evaluation >> of the entropy can be made more numerically stable. >> >> >>> eigenvec.trr file is also attached. >>> >>> Mail with the attached files >>> "Needs moderator approval because of the file size." >>> >>> >>> Thanks, >>> Tarak >>> >>> >>> On Thu, May 15, 2014 at 2:19 PM, tarak karmakar <tarak20...@gmail.com >>> >wrote: >>> >>> | Is your system far out of equilibrium? >>>> >>>> It should not be as I have simulated the system for 100ns. Few >>>> properties >>>> to check to get the feelings of equilibrium are looking good. >>>> >>>> |Do the eigenvalues look reasonable? >>>> Yes. File is attached. >>>> >>>> | Or is your system very large? >>>> My system is a protein dimer. >>>> >>>> It creates the following covariance matrix (c-alpha atoms) >>>> >>>> Constructing covariance matrix (1386x1386) >>>> >>>> eigenvec.trr file is also attached. >>>> >>>> >>>> Thanks, >>>> Tarak >>>> >>>> >>>> >>>> >>>> On Thu, May 15, 2014 at 2:06 PM, David van der Spoel < >>>> sp...@xray.bmc.uu.se >>>> >>>>> wrote: >>>>> >>>> >>>> On 2014-05-15 10:30, tarak karmakar wrote: >>>>> >>>>> Dear Sir, >>>>>> Thanks for the quick reply. >>>>>> Using g_anaeig_d I've got both of them as 'nan' >>>>>> >>>>>> The Entropy due to the Quasi Harmonic approximation is -nan J/mol K >>>>>> The Entropy due to the Schlitter formula is nan J/mol K >>>>>> >>>>>> >>>>> Interesting. Is your system far out of equilibrium? Do the eigenvalues >>>>> look reasonable? Or is your system very large? >>>>> >>>>> >>>>> Tarak >>>>>> >>>>>> >>>>>> On Thu, May 15, 2014 at 1:30 PM, David van der Spoel >>>>>> <sp...@xray.bmc.uu.se>wrote: >>>>>> >>>>>> On 2014-05-15 09:07, tarak karmakar wrote: >>>>>> >>>>>>> >>>>>>> Dear All, >>>>>>> >>>>>>>> I wanted to calculate configurational entropy of a protein by using >>>>>>>> g_covar >>>>>>>> and g_anaeig as follows, >>>>>>>> >>>>>>>> g_covar -f ../traj.xtc -s ../npt_prod -o eigenval -v eigenvec.trr >>>>>>>> -av >>>>>>>> average.pdb >>>>>>>> g_anaeig -v eigenvec.trr -entropy -temp 300 >>>>>>>> >>>>>>>> I got the following results >>>>>>>> >>>>>>>> The Entropy due to the Quasi Harmonic approximation is 31440.3 >>>>>>>> J/mol K >>>>>>>> The Entropy due to the Schlitter formula is nan J/mol K >>>>>>>> >>>>>>>> I went back to check the eigenvector file by dumping it to a new >>>>>>>> file >>>>>>>> 'dump'. There is no 'nan' indeed. >>>>>>>> >>>>>>>> Could anyone comment on why I'm getting Schlitter entropy as 'nan'? >>>>>>>> >>>>>>>> Thanks and regards, >>>>>>>> Tarak >>>>>>>> >>>>>>>> Numerical error due to large numbers. Try compiling gromacs in >>>>>>>> double >>>>>>>> >>>>>>>> precision. >>>>>>> >>>>>>> -- >>>>>>> David van der Spoel, Ph.D., Professor of Biology >>>>>>> Dept. of Cell & Molec. Biol., Uppsala University. >>>>>>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. >>>>>>> sp...@xray.bmc.uu.se http://folding.bmc.uu.se >>>>>>> -- >>>>>>> Gromacs Users mailing list >>>>>>> >>>>>>> * Please search the archive at http://www.gromacs.org/ >>>>>>> Support/Mailing_Lists/GMX-Users_List before posting! >>>>>>> >>>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >>>>>>> >>>>>>> * For (un)subscribe requests visit >>>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-usersor >>>>>>> send a mail to gmx-users-requ...@gromacs.org. >>>>>>> >>>>>>> >>>>>>> >>>>> -- >>>>> David van der Spoel, Ph.D., Professor of Biology >>>>> Dept. of Cell & Molec. Biol., Uppsala University. >>>>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. >>>>> sp...@xray.bmc.uu.se http://folding.bmc.uu.se >>>>> -- >>>>> Gromacs Users mailing list >>>>> >>>>> * Please search the archive at http://www.gromacs.org/ >>>>> Support/Mailing_Lists/GMX-Users_List before posting! >>>>> >>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >>>>> >>>>> * For (un)subscribe requests visit >>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >>>>> send a mail to gmx-users-requ...@gromacs.org. >>>>> >>>>> >>>> >>>> >>> >>> >> >> -- >> David van der Spoel, Ph.D., Professor of Biology >> Dept. of Cell & Molec. Biol., Uppsala University. >> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. >> sp...@xray.bmc.uu.se http://folding.bmc.uu.se >> -- >> Gromacs Users mailing list >> >> * Please search the archive at http://www.gromacs.org/ >> Support/Mailing_Lists/GMX-Users_List before posting! >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> * For (un)subscribe requests visit >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> send a mail to gmx-users-requ...@gromacs.org. >> > > > > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.