AMBER Archive (2007)

Subject: Re: AMBER: PME and counter ions

From: Lars Skjærven (lars.skjarven_at_biomed.uib.no)
Date: Thu Dec 20 2007 - 09:12:53 CST

Hi !

Thank you all for valuable feedback and tips. A great and
comprehensive post here at the end from Thomas. Highly appreciated.
Thanks :-)

As for the instability of my system, I guess I was a bit inaccurate.
So, to clear this up: the complex contains several subunits. According
to the literature these subunits should stick together in a certain
way. However, we observe that some domains of the subunits drift apart
during the simulation. On the other hand, other domains of the
subunits sticks together as expected. It might be better to call this
unexpected behaviour, and not instability.. Anyway, I have loads of
great tips on how to figure this out.

I am also quite uncertain on my choice of pressure/volume/temperature
regulation(s). This might be a long shoot, but is there a general
recomandation what to use at the certain stages of the MD run
(assuming I am interested in the general dynamics of the system)? The
choice I made looks like this (as noted in earlier post):

1) Heating with NVT to 100K
2) Heating with NPT to 300K
3) Equilibration with NPT
4) Production with NVT and ntt=3

Would it be wiser to use ntt=1 instead? I believe I read in the
archive that ntt=1 would make more sense than ntt=3. However, the
amber manual writes: "Unless you are sure you know what you are doing,
please don't use ntt=1"..

Thanks again.

Lars Skjærven
University of Bergen
Norway

On Dec 20, 2007 7:05 AM, Thomas Cheatham <tec3_at_utah.edu> wrote:
>
> > 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
>
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu