Use of Berendsen barostat for reducing the pressure fluctuations (equilibration)

Dear All,

I am trying to equilibrate a box of water having 1000 molecules. I was able to maintain the energy and temperature at almost constant values but the pressure fluctuations were tremendously high (±2000 atm). I had used fix npt.

To reduce the fluctuations, I used fix press/berendsen after fix npt which reduced the fluctuations to a large extent (±12 atm) but the average value is different from what I gave in npt as well as in press/berendsen command. So, can I say that the system has been equilibratred and start with the production run or it is still not equilibrated?

After running fix press/berendsen and observing the fluctuations, I used fix nve and observed the pressure values. The large pressure fluctuations (same as before ±2000) appeared again.

So I am not able to understand whether I should start the production run immediately after using fix press/berendsen command or try something else.

The code is as below:

Initialization

boundary p p p
units real
atom_style full
read_data data.step5a

Dreiding potential information

pair_style lj/charmm/coul/long 9.0 10.0 10.0
pair_coeff 1 1 0.00 2.058
pair_coeff 1 2 0.00 2.612
pair_coeff 2 2 0.1553 3.166
bond_style harmonic
bond_coeff 1 450 0.9572
angle_style harmonic
angle_coeff 1 55 109.47
kspace_style pppm 0.0001
pair_modify mix arithmetic

neighbor 1 bin

neigh_modify delay 0 every 1

min_style cg
thermo_style custom step etotal temp vol press density

minimize 1e-25 1e-25 500000 1000000

velocity all create 300.0 4928459 rot yes mom yes dist gaussian

fix 1 all nve
thermo 1000
thermo_style custom step etotal temp vol press density
timestep 0.25
run 50000

unfix 1

fix 1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
timestep 0.25
thermo_style custom step etotal temp vol press
thermo 1000
run 200000

unfix 1

fix 1 all nve

fix 2 all press/berendsen iso 1.0 1.0 1000.0
thermo_style custom step etotal pe temp vol press density
thermo 1000

timestep 0.01
run 200000

unfix 1
unfix 2

fix 1 all nve
timestep 0.25
thermo_style custom step etotal temp vol press density
thermo 1000

run 200000

If your system is only 1000 molecules, I am not that surprised you see large pressure fluctuations. You might be able to reduce them by playing with the damping time, but it might lead to an unstable simulation. I do think the average pressures of Berendsen and Nosé-Hoover should be the same, but Berendsen does not sample the correct distribution so the fluctuations can be different.

Adding to what Stefan said, using berendsen thermo/baro stat is not recommended as they are giving you an NVT/NPT ensemble. Having said that, with the state of the water you specified (300K, 1bar) pressure fluctuation can be quite large as it’s an incompressible system, have looked at the density values? are they in the range for the water model you’re using?

kasra.

correction: “They are NOT giving you an NVT/NPT ensemble”

Yes the density values are in the range for the model of water used for the simulation.

Regards,
Ajinkya Sarode,
PhD- Mechanical Engineering,
IIT Gandhinagar.

Dear All,

I am trying to equilibrate a box of water having 1000 molecules. I was able
to maintain the energy and temperature at almost constant values but the
pressure fluctuations were tremendously high (+-2000 atm). I had used fix
npt.

To reduce the fluctuations, I used fix press/berendsen after fix npt which
reduced the fluctuations to a large extent (+-12 atm) but the average value
is different from what I gave in npt as well as in press/berendsen command.
So, can I say that the system has been equilibratred and start with the
production run or it is still not equilibrated?

equilibration and magnitude of fluctuations are two completely
different things. you will only get a correct statistical mechanical
ensemble, if the pressure coupling (or temperature coupling for that
matter) is *weak*. thus if you dampen fluctuations noticeably, then
you modify the distribution of states and will get a wrong ensemble.

your equilibration and simulation times, however, are horribly short
(a few 10s of picoseconds), depending on how far you were off
equilibrium in your initial conditions, you may still be far off. as
for how to determine whether you have reached equilibrium, you should
consult your text book on MD. it should be discussed there (or you
need to get a better one). for this short trajectories, you may not
get well converged data.

axel.