[gmx-users] Schlitter Configurational entropy 'NaN'
tarak karmakar
tarak20489 at gmail.com
Thu May 15 11:40:18 CEST 2014
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-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.
>
--
*Tarak KarmakarMolecular Simulation Lab.Chemistry and Physics of Materials
UnitJawaharlal Nehru Centre for Advanced Scientific Research Jakkur P. O.
Bangalore - 560 064Karnataka, INDIAPh. (lab) : +91-80-22082809 *
More information about the gromacs.org_gmx-users
mailing list