# Input file for DPD_Shockwave Bubble Test



label		rerun

variable	t index 1 0
variable	R index 50 0
variable	d index 16 0
variable	k index 800 0
variable	latden equal 2.7
variable	dens index 3 0
variable	q index 81 7885 54343 689 414 547 987 368 13 756 0
variable	h index 10 20 0
variable	m equal -55-$h
variable	j equal $m-21.65
variable	v equal $j-$h
variable	g equal $v-50
variable	b equal round(-1*$v)
variable	p index 1 0 
variable	c equal $p*100
variable	PC equal $p*25
log		DPD3h_S$k_$bB$h_P$c_R$Rd$d.lammps
variable	w equal 1
variable	l equal 0.7
variable	RR equal $R+11*$l
variable	Rr equal $R-4*$l
variable	RRw equal $R+10*$l
variable	Rrw equal $R-2*$l
variable	u equal $R+10
variable	o equal $u+1
variable	move equal ramp($g,$v)
variable	e equal $g+1
# ------------------------ INITIALIZATION ----------------------------
units 		lj
dimension	3
boundary	p	p	f
atom_style	full 
read_data	DPD-sphere$Rdd$d.data
change_box	all z final $g 100
mass		1 1
mass		2 1
mass		3 1
mass		4 1
mass		5 1
mass		6 1

communicate single vel yes
#comm_modify	vel yes
# ------------------------ FORCE FIELDS ------------------------------
pair_style      dpd 1.0 1.0 98$q

pair_coeff      1 1 ${PC} 4.5
pair_coeff      1 2 ${PC} 4.5
pair_coeff      2 2 ${PC} 4.5
pair_coeff      1 3 100 4.5
pair_coeff      2 3 100 4.5
pair_coeff      3 3 ${PC} 4.5
pair_coeff      1 4 25 4.5
pair_coeff      2 4 25 4.5
pair_coeff      3 4 100 4.5
pair_coeff      4 4 25 4.5
pair_coeff      1 5 150 4.5
pair_coeff      2 5 150 4.5
pair_coeff      3 5 150 4.5
pair_coeff      4 5 150 4.5
pair_coeff	5 5 0 4.5
pair_coeff      1 6 25 4.5
pair_coeff      2 6 25 4.5
pair_coeff      3 6 100 4.5
pair_coeff      4 6 25 4.5
pair_coeff      5 6 150 4.5
pair_coeff	6 6 25 4.5

bond_style      harmonic
bond_coeff	* 100 $l

angle_style     harmonic

angle_coeff	* 3 180
angle_coeff	2 2.25 120
angle_coeff	3 1.5 120
angle_coeff	4 2.25 120

dihedral_style	none

neighbor	2.0 bin
neigh_modify    check yes every 1 delay 0 page 1000000 one 100000

# ---------------------- SPATIAL DEFINITION --------------------------
region		box block INF INF INF INF INF INF units box
region		slipidi sphere 0 0 0 ${RR} units box side in
region		slipido sphere 0 0 0 ${Rr} units box side out
region		swateri sphere 0 0 0 ${Rrw} units box side in
region		swatero sphere 0 0 0 ${RRw} units box side out
region		lipid intersect 2 slipidi slipido
region		water union 2 swateri swatero
group		lipid region lipid

lattice		fcc ${latden} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
#create_atoms	4 region water
create_atoms	4 region swatero
create_atoms	6 region swateri
group		4 type 4
group		6 type 6
group		water type 4 6

group		sys union 6 lipid

set		type 1 charge 0.0
set		type 2 charge 0.0
set		type 3 charge 0.0
set		type 4 charge 0.0
# ------------------------- SETTINGS ---------------------------------
compute 	peatom all pe/atom
compute		peall all pe
compute 	keatom all ke/atom
compute		keall all ke
compute		temp all temp
compute		astress all stress/atom

compute		lkeall lipid reduce sum c_keatom
compute		lpeall lipid reduce sum c_peatom

compute		x all property/atom x
compute		xm all reduce ave c_x
compute		y all property/atom y
compute		ym all reduce ave c_y
compute		z all property/atom z
compute		zm all reduce ave c_z
compute		xl all reduce min c_x
compute		yl all reduce min c_y
compute		zl all reduce min c_z
compute		xh all reduce max c_x
compute		yh all reduce max c_y
compute		zh all reduce max c_z
compute		vx all property/atom vx
compute		vy all property/atom vy
compute		vz all property/atom vz

compute		lpx lipid property/atom x
compute		lpxm lipid reduce ave c_lpx
compute		lpy lipid property/atom y
compute		lpym lipid reduce ave c_lpy
compute		lpz lipid property/atom z
compute		lpzm lipid reduce ave c_lpz
variable	lpdiffx atom (c_lpx-c_lpxm)^2
variable	lpdiffy atom (c_lpy-c_lpym)^2
variable	lpdiffz atom (c_lpz-c_lpzm)^2
compute		lpxx lipid reduce sum v_lpdiffx
compute		lpyy lipid reduce sum v_lpdiffy
compute		lpzz lipid reduce sum v_lpdiffz

