[lammps-users] NPT pressure not stable

Hello Dr. Kohlmeyer,

I have been working with the Li-Al-O system. My simulation consists of minimization followed by NPT runs. The temperature in NPT (T_start = T_stop = 300) is stable at around 300 towards the end of the simulation but the pressure (Pstart = Pstop = 0) is not reaching zero. It is either a large positive or large negative (-7000 to +6000) value even during a long period of time. Am I missing something in the code?

units metal
dimension 3
boundary p p p
atom_style charge
read_data final2.data

variable mLi equal 6.941
variable mAl equal 26.9815
variable mO equal 15.9994
mass 1 {mLi} mass 2 {mAl}
mass 3 ${mO}

group Li type 1
group Al type 2
group O type 3

set type 1 charge 0.7
set type 2 charge 1.5
set type 3 charge -1.1

pair_style hybrid/overlay table linear 5000 coul/long 10.0
pair_coeff 1 1 table LiAlO2_TOT2009_WithoutCoulomb.table LiLi
pair_coeff 1 2 table LiAlO2_TOT2009_WithoutCoulomb.table LiAl
pair_coeff 1 3 table LiAlO2_TOT2009_WithoutCoulomb.table LiO
pair_coeff 2 2 table LiAlO2_TOT2009_WithoutCoulomb.table AlAl
pair_coeff 2 3 table LiAlO2_TOT2009_WithoutCoulomb.table AlO
pair_coeff 3 3 table LiAlO2_TOT2009_WithoutCoulomb.table OO
pair_coeff * * coul/long
kspace_style pppm 1.0e-4

timestep 0.0001
compute pe all pe/atom

neighbor 1.0 bin
neigh_modify every 5 delay 0 check no

thermo 20
thermo_style custom step dt time temp pe lx ly lz pxx pyy pzz
thermo_modify norm no

min_style cg
minimize 0 0 1000 10000
write_dump all custom dump.end.min id type x y z modify sort id

fix 1 all npt temp 300 300 1 iso 0 0 1 drag 1
thermo 100
thermo_style custom step temp press ke pe etotal vol density
run 1000000
write_dump all custom dump.end id type x y z modify sort id

Sincerely,
Ankit

1 Like

This topic has been discussed on this mailing list many, many times. Please study the mailing list archives.

That is the wrong method. There can be a large hysteresis, even for macroscopic experiments. There are discussions on this topic and using the coexistence method also in the mailing list archives.

Thanks Dr. Kohlmeyer.
I read a few of the mailing lists regarding this problem and learnt that one can check if the average pressure in the last few runs remains zero. I was able to get a zero average pressure taking the average of the last 1000 runs with 0.001 ps timestep.
But currently I am facing another problem while trying to calculate melting temperature via cooling and heating. When I ramp down the temperature from 2500 to 300 K, my system temperature is overshooting to an absurd value of 20000-25000. Is there anything I can rectify in this section?

units metal
dimension 3
boundary p p p
atom_style charge
#atom_modify map hash

read_data final2.data

variable myT equal 573
variable mLi equal 6.941
variable mAl equal 26.9815
variable mO equal 15.9994
mass 1 {mLi} mass 2 {mAl}
mass 3 ${mO}

group Li type 1
group Al type 2
group O type 3

set type 1 charge 1.75
set type 2 charge 5.25
set type 3 charge -3.5

#Li-Li
pair_style hybrid/overlay buck/coul/long 13 5 zbl 0.5 2
pair_coeff 1 1 buck/coul/long 0 0.285710 0
pair_coeff 1 1 zbl 3 3

Yes I read about hysteresis in this method. But I just wanted to get an approximate range of the melting point.
The problem I am facing is that the temperature is overshooting to 25K or 30K while ramping down from 2500 to 300 Kelvin. Is it because of the potential? The density that I obtained using this potential was very close to the experimental value.

If your model is even somewhat representative of real physics (and that’s not guaranteed), then the process of freezing your system should release a large amount of latent heat (enthalpy of fusion).

How much heat is released, and over what time period? Should you expect your thermostat to be able to remove the heat that quickly? What happens to your system if it can’t?

Do you know of a prior paper which uses a similar method to determine the melting point of a similar system?

Cheers,
Shern

Add-on question: is the force field parameterized to correctly model the liquid phase?