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

                              msg="How may steps to run"

#load force field

include "loadff.inc"

#load MSF built previously

loadmsf file="step2_solvbox.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

#heat up using Berendsen thermostating

ensemble NVT tcoup=Berendsen ttcoup=0.2 Tinit=100 Tref=500

#build MD core

build mdcore itype=vverlet timestep=0.002 nprint=100 \

     trjfile="step3_heatup.dcd" frqtrj=100 trjfmt=charmm \

     enefile="step3_heatup.ene" frqene=10 \

     rstfile="step3_heatup.rst" frqrst=100000 \

     flushene=true \

     randseed=314159 randtype=old \

     frqresort=20000

#run MD

runmd nstepmore=2500

#save PDB

savecrd fmt=pdb file="step3_heatup.pdb"

#change thermostat to Nose-Hoover

#atom lngfric value=5 select="hydrogen"

#atom lngfric value=10 select="! hydrogen"

ensemble tcoup=nosehoover tmass=1000

#change output files name

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

     enefile="step3_melting.ene" frqene=100 \

     rstfile="step3_melting.rst" \

     nstepmore=$nstep

#save MD core

savemds file="step3_melting.mds"

#run MD

runmd

#save mds msf and pdb

savemds file="step3_melted.mds"

savemsf file="step3_melted.msf"

savecrd fmt=pdb file="step3_melted.pdb"