[gmx-users] pH

[baptista at itqb.unl.pt]

Dear Osmany,

>  Hi
> Two questions.I have my system(for example a peptide in a box of water)
> How can I calculate the pH of the sysytem and

First of all, I sincerely doubt that you really want to *compute* the pH.
The reason can be made more clear if I make a parallel with a more
familiar property, the temperature. (I'll ignore volume/pressure stuff in
what follows.) The ideal way to perform a simulation of your solvated
peptide is to use a constant-temperature algorithm, which directly
reflects the heat-bath effect of the surrounding solution. Imagine,
though, that there were no constant-temperature algorithms around and all
you had was a plain constant-energy algorithm. Of course, you could
perform a constant-energy simulation and *compute* the energy afterward
(from the average kinetic energy), but the only reason to do so would be
to try finding out an energy value giving the temperature you wanted for
your system: with it you could then make the simulation at the intended
temperature (note that, in practice, the energy value alone is useless:
eg, if you want to relate your simulations with experiments done in your
peptide at 300K, you have no idea what was the typical energy of a
"peptide+water box" in that experiment). But the constant-temperature
algorithm is the ideal way of performing the simulation, since you don't
need to waste time trying out different energy values. Furthermore, for
the small (non-macroscopic) systems used in simulations you have an
additional, extremely important reason for preferring a
constant-temperature algorithm: a constant-energy algorithm neglects
energy fluctuations that may be crucial for the system. The fluctuations
of all "mechanical" properties (phase functions) in a constant-temperature
system are larger than in a constant-energy one, which, eg, helps the
system to overcome local energy barriers, not to mention the entropic role
of fluctuations. Theoretically, this may be stated by saying that you
should use a proper ensemble when getting fluctuations. Thus, the natural
thermodynamic parameter for your system is the temperature, not the
energy, both from experimental and theoretical viewpoints. Thus, you want
to *impose* the temperature in your simulation, not *compute* it.

The situation concerning pH is identical. A pH-buffered solution works
both as a heat bath and as a H+ bath. Thus, the ideal situation would be
to have a constant-pH algorithm to perform the MD simulations. However,
conventional MD algorithms cannot do that, since they're constant-H+
algorithms. Furthermore, if you use a MM force field, you will never get
protonation/deprotonation of sites in your solute: free H+ will just
wander around, while titrable sites are stuck with their initial
protonation states; if you want to overcome this you have to use a non-MM
Hamiltonian which explicitly allows for proton transfer. In any case,
regardless of the Hamiltonian used, the amount of H+ would never change.
Of course, for a given amount of H+, you may wander to which pH value this
corresponds and even try to *compute* it (eg, using some variant of
Widom's particle insertion method), but, as in the energy/temperature
case, that is not the interesting thing to do. The interesting thing would
be to devise a proper algorithm for constant-pH MD and then *impose* the
intended pH value in your simulation. Which lead to your second
question...

> How can I do a simulation of this system whith a given pH.
> For example I want a simulation of the system whith pH=4 and another
> whith
> pH=9..

The simplest thing to do is to stick to conventional MD and just select
the particular amount of H+ that is "typical" for your system. As noted in
other replies, you can neglect free H+: at most pH values (even at pH=4)
its concentration is orders of magnitude below that of other counterions
(eg, Na+, Cl-). Furthermore, as noted above, a pure MM Hamiltonian will be
unable to move protons between titrable sites (eg, unlike EVB), so you
actually need to select a particular configuration of protonation states,
ie, you need to choose the protonation state of each titrable site in your
molecule. As also noted in some of the previous replies, you can get a
good guess by performing a standard pKa calculation with your starting
structure; these are standard procedures, mostly based on continuum
electrostatics for computing protonation free energies (eg, using MEAD,
UHBD, DELPHI, APBS) and Monte Carlo to perform sampling of protonation
states (eg, using REDTI, PETIT). Unfortunately, as the conformation of the
solute changes along the MD simulation, the protonation states may become
inadequate, due to the strong protonation-conformation coupling that
exists in many cases. Several solutions have been proposed (including by
myself), but most of them are more or less heuristic attempts. The truly
satisfactory solution would be to have a constant-pH MD method, and some
have been proposed in the last years. I have discussed that issue recently
in this list, so please have a look at
http://www.gromacs.org/pipermail/gmx-users/2003-April/031920.html. (Note
that there is a serious theoretical problem with a method by Phil
Hunenberger, whose paper was linked to by a previous reply; see my
previous message.) Unfortunately, constant-pH methods are recent and still
under development and/or testing. Hopefully, they will become standard
methods in the near future. Maybe you could one of these days just specify
"pH = 7.0" in the GROMACS .mdp file and see protonation states changing
during the MD run! :-) Until then, the best solution is probably to one
mentioned above: make an initial good guess with a standard pKa
calculation, and eventually check it later with snapshots from the MD run.

> Thanks...
>

Regards,
Antonio

-- 
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica, EAN, ITQB II, Piso 6, Apartado 127
2781-901 Oeiras, Portugal
phone: +351-214469619(NEW!)   email: baptista at itqb.unl.pt
fax:   +351-214411277         WWW:   http://www.itqb.unl.pt/~baptista
--------------------------------------------------------------------------