LAMMPS RDF Size Dependence

LAMMPS LAMMPS General Discussion
1

Lammps community,
I am computing the RDF for a periodic 3D fcc lattice using Lennard-Jones. I am seeing a size-dependent effect on my RDF results, and I am wondering why. For a periodic box, the formulation of g(r) should not depend on the number of atoms in the system, yet when I increase nx, ny, and nz in my script (seen below), each g(r) value decreases significantly. In testing this problem, the values eventually converged fully at ~1000000 atoms, but I don’t see why they should change at all with the number of atoms.

# Simulation settings
package gpu 1
boundary p p p
atom_style atomic
units metal

# General model information
variable nx equal 1 # number of cells in x direction
variable ny equal 1 # number of cells in y direction
variable nz equal 1 # number of cells in z direction

#Create Atoms
lattice fcc 1.0
region box block 0 ${nx} 0 ${ny} 0 ${nz} units lattice
create_box 1 box
create_atoms 1 box
mass 1 1.0

# Potential
variable rcut equal 2.5
variable padding equal 2*${rcut}
variable sigAA equal 1.0
variable epsAA equal 1.0
pair_style lj/cut/gpu ${rcut}
pair_coeff 1 1 ${epsAA} ${sigAA}
pair_modify shift yes

# Compute radial distribution (RDF) function
comm_modify cutoff ${padding}
compute myRDF all rdf 100 1 1 cutoff ${rcut}
fix fix_rdf all ave/time 1 1 1 c_myRDF[*] file out.rdf mode vector
run 0

I also computed the RDF separately by importing the dump file into OVITO, and did not see a size-dependent effect here – increasing the number of atoms did not have an effect on the g(r) values computed in OVITO (computed with the same number of bins and cutoff).

Any insight into this issue within lammps would be appreciated.

Thank you,
Chloe

2

Please always report which version of LAMMPS you are using and what platform you are running on. What values of nx/ny/nz have you tested and are you comparing?

There should be a small difference in peak heights due to the system size because of a finite size effect: the “pair” of each atom with itself is not included in the g(r) calculation and thus the normalization must take that into account. Since the number of pairs with itself scales O(N), but the total number of pairs O(N**2) the effect gets smaller as the system gets larger. If you don’t do this correction for the finite size of the simulation cell, the g(r) value for r->oo will not become exactly 1.0 but something like N*(N-1)/N**2.

Since this effect is small for systems with many atoms not everybody notices the difference and implements the correction.

For more information, you can have a look into some text books on that matter. I don’t recall whether I found this in “Hansen/McDonald, Therory of Simple Liquids” or “Gray/Gubbins, Theory of Molecular Fluids Vol.1”. It has been a long time (~25 years) since I first worked on something where this mattered: Long-range Structures in Bulk Water. A Molecular Dynamics Study