# do NPT equilibration

# author: mengen

#

# usage: ./modar step7_equil_npt.inp nstep=10000

#        do nstep NPT equilibration

#        50ps is enough

#    

if(! -var nstep) {

  nstep=5000

  echo -hs " no nstep specified"

  echo -hs "   $nstep steps will be"

  echo -hsn "   or tell it now:"

  n=gettext()

  if($n>1000 && $n<100000) nstep=$n

}

#load force field

include "loadff.inc"

#load md system

loadmsf file="step6_heatup.msf"

tell msf_info

#load box setup

include "pbc.inc"

#shake setup

shake enable=true 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

#load NPT coupling method

include "nptc.inc"

#run MD

build mdcore timestep=0.002 nstep=100 nprint=10 \

     trjfile="step8_equil_npt.trj" frqtrj=100 trjfmt=movar \

     rstfile="step8_equil_npt.rst" frqrst=100 \

     enefile="step8_equil_npt.ene" frqene=10 \

     randseed=314159 randtype=old \

     frqresort=20000

#load last coordinates and velocities

loadstate file="step7_equil_nvt.rst" coordvelonly=true

#save MD core for more NPT equilibrating run

savesetup file="step8_equil_npt.mds"

echo -hs "|  running $nstep steps NPT equilibration"

echo -hs "|  it may take a while, please check back later"

runmd nstep=$nstep

echo -hs "|  output files:"

echo -hs "|    step8_equil_npt.mds"

echo -hs "|    step8_equil_npt.rst"

echo -hs "|    step8_equil_npt.trj"

echo -hs "|    step8_equil_npt.ene"

echo -hs "| "

echo -hs "|  please check energies profile step8_equil_npt.ene to make sure system is well equilibrated"

echo -hs "|  if not, please do more NPT equilibration by command:"

echo -hs "|    ./modar step8_equil_npt_more.inp -nstepmore=5000 -o step8_more_1.log"

echo -hs "|  if ok, please go final step to product sampling:"

echo -hs "|    ./modar -nobak step9_product.inp nstep=100000"

echo -hs "| "