Previous Thread
Next Thread
Print Thread
#19449 12/01/08 12:18 PM
Joined: Sep 2004
Posts: 57
B
Baloo Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Sep 2004
Posts: 57
Hi,

I am looking for the contact between two atoms of same molecule and using coor contact command. The output gives me a table in the following format..

Quote:


I-atom J-atom <occupancy> <time> # events






under <time> I get the average time of that particular contact upto one decimal point. Is there any way, I can get more digits after decimal.

Best
-Baloo

Baloo #19450 12/01/08 12:37 PM
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
In source/manip/hbanal.src change the following line:
908 FORMAT(5X,4(1X,A),' - ',4(1X,A),F7.3,4X,F7.1,I8)

The "1" in "F7.1" indicates the number of decimals you want.

Then you install CHARMM again using the same install.com commmand as when CHARMM was installed last time; do not remove any files in the installation tree.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #19451 12/01/08 10:29 PM
Joined: Sep 2004
Posts: 57
B
Baloo Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Sep 2004
Posts: 57
Dear Prof. Lennart,

Thanks, it worked. Can I also increase the digits after decimal for occupancy as well.

Best
-Baloo

Baloo #19452 12/02/08 06:55 AM
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
In principle yes, but since you ask perhaps not. You should perhaps prepare careful arguments if a future referee asks about the statistical significance in occupanices reported with more then three decimal places...


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #19453 12/02/08 09:42 AM
Joined: Sep 2004
Posts: 57
B
Baloo Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Sep 2004
Posts: 57
The reason for asking is that in some cases I get <occupancy> of 0.000. Lets look at the following example.



Quote:


<occupancy> <time> # events
0.000 0.2341 88
0.000 0.3064 47





Both have <occupancy> of 0.000 but if I go upto 4 decimal places then I would be having 0.0004 and 0.0003. I have MD of 49ns. Someone might say that there is only the difference of 0.0001 in <occupancy> but I believe to see some difference even at 4th decimal place is also not bad than having <occupancy> difference of 0.000 with three decimal places.

Best
-Baloo

Baloo #19454 12/02/08 01:41 PM
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
Seeing a difference that is just a spurious result due to insufficient data is not the best thing either. This is your decision to make, but you should use a sound statistical argument. If you repeat the analysis over 5 blocks of 9ns, how much do these numbers vary?


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #19455 12/03/08 11:27 AM
Joined: Sep 2004
Posts: 57
B
Baloo Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Sep 2004
Posts: 57
Dear Prof. Lennart,

Thanks for clarification. However, on carefull inspection of the coor contact outputs. I found one discrepancy. Sum of the number of events of contacts in the output file is different from the number of lines (which should reflect the total number of events) in the IUNIt file, showing more detailed picture of contacts. Even the output summary for any pair indicating number of events of that particular contact is different than the number of lines (which should reflect the number of events of that particular atom-pair) in the IUNIt file.

Any clue, why so?

Best
-Baloo

Baloo #19456 12/03/08 03:11 PM
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
That is strange. There is nothing in the code immediately suggesting how such differences could occur, and in a quick test the various counts match for me.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #19457 12/03/08 09:04 PM
Joined: Sep 2004
Posts: 57
B
Baloo Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Sep 2004
Posts: 57
Things are becoming clear but still some mysteries are unsolved. I believe that it is my mistake that I didn't make it clear that I am trying to figure out the hydrogen bridges using the command..

Quote:


set CUTHB 4.0 ! hydrogen bond cutoff

define OXY SELECT (TYPE O*) .and. .not. SEGID WAT END !solute Oxygens

coor contact VERBOSE Bridge TIP3 CUTHB @CUTHB TCUT 0 IUNIt 21 select OXY end select OXY end FIRST 20 NUNIT 1 NSKIP ?skip BEGIN ?start STOP ?nstep






For the sake of clarity, I considered only the 500ps of MD with much longer hydrogen bond cutoff (I am using 50 Ang. in this test case). The output summary looks like..

xxxxxxxxxxxxxxxxxxx--Output_summary

2500 CORD RECORDS READ FROM 1 UNITS STARTING WITH UNIT 20
RUNNING FROM STEP24750100 TO 25000000 SKIPPING 100 STEPS BETWEEN RECORDS
Time step was 4.0909659E-02 AKMA time units.

Analysis of TIP3 close contact bridges, using cutoff distance= 50.00 and angle= 999.00
Total time= 500.0ps. Resolution= 0.20ps
Occupancy cutoff= 0.00 Lifetime cutoff= 0.0ps

