*FILENAME: hbonds.inp
*PURPOSE: compute static and time-dependent hydrogen bonds from trajectory
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
*MODIFIED: Mike McCallum UOP Oct. 2004
*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set rtf /usr/local/charmm/toppar/top_all22_prot_cmap
!!
set prm /usr/local/charmm/toppar/par_all22_prot_cmap
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open unit 1 card read name  -  
@rtf.inp
read RTF card unit 1  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open unit 2 card read name -  
@prm.inp
read PARA card unit 2  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
!!
!! Input files/etc
if @?input eq 0 goto endoffile
if @?begin eq 0 goto endoffile
open unit 4 card read name @input.psf
read PSF card unit 4
!!

open unit 1 read card name @input.pdb
read coor pdb unit 1
close unit 1

! 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 end select type O end

! now get the time-dependent hydrogen bonds; 
! the output can be voluminous, especially
! with the verbose option

open unit 51 read unform name @input.dcd
traj query unit 51

!open unit 52 read unform name a10-a-stt2.dcd
!! any other files opened in sequence:
!!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)

! Note that by default, HBOND does not specificy a threshold
! angle (cutha).  This is too generous.  I reccomend 160 deg (corresponding
! to a VMD angle of 20). The default distance (O -- H, not O---N)
! distance is 2.4 A (cuthb).
!
! Representative CHARMM output (in the log file):
!
! Analysis of hydrogen bonds using cutoff distance=    2.40 and angle= 160.00
! Total time=    400.0ps.  Resolution=   1.00ps
! Occupancy cutoff= 0.00  Lifetime cutoff=   0.0ps
!
!     I-atom                 J-atom          <occupancy>     <time>  # events
! ---------------------------------------------------------------------------
!      PRO0 ALA  3    HN   -  PRO0 ALA  1    O     0.002        1.0       1
!      PRO0 ALA  4    HN   -  PRO0 ALA  1    O     0.030        1.1      11
!      PRO0 ALA  5    HN   -  PRO0 ALA  1    O     0.198        1.8      43
!      PRO0 ALA  5    HN   -  PRO0 ALA  2    O     0.032        1.1      12
!                PRO0 ALA  5    HN       Summary:  0.230        1.5      55
!
!
! The SKIP is also important.  I have found a SKIP of the file interval works best.
! The charmm log will tell you this if you use a too-small value:
!
! *****  WARNING  ***** SKIP=  100 WAS NOT A MULTIPLE OF THE FILE INTERVAL=  500
! IT HAS BEEN RESET TO  500
!
! Experiment with this so you get a high density of time points.  (in other words,
! At
! some (low) value of SKIP, you will not get any more points.  NOTE
! this will allow very short HB lifetimes!
!
open unit 21 write card name @input-hbond.out
coor HBOND cuthb 2.4 cutha 160 verbose -
select type HN end select type O end -
firstu 51 nunit 1 begin @begin skip 50 IUNIT 21
close unit 21

! note that the trajectory is not closed between these commands, and can be reused
! protein - solvent
!!! deleted:
!!!coor hbond verbose select segid prot end select segid wat end -
!!!firstu 51 nunit 2 begin 10500 skip 2500

! protein - water - ligand
!!! deleted:
!!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

label endoffile
system "echo 'Must specify pdb filename'"

stop