* calculate base pair interactions between the bases
*

faster on

! identifiers
set id1 rr_sol
set id2 rr
set top $CHARMMDATA 
set psfcrd ../output !directory where psf and crd files exist
set trjdir ../output ! directory for the trajectories

! 1: topology
! 2: parameter
set 1 @top/top_all27_na.rtf
set 2 @top/par_all27_na.prm
set 3 12.0 ! cutim
set 4 12.0 ! cutnb
set 5 8.0 ! ctonnb
set 6 10.0 ! ctofnb
set 7 fshift
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 10 read form name @psfcrd/@id1.psf
read psf cards unit 10

open unit 10 read form name @psfcrd/@id1.crd
read coor cards unit 10
coor copy comp

open unit 11 write form name @id2_a2u9.inter
open unit 12 write form name @id2_a3u8.inter
open unit 13 write form name @id2_g4c7.inter
open unit 14 write form name @id2_a5u6.inter
open unit 15 write form name @id2_g6c5.inter
open unit 16 write form name @id2_a7u4.inter
open unit 17 write form name @id2_a8u3.inter
open unit 18 write form name @id2_g9c2.inter

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

update inbfrq 1 imgfrq 1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5 !bycube

set t 1.0  ! time

set begin 100
set stop 5000
set skip 100
calc time = @begin * 0.002
calc stime = @skip * 0.002
calc time_max = @nstep * 0.002

set trajs 1
set traji 1
set trajf 3 
set unit 21
set trajn 0

label trj_loop

open unit @unit unform read name @trjdir/@id2_prod_@trajs.trj

calc trajs = @trajs + @traji
calc trajn = @trajn + @traji
calc unit = @unit + 1

if @trajs le @trajf goto trj_loop

traj iread 21 nread @trajn nskip @skip begin @begin stop @stop

label rms_loop

traj read

!A2U9
inter sele segid RNA1 .and. resid 2 .and. abase show end -
      sele segid RNA2 .and. resid 9 .and. ubase show end

write title unit 11
* @time ?ener ?elec ?vdw
*

!A3U8
inter sele segid RNA1 .and. resid 3 .and. abase show end -
      sele segid RNA2 .and. resid 8 .and. ubase show end

write title unit 12
* @time ?ener ?elec ?vdw
*

!G4C7
inter sele segid RNA1 .and. resid 4 .and. gbase show end -
      sele segid RNA2 .and. resid 7 .and. cbase show end

write title unit 13
* @time ?ener ?elec ?vdw
*

!A5U6
inter sele segid RNA1 .and. resid 5 .and. abase show end -
      sele segid RNA2 .and. resid 6 .and. ubase show end

write title unit 14
* @time ?ener ?elec ?vdw
*

!G6C5
inter sele segid RNA1 .and. resid 6 .and. gbase show end -
      sele segid RNA2 .and. resid 5 .and. cbase show end

write title unit 15
* @time ?ener ?elec ?vdw
*

!A7U4
inter sele segid RNA1 .and. resid 7 .and. abase show end -
      sele segid RNA2 .and. resid 4 .and. ubase show end

write title unit 16
* @time ?ener ?elec ?vdw
*

!A8U3
inter sele segid RNA1 .and. resid 8 .and. abase show end -
      sele segid RNA2 .and. resid 3 .and. ubase show end

write title unit 17
* @time ?ener ?elec ?vdw
*

!G9C2
inter sele segid RNA1 .and. resid 9 .and. gbase show end -
      sele segid RNA2 .and. resid 2 .and. cbase show end

write title unit 18
* @time ?ener ?elec ?vdw
*

calc time = @time + @stime
if time le @time_max goto rms_loop

stop