# molecular mechanics simulation of 3d indentation into MgO

dimension     3
boundary      p p s

units         metal
atom_style    charge

neighbor      3.0 bin
neigh_modify  delay 5

# create geometry: Mg O

lattice       custom 4.217 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 &
              basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.0 0.5 basis 0.0 0.5 0.5 &
              basis 0.5 0.5 0.5 basis 0.0 0.0 0.5 basis 0.0 0.5 0.0 basis 0.5 0.0 0.0
region        box block 0 42.17 0 42.17 0 42.17 units box
create_box    2 box
create_atoms  2 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 &
                    basis 5 2 basis 6 2 basis 7 2 basis 8 2

mass	      1   24.305
mass	      2   16.0

#charges for Mg and O 
group    magnesiumatoms type  1
group    oxygenatoms type  2

set group magnesiumatoms charge 2.0
set group oxygenatoms charge -2.0

# hybrid potential for MgO

pair_style   buck/coul/long 12.0
pair_coeff   1 2 buck/coul/long 2143.768 0.26734 0.0 12.0
pair_coeff   2 2 buck/coul/long 25.31 0.6937 32.32 16.0
pair_coeff   1 1 buck/coul/long 0.0 1.0 0.0  12.0

#kspace_style   pppm 1.0e-4
#kspace_modify slab 3.0

# define groups
region      sunatracte block 0 42.17 0 42.17 0 21.085 units box
group	      bottom region sunatracte
#group	    mobile subtract all bottom

# fix the substrate (bottom) in space

fix	      1 bottom setforce 0.0 0.0 0.0 

# minimize the energy of the system before indentation

thermo        10000000000

min_style     cg         
min_modify    dmax 0.15
minimize      1.0e-10 1.0e-6 30000 30000

write_restart minimization.restart

reset_timestep 0

# compute various values
compute       1 all centro/atom fcc

dump          1 all cfg 100000000 atom*.cfg id type xs ys zs c_1
dump_modify   1 element Mg O

# define varaibles
variable      k equal 1234.689846
variable      k1 equal 246.9379692
variable      inix equal 21.085
variable      iniy equal 21.085
variable      iniz equal 52.17
variable      inc equal 0.2
variable      dec equal 0.1

# minimize the energy of the system with indenter

fix           4 all indent $k sphere ${inix} ${iniy} ${iniz} 10.0 units box
fix_modify    4 energy yes

thermo_style  custom step pe etotal ecoul epair ebond emol press vol f_4 f_4[1] f_4[2] f_4[3]

min_style     cg        
min_modify    dmax 0.15
minimize      1.0e-10 1.0e-6 30000 30000

# indentation process

variable      i loop 1000
label         loading
print         "i = $i"
if            "$i > 200" then "jump creatmgostr.in loading_ends"
variable      height equal "(v_iniz-v_inc*v_i)"

fix	      4 all indent ${k1} sphere ${inix} ${iniy} v_height 10.0 units box
fix_modify    4 energy yes

min_style     cg
min_modify    dmax 0.15
minimize      1.0e-10 1.0e-6 30000 30000

next          i
jump          creatmgostr.in loading

label         loading_ends

write_restart loading.restart

# unloading starts

variable      j loop 1000
label         unloading
print         "j = $j"
if            "$j > 200" then "jump creatmgostr.in unloading_ends"
variable      height equal "(v_iniz-v_inc*200+v_dec*v_j)"

fix	      4 all indent ${k1} sphere ${inix} ${iniy}  v_height 10.0  units box
fix_modify    4 energy yes

min_style     cg
min_modify    dmax 0.15
minimize      1.0e-10 1.0e-6 30000 30000

next          j
jump          creatmgostr.in unloading

label         unloading_ends