AMBER Archive (2007)

Subject: Re: AMBER: PME and counter ions

From: Thomas Cheatham (tec3_at_utah.edu)
Date: Thu Dec 20 2007 - 00:05:41 CST

> What sort of instability do you see? It could be that you should equilibrate
> longer with the restraints still in place. But it could also be that there is
> some feature of the starting coordinates that will cause "trouble" no matter

...just to add a bit more to this thread (and what will turn out to be a
truly long winded response with effort I am sure that my peers will
remind me should be targeted elsewhere)...

The issues of net-charge are special and counter-intuitive. As Chris Moth
mentioned, it is very logical to expect that carrying a net-charge in a
periodic system is unphysical and ridiculous or unrealistic. Yet, in
practice we can calculate free energies of ionic hydration (changing a
net-charge of an ion) that are surprisingly accurate even with only 64
waters surrounding the ion in a periodic cell!

In the late 90's, there were a string of papers from G. Hummer, R. Levy,
B.R. Brooks, and many others that showed and discussed this. The key, and
why this works, centers on how the Ewald sum is commonly calculated. A
crucial aspect of this, which is less intuitive, is the fact that although
the periodic energies should diverge (by a constant), the forces (which we
use for MD) are just fine. Short summary: With virtually all of the Ewald
implementations in common usage, you can run with a net-charge and be OK
(or equivalently, no one yet has definitively shown serious artifact in
simulation running with a net-charge compared to one without; a classic is
the Bader & Chandler 1992 paper looking at the PMF of two charged ions
approaching... Cutoffs led to serious artifact, yet Ewald, even with the
true periodicity, led to the results consistent with expectations for two
isolated ions approaching in water without periodicity; I am still trying
to understand this as it is completely counterintuitive).

Along with these ion-charging simulations, there have been a number of
studies regarding the influence of true periodicity and claims that this
leads to all sorts of serious artifacts. Papers by Hunenberger, McCammon
and others. In fact, if you look at continuum treatments comparing
"periodic" to non-periodic, significant artifacts are clearly present.
However, simulation in explicit solvent has not been able to reproduce
these "intuitive" and expected artifacts (seen with continuum
calculations). Examples include papers by Villareal et al. on peptides
trying to look at the Weber et al. artifacts, Smith & Pettitt on
dipole rotation, and others.

Why do we see no or few artifacts? Why do the energies not diverge?

Two issues; the boundary in the infinite limit (so-called tin-foil
boundary conditions) and the fact that we omit the k=0 term in the
reciprocal sum. Others can please correct me if I am wrong or missing
something here...

Most all Ewald implementations omit the zeroth order term in the
reciprocal sum. This leads to the idea that there is a net-neutralizing
plasma on the system that prevents the energy in the net-charge system
from diverging. In fact, we do not construct or add a special term that
is this ficticious net-neutralizing plasma. It is implicitly present due
to the omission of this term in the reciprocal sum. Complete descriptions
are provided by DeLeeuw, Perram and Smith and also Pollock. The second
issue relates to the Ewald sum in the limit. Effectively, most Ewald
implementions turn a conditionally convergent series (i.e. one whose
results depend on how you do the summation, for example as increasing
spheres vs. slabs) into a sum of two absolutely converging series (direct
space and reciprocal sum). We'll come back to this...

Now based on my limited understanding, if you run with a net-charge, no
matter how large, the periodic energies do not diverge, the forces are OK,
and you'd be hard pressed to find serious issues or artifact. This is
what Case was stating. A slight divergence is that AMBER also has the
ability to smear the net-charge-- this is talked about in the archives--
but this option is not really used in practice anymore. If you did use
this option, in fact, the effect on the system is also relatively modest,
for example if you have your net +250 charge, you can simply subtract
250/#atoms from the charge on each atom. If there are 25,000 atoms in the
system, this modifies each charge by ~0.01 which is likely far less than
the accuracy of the initial charge fit, but I digress... Both of these
options, i.e. smearing or the default (i.e. assuming there is this
nebulous net-neutralizing plasma) "work" and so far no one has shown
either are suspect (although clearly the smearing is less accurate as the
charge increases and the system size decreases).

