[gmx-developers] Reaction Filed crash!

Sat Dec 17 21:47:51 CET 2011

On Fri, 16 Dec 2011, David van der Spoel wrote:

> On 2011-12-16 14:48, Miguel Machuqueiro wrote:
>> 
>> Hi again,
>> 
>> It seems that your test was too short for you to see the kind of crashes
>> we have been observing.
>> Right after the first tests we also decided to move to the "v-rescale"
>> thermostat. All my latest tests have been done with this and sometimes
>> with both (to make sure it was not affecting the results).
>> 
>> A new test using a nstlist=4 and a tau_t=0.05 resulted in a crash after
>> 2.07 ns. The system distorted on the main chain's NH of the last residue
>> of HEWL.
>
> I realize this is not what you want to hear but here goes.
>
> Your cut-off / electrostatics setup is extremely bad for energy conservation.
>
> First, with the exact same mdp file, but T-coupling turned off, the 
> temperature increases linearly by 100 K in 25 ps.
>
> Secondly, with reaction-field-zero and epsilon_RF = 0 I get exactly the same 
> result - this could point to a bug actually. In this case grompp issues many 
> notes and warnings, e.g.:
>
> WARNING 2 [file Lyso_50ns.mdp]:
>  The sum of the two largest charge group radii (0.512292) is larger than 
> rlistlong (1.400000) - rcoulomb (1.400000)
>
> NOTE 3 [file Lyso.top, line 8414]:
>  System has non-zero total charge: 8.999999
>
> If I were to referee a paper based on these premises I would require that 
> authors redo the work with a sane setup containing at least:
>
> - counter ions
> - PME
>
> Reaction fields are sometimes acceptable for dipolar systems if the cut-off 
> is quite a bit larger than the range of ordering in the system. You can guess 
> that water will be ordered significantly around a protein with a charge of +9 
> and no counter-ions.
>
> Therefore, if you do add the counter-ions and use PME you can reduce your 
> cut-off somewhat, use a single range cut-off and dispersion corrections, and 
> probably get the same simulation efficiency. In addition you will be able to 
> publish your results :).

Dear David

When discussing these issues, I think it's important to keep a clear 
distinction between the sampling scheme (which includes the integration 
algorithm) and the Hamiltonian effectively being used (which includes the 
long-range treatment of electrostatics).

The fact remains that the simulations using exactly the same setup (hence 
using the same Hamiltonian) don't crash with previous GROMACS versions. 
So, some recent changes in the sampling scheme have somehow created the 
problem (perhaps the new Trotter scheme for the twin-range interactions).

Concerning the Hamiltonian being used (including RF vs PME), the setup 
that Miguel and I are using didn't came out of the blue -- the suggestions 
you are making were actually previously tested and published by us and 
others. In particular, we compared constant-pH MD simulations using PME 
and RF with and without counterions, with negligible differences in terms 
of structure or (de)protonation. Sure, the cutoff-induced heating in 
constant-energy simulations (specially with highly charged systems) has 
been well-known for a long time, but that doesn't make it a mortal sin in 
a constant-temperature simulation where the properties of interest were 
tested and found to be unaffected.

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. Therefore, for now I prefer to bet on a 
method that inconsistently superimposes continuum concepts on an 
atomic-level description but nonetheless keeps an isotropic physical model 
akin of the solvated protein we are trying to simulate, rather than to bet 
on a method that sticks to the atomic-level description but whose 
underlying physical model is a highly-concentrated protein lattice with 
hardly any resemblance to a solvated protein. We all know that the 
long-range nature of electrostatic interactions is a bitch, but continuum 
models have proven quite useful dealing with it before (eg, Debye-Hückel 
theory avoids the diverging virial expansion terms for electrolytes), so 
maybe they can help us again in devising a decent reaction field for 
solvated molecules. I also know that most people are betting on the other 
horse, but on scientific issues I'd rather be a Bayesian than a 
frequentist. Time may prove me wrong, of course...

In short, we think that our approach is sound for the kind of studies we 
are doing, this being the reason why we are trying to convince you GROMACS 
guys to not give up from keeping the RF method working. GROMACS is a great 
tool that you offered us all and I think it would be a pitty if the 
possibility of using a RF treatment would go away. And, of course, we 
would be glad to answer any of your concerns if you happen to referee one 
of our papers... :)

Thank you for all the suggestions and help.

Best,
Antonio

