Previous Thread
Next Thread
Print Thread
Joined: Dec 2011
Posts: 6
P
Forum Member
OP Offline
Forum Member
P
Joined: Dec 2011
Posts: 6
Dear all,
I calculated several 1 microsecond MD simulations for an homodimer in a water box. The structure of the Protein is stable along the trajectories excepted its N-terminal region (that we previously demonstrated by NMR relaxation).

I now try to calculate the configurational entropy of each monomer from these simulations. For this, i first created two trajectories to remove rotation and translation and water molecules) for each monomer (one such trajectory for monomer A and one for monomer B). Then i calculate the configurationnal entropy of each monomer from the covariance matrix (using the coor cova entropy command, andricioaei & Karplus J. Chem Phys. 2001, 115, 6289).

The exact command i used (for monomer A) is:

coor cova firstu 24 nunit 1 begin @startf skip @skips stop @stope -
sele segid @sega end sele segid @sega end -
entropy unit -1 temp 293 resi

I wrote a script to calculate the configurationnal entropy for increased length of simulation, but apparently the total entropy does not converged (results below).

##########################################################
# STOPE(STEPS) TIMEMAX(NS) SA SB
##########################################################
20000000 40.00 4.54206733 8.67930683
30000000 60.00 4.60703434 9.37740553
40000000 80.00 8.01340579 5.40335717
50000000 100.00 4.51616798 5.27548147
60000000 120.00 4.65081627 10.1002701
70000000 140.00 8.45059191 5.57077596
80000000 160.00 8.75707062 5.62990561
90000000 180.00 5.07571439 5.56165663
100000000 200.00 5.01090967 10.4110834
110000000 220.00 9.02396883 5.61459658
120000000 240.00 5.0635709 10.6797375
130000000 260.00 9.19480905 10.7812281
140000000 280.00 5.15129031 10.8492845
150000000 300.00 9.40821766 10.8243819
160000000 320.00 5.1931215 10.7860315
170000000 340.00 5.22793727 10.7364907
180000000 360.00 9.46028759 5.70225596
190000000 380.00 5.26398717 10.6442232
200000000 400.00 5.3682301 5.67329256
210000000 420.00 9.5172872 10.5949072
220000000 440.00 9.50143223 10.5784023
230000000 460.00 5.21329563 10.588253
240000000 480.00 9.53142293 10.5829299
250000000 500.00 9.51825395 10.5620274
260000000 520.00 9.51029262 5.74673113
270000000 540.00 5.21883109 5.81399319
280000000 560.00 9.5528929 5.90840359
290000000 580.00 9.58950676 10.5252618
300000000 600.00 9.60300341 5.82241847
310000000 620.00 5.22301875 10.5005907
320000000 640.00 5.33463332 10.5064215
330000000 660.00 5.32644774 5.75521376
340000000 680.00 9.62913189 5.87716752
350000000 700.00 9.63007295 5.69807031
360000000 720.00 9.62628892 10.5181028
370000000 740.00 9.61467287 5.8846417
380000000 760.00 9.60119472 10.5287241
390000000 780.00 5.22981812 10.5192489
400000000 800.00 9.5703358 5.87404549
410000000 820.00 5.27369269 5.8547552
420000000 840.00 9.53539689 10.5070513
430000000 860.00 5.33562448 10.5034586
440000000 880.00 9.50372638 10.5079842
450000000 900.00 9.48696289 10.5191612
460000000 920.00 5.21505484 10.5341277
470000000 940.00 5.262273 10.5288986
480000000 960.00 5.20083735 10.5473411
490000000 980.00 5.22263514 5.87958278
500000000 ****** 5.21326949 6.03540297


I wonder if someone tried something like this and what can explain this observation. Alternatively a suggestion to calculate configurationnal entropy from MD trajectories would be appreciated.
Best regards.
Philippe.

Joined: Sep 2003
Posts: 4,850
Likes: 6
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 4,850
Likes: 6
Configuration entropy converges relatively slowly. Are the entropy values on the last line of your table calculated for 1microsecond, or just for the 20 ns from 980 to 1000ns?
There is something strange about your data. Each entropy column seems to contain two distinct data sets, one which converges to 9.5 and one to 5.3 (for SA).


I use the QUASiharmonic command within the VIBRAN molecule, roughly:
read rtf card name top_all36_prot.rtf (open and close are not necessary)
read para card name par_all36_prot.rtf
read psf card name my_prot.psf
open unit @fu read file name my.trj
set T 300
coor dyna firstu @fu nunit @nu sele @solute end

calc nmodes = 3 * ?NSEL
bomlev -2
vibran nmode @nmodes
open unit 11 write form name @s/quasi-@j.txt
outu 11
echu 11

ECHO protein
quasi firstu @fu nunit @nu begin @begin stop @stop -
sele @prot end -
thermo resi temp @T

outu 6

end


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Joined: Dec 2011
Posts: 6
P
Forum Member
OP Offline
Forum Member
P
Joined: Dec 2011
Posts: 6
Thank you for your answer.
Yes it is exactly why i am surprised. It seems that, for each monomer, there are two S values depending on the length of the simulation and when plotting these values as a function of the length of the simulation, each of these two values seems to converge to a plateau differing by about a factor 2 (see figure attached).
The values of S are calculated for the total range meaning for instance for a stopmax of 1000ns, the set of frames are extracted between 20ns to 1000ns each 20ps. So i was expecting that S converge to a single value.
Since my first post, i tried to build a traj for one monomer and to work with it but i found the same 'two values" behaviour.
So i will stop with this "coor cova" approach and try to do it with quasi.

Thank you again.
best regards.
Philippe.

Attached Images
Attached PDF document
fig_entrop.pdf (21.76 KB, 42 downloads)

Moderated by  BRBrooks, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.008s Queries: 21 (0.004s) Memory: 0.7422 MB (Peak: 0.8078 MB) Data Comp: Off Server Time: 2022-09-27 01:42:49 UTC
Valid HTML 5 and Valid CSS