(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 Ross Walker 2006

pKa Calculations using Thermodynamic Integration - SECTION 1

pKa Calculations Using Thermodynamic Integration

By Ross Walker & Mike Crowley

1) Setup the prmtop and inpcrd files for the AspH and Asp- models.

We will the calculations to be run using Amber's IGB=5 model and the FF94 force field. Newer force fields are available and maybe more appropriate but for purposes of this tutorial we will use the same methodology as they did in the paper. To run thermodynamic integration in sander, you need to create two individual prmtops representing the two ends of the thermodynamic cycle. The requirement is that these two prmtop files have exactly the same number of atoms and atom ordering. This means that we cannot simply create ACE-ASH-NME and ACE-ASP-NME but instead we need to follow the approach they used in the paper which is to create two ACE-ASH-NME systems but in the second case we zero the charge on the proton and redistribute the charges on the other atoms in line with what they are defined as for the ASH residue.

One problem we would encounter with this approach is how to make the VDW terms involving the proton disappear. Fortunately in the FF94 force field all protons attached to oxygen atoms were parameterised with a van der Waals radius of zero. This means that there is no VDW interaction with the hydroxyl proton and thus zeroing it's charge removes all of the non-bond interactions with this proton. We don't need to worry about removing the bonds, angles and dihedral terms as these are all internal to the residue and so the average properties of these terms should be identical in the model  and protein systems. There is also no contribution from these terms to the derivative of the potential with respect to lambda.

We will thus setup two ACE-ASH-NME systems with the following charges, as shown in figure 2 of the paper. Note for the work in the paper they made the assumption that only the charges local to the hydroxyl proton change on deprotonation so for the ASP residue they set the two oxygen atom charges to match the FF94 force field values but then folded the changes in all the other atoms into the C(gamma) atom. Thus the charges to be used for the two systems are:

The top charges are the same as those defined in the force field for ASH. You can check this yourself in Xleap if you want.

I have created a simple leap script to do the setup for us. Here it is:

> $AMBERHOME/exe/tleap -s -f create_ash.leaprc

You should check the output to make sure it is working correctly. I.e. does the ASPH model have a charge of -1.0? etc.

This will create: asph_model.prmtop, asph_model.incprd, asp_model.prmtop and asp_model.inpcrd

The next step is to run the thermodynamic integration calculations on this model compound to find delta G(model).

---

CLICK HERE TO GO TO SECTION 2

---

(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 Ross Walker 2006