Problem with Berendsen temperature and pressure control

Hi, all

I try to equilibrate the simulation cell with Berendsen temperature and pressure control. Then I use command like :

Hi, all

      I try to equilibrate the simulation cell with Berendsen temperature
and pressure control. Then I use command like :

-----------------

timestep 0.005

velocity all create 600 376137 dist Gaussian

fix 1 all nve

fix 2 all temp/berendsen 300.0 300.0 0.5

fix 3 all press/berendsen iso 0.0 0.0 5.0

run 30000

----------------

However, the results show that the cell expanded dramatically during
equilibration, which is not expected. Need your help to found the mistake

difficult to say with such limited information.

the mistake is likely due to a combination of three issues that are
not suitable handled:
- your initial geometry is bad or the density of your initial geometry
is too high or your potential parameters not suitable or simply
incorrect leading to very high initial pressure and potential energy
- you are using fix press/berendsen right away (better to run for a
while without volume adjustments and turn it on later, especially)
- you have too short a relaxation time (when starting from a very
high-pressure initial state, the simulation cell will expand too
quickly and then it will take a very long time until you recover an
equilibrium volume/density. this is a deficiency of the berendsen
algorithm by itself and due to the fact that expansion has no barrier
while compression requires atoms rearrangements)).

axel.

Also, what units are you using? 5.0 for press damping
would be way too short for real, but OK for LJ or metal.
You can also monitor thermo output (say every 5 steps
at the beginning) and see if the pressure is dramatically
high. If so it will push your box to be very large. And
it won’t stop until the pressure drops below your target.

Steve

I am using metal units with timestep 0.005

Here are the thermo output

Step Temp PotEng Pxx Pyy Pzz

0 600 -2105908.8 -27644.34 -29965.759 -45951.097

5 491.91554 0 -126.7345 -143.39485 -131.0186

10 482.50976 0 -133.10779 -150.60596 -137.60733

15 473.56497 0 -140.38791 -158.84311 -145.13354

20 465.05856 0 -148.76425 -168.32059 -153.79303

25 458.08825 -12.33533 -155.96773 -178.7153 -159.78392

30 29275.131 4493.3669 51135.045 7635.3225 47767.219

The potential energy turn to 0 just after the simulation begin, but it’s strange that the system can be equilibrated with Nose-Hoover thermostat and barostat(NPT), and without any expansion. here is the corresponding output.

Step Temp PotEng Pxx Pyy Pzz

0 600 -2105908.8 -27644.34 -29965.759 -45951.097

5 256.20591 -2086877.6 -5250.8517 -16162.74 -29422.906

10 243.97896 -2086188.4 -1162.3909 -30646.344 -22442.886

15 297.09628 -2089104.3 -9097.0808 -30328.588 -31436.992

20 303.07141 -2089498 -11048.486 -27708.848 -26192.648

25 317.21551 -2090339.8 -12832.395 -26882.382 -23973.068

30 294.44988 -2089128.1 -11712.894 -26817.628 -23501.67

35 306.03206 -2089825.8 -6337.311 -31135.224 -19820.944

40 292.39328 -2089109.9 -1245.6545 -22868.901 -19019.411

45 312.48164 -2090250 1698.4996 -15156.587 -16029.235

50 311.66776 -2090221.5 6118.1326 -10041.203 -12482.688

55 309.27726 -2090103.4 8666.5804 -11529.919 -13431.579

60 306.80278 -2089982.2 12174.344 -13811.524 -12095.79

65 305.64347 -2089930.9 11700.919 -11511.043 -12280.999

70 310.20745 -2090195.4 8621.3057 -9748.1374 -10872.36

75 308.75486 -2090129.2 5914.8068 -13052.945 -9548.9892

80 302.87402 -2089817.4 3816.4902 -17024.934 -10042.464

85 308.34288 -2090138.7 4032.2649 -18706.622 -9436.3035

90 301.17041 -2089764.8 5655.9022 -13865.015 -8308.9963

95 300.33894 -2089736.6 8181.1731 -9355.3081 -7089.3513

100 309.64968 -2090259.2 10779.477 -9907.8953 -6185.4319

105 311.44141 -2090368.1 13612.33 -10327.685 -6178.7306

110 310.00855 -2090307 16794.254 -7548.9427 -5214.3328

I am using metal units with timestep 0.005

Here are the thermo output

Step Temp PotEng Pxx Pyy Pzz

       0 600 -2105908.8 -27644.34 -29965.759 -45951.097

       5 491.91554 0 -126.7345 -143.39485 -131.0186

      10 482.50976 0 -133.10779 -150.60596 -137.60733

      15 473.56497 0 -140.38791 -158.84311 -145.13354

      20 465.05856 0 -148.76425 -168.32059 -153.79303

      25 458.08825 -12.33533 -155.96773 -178.7153 -159.78392

      30 29275.131 4493.3669 51135.045 7635.3225 47767.219

The potential energy turn to 0 just after the simulation begin, but it’s
strange that the system can be equilibrated with Nose-Hoover thermostat and
barostat(NPT), and without any expansion. here is the corresponding output.

this is not strange at all. as i already mentioned, there is a
systematic deficiency in berendsen algorithms, and your initial
configuration and settings exposes it.
nose-hoover uses a different mechanism than berendsen to manipulate
temperature and volume of your system, that is much more indirect and
those less likely to cause unphysical behavior. judging from the
output below, i still think that there is something not quite right.
your initial configuration and settings seems to indicate an extremely
incompressible medium and the pressure values hint at using an
isotropic volume relaxation for a system, that is anisotropic.

as mentioned many times before, MD follows the GI-GO principle (=
garbage in, garbage out). in other words, if you feed LAMMPS
unreasonable or unsuitable input, you must not expect reasonable
output. computer software does not reason, it blindly follows what you
ask it to do.

axel.

Something very bad happened to the barostat during the first 5 timesteps. Nothing after this point can possibly be right.

0 600 -2105908.8 -27644.34 -29965.759 -45951.097
5 491.91554 0 -126.7345 -143.39485 -131.0186

Clearly, you started the system out with the volume somewhat too large (P ~ -30kbar). After 5 steps, the volume has become huge (small negative pressure, PE=0). In between, the Berendsen barostat probably made the volume extremely small, which caused the subsequent explosion. In order to solve the problem, you need to bring the system to its equilibrium volume gently. Here are some ways to do this:

  1. Increase value of Pdamp by a lot
  2. Use the modulus keyword to enter the correct bulk modulus. (default value is 10.0 bar in metal units, which is probably many orders of magnitude too soft)
  3. Ramp the target pressure from the current value to zero over the duration of the simulaiton

Putting that all together, I suggest:

fix 3 all press/berendsen iso -30000 0.0 5.0 modulus 1000.0

Aidan