compute		wx 6 property/atom x
compute		wxm 6 reduce ave c_wx
compute		wy 6 property/atom y
compute		wym 6 reduce ave c_wy
compute		wz 6 property/atom z
compute		wzm 6 reduce ave c_wz
variable	wdiffx atom (c_wx-c_wxm)^2
variable	wdiffy atom (c_wy-c_wym)^2
variable	wdiffz atom (c_wz-c_wzm)^2
compute		wxx 6 reduce sum v_wdiffx
compute		wyy 6 reduce sum v_wdiffy
compute		wzz 6 reduce sum v_wdiffz

variable	nall equal count(all)
variable	nlip equal count(lipid)
variable	nw equal count(6)
variable	peall equal c_peall
variable	keall equal c_keall
variable	lpeall equal c_lpeall
variable	lkeall equal c_lkeall
variable	peatom atom c_peatom
variable	keatom atom c_keatom
variable	lstress equal c_lstress
variable	wstress equal c_wstress

variable	Rgl equal sqrt((1/v_nlip)*(c_lpxx+c_lpyy+c_lpzz))
variable	Rgw equal sqrt((1/v_nw)*(c_wxx+c_wyy+c_wzz))

variable	lstr atom sqrt(((c_astress[1]-c_astress[2])^2+(c_astress[2]-c_astress[3])^2+(c_astress[1]-c_astress[3])^2+6*(c_astress[4]^2+c_astress[5]^2+c_astress[6]^2))/2)

compute		lstress lipid reduce sum v_lstr
compute		wstress all reduce sum v_lstr
######################################
timestep 	0.01
velocity	all create 1.0 12$q mom yes rot yes

dump 		1 sys cfg 1 DPD3h_S$k_$bB$h_P$c_R$Rd$d_*.cfg id type xs ys zs c_astress[1] c_astress[2] c_astress[3] &
		c_astress[4] c_astress[5] c_astress[6] v_lstr c_peatom c_keatom c_vx c_vy c_vz type mol
dump_modify 	1 element C O SI H C K
dump_modify	1 pad 5

dump 		2 all cfg 1000 DPD3h_S$k_$bB$h_P$c_R$Rd$d_w*.cfg id type xs ys zs c_astress[1] c_astress[2] c_astress[3] &
		c_astress[4] c_astress[5] c_astress[6] v_lstr c_peatom c_keatom c_vx c_vy c_vz type mol
dump_modify 	2 element C O SI H C K
dump_modify	2 pad 5
thermo_style	custom step lz press v_nall
thermo_modify	lost warn

run		0

thermo		10
fix		out all ave/time 1 $w $w v_keall v_peall v_lkeall v_lpeall v_lstress v_wstress &
		v_Rgl v_Rgw file DPD3h_S$k_$bB$h_P$c_R$Rd$d.txt
fix		1 all nve

change_box	all boundary p p m
region		zlr block INF INF INF INF INF $e units box
region		zhr block INF INF INF INF 99 INF units box
region		delete union 2 zlr zhr
delete_atoms	region delete
fix		piston all wall/reflect zlo v_move zhi EDGE
dump_modify	1 every 1000
run		$k



region		bubble2 sphere -12.5 0 $j $h units box side in
region		bubble3 sphere 12.5 0 $j $h units box side in
region		bubbles union 2 bubble2 bubble3
group		bubble2 region bubbles

group		delete2 intersect bubble2 4

delete_atoms	group delete2


unfix		piston
fix		piston all wall/reflect zlo $v zhi EDGE


dump_modify	1 every 1
dump_modify	2 every 1
run		1
dump_modify	1 every 100
dump_modify	2 every 25
run		250

region		bubble1 sphere 0 0 $m $h units box side in
group		bubble1 region bubble1
group		delete1 intersect bubble1 4
delete_atoms	group delete1

dump_modify	2 every 1
run		1
dump_modify	1 every 100
dump_modify	2 every 50
run		300
dump_modify	1 every 100
dump_modify	2 every 200
#run		1000
dump_modify	1 every 200
dump_modify	2 every 6000
#run		6000
######################################
# SIMULATION DONE
clear
#jump 		DPD-SB2.in
next		q
print 		"Run Done"
next		h
if		$h!=0 then "jump DPD-SB3.in rerun"
variable 	h delete
next		p
if		$p!=0 then "jump DPD-SB3.in rerun"
variable 	p delete
next		k
if		$k!=0 then "jump DPD-SB3.in rerun"
variable 	k delete
print		"Finally Finished"

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