Simulation halts at low temperature in NPT

LAMMPS LAMMPS Mailing List Mirror
1

Hi William, Axel,

Unfortunately I run into problems with the simulation halting after just a few timesteps and the volume either exploding or collapsing (and energies exploding).

By some coincidence, I spent the last days investigating the exact same phenomenon doing NPT in metallic systems (B~140GPa) with MEAM potential. Cell geometry in the reference structure (which is always a stable state in MEAM) was fully relaxed to press->0 using minimize and box/relax before.

The problem seems to come from the Nose-Hoover dynamics itself, which do not work at low temperatures. The barostat mass is taken as kT/w^2 (see also MTK1994, Section II.F). At low temperatures, that mass vanishes and the barostat has no inertia. Therefore, even the also very small thermal pressure is enough to fully "explode" an otherwise stable system unless one uses absurdly small dt.

For example, consider this test:
- take relaxed structure (Etot = -3057.8151 eV)
- equilibrate at 30K, p=1bar
- step down to 25, 20, 15, 1, 5, 1K
- record ensemble averages of quantities in stable state after each step

== With temp/berendsen, press/berendsen, nve integrator:
T p H
30.00 0.90 -3051.10
25.00 0.99 -3052.22
20.00 0.96 -3053.34
15.00 1.02 -3054.46
10.00 1.03 -3055.58
5.00 1.01 -3056.69
1.00 1.01 -3057.59

== With npt:
T p H
30.01 1.08 -3051.10
25.01 1.20 -3052.22
20.03 1.22 -3053.33
15.02 1.24 -3054.45
10.00 11.60 -3055.54
  5.51 4387.24 -62.44
   # simulation aborts due to broken numerics

I can't really think of anything that could be done about that, except manually scaling the barostat time constants for low temperatures to "fake" a larger inertia... so I tried doing that with a bit of guessed parameters:

   tauscaled = taubarostat * (1 + exp((25 - Tnext) * 0.15))
   lmp.command(f"fix int all npt temp {Tnext} {Tnext} {tauthermostat} iso {press} {press} {tauscaled} nreset 2000")

== With npt with scaled tau_press:
T p H
30.00 1.01 -3051.10
25.02 1.06 -3052.22
20.02 1.08 -3053.34
15.08 1.25 -3054.44
10.02 0.94 -3055.57
5.00 1.28 -3056.69
1.00 0.94 -3057.59

It's a bit arbitrary, but it works and I don't have to give up the NPT ensemble. Maybe you can do something similar for your system?

Best,
Sebastian

2

Hi Axel, Sebastian and everyone else,

Thanks Axel for catching the wrong geometry! I (recklessly perhaps) reused part of an earlier mailers geometry for the NaCl structure and noticed that the Cl atoms in the basis weren't at (a/2 a/2 a/2) but I wasn't sure if it would make a difference. Your explanation for what could occur makes a lot of sense.

Sebastian, your damping constant scaling looks intriguing.

I tried the following code to run after the cooling phase (dt = 1.0 initially):

fix 1 all temp/berendsen 0.01 0.01 \(100\.0\*dt\) \# Thermostat keeps temp at 0\.01 K fix 2 all press/berendsen iso 1\.0 500000\.0 (1000.0*dt) # Raise pressure from 1kbar to 500kbar
fix 3 all nve

run 1200000

Unfortunately I get an error message instantly upon starting the compression phase (step 10k + 1):

ERROR on proc 0: Out of range atoms - cannot compute PPPM (../pppm.cpp:1940)
Last command: run 1200000

According to the following quote of the error page this might be caused due to too high a timestep or that it might be fixed by changing the neigh_modify command:

"Out of range atoms - cannot compute PPPM
One or more atoms are attempting to map their charge to a PPPM grid point that is not owned by a processor. This is likely for one of two reasons, both of them bad. First, it may mean that an atom near the boundary of a processor’s sub-domain has moved more than 1/2 the neighbor skin distance without neighbor lists being rebuilt and atoms being migrated to new processors. This also means you may be missing pairwise interactions that need to be computed. The solution is to change the re-neighboring criteria via the neigh_modify command. The safest settings are “delay 0 every 1 check yes”. Second, it may mean that an atom has moved far outside a processor’s sub-domain or even the entire simulation box. This indicates bad physics, e.g. due to highly overlapping atoms, too large a timestep, etc."

I decreased the timestep from 1.0 to 0.1 but still got the same error so I added
neigh_modify delay 0 every 1 check yes

keeping the neighbor setting at default, I'm not sure this should change anything however since afaik the above command is the default (?). In any case I still got the same error.

Thinking that the issue might arise due to a too fast cooling phase I increased the length from 10k steps to 100k steps, no change however. There was also no change if I swapped out fix temp/berendsen with fix langevin. I tested changing the fix npt in the cooling phase with a temp/berendsen + press/berendsen + fix nve-combo (cooling from 10 K to 0.01 K while either increasing pressure to 1kbar or keeping it constant at 1bar) but this caused the simulation to crash at step 1 instead of step 10k + 1. I guess my potential could be the culprit but as mentioned earlier it doesn't look like it should.

I'm at a loss for why these combination of fixes crashes the simulation instantly when applied. I've tried searching the mailing list for this error but didn't find any other solutions than the ones I tried above.
Changing the kspace_style from pppm to ewald causes a similar hang (with berendsen thermo-/barostats) as the npt fix after the third or fourth step:

Step TotEng E_pair E_long Enthalpy Temp Press Volume
       0 -47391.457 -47543.777 -36963.937 -47716.572 100 -1941.5305 11481.993
       1 -47392.759 -47543.651 -36963.91 -47718.53 99.063181 -2415.6346 9247.0868
       2 -46246.597 -46397.105 -37086.479 -33894.261 98.810861 9015.7812 93944.172
       3 -26003.574 -26153.164 -25980.669 -34629.404 98.208061 -17001.422 34788.826

Does anyone have any ideas?

I will try Sebastians solution in the meantime.

Best regards,
William

3

Hi again,

I just wanted to give an update that using Sebastians method of the scaled Pdamp was successful.

These were the lines of code in question, in my input file:

variable dt equal 1.0
variable tauscaled equal \(1000\.0\*dt\)\*\(1\+exp\(\(25\-{epsilon})*0.15))
variable Tstart equal 10.0
variable epsilon equal 0.01
variable Pstart equal 1.0
variable Pend equal 500000.0

fix cooling all npt temp \{Tstart\} {epsilon} \(100\.0\*dt\) iso 1\.0 1\.0 (1000.0*dt)
run 10000
unfix cooling
fix compression all npt temp \{epsilon\} {epsilon} \(100\.0\*dt\) iso {Pstart} \{Pend\} {tauscaled}
run 1200000

Because my "Tnext"(epsilon = 0.01) variable never changed I could've simply used a fixed Pdamp constant (of the order 10^4, I haven't tried different Pdamp in this range). It was surprising, to my naive understanding, how well the compression proceeded even with such a high Pdamp (P,V,H increased smoothly).

Thank you so much, Axel and Sebastian, for all the help!

Best regards,
William

1 Like