*  ------------------------------
*      CHARMM setup script
*  ------------------------------
*


set sim aa
set ff all36_na
set ffname c36
set input "aa_charmm.pdb"
set inputions "aa_ions.pdb"
set outvac "aa-c36-vac.pdb"
set outcio "aa-c36-cio.pdb"

set comment "Initial ApA model"

!
! number of waters and guess on initial box size
!
set waters = 1075
set box = 36.00


set finish = 0

!
!  we need to build a box bigger than we need and cut it down to
!  appropriate size.  This is the "pad".
!

set pad = 13.0
set solutecut 2.5

bomblev -3

!
!  read in force field and topology information
!
open unit 1 read card name top_@ff.rtf
read rtf card unit 1
close unit 1

open unit 1 read card name par_@ff.prm
read param card unit 1 flex
close unit 1

!
!  create rna ("rna")
!

read sequence card
2
ADE ADE
generate rna first 5ter last 3ter setup warn

!
!  CHARMM by default builds RNA; to get DNA you patch or modify the
!  RNA using the commands below...  Commented out at present.
!
!patch dog5 rna 1
!patch doa  rna 2
!patch doc  rna 3
!patch doc3 rna 4

open read card unit 1 name @input
read coordinates pdb unit 1 select segid rna end
close unit 1

ic param
ic build

hbuild
autogenerate angles dihedrals

coor orient mass select all end
coor stat sele all end

!
!  Write out the in-vacuo files needed
!
open unit 1 write card name @sim_@ffname_vac.psf
write psf unit 1 card
* @sim @comment
*
close unit 1

open unit 1 write card name @sim_@ffname_vac.pdb
coor write unit 1 pdb select all end
* @sim @comment
*
close unit 1

open unit 1 write card name @sim_@ffname_vac.cor
coor write unit 1 card select all end
* @sim @comment
*
close unit 1



!
!  create Na+ counterions ("ions")
!


read sequence card
11
NA   NA   CLA  NA   CLA  NA   CLA  NA   CLA  NA   CLA
generate ions setup first none last none

open read card unit 1 name @inputions
read coordinates pdb unit 1 select segid ions end
close unit 1



coordinate orient mass select all end
coordinate statistics select all end


!
!  MINIMIZE ION POSITIONS
!
!  ...change the ion "type" from IP (31) to IB (32)
!
!scalar type set 32 select segid ions end
!
!cons fix select segid rna .and. .not. type H* end
!
!mini abnr nstep 1000
!
!  ...change the ion "type" back!
!
!scalar type set 31 select segid ions end
!
!open unit 1 write card name @outcio
!write coor pdb unit 1 select all end
!* B-DNA @sim
!* ion positions gas-phase minimized
!*
!close unit 1
!
!cons fix select none end


open unit 1 write card name @sim_@ffname_cio.psf
write psf unit 1 card
* @sim @comment
* 1 Na+ + 5 NaCl ions
*
close unit 1

open unit 1 write card name @sim_@ffname_cio.pdb
coor write unit 1 pdb select all end
* @sim @comment
* 1 Na+ + 5 NaCl ions
*
close unit 1

open unit 1 write card name @sim_@ffname_cio.cor
coor write unit 1 card select all end
* @sim @comment
* 3 Na+ ions
*
close unit 1


define SOLUTE select .not. segid wat end


set boxlen 31.1661028
read sequence wat 1000
generate wat noangle nodihedral

open read formatted unit 1 name wat1000.cor.tip3p
read coordinates card unit 1 append select segid wat end
close unit 1

!
!  Make the box pretty big for now...
!
generate watx duplicate wat
coor duplicate select segid wat end select segid watx end
coor translate xdir @boxlen select segid watx end
join wat watx renumber

generate waty duplicate wat
coor duplicate select segid wat end select segid waty end
coor translate ydir @boxlen select segid waty end
join wat waty renumber

generate watz duplicate wat
coor duplicate select segid wat end select segid watz end
coor translate zdir @boxlen select segid watz end
join wat watz renumber

define junk select .not. init show end

coordinate orient mass norotate select segid wat end
coordinate statistics select all end

calc xbox ?XMAX - ?XMIN
calc ybox ?YMAX - ?YMIN
calc zbox ?ZMAX - ?ZMIN

set boxmax
let boxmax = max @xbox @ybox
let boxmax = max @boxmax @zbox

calc cutimg = 1.0

set boxincr 0.002
calc halfbox @box / 2.0 + 2.0

!if box .gt. @maxbox stop

delete atom select .byres. (atom WAT * O .and. (property x .lt. -@halfbox )) end
delete atom select .byres. (atom WAT * O .and. (property x .gt.  @halfbox )) end
delete atom select .byres. (atom WAT * O .and. (property y .lt. -@halfbox )) end
delete atom select .byres. (atom WAT * O .and. (property y .gt.  @halfbox )) end
delete atom select .byres. (atom WAT * O .and. (property z .lt. -@halfbox )) end
delete atom select .byres. (atom WAT * O .and. (property z .gt.  @halfbox )) end

join wat renumber

label loop2


set angle = 109.4712206344907

crystal define octahedral @box @box @box @angle @angle @angle
crystal build cutoff 14.0

coordinate copy comp

image byres select segid wat end
update inbfrq 0 ihbfrq 0

coordinate diff comp
coordinate dist weigh comp

delete atom select .byres. (atom WAT * O .and. (property wcomp .gt. 0.001 )) end

join wat renumber

crystal free
incr box by -@boxincr

!
!  fudge factor for waters
!    set the # value to be total of non-water residues in the system
!
calc numwat ?NRES - 13
if @waters .lt. @numwat goto loop2
incr box by @boxincr



crystal define octa @box @box @box @angle @angle @angle
crystal build cutoff @cutimg

image byresidue select .not. SOLUTE end
coor copy comp
update cutim 1.0
coor diff comp
coor dist comp weig
delete atom sele .byres. prop wcomp .gt. 0.0 end




open unit 1 write card name @sim_@ffname.psf
write psf unit 1 card
* @sim @comment
* ?XTLA ?XTLB ?XTLC ?XTLALPHA ?XTLGAMMA ?XTLBETA
* 3 Na+ ions + @numwat waters
*
close unit 1

open unit 1 write card name @sim_@ffname.cor
coor write unit 1 card select all end
* @sim @comment
* ?XTLA ?XTLB ?XTLC ?XTLALPHA ?XTLGAMMA ?XTLBETA
* 3 Na+ ions + @numwat waters
*
close unit 1

open unit 1 write card name @sim_@ffname.pdb
coor write unit 1 pdb select all end
* @sim @comment
* ?XTLA ?XTLB ?XTLC ?XTLALPHA ?XTLGAMMA ?XTLBETA
* 3 Na+ ions + @numwat waters
*
close unit 1