Re: [AMBER] Explicit Solvent Constant pH MD

From: Jason Swails <>
Date: Tue, 30 Jun 2015 23:11:43 -0400

On Tue, Jun 30, 2015 at 4:53 PM, Rahul Ramesh <> wrote:

> Hello
> I am using AMBER's constant pH MD to study the effect of drug on a polymer
> having carboxylic acid side chain. To start off with , I am currently
> trying to simulate a monomer which has a cellulose ring with a COOH group
> attached. The simulation is performed using amber 14 explicit solvent
> constant pH MD following Jason Swails tutorial. Since my system does not
> resemble any protein residue , I modified the * *code to
> incorporate my monomer. The force field used is AMBER 99 SB. The cpin file
> generated looks fine. The reference energy for the protonated deprotonated
> transition was obtained by performing a Thermodynamic Integration (TI) . As
> a test case , I set the solvph equal to the pKa of the carboxylate side
> group. While performing MD , there are no transitions between the
> protonated and deprotonated states of the monomer. The system is always in
> the protonated state for 2 ns.
> I have four questions regarding this ,
> 1) As of now , I have only one titratable residue ( monomer has 35 atoms
> totally ) which is initially in its protonated state , I would like to know
> if ions need to be added in the system while doing explicit solvent
> constant pH MD. If so how many ?

​I would add a "typical" amount. The ions don't play any direct role in
the protonation state change evaluations -- only through its influence of
the system dynamics. So the exact number in this case is not terribly

> 2) Is there any reason why transitions between protonated and deprotonated
> wouldn't happen? Is the amber 14 constant pH MD tuned only to titrate
> proteins?

​No, the code is completely general and will work for literally any
titratable residue definition that is properly defined in If
you are not getting transitions for the model compound when titrating at pH
= pKa, then your reference energy is wrong. It's really as simple as that
(which is nice). But there are so many ways that you can get the reference
energy wrong that it's impossible to guess exactly which one you happened
to do here.

> 3) I would like to know the units of the reference energy which are given
> ​​
> in the python code * *


4) Suppose I want to perform this simulation in implicit solvent just in
> case explicit solvent doesn't work, how can I obtain the partial charge
> information for my polymer accounting for implicit solvents

​The residue definition for implicit and explicit solvent are exactly the
same. And the reference energies are derived similarly and are often very
similar. Both use TI in GB as a starting point (which gets you very
close), and then uses a titration to compute the pKa of the model
compound. The shift between the actual computed pKa and the target pKa you
get from experiment can be used to tweak the reference energy to its
correct value.

The fact that TI and the procedure described above agree so well is
necessary, but not necessarily sufficient to show that the method is
implemented correctly. There are some disparities between the TI and
reference energies in some instances; for carboxylates there is an entropic
factor that will make the TI energy differ from the reference energy (there
are 2 sets of 2 degenerate proton positions on each oxygen for the
protonated state and only a single deprotonated state). But for histidine,
the TI and reference energies are spot-on.

One thing you need to take care with is that your TI has to be computed
with *exactly* the same settings (same solvent radii, same charge sets,
same residue with the same dummy atoms in the same places, same GB model,
... etc.). Any difference between these two will make your reference
energy wrong by a potentially very large amount. This is actually why the
implicit and explicit reference energies for AS4 and GL4 are so different
-- the solvent radii in the explicit definitions for the carboxylate
oxygens are smaller than in implicit solvent, so the reference energy is
therefore very different.

​In my opinion, what you are trying to do is a challenging, long-term
problem. I spent the better part of 1 to 2 years studying the method, the
code, and experimenting with different ideas before I understood the model
well enough to extend it the way I did to explicit solvent and to compute
reference energies and define titratable residues reliably and efficiently.
 (Actually, the ParmEd program was born from the collection of the tools I
wrote specifically to try out my ideas for constant pH MD).

If you are prepared to invest the necessary time to learn how to do this,
the best advice that I can give is step back from your current residue you
are trying to define. Instead, start with with something you *already*
know the answer for. Make a copy of (and the cpinutils
package). Delete the AS4 or GL4 definition in there, and try to create the
definition yourself. Once you have a recipe that you can use to reproduce
the same AS4 or GL4 definition as mine, then take that recipe and apply it
to your residue. Even better if you can take that recipe and turn it into
a generalized script that people can use for their *own* titratable

The concept of defining a titratable residue is itself simple and
straightforward. But this is an excellent example of the phrase "the
devil's in the details". There are a *lot* of details that are omitted
when people talk about defining references and reference energies (simply
because it would be way too verbose to include them all).

I hope this helped,

P.S., you are welcome to ask questions to the list as problems arise, but I
doubt there are very many people that will know the details well enough to
provide much detailed help. I will try and assist as time permits, but
this will likely require a significant time investment on your part.

Jason M. Swails
Rutgers University
Postdoctoral Researcher
AMBER mailing list
Received on Tue Jun 30 2015 - 20:30:02 PDT
Custom Search