* calculate various DNA hydration numbers
*

faster on

! identifiers
set id1 ecor_sod
set id2 test
set range 500_1000

set trajf 950  ! final traj identifier
set traji 50   ! traj increment value
set file1 500    ! initial traj file

set skip  500
set begin 250500   ! initial time frame
set begin2 250500
set stop 500000 ! final time frame

calc nfile = (@stop - @begin + @skip)/@skip

!total number of nucleotides in one strand
set restot 12

set nf 50    ! number of frames per analysis section

bomlev -2

! using the real space trunction scheme without the PME correction

! 1: topology
! 2: parameter
! 3: cutim
! 4: cutnb
! 5: ctonnb
! 6: ctofnb
! 7: electrostatic cutoff regime (shift or switch)
! 8: electrostatic pairing method (atom or group)
! 9: VDW pairing method (vatom or vgroup)
! VDW cutoff regime is vwsitch
set 1 top_all27_na.rtf
set 2 par_all27_na.prm
set 3 4.0 ! cutim
set 4 4.0 ! cutnb
set 5 4.0 ! ctonnb
set 6 4.0 ! ctofnb
set 7 switch
set 8 atom
set 9 vatom

! read topology and parameters
open unit 9 read form name @1
read rtf card unit 9

open unit 9 read form name @2
read para card unit 9

open unit 11 read form name @id1_@id2.psf
read psf cards unit 11

! We read coordinates of ecor1 to allow correl to work
open read unit 12 card name @id1_@id2_b.crd
read coor card unit 12

! cell parameters
set a 60.8529
set b 36.9108
set c 36.9108

crystal define orthorhombic @a @b @c 90. 90. 90.

crystal build noper 0 cutoff @3

!open output files (interaction energies of each window)
!open unit 91 write form name solv_dna_total_@id2_@range.hyd
!open unit 92 write form name solv_dna_elec_@id2_@range.hyd
!open unit 93 write form name solv_dna_vdw_@id2_@range.hyd
open unit 94 write form name solv_dna_bas_@id2_@range.hyd
open unit 95 write form name solv_dna_bkb_@id2_@range.hyd

! SET UP ATOM DEFINITIONS
define GBASE  sele resn gua .and. (type N1 .or. type C2 .or. type N2 -
.or. type N3 .or. type H8 -
.or. type C4 .or. type C5 .or. type C6 .or. type O6 .or. type N7 -
.or. type C8 .or. type N9 .or. type H21 .or. type H22 .or. type H1) end

define ABASE  sele resn ade .and. (type N1 .or. type C2 .or. type H2 -
.or. type N3 .or. type H8 -
.or. type C4 .or. type C5 .or. type C6 .or. type N7 -
.or. type C8 .or. type N9 .or. type H61 .or. type H62 .or. type N6) end

define CBASE  sele resn cyt .and. (type N1 .or. type C2 -
.or. type N3 .or. type C4 .or. type C5 .or. type C6 .or. type O2 -
.or. type N4 .or. type H41 .or. type H42 -
.or. type H5 .or. type H6) end

define TBASE  sele resn thy .and. (type N1 .or. type C2 -
.or. type N3 .or. type C4 .or. type C5 .or. type C6 .or. type O2 -
.or. type O4 .or. type C5M .or. type H51 .or. type H52 .or. type H53 -
.or. type H6 .or. type H3) end

define UBASE  sele resn ura .and. (type N1 .or. type C2 -
.or. type N3 .or. type C4 .or. type C5 .or. type C6 .or. type O2 -
.or. type O4 .or. type H3 -
.or. type H5 .or. type H6) end

!ALL solute ATOMS
define dna sele -
         segi dna1:dna2 .and. .not. -
         (hydr .or. resi 1 .or. resi @restot) show end

!define st1 sele segi dna1 end
!define st2 sele segi dna2 end

! BASES 
!define base sele -
!        ((segi dna1:dna2) .and. -
!        (gbase .or. cbase .or. ubase .or. abase .or. tbase)) -
!	end

!BACKBONE 
!define bkb sele -
!        (segi dna1:dna2) .and. -
!        (type P .or. type O1P .or. type O2P .or. -
!        type %%' .or. type %%'' .or. -
!        type %%T) show end

define pur sele -
      (gbase .or. abase) .and. -
      .not. (resi 1 .or. resi @restot) show end

define pyr sele -
      (cbase .or. tbase .or. ubase) .and. -
      .not. (resi 1 .or. resi @restot) show end

define atom_c8 sele -
      type c8 .and. (gbase .or. abase) .and. -
      .not. (resi 1 .or. resi @restot) show end

define atom_c6 sele -
      type c6 .and. (cbase .or. tbase .or. ubase) .and. -
      .not. (resi 1 .or. resi @restot) show end

