C-rescale barostat

User discussions
, , ,
1

GROMACS version:2023.2
GROMACS modification: No

Hello,

I am working on a system with 700 glycerol molecules using OPLS FF. With the time-step of 1 fs, I am running NPT simulation using C-rescale (tau_p= 3ps) barostat with a reference pressure of 1 bar and NH (tau_t=0.1) thermostat with a reference temperature of 300 K.

Here is the
C-r.mdp (3.4 KB)
file for reference.

Reading this post from GitLab, where @hess
explained the approach to choose tau values https://gitlab.com/gromacs/gromacs/-/merge_requests/2624 and the thread of this discussion, I understood that tau_t has to be a smaller value compared to tau_p.

The problem with my simulation is that once the system reaches the reference pressure, it attains the desired density but on running the simulation for a longer time using the same barostat (C-rescale) to reach equilibrium, the density decreases. I am assuming that the configuration of the system is not at equilibrium and is showing this behaviour to be at equilibrium but is this normal?

I have tried testing my hypothesis with different values of tau_p and tau_t to see the impact of oscillation but the result remains persistent.

@hess I tested the simulation with tau_t= 1ps and tau=5 ps for V-rescale thermostat and C-rescale barostat (as per this thread) but the density hiked up drastically and then decreased eventually.

Could you please suggest to me if my approach is wrong for using the barostat?

1 Like
2

You are now generating velocities in your NPT stage. Have you run an NVT equilibration before? I would recommend doing that and set gen-vel = no. I think tau_t=1ps and tau_p=5ps are reasonable. But you might want to try running it in separate stages, e.g.:
NVT with tau-t=0.1 (perhaps 0.5 - 1 ns)
NPT with tau-t=0.2 and tau-p=1 (perhaps 1-2 ns)
NPT with tau-t=1 and tau-p=5 (2-50 ns, depending on when it is stable)

By the way, I just noticed from your mdp file that you are writing uncompressed and compressed trajectory output very often. Do you even take a look at force output at all? Do you use velocity output? Unless you’ve got very large storage capacity, I’d recommend something along the lines of:

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 0
nstvout                  = 0
nstfout                  = 0
; Output frequency for energies to log file and energy file
nstlog                   = 10000
nstenergy                = 100000
; Output frequency and precision for xtc file
nstxtcout                = 100000

But that is just my recommendation. There can be all sorts of valid reasons for outputting more often than that, as well as uncompressed positions and velocities (e.g. if you are generating starting configurations for umbrella sampling). But during equilibration you usually only need such detailed output if you need to analyze crashes.

3

Dynamics of pure glycerol is extremely slow; it’s viscosity is a factor 140 higher than that of water. How fast things equilibrate will depend on how close to equilibrium your initial configuration is.