berendsen barostat issue

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I am having some issues implementing the berendsen barostat in my glycine water system. My initial attempt consisted of the following fixes

fix 1 all nve
fix 2 all langevin 270.0 270.0 2000.0 699483
fix 3 all shake 0.000001 20 100000 b 1 3 5 a 4 #constrain N-H C-H and O-H
fix 4 all press/berendsen iso 1.0 1.0 1000.0

However, within five time steps I actually got a negative volume and quickly a pppm atom out of range error popped up. The volume of the system doubled between steps 1 and 2. Clearly the dynamcs were very bad.

I did some digging on the email list, and found discussions where people said to first equilibrate in NVT. So I tried the following

fix 1 all nve
fix 2 all langevin 270.0 270.0 2000.0 699483
fix 3 all shake 0.000001 20 100000 b 1 3 5 a 4 #constrain N-H C-H and O-H
run 10000

fix 4 all press/berendsen iso 1.0 1.0 1000.0

run 10000

Now I know 10000 isn’t long enough from a true equilibration, but I think it should be long enough to get out and kinks that could produce instability numerics. However, again this blew up within 5 to 10 steps as before.

I then tried

fix 1 all nve
fix 2 all langevin 270.0 270.0 2000.0 699483
fix 3 all shake 0.000001 20 100000 b 1 3 5 a 4 #constrain N-H C-H and O-H
fix 4 all press/berendsen iso 1.0 1.0 1000000.0

Now this produces a stable numerical integration. However, from the examples I could find on the mailing list, my own experience, and the lammps doc example, that seems like an outrageously large pdamp. However, the volume of the system is still changing by 50A^3 between time steps, which is reasonable. Could anyone comment here? If I think about my knowledge of the berendsen algorithm, it looks something like

mu = 1 -(dt/tau)(ptarget-pressure)/nktv2p
mu = mu**(1.0d0/3.0d0)
box = boxmu

We would then also scale the particle coordinates. Am I incorrect to interpret the lammps pdamp parameter as tau in the above equation?

Also of note, fix 1 all npt temp 270.0 270.0 100.0 iso 1.0 1.0 1000.0 produces a stable numerical integration.

Input file and data file are attached.

glycine-water.in (2.62 KB)

data.glycine-test (801 KB)

2

First, turn off langevin and shake, just run nve with press/berendsen.
Print out thermo every tilmestep and compare the actual pressure
to the target pressure. Unless they are widely different with a bad
damping constant,
you should not get the a negative volume in 5 time steps.

Steve

3

Thanks for your response, Steve.

I have tried that already, however, to no avail. The starting pressure, at least from my experience, is a pretty reasonable pressure value of -2605.2055 atm. Thats what has me so perplexed.

4

If you print the pressure and volume for each of the first

5 steps, along with your damping constant and target pressure,

then it should be clear why the volume is changing as it is.

The Berendsen formula is fairly simple.

Steve

5

Maybe you could try to use the modulus keyword for press/berendsen. The Berendsen barostat also requires to know the (approximate) bulk modulus of your system: the default in LAMMPS is 10 pressure units. If you’re using real units, that leads to a modulus of 10 atm, which is way smaller (four orders of amgnitude) than what you’d want for water. Consequently, your volume changes will also be too large. Making the modulus larger actually corresponds with using a large pdamp value, consistent with your observation that a much larger pdamp value gives a more reasonable behavior.

Kristof

6

If Kristof’s explanation is not the answer, then you need to post

your press/berendsen command (and units) and the thermo

output for the first few steps (every step), that I asked for

previously.

Steve