GCMC with region keyword stops after couple of steps

Dear users and developers,
I am using GCMC-NVT in LAMMPS to study surface adsorption.
The simulation runs without problem. But when I activate the region keyword with GCMC it gets stuck after couple of GCMC trials.

The surface (~ 52 x 52 nm^2) faces Z direction with a large vacuum extending 1.5 — 5 nm .

The adsorbate is roughly the size of a cyclohexane.

I like to limit the GCMC insertion in the vicinity of surface. Otherwise, after 10^6 steps, GCMC fills the vacuum with adsorbates.
So, I have to either do GCMC with region keyword or shrink the vacuum to ~1nm.

With the second option LAMMPS fails with the “Fix gcmc put atom outside box” error message as I use “boundary p p f” for the slab.

With the region keyword (with large vacuum region) it gets stuck without any error message. The documentation says with the region keyword LAMMPS tries 1000 insertions before moving on to the next step. I think this is causing the problem. I’m wondering if you know a way to overcome this issue.

my LAMMPS version is 12Dec18.
Following is my input.

Thanks for your time.
Sal

units real
boundary p p p
neighbor 2.0 bin
neigh_modify delay 1 every 1 check yes
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic
pair_style lj/charmm/coul/long 10 12
pair_modify mix arithmetic
kspace_style pppm 1e-5
kspace_modify order 6
special_bonds charmm

read_data res.data
molecule adsmol ads.txt
group ads type 1:12

region vcm block -27.93 27.93 -26.8535 26.8535 15.0 25.0 side in
velocity all create 300 12345678 dist uniform
variable tfac equal 6.0/3.0
variable mu index -10.0
variable disp index 1.0

compute mdtemp all temp
compute_modify mdtemp dynamic/dof yes
fix md all nvt temp 300 300 100
fix_modify md temp mdtemp
fix mygcmc ads gcmc 100 10 10 0 54341 300 {mu} {disp} mol adsmol &
tfac_insert ${tfac} group ads region vcm

variable adsN atom “type==5”
group adsN dynamic all var adsN
variable nads equal count(adsN)

thermo_style custom step temp press pe ke density atoms v_nads
thermo 100
timestep 1.0
run 1000

Hi everyone,
To provide more information about this issue, following is the last couple of steps in the output that has been hanging over 24 hours! I don’t understand why GCMC doesn’t move on to the next step when it can’t find a place to insert the molecule. This problem doesn’t happen without the region keyword even when the vacuum is completely filled with adsorbates after millions of steps.

Thanks for you time.
Sal

Step Temp Press PotEng KinEng Density Atoms v_nads.
.
.

6200 363.87641 1180.213 -450763.77 6154.284 1.5854875 6893 43
6300 361.33981 2195.2608 -450699.47 6111.3822 1.5854875 6893 43
6400 361.1288 944.52535 -450562.07 6107.8134 1.5854875 6893 43
6500 376.34145 1085.2092 -450705.71 6365.1066 1.5854875 6893 43

Hi Sal,

About your first option.
If your adsorbate has affinity for your surface, you should see adsorptions at chemical potentials below the chemical potential thats corresponds to the equilibrium condition between a bulk condensed (simulated) phase of your adsorbate and gaseous ( infinity reservoir at fix chemical potential phase ) . So my question… Do you know the equilibrium chemical potential of a bulk phase of your adsorbate, you should get this value in order to choose the chemical potential of your adsorption simulations in a more rational way, with that you might observed adsorption in your surface with filling the box. This is the same as knowing the vapor pressure of adsorbate at the simulation temperature.

For the second option. Maybe you could attached a simple example to reproduce the error message.

Hopes this help
Good luck

Hi everyone,
To provide more information about this issue, following is the last couple of steps in the output that has been hanging over 24 hours! I don’t understand why GCMC doesn’t move on to the next step when it can’t find a place to insert the molecule. This problem doesn’t happen without the region keyword even when the vacuum is completely filled with adsorbates after millions of steps.

Thanks for you time.
Sal

Step Temp Press PotEng KinEng Density Atoms v_nads.
.
.

6200 363.87641 1180.213 -450763.77 6154.284 1.5854875 6893 43
6300 361.33981 2195.2608 -450699.47 6111.3822 1.5854875 6893 43
6400 361.1288 944.52535 -450562.07 6107.8134 1.5854875 6893 43
6500 376.34145 1085.2092 -450705.71 6365.1066 1.5854875 6893 43

Thanks for your suggestion Matias.

Unfortunately, moderator doesn't allow me to send my files as their
size is 2mb and Gmail blocks compressed files.

I did calculate the excess chemical potential with FEP method on the
bulk adsorbate molecules. The self solvation free energy for my
adsorbate is ~40 kcal/mol.

Hi Sal,

About your first option.
If your adsorbate has affinity for your surface, you should see adsorptions
at chemical potentials below the chemical potential thats corresponds to
the equilibrium condition between a bulk condensed (simulated) phase of
your adsorbate and gaseous ( infinity reservoir at fix chemical potential
phase ) . So my question... Do you know the equilibrium chemical potential
of a bulk phase of your adsorbate, you should get this value in order to
choose the chemical potential of your adsorption simulations in a more
rational way, with that you might observed adsorption in your surface with
filling the box. This is the same as knowing the vapor pressure of
adsorbate at the simulation temperature.

For the second option. Maybe you could attached a simple example to
reproduce the error message.

Hopes this help
Good luck