units           metal
variable a equal 3.59571

boundary        p p p

############ build base block ########################################
lattice         fcc $a
region          box block -3 3 -3 3 -5 0

create_box      2 box
create_atoms    1 box


pair_style hybrid eam/alloy/opt lj/cut 1.0 
pair_coeff * * eam/alloy/opt ./potential Cu NULL 

pair_coeff 2 2 none
pair_coeff 1 2 lj/cut 0.0001 2 3.0

mass 2 0.000002

thermo 1
thermo_style custom step temp pe press pxx pyy pzz 
minimize 0.0 1.0e-2 40000 40000
fix 1 all box/relax  aniso 0.0
minimize 1.0e-12 1.0e-8 40000 40000
unfix 1

#####################################################################

######################## therm base atoms ########################### 
############### this is just to stop wild movement of
############### the slab atoms, which would cause
############### misleading behaviour
############### thermalise first with PBC, then create void in Z

thermo 1000
thermo_style custom step temp pe press pxx pyy pzz 

###### thermalise under full periodic boundaries ####################
fix 1 all npt temp 300 300 0.1 aniso 0.0 0.0 8
velocity all create 300 4572
run 2000
velocity all create 300 4572
fix 1 all npt temp 300 300 0.1 aniso 0.0 0.0 2
run 850
unfix 1

###### store box variables for indenter ########
variable lid equal zhi
variable low equal zlo
variable upzero equal ${lid}+5.0
variable upone equal ${lid}+10.8
variable uptwo equal ${upone}+12.0
variable down equal ${low}-5.0
variable ceil equal ${lid}+50
variable btop equal ${low}+3.0

###### thermalise again with PBC removed in Z #########################
change_box all z final ${down} ${ceil} boundary p p f units box
fix 1 all npt temp 300 300 0.5 x 0.0 0.0 1 y 0.0 0.0 1
velocity all create 300 4572
run 170
velocity all create 300 4572
run 1000
unfix 1

###### set base group to stop drift ###################################
region base block EDGE EDGE EDGE EDGE EDGE ${btop} units box
group base region base
group mobile subtract all base

###### build indenter #################################################
################## do this by creating a region of
################## two sets of offset planes

lattice sc 1
region top block EDGE EDGE EDGE EDGE ${upzero} EDGE units box

####################### berkovich #####################################
region p1 plane 1 1 ${upone} -0.417074 0.0      0.908872 side in units box 
region p2 plane 1 1 ${upone} 0.208537 -0.361197 0.908872 side in units box 
region p3 plane 1 1 ${upone} 0.208537  0.361197 0.908872 side in units box 
region ind_outer intersect 3 p1 p2 p3 

region p1a plane 1 1 ${uptwo} -0.417074 0.0	 0.908872 side in units box 
region p2a plane 1 1 ${uptwo} 0.208537 -0.361197 0.908872 side in units box 
region p3a plane 1 1 ${uptwo} 0.208537  0.361197 0.908872 side in units box 
region ind_inner intersect 3 p1a p2a p3a side out

region indenter intersect 2 ind_outer ind_inner

create_atoms 2 region indenter
group indenter region indenter 
displace_atoms indenter move 0 0 -8.183 units box

group slab type 1

velocity base set 0 0 0
reset_timestep 0

compute f indenter group/group slab

variable zpos equal bound(indenter,zmin)
variable ci equal 2.817-v_zpos

################## you can vary this number 
################## to change the target load.
variable tl equal 10

################## this variable sets the 
################## indenter velocity (via fix 7)
################## according to the load
variable iv equal (c_f[3]-v_tl)+(0.0001*(c_f[3]-v_tl)*(c_f[3]-v_tl)*(c_f[3]-v_tl))
################## to prove a point
################## you can set the velocity as a constant
################## by commenting the line above
################## and uncommenting this line below
################## (you will get a 'lost atoms' warning)
#variable iv equal -0.5

################## only integrate the 'slab' atoms
################## because 'fix move' updates indenter positions
fix 4 slab nve

fix 5 base setforce 0 0 0
fix 6 indenter setforce 0 0 0
fix 7 indenter move variable NULL NULL NULL NULL NULL v_iv

thermo 1000
thermo_style custom step press v_ci c_f c_f[3] v_iv 
thermo_modify lost warn
neighbor 7.0 bin
neigh_modify exclude type 2 2 every 1 delay 2 check yes one 4000

run 100000