Using fix gcmc

Dear Lammps-users,

Recently I am trying to use fix gcmc to simulate Lithiation in amorphous Silicon bulk. However what I find is that the temperature I get from the thermal output is extremely high initially, much higher than what I assigned in fix gcmc command (300 K), and then falls close to 0 K. I have attached a log file of the same. I am using lammps version from August 10, 2015.
The changes I tried are:

i. Used fix nve/limit to impose a limit on the maximum distance an atom can move

ii. Tried both by using and not using fix nvt with GCMC

When I dont use fix NVT with gcmc, of course I achieve the desired temperature but I lose the included atoms immediately when I try NPT equilibration after this.

iii. Changed the valus of chemical potential

iv. changed the fix gcmc settings to change the speed at which atoms are inserted, i.e. run it more or less frequently

v. Tried with the latest Lammps version.

I know fix gcmc was designed to study distribution and diffusion of gases in porous media, but is it not possible to simulate the effect of Lithiation in bullk amorphous silicon. So should I be using fix deposit for the same. Any suggestion would be great.

I am using MEAM potential.

Here is the input script for nvt+gcmc

units metal
atom_style atomic
boundary p p p

atom_style atomic
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

read_data amor.si_box

pair_style meam
pair_coeff * * library.meam Si Li LiSi.meam Si Li

compute_modify thermo_temp dynamic yes

group g_Si type 1
group g_Li type 2
region lireg block 0 30 0 30 12 30
group lithiation region lireg

variable chempot equal -2.14

compute tt all temp
compute_modify tt dynamic yes extra 0

dump 2 all atom 1000 inserttestmeamnpt

velocity all create 300 481516

timestep 0.001

variable nM loop 150
label loop
variable b equal v_chempot+0.14*v_nM

if “${nM} > 1” then “unfix f_npt”

fix flith all gcmc 500 20 20 2 29494 300.0 $b 0.01

thermo_style custom step atoms f_flith[1] f_flith[2] f_flith[3] f_flith[4] f_flith[5] &
f_flith[6] temp pe ke etotal press vol v_chempot

thermo 500
thermo_modify temp tt lost/bond ignore lost warn

#fix 7 all nve/limit 0.1
fix f_nvt all nvt temp 300 300 0.1
run 5000
next nM
unfix flith
unfix f_nvt

thermo_style custom step atoms temp pe ke etotal press vol
thermo 1000
thermo_modify temp tt lost/bond ignore lost warn

variable myP equal press
variable myT equal temp
fix f_npt all npt temp {myT} 300.0 0.1 iso {myP} 0.0 1

run 100000
jump gcmcex4 loop

log.lammps (2.26 KB)

Here are some observations:

  1. The log file and input script that you provide do not seem to be from the same simulation
  2. The input script contains some strange logic, such as putting the “next” command in the middle of the loop body.
  3. The first reported temperature after initialization is 40,000K.
  4. the simulation is already garbage after 500 steps, but probably much sooner than that.

I suspect your troubles are caused by insertion of an Li atom at a very small separation from a Si atom. This position may correspond to an unphysical region of the potential with a negative divergence as R goes to zero. For MD, this configuration will never appear, because it is separated from physical configurations by a large energy barrier. You can understand this better by performing a simulation with just 1 Si atom at a low chemical potential and checking the energy after every attempted insertion.

The solution would be to fill in the unphysical region with a repulsive well. This could be done by modifying the pair_meam.cpp code, or by using pair_style hybrid/overlay meam zbl

Aidan

Thanks a lot Aiden…

Also I had a query… According to the errors manual:

*Cannot do GCMC on atoms in atom_modify first group:* This is a restriction due to the way atoms are organized in a list to enable the atom_modify first command.
Again, in case of using read_restart, i understand we have to turn off atom sorting using the command atom_modify 0 0.
So, is there a problem with using read_restart and fix gcmc? I tried using but though the fix output shows that atoms are being inserted there is actually no physical presence of the inserted atoms! I am my script here and I attach the logfile .. (i double checked.. its the correct log file!)

read_restart restartnew2.MD_260000

change_box all z final 0 100 boundary p p f
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
atom_modify sort 0 0

variable chempot equal -2000
variable density equal atoms/vol
variable specdensity equal v_density1.39156E-251e27
variable chempotkcal equal -2
variable x equal xlo

variable tmpx equal “xlo”
variable tmpy equal “xhi”
variable tmpz equal “lz”
variable Lx equal {tmpx} variable Ly equal {tmpy}
variable Lz equal {tmpz} print "Initial Lengths, Lx = {Lx}, Ly = {Ly}, Lz = {Lz}"

region r_surf prism 0 33 0 33 17 34 0 0 0

region r_fix prism INF INF INF INF 0 4 0 0 0

compute tt all temp
compute_modify thermo_temp dynamic yes

group g_Si type 1

#timestep 0.001

group g_depo type 2

group g_surf region r_surf

group g_fix region r_fix

#timestep 0.001

group mobile subtract all g_fix

dump mydump all atom 100 tmpgcmc.deposit

compute add g_depo temp
compute blockt mobile temp
compute subt g_fix temp
compute_modify add dynamic yes extra 0
compute_modify blockt dynamic yes extra 0
compute_modify subt dynamic yes extra 0

fix flith all gcmc 500 20 10 2 29494 300.0 -2.0 0.01 region r_surf

thermo_style custom step atoms f_flith[1] f_flith[2] f_flith[3] f_flith[4] f_flith[5] f_flith[6] temp epair etotal vol press
thermo 100
thermo_modify lost/bond ignore lost warn

run 50000
unfix flith

fix f_nvtall all nvt temp 300 300 100

thermo_style custom step atoms temp epair etotal press
thermo 1000
thermo_modify lost/bond ignore lost warn
run 150000

log.lammps (91.9 KB)

restartnew2.MD_260000 (150 KB)

The error you mention about atom_modify first group has

nothing to do with atom_modify sort, which is a different option.

There will be no issue with a first group being defined

if you didn’t use the atom_modify first command.

Steve