Constant pH Molecular Dynamics Example - Calculating pKas for titratable side chains in HEWL

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Jason Swails & T. Dwight McGee Jr 2013

Section 3: Running the Production Simulations

Preparing the input files:

This section describes running the actual production simulations we will use to simulate our system. All of our results will be calculated from the data generated by the simulations here.

The simulations that should run in this phase depend on what you are trying to learn from your simulations. If you are interested only in the coupling between protein dynamics and protonation states in a particular pH environment, then you can run simulations at the single pH of interest. If you want to calculate pK_a values, on the other hand, you will need to calculate the pK_a based on a titration curve calculated over a pH range that spans the inflection points of your titratable residues.

The best pH values to run your simulations with depends on your system, what you want to learn, and the residues you are titrating. In this tutorial, we are interested in calculating pK_a values for the acidic-range residues of HEWL, as well as determining some coupling between protonation state and structure.

Therefore, we will run simulations at pH values of 0 through 7 with 1 pH-unit intervals.

The template input file for the production dynamics is the same as the equilibration stage from Section 2 (shown below)

Template mdin

Implicit solvent constant pH molecular dynamics &cntrl imin=0, irest=1, ntx=5, ntpr=1000, ntwx=1000, nstlim=1000000, dt=0.002, ntt=3, tempi=300, temp0=300, tautp=2.0, ig=-1, ntp=0, ntc=2, ntf=2, cut=30, ntb=0, igb=2, saltcon=0.1, nrespa=1, tol=0.000001, icnstph=1, solvph=<pH>, ntcnstph=5, gamma_ln=5.0, ntwr=10000, ioutfm=1, /

The input files we will use in our simulations have <pH> replaced with the target pH. A tarball with all 8 mdin files for pH values 0 through 7 can be downloaded here: mdins.tar.

Unpack the mdin files into the current directory with the command

tar xvf mdins.tar

Running the simulations

This simulation is identical to the equilibration dynamics we ran in Section 2, so we will use a cluster to run the calculations here as well. The template job file for all 8 simulations is shown below.

Template PBS script

#!/bin/sh #PBS -A UT-ABCD1234 #PBS -l walltime=10:00:00 #PBS -l nodes=2:ppn=12 #PBS -N CpHMD_Tutorial #PBS -o pbs.out #PBS -j oe #PBS -m abe #PBS -M fake.email@gmail.com # Change to the directory we submitted the job from cd $PBS_O_WORKDIR # Source the resource script necessary to run Amber on NICS Keeneland source ~/keeneland.setup # Run Amber on every requested node and core using Amber-installed mpich2 $AMBERHOME/bin/mpiexec -f $PBS_NODEFILE \ sander.MPI -O -i pH_<pH>.mdin -p 4LYT.parm7 -c 4LYT.equil.rst7 \ -cpin 4LYT.equil.cpin -o <pH>/4LYT.md1.mdout \ -cpout <pH>/4LYT.md1.cpout -r <pH>/4LYT.md1.rst7 \ -x <pH>/4LYT.md1.nc -cprestrt <pH>/4LYT.md1.cpin

To organize the output files by pH, I created a separate directory with the pH value as the directory name. Inside these directories, the file names are all the same. A tarball with the output files can be downloaded here: production_output.tar.bz2

Some important things to notice about this simulation:

The same equilibration starting structure is used for all of the simulations (the one generated in the previous section).
The cpin file was the cprestrt file generated by equilibration stage. This will set the "initial" protonation states for each titratable residue to the same state that the equilibration simulation ended in. Any subsequent simulation should be started with the cprestrt file from the previous run.

Click here to go to section 4

Click here to go back to section 2