CHARMM Development Project Forums CHARMM Community Script Archive hbonds.inp

*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

Hi.
Which column in the psf specifies donor/acceptor status?
Thanks.

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

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

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.

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.

Check for DONOr and ACCEptor statements in the RESIdues of interest in your RTF. What is so difficult about this?

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.

Thank you all. I know better now.

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