The boundary conditions are a little more subtle and you can read about
"tin-foil", "conducting", or setting up boundary conditions at the
appropriate dielectric. If you do nothing, i.e. you throw out the
conditionally convergent part of the summation, you get tin-foil by
default. This corresponds to an infinite dielectric surrounding the
infinite periodic cell in the Ewald summation limit. Almost all Ewald
implementations (AMBER, CHARMM, NAMD) use these conditions since it is
easy, i.e. we do not need an extra term that takes care of this
conditional convergence. If you do want to consider this extra term, you
get conducting boundary conditions. This leads to a term that is
proportional to the dipole on the unit cell; see articles by Roberts and
Schnitker, ~1994-1995. This term may be important for small molecule
crystals (but read the Boresch papers). What I actually think may be the
most reasonable way--but what few people implement including us in the
AMBER code--is to think that in the limit, the dielectric should be
similar to what my system's dielectric is. This can be constructed
through reaction fields, etc. See work by Boresch and Steinhauser (1997).

So, now that I've told you more than you want to know about boundary
conditions and net-charge corrections, lets get to the heart of your
problem and try to decipher the very appropriate comments of Dave Case and
others. What you are seeing is likely not due to the 250 Na+ ions or
net-neutralizing plasmas (unless of course an explicit ion is in a very
unhappy place leading to local instability due to high forces that blow
the system up).

(1) what is stability?

When you say your simulation is unstable, this can mean many things. We
cannot interpret this. It could be SHAKE failures, it could be the
protein moves, it could be that the system blows up, it could be that
force constants are too high on some degrees of freedom, it could be
inaccuracy in the force field, it could be a poor initial model, it
could be something with imaging/wrapping, etc. Without more information
we cannot provide useful insight.

In general, we parse simulations into the "equilibration" and "production"
phases. Equilibration is how long we run the system such that all the
anomalous forces are gone and the initial structure is happy and unbiased
by its surroundings. The way many of us do this is to very weakly
restrain the structure of interest (aka protein) and then run as long as
possible simulations to relax the water, ions and surroundings. This
requires careful attention to issues of pressure, density, etc. If this
is done correctly, and the simulation still blows up, this is indication
of a more serious underlying problem (but we cannot tell you what it is
without further information). Serious problems usually means poor initial
structure, overlap leading to large initial forces, etc. etc.

(2) initial ion placement

LEaP, etc. use very crude electrostatics to place ions; these may or may
not be representative of reality. What I prefer to do (and this is
subjective) is (a) add solvent as suggested by Case, (b) add ions, and (c)
randomize the initial ion positions away from the solute of interest (aka
protein).

Then I run a long equilibration (at least 1-2 ns which is short by current
standards), allowing the ions and solvent to relax. While this is being
done, I weakly restrain the protein. Force constants of 0.5 kcal/mol-AA
are sufficient. Then I monitor density, pressure, etc. to see that the
system has relaxed.

Now back to the question of the "plasma" vs. explicit ions. If I had good
ion parameters, I would automatically want to include these as they are
more realistic (i.e. nature doesn't have this plasma) and I can track them
to look at specific ion interaction, lifetimes, etc. But, this is also
why I initially place them away from the solute and equilibrate as much as
I can (so I am not biased by the initial placement; this same argument
holds for crystal waters; normally I omit crystal waters and hope that the
simulation will spontaneously reproduce them). A way to randomize the
initial ion positions is with ptraj:

randomizeions :K+,Cl- around :R* by 6.0 overlap 4.0 noimage

The command above will move K+ and Cl- ions such that they are more than 6
angstroms away from :R* (an RNA) and move than 4 A from each other.

If I were to do a quick and dirty run (since I have not yet seen
a compelling argument that the plasma is "worse", I don't
have to worry about bad ion parameters that spontaneously
crystallize, etc; see work by Auffinger, Papoain, Pappu, and I don't
have to worry about the initial placement, ...), plasma (i.e. no added
ions) may be the way to go. However, if you look at my papers, I always
add ions since these are more real. Go figure?

--tom
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu