# author: mengen

#

# usage: ./modar step5_gen_large_box.inp npair=? nstepmini=1000

#     or ./modar step5_gen_large_box.inp boxtype=cubic A=? B=? C=? nstepmini=1000

#

#load force field

include "loadff.inc"

anion="bmi"

cation="pf6"

if(! -var nstepmini) nstepmini=200

if(!-var boxtype && !-var npair) {

  :retry

  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        "

  gettext ichoice default=1 min=1 max=5 \

                  msg="please make a choice:"

  if($ichoice=="1") boxtype="cubic"

  else if($ichoice=="2") boxtype="tetragonal"

  else if($ichoice=="3") boxtype="rectangular"

  else if($ichoice=="4") boxtype="orthorhombic"

  else if($ichoice=="5") boxtype="octahidralt"

  else goto retry

}

#load raw solvent box unit prepared

loadmsf file="step4_equil_npt.msf"

tell box

tell density

lunit=$boxA

x=$lunit/2

atom trans vector=($x,$x,$x) select="all"

#get box size if no npair specified

if(!-var npair) {

 if(!-var a) getvar a default=48 min=20 max=1000 \

                     msg="please tell box leng a"

 if(!-var b) getvar b default=48 min=20 max=1000 \

                     msg="please tell box leng b"

 if(!-var c) getvar c default=48 min=20 max=1000 \

                     msg="please tell box leng c"

}

#compute box size for number pairs case

if(-var npair) {

 tell msf_info select="all"

 tell density  select="all"

 boxtype="cubic"

 v=$volume*$npair*2/$nres

 a=pow($v,1/3)

 b=$a

 c=$a

 echo -hs " number of pairs: $npair"

 echo -hs " the box will be cubic with size:"

 echo -hs "      $a x $a x $a"

}

#process boxtype

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

  boxl=$a

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

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

  a=$boxl

  b=$boxl

  c=$boxl

}

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

  boxl=$a

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

  a=$boxl

  b=$boxl

}

else if($boxtype .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 boxtype: $boxtype"

  goto tellboxtype

}

if(!-var npair) {

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

  tell box

  tell msf_info

  npair=int($nres/2*$boxvol/$lunit/$lunit/$lunit+0.5)

}

#compute number of replication for each dimmension

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

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

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

#do replication

echo -hs "| building box ..."

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

#move center to origin

cx=$a/2

cy=$b/2

cz=$c/2

atom trans vector=(-$cx,-$cy,-$cz) select=all

#setup crystal/pbc box

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

tell boxsize

echo "crystal type=$boxtype 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 anions outside the box

tell nres select="resn $anion"

while($nres>$npair) {

  tell faroutboxres select="resn $anion"

  if(-z $resexpr) tell farres select="resn $anion"

  delete atom select="residue $resexpr"

  tell nres select="resn $anion"

}

#remove cations outside the box

tell nres select="resn $cation"

while($nres>$npair) {

  tell faroutboxres select="resn $cation"

  if(-z $resexpr) tell farres select="resn $cation"

  delete atom select="residue $resexpr"

  tell nres select="resn $cation"

}

#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=10.0 nbcutoff=8.0 swcutoff=7.0 eps=1.0 \

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

#do minimization

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

nbsoft 1.0

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

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

nbsoft 0

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

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

#save msf and pdb

savemsf file="step5_largebox.msf"

savepdb fmt=pdb file="step5_largebox.pdb"

echo -hs "| "

echo -hs "|  output files:"

echo -hs "|    step5_largebox.msf"

echo -hs "|    step5_largebox.pdb"

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

echo -hs "| "