[gmx-developers] Reaction Filed crash!

[David van der Spoel]

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 :).


>
> 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