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  		5
variable pistontype equal   6
variable NAtype equal   	4

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

read_data pristine.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}  # all carbon atoms
group ox        type ${Otype}
group hy        type ${Htype}
group na        type ${NAtype}
group cl        type ${CLtype}

# 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

region rigidmembrane1_zone cylinder z 15.0 15.0 15.0 INF INF
group in_rigidmembrane_1zone region rigidmembrane1_zone
group membrane1_free subtract membrane1_atoms in_rigidmembrane_1zone

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

group thermostat_target union saltwater membrane1_free piston1 piston2

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 C H
# Lennard-Jones: sigmaAB = (1/2)(sigmaAA + sigmaBB), epsilonAB = (epsilonAAepsilonBB)1/2
# Syntax: pair_coeff atom_i atom_j epsilon sigma
pair_coeff ${Ctype} ${CLtype} lj/cut/tip4p/long 0.031702 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.12613 3.2793 # C-O
pair_coeff ${Ctype} ${pistontype} lj/cut/tip4p/long 0.019208 3.3749 # C-PIS
pair_coeff ${Ctype} ${NAtype} lj/cut/tip4p/long 0.12027 2.8293 # C-Na
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.046549 4.1617 # Cl-O
pair_coeff ${CLtype} ${pistontype} lj/cut/tip4p/long 0.0070888 4.2572 # CL-PIS
pair_coeff ${CLtype} ${NAtype} lj/cut/tip4p/long 0.044388 3.7117 # Cl-Na
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-PIS
pair_coeff ${Htype} ${NAtype} lj/cut/tip4p/long 0 0 # H-Na
pair_coeff ${Otype} ${Otype} lj/cut/tip4p/long 0.1852 3.1589 # O-O
pair_coeff ${Otype} ${pistontype} lj/cut/tip4p/long 0.028203 3.2545 # O-PIS
pair_coeff ${Otype} ${NAtype} lj/cut/tip4p/long 0.1766 2.7089 # O-Na
pair_coeff ${pistontype} ${pistontype} lj/cut/tip4p/long 0.004295 3.35 # PIS-PIS
pair_coeff ${pistontype} ${NAtype} lj/cut/tip4p/long 0.026894 2.8045 # PIS-Na
pair_coeff ${NAtype} ${NAtype} lj/cut/tip4p/long 0.1684 2.2589 # Na-Na

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 ----------------
fix pistonfreeze bothpistons setforce 0.0 0.0 0.0
compute selective_thermostat thermostat_target temp
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 equilib.lammpstrj
run     50000 

write_restart pristine.*.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
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 2500000 dynamics.restart

dump    fulldump     all atom 1000 dynamics.lammpstrj
run     10000