# heat up MD system through Berendsen thermostat

# author: mengen

#

# usage: ./modar step6_heat.inp Ttar=300 Tinit=50 nstep=5000

#        Ttar is the target temperature

#        Tinit is the initial temperature

#        tTar default value is 300 k

#        nstep default is 1000

#

if(! -var ttar) {

  ttar=300

  echo -hs " no target temperature specified"

  echo -hs "   default value $ttar K will be used"

  echo -hsn "   or tell it now:"

  t=gettext()

  if($t>10 && $t<100000) ttar=$t

}

if(! -var tinit) {

  tinit=50

  echo -hs " no initial temperature specified"

  echo -hs "   default value $tinit K will be used"

  echo -hsn "   or tell it now:"

  t=gettext()

  if($t>10 && $t<1000) tinit=$t

}

if(! -var nstep) {

  nstep=5000

  echo -hs " no nstep for tempering specified"

  echo -hs "   $nstep steps will be"

  echo -hsn "   or tell it now:"

  n=gettext()

  if($n>100 && $n<10000) nstep=$n

}

ttcoup=0.1*($nstep/1000)

if(! -var rand) {

  rand=false

  echo -hs " generate initial velocity randomly or not"

  echo -hs "   default is not"

  echo -hsn "   or say yes here to change it:"

  asw=gettext()

  if($asw=="yes" || asw=="y") {

    rand=true

    echo -hs " initial velocities will be generated randomly"

  }

}

echo -hs "|  target temperatue: $ttar"

echo -hs "|  initial temperature: $tinit"

echo -hs "|  nstep: $nstep"

echo -hs "|  using Berendsen tempering"

if($rand=="true" && !-var iseed) {

  echo -hs "|  enable rand state"

  echo -hs "|  result cannot reproduce"

  iseed=int(randsys(100000))

}

else {

  if(!-var iseed) iseed=314159

}

echo -hs "|  random seed: $iseed"

#load force field

include "loadff.inc"

#load md system with minimized state

loadmsf file="step5_minimized.msf"

#load PBC crystal setup

include "pbc.inc"

tell msf_info

#setup image groups

imggroup method=resid sorted=true select="resname tip* pot cl* sod"

#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=6

#shake setup

:shake

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

#center group

putcenter select="segn PRO* RNA* DNA*"

#setup friction cooefficient for langevin may be used later

atom lngfric value=5.0 select="hydrogen"

atom lngfric value=10.0 select="! hydrogen"

#setup Berendsen thermostating

ensemble nvt tcoup=Berendsen ttcoup=$ttcoup tref=$ttar tinit=$tinit

#build MD core

build mdcore timestep=0.002 nprint=10 \

     trjfile="step6_heatup.trj" frqtrj=100 trjfmt=modar \

     rstfile="step6_heatup.rst" frqrst=100 \

     enefile="step6_heatup.ene" frqene=10 \

     randseed=$iseed randtype=old \

     frqresort=20000

runmd nstep=$nstep

savemsf file="step6_heatup.msf"

savemds file="step6_heatup.mds"

savepdb file="step6_heatup.pdb"

echo -hs "|  output files:"

echo -hs "|    step6_heatup.mds   the binary package of MD setup"

echo -hs "|    step6_heatup.rst   the MD state for restarting run"

echo -hs "|    step6_heatup.trj   the trajectory"

echo -hs "|    step6_heatup.ene   the energies profile"

echo -hs "|    please use viewene to check "step6_heatup.ene""

echo -hs "|    please use movar to view the trajectory"

echo -hs "| "

echo -hs "|  it's ready to go next step to do a short NVT equilibration"

echo -hs "|    ./modar -nobak step7_equil_nvt.inp nstep=2000"

echo -hs "| "