#! /bin/bash

#rlc# # for /bin/sh use 1>fname 2>&1 in place of >&

#INIT {{{
init_structure () {
echo "Starting pdb2gmx"
pdb2gmx -ignh -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top << __EOD_pdb2gmx__
9
1
__EOD_pdb2gmx__
  echo "pdb2gmx finished"

# GENBOX ###
  editconf -f ${MOL}.gro -o ${MOL}.gro -d 1.0 1> output.genbox  2>&1
  genbox -cp ${MOL}.gro -cs -o ${MOL}_b4ions.gro -p ${MOL}.top 1> output.genbox  2>&1
  echo "editconf and genbox finished"

#em.mdp before adding ions{{{
cat <<__EOD_em__ > em.mdp
;define              =  -DFLEX_SPC
;constraints         =  none
;       Energy minimizing stuff
;
integrator          =  steep
emtol               =  1000.0
emstep              =  0.01
nsteps              =  50000
dt                  =  0.002    ; ps !
nstlist             =  1
ns_type             =  grid
coulombtype         =  pme
vdwtype             =  cutoff
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
pbc                 =  xyz
;
; Energy groups
;energygrps          = protein water_and_ions
__EOD_em__
#}}}

#run grompp to determine total charge and then run genion and re-run grompp and continue
  grompp -f em.mdp -c ${MOL}_b4ions.gro -p ${MOL}.top -o ${MOL}_b4ions.tpr 1> output.grompp_ions 2>&1
}
#}}}

#EM {{{
energy_minimization () {
make_ndx -f ${MOL}_ions << __EOD_make_ndx__
q
__EOD_make_ndx__

#em.mdp{{{
cat <<__EOD_em__ > em.mdp
;define              =  -DFLEX_SPC
;constraints         =  none
;       Energy minimizing stuff
;
integrator          =  steep
emtol               =  1000.0
emstep              =  0.01
nsteps              =  50000
dt                  =  0.002    ; ps !
nstlist             =  1
ns_type             =  grid
coulombtype         =  pme
vdwtype             =  cutoff
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
pbc                 =  xyz
;
; Energy groups
energygrps          = protein water_and_ions
__EOD_em__
#}}}
  grompp -n index.ndx -f em.mdp -c ${MOL}_ions.gro -p ${MOL}.top -o ${MOL}_em.tpr 1> output.grompp_em 2>&1
  echo "starting energy minimisation mdrun..."
  #mdrun -nice 19 -s ${MOL}_em -e ${MOL}_ener_em -o ${MOL}_em -c ${MOL}_em -v 1 -g ${MOL}_mdrun_em.log> output.mdrun_em  2>&1
  mdrun -nice 19 -deffnm ${MOL}_em -v 1 > output.mdrun_em  2>&1
  echo "em mdrun finished"
}
#}}}

#prNVT {{{
position_restrained_md_nvt () {
#prNVT.mdp{{{
cat <<__EOD_pr__ > prNVT.mdp
title               = NVT equilibration
continuation        = no
define              =  -DPOSRES
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002	; ps !
nsteps              =  50000	; total 100.0 ps.
nstcomm             =  10
comm_mode           =  Linear
comm_grps           =  System
nstfout             =  0
nstxout             =  100
nstvout             =  100
nstlog              =  100
nstenergy           =  100
ns_type             =  grid
nstlist             =  5
coulombtype         =  pme
vdwtype             =  cutoff
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.4
pbc                 =  xyz
fourierspacing      = 0.135
; Berendsen temperature coupling is on in two groups, solvent and ions together, all protein together
; use v-rescale as it is better than Berendsen
tcoupl              =  v-rescale
tc-grps		    =  Water_and_ions protein
tau_t               =  0.1 0.1
ref_t               =  298 298
; Pressure coupling is not on
pcoupl              =  no
dispcorr            =  enerpres
; Generate velocites is on at 298 K.
gen_vel             =  yes
gen_temp            =  298
gen_seed            =  -1
; Energy groups
energygrps          = water_and_ions protein
__EOD_pr__
#}}}
  grompp -n index.ndx -f prNVT.mdp -c ${MOL}_em.gro -p ${MOL}.top -o ${MOL}_prNVT.tpr 1> output.grompp_prNVT 2>&1
  echo "grompp for pr NVT mdrun finished"

  echo "starting pr NVT mdrun..."
  mpirun -np ${NUM_PROCS} mdrun_mpi -nice 19 -deffnm ${MOL}_prNVT > output.mdrun_prNVT 2>&1 < /dev/null

  echo "pr NVT mdrun finished"
}
#}}}

