# generate a water box

# author: mengen

#

# usage: ./modar step2_gen_solvbox.inp type=cubic A=? B=? C=? nstepmini=1000

#        type must be presented, otherwise job will be aborted

#        A,B,C if not specified, will be assigned according to the size of solute

#        and do nstepmini minimization for the solvent box, default value=100

#

#load solute_size.inc genetated by step1_build_solute.inp

include "solute_size.inc"

if(! -var nstepmini) nstepmini=100

if(!-var type) {

  :retrytype

  type="cubic"

  echo -hs " no box type specified"

  echo -hs " available box types:"

  echo -hs "   1 for cubic           a=b=c        "

  echo -hs "   2 for tetragonal      a=b          "

  echo -hs "   3 for rectangular     a/A=b/B=c/C  "

  echo -hs "   4 for orthorhombic    a,b,c        "

  echo -hs "   5 for octahidralt     a=b=c        "

  echo -hsn "please make a choice in above:"

  input=gettext()

  if($input=="1") type="cubic"

  else if($input=="2") type="tetragonal"

  else if($input=="3") type="rectangular"

  else if($input=="4") type="orthorhombic"

  else if($input=="5") type="octahidralt"

  else                 goto retrytype

}

#assign A, B, C if not assigned

aa=48

ba=48

ca=48

if(-var xext) aa=int($xext)+2*10+1

if(-var yext) ba=int($yext)+2*10+1

if(-var zext) ca=int($zext)+2*10+1

if(! -var a) {

  gettext a default=$aa min=$aa max=1000 \

            msg="no box length A specified"

  if(! -var b) {

    gettext b default=$ba min=$ba max=1000 \

            msg="no box length B specified"

  }

  if(! -var c) {

    gettext c default=$ca min=$ca max=1000 \

            msg="no box length C specified"

  }

}

if(!-var b) b=$a

if(!-var c) c=$a

if($a<1) $a=$aa

if($b<1) $b=$ba

if($c<1) $c=$ca

#process type

if($type .head. "cubi") { #cubic

  boxl=$a

  if($boxl<$b) boxl=$b

  if($boxl<$c) boxl=$c

  a=$boxl

  b=$boxl

  c=$boxl

}

else if($type .head. "tetr") { #tetragonal

  boxl=$a

  if($boxl<$b) boxl=$b

  a=$boxl

  b=$boxl

}

else if($type .head. "octa") { #octahidralt

  boxl=$a

  if($boxl<$b) boxl=$b

  if($boxl<$c) boxl=$c

  a=$boxl

  b=$boxl

  c=$boxl

}

else{

  echo -hs "|  ** unknown type: $type"

  goto telltype

}

#load force field

include "loadff.inc"

#load raw solvent box unit prepared

loadmsf file="waterbox.msf"

tell box

lunit=$boxA

#compute number of replication for each dimmension

nx=int(($a+8)/$lunit+1)

ny=int(($b+8)/$lunit+1)

nz=int(($c+8)/$lunit+1)

#do replication

echo -hs "| building box ..."

replicate nx=$nx ny=$ny nz=$nz stride=($lunit,$lunit,$lunit)

#move center to origin

atom transto pos=(0,0,0) select=all

tell center select=all

#setup crystal/pbc box

crystal type=$type a=$a b=$b c=$c

tell boxsize

echo "crystal type=$type a=$a b=$b c=$c" > pbc.inc

tell bestfft

echo "fftx=$fftx" >> pbc.inc

echo "ffty=$ffty" >> pbc.inc

echo "fftz=$fftz" >> pbc.inc

#remove residues outside the box

removeresidueoutsidebox ext=4 select="all"

#remove residues overlapped

resolveresidueoverlap cutoff=2.2 select="name OH2"

#resolveresidueoverlap cutoff=1.0 select="all"

#check bad overlap

tell contactsamong cutoff=1.5 select="name O*"

#checkmsf

checkmsf

tell msf_info

tell density

echo -hs "| "

echo -hs "|  Summary of solvent box:"

echo -hs "|    box type: $boxtype"

echo -hs "|    size in angstrom: $a x $b x $c"

echo -hs "|    number of residue: $nres"

echo -hs "|    number of atoms: $natom"

echo -hs "|    density: $density"

echo -hs "|    dimmension of FFT suggested: $fftx x $ffty x $fftz"

#resort sequence

atom resortseq  start=1 select=all

tell center select=all

tell size   select=all

#image groups

imggroup method=resid sorted=true select="all"

#nonbond setup

nonbond type=pme nblcutoff=14.0 nbcutoff=10.0 swcutoff=8.0 eps=1.0 \

        beta=0.34 ftx=$fftx fty=$ffty ftz=$fftz bsorder=6

#do minimization for solvent

echo -hs "|  running $nstepmini steps minimization for the solvent box, please wait a moment"

minimize type=SDfuse nstep=$nstepmini nprint=10 ftol=1e-5

minimize type=SD     nstep=$nstepmini nprint=10 ftol=1e-5

tell pot

echo -hs "| total potential energy: $pot kcal/mol"

#save msf and pdb

savemsf file="step2_solvbox.msf"

savepdb fmt=pdb file="step2_solvbox.pdb"

echo -hs "| "

echo -hs "|  output files:"

echo -hs "|    step2_solvbox.msf"

echo -hs "|    step2_solvbox.pdb"

echo -hs "|    pbc.inc           the crystal unit type and size"

echo -hs "| "

echo -hs "|  it's ready to go next step to add countour ions by command:"

echo -hs "|    ./modar -nobak step3_add_ions.inp"

echo -hs "| "