boundary f f p 
units real
atom_style full
bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long 15
kspace_style pppm/electrode 1e-7
kspace_modify wire 3 

variable unit_length equal 22.2404 
variable unit_water_num equal 1000
variable unit_ion_num equal 10
variable unit_rep equal 1


variable elect_dist equal 16.28
variable Lx equal ${unit_rep}*${unit_length}
variable displace equal ${elect_dist}+${Lx}
variable Ltop equal ${elect_dist}+${Lx}
variable Lend equal ${Ltop}+${elect_dist}
variable Llimit equal ${Lend}+1

variable Lside equal 3.663e+01
variable creation_low equal 1
variable creation_high equal ${Lside}-1

region box block 0.0 ${Llimit} 0.0 ${Lside} 0.0 ${Lside}

create_box 6 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 &
extra/angle/per/atom 1 &
extra/special/per/atom 5

region bot block 0 ${elect_dist} EDGE EDGE EDGE EDGE

print ${elect_dist}
print ${Lx}
print ${Ltop}
print ${Lend}
print ${Llimit}
include SPCFw-NaCl.ff
molecule water H2O-SPCFw.mol 

lattice	fcc 4.07 # Gold electrode settings

create_atoms 6 region bot units lattice #displace to top
group top_electrode type 6 
displace_atoms top_electrode move ${displace} 0 0 units box

variable x_min_top equal bound(top_electrode,xmin)
variable x_max_top equal bound(top_electrode,xmax)

print "X boundaries of top electrode ${x_min_top} ${x_max_top}"

create_atoms 5 region bot units lattice
group bot_electrode type 5
variable x_min_bot equal bound(bot_electrode,xmin)
variable x_max_bot equal bound(bot_electrode,xmax)

print "X boundaries of bot electrode ${x_min_bot} ${x_max_bot}"

region bot_electrode block ${x_min_bot} ${x_max_bot} EDGE EDGE EDGE EDGE
region top_electrode block ${x_min_top} ${x_max_top} EDGE EDGE EDGE EDGE

variable x0_water equal ${x_max_bot}
variable x1_water equal ${x_min_top}



lattice none 1.0
region creation block ${x0_water} ${x1_water} ${creation_low} ${creation_high} ${creation_low} ${creation_high}


variable Nwater equal ${unit_rep}*${unit_water_num}
variable Nna equal ${unit_rep}*${unit_ion_num}
variable Ncl equal ${unit_rep}*${unit_ion_num}

create_atoms 0 random ${Nwater} 311121 creation mol water 13131 overlap 1.0 maxtry 1000
create_atoms 3 random ${Nna} 311121 creation overlap 1.0 maxtry 1000
create_atoms 4 random ${Ncl} 311121 creation overlap 1.0 maxtry 1000
set type 3 charge 1
set type 4 charge -1
write_data init.data

group electrolyte type 1 2 3 4
group water type 1 2
variable ke equal ke
variable evdwl equal evdwl
variable ecoul equal ecoul
variable elong equal elong
variable etotal equal etotal

compute ctemp electrolyte temp
fix etotal_out all ave/time 100 1 100 v_ke v_etotal v_evdwl v_ecoul v_elong file min.energy

thermo 100
thermo_modify temp ctemp
thermo_style custom step temp ke pe v_etotal v_evdwl v_ecoul v_elong


timestep 1e-2

fix mywall electrolyte wall/reflect ylo EDGE yhi EDGE
fix mynvt electrolyte nvt temp 1e-5 1e-5 0.1 


dump mindump all custom 1000 min.lammpstrj id mol type q x y z vz vy vz

run 20000

unfix mynvt
timestep 1
fix mynvt electrolyte nvt temp 1e-5 1e-5 1
run 50000
write_data min.data