Fix GCMC loosing atoms

Carlos,
I have attached a complete input script and data file for a system that
fails with a lost atoms error. As far as testing with other codes, I have
used an old in-house MC code without issues for this system at higher
concentrations, where the lammps GCMC fix frequently runs into lost atom
problems. Most of the low concentration systems complete without issue.

Thanks,
Chris

variable t equal 300.00
variable temp string "300"
variable seed equal 31798.000000
variable phi string "q"
variable numH equal 200
variable mu equal -2.4000e+00
variable bulkH equal 3.0000e-05
variable e equal -0.01
units metal
boundary p f p
atom_style atomic
### set up simulation domain
read_data phi_\{phi\}\.data change\_box all x scale 1\.004347451179 y scale 1\.004347451179 z scale 1\.004347451179 variable rxlo equal xlo\+0\.0001 variable rxhi equal xhi\-0\.0001 variable rzlo equal zlo\+0\.0001 variable rzhi equal zhi\-0\.0001 region active block {rxlo} \{rxhi\} 47\.0000 107\.0000 {rzlo}
${rzhi} side in
group hydrogen type 4
group nickel type 1 2 3
# EAM potential -- Ni
pair_style eam/alloy
pair_coeff * * /Users/cjobrie/lammps-pot/NiH_jea.eam.alloy Ni Ni Ni H
mass * 58.710
mass 4 1.0080
#re-neighbor
neighbor 1.5 bin
neigh_modify every 5 delay 0
#THERMALLY EQUILIBRATE
thermo 0
log logqequ.lammps
velocity all create t {seed} mom yes rot yes dist gaussian
fix 7 all nvt temp $t t 10\.0 timestep 0\.001 run 100000 reset\_timestep 0 \#GCMC SETTINGS \# fID gID gcmc Nmdsteps Xchange Moves attype seed T mu maxtrans options fix gcmc hydrogen gcmc 100 {numH} 0 4 \{seed\} {t}.0 \{mu\} 0\.01 region active full\_energy compute\_modify thermo\_temp dynamic yes variable nh equal count\(hydrogen\) fix avenh all ave/time 50 100 5000 v\_nh ave one thermo 100000 thermo\_style custom step temp pe f\_avenh thermo\_modify format float %24\.16g restart 100000 restarteq\.{phi}.*
timestep 0.001
#re-neighbor
neighbor 1.5 bin
neigh_modify every 1 delay 0
run 3000000
#GCMC AVERAGING RUN
unfix avenh
variable cex equal "(count(hydrogen) -
count(nickel,active)*v_bulkH) / (lx*lz)"

unfix avenh
variable cex equal "(count(hydrogen) -
count(nickel,active)*v_bulkH) / (lx*lz)"
variable nh equal count(hydrogen)
variable tau equal (pxx+pzz)/2
compute PE all pe/atom
compute mype all pe
fix avetau all ave/time 2 250 1000 v_tau start 1000001 ave
running
fix avetemp all ave/time 2 250 1000 c_thermo_temp start
1000001 ave one
fix avecex all ave/time 2 250 1000 v_cex start 1000001 ave
running
fix avepe all ave/time 2 250 1000 c_mype start 1000001 ave
running
log logqave.lammps
dump 4 all custom 1000 dump.phi_\{phi\}\.\* id type x y z c\_PE dump\_modify 4 sort id pad 7 restart 20000 restart\_{phi}
thermo 1000
thermo_style custom step f_avetemp pxx pzz f_avetau f_avepe v_nh v_cex
f_avecex
thermo_modify format float %24.16g
timestep 0.001
run 200000

phi_q.data (276 KB)

Chris,

First, please confirm that this problem exists with the latest changes
to fix gcmc.