I-atom J-atom <occupancy> <time> # events
---------------------------------------------------------------------------
LK13 LK13 0 O1 - LK13 LK13 0 O2 19.064 18.8008 507
LK13 LK13 0 O1 - LK13 LK13 0 O3 17.801 17.5550 507
LK13 LK13 0 O1 - LK13 LK13 0 O4 20.105 19.8276 507
LK13 LK13 0 O1 - LK13 LK13 0 O5 20.650 20.3653 507
LK13 LK13 0 O1 - LK13 LK13 0 O6 21.656 21.3566 507
LK13 LK13 0 O1 Summary: 99.276 19.6 2535
LK13 LK13 0 O2 - LK13 LK13 0 O1 19.064 18.8008 507
LK13 LK13 0 O2 - LK13 LK13 0 O3 21.025 20.7349 507
LK13 LK13 0 O2 - LK13 LK13 0 O4 21.825 21.5239 507
LK13 LK13 0 O2 - LK13 LK13 0 O5 21.888 21.5862 507
LK13 LK13 0 O2 - LK13 LK13 0 O6 22.781 22.4667 507
LK13 LK13 0 O2 Summary:106.584 21.0 2535
LK13 LK13 0 O3 - LK13 LK13 0 O1 17.801 17.5550 507
LK13 LK13 0 O3 - LK13 LK13 0 O2 21.025 20.7349 507
LK13 LK13 0 O3 - LK13 LK13 0 O4 22.172 21.8659 507
LK13 LK13 0 O3 - LK13 LK13 0 O5 21.850 21.5479 507
LK13 LK13 0 O3 - LK13 LK13 0 O6 22.505 22.1941 507
LK13 LK13 0 O3 Summary:105.352 20.8 2535
LK13 LK13 0 O4 - LK13 LK13 0 O1 20.105 19.8276 507
LK13 LK13 0 O4 - LK13 LK13 0 O2 21.825 21.5239 507
LK13 LK13 0 O4 - LK13 LK13 0 O3 22.172 21.8659 507
LK13 LK13 0 O4 - LK13 LK13 0 O5 25.494 25.1424 507
LK13 LK13 0 O4 - LK13 LK13 0 O6 25.882 25.5243 507
LK13 LK13 0 O4 Summary:115.478 22.8 2535
LK13 LK13 0 O5 - LK13 LK13 0 O1 20.650 20.3653 507
LK13 LK13 0 O5 - LK13 LK13 0 O2 21.888 21.5862 507
LK13 LK13 0 O5 - LK13 LK13 0 O3 21.850 21.5479 507
LK13 LK13 0 O5 - LK13 LK13 0 O4 25.494 25.1424 507
LK13 LK13 0 O5 - LK13 LK13 0 O6 0.805 0.7921 508
LK13 LK13 0 O5 Summary: 90.688 17.9 2536
LK13 LK13 0 O6 - LK13 LK13 0 O1 21.656 21.3566 507
LK13 LK13 0 O6 - LK13 LK13 0 O2 22.781 22.4667 507
LK13 LK13 0 O6 - LK13 LK13 0 O3 22.505 22.1941 507
LK13 LK13 0 O6 - LK13 LK13 0 O4 25.882 25.5243 507
LK13 LK13 0 O6 - LK13 LK13 0 O5 0.805 0.7921 508
LK13 LK13 0 O6 Summary: 93.628 18.5 2536
-----------------------------------------------------
OVERALL
Average number (over ISEL!)=101.834
<lifetime>= 20.09ps
Total number/frame= 611.0

xxxxxxxxxxxxxxxxxxx

The output file with detailed hydrogen bridge statistics looks like..

xxxxxxxxxxxxxxxxxxx--Output_detailed

