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