#####################################################
#                                                   #
# CaF2 Nanoindentation LAMMPS INPUT SCRIPT          #
# Potential used: Buckingham Potential              #
# Filename: in.CaF2_BUCK_Nanoindentation.txt        #
# Author: Janel Chua, 2017                          #
#                                                   #
# Execute the script through:                       #
# lmp_exe < in.CaF2_BUCK_Nanoindentation.txt        #
#                                                   #
#####################################################


# ------------------------ INITIALIZATION ----------------------------

variable 	fname index CaF2[100]10x10x10.txt
neighbor	3.0 bin
neigh_modify	delay 0 every 1 check yes

# ----------------------- ATOM DEFINITION ----------------------------

units 	metal # ps, Angstrom
echo 		both
dimension	3
boundary	p	p	p
atom_style 	charge
read_data	${fname}
replicate	1 1 1

# ------------------------ FORCE FIELDS ------------------------------

pair_style buck/coul/long 12 # cutoff
kspace_style pppm 1.0e-4
pair_coeff 1 1 0.00   1   0.00 
pair_coeff 1 2 1717.441   0.287   0.102
pair_coeff 2 2 2058.994  0.252   16.703 

# potential coeffs obtained from high-temp phase diagram of CaF2 from classical atomistic simulations
# rho cannot be 0

#------------------------- SETTINGS ---------------------------------

compute csym all centro/atom fcc
compute peratom all pe/atom
compute eng all pe/atom
compute eatoms all reduce sum c_eng

# ------------------------ Energy Minimization -----------------------

reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
dump 		1 all atom 1 Min.lammpstrj 
thermo 10 
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
min_style cg 
minimize 1e-25 1e-25 10000 10000 
unfix 1
undump 	1
timestep	0.005

# ------------------------ Equilibration -----------------------

reset_timestep	0
velocity all create 300 12345 mom yes rot no # specifies the thermal velocity
fix 1 all npt 	temp 300 300 1 iso 1 1 1 drag 0.0  

dump 		2 all atom 250 Equilibrium.lammpstrj 

# Set thermo output
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp

# Run for at least 10 picosecond (assuming 1 fs timestep)
run 10000 # This is 20 femtosecond # you have to equilibrate it till the structure stops moving
unfix	1
undump	2
timestep	0.000001

# ------------------------rigid boundary ------------------------

region     1 block 0 55 0 7 0 55 
group      anvil1 region 1
region     2 block 0 7 0 55 0 55 
group      anvil2 region 2
region     3 block 48 55 0 55 0 55 
group      anvil3 region 3
group      anvil union anvil1 anvil2 anvil3
group      mobile subtract all anvil

# ------------------------ Nanoindentation -----------------------

reset_timestep	0
dump       1 all atom 1000 dump.indent83d_unload.dump.lammpstrj
log log5050_indent8_3d_voidless_unload.dat

# initial velocities
compute	   new mobile temp
velocity   mobile create 300 482748 temp new
fix	   1 mobile  nvt temp 300.0  300.0 0.05
#fix	   1 all nve
fix	   2 anvil setforce 0.0 0.0 0.0

# assigning velocity to the indenter in y direction/loading direction
variable   y equal "90-step*dt*0.1" # 135+8
print      "y is $y" 

# indenter position and radius at onset of loading
fix	   	4 mobile indent 1000.0 sphere 25 v_y 25 20.0 units box
thermo	   	100
thermo_style    custom step temp  v_y f_4[1] f_4[2] f_4[3] 
run	   	100000