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