[gmx-developers] Reaction Filed crash!

David van der Spoel spoel at xray.bmc.uu.se
Sun Dec 18 16:29:30 CET 2011

On 12/18/11 1:27 AM, baptista at itqb.unl.pt wrote:
> On Sat, 17 Dec 2011, Berk Hess wrote:
>> I will try to clarify the situation a bit.
>> Pre-4.5 there were 4 approximation/errors in this typical "Gromos" setup
>> and a complicating factor:
>> 1. Non-buffered pair-lists reused for nstlist-1 steps: no real
>> Hamiltonian
>> 2. Reaction-field with epsilon_rf<infinity: discontinuous energy and
>> force at cut-off
>> 3. Twin-range interactions (or multiple time stepping) with 5*2=10 fs
>> for the forces
>> beyond 0.9 nm, but the reorientation of water only allows time steps
>> up to 5 fs
>> 4. A non-reversible multiple time stepping scheme
>> 5. A thermostat with strong coupling (pre-4.0 also Berendsen: one
>> other error)
>> Factors 1-4 cause serious integration errors in the long-range
>> electrostatics
>> which all lead to serious heating of the system, but different for
>> water and protein.
>> Factor 5 couples off this heat, but with such strong heating you no
>> longer know
>> what happening.
>> In version 4.5 I fixed factor 4. This takes out one source of error
>> and heating.
>> But as different errors affect different aspects/components
>> differently, it's very
>> difficult to tell what will happen.
>> Factors 1-3 can be fixed at a high computational cost (although 4.6
>> might help here).
>> My hypothesis is that factor 4 caused massive heating everywhere and
>> combined
>> with a strongly coupled thermostat, all effects a seriously dampened,
>> including those of 1-3.
>> Fixing factor 4 gives less heating, thus less thermostat coupling and
>> possibly less
>> uniform heating between water and protein (due to effects of 1-3 being
>> damped less),
>> causing hydrogens in the protein to overheat. But as there are so many
>> complicating factors,
>> it's difficult to find out what's really happening.
>> When I fixed factors 1-3, by different mdp settings, I get perfect
>> energy conservation in 4.5,
>> but that's at a high computational cost.
>> Cheers,
>> Berk
> Thanks for clarifying the different factors involved. I see what you
> mean. The relation between these factors is intricate and difficult to
> antecipate, and trying to fix one of them ended up making things
> unstable when using some setups (such as ours).

As I pointed out earlier, due to large inherent heating in your setup, 
the system becomes unstable, and one can delay the "fatal fluctuation" 
but not indefinitely.

> You mention that 4.6 might help with factors 1-3. Does that mean that
> you are making further changes? Assuming that we decide to stick with
> our current setup, would you then advise us to stay with 4.0 for now
> and wait for 4.6 to be released?

If Berk has implemented these things correctly (which he usually does) 
4.6 will make it more efficient to run with perfect energy conservation, 
but the inherent problems in your setup remain.

What most likely "saved" you in previous versions, is the twin-range 
scheme (1) where the forces from far-away molecules are kept constant 
even though atoms may  have  moved out of the cut-off or into the 
cut-off.  In addition, the RF with epsilon < infinity in the original 
implementation has a minimum in the energy slightly beyond the cut-off. 
This may also improve stability.

In summary: energy conservation is your friend!

> Cheers,
> Antonio
>> On 12/17/11 22:12 , Shirts, Michael (mrs5pt) wrote:
>>> A couple thoughts:
>>> 1. if something was stable before (even if it was not best
>>> practices), then
>>> I agree that it should be stable in the current code, unless the
>>> previous
>>> algorithm was wrong to begin with. I haven't looked into the details
>>> best
>>> to know what the problem is yet.
>>> 2. Massive heating damped with temperature control will drive you
>>> away from
>>> an equilibrium to a steady state where some properties become
>>> incorrect even
>>> if others are correct. So you do have to be very careful that the
>>> properties you looked at before are the same ones. Where exactly the
>>> failure point is is very hard to say very clear information about the
>>> specific system.
>>> 3.
>>>> We can also discuss if, as a general rule, using a continuum reaction
>>>> field is better or worse than using a lattice method such as PME,
>>>> but that
>>>> is a different issue. Like many other people (Wilfred van Gunsterem,
>>>> Arieh
>>>> Warshel, Alan Mark, Philippe Hünenberger, etc), I'm really not
>>>> convinced
>>>> that lattice methods offer any real advantage, because I find the
>>>> published evidences for either their benefits or their artifacs to be
>>>> rather weak or contradictory.
>>> For HOMOGENEOUS systems, then I agree; the configurations sampled can be
>>> achieved by cheaper methods than PME -- basically, the configurations
>>> are
>>> dominated by short-range effects, and beyond a certain range, all
>>> that the
>>> extra interactions do is affect the overall energy, which can be
>>> calculated
>>> pretty well by many methods (including reaction field).
>>> When lattice methods are needed is inhomogeneous systems, such as
>>> membranes
>>> and interfaces, because the long range order affects the configurations
>>> sampled, and RF and other continuum methods simply can't handle those
>>> well;
>>> properties will depend on cutoffs, etc.
>>> Best,
>>> ~~~~~~~~~~~~
>>> Michael Shirts
>>> Assistant Professor
>>> Department of Chemical Engineering
>>> University of Virginia
>>> michael.shirts at virginia.edu
>>> (434)-243-1821
>> --
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the www
>> interface or send it to gmx-developers-request at gromacs.org.

David van der Spoel, PhD, Professor of Biology
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se

More information about the gromacs.org_gmx-developers mailing list