>
>
>> 
>> If you care to take a look at my system, here goes the files needed to
>> reproduce my crashes:
>> 
>> http://dl.dropbox.com/u/3245083/Lyso_RF-test_GMX-4.5.x.rar
>> 
>> Once again, thank you for all the effort.
>> Cheers,
>> Miguel
>> 
>> 
>> On 13-12-2011 17:02, Berk Hess wrote:
>>>  Hi,
>>>
>>>  I did one more check. I can run a 128 residue protein for 100 ps with
>>>  both nstlist=4 and nstlist=5 without warnings. But I do see that the
>>>  protein is much more above the target temperature than the water.
>>>
>>>  One difference I did notice is that you are using a Berendsen
>>>  thermostat. I would always use the Bussi thermostat (v-rescale)
>>>  instead. With that thermostat you can make tau_t arbitrarily small.
>>>  So I would try tau_t=0.05
>>>
>>>  Cheers,
>>>
>>>  Berk
>>>
>>>  On 12/12/2011 09:27 PM, Miguel Machuqueiro wrote:
>>> >
>>> > Hi Berk, and all the devel guys.
>>> >
>>> > Once again, thank you for not giving up on us.
>>> >
>>> > It took me a while to do the tests to check your suggestions and my
>>> > results are not very famous.
>>> >
>>> > A decrease in the nstlist from 5 to 4 was not able to provide
>>> > stable MD simulations in my HEWL system. All attempts (several by
>>> > now) ended up in lincs crashed with several stepxxxxxx.pdb files
>>> > being generated. I also tried simulations using nstlist values of
>>> > "2" and "1" to be sure: - The system with a value of "2" held on
>>> > for a while, but eventually crashed after ~15 ns. I reproduced this
>>> > crash in 3 replicates. - The system with a nstlist of "1" went
>>> > through the 50 ns test without any problems. Unfortunately, this
>>> > simulation was very slow; 4 times slower than PME (remember that
>>> > nstlist=5 gives RF simulations much faster than PME).
>>> >
>>> > The idea that RF in previous versions of Gromacs was stable
>>> > probably due to error cancellation raises in me several concerns.
>>> > Unfortunately, with the correction in the current implementations,
>>> > I am not able to use the RF option for the treatment of the long
>>> > range electrostatics.
>>> >
>>> > If you are interested, I could send you the files used to run my
>>> > tests. Even though we modified our code in order to use GMX version
>>> > 4.0.7 (in use at the moment), we are still very interested in using
>>> > the latest and faster versions of Gromacs.
>>> >
>>> > Thank you and all the guys doing the hard work! Miguel
>>> >
>>> >
>>> >
>>> > On 27-11-2011 14:59, Berk Hess wrote:
>>> >> Hi,
>>> >>
>>> >> I found a solution for your problem!
>>> >>
>>> >> The basic issue is that a time step for the long range forces of
>>> >> 5*2 fs = 10 fs is too long in a system where charges move fast
>>> >> (probably mainly the water hydrogens). With the old,
>>> >> non-reversible, scheme there must be some cancellation of
>>> >> errors, which makes the integration stable (although still very
>>> >> inaccurate). With the new, reversible, integration scheme the
>>> >> hydrogens in the protein heat up by a factor 1.5 with a
>>> >> T-coupling strength of tau_t=0.1 (worse than the old scheme)
>>> >> and they become unstable. Changing the long time step to 8 fs
>>> >> (nstlist=4 with dt=2fs) reduces the energy drift by a factor of
>>> >> 4 and makes protein simulations stable. With this 8 fs the energy
>>> >> drift, for my protein test system, with the new scheme is a
>>> >> factor 3 lower than with the old scheme.
>>> >>
>>> >> The proper solution is using PME. But you said you couldn't use
>>> >> that for whatever reason. I hope that that is not because you
>>> >> need pair forces. Note that the reaction field force is NOT a
>>> >> pair force, it is a collective effect of all charges in the
>>> >> cut-off sphere, very similar to PME.
>>> >>
>>> >> Cheers,
>>> >>
>>> >> Berk
>>> >>
>>> >> On 2011-11-03 18:52, Miguel Machuqueiro wrote:
>>> 
>> 
>> 
>> --
>> ============================================
>> Miguel Machuqueiro
>> Department of Chemistry and Biochemistry
>> Faculty of Sciences, University of Lisbon
>> Campo Grande, Edifício C8 (sala 8.5.47)
>> 1749-016 Lisboa, Portugal
>> Tel. : +351 217500112 (int.ext.28547)
>> Mobile: +351 967562285
>> E-mail: machuque at fc.ul.pt
>> www1: http://webpages.fc.ul.pt/~mamachuqueiro
>> www2: http://intheochem.fc.ul.pt
>> ______________________________________________
>> 
>> 
>> 
>> 
>
>
> -- 
> David van der Spoel, Ph.D., Professor of Biology
> Dept. of Cell & Molec. Biol., Uppsala University.
> Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
> -- 
> 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.
>

-- 
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619         email: baptista at itqb.unl.pt
fax:   +351-214411277         WWW:   http://www.itqb.unl.pt/~baptista
--------------------------------------------------------------------------