Simulation halts at low temperature in NPT

Hi everyone,

I’m currently trying to pressurize a NaCl crystal in an NPT ensemble at close to 0 K (I’ve used T = 0.1 K) from p = 1 bar to 600 kbar to record and plot the enthalpy vs pressure.

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

At first I tried to simply start at T = 0.1 K and increase the pressure 1 bar per fs (my timestep is 1 in real units) for 600k steps but after the simulation halted after just a few timesteps I guessed that this was due to the large pressure change per timestep so I ran the simulation again with 1.2 million time steps but had the same thing happen. I then tested different temperatures and noticed that I could run the simulation at temperatures above T = 5 K so I cooled the system from T = 10 K to T = 0.1 K for 10k timesteps and then ran it for 10k timesteps at T = 0.1 K and p = 1bar (not sure this was necessary or made sense) and then increased pressure from p = 1 bar to p = 600 kbar for 1.2 million steps. Unfortunately the system hangs at time step 5 of the equilibration phase.

Here is my input file (Lammps version is: 64-bit 15Apr2020)

/// DEFINITIONS ///

dimension 3
units real
boundary p p p
atom_style charge

/// CRYSTAL STRUCTURE ///

Na atoms

lattice fcc 5.646 origin 0 0 0 # Lattice for Na system
region box block 0 4 0 4 0 4
create_box 2 box
create_atoms 1 box

Cl atoms

lattice fcc 5.646 origin 0.5 0 0 # Lattice for Cl system
create_atoms 2 box

Other data

mass 1 28.990
mass 2 35.453
set type 1 charge +1.0
set type 2 charge -1.0

/// POTENTIALS ///

pair_style born/coul/long 7.6 15.0
pair_coeff 1 1 9768.28 0.317 0 24.187 -11.52 # Na-Na interactions
pair_coeff 1 2 28937.74 0.317 0 161.21 -161.21 # Na-Cl interactions
pair_coeff 2 2 80367.73 0.317 0 1669.6 -3353.63 # Cl-Cl interactions

kspace_style ewald 1.0E-5

variable dt equal 1
variable epsilon equal 0.1
variable Pstart equal 1.0
variable Pend equal 500000

timestep ${dt}

Energy minimization

minimize 1.0e-10 1.0e-10 100000 100000

Simulation

velocity all create 10.0 12345 mom yes rot no dist gaussian # Take initial velocity from gaussian distribution at T = 10K
fix cooling all npt temp 10.0 {epsilon} (100.0dt) iso {Pstart} {Pstart} $(100.0dt) # Cooling at 1bar from 10 K to 0.1 K

thermo 1

thermo_style custom step etotal enthalpy temp press vol

#dump 1 all custom 10 NaCl_B1.out id xs ys zs

run 10000

unfix cooling

fix equilibration all npt temp {epsilon} {epsilon} (100.0*dt) iso {Pstart} {Pstart} (100.0*dt) # Equilibrate for 10k steps at 0.1 K 1bar
run 10000

unfix equilibration
fix pressurization all npt temp {epsilon} {epsilon} (100.0*dt) iso {Pstart} {Pend} (100.0*dt) # pressurize for 1.2mil steps from 1bar to 500kbar

run 12000000

write_restart B1B2.equil

Here is the log thermo_style output of the 4 last timesteps:

Step TotEng Enthalpy Temp Press Volume
10001 -47516.781 -47523.273 0.29009039 -39.072283 11392.761
10002 -47514.678 -47077.577 0.29247004 2656.6289 11281.704
10003 -42338.089 -54381.349 0.19064583 -38484.138 21457.861
10004 -9.4308335e+21 -3.4579103e+22 0 -3.9915599e+30 0.00043200578

My potentials gave accurate (within 2%) densities at 300 K and 1 bar but this is of course only one test of validity and due to the repulsive part being an exponential, if particles have more than 44 eV in energy they (perhaps) can overcome the repulsion and fall into the singularity due to the (negative) VdW-terms?

I read earlier in the mail list that to approximate absolute zero you need to set temperature to a low non-zero value, maybe my temperature is too low?

Hopefully I didn’t leave out any crucial information, if so I apologize in advance. I’m still very new to this mail list and Lammps in general!

Best regards,

William Green

Doing simulations at very high pressure can be tricky. You have to be very careful with the choice of potential functional forms and parameters as well as simulation methods and settings.

