# Film of 512 water molecules, with fix evaporate.

# Input variables.
variable run    string water02   # output log
variable ts     equal  1         # timestep
variable cutoff equal  12        # cutoff for real-space interactions
variable ti     equal  300       # initial temperature
variable tf     equal  400       # equilibrium temperature
variable p      equal  1.        # equilibrium pressure
variable skin   equal  2.        # Skin distance for the neighbour list.
variable cl     equal  4         # correlation length for averaging
variable s      equal  50        # sample interval for averaging
variable prod   equal  50000     # Production steps
variable weps   equal  0.4       # Wall epsilon
variable wsigma equal  3.15      # Wall sigma

# PBC 
boundary p p f

# Force field setup.
units           real
atom_style      full
bond_style      harmonic
angle_style     harmonic
pair_style      lj/cut/coul/cut ${cutoff} 
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 angle yes dihedral yes

# Neighbour list.
neighbor        ${skin} bin
neigh_modify    check yes

# Create a random seed.
variable r_tmp equal random(10000,1000000,48519)
variable r format r_tmp %.0f

# Extra stuff
processors * * 1
neigh_modify delay 0 every 10 check yes

# ----------------- Atom Definition Section -----------------
read_data "water02.data"

# ----------------- Settings Section -----------------
mass 1  15.9994
mass 2   1.0079

pair_coeff  1 1 0.102 3.188
pair_coeff  2 2 0.0 0.0

bond_coeff  1  450.0 0.9572
angle_coeff 1  55.0 104.52

# Walls along the non-periodic direction.
fix ref all wall/reflect zhi EDGE
fix bot all wall/lj126 zlo EDGE ${weps} ${wsigma} ${cutoff}


# Apply the SHAKE algorithm to hydrogen atoms.
fix SHK all shake .0001 10 0 m 1.0079 a 1

# Derived variables.
# thermo interval.
variable d       equal ${cl}*$s
# dump interval ("sample interval" frames).
variable dcycle  equal ${prod}/$s
variable tcouple equal ${ts}*100
variable pcouple equal ${ts}*1000

# Output.
thermo          $d
thermo_style custom step etotal evdwl ecoul elong ebond eangle atoms &
ke pe temp press vol density cpu
thermo_modify flush yes
compute_modify thermo_temp dynamic/dof yes

# Trajectory.
dump TRJ all dcd ${dcycle} ${run}.dcd
dump_modify TRJ unwrap yes

# Set the velocity to the final temperature.
if "$(is_defined(variable,ti))" then "variable t0 equal ${ti}" else "variable t0 equal ${tf}"
velocity all create ${t0} $r dist gaussian mom yes rot yes

# Heating the system.
fix NVT all nvt temp ${ti} ${tf} ${tcouple}
fix_modify NVT temp thermo_temp
run ${prod}
unfix NVT

# Evaporation
region evap_space block EDGE EDGE EDGE EDGE 35 EDGE
fix EVP all evaporate $(v_d*10) 3 evap_space $r molecule yes
fix NVT all nvt temp ${tf} ${tf} ${tcouple}
fix_modify NVT temp thermo_temp
run $(v_prod*2)

# New equilibration
unfix EVP
run ${prod}