#get density and box size from user if not specified

if(! -var msf)     getvar msf checkfile=true \

                       msg="Please tell the msf file of basic unit"

if(! -var density) getvar density msg="no density specified" \

                                  default=1.0 min=1e-6 max=100

if(! -var a)       getvar a msg="no box length specified"\

                            default=32 min=18 max=1000

#load force field

include "loadff.inc"

#load MSF built previously

loadmsf file="$msf"

#reorient molecule and move center to origin

atom reorient select="all"

#get atoms count for the single unit

tell natom select="all"

natomu=$natom

#compute the unit extension, assume it's in cubic shape

tell mass select="all"

vol=$totmass/0.602214199/$density

lunit=pow($vol,1/3)

#compute number of replication for each dimmension

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

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

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

#do replication

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

#make random orient for each single unit

if(-var rotate) {

 iu=0

 nu=$nx*$ny*$nz

 while($iu<$nu) {

   ias=$iu*$natomu

   iae=$ias+$natomu-1

   tell center select="atomidx $ias:$iae"

   xv=rand(1.0)

   yv=rand(1.0)

   zv=sqrt(1-$xv*$xv-$yv*$yv)

   th=rand(360)

   atom rotate origin=($xce,$yce,$zce) vector=($xv,$yv,$zv) theta=$th \

        select="atomidx $ias:$iae"

   iu=$iu+1

 }

}

#the actual box size

a=$nx*$lunit

#move center to origin

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

#setup crystal/pbc box

pbcbox type=cubic a=$a

#save PBC setup for further simulation

tell box

echo "crystal type=$boxtype a=$boxa" > pbc.inc

tell bestfft

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

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

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

#image groups

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

#nonbond setup

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

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

tell msf_info

echo -hs "number of residues: $nres"

echo -hs "number of atoms: $natom"

tell box

echo -hs "box type: $boxtype"

echo -hs "box size: $boxa x $boxb x $boxc"

tell density

echo -hs "density: $density"

#do minimization for solvent

echo -hs "running minimization ..."

minimize type=sdfuse nstep=200 nprint=10 ftol=1e-5 \

         frqene=1 enefile=step2_solvbox.ene

minimize type=mc     nstep=200 nprint=10 ftol=1e-5 \

         frqene=1 enefile=step2_solvbox.ene

minimize type=sd     nstep=200 nprint=10 ftol=1e-5 \

         frqene=1 enefile=step2_solvbox.ene

#save msf and crd

savemsf file="step2_solvbox.msf"

savecrd fmt=pdb file="step2_solvbox.pdb"

tell density

tell msf_info