[gmx-users] Schlitter Configurational entropy 'NaN'
David van der Spoel
spoel at xray.bmc.uu.se
Thu May 15 19:20:33 CEST 2014
On 2014-05-15 18:41, tarak karmakar wrote:
> Hi Justin,
> That's true and it is a C code and 0 index loop.
> Then the strange jump 999 to 1000 came in to the picture.
> But what about dealing with very small numbers?
What does g_covar print?
It seems the problem might be there, because the jump by 10 orders of
magnitude in eigenvalue does not seem reasonable.
>
>
> Tarak
>
>
> On Thu, May 15, 2014 at 9:30 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> 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-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.
>>>>
>>>>
>> --
>> ==================================================
>>
>> 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
>>
>> ==================================================
>>
>> --
>> 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
More information about the gromacs.org_gmx-users
mailing list