* OVERLAY RNA WITH A BOX CONTAINING
* WATER AND SODIUM IONS
* THE RESULTANT SYSTEM IS NEUTRAL!!!
*

!Input data
!1) box dimensions, xdim, ydim, zdim
!2) identifiers, id1, id2
!3) topology filename
!4) parameter filename
!5) DNA coordinate filename (and the corresponding sequence
!   and DNA patch information)
!6) waterbox coordinate filename
bomlev -1

!set box dimension
set xdim 28
set ydim 21
set zdim 19

! identifiers
set id1 rr_gen
set id2 rr_sol
set top $CHARMMDATA
set out output

! 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/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 shift
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 20 read form name @out/@id1.psf
read psf cards unit 20

open unit 20 read form name @out/@id1.crd
read coor card unit 20

! info on size of structure; use this information when
! creating waterbox
coor stat
coor orient sele segid rna1:rna2 .and. .not. hydrogen show end
coor stat

! read pre-equilibrated sodium/water box
open unit 30 read form name @out/wat_sod_@xdim_@ydim_@zdim.crd
!open unit 30 read form name wat_30_18_18.crd
read sequence coor unit 30

generate solv first none last none noangle nodihedral

rewind unit 30
read coor card unit 30 append

! delete waters overlapping with the solute
delete atom sele .byres. (segid solv .and. (.not. hydrogen) .and. -
(((segid rna1 .or. segid rna2) .and. .not. hydr) .around. 1.8)) end

coor stat sele .not. hydrogen end
coor stat sele resn tip3 end
coor stat sele resn sod end
set nsod ?nsel
!coor stat sele resn mg show end
!coor stat sele resn cl end
coor stat sele type p* end
set npho ?nsel
coor stat sele type n* end
coor stat sele type h* end
coor stat sele type o* end
coor stat sele type c* .and. .not. resn cl end

!define term to specify selection of all solute atoms
define solute sele segid rna1 .or. segid rna2 end

!determine if sodiums need to be deleted or added
set charge ?cgtot
if charge eq 0 goto write_crd
if charge lt 0 goto sod_add

!if charge gt 0.0 delete sodiums to create neutral system

!calculate final number of sodiums required for neutral system
calc nsodfinal = @nsod - @charge

set dist 6.0

label sod_del_loop

coor stat sele resn sod .and. -
(solute .around. @dist) end

! this statement requires that each phosphorous
! atom correspond to a -1 unit charge
if nsodfinal eq ?nsel goto del_sod

calc dist = @dist + 0.01

if dist le 15.0 goto sod_del_loop

! exit if proper number of sodiums to neutralize
! system is not found. when this occurs narrow
! the range of distances and decrease dist
! increment to 0.01

stop

label del_sod

!delete the sodiums furthest from the DNA
!
delete atom sele resn sod .and. -
.not. (solute .around. @dist) end

coor stat sele resn sod end
coor stat sele type p* end

goto write_crd

label sod_add

calc nadd = abs ( @charge )

! add sodiums to make the system neutral
! generate the sodiums
! 
read sequence sod 1
generate sod noangle nodihedral

read coordinate card append
* sodium starting coordinates 
*
    1
    1    1 SOD  SOD    0.00000   0.00000   0.00000 SOD  1      0.00000

!use replica to create required number of sodiums
!
calc naddm1 = @nadd - 1
replica A nreplica @naddm1 -
  sele segid sod end

set count 1
label join_loop

join sod A@count renumber 

calc count = @count + 1

if count le @naddm1 goto join_loop

coor print sele resn sod end

! loop over the sodiums, place in random positions in the
! box, checking to avoid overlap with the RNA and close
! overlaps with water

set r 1
set i 87654

random uniform scale 1. offset 0.0 iseed @i

label loop

goto skip_redo

label redo_resi
! place ion back at the origin and generate new coordinates
coor trans xdir -@xcoor ydir -@ycoor zdir -@zcoor sele segid sod .and. resi @r end

label skip_redo
! x coordinate
set xcoor ?rand
calc xcoor = ( @xcoor * @xdim ) - ( @xdim / 2.0 )
! y coordinate
set ycoor ?rand
calc ycoor = ( @ycoor * @ydim ) - ( @ydim / 2.0 )
! z coordinate
set zcoor ?rand
calc zcoor = ( @zcoor * @zdim ) - ( @zdim / 2.0 )

coor trans xdir @xcoor ydir @ycoor zdir @zcoor sele segid sod .and. resi @r end

! check for close interaction with the tetraplex
! 3.0 selected to avoid first solvation shell
coor mindist -
 sele segid sod .and. resi @r end -
 sele solute .and. .not. hydrogen end
if ?mind lt 3.0 goto redo_resi

! check for extremely close interaction with water oxygen
!
coor mindist -
 sele segid sod .and. resi @r end -
 sele type oh2 end
if ?mind lt 0.5 goto redo_resi

! check for close interaction with other sodiums
!
coor mindist -
 sele segid sod .and. resi @r end -
 sele segid sod .and. .not. resi @r end
if ?mind lt 3.0 goto redo_resi

incr r by 1
if r le @nadd goto loop

coor print sele segid sod end

label skip_ion_add

label write_crd

!write coordinates, PSF and IC tables
 
open unit 20 write form name @out/@id2.crd
write coor cards unit 20
* @1
* @2
*

open unit 20 write card name @out/@id2.pdb
write coor pdb unit 20
* @1
* @2
*


open unit 20 write form name @out/@id2.psf
write psf cards unit 20
* @1
* @2
*

!check total charge
set charge ?cgtot

coor stat sele .not. hydrogen end
coor stat sele resn tip3 end
coor stat sele resn sod end
!coor stat sele resn cl end
!coor stat sele resn mg end
coor stat sele type p* end
coor stat sele type n* end
coor stat sele type h* end
coor stat sele type o* end
coor stat sele type c* .and. .not. resn cl end
coor stat sele segid rna1:rna2 end
coor stat sele all end

stop