if(! -var nstep) getvar nstep default=500000 min=0 \

                        msg="how many steps to run"

#load force field

include "loadff.inc"

#load MSF built previously

loadmsf file="step3_melted.msf"

#load MD core

loadmdcore file="step3_melting.mds"

loadstate  file="step3_melting.rst"

#aneal to 300K through Berendsen thermostating

ensemble NVT tcoup=Berendsen ttcoup=0.2 tref=300.0

#change output files name

mdcp trjfile="step4_annealing.trj" frqtrj=100 trjfmt=movar \

     enefile="step4_annealing.ene" frqene=10 \

     rstfile="step4_annealing.rst"

#run MD

runmd nstepmore=2500

#save PDB

savecrd fmt=pdb file="step4_annealed.pdb"

#change thermostat to Nose-Hoover

ensemble tcoup=nosehoover tmass=1000

#change output files name

mdcp trjfile="step4_equil_nvt.trj" frqtrj=500 trjfmt=movar \

     enefile="step4_equil_nvt.ene" frqene=100 \

     rstfile="step4_equil_nvt.rst" \

     nstepmore=$nstep

#save MD core

savemds file="step4_equil_nvt.mds"

if($nstep==0) {

  echo -hs "please run MD by command"

  echo -hs " ./modar step4_equil.mds nstepmore=1000000"

  stop

}

#run MD

runmd