# ------------------------ INITIALIZATION ----------------------------
clear
units		real
dimension	3
newton      off 
boundary	p	p	f

neighbor 2.0 bin
neigh_modify every 1 delay 3

atom_style	full
bond_style	harmonic
angle_style	harmonic


pair_style	lj/cut/coul/long 10.0
kspace_style        pppm 1.0e-6
kspace_modify       slab 3.0

# ----------------------- ATOM DEFINITION ----------------------------
read_data pore0.2full.lmpdat
print "READ DATA DONE"



pair_coeff 1 1 	0.1499 4.0447 # Cl-Cl
pair_coeff 2 2 	0.0460 0.4# H-H
pair_coeff 3 3 	0.04687 2.4299 # Na-Na
pair_coeff 4 4 	0.1521 3.15070000 # O-O
pair_coeff 5 5 	0.07000000 3.550000 # C-C
pair_coeff 6 6 	0.07000000 3.550000 # C-C

bond_coeff 1 450 0.9572 # O-H
angle_coeff 1 	55 104.52

print "All CO-EFFICIENTS DONE"


group G_chlorine 	type 1
group G_hydrogen 	type 2
group G_sodium 		type 3
group G_oxygen 		type 4
group lowerE 		type 5			#lower_electrode
group upperE 		type 6			#upper_electrode

group water 		type 2 4
group salt		type 1 3
group All_liquid 	type 1 2 3 4
group All_electrode 	type 5 6


print "Grouping DONE"

#----------------------------minimize energy------------------------------

thermo 100	
thermo_style multi
neigh_modify 	exclude group lowerE lowerE
neigh_modify 	exclude group upperE upperE

min_style cg
minimize 1e-4 1e-4 2500 10000
print "Minimization DONE"
write_data	Min.dat

#--------------------------Simulation--------------------------------------


run_style       verlet
timestep 	1

velocity All_electrode set 0.0 0.0 0.0

fix zwalls All_liquid wall/reflect zlo EDGE zhi EDGE

dump dump1 all custom 100 yt.dump.* id type x y z ix iy iz vx vy vz

fix 	myforce All_electrode setforce 0.0 0.0 0.0
fix     myshake water shake 1.0e-4 10 1 b 1 a 1

compute comp_oxygen G_oxygen chunk/atom bin/1d z lower 2 units box
fix 	fix_oxygen G_oxygen ave/chunk 1 1 1000 comp_oxygen density/number norm all file density_Oxygen.profile

compute comp_hydrogen G_hydrogen chunk/atom bin/1d z lower 2 units box
fix 	fix_hydrogen G_hydrogen ave/chunk 1 1 1000 comp_hydrogen density/number norm all file density_hydrogen.profile

compute comp_sodium G_sodium chunk/atom bin/1d z lower 2 units box
fix 	fix_sodium G_sodium ave/chunk 1 1 1000 comp_sodium density/number norm all file density_sodium.profile

compute comp_chlorine G_chlorine chunk/atom bin/1d z lower 2 units box
fix 	fix_chlorine G_chlorine ave/chunk 1 1 1000 comp_chlorine density/number norm all file density_chlorine.profile

compute comp_water water chunk/atom bin/1d z lower 2 units box
fix 	fix_water water ave/chunk 1 1 1000 comp_water density/number norm all file density_water.profile

compute comp_salt salt chunk/atom bin/1d z lower 2 units box
fix 	fix_salt salt ave/chunk 1 1 1000 comp_salt density/number norm all file density_salt.profile

compute comp_All_liquid All_liquid chunk/atom bin/1d z lower 2 units box
fix 	fix_All_liquid All_liquid ave/chunk 1 1 1000 comp_All_liquid density/number norm all file density_All_liquid.profile



fix 1 All_liquid nvt temp 298.0 373.0 100.0
run 10000
unfix 1

fix 2 All_liquid nvt temp 373.0 373.0 100.0
run 10000
unfix 2

fix 3 All_liquid nvt temp 373.0 298.0 100.0
run 10000
unfix 3

fix 4 All_liquid nvt temp 298.0 298.0 100.0
run 100000
unfix 4



######################################
# SIMULATION DONE
print "All done"