Project

General

Profile

Bug #1496

Langevin sd integrator: wrong average temperature

Added by Andrey Frolov over 3 years ago. Updated over 3 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Category:
mdrun
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Dear GMX-developers,

I perform MD simulation of one point on the phase space of supercritical CO2. I use NVT ensemble with the reference T=308.15K, number molecules 548, the box size is 3.6^3 nm^3. I use different thermostats: Berendsen, v-rescale, Nose-Hoover and Langevin (sd). The problem is that "sd" integrator gives a slightly elevated temperature (by ~ 3K) compared to the reference one. Also, the potential energy is higher for simulations with "sd" integrator. This strongly depends on the integration timestep, the smaller the timestep the lower the deviation. However, it seems to be independent of tau_p: tau_p = 2 and 5 ps produce similar elevated average temperature.

I should say that for the simulations of bulk tip3p water at ambient conditions I observe similar behavior (data not shown): the temperature tends to be slightly elevated with "sd" integrator (less than by ~ 1K, comparable to the statistical error)

Total simulation time for each parameter set is 5 ns, the averages for the last 4 ns:

integrator,tcoupl,tau_t[ps],dt[ps] Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
md,ber,2.0,dt=0.002 Temperature 308.2 0.0 3.4 0.0 (K)
md,nose,2.0,dt=0.002 Temperature 308.1 0.0 6.1 0.0 (K)
md,vrescale,2.0,dt=0.002 Temperature 308.1 0.1 5.9 0.6 (K)
sd,no,2,dt=0.0005 Temperature 308.4 0.2 6.2 -0.3 (K)
sd,no,2,dt=0.001 Temperature 308.8 0.1 6.2 -0.4 (K)
sd,no,2,dt=0.002 Temperature 311.5 0.1 6.1 0.6 (K)

gromacs_4.6.5 (double precision).
Also, I should note that the compressibility of the fluid is about an order higher than that of water. (~ 7.8e-4 bar^-1).

Please, see the tpr,top,gro,mdp files for the case of: sd,no,2,dt=0.002
Please, note that there are no constraints in these simulations.

Thank you very much.
Andrey Frolov
----
Junior researcher
Institute of solution chemistry
Russian Academy of Sciences

nvteq_sd_taut2.tpr (79.7 KB) nvteq_sd_taut2.tpr Andrey Frolov, 05/06/2014 10:27 PM
nvteq.mdp (6.38 KB) nvteq.mdp Andrey Frolov, 05/06/2014 10:34 PM
topol.top (1.48 KB) topol.top Andrey Frolov, 05/06/2014 10:34 PM
nvteq.gro (111 KB) nvteq.gro Andrey Frolov, 05/06/2014 10:34 PM

Related issues

Related to GROMACS - Bug #1418: Parrinello-Rahman barostat: temperature and potential energy depend on tau_pClosed2014-01-08

Associated revisions

Revision b8a3738b (diff)
Added by Nicolae goga over 3 years ago

Default SD integrator is now sd1 by Nicu Goga

The default SD integrator is now sd1, sd2 is deprecated. grompp
and mdrun both warn about sd2, and recommend sd.

The new sd1 integrator was developed by Herman Berendsen and
implemented by Nicu Goga.
Without constraints the sd1 integrator is nearly unchanged.
With constraints the new sd1 integrator is now as accurate as normal
MD with a leap-frog integrator, whereas the sd2 integrator will give
a slightly too high temperature.

Fixes #1496

Change-Id: I02b4537b05742ed499a84b625c4d4bf8994c0304

History

#1 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1496.
Uploader: Berk Hess ()
Change-Id: I02b4537b05742ed499a84b625c4d4bf8994c0304
Gerrit URL: https://gerrit.gromacs.org/3419

#2 Updated by Berk Hess over 3 years ago

Szilard Pall noticed a temperature issue with the SD integrator. This is probably the same issue. Could you try the new sd integrator that I put in patch:
https://gerrit.gromacs.org/#/c/3419/
and report back if this fixes the issue or not?

#3 Updated by Berk Hess over 3 years ago

I think the issue is that the angle, and possibly also the bond, potential is too stiff for a 2 fs time step. Plain MD will integrate such harmonic potential "too" well, and you will get perfect energy conservation, but your effective temperature will be higher than what the kinetic energy says. SD does not have this issue and will show you the real temperature.

#4 Updated by Andrey Frolov over 3 years ago

Dear Berk,

I am afraid I do not fully understand your argument. Please, correct me if I understood it wrongly.

The integration timestep is too big for MD to sample correctly the bonded terms of CO2: bending and stretching. The integration errors should result in a continuous increase of kinetic energy of atoms (right?) and the temperature will start to rise. The thermostats will react to this temperature rise in a different way.

In Langevin dynamics the strength of thermostat is not related to the temperature deviation from its reference value and is determined only by the friction coefficient. Therefore, if the friction coefficient is not extremely big, then it might happen that the strength of the thermostat is not high enough to remove the artificial heating associated with integration errors. From one side there is a constant increase of kinetic energy due to numerical errors, on the other side there is a constant decrease of kinetic energy due to exchange with an external bath via the stochastic term of Langevin equation. At some point a stationary state is achieved at a slightly elevated temperature.
I think that Anderson thermostat would suffer from the same problem, right?

