Fix GCMC for pure LJ bulk fluid

Hi Lammps users and developers,

I am trying simulate pure LJ fluid with fix GCMC. I had a couple of questions:

  1. In lammps , lj units mu* = mu/epsilon ?

  2. According to Johnson et. al., Mol. Phys, 78(3),1993 , rho* = 0.001, T* = 0.75 corresponds to a mu* = -5.19
    I am using it to validate my gcmc simulation in lammps in 10 X 10 X 10 (sigma)^3 box with a 50% creation/destruction moves and rest being displacement moves.

The densities I am getting are around O(0.6). I request you to please look at my setup in the script below.

echo screen
units lj
atom_style atomic
boundary p p p

region box block 0 10 0 10 0 10

create_box 1 box

mass 1 1.0

pair_style lj/cut 2.5

pair_modify shift yes

pair_coeff 1 1 1.0 1.0 2.5

neighbor 0.3 bin

neigh_modify every 1 delay 0 check no

variable vol equal 101010

variable dens equal atoms/v_vol

variable chempot equal -5.19

compute_modify thermo_temp dynamic yes

fix fixgcmc all gcmc 1 50 100 1 123213 0.75 ${chempot} 1

thermo_style custom step temp press pe ke atoms v_dens v_chempot

thermo 100

dump m0 all movie 1 test.mpg type type size 1280 720 zoom 1.5

run 1000000

Thanks for help in advance.

Thanks & Regards

Ashutosh Rathi

Graduate StudentDepartment of Chemical Engineering
University of Massachusetts, Amherst

“Success starts with being different and ends with being superb”

In LAMMPS, mu and epsilon are whatever you specify them to be. The choice of units via the units command primarily affects the values used for physical constants, specifically Boltzmann’s constant. If you specify lj units, then k=1, which is probably what you want here. I can not speak to how Johnson et al defined, but it looks like your script is following the commonly used conventions for LJ. Things to worry about are:

  1. Johnson is probably generating a vapor phase, but you are initializing you system with the default lattice “lattice none 1.0”, which for lj units means a reduced density rho*=1.0, so you may be getting stuck in a metastable solid or liquid phase
  2. The LJ model is not uniquely defined by sigma and epsilon. You also need to worry about rcut, shift yes/no, tail yes/no.


Dear Aidan,

I am using 7-Dec-2015 version of lammps.

In regards to the initialization of system, I am a bit confused on this as to how “lattice none 1.0” would affect my system when I did not use it to create any atoms initially. I start gcmc with an empty system. So if I understand correctly, position of atoms will be in accordance with gcmc criteria and lattice should not have any effect on it.

Also, if you could share any sample/test scripts for atomic lj fluids using fix gcmc , it would be a great help.

Thanks & Regards
Ashutosh Rathi

If you are starting from an empty system, then yes, metastability should not be a problem. However, you still have to worry about point 2, the definition of LJ. Here is another suggestion. Reduce the value of the chemical potential sufficiently to generate a stable vapor. I also noticed this comment in the doc page:

If LJ units are used, note that a value of 0.18292026 is used by this fix as the reduced value for Planck’s constant. This value was derived from LJ parameters for argon, where h* = h/sqrt(sigma^2 * epsilon * mass), sigma = 3.429 angstroms, epsilon/k = 121.85 K, and mass = 39.948 amu.

The point is, there is no unique or even widely accepted reference point for the zero of chemical potential. Depending on the application area, different choices are made. Here is what the code actually does:

beta = 1.0/(force->boltzreservoir_temperature);
double lambda = sqrt(force->hplanck
zz = exp(beta