|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
*FILENAME: hbonds.inp *PURPOSE: compute static and time-dependent hydrogen bonds from trajectory *AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003 * !unix environment variable CHM_HOME points to CHARMM installation directory read rtf card name $CHM_HOME/toppar/top_all22_prot.inp read para card name $CHM_HOME/toppar/par_all22_prot.inp ! we assume there are three segments in this psf ! prot = a protein ! lig = a bound ligand ! wat = solvent waer read psf card name my_psf.psf read coor card name my_structure.crd ! hydrogen bonding patterns in initial structure ! the "coor hbond" depends on DONOr/ACCEptor information in PSF to ! identify possible hbonding partners; see corman.doc for details ! we use strictly geometric criteria on the H...A distance, and optionally ! on the D-H...A angle to define the hydrogen bond ! first look at the backbone-backbone hydrogen bonds coor hbond verbose select type HN .or. type O end select type HN .or. type O end ! protein-solvent !!NB your psf MUST have DONOr/ACCEptor specifications ! for the water molecules for this to work; ! if any of the following lines are missing from RESI TIP3 ! in your RTF (and hence in the PSF), add missing bits to RTF and ! regenerate the PSF ! ACCEPTOR OH2 ! DONOR H1 OH2 ! DONOR H2 OH2 coor hbond verbose select segid prot end select segid wat end
! now get the timedependent hydrogen bonds; the output can be voluminous, especially ! with the verbose option open unit 51 read unform name my_traj1.cor open unit 52 read unform name my_traj2.cor ! the SKIP 2500 (=5ps) means that for a hydrogen bond to be counted as broken ! in the calculation of average lifetime it has to be broken for > 5ps (otherwise ! we are unlikely to detect the break)
! backbone-backbone coor hbond verbose select type HN .or. type O end select type HN .or. type O end - firstu 51 nunit 2 begin 10500 skip 2500
! note that the trajectory is not closed between these commands, and can be reused ! protein - solvent coor hbond verbose select segid prot end select segid wat end - firstu 51 nunit 2 begin 10500 skip 2500
! protein - water - ligand coor hbond verbose select segid prot end select segid lig end - bridge tip3 firstu 51 nunit 2 begin 10500 skip 2500
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: May 2004
Posts: 222
Forum Member
|
Forum Member
Joined: May 2004
Posts: 222 |
Hi. Which column in the psf specifies donor/acceptor status? Thanks.
CHARMM 30b1 driven by
1/ Xeon (32 bits)
2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel
3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum.
4/ Gracious support from the forum.
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
It's not a column, but 2 sections, a donor list and an acceptor list; the list is a sequence of integers, which is the atom number index defined in the first section of the PSF with the charges, etc.
Write out a formatted PSF file and have a look for yourself ...
Rick Venable computational chemist
|
|
|
|
Joined: May 2004
Posts: 222
Forum Member
|
Forum Member
Joined: May 2004
Posts: 222 |
Just to make sure the explanation made is understood, i 'll use an example.
so in the example below, one always read along a row with a donor index followed by an acceptor index, donor then acceptor kind of fashion.
so atom 1 and 5, 2 and 1, 3 and 1 are examples of donor/acceptor pairs?
30279 !NBOND: bonds 1 5 2 1 3 1 4 1 5 6 7 5 7 8 7 9 10 7 10 11 10 12 13 10 13 14 13 15 16 13 16 17 16 18 19 16 19 20 19 21 19 22 23 5 23 25 24 23 25 26 25 27 27 28 29 27 29 30 31 29 31 32 31 33 31 34 35 29 35 36 35 37 35 38 39 27 39 41 40 39 41 42 41 43 43 44 45 43 45 46 45 47 48 45 49 48 49 50 51 49 51 52 53 51
CHARMM 30b1 driven by
1/ Xeon (32 bits)
2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel
3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum.
4/ Gracious support from the forum.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
NO, this is not correct (and why are you looking at the list of bonds?).
Add DONOr and ACCEptor statements as needed to your RTF and regenerate the PSF - do not play directly with the PSF.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: May 2004
Posts: 222
Forum Member
|
Forum Member
Joined: May 2004
Posts: 222 |
I am looking just to make sure i have the correct output.
I ... read parm read rtf read psf read coor
write new.psf write new.coor
This is how i regenerate psf. So i want to know if the new psf is what i want. Hence i need to know where are the donor acceptor sections.
Help.
CHARMM 30b1 driven by
1/ Xeon (32 bits)
2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel
3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum.
4/ Gracious support from the forum.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Check for DONOr and ACCEptor statements in the RESIdues of interest in your RTF. What is so difficult about this?
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
The HBOND sections of the PSF begin with lines containing the strings "!NDON: donors" and likewise "!NACC: acceptors", and are listed in pairs as donor,donor-antecedent or acceptor,acceptor-antecedent. It is, however, mostly irrelevant as the HBOND term of the potential isn't used except for analysis, and for that it is trivial to use the DONOR and ACCEPTOR commands to define them on the fly. You don't need a new PSF for that.
Rick Venable computational chemist
|
|
|
|
Joined: May 2004
Posts: 222
Forum Member
|
Forum Member
Joined: May 2004
Posts: 222 |
Thank you all. I know better now. 
CHARMM 30b1 driven by
1/ Xeon (32 bits)
2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel
3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum.
4/ Gracious support from the forum.
|
|
|
|
Joined: Mar 2005
Posts: 42
Forum Member
|
Forum Member
Joined: Mar 2005
Posts: 42 |
Hello!!
I just have a questions about the output of the script. I have a implicit solvent trayectory of 1fs timestep, so I commented the water H-bond part, and I changed the coor hbond setup as follows:
!backbone-backbone coor hbond verbose select type HN .or. type O end select type HN .or. type O end - first 51 nunit 1 begin 100 skip 100 stop 200000
Part of the output is :
READING TRAJECTORY FROM UNIT 51 NUMBER OF COORDINATE SETS IN FILE: 2000 NUMBER OF PREVIOUS DYNAMICS STEPS: 100 FREQUENCY FOR SAVING COORDINATES: 100 NUMBER OF STEPS FOR CREATION RUN: 200000 ... I-atom J-atom (Bridge) Lifetime Endtime PROT 37 LYS HN - PROT 33 ALA O 0.60 1.00 PROT 43 LYS O - PROT 46 ALA HN 0.20 1.10 PROT 51 LEU HN - PROT 51 LEU O 0.20 0.40
I am not sure of understanding the 'lifetime' and 'endtime'. I saved structures every 100. I interpret the first line hbond (37-33) as existing during the time window 0.4-1.0, and in the second case (43-46) during 0.90-1.10, Is that correct?? I also noted cases as the third line were a Hbond is considered beetween C and N of the same residue. I read that one can define CUTHA = 180 to allow only linear Hbonds, nevertheless, the angles of the first structure (showed earlier in the output) range from 170 to 90 degrees, being the average angle ~160. How can i avoid strange angles without being too restrictive?
Thanks
Last edited by Leonardo; 04/23/05 02:09 AM.
Leonardo Sepulveda DurĂ¡n
Universidad de Chile
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Your interpretation of the lifetime and endtime columns is correct. I do not think that you have to worry about the angular dependence - in practice it works very well with just a distance dependent cutoff on the H-Acceptor distance (see J. Am. Chem. Soc., 1992. 114: p. 4028-4035). It may be that more distorted hydrogen bonds do actually occur in implicit solvent simulations (I have not looked at this, so it is just a speculation). Note that with a 100 step resolution you may underestimate the hydrogen bond lifetimes; we normally use a resolution of 1-5 ps.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: May 2004
Posts: 222
Forum Member
|
Forum Member
Joined: May 2004
Posts: 222 |
In determining prot-water-ligand bridges, i've specified 'tip3 .and. resi 136' instead of 'TIP3' to make it more specific. There is no error in the output but i cant be sure it is correct. Can someone confirm? The command are as followed: coor hbond verbose cuthb @CUTHB cutha @CUTHA select ONE end select TWO end- bridge tip3 .and. resi 136 firstu 51 nunit @N skip @SKIP
CHARMM 30b1 driven by
1/ Xeon (32 bits)
2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel
3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum.
4/ Gracious support from the forum.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Incorrect, but you can achieve this effect by renaming tip3 136: rename resn t136 select atom wat 136 * end coor hbond ... bridge t136 ....
It is typically not necessary to define a cutoff for the angle.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jan 2004
Posts: 91
Forum Member
|
Forum Member
Joined: Jan 2004
Posts: 91 |
hi,
i would like to know how to interpret the results of lifetime and endtime. as i read the post on hbonds, it is still not clear with the output i have.
######################################################
I-atom J-atom (Bridge) Lifetime Endtime LEC 201 ASP OD1 - MIC 13 BOG HO4 60.00 70.00 LEC 201 ASP OD2 - MIC 13 BOG HO3 140.00 150.00 LEC 201 ASP OD2 - MIC 13 BOG HO3 120.00 280.00 LEC 201 ASP OD1 - MIC 13 BOG HO4 230.00 310.00 LEC 173 LEU O - MIC 20 BOG HO3 70.00 410.00 LEC 195 ASP OD2 - MIC 8 BOG HO4 90.00 410.00 LEC 195 ASP OD2 - MIC 8 BOG HO4 100.00 520.00 LEC 145 LYS HZ2 - MIC 16 BOG O4 90.00 740.00 LEC 195 ASP OD1 - MIC 8 BOG HO4 150.00 770.00 LEC 108 GLU OE2 - MIC 16 BOG HO2 90.00 800.00 LEC 129 ASP OD1 - MIC 16 BOG HO4 150.00 800.00 LEC 145 LYS HZ2 - MIC 16 BOG O2 190.00 800.00
I-atom J-atom (Bridge) Lifetime Endtime LEC 100 TYR O - MIC 16 BOG O4 WAT1 7106 250.00 260.00 LEC 174 THR O - MIC 20 BOG O6 WAT1 6496 60.00 470.00 LEC 192 THR OG1 - MIC 8 BOG HO6 WAT1 5614 60.00 520.00
#########################################################
as seen in the second and third row the protein and lipid atoms are the same except their two different lifetime and endtimes. which one should i consider in this case. in general:--
the endtime --- is it the actual residence time b/n those atoms?
the lifetime --- is it the duration of active hbond ?
thank you
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
In this section of the verbose-mode output the endtime is the timepoint at which this particular instance of the bond between these two atoms was broken. The lifetime is the duration of this particular instance of the hydrogen bond or bridge. Thus the information in your second and third rows means that there was a hydrogen bond between these atoms in the two intervals (10,150) and (160,280).
You are probably better off looking at the more compact summary information at the end of the output.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jan 2006
Posts: 22
Forum Member
|
Forum Member
Joined: Jan 2006
Posts: 22 |
Hello,
I also have questions of the output od hydrogen bond analysis even I've read the above posted but it's still unclear to me.
In protein-proetin output section, I have:
Analysis of hydrogen bonds using cutoff distance= 2.40 and angle= 999.00 I-atom J-atom r(H-A) [A] D-H..A [deg] -------------------------------------------------------------------- PROT GLY 4 O - PROT LYS 56 HN 1.90 168.7 PROT GLY 4 O - PROT THR 57 HN 1.86 157.8 PROT TRP 6 HN - PROT THR 57 O 1.88 169.8 PROT TRP 6 O - PROT HSD 59 HN 1.99 163.5 PROT GLN 7 HN - PROT VAL 22 O 1.79 167.7 PROT GLN 7 O - PROT VAL 22 HN 1.89 170.1 PROT LEU 8 HN - PROT HSD 59 O 1.77 175.0 PROT THR 11 HN - PROT ILE 18 O 1.86 165.5 PROT THR 11 O - PROT ILE 18 HN 1.91 159.1
and the lifetime section is:
I-atom J-atom (Bridge) Lifetime Endtime PROT 2 SER HN - PROT 2 SER O 5.00 10.00 PROT 2 SER O - PROT 2 SER HN 5.00 10.00 PROT 5 ILE O - PROT 24 VAL HN 5.00 10.00 PROT 17 VAL HN - PROT 17 VAL O 5.00 10.00 PROT 17 VAL O - PROT 17 VAL HN 5.00 10.00 PROT 24 VAL HN - PROT 5 ILE O 5.00 10.00 PROT 41 GLU O - PROT 45 PHE HN 5.00 10.00 PROT 45 PHE HN - PROT 41 GLU O 5.00 10.00
My questions are: 1. For lifetime section, Does the first line mean that the bridge between HN of SER2 and O of the same residue occurs within 5ps and it will be broken after 10 ps? I just want to make sure whether I'm correctly understand.
2. Why the lines of the first section (protein-protein output section in this case) are not match with the lines of the lifetime section? For example,
If the protein-protein section starts with
PROT GLY 4 O - PROT LYS 56 HN 1.90 168.7 PROT GLY 4 O - PROT THR 57 HN 1.86 157.8 PROT TRP 6 HN - PROT THR 57 O 1.88 169.8
The lifetime section should also start with:
I-atom J-atom (Bridge) Lifetime Endtime
PROT GLY 4 O - PROT LYS 56 HN ... .... PROT GLY 4 O - PROT THR 57 HN ... .... PROT TRP 6 HN - PROT THR 57 O ... ....
Thanks in advance.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
1/ Yes, your understanding is correct. 2/ You are apparently looking at two different COOR HBOND outputs. The one you call protein-protein is from a static analysis of a single coordinate set; a similar output, but with hydrogen bond lifetime data, is printed from the analysis of a trajectory. Your other output, which you call lifetime section, is simply a chronological record of all hydrogen bond forming/breaking events in the trajectory, as such it is time-ordered; it is usually not necessary to look at this output. The summary information, giving the average occupancies and lifetimes is more informative and easier to grasp.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jan 2006
Posts: 22
Forum Member
|
Forum Member
Joined: Jan 2006
Posts: 22 |
Hello,
Thank you Prof. Lennart. However I still have questions.
I've looked at the summary section and I've found something like this:
I-atom J-atom <occupancy> <time> # events --------------------------------------------------------------------------- PROT SER 2 HN - PROT SER 2 O 0.44 18.3 6 PROT SER 2 O - PROT SER 2 HN 0.44 18.3 6 PROT SER 2 O - PROT GLY 4 HN 0.06 7.5 2 PROT SER 2 O - PROT ILE 5 HN 0.28 7.8 9 PROT SER 2 O Summary: 0.78 11.2 17 PROT GLY 4 HN - PROT SER 2 O 0.06 7.5 2 PROT GLY 4 O - PROT LYS 56 HN 1.00 250.0 1 PROT GLY 4 O - PROT THR 57 HN 0.98 122.5 2 PROT GLY 4 O Summary: 1.98 186.2 3 PROT ILE 5 HN - PROT SER 2 O 0.28 7.8 9 PROT ILE 5 O - PROT VAL 24 HN 0.08 5.0 4 PROT TRP 6 HN - PROT TRP 6 O 0.10 5.0 5 PROT TRP 6 HN - PROT THR 57 O 1.00 250.0 1 PROT TRP 6 HN Summary: 1.10 127.5 6
My questions are:
1. what's the meaning of time and events? Are the time and event refer to the duration of h-bond (ps) and a numbers of time that we can find this h-bond, respectively? And what's the maximum value of "summary of occupancy" per residue?
2. In the last section of this summary, I've found:
OVERALL Average number (over ISEL!)= 0.6 <lifetime>= 83.52ps Total number/frame= 175.3
Is the "lifetime" refer to the total lifetime or what?
Thanks in advance.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
<...> means average. The symbol # means "number". So "<time>" means "average [life]time" (in ps) and "# events" means "number [of] events", ie how many times this hydrogen bond was formed. I do not know exactly what the maximum number of the occupancy would be since it depends on how many simultaneous hydrogen bonds that particular atom can form; the summary is not per residue. The final numbers are the averages over all selected donors/acceptors over the whole trajectory. Check corman.doc.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jan 2006
Posts: 22
Forum Member
|
Forum Member
Joined: Jan 2006
Posts: 22 |
Thank you to Prof. Lennart for your clarification.
For my question below, I'm not sure whether I should post it here or should be in another section.
I just have problem on h-bond analysis when reading more than 1 trajectory. For example, I used charmm 29b2 to simulate the first trajectory and used c32b1 for the other trajectories (because we just have c32b1).
I can not use neither c29b2 nor c32b1 to analyse all trajectories produced by both c29b2 and c32b1.The header produced by c29b2 and c32b1 does not match each other. However, if I use c29b2 to analyse only trajectory produced by c29b2 or vice versa for c32b1, it works well. I'd like to know that are there any possiblities to read trajectories produced by different versions of charmm wihtout any configtion.
Thanks in advance.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Your CHARMM versions probably have been compiled with different settings for little_endian vs big_endian representation of binary data in unformatted files; for c32b1 on gnu/linux systems the default is big_endian with the intel and pgi compilers. It is easy to modify install.com to get the native format /little_endian/ to be used instead. If your system uses the intel compiler this behavior can be changed at runtime without reinstallation; you can also convert the files you already have. Search the forums, or look at some recent posts on format of dcd files.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jan 2006
Posts: 22
Forum Member
|
Forum Member
Joined: Jan 2006
Posts: 22 |
Thank you very much Prof. Lennart. I've checked the format dcd file, as you mentioned. In my case c29b2 writes little_endian and c32b1 writes big_endian. After I changed dcd file (from c32) in big_endian format to little_endian (or little, c29, to big) but still can not analyse both trajectories together.
PROT 110 VAL HN - PROT 106 ILE O 5.00 250.00 PROT 111 ARG O - PROT 114 ALA HN 5.00 250.00 PROT 114 ALA HN - PROT 111 ARG O 5.00 250.00 ***** ERROR ***** HEADERS DO NOT MATCH Ă´^A^@^@ CORD VELD *** LEVEL -1 WARNING *** BOMLEV IS -3
I have no idea how to deal with the header. Please help. Thank you.
|
|
|
|
Joined: Nov 2005
Posts: 51
Forum Member
|
Forum Member
Joined: Nov 2005
Posts: 51 |
Dear CHARMMERS,
From the hbond analysis output, I have seen some part of the output which looks like this: Analysis of hydrogen bonds using cutoff distance= 2.80 and angle= 999.00 Total time= 600.0ps. Resolution= 5.00ps Occupancy cutoff= 0.00 Lifetime cutoff= 0.0ps
I-atom J-atom <occupancy> <time> # events --------------------------------------------------------------------------- PROT SER 2 O - PROT GLY 4 HN 0.62 11.9 31 PROT SER 2 O - PROT ILE 5 HN 0.83 29.1 17 PROT SER 2 O Summary: 1.44 20.5 48 PROT PRO 3 O - PROT ILE 5 HN 0.02 5.0 2 PROT GLY 4 HN - PROT SER 2 O 0.62 11.9 31 PROT GLY 4 O - PROT LYS 56 HN 1.00 600.0 1 PROT GLY 4 O - PROT THR 57 HN 1.00 600.0 1 PROT GLY 4 O Summary: 2.00 600.0 2
If I would like to the find % occupation of hbond , for example, % hbond occupation of the O-Ser2 and other residues, how can I obtain it?
I am not sure if I can use "1.44*100" because it's higher than 100%. If I think that 100% means this kind of hbond occurs entire the simulation time, how can I calculate the % hbond occupation? Any advice would be appreciated. Thank you.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
The carbonyl of Ser 2 apparently hydrogen bonds to the amides of Gly4 and Ile5, partly at the same time. 1.44*100 is correct if you want to express this total occpancy as a percentage.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Nov 2005
Posts: 51
Forum Member
|
Forum Member
Joined: Nov 2005
Posts: 51 |
Thank you Prof. Lennart. From your reply, it means that % hbond occupation can higher than 100%, is my understanding correct? First I though that the highest value could be 100% but I am not sure.
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
The maximum fractional occupancy is 1.0 when there is a 1:1 correspondence between donor and acceptor, but a good H-bond acceptor can have 2 donors, and the fractional occupancy can exceed 1.0 in that case.
Rick Venable computational chemist
|
|
|
|
Joined: Oct 2003
Posts: 39
Forum Member
|
Forum Member
Joined: Oct 2003
Posts: 39 |
HI,
I will do this by hand.
Last edited by zhangxd68; 02/15/07 05:22 AM.
|
|
|
|
Joined: Feb 2006
Posts: 4
Forum Member
|
Forum Member
Joined: Feb 2006
Posts: 4 |
Hi,
1) When i was doing analysis for a trajectory for a 20ns simulation with 2000 frames with: "coor hbond verbose sele monitor end sele .not. monitor end firstu 20 nunits 1 begin 10000 skip 10000"
the following was observed in the output: Analysis of hydrogen bonds using cutoff distance= 2.40 and angle= 999.00 Total time= 977.8ps. Resolution= 0.49ps Occupancy cutoff= 0.00 Lifetime cutoff= 0.0ps
>> why only 977.8 ps?
2) the spec "skip" and "nskip" for the traj-spec is same?
regards, Bryan
|
|
|
|
Joined: Jul 2004
Posts: 165
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 165 |
Dear Prof. Nilsson, I would like to analyze the hydrogen bond between lipid and peptide to compare with the literature. The criteria that they use is the hydrogen bond said to exist if donor -acceptor distance is < 2.5 angstrom and the donor-hydrogen..acceptor hydrogen bond form angle > 145 degree. to set that, I would say cuthb 2.5, but to set the angle I could not simply say cutha 145, correct? how could I reach that angle? thank you
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Why not? From corman.doc: "[hydrogen bond exists if ] the angle formed by D-H..A is greater than this CUTAngle (in degrees, 180 is a linear H-bond)".
Note that adding the angular criterion makes almost no difference in practice.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jul 2004
Posts: 165
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 165 |
perhaps it is a naive question. Quote:
It is typically not necessary to define a cutoff for the angle.
Since there are many difference criteria have been used in the lipid simulation, is this criteria is also reasonable for the lipid system. I am not so sure becuse we are here mainly focused on the peptide/protein systems. thank you very much
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
The forces between atoms are computed in CHARMM w/o regard to the molecule type; for H-bonds, the partial charges and VDW radii play a major role, esp. in the close approach distance.
If you are comparing to another simulation result, you should probably try to use the same distance and angle criteria.
Otherwise, I suggest trying the COOR HBOND defaults; we've done that in a publication to evaluate lipid:water and lipid:sugar H-bonds.
Making sure you've properly accounted for image interactions may be a more important issue, esp. if the peptide reaches the edge of the simulation cell during the simulation.
|
|
|
|
Joined: Mar 2007
Posts: 23
Forum Member
|
Forum Member
Joined: Mar 2007
Posts: 23 |
Hi, I used the code listed here for estimating the number of hydrogen bonds and lifetime. I am wondering is there any way to get the average distance between the donors and acceptors as well ?
Thanks a lot Diwakar
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
As of CHARMM version 32 the COOR HBONd command can also give the distance distribution, from which you should easily be able to compute the average. From corman.doc: "Distance and lifetime histograms can be computed for all (putative) hydrogen bonds encountered in the analysis"
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Mar 2007
Posts: 23
Forum Member
|
Forum Member
Joined: Mar 2007
Posts: 23 |
Thanks a lot for the help. I have a couple more questions about calculating the number and lifetime of hydrogen bonds. 1) Is it possible to exclude the intra-molecular hbonds? I calculated the hbonds by selecting each molecule separately and then counting the hbonds of this molecule with other molecule. This process excludes the intra-molecular bonds but I was wondering if there is a better way of doing it.
2) How can I get the number of hbonds as a function of time?
Regards, Diwakar
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
1/ You did the right thing, there is no other way to do it. Using the two selections is not all that difficult, is it? 2/ Use a COOR HBOND (without keyword VERBose) in a CHARMM loop and write the number of hydrogen bonds (available as a CHARMM substitution variable, see subst.doc) to a file.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Mar 2007
Posts: 23
Forum Member
|
Forum Member
Joined: Mar 2007
Posts: 23 |
I am trying to get water water hydrogen bonds. I have added the donor list to TIP3 water model as described in the hbonds.inp. I get the same psf without any extra donors on regenerating the psf using this code:
read rtf card name top_all27.rtf read para card name par_all27.prm
read psf card name @pdbfile@qual.psf read coor card name @pdbfile@qual.crd
write psf card name @pdbfile@qual-water-donor.psf write coor card name @pdbfile@qual-water-donor.crd
I want to know if anything more needs to be done to regenrate the psf?
Thanks Diwakar
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
Yes, you must use the GENERATE command instead of READ PSF
An alternative is the DONOR and ACCEPTOR commands; they seem to work as long as there only one of each in a given CHARMM input.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
DONOr/ACCEptor ADD have some issues, and do not work properly in older versions of CHARMM. It is really very easy to regenerate the PSF so I strongly recommend that route; if your water segment is the last one all it takes is: read rtf card name ... ! the new RTF w/ DONO, ACCE for TIP3 read param card name .. read psf card name old.psf delet atom sele segid wat read seque tip3 345 gene wat noangle nodihe write psf card name new.psf * PSF for ... with DONO/ACCE for TIP3 *
If water residue numbering in the original PSF was not sequential one may change it to sequential with the JOIN command (struct.doc).
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
A good point; you don't have to start over from the beginning of your model building. What's most important is that the number and ordering of atom types in the PSF does not change. Reading the old PSF and deleting the water segment alone won't work if additional segments were added after the water, such as ions. In that case, you'd need to also delete the segments that follow the water segment in the PSF, and then add them back via GENERATE, in the same order as was done initially.
|
|
|
|
Joined: Jul 2004
Posts: 165
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 165 |
Dear all,
Sorry if this is a naive question, but I feel a little bit confused about ?NHBOND and ?AVNOHB.
I want to calculate the total number of hydrogen bond between lipid and water and average number of hydrogen bonds per lipid.
Could I just simply say
coor hbond sele resn lip end - sele resn tip3 end firstu 18 nunit 1 skip 5000 echo ?NHBOND to get the total number of hydrogen bond?
However, I am not sure if I can apply just ?avnohb to get the correct value for number of hydrogen bond per residue. In addition, is it possible with the charmm command to get the standard deviation for that?
thank you
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Yes, ?NHBOND gives the total number The average for ?avnohb is however over selected atoms, not residues (but the conversion should be easy. No, there is no additional statistics available; it is straightforward to output the numbers of interest for relevant timeframes using a CHARMM loop and post-process that with your favorite graphics/spreadsheet/databases/whatever program. Plese not that in this case the VERBose keyword for CORR HBON should not be used.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Mar 2007
Posts: 23
Forum Member
|
Forum Member
Joined: Mar 2007
Posts: 23 |
Sorry for the late response. I tried to generate the new psf file based on the suggestions by lennart and rmv. The system consists of the solute molecules in a water box. I started from a water box and then I add the solute molecules randomly. After the addition of solute molecules, I remove overlapping water molecules. If I setup the system again using the same procedure, the solute molecules are placed randomly again. Therefore, it is not possible to generate the same system again. However, I found a quick way to add water donors in the psf. I edited the donor section of the psf using awk and it works. I need to change the Ndonors and the list of the donor atoms only.
Thanks a lot for all your suggestions.
Diwakar
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Hand editing of the PSF should not be necessary. Note that for regenerating the PSF you do not need to perform any coordinate manipulations, all you need to know is how many molecules you have of a certain type, not where they are. Let's say you start with 100 water molecules, then add 30 solute molecules and delete 50 water molecules. To regenerate this PSF you can either read seque tip3 50 gene wat noangle nodihe read seque myml 30 gene solu
Sequences could also be read from coordinate files (one for each segment) prepared from the final state of your inital setup.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Jul 2004
Posts: 165
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 165 |
could you please give me a hint how to get the hydrogen bond autocorrelation function of i.e. N–H---O=C atoms ?
Thank you in advance.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Not directly implemented, but there are a couple of thing you can try.
1. If you are interested in specific donor-acceptor pairs you can use CORRel to extract the timeseries of the distance between the atoms, and the process that with your own program.
2. The MRD option of coor anal (corman.doc) does calculate the intermittent correlation function (ie, no requirement that the "donor"/"acceptor" are within cutoff distance at all times). Something like this (not tested): COOR ANAL SELE first-atom(s) END SITE SELE second-atom END - FIRSTU 51 NUNIT 2 NCORS 100 IMRD 21 RRES 3.0
This should output the correlation function in the fourth column for atom(s) in the first selection around the atom in the second selection.
Note that correlation times extracted from these correlation functions are not the same thing (and possiblby more approximate) as the (exact) averages calculated by COOR HBONd. You can also get lifetime distributions from COOR HBONd, which should be related.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
1. If you are interested in specific donor-acceptor pairs you can use CORRel to extract the timeseries of the distance between the atoms, and the process that with your own program.
It should be noted that CORREL can also easily calculate an autocorrelation function from a distance time series; there's no need to write a program for this.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Well, in this case it is usually not <r(t)r(t+T)> that is sought (not even <dr(t)dr(t+T)>), but rather something like <h(t)h(t+T)>, where h=1 if r<rcut, and 0 otherwise.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
It appears that one of the pieces needed is present in c33 and later, but was never documented. What's missing from correl.doc is a description of the following MANTime command--
STATE realmin realmax
Returns a state function in Q(t), with 1.0 for values between realmin and realmax, and 0.0 for values outside the range.
|
|
|
|
Joined: Nov 2019
Posts: 52
Forum Member
|
Forum Member
Joined: Nov 2019
Posts: 52 |
Hello!
! ACCEPTOR OH2 ! DONOR H1 OH2 ! DONOR H2 OH2
Regarding with adding the donor and acceptor list for water, I'm wondering whether this also applies to SWM4 drude model. There is only ACCEPTOR OH2 in original SWM4.
RESI SWM4 0.000 ! SWM4-NDP water GROUP ATOM OH2 ODW 0.00000 TYPE DOH2 ALPHA -0.97825258 THOLE 1.3 ATOM OM LPDW -1.11466 ATOM H1 HDW 0.55733 ATOM H2 HDW 0.55733 BOND OH2 H1 BOND OH2 H2 BOND OH2 OM BOND H1 H2 ! for SHAKE ANGLE H1 OH2 H2 ACCEPTOR OH2
Last edited by Allen_123; 06/25/20 09:56 AM.
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
Yes. coor hbond will only recognize hbonds between atoms that have been designated as DONOrs and ACCEptors.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
|