Project

General

Profile

Bug #2147

Parrinello-Rahman barostat not properly working

Added by Manuel Martinez over 2 years ago. Updated over 2 years ago.

Status:
Feedback wanted
Priority:
Normal
Assignee:
-
Category:
-
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Hi! I am posting here after widely looking for a solution by myself but not finding one. Once I get a successful equilibration using Berendsen barostat and v-Rescale termostat (p=1 bar, T=300 K), I turned to P-R and N-H couplings for data collection. While N-H thermo works fine, the Parrinello-Rahman barostat provides a greater pressure than the one desired (p~1.5 bar instead of p=1 bar). Is this a bug for recent Gromacs versions? Or is just a setting parameter for mdp file I can't figure out? Thanks in advance.

Manuel

agua.mdp (11.5 KB) agua.mdp Manuel Martinez, 03/23/2017 10:19 PM
agua.top (1.01 KB) agua.top Manuel Martinez, 03/23/2017 10:19 PM
conf_relaj.g96 (238 KB) conf_relaj.g96 Manuel Martinez, 03/23/2017 10:19 PM
grompp.mdp (3.07 KB) grompp.mdp Manuel Martinez, 03/23/2017 10:26 PM
nacl_aq.top (1.45 KB) nacl_aq.top Jeff Thompson, 04/15/2017 04:04 PM
sd_4.6.5.mdp (956 Bytes) sd_4.6.5.mdp Jeff Thompson, 04/15/2017 04:04 PM
sd_5.1.2.mdp (960 Bytes) sd_5.1.2.mdp Jeff Thompson, 04/15/2017 04:04 PM
nacl_aq.gro (1.27 MB) nacl_aq.gro Jeff Thompson, 04/15/2017 04:04 PM

History

#2 Updated by Manuel Martinez over 2 years ago

I used Gromacs 5.0.4.

#3 Updated by Mark Abraham over 2 years ago

It can take a long simulation to make a converged measurement of pressure, given its fluctuations. Can you describe/share how you made your observations?

#4 Updated by Jeff Thompson over 2 years ago

It does seem that there is a difference in how Parinello-Rahman pressure coupling is being implemented in version 5.1.2 vs. 4.6.5.

The attached .mdp files run 12 ns stochastic dynamics at 298.15 K and 1 bar with Nose-Hoover temperature coupling and Parinello-Rahman pressure coupling. They are identical except for one line accounting for the renaming of 'verlet_buffer_drift' to 'verlet_buffer_tolerance' in version 5*. I've also included my topology and configuration files (for a model of an aqueous NaCl solution).

With version 4.6.5, I get an average density of 1.02258(4) g/cm^3 (discarding the first 2 ns as equilibration).

With version 5.1.2, the density drops by ~1% to 1.01180(4) g/cm^3.

I haven't yet had a chance to pick through the code to see what might've changed between the two versions, but it seemed worth sharing this observation here.

#5 Updated by Mark Abraham over 2 years ago

Jeff Thompson wrote:

It does seem that there is a difference in how Parinello-Rahman pressure coupling is being implemented in version 5.1.2 vs. 4.6.5.

The attached .mdp files run 12 ns stochastic dynamics at 298.15 K and 1 bar with Nose-Hoover temperature coupling and Parinello-Rahman pressure coupling. They are identical except for one line accounting for the renaming of 'verlet_buffer_drift' to 'verlet_buffer_tolerance' in version 5*. I've also included my topology and configuration files (for a model of an aqueous NaCl solution).

With version 4.6.5, I get an average density of 1.02258(4) g/cm^3 (discarding the first 2 ns as equilibration).

With version 5.1.2, the density drops by ~1% to 1.01180(4) g/cm^3.

I haven't yet had a chance to pick through the code to see what might've changed between the two versions, but it seemed worth sharing this observation here.

There's not enough information here to even speculate that the coupling algorithm is the issue. Or whether the issue is a fixed bug or an introduced one. A different force field parameter could produce this. What does gmx check report for your tpr files?

#6 Updated by Jeff Thompson over 2 years ago

Mark Abraham wrote:

Jeff Thompson wrote:

It does seem that there is a difference in how Parinello-Rahman pressure coupling is being implemented in version 5.1.2 vs. 4.6.5.

The attached .mdp files run 12 ns stochastic dynamics at 298.15 K and 1 bar with Nose-Hoover temperature coupling and Parinello-Rahman pressure coupling. They are identical except for one line accounting for the renaming of 'verlet_buffer_drift' to 'verlet_buffer_tolerance' in version 5*. I've also included my topology and configuration files (for a model of an aqueous NaCl solution).

With version 4.6.5, I get an average density of 1.02258(4) g/cm^3 (discarding the first 2 ns as equilibration).

With version 5.1.2, the density drops by ~1% to 1.01180(4) g/cm^3.

I haven't yet had a chance to pick through the code to see what might've changed between the two versions, but it seemed worth sharing this observation here.

There's not enough information here to even speculate that the coupling algorithm is the issue. Or whether the issue is a fixed bug or an introduced one. A different force field parameter could produce this. What does gmx check report for your tpr files?

You are absolutely right--definitely a lapse in logic on my part. It seemed there was little reason to suspect the force field parameters, as the top file was self-contained. But gmx check alerted me to the fact that the default sd integrator changed between the two versions. I'm checking what happens if I revert to sd2 in v5.1.2.

In addition, some data from NVT simulations suggest that it may be the pressure calculation, not the coupling algorithm, that changed. I realize that the different behavior could represent a fixed bug rather than an introduced one. I have some digging to do to figure out which one it is.

#7 Updated by Berk Hess over 2 years ago

  • Status changed from New to Feedback wanted
  • Priority changed from High to Normal

Did you find the cause of the difference?

#8 Updated by Erik Lindahl over 2 years ago

Unless we get any feedback the next few days, we'll close this bug.

#9 Updated by Jeff Thompson over 2 years ago

I haven't yet found the source of the difference, though, on empirical grounds, it appears to be unrelated to the choice of integrator (sd1 vs. sd2).

The ~1% difference in density remains problematic from the standpoint of my particular research problem. But this problem is a bit tangential to my work at the moment, and I can't commit to trying to resolve this issue in the short-term.

Since it is still not clear whether the issue is indicative of a bug or a bugfix, and given the lack of (1) continued feedback from the OP and (2) demonstration of a direct connection between my issue and that of the OP, it may make sense to close this bug.

#10 Updated by Manuel Martinez over 2 years ago

Hi everybody

I very sorry for not writing before and the lack of feedback; I was not aware of this discussion and at the same time I consider not having new valuable news. If still useful, I can say these observations:

I also found the same behavior as Jeff about density. After running several simulations, the anomalous behavior of Parrinello-Rahman barostat seems to be related to the following mdp parameters, when using Verlet cut-off scheme: nstpcouple, tau_p, nstcalcenergy and nstenergy (which in some way is obvious :D). I tried different reasonable combinations for these options observing notorious differences but not obtaining satisfactory pressure coupling for 1 bar (NH thermal coupling works fine). The "best run" I got used these:

nstlist = 50
nstpcouple = 10
tau_p = 2.5
nstcalcenergy = 50
nstenergy = 50
nstcomm = 50

I guess there is an optimal set of values for them but it can not be straightforward to find. For my current research I can use Berendsen barostat as an alternative but I will be attentive if more information arises. Thanks again.

Regards, Manuel

Also available in: Atom PDF