* generate 5'-cgcgaattcgcg-3' dna duplex and overlay with
* a waterbox that includes sodiums.  number of sodiums
* the optimized to yield a neutral system.  
*

!this version translate and rotates the DNA so it is
!protruding from the waterbox, allowing for a test
!of image generation by CRYSTAL

!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 60.
set ydim 36.
set zdim 36.

! identifiers
set id1 ecor_sod
set id2 test

! 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 14.0 ! cutim
set 4 14.0 ! cutnb
set 5 10.0 ! ctonnb
set 6 12.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

! generate segid dna1
read sequence card
* cgcgaattcgcg
*
12
cyt gua cyt gua ade ade thy thy cyt gua cyt gua

generate dna1 first 5ter last 3ter setup warn

!set variable for later use
set nstrand1 = ?nres

! generate segid dna2
read sequence card
* cgcgaattcgcg
*
12
cyt gua cyt gua ade ade thy thy cyt gua cyt gua

generate dna2 first 5ter last 3ter setup warn

!set variable for later use
calc nstrand2 = ?nres - @nstrand1

!change ribose to deoxy
!patch for deoxy, 1 for pyrimidine, 2 for purine
!first strand
set nbase = 1
label patchsegid1
  set pres = deo1
   define thisresidue select segid dna1 .and. resi @nbase end
    if ?selresn eq GUA set pres = deo2
    if ?selresn eq ADE set pres = deo2
   patch @pres dna1 @nbase setup warn
  calc nbase = @nbase + 1
if nbase le @nstrand1 goto patchsegid1

!second strand
set nbase = 1
label patchsegid2
  set pres = deo1
   define thisresidue select segid dna2 .and. resi @nbase end
    if ?selresn eq GUA set pres = deo2
    if ?selresn eq ADE set pres = deo2
   patch @pres dna2 @nbase setup warn
  calc nbase = @nbase + 1
if nbase le @nstrand2 goto patchsegid2

autogenerate angle dihedral

open unit 20 read form name deo_1a.crd
read coor card unit 20

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

!translate DNA to test image building when DNA is
!partially out of the box
coor trans xdir 0.0 ydir 1.0 zdir 0.0 dist 10.0

!rotate DNA about long (x-axis)
coor axis sele atom dna1 1 c1' show end -
          sele atom dna2 1 c1' show end
coor rotate axis phi 90. sele all end

! read pre-equilibrated sodium/water box
open unit 30 read form name wat_sod_30_18_18.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 dna1:dna2) .and. .not. hydr) .around. 1.6)) 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 dna1 .or. segid dna2 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.1

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 DNA 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 dna
! 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
 
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

! 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

!turn on image centering of solvent
image byres xcen 0.0 ycen 0.0 zcen 0.0 -
     sele resn tip3 .or. resn sod .or. resn cl end

!kappa should be optimized based on the value of
!kappa at which the 1st derivative of the energy
!with respect to kappa goes to 0.0
set kap 0.36
set ord 6

!fftx dimensions should correspond to 1 A or less; must
!be power of 2, 3 or 5
update inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa @kap fftx 64 ffty 32 fftz 32 order @ord -
@7 @8 @9 vswitch cutimg @3 cutnb @4 ctofnb @6 ctonnb @5 imall


open unit 20 write form name @id1_@id2_rotate.crd
write coor cards unit 20 
* 5'-CGCGAATTCGCG-3'  generated coordinates, all-hydrogen
*    ||||||||||||     in a 60x36x36 cube of sodium/water
* 3'-GCGCTTAAGCGC-5'  1.6 A atom delete
* with images
*


open unit 20 write form name @id1_@id2_rotate_image.crd
write coor cards unit 20 image
* 5'-CGCGAATTCGCG-3'  generated coordinates, all-hydrogen
*    ||||||||||||     in a 60x36x36 cube of sodium/water
* 3'-GCGCTTAAGCGC-5'  1.6 A atom delete
* with images
*


stop