# =================== Variable_settings =======================

variable	thermo_step		equal	10 
variable	rl_run_step		equal	500 
variable	sp_run_step		equal	500 
variable	sp_num_dump		equal	20 

variable	c_velocity		equal	1 

variable	v_element	string	'c'

# =================== Model_settings ==========================

newton		on
units		metal
boundary	p p p

atom_style	full

read_data	cell-diamond.data

replicate	10 10 30

# =========== group =============

region c1 block INF INF INF INF INF 6.0
region c2 block INF INF INF INF 6.0 $((zhi-zlo)*2/3)
region c3 block INF INF INF INF $((zhi-zlo)*2/3) $((zhi-zlo)*2/3+10.0)
region c4 block INF INF INF INF $((zhi-zlo)*2/3+10.0) INF
region c5 block INF INF INF INF $((zhi-zlo)*2/3+40.0) INF
region c6 block INF INF INF INF $(zhi-10.0) INF

region p1 plane $((xhi-xlo)/2+xlo) $((yhi-ylo)/2+ylo) $((zhi-zlo)*2/3+10.0) -1  0 1.8 side in
region p2 plane $((xhi-xlo)/2+xlo) $((yhi-ylo)/2+ylo) $((zhi-zlo)*2/3+10.0)  1  0 1.8 side in
region p3 plane $((xhi-xlo)/2+xlo) $((yhi-ylo)/2+ylo) $((zhi-zlo)*2/3+10.0)  0  1 1.8 side in
region p4 plane $((xhi-xlo)/2+xlo) $((yhi-ylo)/2+ylo) $((zhi-zlo)*2/3+10.0)  0 -1 1.8 side in
region p5 plane $((xhi-xlo)/2+xlo) $((yhi-ylo)/2+ylo) $(zhi) 0 0 1 side out
region p_top intersect 5 p1 p2 p3 p4 p5

group c1 region c1
group c2 region c2
group c3 region c3
group c4 region c4
group c5 region c5
group c6 region c6

group p1 region p1
group p2 region p2
group p3 region p3
group p4 region p4
group p5 region p5
group p_top region p_top

group c4_del subtract c4 p_top

delete_atoms group c3
delete_atoms group c4_del
delete_atoms group c5
delete_atoms group c6

# write_dump p1 custom dump.p1.lammpstrj id element x y z
# write_dump p2 custom dump.p2.lammpstrj id element x y z
# write_dump p3 custom dump.p3.lammpstrj id element x y z
# write_dump p4 custom dump.p4.lammpstrj id element x y z
# write_dump p5 custom dump.p5.lammpstrj id element x y z
write_dump p_top custom dump.p_top.lammpstrj id element x y z
write_dump all custom dump.all.lammpstrj id element x y z

group			mobile union p_top
group			base   union c1 c2 c3

# =================== Force_field_settings ====================

pair_style		tersoff 
pair_coeff		* * SiC.tersoff C 

# =================== Thermo_style settings ===================

compute			f_c0 p_top reduce sum fx fy fz 

variable 		f_all equal sqrt(c_f_c0[1]*c_f_c0[1]+c_f_c0[2]*c_f_c0[2]+c_f_c0[3]*c_f_c0[3])
variable 		f_z   equal sqrt(c_f_c0[3]*c_f_c0[3])

compute 		mobilePress mobile stress/atom NULL

compute 		coord_atom all coord/atom cutoff 2.1 *

thermo_style	custom step cpu cpuremain temp etotal c_f_c0[*] v_f_all
thermo			${thermo_step}
thermo_modify	flush yes

# =================== Dump settings ===========================

dump			strj	all		custom	$(v_sp_run_step/v_sp_num_dump) all.lammpstrj id element x y z fx fy fz
dump_modify		strj	flush	yes		element ${v_element} sort id

dump			strj2	c3		custom	$(v_sp_run_step/v_sp_num_dump) c3.lammpstrj id element x y z fx fy fz c_coord_atom[*]
dump_modify		strj2	flush	yes		element ${v_element} sort id

dump			strj3	p_top		custom	$(v_sp_run_step/v_sp_num_dump) p_top.lammpstrj id element x y z fx fy fz
dump_modify		strj3	flush	yes		element ${v_element} sort id

dump 			strj4	mobile custom $(v_sp_run_step/v_sp_num_dump) mobilePress.lammpstrj id element x y z  c_mobilePress[*]
dump_modify		strj4	flush	yes		element ${v_element} sort id

# shell mkdir img

# dump			image_x	all		image	$(v_sp_run_step/v_sp_num_dump) ./img/img_x.*.jpg type type view 90 90  size 1000 1400 zoom 2.0 axes yes 0.1 0.01 box yes 0.003
# dump			image_y	all		image	$(v_sp_run_step/v_sp_num_dump) ./img/img_y.*.jpg type type view 90  0  size 1000 1400 zoom 2.0 axes yes 0.1 0.01 box yes 0.003

# =================== Minimize run ============================

fix			freeze	mobile	setforce 0.0 0.0 0.0
min_style	cg
minimize	1.0e-4 1.0e-6 100 1000
unfix		freeze

reset_timestep	0

# =================== Relaxed run ============================

fix			npt_1			base npt temp 300.0 300.0 100.0 iso 0.0 0.0 1000.0
run			${rl_run_step}
unfix		npt_1

# ================ inject run ==========================

velocity	all	set 0.0 0.0 0.0

fix			f_nve_base		base nvt temp 300.0 300.0 100.0

fix			f_static 		c1 setforce 0.0 0.0 0.0

fix			f_rgd_up		mobile rigid/nvt group 1 mobile temp 300.0 300.0 100.0 torque * off off off

fix 		m_v 			mobile move linear 0.0 0.0 $(-v_c_velocity)

timestep	0.005

run			${sp_run_step}

unfix		m_v

# =================== deinject run ============================

fix 		m_v 			mobile move linear 0.0 0.0 $(v_c_velocity)

run			${sp_run_step}

unfix		m_v