Some suggestions:

  • use a small timestep and/or use fix dt/reset. you are in the steep repulsive part of the potential and thus forces can change rapidly with small displacements and thus cause instability of the time integration algorithm. fix nve/limit can also help avoiding instabilities while the system is still relaxing
  • use fix press/berendsen for pressure control and fix langevin (or csld, csvr, temp/berendsen, gld, gle etc.) instead of fix npt. nose-hoover thermostats and barostats are crucial for correct sampling of statistical mechanical ensembles in equilibrium, but not good at consistent transitions and to apply a suitable damping to avoid feedback oscillations which may result in crashes. if you just want to achieve a specific density, it may be simpler to simply use fix deform.
  • carefully check the shape of your potential at typical distances as you decrease the volume.
  • be mindful of the speed with which you transform your system. if you compress too fast, you may trigger instabilities quickly.

axel.

Hi Axel,

Thanks for all the useful information!

I added berendsen thermostats and barostats and a nve integrator like this (these are the changes from my previous input-file):

variable dt equal 1
variable epsilon equal 0.1
variable Pstart equal 1.0

timestep ${dt}

velocity all create ${epsilon} 12345 mom yes rot no dist gaussian # Take initial velocity from gaussian distribution at T = 0.1K

thermo 1

thermo_style custom step etotal lx ly lz temp press vol

fix equilibrationP all press/berendsen iso {Pstart} {Pstart} (1000.0*dt) # Equilibrate for 10k steps at T = 0.1 K p = 1bar fix equilibrationT all temp/berendsen {epsilon} {epsilon} (100.0*dt)
fix integrator all nve

run 10000

I tried to isolate the problem but having the system keep 0.1 K temperature and 1.0 bar pressure so that your last point about the speed of which the system is transformed was eliminated.

Unfortunately it seems that no matter what I tried (I tried a fix langevin instead of temp/berendsen as well) the simulation still hangs (without error) at time step 4. I tried different timesteps as well (0.1, 0.01 fs) but neither made the integration work.

the log for the 4 timesteps was as follows:

Step TotEng Lx Ly Lz Temp Press Volume
1 -47531.291 22.56 22.56 22.56 0.1 -2347.1748 11481.993
2 -47531.291 20.63451 20.63451 20.63451 0.099919567 -3067.4653 8785.824
3 -45662.811 49.494892 49.494892 49.494892 0.099512195 9275.4674 121249.83
4 -24017.835 40.456109 40.456109 40.456109 0.099189723 -8309.8744 66214.383

I must have used the wrong syntax or misunderstood the documentation because even setting temperature to 100 K will not let the integration proceed beyond step 4. Just to make sure I wasn’t going crazy I again tested a fix npt and it proceeds as long as temperature is at or above 2 K. I thought that fix nve could be used alongside the thermostat + barostat to integrate the atom trajectories/positions?

I also plotted the interpolated tables of the potential and as long as the separation is above 1 Å the force is repulsive. Since I’m simulating a NaCl structure (where if I’m not mistaken the nearest neighbor distance is a/2) this should imply that the lattice constant cannot be lower than 2 Å or the volume of the 4x4x4 cell system cannot be lower than 512 Å^3. As you can see the volume was far greater than this at all time steps. Of course, just because the atoms aren’t collapsing on each other doesn’t mean that the distance won’t be a problem, but as you said the forces might rapidly change with a small change in separation, but when plotting the force it looks very smooth (almost linearly horizontal) at a > 3 Å.

​Any further feedback would be greatly appreciated.

Best regards,

William Green

William,

I think a fundamental problem is your very bad initial geometry. because of moving the origin by 0.5 0.0 0.0 and not 0.5 0.5 0.5
you are placing the Cl anions in the same sheet as the Na cations and squeeze them in between the Na atoms.
that will be meta stable and thus minimization does not really help to resolve this, but once you start dynamics that breaks the symmetry and thus cancellation of the forces, you will get huge displacements and the symptoms you see.

check out the following minimal test input instead:

dimension 3
units real
boundary p p p
atom_style charge

lattice fcc 5.646 origin 0 0 0
region box block 0 4 0 4 0 4
create_box 2 box
create_atoms 1 box
lattice fcc 5.646 origin 0.5 0.5 0.5
create_atoms 2 box

mass 1 28.990
mass 2 35.453
set type 1 charge +1.0
set type 2 charge -1.0

pair_style born/coul/long 8.0
pair_coeff 1 1 9768.28 0.317 0 24.187 -11.52
pair_coeff 1 2 28937.74 0.317 0 161.21 -161.21
pair_coeff 2 2 80367.73 0.317 0 1669.6 -3353.63

kspace_style pppm 1.0e-5

timestep 1.0
minimize 1.0e-10 1.0e-10 100000 100000

reset_timestep 0

velocity all create 10.0 12345 mom yes dist gaussian
fix cooling all npt temp 10.0 0.01 100.0 iso 1.0 1000.0 1000.0 # Cooling at 1bar from 10 K to 0.1 K

thermo 100
thermo_style custom step etotal epair elong temp press vol

run 10000

axel.