# ==============================
# Source the necessary Tcl files
# ==============================
# This depends on where and how NAMD is installed on your computer.
# Relative to the main NAMD directory, the appropriate file is
# /lib/namdcph/namdcph.tcl
#
source ../../lib/namdcph/namdcph.tcl

# ==================================
# These keywords all follow as usual
#
set topo_dir "../topology"
structure      $topo_dir/bbl.hmr.psf
coordinates    $topo_dir/bbl.pdb
temperature    298.0
extendedSystem $topo_dir/bbl.xsc

wrapWater      on
wrapAll        on
outputEnergies 5000
DCDFreq        5000

timestep           2.0
fullElectFrequency 1
rigidBonds         ALL
langevin           on
langevinTemp       298.0
langevinDamping    1.0
langevinHydrogen   no

switching         on
VDWForceSwitching on
LJCorrection      on
switchDist        10.0
cutoff            12.0
pairlistDist      14.0
exclude           scaled1-4
1-4scaling        1.0
PME               on
PMEGridSpacing    1.0

# ===============================================
# Begin constant-pH MD keywords and modifications
#
# Load force field files as usual, but add constant-pH specific parameters
set toppar_dir "../../toppar"
paratypecharmm on
parameters $toppar_dir/par_all36_prot.prm
parameters $toppar_dir/par_cph36_prot.prm
parameters $toppar_dir/par_all36_solvent.prm
# Load constant-pH specific topology files 
cphConfigFile $toppar_dir/conf_cph36_prot.json
topology $toppar_dir/top_cph36_prot.rtf

# We will be running multiple pH values sorted into their own directories, but
# otherwise using the same naming scheme.
set pHList [list 2.8 3.2 3.6 4.0 4.4 4.8 5.2 5.6 6.0]
set pH [lindex $pHList [myReplica]]
pH $pH
outputname        $pH/bbl_prod0
stdout            $pH/bbl_prod0.log
cphMDBasename     $pH/namdcph.md
cphSwitchBasename $pH/namdcph.sw

# Experimental values from Arbely, et al.
# For clarity, each residue is sorted by type with a separate command:
#
# ASP
cphSetResiduepKai PROA:129:ASP 3.9 PROA:145:ASP 3.7 PROA:162:ASP 3.2
# GLU
cphSetResiduepKai PROA:141:GLU 4.5 PROA:161:GLU 3.7 PROA:164:GLU 4.5
# HIS - Note that only experimental macroscopic pKas are reported. We assume
# that the relative microscopic pKa is the same as for the reference values and
# then adjust to match the experimental values (H142 - 6.5, H166 - 5.4).
cphSetResiduepKai PROA:142:HIS {6.6 7.1} PROA:166:HIS {5.5 6.0}

# With the current settings this implies 1 ps between switching attempts,
# which will be 15 ps in length.
#
cphNumstepsPerSwitch 7500
# Minimize after randomizing protonation states at start.
cphNumMinSteps 100
cphRun 500 5
# Don't exit until all simulations have finished.
replicaBarrier
exit