atom_style  full
units   real
dimension 3
boundary p p p

variable Ctype equal  	 	1
variable CLtype equal  		2
variable Htype equal   		3
variable Otype equal  		4
variable pistontype equal   5
variable NAtype equal   	6
variable zCotype equal 		7
variable zChtype equal 		8
variable zOtype equal 		9
variable zHctype equal 		10
variable zHotype equal 		11

variable seed equal 		1
variable pressure equal 	0.05
variable 1atm equal 		-0.000041
variable waterbond equal    1
variable waterangle equal   1

read_data single_system.dat
variable natoms equal "count(all)"

# Define interaction parameters
pair_style hybrid lj/cut/tip4p/long ${Otype} ${Htype} ${waterbond} ${waterangle}  0.1546 13.0 airebo 3.0 0 1
pair_modify tail yes
bond_style harmonic
angle_style harmonic
dihedral_style none

# ------------- Regions and Groups ---------------- 
group carbon    type ${Ctype} ${pistontype}
group ox        type ${Otype}
group hy        type ${Htype}
group na        type ${NAtype}
group cl        type ${CLtype}
group c_hydrogen type ${zChtype}
group h_hydrogen type ${zHctype}
group c_hydroxyl type ${zCotype}
group o_hydroxyl type ${zOtype}
group h_hydroxyl type ${zHotype}
group hydroxyl_group union c_hydroxyl o_hydroxyl h_hydroxyl
group hydrogen_group union c_hydrogen h_hydrogen
group functional_group union hydrogen_group hydroxyl_group

# Defining the positions of all four carbon planes
variable zpiston1 equal -50
variable zmembrane1 equal 0
variable zpiston2 equal "20"

# Defining the piston and the membrane carbons
variable zmin equal ${zpiston1}-1.0
variable zmax equal ${zpiston1}+1.0
region piston1zone   block INF INF INF INF ${zmin} ${zmax} units box
variable zmin equal ${zmembrane1}-1.0
variable zmax equal ${zmembrane1}+1.0
region  membrane1zone  block INF INF INF INF ${zmin} ${zmax} units box
variable zmin equal ${zpiston2}-1.0
variable zmax equal ${zpiston2}+1.0
region piston2zone   block INF INF INF INF ${zmin} ${zmax} units box
group piston1            region piston1zone             # whole piston1 (carbon)
group piston2          region piston2zone               # whole piston1 (carbon)
group bothpistons          union piston1 piston2
group membrane1_atoms region membrane1zone 
group totalmembrane     union membrane1_atoms
group not_hy subtract all hy

# defining the subset of each membrane that is allowed to move
variable free_radius equal 16.5
region freemembrane1_zone cylinder z 15.0 15.0 ${free_radius} INF INF
group in_freemembrane_1zone region freemembrane1_zone
group membrane1_free intersect in_freemembrane_1zone membrane1_atoms
group membrane1_rigid subtract membrane1_atoms membrane1_free

# defining water groups
group water     union ox hy
group saltwater union water na cl
group notsaltwater subtract all saltwater

# Defining the group of atoms that go into the Nose-Hoover thermostat
group thermostat_target union saltwater membrane1_free piston1 piston2

# setting hydrogen and hydroxyl charges
set group c_hydrogen charge -0.115
set group h_hydrogen charge 0.115
set group c_hydroxyl charge 0.2
set group o_hydroxyl charge -0.64
set group h_hydroxyl charge 0.44

# Setting water charges
set group ox charge -1.0484
set group hy charge  0.5242
set group na charge  1.
set group cl charge -1.

