Previous Thread
Next Thread
Print Thread
Page 1 of 6 1 2 3 4 5 6
#243 10/23/03 09:33 AM
Joined: Sep 2003
Posts: 4,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
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
lennart #244 03/04/05 08:26 AM
Joined: May 2004
Posts: 222
C
Forum Member
Offline
Forum Member
C
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.
chuan #245 03/04/05 05:33 PM
Joined: Sep 2003
Posts: 8,650
Likes: 26
rmv Online Content
Forum Member
Online Content
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

rmv #246 03/15/05 07:11 AM
Joined: May 2004
Posts: 222
C
Forum Member
Offline
Forum Member
C
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.
chuan #247 03/15/05 07:36 AM
Joined: Sep 2003
Posts: 4,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
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
lennart #248 03/15/05 09:31 AM
Joined: May 2004
Posts: 222
C
Forum Member
Offline
Forum Member
C
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.
chuan #249 03/15/05 12:49 PM
Joined: Sep 2003
Posts: 4,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
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
chuan #250 03/15/05 08:25 PM
Joined: Sep 2003
Posts: 8,650
Likes: 26
rmv Online Content
Forum Member
Online Content
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

chuan #251 03/16/05 01:40 AM
Joined: May 2004
Posts: 222
C
Forum Member
Offline
Forum Member
C
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.
lennart #252 04/23/05 01:53 AM
Joined: Mar 2005
Posts: 42
L
Forum Member
Offline
Forum Member
L
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
Page 1 of 6 1 2 3 4 5 6

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~deb10u3 Page Time: 0.019s Queries: 35 (0.013s) Memory: 0.7878 MB (Peak: 0.8839 MB) Data Comp: Off Server Time: 2023-06-05 19:44:47 UTC
Valid HTML 5 and Valid CSS