Simulating liquid Argon in lbulk

Hi ,

I am trying to simulate liquid argon in bulk using NVT. I’m following the old paper of Rahman et. al (Phys.Rev, Vol 136,1964) with a temperature 94.4K and density as 1.374 g/cm^3. However after around 5000ps the system is showing the pressure of 150-200 bar, whereas according to NIST the pressure at this state point should be 73 bar. Why the system is not equilibrating to the accurate pressure? I kept the density same as 1.374 g/cm^3

Here is my script:

#-----------initialize simulation--------------
units metal
dimension 3
boundary p p p
atom_style sphere
atom_modify map array

#-----------create atoms--------------------

region simbox block -17.3893 17.3893 -17.3893 17.3893 -17.3893 17.3893 side in units box
create_box 1 simbox
create_atoms 1 random 864 125698 NULL

set type 1 mass 39.95

set type 1 diameter 3.4

#--------------grouping---------------------------------

group Ar type 1

---------- Define Interatomic Potential ---------------------

pair_style lj/cut 10.5
pair_coeff * * 0.0103 3.4 10.5
pair_modify shift yes

#-------------Neighbour list-------------------------------------

neighbor 2 bin
neigh_modify delay 10 check yes

----------------minimization of solvent alone---------------------

fix 1 Ar setforce 0.0 0.0 0.0

thermo 100
thermo_style custom step temp pe ke etotal press vol enthalpy

min_style cg
minimize 1e-25 1e-25 100000 1000000
unfix 1

#-------------- minimization of solution-----------------------

thermo 100
thermo_style custom step temp pe ke etotal press vol enthalpy

min_style cg
minimize 1e-25 1e-25 100000 1000000

print “MINIMIZATION DONE”

#write_restart solvated_after.minimize

#----------------equilibration--------------------------------

timestep 0.0001
velocity all create 94.4 123475 dist gaussian
fix 3 all nvt temp 94.4 94.4 0.01
fix 4 Ar momentum 100 linear 1 1 1
thermo 1000
thermo_style custom step time temp pe ke etotal press vol density
dump 6 all xyz 200000 dump_Ar_liquid_500.xyz

log lognew.txt

#----------------- Run-------------------

run_style verlet
run 100000000

Thanks in advance,

Kanka Ghosh

Dept. of Physics

IIT Madras

India

Maybe you need tail corrections on the LJ potential?

E.g. pair_modify tail yes

Steve

Thanks Steve. I’ll do that check. But as I was using pair_modify shift yes to shift the LJ potential to zero at cutoff, I couldn’t use that. Can you tell any other possibility for the discrepancy? As the calculation for liquid Argon seems trivial.

Thanks,

Kanka

It’s actually not trivial. You have to know precisely how the

earlier calculation was done (e.g. tail corrections, cutoffs, shifting, etc). And you

have to figure out the corresponding settings in LAMMPS,

if you want any hope of getting the same answers.

Steve

Dear Kanka

I get the desired numbers (with acceptable accuracy) with an input file which is available at the end of this email. I set T=94.4 and P=73 bar by using NPT ensemble to get the corresponding density.

a few suggestions:

1- use a higher cutoff

2- check the epsilon and sigma (and also the unit conversions)

3- use a NPT first

4- check if your system is solid or liquid!! (by a visualization software like VMD).

Input file:

Thanks Mostafa and Steve!!! Thanks for your help.

Hello Mostafa,

Thanks for your concern.Looking at your output file it’s clear that your system is equilibrating at the desired density. But you simulated NPT system with target P=73 bar (or atm in real unit), whereas your output pressure values are not very close to 73. Have you noticed that??

Thanks,

Kanka

Hello Mostafa,
   Thanks for your concern.Looking at your output file it's clear that
your system is equilibrating at the desired density. But you simulated NPT
system with target P=73 bar (or atm in real unit), whereas your output
pressure values are not very close to 73. Have you noticed that??

​kanka,

please have a look at an MD text book again.​

what LAMMPS output in the standard thermodynamic output is the
*​instantaneos* pressure. in condensed systems, this will *always*
fluctuate a lot, and thus the variations for low pressure targets will be
large since there is not much compressibility. the only pressure value that
you can compare to, is an averaged pressure over a sufficiently long chunk
of a trajectory taken after it is properly equilibrated.

​axel.​

Dear Kanka

Yes I noticed that and as Axel said it is ok.