!Ribo
!define ribose sele -
!      ((type %%' .or. type %%'') .and. -
!      .not. (hydrogen .or. type O5' .or. type O3')) .or. -
!      (segid dna1 .and. resi 1 .and. (type o5' .or. type h5t)) .or. -
!      (segid dna1 .and. resi @restot .and. (type o3' .or. type h3t)) .or. -
!      (segid dna2 .and. resi 1 .and. (type o5' .or. type h5t)) .or. -
!      (segid dna2 .and. resi @restot .and. (type o3' .or. type h3t)) show end

define sugar_o sele -
      (type o4') .and. .not. -
      ((segid dna1 .and. resi 1 .and.  type o4') .or. -
      (segid dna1 .and. resi @restot .and. type o4') .or. -
      (segid dna2 .and. resi 1 .and. type o4') .or. -
      (segid dna2 .and. resi @restot .and. type o4')) show end

define ester_o sele - 
      (type o3' .or. type o5') .and. .not. -
      ((segid dna1 .and. resi 1 .and. (type o5')) .or. -
      (segid dna1 .and. resi @restot .and. (type o3')) .or. -
      (segid dna2 .and. resi 1 .and. (type o5')) .or. -
      (segid dna2 .and. resi @restot .and. (type o3'))) show end

define anionic_o sele -
      (type o1p .or. type o2p) show end

!phosphate atom
!define phosphate sele -
!      (type P .or. type O1P .or. type O2P -
!      .or. type O5' .or. type O3') - 
!      .and. .not. ribose show end

calc restotn1 = @restot - 1

define major sele -
                  resi 2:@restotn1 .and. -
                  ((resn gua .and. (type n7 .or. type o6)) .or. -
                  (resn cyt .and. type n4) .or. -
                  (resn ade .and. (type n7 .or. type n6)) .or. -
                  (resn thy .and. type o4)) show end

define minor sele -
                  resi 2:@restotn1 .and. -
                  ((resn gua .and. (type n3 .or. type n2)) .or. -
                  (resn cyt .and. type o2) .or. -
                  (resn ade .and. (type n3 .or. type c2)) .or. -
                  (resn thy .and. type o2)) show end

set unit 30
set time @file1 
set timeinc @traji 
set timef @trajf 

label open_traj

open unit @unit read unform name @id1_@id2_@time.trj
incr unit by 1
incr time by @timeinc

if time le @timef goto open_traj

! determine averages over 500 ps blocks (every 10 ps) and average.  therefore
! omit determination of std dev.

set hyddna 0
set hydpyr 0
set hydpur 0
set hydc8 0
set hydc6 0
set hydmaj 0
set hydmin 0
set hydo4 0
set hydester 0
set hydanionic 0

! detemine number of atoms for the hydration determination
! for normalization purposes
coor stat -
sele dna end
set natomdna ?nsel

coor stat -
sele pyr end
set natompyr ?nsel

coor stat -
sele pur end
set natompur ?nsel

coor stat -
sele atom_c8 end
set natomc8 ?nsel

coor stat -
sele atom_c6 end
set natomc6 ?nsel

! major groove excluding terminal base pairs
coor stat -
sele major end
set natommaj ?nsel

coor stat -
sele minor end
set natommin ?nsel

coor stat -
sele sugar_o end
set natomo4 ?nsel

coor stat -
sele ester_o end
set natomester ?nsel

coor stat -
sele anionic_o end
set natomanionic ?nsel

write title unit 94
* natom @natompyr @natompur @natomc8 @natomc6 @natommin @natommaj
*

write title unit 95
* natom @natomdna @natomo4 @natomester @natomanionic
*

!counters for the inner loop (read_traj)
set i 1		!step counter
!set time
set ti 1.0	!time counter

set stepi @begin
set stepf @stop

!loop thru the traj file
traj iread 30 nread 10 begin @stepi stop @stepf skip @skip

label read_traj

traj read

! check need for update
update inbfrq -1 imgfrq -1 ihbfrq 0 -
@7 @8 @9 vswitch cutimg @3 cutnb @4 ctofnb @6 ctonnb @5 !bycube

! detemine hydration numbers
coor stat -
sele type oh2 .and. -
(dna .around. 3.50) end
calc hyddna = @hyddna + ?nsel

coor stat -
sele type oh2 .and. -
(pyr .around. 3.50) end
calc hydpyr = @hydpyr + ?nsel

coor stat -
sele type oh2 .and. -
(pur .around. 3.50) end
calc hydpur = @hydpur + ?nsel

coor stat -
sele type oh2 .and. -
(atom_c8 .around. 3.50) end
calc hydc8 = @hydc8 + ?nsel

coor stat -
sele type oh2 .and. -
(atom_c6 .around. 3.50) end
calc hydc6 = @hydc6 + ?nsel

coor stat -
sele type oh2 .and. -
(major .around. 3.50) end
calc hydmaj = @hydmaj + ?nsel

coor stat -
sele type oh2 .and. -
(minor .around. 3.50) end
calc hydmin = @hydmin + ?nsel

coor stat -
sele type oh2 .and. -
(sugar_o .around. 3.50) end
calc hydo4 = @hydo4 + ?nsel

coor stat -
sele type oh2 .and. -
(ester_o .around. 3.50) end
calc hydester = @hydester + ?nsel

coor stat -
sele type oh2 .and. -
(anionic_o .around. 3.50) end
calc hydanionic = @hydanionic + ?nsel

!stop

incr i by 1 

prnlev -5
if i le @nf goto read_traj
prnlev 5

! calculate statistics and write data
! total energies

calc hyddna = @hyddna/@nf/@natomdna
calc hydpyr = @hydpyr/@nf/@natompyr
calc hydpur = @hydpur/@nf/@natompur
calc hydc8 = @hydc8/@nf/@natomc8
calc hydc6 = @hydc6/@nf/@natomc6
calc hydmin = @hydmin/@nf/@natommin
calc hydmaj = @hydmaj/@nf/@natommaj
calc hydo4 = @hydo4/@nf/@natomo4
calc hydester = @hydester/@nf/@natomester
calc hydanionic = @hydanionic/@nf/@natomanionic

write title unit 94
* @range @hydpyr @hydpur @hydc8 @hydc6 @hydmin @hydmaj
*

write title unit 95
* @range @hyddna @hydo4 @hydester @hydanionic
*

stop