[lammps-users] Buckingham potential: fluid EoS comparison with Monte Carlo results

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,

I tried to do rather a simple thing: to reproduce the equation of state of the
fluid of particles interacting via the Buckingham (exp-6) potential and to
compare the results with the data from the literature (Fried&Howard
J.Chem.Phys. Vol.109 p.7338, http://dx.doi.org/10.1063/1.476520 ).
But there is a large discrepancy in Lammps and MC results!

For example, if we consider the EoS in the dimensionless parameters Z=PV/NkT
and U=Uex/NkT

for:
alpha=15.5 (the conventional Buckingham pot. paramenter)
T*=100 (temperature in reduced units)
rho*=5.8025 (density in reduced units)

we have:
from Lammps equilibrium NVT trajectory Z=1.4 U=0.1
the MC accurate results Z=43.2 U=12.7

I checked the effects of the cut-off radius, the averaging time, the system
size. None of them has the strength enough to converge the results.

Below I attach the sample input script.

The problem might be not directly connected with Lammps, but I feel it can.

I think I should check this via the simple MD code,
but before doing this I hope that maybe the community can point out
if I miss some important aspects of the Lammps application in the case.

Thank you.

Vladimir Stegailov

# LAMMPS 12Apr2006 tested
# Buckingham

variable x index 10
variable y index 10
variable z index 10

variable f equal 15.5 #alpha
variable t equal 100 #T* (epsilon = 1)
variable d equal 5.8025 #rho* (r_m = 1)

variable l equal pow(div(1,$d),div(1,3)) #sc lattice constant
variable o equal 3.0 #cut-off distance
variable a equal mult(div(6,sub($f,6)),exp($f)) #A(LAMMPS exp-6 param.)
variable r equal div(1,$f) #rho (LAMMPS exp-6 param.)
variable c equal div($f,sub($f,6)) #C (LAMMPS exp-6 param.)

#reduced pressure Z=PV/NkT
variable p equal div(mult(press,vol),mult(atoms,temp))
#reduced excess energy Uex=dU/NkT
variable u equal div(pe,mult(atoms,temp))

print lattice const=$l A_buck=$a rho_buck=$r C_buck=$c

units lj
atom_style atomic

lattice sc $l
region box block 0 $x 0 $y 0 $z units lattice
create_box 1 box
create_atoms 1
mass 1 1.0

velocity all create $t 35058 dist gaussian mom yes rot yes

pair_style buck $o
pair_coeff 1 1 $a $r $c

timestep 0.001

thermo 50
thermo_modify norm no
fix 0 all print 50 "Z= $p Uex= $u"

#preparation (in order to melt the initial lattice)
fix 1 all nve
run 1000
unfix 1

#production run
velocity all create $t 35058 dist gaussian mom yes rot yes
fix 1 all nvt $t $t 0.1
run 10000

2

Vladimir,

You're not simulating at the correct density. Try the
script below. The answers I'm getting when I use the
correct density give much better agreement with the
published MC results. The essential change is to use
"lattice sc $d" instead of "lattice sc $l". See also:

http://www.cs.sandia.gov/~sjplimp/lammps/doc/lattice.html

Paul

# LAMMPS 12Apr2006 tested
# Buckingham

variable x index 10
variable y index 10
variable z index 10

variable f equal 15.5
#alpha
variable t equal 100 #T*
(epsilon = 1)
variable d equal 5.8025
#rho* (r_m = 1)

variable l equal pow(div(1000,$d),div(1,3))
#sc lattice constant
variable o equal 3.0
#cut-off distance
variable a equal mult(div(6,sub($f,6)),exp($f))
#A(LAMMPS exp-6 param.)
variable r equal div(1,$f) #rho
(LAMMPS exp-6 param.)
variable c equal div($f,sub($f,6)) #C
(LAMMPS exp-6 param.)

#reduced pressure Z=PV/NkT
variable p equal
div(mult(press,vol),mult(atoms,temp))
#reduced excess energy Uex=dU/NkT
variable u equal div(pe,mult(atoms,temp))

print lattice const=$l A_buck=$a rho_buck=$r
C_buck=$c

units lj
atom_style atomic

lattice sc $d
region box block 0 $x 0 $y 0 $z units lattice
create_box 1 box
create_atoms 1
mass 1 1.0

velocity all create $t 35058 dist gaussian
mom yes rot yes

pair_style buck $o
pair_coeff 1 1 $a $r $c

timestep 0.001

variable h equal div(atoms,vol)

thermo 50
thermo_modify norm no
fix 0 all print 50 "Z= $p Uex= $u
rho=$h"

#preparation (in order to melt the initial lattice)
fix 1 all nve
run 1000
unfix 1

#production run
velocity all create $t 35058 dist gaussian
mom yes rot yes
fix 1 all nvt $t $t 0.1
run 10000