[gmx-users] Schlitter Configurational entropy 'NaN'
Justin Lemkul
jalemkul at vt.edu
Thu May 15 18:02:31 CEST 2014
On 5/15/14, 11:51 AM, tarak karmakar wrote:
> Part of the eigenval.xvg file
>
> 995 3.02979e-06
> 996 2.98051e-06
> 997 2.93131e-06
> 998 2.75594e-06
> 999 2.69955e-06
> 1000 2.56434e-06
> 1001 2.66401e-16
> 1002 1.43252e-16
> 1003 1.0735e-16
> 1004 8.0728e-17
> 1005 4.3465e-17
> 1006 1.60571e-17
> 1007 1.49006e-17
> 1008 1.44694e-17
> 1009 1.07083e-17
> 1010 9.77564e-18
> 1011 8.25739e-18
> 1012 8.18473e-18
> 1013 6.85267e-18
> 1014 6.79135e-18
> 1015 5.7498e-18
> 1016 5.49917e-18
> 1017 4.84525e-18
> 1018 4.15989e-18
> 1019 4.06407e-18
> 1020 3.70677e-18
> .............................
> .............................
> 1362 -3.16429e-18
> 1363 -3.40673e-18
> 1364 -3.85766e-18
> 1365 -4.16796e-18
> 1366 -4.16839e-18
> 1367 -4.17689e-18
> 1368 -4.48977e-18
> 1369 -4.54731e-18
> 1370 -4.58794e-18
> 1371 -4.86076e-18
> 1372 -5.05864e-18
> 1373 -5.73908e-18
> 1374 -6.1385e-18
> 1375 -6.98907e-18
> 1376 -7.5897e-18
> 1377 -1.23406e-17
> 1378 -1.24701e-17
> 1379 -1.41182e-17
> 1380 -1.47709e-17
> 1381 -1.50144e-17
> 1382 -5.73596e-17
> 1383 -1.03466e-16
> 1384 -1.277e-16
> 1385 -1.64873e-16
> 1386 -3.2657e-16
>
> These eigenvalues are very less compared to the below 1000 numbers.
> I see there is no strange jump from 999 to 1000 but from 1000 to 1001.
>
The array in the code is zero-indexed, so this is actually the location of the
problem.
-Justin
> Lambda values are considered to be zero in these cases and thus
> w = sqrt(BOLTZMANN*temp/lambda)/NANO = infinity
> &
> dd = 1+kteh*eigval[i];
> deter += log(dd)
> is failed.
>
> 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
>
> So, there is precision problem in the code; 10^-16 is taken as zero!!
>
> Tarak
>
>
>
>
> On Thu, May 15, 2014 at 5:23 PM, David van der Spoel
> <spoel at xray.bmc.uu.se>wrote:
>
>> On 2014-05-15 11:40, tarak karmakar 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)?
>>>
>> Check your eigval.xvg file, what is the value at i = 1000?
>> It seems that lambda is the problem, hence I suspect the eigval[1000] to
>> be strange.
>>
>>
>>> Tarak
>>>
>>>
>>>
>>> On Thu, May 15, 2014 at 2:54 PM, David van der Spoel
>>> <spoel at 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 <tarak20489 at 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 <
>>>>>> spoel at 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
>>>>>>>> <spoel at 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.
>>>>>>>>> spoel at 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-request at 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.
>>>>>>> spoel at 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-request at 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.
>>>> spoel at 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-request at 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.
>> spoel at 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-request at gromacs.org.
>>
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list