|
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
|
|
|
|
|