# [gmx-users] Schlitter Configurational entropy 'NaN'

tarak karmakar tarak20489 at gmail.com
Thu May 15 12:02:28 CEST 2014

```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 <tarak20489 at 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 <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.
>>
>
>
>
>
```