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

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.

# 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

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

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