if(! -var nstep) getvar nstep default=50000 min=1000 \

                              msg="How may steps to run"

#load force field

include "loadff.inc"

#load MSF built previously

loadmsf file="step7_two_phases.msf"

#setup crystal/pbc box

include "pbc.inc"

#image groups

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

#constraint

shake type=bondH tol=1e-6 scale=1.0 maxite=500 ref=para

#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

#restraint for the boundary

restrain zrange applyto=residue kforce=5 startpos=-24 endpos=24 select="segn ILC"

restrain zblock applyto=residue kforce=5 startpos=-20 endpos=20 select="resn TIP*"

#langevin-piston NPT coupling

tell tmass #print total mass, and put total mass number to $totmass

set pmass = int($totmass/50.0)

ensemble npt pmass=$pmass pgamma=20.0 \

         pxx=1.0 pyy=1.0 pzz=1.0 \

         tbath=300.0 tref=300.0 tmass=2000

#build MD core

build mdcore timestep=0.002 nprint=250 \

     trjfile="step8_equil_npt.trj" frqtrj=250 \

     enefile="step8_equil_npt.ene" frqene=10 \

     rstfile="step8_equil_npt.rst" frqrst=100000 \

     flushene=true \

     randseed=314159 randtype=old \

     frqresort=20000

#save MD core

savemds file="step8_equil_npt.mds"

if($nstep==0){

  echo -hs "please run MD by command"

  echo -hs " ./modar step8_equil_npt.mds nstepmore=50000"

  stop

}

#run MD

runmd nstep=$nstep

restrain clear

savemsf file="step8_equil_npt.msf"