# ------------- Coefficients ---------------- 
# AIREBO:
pair_coeff * * airebo CH.airebo_real C NULL NULL NULL NULL NULL NULL C NULL H NULL
# Lennard-Jones: sigmaAB = (1/2)(sigmaAA + sigmaBB), epsilonAB = (epsilonAA*epsilonBB)^1/2
# Syntax: pair_coeff atom_i atom_j epsilon sigma
pair_coeff ${Ctype} ${CLtype} lj/cut/tip4p/long 0.0317022 4.2821 # C-CL
pair_coeff ${Ctype} ${Htype} lj/cut/tip4p/long 0 0 # C-H
pair_coeff ${Ctype} ${Otype} lj/cut/tip4p/long 0.118238 3.28203 # C-O
pair_coeff ${Ctype} ${pistontype} lj/cut/tip4p/long 0.0859 3.3997 # C-piston
pair_coeff ${Ctype} ${NAtype} lj/cut/tip4p/long 0.120273 2.8293 # C-NA
pair_coeff ${Ctype} ${zCotype} lj/cut/tip4p/long 0.0777095 3.47485 # C-zCo
pair_coeff ${Ctype} ${zOtype} lj/cut/tip4p/long 0.115388 3.23485 # C-zO
pair_coeff ${Ctype} ${zHotype} lj/cut/tip4p/long 0 0 # C-zHo
pair_coeff ${CLtype} ${CLtype} lj/cut/tip4p/long 0.0117 5.1645 # CL-CL
pair_coeff ${CLtype} ${Htype} lj/cut/tip4p/long 0 0 # CL-H
pair_coeff ${CLtype} ${Otype} lj/cut/tip4p/long 0.0436369 4.16443 # CL-O
pair_coeff ${CLtype} ${pistontype} lj/cut/tip4p/long 0.0317022 4.2821 # CL-piston
pair_coeff ${CLtype} ${NAtype} lj/cut/tip4p/long 0.0443878 3.7117 # CL-NA
pair_coeff ${CLtype} ${zCotype} lj/cut/tip4p/long 0.0286794 4.35725 # CL-zCo
pair_coeff ${CLtype} ${zChtype} lj/cut/tip4p/long 0.0231991 4.07475 # CL-zCh
pair_coeff ${CLtype} ${zOtype} lj/cut/tip4p/long 0.0425852 4.11725 # CL-zO
pair_coeff ${CLtype} ${zHctype} lj/cut/tip4p/long 0.0187662 3.79225 # CL-zHc
pair_coeff ${CLtype} ${zHotype} lj/cut/tip4p/long 0 0 # CL-zHo
pair_coeff ${Htype} ${Htype} lj/cut/tip4p/long 0 0 # H-H
pair_coeff ${Htype} ${Otype} lj/cut/tip4p/long 0 0 # H-O
pair_coeff ${Htype} ${pistontype} lj/cut/tip4p/long 0 0 # H-piston
pair_coeff ${Htype} ${NAtype} lj/cut/tip4p/long 0 0 # H-NA
pair_coeff ${Htype} ${zCotype} lj/cut/tip4p/long 0 0 # H-zCo
pair_coeff ${Htype} ${zChtype} lj/cut/tip4p/long 0 0 # H-zCh
pair_coeff ${Htype} ${zOtype} lj/cut/tip4p/long 0 0 # H-zO
pair_coeff ${Htype} ${zHctype} lj/cut/tip4p/long 0 0 # H-zHc
pair_coeff ${Htype} ${zHotype} lj/cut/tip4p/long 0 0 # H-zHo
pair_coeff ${Otype} ${Otype} lj/cut/tip4p/long 0.16275 3.16435 # O-O
pair_coeff ${Otype} ${pistontype} lj/cut/tip4p/long 0.118238 3.28203 # O-piston
pair_coeff ${Otype} ${NAtype} lj/cut/tip4p/long 0.165551 2.71163 # O-NA
pair_coeff ${Otype} ${zCotype} lj/cut/tip4p/long 0.106964 3.35717 # O-zCo
pair_coeff ${Otype} ${zChtype} lj/cut/tip4p/long 0.0865246 3.07468 # O-zCh
pair_coeff ${Otype} ${zOtype} lj/cut/tip4p/long 0.158828 3.11718 # O-zO
pair_coeff ${Otype} ${zHctype} lj/cut/tip4p/long 0.0699912 2.79218 # O-zHc
pair_coeff ${Otype} ${zHotype} lj/cut/tip4p/long 0 0 # O-zHo
pair_coeff ${pistontype} ${pistontype} lj/cut/tip4p/long 0.0859 3.3997 # piston-piston
pair_coeff ${pistontype} ${NAtype} lj/cut/tip4p/long 0.120273 2.8293 # piston-NA
pair_coeff ${pistontype} ${zCotype} lj/cut/tip4p/long 0.0777095 3.47485 # piston-zCo
pair_coeff ${pistontype} ${zChtype} lj/cut/tip4p/long 0.0628602 3.19235 # piston-zCh
pair_coeff ${pistontype} ${zOtype} lj/cut/tip4p/long 0.115388 3.23485 # piston-zO
pair_coeff ${pistontype} ${zHctype} lj/cut/tip4p/long 0.0508487 2.90985 # piston-zHc
pair_coeff ${pistontype} ${zHotype} lj/cut/tip4p/long 0 0 # piston-zHo
pair_coeff ${NAtype} ${NAtype} lj/cut/tip4p/long 0.1684 2.2589 # NA-NA
pair_coeff ${NAtype} ${zCotype} lj/cut/tip4p/long 0.108805 2.90445 # NA-zCo
pair_coeff ${NAtype} ${zChtype} lj/cut/tip4p/long 0.0880136 2.62195 # NA-zCh
pair_coeff ${NAtype} ${zOtype} lj/cut/tip4p/long 0.161561 2.66445 # NA-zO
pair_coeff ${NAtype} ${zHctype} lj/cut/tip4p/long 0.0711958 2.33945 # NA-zHc
pair_coeff ${NAtype} ${zHotype} lj/cut/tip4p/long 0 0 # NA-zHo
pair_coeff ${zCotype} ${zCotype} lj/cut/tip4p/long 0.0703 3.55 # zCo-zCo
pair_coeff ${zCotype} ${zChtype} lj/cut/tip4p/long 0.0568665 3.2675 # zCo-zCh
pair_coeff ${zCotype} ${zOtype} lj/cut/tip4p/long 0.104386 3.31 # zCo-zO
pair_coeff ${zCotype} ${zHctype} lj/cut/tip4p/long 0.0460003 2.985 # zCo-zHc
pair_coeff ${zCotype} ${zHotype} lj/cut/tip4p/long 0 0 # zCo-zHo
pair_coeff ${zChtype} ${zOtype} lj/cut/tip4p/long 0.0844393 3.0275 # zCh-zO
pair_coeff ${zChtype} ${zHotype} lj/cut/tip4p/long 0 0 # zCh-zHo
pair_coeff ${zOtype} ${zOtype} lj/cut/tip4p/long 0.155 3.07 # zO-zO
pair_coeff ${zOtype} ${zHctype} lj/cut/tip4p/long 0.0683045 2.745 # zO-zHc
pair_coeff ${zOtype} ${zHotype} lj/cut/tip4p/long 0 0 # zO-zHo
pair_coeff ${zHctype} ${zHotype} lj/cut/tip4p/long 0 0 # zHc-zHo
pair_coeff ${zHotype} ${zHotype} lj/cut/tip4p/long 0 0 # zHo-zHo

