I am unable to calculate the entropy values either schlitter or andricioaei

[vijay kumar narsapuram]

My system is 10ns simulated production trajectory with bound ligand.

The Andricioaei entropy can not be calculated due to
the presence of an eigenvalue less than or equal to zero
Entropy (Schlitter)   is -nan (J/molK)
done.
I have tried following method to improve the results, but no use..
1)increase in step size.
2)removel of flexible regions.

please suggest some other changes if available.

and i wanted to know the significance of entropy values which we get(schlitter or  andricioaei).which one should i consider for my analysis..?

Thank you

Regards
vijay kumar N

[Nicholas M Glykos]

> Entropy (Schlitter) is -nan (J/molK)

The simplest explanation is that you have an eigenvalue < -1.0 (check the
contents of the file containing the eigenvalues as produced by carma).

Have you removed overall rotations/translations before calculating entropy ?

After excluding ther flexible regions, have you re-removed overall
rotations/translations ?

> get(schlitter or andricioaei).which one should i consider for my
> analysis..?

Schlitter's is more stably calculatable.




--


Nicholas M. Glykos, Department of Molecular Biology
and Genetics, Democritus University of Thrace, University Campus,
Dragana, 68100 Alexandroupolis, Greece, Tel/Fax (office) +302551030620,
Ext.77620, Tel (lab) +302551030615, http://utopia.duth.gr/~glykos/

[vijay kumar narsapuram]

The least value observe to be -0.16927 for eigenvalue

and i have performed removal of rotations/translations after every selection of residues. 

[Nikolaos Glykos]

> The least value observe to be -0.16927 for eigenvalue

Have you excluded the lowest six eigenvalues ? (corresponding to
overall rotations/translations) :

The basic iterations for Schlitter's estimate is

entropy = 0.0;
for ( i=1 ; i <= eigendim -6 ; i++ )
entropy += log1p( CONST2 * Temp * M_E * M_E * eigenvalues[i] );


where

#define CONST2 (double)(0.020614869) /* (AMU*1e-20) /
(HBAR^2 / Kb) */

end M_E is exp(1).


This means that with an eigenvalue of -0.16927, any temperature higher
than
~39 K will result to a NaN.


Now, well-behaved (stable) structures have diminishingly small
eigenvalues near
the end of the table. For example, the analysis of a ~360 residues
trimer gave
these eigenvalues

0.0010729906
0.0010554525
0.0000028583 \
-0.0000126419 |
-0.0000164017 |== Last 6, corresponding to overall rot/trans
-0.0000191536 |
-0.0000244375 |
-0.0000265611 /


If you observe such large eigenvalues either the protein is extremely
flexible (disordered ?) or you have a 'bad' frame in the DCD or any
other
reason that invalidates the harmonic approximation.

If you have reasons to believe that the protein is well-behaved and that
this may be a bug in carma, I'll be happy to look into it, but I'll be
needing
access to your DCD-PSF files.

[Nicholas M Glykos]

> and negative eigenvalue this time started from different frame 12667 but
> not as previously mentioned 13000 frame. can u please specify any
> particular script to remove negative eigen values..?

This doesn't make sense. Eigenvalues are _not_ a function of specific
frames.

How big is the protein ? How many atoms carma reports when it starts ?

[vijay kumar narsapuram]

No of atoms in protein =5288

[Nicholas M Glykos]

> No of atoms in protein =5288

This may also contribute to your problems. Maybe 10 ns worth of simulation
are not sufficient to meaningfully populate the variance-covariance
matrix, especially if the protein is highly flexible. I find it unlikely
that there is a your-protein-specific bug in carma, but you can email me
(directly me, not the list) your carma.XXXX.eigenvalues.dat and
carma.XXXX.eigenvectors.dat files to have a look at them.

[Nicholas M Glykos]

I've seen your files. The negative eigenvalues start very early: already
for the eigenvector 12664 (out of a total of 15864 eigenvectors) they
become negative. I suspect that a possible reason can be seen in the very
large value of the first eigenvalue, indicating the presence of large
amplitude motion. This is also evident from the graph of atomic
fluctuations per atom per eigenvector which show large mobilities for the
first ~350 atoms, and then again from atom ~4140 to the end (especially
after atom 4910). If the problem is that there is indeed large scale
protein motion present in your simulation, then the assumptions of the
quasi-harmonic approximation break down, you can not meaningfully remove
overall rotations/translations, and you must think carefully about whether
it makes structural sense to sub-divide the calculations (maybe on a
per-domain basis ?).

My twocents.

[Nicholas M Glykos]

> But neglecting residues doesn't add any significance to my
> calculations.since, I am calculating other energy values by considering
> all these atoms without neglecting any.
>
> can you please suggest any other alternative for calculating solute
> entropy..?


OK, you asked for it, you got it : from

http://utopia.duth.gr/~glykos/progs/

download the linux executable named 'carma64', chmod it, and install it.
Then use the new executable as usual, _but_ add the flag '-force' in the
command line.

The '-force' flag tells carma to stop the calculation of entropies before
they become undefined.

Carma will emit cautionary messages in big bold red letters (and with
exclamation marks !) so that you know that you'd better know what you are
doing. It will also tell you how many eigenvalues actually used for the
calculation. I have tested the program and it works with Andricioaei's
formula, but I do not have a system for which Schlitter's formula goes to
Nan, so test it and tell me how it goes.

[vijay kumar narsapuram]

ERROR:Unidentified keyword or argument : -force. Abort.

when i used following command i got above error.

./carma64 -v -fit -force -atmid ALLID newdcd-1.dcd newpsf.psf > carma_entropy_after_newdcd-1.out

./carma64 -v -cov -eigen -mass -force -temp 300 -atmid ALLID carma.fitted.dcd newpsf.psf > carma_entropy_after_newdcd-2.out

[Nicholas M Glykos]

> ERROR:Unidentified keyword or argument : -force. Abort.

For some reason you are still using the old executable (have you
downloaded the carma archive instead of the carma64 executable ?).
Try this from the unix shell :

wget http://utopia.duth.gr/~glykos/progs/carma64
chmod 755 carma64

./carma64 -v -cov -eigen -mass -force -temp 300 -atmid ALLID carma.fitted.dcd newpsf.psf > carma_entropy_after_newdcd-2.out

[vijay kumar narsapuram]

After using new executable carma64.
This was the output..

Entropy calculation will ignore negative eigenvalues !

Entropy (Andricioaei) using only  11278 eigenvalues is 57985.630383 (J/molK)
Entropy (Schlitter)   using only  13122 eigenvalues is 59801.236248 (J/molK)
done.

[Nicholas M Glykos]

> After using new executable carma64.
> This was the output..
>
> Entropy calculation will ignore negative eigenvalues !
>
> Entropy (Andricioaei) using only 11278 eigenvalues is 57985.630383 (J/molK)
> Entropy (Schlitter) using only 13122 eigenvalues is 59801.236248 (J/molK)
> done.

OK, looks sensible, ie. (a) more eigenvalues are used with Schlitter's
estimate, and, (b) the two estimates obtained are of comparable magnitude
although Andricioaei's uses almost 2000 eigenvalues fewer than
Schlitter's.

Having said that, I still have my doubts whether this is a reasonable
solution to your problem : If it were for just a couple a slightly
negative eigenvalues near the end of the eigenvalues table, I would not
worry too much. But when using only 13122 out of 15864 eigenvalues, I do
not have a clue what is the magnitude of error introduced in the
calculation. Maybe the quasi-harmonic approximation is not valid for your
problem and a completely different approach is needed ?

Good luck with it,
Nicholas