"7 Jul 2015
Aidan Thompson (Sandia) fixed some issues with the fix gcmc command.
Specifically, three changes have been made that affect how molecules
are inserted and rotated. The first change corrects a problem with the
periodic image flags that affected rotation moves for molecules that
had previous been inserted at a periodic boundary. The second change
modifies the distribution of rotation angles used for molecule
insertions and rotations. The third change eliminates a problem
inserting charged atoms and molecules when full_energy option is not
used."

None of these advertized changes seem relevant to your problem, but
nonetheless it is important to check.

If the problem still exists, please try to construct a small fast
instance. Think of the movie "Gone in Sixty Seconds." I will take a
look at it.

Aidan

Aidan,
I am currently using the 7/22/15 version of lammps. Unfortunately, I don’t
really have an example that instantly yields the error. It is
unpredictable and can occur at any time step. The included input file will
sometimes show the error. Also, I have noted that the presence of the
region that GCMC is confined to does not cause the error, it will
occasionally loose atoms without the region.

variable t equal 300.00
variable temp string "300"
variable seed equal 87248.000000
variable theta string "a"
variable numH equal 200
variable mu equal -2.4000e+00
variable bulkH equal 3.0000e-05
variable e equal -0.01
units metal
boundary p f p
atom_style atomic
### set up simulation domain
read_data theta_\{theta\}\.data variable rxlo equal xlo\+0\.0001 variable rxhi equal xhi\-0\.0001 variable rzlo equal zlo\+0\.0001 variable rzhi equal zhi\-0\.0001 region active block {rxlo} \{rxhi\} 30\.9934 90\.9934 {rzlo}
${rzhi} side in
group hydrogen type 4
group nickel type 1 2 3
# EAM potential -- Ni
pair_style eam/alloy
pair_coeff * * /Users/cjobrie/lammps-pot/NiH_jea.eam.alloy Ni Ni Ni H
mass * 58.710
mass 4 1.0080
#THERMALLY EQUILIBRATE

#re-neighbor
neighbor 1.5 bin
neigh_modify every 5 delay 0
thermo 0
log logaequ.lammps
velocity all create t {seed} mom yes rot yes dist gaussian
fix 7 all nvt temp $t t 10\.0 timestep 0\.001 run 1000 reset\_timestep 0 \#GCMC SETTINGS \# fID gID gcmc Nmdsteps Xchange Moves attype seed T mu maxtrans options fix gcmc hydrogen gcmc 100 {numH} 0 4 \{seed\} {t}.0 ${mu}
0.01 region active full_energy
compute_modify thermo_temp dynamic yes
variable nh equal count(hydrogen)
fix avenh all ave/time 10 100 1000 v_nh ave one
thermo 1000
thermo_style custom step temp pe f_avenh
thermo_modify format float %24.16g
timestep 0.001
#re-neighbor
neighbor 1.5 bin
neigh_modify every 1 delay 0
run 3000000

theta_a.data (66.8 KB)

This is too vague and too big.

First, tell me a few more things:

1. Is the error reproducible i.e. does it happen every time
2. Roughly how long does it take to generate the error

Second, at least make an effort to eliminate parts of the computation
that are not needed to generate the error. The less lines of script I
have to look at, and the fewer ticks of the clock I have to wait, the
more likely I will be to isolate the problem. In fact:

P( time_find_bug < time_run_out of_patience ) ~
exp(-k1*lines_of_script) exp(-k2*time_to_error)
exp(-k3*number_of_atoms)

The value of 1/k1, 1/k2, and 1/k3 are roughly 20, 60secs, and 100, respectively.

Aidan

P( time_find_bug < time_run_out of_patience ) ~
exp(-k1lines_of_script) exp(-k2time_to_error)
exp(-k3*number_of_atoms)
The value of 1/k1, 1/k2, and 1/k3 are roughly 20, 60secs, and 100, respectively.

Nice. Finally, a quantitative formula for patience thresholds.

I suggest a Nature paper broadening this with coefficients

for Internet memes, political scandals, interpersonal relationships, etc.

Steve