I-atom J-atom (Bridge) Lifetime Endtime
LK13 0 LK13 O1 - LK13 0 LK13 O2 WAT 1 6641775.61 500.00
LK13 0 LK13 O1 - LK13 0 LK13 O3 WAT 1 6641144.01 500.00
LK13 0 LK13 O1 - LK13 0 LK13 O4 WAT 1 6642296.21 500.00
LK13 0 LK13 O1 - LK13 0 LK13 O5 WAT 1 6642568.81 500.00
LK13 0 LK13 O1 - LK13 0 LK13 O6 WAT 1 6643071.41 500.00
LK13 0 LK13 O2 - LK13 0 LK13 O1 WAT 1 6641775.61 500.00
LK13 0 LK13 O2 - LK13 0 LK13 O3 WAT 1 6642756.21 500.00
LK13 0 LK13 O2 - LK13 0 LK13 O4 WAT 1 6643156.21 500.00
LK13 0 LK13 O2 - LK13 0 LK13 O5 WAT 1 6643187.81 500.00
LK13 0 LK13 O2 - LK13 0 LK13 O6 WAT 1 6643634.21 500.00
LK13 0 LK13 O3 - LK13 0 LK13 O1 WAT 1 6641144.01 500.00
LK13 0 LK13 O3 - LK13 0 LK13 O2 WAT 1 6642756.21 500.00
LK13 0 LK13 O3 - LK13 0 LK13 O4 WAT 1 6643329.61 500.00
LK13 0 LK13 O3 - LK13 0 LK13 O5 WAT 1 6643168.41 500.00
LK13 0 LK13 O3 - LK13 0 LK13 O6 WAT 1 6643496.01 500.00
LK13 0 LK13 O4 - LK13 0 LK13 O1 WAT 1 6642296.21 500.00
LK13 0 LK13 O4 - LK13 0 LK13 O2 WAT 1 6643156.21 500.00
LK13 0 LK13 O4 - LK13 0 LK13 O3 WAT 1 6643329.61 500.00
LK13 0 LK13 O4 - LK13 0 LK13 O5 WAT 1 6644990.81 500.00
LK13 0 LK13 O4 - LK13 0 LK13 O6 WAT 1 6645184.41 500.00
LK13 0 LK13 O5 - LK13 0 LK13 O1 WAT 1 6642568.81 500.00
LK13 0 LK13 O5 - LK13 0 LK13 O2 WAT 1 6643187.81 500.00
LK13 0 LK13 O5 - LK13 0 LK13 O3 WAT 1 6643168.41 500.00
LK13 0 LK13 O5 - LK13 0 LK13 O4 WAT 1 6644990.81 500.00
LK13 0 LK13 O5 - LK13 0 LK13 O6 WAT 1 6645753.21 500.00
LK13 0 LK13 O6 - LK13 0 LK13 O1 WAT 1 6643071.41 500.00
LK13 0 LK13 O6 - LK13 0 LK13 O2 WAT 1 6643634.21 500.00
LK13 0 LK13 O6 - LK13 0 LK13 O3 WAT 1 6643496.01 500.00
LK13 0 LK13 O6 - LK13 0 LK13 O4 WAT 1 6645184.41 500.00
LK13 0 LK13 O6 - LK13 0 LK13 O5 WAT 1 6645753.21 500.00

xxxxxxxxxxxxxxxxxxx


The reason of using 50 Ang. of hydrogen bond cutoff was that I can see the <occupancy> of 1 in all the atom pairs connected via hydrogen bridge.

Some doubts..

1- In Output_summary, <occupancy> is some kind of percentage and it should be in the range of 0-1. But it has values like 19.064, 17.801, ...

2- Since, very high hydrogen bond cutoff has been given, I thought that every atom pair of the solute must be connected via hydrogen bridge and there should be only 1 event from start to the end irrespective of whether 500ps or 1ns or entire MD trajectory is taken into consideration. But in Output_summary, I get 507 events for every pair.

3- In Output_summary, it is difficult to understand that why every pair has different <occupancy> and <time>, even with high hydrogen bond cutoff. Actually <occupancy> should be identical and same should be true of <time> and events.

4- In Output_detailed, under lifetime column what these numbers like 6641775.61, 6641144.01, ... indicate. I thought the unit of lifetime is 'ps', but the units of these numbers can not be 'ps', as I don't have such long MD simulation.

5- It would be very helpful if some one can comment on the units of lifetime, endtime & <time>. Whether they are SI or AKMA. Sorry to say but units are not very clear to me.

Best
-Baloo

Baloo #19458 12/03/08 09:36 PM
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
Yes, trying to describe what you did in your own words instead of showing actual input/output is usually a very big mistake (this is an understatement).

With bridges things get much more complicated since there can be multiple paths through single bridging molecules, and the combination of CONTact and BRIDge is potentially even more complicated; the results have to be analyzed keeping this in mind. You have managed to pick a particularly unintended set of parameters (very frequent sampling, contact, bridge, very long cutoff). Note that with the CONTact variant it may be incorrect to speak about hydrogen bonds/bridges. There is a built-in limitation in the "event counter" in that the accumulated lifetime should not exceed 65000 analyzed frames, a limit that should normally not be exceeded even today, but your "test case" may violate this, and then all bets are off as to what numbers you will see in the output files.

As to your questions (you probably do not mean "doubts" - a doubt is when you think something is not true):
Times are reported in ps. Distances are in A (not Ang.) (usage.doc)

If a given atom throughout the simulation is within the cutoff distance for CONTact/HBONd with two other atoms the occupancy will be 2.0. Is this strange? Should percentages be restricted to < 100%?


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Moderated by  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~deb10u2 Page Time: 0.014s Queries: 34 (0.010s) Memory: 0.7865 MB (Peak: 0.8754 MB) Data Comp: Off Server Time: 2023-01-27 08:54:24 UTC
Valid HTML 5 and Valid CSS