#prNPT {{{
position_restrained_md_npt () {
#prNPT.mdp{{{
cat <<__EOD_pr__ > prNPT.mdp
title               =  NPT equilibration
define              =  -DPOSRES
continuation        =  yes
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002	; ps !
nsteps              =  50000	; total 100.0 ps.
nstcomm             =  10
comm_mode           =  Linear
comm_grps           =  System
nstfout             =  0
nstxout             =  100
nstvout             =  100
nstlog              =  100
nstenergy           =  100
ns_type             =  grid
nstlist             =  5
coulombtype         =  pme
vdwtype             =  cutoff
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.4
;rvdw_switch         =  1.0
pbc                 =  xyz
fourierspacing      = 0.135
; Berendsen temperature coupling is on in two groups, solvent and ions together, all protein together
; use v-rescale as it is better than Berendsen
Tcoupl              =  v-rescale
tc-grps		    =  protein water_and_ions 
tau_t               =  0.1 0.1
ref_t               =  298 298
; Pressure coupling is not on
pcoupl              =  Parrinello-Rahman
pcoupltype          =  isotropic
tau_p               =  2.0
ref_p               =  1.0
refcoord_scaling    = all
compressibility     =  4.5e-5
dispcorr            =  enerpres
; Generate velocites is on at 298 K.
gen_vel             =  no
gen_temp            =  298
gen_seed            =  -1
; Energy groups
energygrps          = water_and_ions protein
__EOD_pr__
#}}}
  grompp -n index.ndx -f prNPT.mdp -c ${MOL}_prNVT.gro -p ${MOL}.top -o ${MOL}_prNPT.tpr 1> output.grompp_prNPT 2>&1
  echo "grompp for pr NPT mdrun finished"

  echo "starting pr NPT mdrun..."
  mpirun -np ${NUM_PROCS} mdrun_mpi -nice 19 -deffnm ${MOL}_prNPT > output.mdrun_prNPT 2>&1 < /dev/null

  echo "pr NPT mdrun finished"
}
#}}}

#MD{{{
full_md () {

######### MD #########
#md.mdp{{{
cat <<__EOD_md__ > md.mdp
title               =  Full MD run
continuation        =  yes
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002  ; ps !
nsteps              =  10000000  ; total 20000 ps.
nstcomm             =  1
comm_mode           =  Linear
comm_grps           =  System
nstxout             =  10000
nstvout             =  0
nstfout             =  0
nstlist             =  5
nstlog              =  10000
nstenergy           =  10000
ns_type             =  grid
coulombtype         =  pme
vdwtype             =  cutoff
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.4
;rvdw_switch         =  1.0
pbc                 =  xyz
fourierspacing      = 0.135
; Berendsen temperature coupling is on in two groups
; use v-rescale as it is better than Berendsen
Tcoupl              =  v-rescale
tc-grps		    =  water_and_ions protein
tau_t               =  0.1 0.1 
ref_t               =  298 298
; Pressure coupling is not on
pcoupl              =  Parrinello-Rahman
pcoupltype          =  isotropic
tau_p               =  2.0
ref_p               =  1.0
compressibility     =  4.5e-5
; Generate velocites is on at 298 K.
gen_vel             =  no
gen_temp            =  298
gen_seed            =  -1
; Energy groups
energygrps          = water_and_ions protein
__EOD_md__
#}}}
  grompp -n index.ndx -f md.mdp -c ${MOL}_prNPT.gro -p ${MOL}.top -o ${MOL}_md.tpr 1> output.grompp_md 2>&1

  echo "grompp for mdrun finished"

  echo "starting mdrun for full md..."
  mpirun -np ${NUM_PROCS} mdrun_mpi -nice 19 -deffnm ${MOL}_md -v 1 > output.mdrun_md 2>&1 < /dev/null

  echo "full md mdrun finished"
}
#}}}
################################## end of function definitions ########################

# get molecule name from command-line argument
MOL=$1
MOL=`echo $MOL | sed 's/.pdb//'`
export MOL

# set the number of processors
NUM_PROCS=2
export NUM_PROCS

OS=`uname -s`
OSm=`uname -m`
if [[ $OSm == "x86_64" ]]; then
  OS=${OS}_${OSm}
fi
#. /software/gromacs/gromacs-4.5.3/${OS}/bin/GMXRC

init_structure

################# FIX CHARGES: ##########################################################
# Look at output.grompp_ions to to check for a net charge, then
##   add negative chloride ions to neutralize system
#genion -s ${MOL}_b4ions.tpr -o ${MOL}_ions.gro -p ${MOL}.top -g ${MOL}_genion.log -nn 4 -nname CL -seed 999 <<__EOD_GENION__
#13
#EOD_GENION__
##   add positive sodium ions to neutralize system
#genion -s ${MOL}_b4ions.tpr -o ${MOL}_ions.gro -p ${MOL}.top -g ${MOL}_genion.log -np 1 -pname NA -seed 999 << __EOD_GENION__
#13
#__EOD_GENION__
#############################################################################################

#energy_minimization
#position_restrained_md_nvt
#position_restrained_md_npt
#full_md