Questions and problems?

Equilibration

There are two aspects of equilibration - getting things started at all (mainly addressed under the topic blowup), and the finer point of how to arrive at a truly equilibrated system. Here, we deal mainly with the latter topic, presenting a range of opinion.

Note that when a system is solvated by superimposing water boxes and subtracting overlapping waters, constant pressure equilibration must be used to achieve proper density because of vdw voids at the solute and box boundaries. Thus to achieve a 10A box, an initial ~12.5A box is required. Denser initial packing can be achieved in Leap by setting the closeness parameter less than the default of 1.0 - in this case there will be irregular close contacts and lesser voids.

I'm trying to achieve a temperature of 300K, and pressure 1.0 bar. I start with a much-minimized initial structure. So far, I've gone to 30 ps and the temperature has climbed to 500K and pressure shows no signs of stabilizing: it fluctuates wildly between + & - 300 bar! I think I must have some parameters set incorrectly ...

Did you use gradual warming? Even when things are well-minimized, cold-starting at 300K can create some high energies on the 1st step that take a long time to resolve and can e.g. mess up the structure of nucleic acids. The pressure is more or less consistent with the performance of the constant pressure algorithm used in AMBER.

Andrew Pohorille:

There is no universal strategy and some physical insight is needed. Two general points may be worth mentioning though:

  1. probably nothing disorders the system better than thermal energy; however one should remember that this is a double-edged sword - too high temperature may destroy the integrity of the system (or its parts).
  2. there is great virtue to an old (and simple) technique of occasional resampling velocities from the Maxwell-Boltzmann distribution during equilibration. In particular, it will take care of hot and/or cold spots.
Bill Ross:

Nucleic acids seem to need extra care in equilibration. In water, I first minimize the whole system for 50-500 steps to resolve bad steric contacts, then do a constant volume run at high temperature holding the DNA or RNA (but not counterions) fixed in order to disorder the periodicity of the 216-water subboxes. Then I start again without the velocities and warm gradually from 10-300K over 2-20 ps with constant pressure.

An alternative would be to restrain the solute atoms to their initial positions using a harmonic potential, and gradually relax the potential.

George Seibel:

As a general rule, I don't recommend [gradual warming as an equilibration technique]. The model that you start with is most likely not at its global minimum; parts of the system will usually be trapped in local minima that may still be very high in potential energy. In the course of the simulation, these regions may be dislodged, and their PE will be transformed into a high level of kinetic energy, localized in a small area. The system will now have a "hot spot". The Berendsen temperature control algorithm will be of no help in dealing with a hot spot; it will scale the velocities of all particles such that a desired average temperature is maintained. The net result will be to cool the entire system, without significantly altering the ratio between the hottest and coldest parts. As the temperature is raised, all velocities will be scaled up, but the ratio will again not change. If the simulation temperature is 10 K and you get a couple residues at 1000 K, the ratio of hot to cold parts is 100. If instead you start the simulation at 300 K and get a couple of 1000 K residues, the ratio of hot to cold is only 3.3. Based on this reasoning, one might choose to start the simulation at a HIGHER temperature than the target temperature, and scale DOWN instead. One might also periodically stop the simulation and randomly reassign a Maxwellian distribution of velocities. I usually just start at the target temperature and leave it there.

If the system is sufficiently well coupled, energy will transfer throughout the system and the hot and cold spots will work themselves out. The problem is, many systems are not that well coupled, and it's hard to say how long (if ever) it will take for thermal equilibration to occur. Longer nonbonded cutoffs will tend to improve coupling.

Another equilibration issue involves the use of the NPT ensemble. If the system is not reasonably close to a good minimum when you initiate a constant pressure simulation, you can get P-V oscillations. One way to get a decent starting structure for NPT is to start with a short NVT simulation. This is what we do in the Amber water demo. 


George Seibel, SmithKline Beecham
seibelgl@sb.com
Thomas Cheatham, III:

In some sense, the "structure" is not too sensitive to the initial equilibration protocol as long as the water structure is sufficiently relaxed. This happens rather quickly (within ~25-100 ps). See 

	Norberto de Souza and Ornstein, 
	J Biomol Struct Dyn 14, 607-611 (1997)
for a discussion of "equilibration". Note also that whether the entire system is relaxed within the time frame of "equilibration" is debatable; in many simulations, the DNA has not even fully equilibrated on a multi-nanosecond time frame in terms of some properties (such as the agreement between symmetric halves of a duplex, etc) however A->B transitions are seen on a < 1 ns time frame. It all depends what you are looking at and what properties you are trying to represent.

Also see Tom Cheatham's discussion of equilibration in his DNA tutorial.

Using an electrostatic cutoff of 13A I observed a nearly identical dynamic behaviour for the two systems in the 30ps trajectories.

Because the electrostatic interactions are very long ranging, such a short cutoff may not be sufficient to get the mean electrostatics correctly.

Thus, the number of counter ions seems to have a considerable influence on the dynamic behaviour of the protein. ... What would be an appropriate parameter to decide what is the most realistic number of counter ions?

The Maxwell relaxation time is the relevant number: 


	tM = lDH^2 / D
where lDH is the Debye length and D the diffusion constant for the ion. For mixtures of ions see: 

  Joseph B. Hubbard (1987).
  Non-equilibrium theories of electrolyte solutions.
  The Physics and Chemistry of Aqueous Ionic Solutions.
  Dortrecht, D. Reidel. 95--128.
tM is the relaxation time of an electrolyte if you apply a small perturbation. To be sure to get the influence of the electrolyte (co- and counter ions) not completely wrong you have to simulate at least a time range of tM. 

Eberhard von Kitzing, Max-Planck-Institut fuer Medizinische Forschung,
vkitzing@silia.mpimf-heidelberg.mpg.de