kspace_style pppm/tip4p 1.0e-4

# Bond and angles coeffs
# Syntax: angle_coeff N spring_constant theta_0
bond_coeff * 0.0 0.0 # Zero by default
angle_coeff * 0.0 0.0 # Zero by default

# Water
bond_coeff  ${waterbond} 0.0 0.9572 # H2O bond (TIP4P-2005)
angle_coeff ${waterangle} 0.0 104.52  # H2O angle (TIP4P-2005)

neighbor        2.0 bin
neigh_modify	every 1 delay 10 check yes

# ------------- Setup ----------------
# For the equilibration phase, keep the piston rigid.
fix pistonfreeze bothpistons setforce 0.0 0.0 0.0

# Make sure the thermostat is looking at the relevant atoms by defining a new temperature calculation that only looks at the moveable atoms
compute selective_thermostat thermostat_target temp
# And require LAMMPS to use this compute when doing anything related to temperature
thermo_modify temp selective_thermostat

# ------------- Minimization ----------------        
thermo_style one       
thermo 10
run 0

fix freezewater saltwater setforce 0.0 0.0 0.0
minimize 1.0e-4 1.0e-6 1000 1000
unfix freezewater

fix relax all box/relax x 0.0 y 0.0
minimize 1.0e-4 1.0e-6 100 1000
unfix relax
unfix pistonfreeze

# ------------- Equilibration ----------------
fix pistonkeep bothpistons setforce 0.0 0.0 NULL
fix piston1thrust piston1 aveforce NULL NULL 0.01 # squeezing the water between two pistons at about 150 MPa each
fix piston2thrust piston2 aveforce NULL NULL -0.01
fix 1 water shake 1.0e-4 100 0 b ${waterbond} a ${waterangle}
fix NVTequilib thermostat_target nvt temp 300 300 50
timestep 0.5

thermo 1000
dump    equilibdump not_hy atom 100 single.equilib.lammpstrj
run     1000

write_restart single.*.restart
print "Equilibration completed."
unfix NVTequilib
undump equilibdump
    
# ------------- Dynamics --------
unfix piston1thrust
unfix piston2thrust

fix piston1push piston1 aveforce NULL NULL ${pressure}
fix piston2push piston2 aveforce NULL NULL ${1atm}
fix NVT thermostat_target nvt temp 300 300 50

thermo 1000

# Counting the number of waters and ions in the feed region
region feedzone block INF INF INF INF INF ${zmembrane1} units box 
variable num_feed_waters equal "count(ox,feedzone)"
variable num_feed_na equal "count(na,feedzone)"
variable num_feed_cl equal "count(cl,feedzone)"

thermo_style    custom step temp etotal vol press v_num_feed_waters v_num_feed_na v_num_feed_cl 
restart 2500 single.dynamics.restart
dump    fulldump     all atom 1000 single.dynamics.lammpstrj
run     2500