Berendsen or Bussi-Parrinello thermostats will scale the velocities to match the target temperature, and what is important is that the thermostat will react proportionally to the deviation of temperature from the target value (T - T0): the bigger the deviation the stronger is the reaction of thermostat. Therefore, these thermostats can ensure that the average temperature corresponds (or almost corresponds?) to its target value.

Does it sound reasonable?

PS: Szilárd, Thank you for finding the problem with the bad energy conservation! Indeed, I did not notice that mdrun changed the nstlist from 1 to 40. I agree that more warning messages would be helpful in this case. Stdout of mdrun prints only this note, which does not state that it is going to change the nstlist setting:
"
With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
that with the Verlet scheme, nstlist has no effect on the accuracy of
your simulation.
"
The note about the change of nslist appears only in log file:
"Changing nstlist from 1 to 40, rlist from 0.9 to 0.9"

I think that it would be helpful to have this message in mdrun stdout also or not to allow automatic change of nstlist option.

PPS: I think that all other requests for my input files are no longer important, please, let me know if it is not the case.

Thank you.

#5 Updated by Andrey Frolov over 3 years ago

I am sorry for messing up things, I just found that the PS of the previous post is not related to this redmine issue. I repost this to "gerrit".

#6 Updated by Szilárd Páll over 3 years ago

Andrey Frolov wrote:

PS: Szilárd, Thank you for finding the problem with the bad energy conservation! Indeed, I did not notice that mdrun changed the nstlist from 1 to 40. I agree that more warning messages would be helpful in this case. Stdout of mdrun prints only this note, which does not state that it is going to change the nstlist setting:
"
With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
that with the Verlet scheme, nstlist has no effect on the accuracy of
your simulation.
"

That's printed by grompp, not mdrun!

The note about the change of nslist appears only in log file:
"Changing nstlist from 1 to 40, rlist from 0.9 to 0.9"

I think that it would be helpful to have this message in mdrun stdout also or not to allow automatic change of nstlist option.

I'm pretty sure this appears on the stdout too! Additionally, I see no good reason for disabling the automatic nstlist changing by mdrun. Thanks to the automated buffer estimation, for a given tolerance set by verlet-buffer-tolerance, mdrun is free to use any nstlist value that introduces no more drift than specified. One just needs to get used to the fact that with the Verlet scheme you should rely on verlet-buffer-tolerance to control the drift (instead of guesstimating the appropriate nstlist) or set it to -1 for full manual control of rlist and nstlist.

A patch that Berk pushed up yesterday does a minor tweak, though, it disables nstlist changing when nstlist=1 is set in the mdp.

#7 Updated by Berk Hess over 3 years ago

I think the nstlist/rlist message current only appear on stdout with mdrun -v.

Andrey, could you check that the too stiff angle (and/or bond) potential with respect to the time step is indeed the cause of this issue?
I would suggest to run the leap-frog integrator with time steps 0.5, 1.0 and 2.0 fs and check the angle and bond energy. I would expect the angle energy with dt=2 fs to be too high by a percentage similar to the temperature error in SD.

#8 Updated by Andrey Frolov over 3 years ago

Indeed, the bond and angle potentials are too stiff. See the bond and angle potentials below. Energy stabilizes at dt 0.0005 ps.

The increase of Bond and Angle energies at 0.002 ps is ~3 %, while temperature increase ~ 1 %: the values are comparable.

Integrator,tcoupl,tau_t[ps],dt[ps] Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------
sd,no,2,000125 Bond 1405.6 0.9 60.0 -5.1 (kJ/mol)
sd,no,2,00025 Bond 1406.6 0.9 60.5 -0.2 (kJ/mol)
sd,no,2,0005 Bond 1408.6 0.5 60.2 -0.6 (kJ/mol)
sd,no,2,001 Bond 1413.3 1.4 59.8 -2.4 (kJ/mol)
sd,no,2,002 Bond 1441.3 1.1 61.7 -0.5 (kJ/mol)

sd,no,2,000125 Angle 1406.0 1.1 61.0 -2.8 (kJ/mol)
sd,no,2,00025 Angle 1403.3 0.9 60.7 0.3 (kJ/mol)
sd,no,2,0005 Angle 1404.4 2.2 58.4 -6.6 (kJ/mol)
sd,no,2,001 Angle 1414.0 0.9 60.9 -1.6 (kJ/mol)
sd,no,2,002 Angle 1437.9 0.7 61.4 0.3 (kJ/mol)

sd,no,2,000125 LJ -4059.5 1.1 47.3 -5.0 (kJ/mol)
sd,no,2,00025 LJ -4058.6 0.4 47.4 1.2 (kJ/mol)
sd,no,2,0005 LJ -4058.6 0.6 47.4 -0.5 (kJ/mol)
sd,no,2,001 LJ -4058.6 0.5 47.4 0.1 (kJ/mol)
sd,no,2,002 LJ -4057.4 0.7 47.3 1.7 (kJ/mol)

#9 Updated by Berk Hess over 3 years ago

  • Status changed from New to Closed
  • Assignee changed from Szilárd Páll to Berk Hess

I'm closing this issue.
Might issue #1418 be caused by the same problem?

#10 Updated by Andrey Frolov over 3 years ago

Yes, I think that the too issues are caused by the same problem. So, the issue #1418 can be safely closed as well. Thank you very much for your time and help.

#11 Updated by Rossen Apostolov over 3 years ago

  • Related to Bug #1418: Parrinello-Rahman barostat: temperature and potential energy depend on tau_p added

Also available in: Atom PDF