GCMC Simulation for NaCl Using Born-Mayer-Huggins Potential in LAMMPS

LAMMPS LAMMPS General Discussion
1

Subject: Help with GCMC Simulation for NaCl Using Born-Mayer-Huggins Potential in LAMMPS

Hello everyone,

I’m attempting to run a Grand Canonical Monte Carlo (GCMC) simulation in LAMMPS for NaCl using the Born-Mayer-Huggins potential. My goal is to include the following moves in the simulation:

  1. Insertion/deletion of atoms,
  2. Position swapping,
  3. Displacement,
  4. Mutation (identity exchange) moves.

However, my simulation repeatedly fails with the error:
ERROR: Lost atoms: original 61 current 59

I’ve checked the simulation box boundaries and potential parameters, but I’m unsure if there are other aspects specific to GCMC or the Born-Mayer-Huggins potential that could be causing this issue. Has anyone encountered similar problems with atom loss in GCMC simulations, especially with NaCl or other ionic compounds? I’d appreciate any guidance on setup tips or parameters that might help stabilize the system.

Thanks in advance for your help!

# GCMC simulation for NaCl system at 800 K with GCMC statistics output

units real
atom_style charge
boundary p p p  # Periodic boundaries

# Read the final configuration from previous simulation (assuming restart file or data file)

# Define simulation box dimensions
region box block 2.273551041 101.5614489589 2.273551041 101.5614489589 2.273551041 101.5614489589
create_box 2 box

# Define atom types and charges
create_atoms 1 random 50 12345 box  # 50 Na atoms
create_atoms 2 random 50 67890 box  # 50 Cl atoms

set type 1 charge +1.0  # Na charge
set type 2 charge -1.0  # Cl charge


# Define masses (make sure this matches with the system)
mass 1 22.989  # Na mass
mass 2 35.453  # Cl mass

# Pair potentials (Born-Mayer-Huggins or any appropriate potential for your system)
pair_style born/coul/long 12.0
pair_coeff 1 1 6.08106 0.317 2.34 24.1789 -11.5141  # Na+ - Na+
pair_coeff 2 2 3.64817 0.317 3.17 1669.6160 -3353.917  # Cl- - Cl-
pair_coeff 1 2 4.86346 0.317 2.755 161.2047 -200.06409  # Na+ - Cl-

# K-space style for long-range interactions
kspace_style pppm 1.0e-5

# Define the temperature for the GCMC run
variable rdm equal 544645 
variable pre equal 1
variable tem equal 800  # Temperature in Kelvin
variable dis index 0.5 # Displacement:
variable tfac equal 5.0/3.0
variable mu equal -91.7829828 #  Chemical Potential
variable tdamp equal 0.5
# Neighbor settings
neighbor 3.0 bin
neigh_modify every 5 delay 10 check yes one 10000

# Group definitions
group Na type 1 
group Cl type 2 
group NaCl type 1 2

# Fix for temperature control
fix nvt all nvt temp ${tem} ${tem} ${tdamp}

# Fix GCMC for Na and Cl
#fix gcmc_Na Na gcmc 100 50 0 0 ${rdm} ${tem} ${mu} ${dis} tfac_insert ${tfac}   
#fix gcmc_Cl Cl gcmc 100 50 0 0 ${rdm} ${tem} ${mu} ${dis} tfac_insert ${tfac}

#fix gcmc_NaCl all gcmc 100 0 50 0 ${rdm} ${tem} ${mu} ${dis} tfac_insert ${tfac}  full_energy
fix gcmc_NaCl all gcmc 150 1 50 1 ${rdm} ${tem} ${mu} ${dis} tfac_insert ${tfac}  full_energy

#GCMC Steps Between Trials: 150 
#Insertion Attempts: 50 
#Translation Attempts: 50 
#Deletion Attempts: 2

# Time step for the simulation
timestep 1.0  # 1 fs timestep

# Output thermodynamic information every 1000 steps
thermo 1000
#thermo_style custom step temp press density atoms f_gcmc_Na[1] f_gcmc_Na[2] f_gcmc_Na[3] f_gcmc_Na[4] f_gcmc_Na[5] f_gcmc_Na[6] f_gcmc_Na[7] f_gcmc_Na[8] f_gcmc_Cl[1] f_gcmc_Cl[2] f_gcmc_Cl[3] f_gcmc_Cl[4] f_gcmc_Cl[5] f_gcmc_Cl[6] f_gcmc_Cl[7] f_gcmc_Cl[8]

thermo_style custom step temp press density atoms f_gcmc_NaCl[1] f_gcmc_NaCl[2] f_gcmc_NaCl[3] f_gcmc_NaCl[4] f_gcmc_NaCl[5] f_gcmc_NaCl[6] f_gcmc_NaCl[7] f_gcmc_NaCl[8]

# Dump every configuration to file
dump gcmc_dump all custom 100 dump.gcmc.* id type x y z

# Run GCMC for 1,000,000 steps
run 1000000
2

Your create_atoms seems to be intended to create 50 sodium ions and 50 chloride ions, but your pasted log output indicates you have only 61 atoms in the system (and then 59), so something has went wrong very early on.

I would also encourage you to search the literature to see what other researchers have done for similar situations. (By contrast, the LAMMPS forum is largely run by a handful of academic volunteers who may or may not have seen this particular simulation type before.)

In particular I am not sure whether it even makes sense to insert or delete single ions, since such a Monte Carlo move changes the overall charge of the box and the resulting energy change is not physically defined (even though a mathematical result can be reached taking certain assumptions).

3

you right I wil like to delete and insert Na+ and Cl- simultaneously I was woundering how can I do that in lammps

4

Before discussing the technical details, what is the purpose of using GCMC and not just plain MD? When using LAMMPS you have to keep in mind that it was designed for MD and using any kind of MC step is rather inefficient and thus should be used rarely (say every 1000 MD steps or so). Otherwise, it would be advisable to look into using a real MC code which doesn’t suffer from the many limitations of MC steps in LAMMPS.

What is the reason for simulating such a tiny system? In particular in combination with MC steps. This makes almost all MC steps very disruptive and it is very difficult to dissipate and energy from moves that increase the potential energy. Also, for such a small system, you have serious finite size problem with atoms interacting with its own periodic copies.

As I have mentioned in a message yesterday, this type of potential has a problem when atoms come too close, since the non-bonded repulsion is not strong enough to overcome the Coulomb interaction for atoms of opposite charge if they come too close. With MC insertion or displacement move, that is likely to happen eventually and thus can crash or invalidate your simulation. Already for a plain MD run it is important to start from a sufficiently well prepared initial geometry without local high potential energy areas that would lead to fast moving atoms.

It is , unlikely that somebody here has exactly the expertise you need. All we can do is give some general advice, but that would be for the most part off-topic here and more suitable for the Science Talk category. As usual, the advice is to study the published literature and to find a suitably experienced collaborator.

5

Thank you

my objectif was to generate some synsthetic data but I think I will write a custom GCMC code

Thank you for your input

Best,

Conrard T.

6

This doesn’t really answer any of my questions but suit yourself since it very much looks like you are much better off not using LAMMPS for your purposes.

7

I interestd in using GCMC because I assume my system is dominated by frictional or viscous forces and also I want ot have a direct control of Chemical Potential.
I am also interested in calculating the viscosity.