# Toy system for testing "sorted dumps", 2018-08-06:
#
# this system has a slab-geometry containing two hex layers 
# of spheres at zmin and zmax, both layers are charged 
# either +E  (lower layer) or -E (upper layer), generating 
# a field gradient from bottom to top in the neutral but in-
# homogeneous-charged system. Vapor from molecules with dipole 
# moments (SPC water) will possibly align dipoles along the 
# gradient and condense out. Hopefully, the slab/pppm method 
# in LAMMPS with 4A slab buffer will be able to handle such 
# a system ...
#

units           real
atom_style      full
boundary        p p f
lattice         hcp  3.0  

# µ_ex of -11.25 leads to a pressure very close to the 
# saturation pressure of the rigid SPC w/SHAKE model
variable temp    equal 300
variable muex    equal -11.25
variable dspl    equal 0.0
variable zborder equal 0.5

variable xbox    equal 24
variable ybox    equal v_xbox/sqrt(3.0)
variable zbox    equal v_xbox/3.0-2.0
variable zboxeps equal v_zbox+1.0e-5
variable zm3     equal v_zbox-v_zborder

region           full   block  0 ${xbox}  0 ${ybox}  0 ${zboxeps}    units lattice
region           center block  0 ${xbox}  0 ${ybox}  ${zborder} ${zm3} units lattice
create_box   4   full   bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2

# create top- and bottom layers
region           btmlayer  block  0  ${xbox}  0 ${ybox}     0       0    units lattice
region           toplayer  block  0  ${xbox}  0 ${ybox}  ${zbox} ${zbox} units lattice
create_atoms 3   region btmlayer  units lattice
create_atoms 4   region toplayer  units lattice
set type     3   mol 1 charge  0.025
set type     4   mol 2 charge -0.025

# preload molecule for monte-carlo and start up with some molecules
molecule         h2o  ./spc.shake.mol
create_atoms   	 0 random 33 54321 center  mol h2o 464563    units box
group water      type 1 2

# fix the walls
group frozen     type 3 4
fix   walls      frozen   setforce 0.0 0.0 0.0

# fairly standard potential stuff
pair_style     lj/cut/coul/long 10.0 15.0
kspace_style   pppm      1.0e-5
kspace_modify  slab      4
bond_style     harmonic
angle_style    harmonic

# I made those (*-3, *-4) up, I confess
pair_coeff      1 1  0.15535 3.1656
pair_coeff      1 3  0.15535 3.1656
pair_coeff      1 4  0.15535 3.1656
pair_coeff      3 3  0        0       
pair_coeff      3 4  0        0       
pair_coeff      4 4  0        0       
pair_coeff      * 2  0        0       
bond_coeff      1    1000       1.0   
angle_coeff     1     100       109.47

mass 1 15.9994   # OW
mass 2 1.0       # HW
mass 3 24.0      # W+
mass 4 24.0      # W-

# put the configuration files into subdirectory in order
# to keep order in the main directory
shell "mkdir -p config/; rm config/*.*"
dump        conf  all custom  1000 config/d.*.lcd  id mol type element  x y z 

# sort or no sort -- this is the question: here be dragons!
# dump_modify conf  sort  id
dump_modify conf  pad   8
dump_modify conf  pbc   yes
dump_modify conf  element OW HW W+ W-

# - - - - MD settings - - - -
neighbor        2.0 bin
neigh_modify    every 1 delay 1 check yes
velocity       	all create ${temp} 54654
timestep        1.0

# - - - - prepare extra data file for water properties - - - - 
variable   s  equal step
variable   m  equal "count(water) / 3"
variable   d  equal density
variable   p  equal press
variable   t  equal temp
fix        Dens  water print 1000 "$s $m $t $p $d" file density.log title "step nmol temp press density" screen no 

# - - - - correct d.o.f. for rigid molecule - - - - - 
variable    tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)

variable    oxygen atom "type==1"
variable    hydrogen atom "type==2"
group       oxygen dynamic all var oxygen
group       hydrogen dynamic all var hydrogen
variable    nO equal count(oxygen)
variable    nH equal count(hydrogen)
variable	tacc equal f_mc[2]/(f_mc[1]+0.1)
variable	iacc equal f_mc[4]/(f_mc[3]+0.1)
variable	dacc equal f_mc[6]/(f_mc[5]+0.1)
variable	racc equal f_mc[8]/(f_mc[7]+0.1)

# - - - - - set up fixes for simulation - - - - - - 
thermo_style custom step temp v_nO v_nH atoms press density pe ke v_iacc v_dacc v_tacc v_racc
thermo       200
fix          nvt water  nvt temp ${temp} ${temp} 100
fix          shk water  shake 0.0001 50 0  b 1 a 1  mol h2o
fix          mc  water  gcmc 100 100 0 0 54341 ${temp} ${muex} 0 mol h2o  region center tfac_insert ${tfac} shake shk
compute_modify  thermo_temp dynamic/dof yes
compute_modify  nvt_temp dynamic/dof yes

# - - - - go ahead - - - - -
run          40000000