# Simulation with toy surfactant molecules

units		lj
atom_style	full
dimension	3
boundary p p f

# Define a square lattice in a narrow region to create wall particles
lattice		fcc 2.5
region		box block 0 42.0 0 42.0 0 0.1 units box
create_box	4 box bond/types 2 angle/types 2 extra/angle/per/atom 2 extra/bond/per/atom 2 extra/special/per/atom 6
create_atoms	4 box

kspace_style pppm 6.0e-2
kspace_modify slab 3.0

dielectric 40.0

# Set charge for wall particles
group wall type 4

mass 1 1.0
mass 2 1.0
mass 3 1.0
mass 4 1.0

# Change the box size, so that we can add surfactants and solvent above the wall
change_box all z final 0 200.0 units box
region    abovebox block 0 42.0 0 42.0 2.0 195.0 units box

# The geometry of the surfactant molecule is dfined in 'surf.mol'
molecule surfactant surf.mol.lj

# Define potential parameters

pair_style lj/cut/coul/long 1.12246204831 3.0
pair_coeff	1 1   1.0 2.5 2.80615512078 # Surf Head, Surf Head
pair_coeff  1 2   1.0 1.0 # Surf Head, Surf Tail
pair_coeff  1 3   1.0 1.0 # Surf Head, Counterion
pair_coeff  1 4   1.0 1.0 # Surf Head, Wall
pair_coeff  2 2   1.7 1.0 2.5 # Surf Tail, Surf Tail
pair_coeff  2 3   1.0 1.0 # Surf Tail, Counterion
pair_coeff  2 4   1.0 1.0 # Surf Tail, Wall
pair_coeff  3 3   1.0 1.0 # Counterion, Counterion
pair_coeff  3 4   1.0 1.0 # Counterion, Wall
pair_coeff  4 4   0.0 0.0 0.0 # Wall, Wall

# Shift the potential 
pair_modify shift yes 

# Write Restart
restart 10000 rst.1 rst.2 

bond_style harmonic
bond_coeff 1 80.0 1.0
bond_coeff 2 80.0 1.0

angle_style harmonic
angle_coeff 1 20.0 180.0
angle_coeff 2 32.0 180.0

# Variables for computing the necessary number of counterions
variable headchg equal 0.0
variable wallchg equal 0.0
variable counterionchg equal 0.0

# Randomly position surfactant molecules above the wall
create_atoms 0 random 442 3717 abovebox overlap 1.0 mol surfactant 6354 units box
set   type 1 charge -4

# Randomly position surfactant counterions
create_atoms 3 random 1090 3539 abovebox overlap 1.0 units box
set 	type 3 charge 4.0
set		group  wall charge -1.0

# Set velocities
velocity	all create 1.44 7501 loop geom

# Declare wall BC at top of the box
fix top all wall/harmonic  zhi EDGE 10.0 0.0 4.0 units box

# Fluid (mobile) particles
group  fluid type 1 2 3

timestep 0.01
fix 	1 fluid nve/limit 0.1
run		10000
unfix 1

fix   1 fluid langevin 1.0 1.0 10.0 65948
fix		2	fluid	nve

run 1000000000