######################################
# LAMMPS INPUT SCRIPT
# HAp deformation

# ---------- Setup Variables --------------------- 
variable atoms_coordinate string "system.lmp"

variable fname string "friction"
variable potname string "reaxff.txt"

variable relax_time_step equal 0.1
variable initial_temp_relax equal 30

log log.${fname}.txt

######################################
# INITIALIZATION
units real
dimension 3
boundary p p p
atom_style full

######################################
# DATA FILE
read_data ${atoms_coordinate}												
region 1 block INF INF INF 0.07215 INF INF

######################################
#GROUP
group 1 molecule 1
group 2 molecule 2
group 3 molecule 3
group 4 molecule 4
group 5 molecule 5
group 6 molecule 6

group 7 molecule 1 2 3
group 8 molecule 4 5 6
group 9 molecule 1 2 3 4 5
group 10 region 1
######################################
# DEFINE INTERATOMIC POTENTIAL
include ${potname}
special_bonds coul 0.0 1.0 1.0
neighbor 2.0 bin
neigh_modify every 1 delay 1 check yes page 10000000000
neigh_modify one 1000000000 

######################################
# DEFINE COMPUTES 

compute 1 all pe/atom	
										
compute 2 all displace/atom	

compute 3 3 temp
compute 4 3 stress/atom 3
compute 5 3 reduce ave fx fy fz

compute 6 4 temp
compute 7 4 stress/atom 6
compute 8 4 reduce ave fx fy fz

compute 9 2 temp
compute 10 5 temp											
compute 13 9 temp
compute 11 all stress/atom 13
compute 12 all stress/atom 13 virial		#exclude KE
######################################
#MINIMIZE(because use NVE)
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
minimize 1.0e-4 1.0e-6 100 1000
min_style cg

#####################################
#PRESSING

reset_timestep 0
timestep 0.1
velocity all create ${initial_temp_relax} 12345 mom yes rot yes dist gaussian

#FIX_COMMANDS
fix 2 all nvt temp 30 30 $(100*dt) drag 0.2

fix 4 1 setforce 0 0 0
  
fix 5 3 ave/time 1 100 100 c_5[1] c_5[2] c_5[3] file force_HAp_interface.txt
fix 6 4 ave/time 1 100 100 c_8[1] c_8[2] c_8[3] file force_sphere_interface.txt
fix 7 all ave/time 1 100 100 c_3 c_6 c_9 c_10 file temp_ave.txt

#fix 11 6 rigid single
fix 12 6 move linear 0 0 -0.005

# Set thermo output
thermo 100
thermo_style custom cpu step lx ly lz pxx pyy pzz pxy pxz pyz pe ke etotal press vol density temp
thermo_modify lost ignore
thermo_modify lost/bond ignore

#DUMP
dump 1 all cfg 100 punwall_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_11[1] c_11[2] c_11[3] c_11[4] c_11[5] c_11[6] c_12[1] c_12[2] c_12[3] c_12[4] c_12[5] c_12[6]
dump_modify 1 element Ca P O H Ti
dump 2 3 cfg 100 H_i_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_4[1] c_4[2] c_4[3] c_4[5] c_4[6]
dump_modify 2 element Ca P O H Ti
dump 4 4 cfg 100 Ti_i_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_7[1] c_7[2] c_7[3] c_7[5] c_7[6]
dump_modify 4 element Ca P O H Ti
dump 5 8 cfg 100 Ti_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_7[1] c_7[2] c_7[3] c_7[5] c_7[6]
dump_modify 5 element Ca P O H Ti
dump 6 10 cfg 100 cutted_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_11[1] c_11[2] c_11[3] c_11[4] c_11[5] c_11[6] c_12[1] c_12[2] c_12[3] c_12[4] c_12[5] c_12[6]
dump_modify 6 element Ca P O H Ti

#CYCLES
run 1280

unfix 12
undump 1
undump 2
undump 4
undump 5
undump 6
#####################################
#MOVING

timestep ${relax_time_step}

#FIX_COMMANDS

fix 12 6 move linear 0.0005 0 0

# Set thermo output
thermo 1000
thermo_style custom cpu step lx ly lz pxx pyy pzz pxy pxz pyz pe ke etotal press vol density temp
thermo_modify lost ignore
thermo_modify lost/bond ignore

#DUMP
dump 1 all cfg 1000 unwallm_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_11[1] c_11[2] c_11[3] c_11[4] c_11[5] c_11[6] c_12[1] c_12[2] c_12[3] c_12[4] c_12[5] c_12[6]
dump_modify 1 element Ca P O H Ti
dump 2 3 cfg 1000 H_im_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_4[1] c_4[2] c_4[3] c_4[5] c_4[6]
dump_modify 2 element Ca P O H Ti
dump 4 4 cfg 1000 Ti_im_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_7[1] c_7[2] c_7[3] c_7[5] c_7[6]
dump_modify 4 element Ca P O H Ti
dump 5 8 cfg 1000 Tim_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_7[1] c_7[2] c_7[3] c_7[5] c_7[6]
dump_modify 5 element Ca P O H Ti
dump 6 10 cfg 1000 cuttedm_*.cfg mass type xsu ysu zsu q c_1 c_2[1] c_2[2] c_2[3] c_2[4] fx fy fz c_11[1] c_11[2] c_11[3] c_11[4] c_11[5] c_11[6] c_12[1] c_12[2] c_12[3] c_12[4] c_12[5] c_12[6]
dump_modify 6 element Ca P O H Ti

#CYCLES
run 200000

unfix 12

###################################
#TAKE UP

#FIX_COMMANDS
fix 12 6 move linear 0 0 0.005

# Set thermo output
thermo 1000
thermo_style custom cpu step lx ly lz pxx pyy pzz pxy pxz pyz pe ke etotal press vol density temp
thermo_modify lost ignore
thermo_modify lost/bond ignore

#CYCLES
run 100000
######################################
# SIMULATION DONE
print "All done"