Constant pH simulation of BPTI Explicit(Amber)

Creating Input Files:

loadAmberParams frcmod.ionsjc_tip3p(if this is not loaded, there will be error "Could not find vdW (or other) parameters for type: Cl-" )
source leaprc.constph
# Load the PDBs
6pti = loadPDB 6pti_fixed.pdb
# Make the disulfide bonds
bond 6pti.5.SG 6pti.55.SG
bond 6pti.14.SG 6pti.38.SG
bond 6pti.30.SG 6pti.51.SG
# Solvate
solvateBox 6pti TIP3PBOX 10.0
#Add ions(make the salt concentration 0.1M)
addIonsRand 6pti CL 13 NA 7
# Save topology files
saveAmberParm 6pti 6pti.parm7 6pti.rst7
# Quit

tleap -f

Creating Cpin File -p 6pti.parm7 -o 6pti.cpin -op 6pti.mod0.parm7

note: edit line 256 in


-op: Carboxylate residues in explicit solvent simulations require a modified topology file, thus need to print a new parameter file.

System Minimization

Minimization to prepare for constant pH MD
imin = 1,                                   ! Turn on minimization
ncyc=1000,                              ! If ncyc<maxcyc steepest descent before switching to conjugated gradient for remaining
maxcyc = 5000,                       ! Total number of minimization cycles
cut = 8,                                    ! Nonbonded cutoff in angstroms
ntpr=50,                                   ! Print frequency
restraint_wt = 10,                     ! 10 kcal/mol/A**2 restraint force constant
restraintmask = '@CA,C,O,N'     ! Restraints on the backbone atoms only

mpirun -np 12 pmemd.MPI -O -i min.mdin -p 6pti.mod0.parm7 -c 6pti.rst7 -o 6pti.min.mdout -r 6pti.min.rst7 -ref 6pti.rst7 -cpin 6pti.cpin

Heating The System

Explicit solvent constant pH initial heating mdin
imin=0, irest=0, ntx=1,
ntpr=1000, ntwx=1000, nstlim=200000,
dt=0.002, ntt=3, gamma_ln=5.0, ig=-1,
TYPE='TEMP0', ISTEP1=0, ISTEP2=150000,
VALUE1=10.0, VALUE2=300.0,
&wt TYPE='END' /

mpirun -np 12 pmemd.MPI -O -i heat.mdin -c 6pti.min.rst7 -p 6pti.mod0.parm7 -ref 6pti.min.rst7 -cpin 6pti.cpin -o 6pti.heat.mdout -r 6pti.heat.rst7 -x

Equilibrating The System

Explicit solvent constant pH molecular dynamics
imin=0, irest=1, ntx=5,
ntpr=1000, ntwx=1000, nstlim=1000000,
dt=0.002, ntt=3, tempi=300,
temp0=300, gamma_ln=5.0, ig=-1
ntp=1, ntc=2, ntf=2, cut=8,
ntb=2, iwrap=1,ioutfm=1,

mpirun -np 12 pmemd.MPI -O -i equi.mdin -p 6pti.mod0.parm7 -c 6pti.heat.rst7 -cpin 6pti.cpin -o 6pti.equil.mdout -cpout 6pti.equil.cpout -r 6pti.equil.rst7 -x -cprestrt 6pti.equil.cpin

heating is NVT, and equilibrating is NPT, "Another equilibration issue involves the use of the NPT ensemble. If the system is not reasonably close to a good minimum when you initiate a constant pressure simulation, you can get P-V oscillations. One way to get a decent starting structure for NPT is to start with a short NVT simulation. This is what we do in the Amber water demo".   --George Seibel

It is better use CPU code than GPU code for equilibration,since GPU code does not automatically reorganize grid cells and if the density changed too much, GPU code will fail.

Production Run

Explicit solvent constant pH molecular dynamics
imin=0, irest=1, ntx=5,ntxo=2,
ntpr=1000, ntwx=1000, nstlim=1000000,
dt=0.002, ntt=3, tempi=300,
temp0=300, gamma_ln=5.0, ig=-1,
ntc=2, ntf=2, cut=8, iwrap=1,
ioutfm=1, icnstph=2, ntcnstph=100,
solvph=3.0, ntrelax=100, saltcon=0.1,

*icnstph=1 implicit
=2 explicit

for i in $(seq 0 14) ; do sed "s/solvph=3.0/solvph=$i/" md.mdin > pH_$i.mdin; done

for ph in $(seq 0 14); do mpirun -np 2 pmemd.MPI -O -i pH_${ph}.mdin -p 6pti.mod0.parm7 -c 6pti.equil.rst7 -cpin 6pti.cpin -o ${ph}/6pti.md0.mdout -cpout ${ph}/6pti.md0.cpout -r ${ph}/6pti.md0.rst7 -x ${ph}/ -cprestrt ${ph}/6pti.md0.cpin & done

Calculate pKas

for ph in $(seq 0 14); do cphstats -i 6pti.cpin ./${ph}/6pti.md0.cpout -o pH${ph}_calcpka.dat --population pH${ph}_populations.dat; done

This will generate 2 files for each pH: pH<ph>_calcpka.dat & pH<ph>_populations.dat 0_14

Multiple Run

pK1/2 of titrable residues in BPTI
#res name trial1 trail2 trial3 trial4 trial5 ave\pmstd
3 ASP 3.17 3.13 3.34 3.52 3.22 3.27\pm0.14
7 GLU 4.73 4.77 4.89 4.68 4.82 4.78\pm0.07
10 TYR 11.15 10.59 10.35 10.07 10.63 10.56\pm0.36
15 LYS 10.23 10.14 10.15 9.94 9.93 10.02\pm0.04
21 TYR 11.60 11.47 11.51 11.54 11.60 11.54\pm0.05
23 TYR 12.08 11.68 11.73 12.08 11.65 11.51\pm0.17
26 LYS 9.74 10.26 10.28 10.25 10.26 10.16\pm0.21
35 TYR 11.66 12.08 12.47 12.49 12.17 11.27\pm0.43
41 LYS 9.79 9.74 9.74 10.09 10.12 9.90\pm0.17
46 LYS 9.52 9.54 9.84 9.66 9.54 9.62\pm0.12
49 GLU 3.67 3.92 3.50 3.76 3.95 3.76\pm0.17
50 ASP 2.74 2.71 3.62 3.33 3.33 3.15\pm0.36