# do minimization

# author: mengen

#

# usage: ./modar step5_do_minimization.inp [nstep1=500] [nstep2=500] [nstep3=500] [nstep4=500]

#

#load force field

include "loadff.inc"

#load md system

loadmsf file="step4_all.msf"

#load PBC crystal setup

include "pbc.inc"

tell msf_info

#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

echo -hs " will do 4 stages minimization"

echo -hs "   stage 1. only water molecules, others freezed"

echo -hs "   stage 2. only water and ions, others freezed"

echo -hs "   stage 3. every thing except the main chain atoms"

echo -hs "   stage 4. whole system"

if(! -var nstep1) {

  nstep1=200

  echo -hs " number cycles of stage1 not specified"

  echo -hs "   default value $nstep1 will be used"

  echo -hsn "   or tell it now:"

  nstep=gettext()

  if($nstep>0) nstep1=$nstep

}

if(! -var nstep2) {

  nstep2=200

  echo -hs " number cycles of stage2 not specified"

  echo -hs "   default value $nstep2 will be used"

  echo -hsn "   or tell it now:"

  nstep=gettext()

  if($nstep>0) nstep2=$nstep

}

if(! -var nstep3) {

  nstep3=200

  echo -hs " number cycles of stage3 not specified"

  echo -hs "   default value $nstep3 will be used"

  echo -hsn "   or tell it now:"

  nstep=gettext()

  if($nstep>0) nstep3=$nstep

}

if(! -var nstep4) {

  nstep4=200

  echo -hs " number cycles of stage4 not specified"

  echo -hs "   default value $nstep4 will be used"

  echo -hsn "   or tell it now:"

  nstep=gettext()

  if($nstep>0) nstep4=$nstep

}

#do minimization for solvent

echo -hs "|  running $nstep1 steps minimization for solvent, please wait a moment"

atom fix select="! resname tip*"

domin method=MC nstep=$nstep1 nprint=10 ftol=1e-5 iseed=314159 \

      trjfile="step5_mini1.trj" frqtrj=10 \

      enefile="step5_mini1.ene" frqene=1

tell pot

echo -hs "| total potential energy is: $pot kcal/mol"

#do minimization for solvent with ions

echo -hs "|  running $nstep2 steps minimization for solvent and ions"

atom fix select="! resname tip* pot cl* sod"

domin method=SDfuse nstep=$nstep2 nprint=10 ftol=1e-5 fuse=0.2 \

      trjfile="step5_mini2.trj" frqtrj=10 \

      enefile="step5_mini2.ene" frqene=1

tell pot

echo -hs "| total potential energy is: $pot kcal/mol"

#do minimization for side chain, solvent and ions

defcls name="mainchain" select="(name CA N O C && segn PRO*) \

                || (segn DNA* RNA* && name C3' C4' C5' P)"

if($nasel>0) {

  echo -hs "|  running $nstep3 steps minimization for solvent, ions and side chains"

  atom fix select="mainchain"

  domin method=steep nstep=$nstep3 nprint=10 ftol=1e-5 fuse=0.2 \

       trjfile="step5_mini3.trj" frqtrj=10 \

       enefile="step5_mini3.ene" frqene=1

  tell pot

  echo -hs "| total potential energy is: $pot kcal/mol"

}

echo -hs "|  running $nstep4 steps minimization for whole system"

#do minimization for all

atom fix select="none"

domin method=SD nstep=$nstep4 nprint=10 ftol=1e-5 fuse=0.2 \

      trjfile="step5_mini4.trj" frqtrj=10 \

      enefile="step5_mini4.ene" frqene=1

tell pot

echo -hs "| total potential energy is: $pot kcal/mol"

#save msf and pdb

savemsf file="step5_minimized.msf"

savepdb fmt=pdb file="step5_minimized.pdb"

tell logfilename

echo -hs "|  output files:"

echo -hs "|    step5_minimized.mtf"

echo -hs "|    step5_minimized.pdb"

if(-var slogfn) {

  exec "grep minstep $slogfn | grep -v grep | awk '{print \$4}' > minipot.txt"

  echo -hs "|  please plot potential profile "minipot.txt" to check minimization result"

  echo -hs "|  for example: ./xmgr minipot.txt"

}

echo -hs "|  it's ready to go next step to heatup system by command:"

echo -hs "|    ./modar -nobak step6_heatup.inp"

echo -hs "| "