atom_style	full
units	real
dimension	3
boundary	p p f
neighbor	2.0 bin
# neigh_modify	every 1 delay 0 check yes

variable Ctype equal	1
variable Otype equal	2
variable Htype equal	3
variable Natype equal	4
variable Cltype equal	5

variable waterbond equal 1
variable waterangle equal 1

# Define interaction parameters
pair_style lj/cut/coul/long 10.0
#pair_modify tail yes
kspace_style pppm 1.0e-4
kspace_modify slab 3.0 #define empty space in z direction for non-periodicity

bond_style harmonic
angle_style harmonic
#dihedral_style harmonic
#improper_style harmonic


read_data	geometry.data

#Coefficients

pair_coeff ${Htype} ${Htype} 0.0 0.0 #H-H
pair_coeff ${Otype} ${Htype} 0.0 0.0 # H-O
pair_coeff ${Otype} ${Otype} 0.15535 3.166 # O-O
pair_coeff ${Ctype} ${Ctype} 0.07029 3.55 #C-C
pair_coeff ${Ctype} ${Htype} 0.0 0.0 # C-H
pair_coeff ${Ctype} ${Otype} 0.104497 3.358 # C-O

pair_coeff ${Otype} ${Natype} 0.14211 2.661 # O-NA
pair_coeff ${Otype} ${Cltype} 0.12464 3.7435 # CL-O
pair_coeff ${Ctype} ${Natype}  0.09559 2.853 # C-NA
pair_coeff ${Ctype} ${Cltype} 0.08384 3.9355 # C-CL
pair_coeff ${Natype} ${Cltype} 0.11402 3.2385 # CL-NA

pair_coeff ${Htype} ${Natype}  0.0 0.0 # H-NA
pair_coeff ${Htype} ${Cltype}  0.0 0.0 # CL-H

pair_coeff ${Natype} ${Natype}  0.130 2.156 # NA-NA
pair_coeff ${Cltype} ${Cltype} 0.1 4.321 # CL-CL

#groups
group carbon	type ${Ctype}
group ox        type ${Otype}
group hy        type ${Htype}
group Na        type ${Natype}
group Cl        type ${Cltype}
group all_water	union hy ox

# Setting charges
set group ox charge -0.8476  #SPE/C
set group hy charge  0.4238
set group Na charge  1.0
set group Cl charge  -1.0
set group carbon charge 0.0

variable membranel equal -1.77241 #membrane z positions
variable membraner equal 22.2980

#group membrane
variable zmin	equal ${membranel}-0.5
variable zmax	equal ${membraner}+0.5
region	membraneregion	block		INF INF INF INF ${zmin} ${zmax} units box
group		membrane	region	membraneregion


#new regions for pressure calculation
region left_box block 5.0 25.0 5.0 25.0 -25.0 -5.0 units box
region right_box block 5.0 25.0 5.0 25.0 26.0 46.0 units box

group left_water region left_box 
group right_water region right_box 

#bond,angle,dihedral coefficients
bond_coeff * 0.0 0.0
angle_coeff * 0.0 0.0 
#dihedral_coeff * 0.0 1 0 
#improper_coeff * 0.0 0 

bond_coeff  ${waterbond} 0.0 1.0 # H2O bond (SPE/C)
angle_coeff ${waterangle} 0.0 109.47 # H2O angle (SPE/C)


fix wall_bound all wall/reflect zlo EDGE zhi EDGE 

fix membrane_freeze membrane setforce 0.0 0.0 0.0  #freeze membrane

#minimize
#fix minz all box/relax x 1.0 y 1.0 vmax 0.001
#minimize 1e-4 1e-6 1000 1000
#unfix minz

fix water_shak all_water shake 1.0e-4 500 0 b ${waterbond} a ${waterangle}

#equilibration using NPT
fix fix_npt all_water npt temp 300 300 100 x 1.0 1.0 1000.0 y 1.0 1.0 1000.0 couple none

#how to use npt to maintain the same pressure in two separate #regions(pure water region and saltwater region), do we need two #npt fixes, one for each region??
#If so, how to define the pressure fo the two npt commands for #the two separate regions?


compute     mobile all_water temp
velocity    all_water create 300.0 12345678 dist gaussian

variable voll equal 8000.0  #volume of the small regions taken

#compute pressure in a small box in the pure water region
compute allperatom all stress/atom mobile 
compute allp all reduce/region left_box sum c_allperatom[1] c_allperatom[2] c_allperatom[3]
variable ppvalue equal -(c_allp[1]+c_allp[2]+c_allp[3])/(3*${voll})   

#compute pressure in a small box in the salt water region
compute allperatom2 all stress/atom mobile 
compute allp2 all reduce/region right_box sum c_allperatom2[1] c_allperatom2[2] c_allperatom2[3]
variable ppvalue2 equal -(c_allp2[1]+c_allp2[2]+c_allp2[3])/(3*${voll})



thermo 1000
thermo_style custom step temp density press v_ppvalue v_ppvalue2
thermo_modify temp mobile flush yes
timestep 1.0

dump            equil1 all custom 5000 dump_equil.lammpstrj id type x y z

run             500000
undump   equil1

write_restart   file_equil.restart

unfix fix_npt

#final dynamics in NVT
fix fix_nvt all_water nvt temp 300 300 100

dump            final all custom 1000 dump_dynamics.lammpstrj id type x y z
run             100000
undump final

write_restart   file_dynamics.restart