* calculate various Solvent-DNA interaction energies using
* 10 A truncation cutoff. 
*

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

set nf 50    ! number of frames per analysis section

! set variables based on particular sequence
set restot 12
set first1 1
calc first2 = @restot + 1
calc last1 = @restot
calc last2 = @restot * 2

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 8.0 ! cutim
set 4 8.0 ! cutnb
set 5 8.0 ! ctonnb
set 6 8.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_@dist_@runid.inter
open unit 92 write form name solv_dna_elec_@dist_@runid.inter
open unit 93 write form name solv_dna_vdw_@dist_@runid.inter
open unit 94 write form name solv_dna_bas_@dist_@runid.hyd
open unit 95 write form name solv_dna_bkb_@dist_@runid.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 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

!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 16 .and. (type o3' .or. type h3t)) .or. -
      (segid dna2 .and. resi 1 .and. (type o5' .or. type h5t)) .or. -
      (segid dna2 .and. resi 16 .and. (type o3' .or. type h3t)) 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

define major sele (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 (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 0
set timeinc 50
set timef 450

label open_traj

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

if time le @timef goto open_traj

! determine averages over 100 ps blocks and average.  therefore
! omit determination of std dev.

set ener1 0
set elec1 0
set vdw1 0
set int1 0
set ener2 0
set elec2 0
set vdw2 0
set int2 0
set ener3 0
set elec3 0
set vdw3 0
set int3 0
set ener4 0
set elec4 0
set vdw4 0
set int4 0
set ener5 0
set elec5 0
set vdw5 0
set int5 0
set ener6 0 !no internal terms
set elec6 0
set vdw6 0
set ener7 0
set elec7 0
set vdw7 0
set ener8 0
set elec8 0
set vdw8 0

set hyddna 0
set hydpho 0
set hydrib 0
set hydmaj 0
set hydmin 0
set hydbas 0
set hydbkb 0

! detemine number of atoms for the hydration determination
! for normalization purposes
! full dna
coor stat -
sele segid dna1:dna2 .and. .not. hydrogen show end
set natomdna ?nsel

! full backbone
coor stat -
sele bkb show end
set natombkb ?nsel

! minor groove excluding terminal base pairs
coor stat -
sele minor .and. .not. (ires @first1 .or. ires @last1 .or. ires @first2 .or. ires @last2) -
show end
set natommin ?nsel

! major groove excluding terminal base pairs
coor stat -
sele major .and. .not. (ires @first1 .or. ires @last1 .or. ires @first2 .or. ires @last2) -
show end
set natommaj ?nsel

! phosphate
coor stat -
sele phosphate show end
set natompho ?nsel

! ribose 
coor stat -
sele ribose show end
set natomrib ?nsel

! all bases except termini 
coor stat -
sele ((abase .or. gbase .or. ubase .or. tbase .or. cbase) .and. -
.not. (hydrogen .or. ires @first1 .or. ires @last1 .or. ires @first2 .or. ires @last2)) -
show end
set natombas ?nsel

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

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

label read_traj

traj read

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

!solv_end
inter -
sele dna end -
sele resn tip3 .or. resn sod end 

calc ener1 = @ener1 + ?ener
calc elec1 = @elec1 + ?elec + ?imel
calc vdw1 = @vdw1 + ?vdw + ?imnb

!solv_strand 1
inter -
sele st1 end -
sele resn tip3 .or. resn sod end

calc ener2 = @ener2 + ?ener
calc elec2 = @elec2 + ?elec + ?imel
calc vdw2 = @vdw2 + ?vdw + ?imnb

!solv_strand 2
inter -
sele st2 end -
sele resn tip3 .or. resn sod end

calc ener3 = @ener3 + ?ener
calc elec3 = @elec3 + ?elec + ?imel
calc vdw3 = @vdw3 + ?vdw + ?imnb

!solv_base 
inter -
sele base end -
sele resn tip3 .or. resn sod end

calc ener4 = @ener4 + ?ener
calc elec4 = @elec4 + ?elec + ?imel
calc vdw4 = @vdw4 + ?vdw + ?imnb

!solv_bkb
inter -
sele bkb end -
sele resn tip3 .or. resn sod end

calc ener5 = @ener5 + ?ener
calc elec5 = @elec5 + ?elec + ?imel
calc vdw5 = @vdw5 + ?vdw + ?imnb

!solv_ribose
inter -
sele ribose end -
sele resn tip3 .or. resn sod end

calc ener6 = @ener6 + ?ener
calc elec6 = @elec6 + ?elec + ?imel
calc vdw6 = @vdw6 + ?vdw + ?imnb

!solv_phosphate
inter -
sele phosphate end -
sele resn tip3 .or. resn sod end

calc ener7 = @ener7 + ?ener
calc elec7 = @elec7 + ?elec + ?imel
calc vdw7 = @vdw7 + ?vdw + ?imnb

! detemine hydration numbers
! full dna
coor stat -
sele type oh2 .and. -
((segid dna1:dna2 .and. .not. hydrogen) .around. 3.50) end

calc hyddna = @hyddna + ?nsel

! full backbone
coor stat -
sele type oh2 .and. -
(bkb .around. 3.50) end

calc hydbkb = @hydbkb + ?nsel

! minor groove excluding terminal base pairs
coor stat -
sele type oh2 .and. -
((minor .and. .not. (ires @first1 .or. ires @last1 .or. ires @first2 .or. ires @last2)) -
 .around. 3.50) end

calc hydmin = @hydmin + ?nsel

! major groove excluding terminal base pairs
coor stat -
sele type oh2 .and. -
((major .and. .not. (ires @first1 .or. ires @last1 .or. ires @first2 .or. ires @last2)) -
 .around. 3.50) end

calc hydmaj = @hydmaj + ?nsel

! phosphate
coor stat -
sele type oh2 .and. -
(phosphate .around. 3.50) end

calc hydpho = @hydpho + ?nsel

! ribose 
coor stat -
sele type oh2 .and. -
(ribose .around. 3.50) end

calc hydrib = @hydrib + ?nsel

! all bases except termini 
coor stat -
sele type oh2 .and. -
(((abase .or. gbase .or. ubase .or. tbase .or. cbase) .and. -
.not. (hydrogen .or. ires @first1 .or. ires @last1 .or. ires @first2 .or. ires @last2)) -
 .around. 3.50) end

calc hydbas = @hydbas + ?nsel

!stop

incr i by 1 

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

! calculate statistics and write data
! total energies

! solv_dna
set tave1 @ener1        !mean
divi tave1 by @nf

! solv_strand 1
set tave2 @ener2        !mean
divi tave2 by @nf

! solv_strand 2
set tave3 @ener3        !mean
divi tave3 by @nf

! solv_base
set tave4 @ener4        !mean
divi tave4 by @nf

! solv_bkb
set tave5 @ener5        !mean
divi tave5 by @nf

! solv_ribose
set tave6 @ener6        !mean
divi tave6 by @nf

! solv_phosphate
set tave7 @ener7        !mean
divi tave7 by @nf

!electrostatic energy

! solv_dna
set eave1 @elec1        !mean
divi eave1 by @nf

! solv_strand1
set eave2 @elec2        !mean
divi eave2 by @nf

! solv_strand2
set eave3 @elec3        !mean
divi eave3 by @nf

! solv_base
set eave4 @elec4        !mean
divi eave4 by @nf

! solv_bkb
set eave5 @elec5        !mean
divi eave5 by @nf

! solv_ribose
set eave6 @elec6        !mean
divi eave6 by @nf

! solv_phosphate
set eave7 @elec7        !mean
divi eave7 by @nf

!vdw energy

! solv_dna
set vave1 @vdw1        !mean
divi vave1 by @nf

! solv_strand 1
set vave2 @vdw2        !mean
divi vave2 by @nf

! solv_strand 2
set vave3 @vdw3        !mean
divi vave3 by @nf

! solv_base
set vave4 @vdw4        !mean
divi vave4 by @nf

! solv_bkb
set vave5 @vdw5        !mean
divi vave5 by @nf

! solv_ribose
set vave6 @vdw6        !mean
divi vave6 by @nf

! solv_phosphate
set vave7 @vdw7        !mean
divi vave7 by @nf

! convert c to distance in A
divi c by 10

write title unit 91
* @runid @tave1 @tave2 @tave3 @tave4 @tave5 @tave6 @tave7
*

write title unit 92
* @runid @eave1 @eave2 @eave3 @eave4 @eave5 @eave6 @eave7
*

write title unit 93
* @runid @vave1 @vave2 @vave3 @vave4 @vave5 @vave6 @vave7
*

calc hyddna = @hyddna/@nf/@natomdna
calc hydbkb = @hydbkd/@nf/@natombkb
calc hydmin = @hydmin/@nf/@natommin
calc hydmaj = @hydmaj/@nf/@natommaj
calc hydpho = @hydpho/@nf/@natompho
calc hydrib = @hydrib/@nf/@natomrib
calc hydbas = @hyddna/@nf/@natombas

write title unit 94
* @runid @hyddna @hydbas @hydmin @hydmaj
*

write title unit 95
* @runid @hydbkb @hydpho @hydrib
*

stop