*  Generate coordinates of a "box" of TIP3P water
*  of specified dimensions based on tip32.crd
*  The resulting waterbox: watertmp.crd is then
*  trimmed to obtain a waterbox of the desired size.  
*  Note that the dimensions of watertmp.crd must
*  be significantly larger than the final waterbox
*  to avoid problems with the "rounded" corners in
*  waterbox.crd.   
*  May also be used to generate a waterball (see below)
*  adm jr.
*
        
! all images with one or more atoms within the following distance
! of the primary atoms will be generated.  This is what controls
! the temporary box size.  Note that this value should be
! ca. one half that of the maximum dimension of the final box
set q 28.0
set out output

set top $CHARMMDATA

!set box dimensions: set following values to half of
!the final box dimensions
set x 28
set y 21
set z 19

!set variables to create waterbox without ions
!set isod 0
!set id wat

!uncomment the following to create waterbox that
!contains sodium
set isod 1
set id wat_sod

! 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 @q ! cutim
set 4 @q ! cutnb
set 5 @q ! ctonnb
decr 5 by 2
set 6 @q ! ctofnb
decr 6 by 1
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

if isod eq 1 goto sod_box

open unit 20 read card name TIP32.CRD
read sequence coor unit 20

generate SOLV setup warn noangle nodihedral

rewind unit 20
read coor cards unit 20

! create images via CRYSTAL
! the unit cell parameters correspond to the TIP32.CRD waterbox
crystal define cubic 18.856 18.856 18.856 90.00 90.00 90.00

! build the crystal including all images within the cutoff
crystal build Noper 0 cutoff @3

goto wat_box

!section to create waterbox that includes sodium ions
label sod_box


!read coordinates and generate
open unit 20 read card name wat_sod_1b.crd
read sequence coor unit 20

generate SOLV setup warn noangle nodihedral

rewind unit 20
read coor cards unit 20

! create images via CRYSTAL
! the unit cell parameters correspond to the wat_sod_1b.crd waterbox
crystal define cubic 18.65347 18.65347 18.65347 90.00 90.00 90.00

! build the crystal including all images within the cutoff
crystal build Noper 0 cutoff @3

label wat_box

! update required to build the images
! IMALL specifies that all images be built, not
! just the unique images
update ihbfr 0 inbfr 1 cutnb @3 cutim @3 IMALL

COOR STAT image sele all end

!write coordinates, including image atoms, to a temporary
!file to be read later by the script
open unit 20 write form name watertmp.crd 
write coor images cards unit 20 
close unit 20

!read temporary coordinate file and trim waterbox to appropriate size

!must delete all atoms prior to generating new system
delete atom sele all end
bomlev -2

!the following error is assocaited with format
!problems and can (usually) be ignored
!
!      ***** LEVEL -1 WARNING FROM <ENCODI> *****
!      ***** STRING TOO SMALL
!      ******************************************
!      BOMLEV ( -2) IS NOT REACHED. WRNLEV IS  5
!
wrnlev -2

open unit 20 read form name watertmp.crd
read sequence coor unit 20

generate SOLV setup warn noangle nodihedral

rewind unit 20
read coor cards unit 20

coor stat sele .not. hydrogen end

!uncomment following to create a water sphere
!goto shpere

scalar wmain = x
delete atom sele .not. .byres. -
((type oh2 .or. type sod) .and. -
(prop abs wmain .lt. @x)) end

scalar wmain = y
delete atom sele .not. .byres. -
((type oh2 .or. type sod) .and. -
(prop abs wmain .lt. @y)) end

scalar wmain = z
delete atom sele .not.  .byres. -
((type oh2 .or. type sod) .and. -
(prop abs wmain .lt. @z)) end

! the following commands may be used to offset the
! center of the box
!coor trans xdir -4.0 ydir 0.0 zdir 0.0
coor stat sele .not. hydrogen end

open unit 20 write card name @out/@id_@x_@y_@z.crd
write coor cards unit 20
* Waterbox of dimensions @xx@yx@z A cube prepared from a
* prequilibrated waterbox or sodium/waterbox
* ?nsel waters/sodiums
*

open unit 20 write card name @out/@id_@x_@y_@z.pdb
write coor pdb unit 20
* Waterbox of dimensions @xx@yx@z A cube prepared from a
* prequilibrated waterbox or sodium/waterbox
* ?nsel waters/sodiums
*

!delete temporary waterbox, this command fails on some systems
system "rm watertmp.crd"

stop

label sphere

! to prepare a sphere of water centered around the origin
set j 18  ! radius of sphere
delete atom sele .not. (.byres. (point 0.0 0.0 0.0 cut @j .and. type oh2)) end

coor stat
coor stat sele .not. hydrogen end

open unit 20 write card name waterball_@j.crd  
write coor cards unit 20
* @j A sphere of water prepared from a
* prequilibrated waterbox or sodium/waterbox
* ?nsel waters/sodiums
*

!delete temporary waterbox, this command fails on some systems
system "rm watertmp.crd"

stop