#
# Runs LAMMPS simulation electrolyte solutions
# echo            screen
variable        dcdfreq index 10000 
variable        thermofreq index 100
variable        rand equal 990416
variable        nlayers equal 100
variable        nsteps index 10000000
variable        temp index 250
variable        deltat equal 20
variable        Thot equal (v_temp+v_deltat)
variable        Tcold equal (v_temp-v_deltat)
variable        dofh  equal 1.5925
variable        dofo  equal 2.8151
variable        seedH equal 35543
variable        seedC equal 3356
variable        unused equal 69420
variable        dz equal 4.0
variable        inputname  index startnemd
variable        outputname string nemd1.TH${Thot}.TC${Tcold}
# variable        outputname index run3-113

# Define simulation conditios
units           real
boundary        p p p
atom_style      full

# Forcefield
pair_style      lj/long/tip4p/long long long  1 2 1 1 0.1546 10 10
pair_style      goop
pair_modify     mix arithmetic
kspace_style    pppm/disp/tip4p 1.0e-5
kspace_modify   force/disp/real 0.0001
kspace_modify   force/disp/kspace 0.0002


# Reads in previous configurations
read_data       ${inputname}.data


# Intramolecular interactions
bond_style      hybrid harmonic # OW-HW
bond_coeff      1 harmonic 5.102 0.9572
angle_style     hybrid harmonic # HW-OW-HW
angle_coeff     1 harmonic 2.198 104.52

# Intermolecular interactions: Define pair coefficients
pair_coeff 1 1 0.1852 3.1589 # O O TIP4P/2005 x
pair_coeff 1 2 0.0000 0.0000 # O H TIP4P/2005 x
pair_coeff 2 2 0.0000 0.0000 # H H TIP4P/2005 x

# Define Masses
mass 1 15.9994 # OW
mass 2 1.00794 # HW
# Define groups
group           water type 1 2
group           oxygen type 1
group           hydrogen type 2


#SIMULATION
neighbor        2 bin
neigh_modify    delay 5 every 1

# Thermostats
# Regions
#set up the  at the left end of the box
variable lx equal lx
variable ly equal ly
variable lz equal lz


variable  cbl equal $(-v_lz)/2.0+${dz}
variable  cbh equal ${lz}/2.0-${dz}
region cold1 block INF INF INF INF INF  ${cbl} units box
region cold2 block INF INF INF INF ${cbh} INF   units box

#set up the "cold region" in the middle of the box
#region hot_region block INF INF INF INF ${cbl} ${cbh} units box
region hot_region block INF INF INF INF -${dz} ${dz} units box

#cold regions  unification
region cold_region union 2 cold1 cold2 units box

#Apply thermostats
compute T_hot water temp/region hot_region
compute T_cold water temp/region cold_region
fix thermostatH water temp/csvr $(v_Thot*2/3) $(v_Thot*2/3) $(500*dt) $(v_seedH)
fix thermostatC water temp/csvr $(v_Tcold*2/3) $(v_Tcold*2/3) $(500*dt) $(v_seedC)
fix_modify thermostatH temp T_hot
fix_modify thermostatC temp T_cold

fix cons_mom all momentum 1 linear 1 1 1
compute mymom all momentum
compute mymomwater water momentum


variable function equal 12*sqrt(-1260) # Just for testing
variable pi equal PI

# Reset timestep
reset_timestep 0

timestep        1

fix 2 all nve what 

log ${outputname}.log

fix shakew water rattle 0.00001 200 0 b 1 a 1 # FB Shake Water

fix shakew2 water raatle  0.00001 200 0 b 1 a 1 # FB Shake Water

compute mom all mumentum
thermo          ${thermofreq}

thermo_style    custom step time etotal pe ke temp press pxx pyy pzz vol lx ly lz atoms density f_thermostatH f_thermostatC c_T_hot c_T_cold  c_mymom[*] c_mymomwater[*] f_imaginary
# dump            1 all dcd ${dcdfreq} ${outputname}.dcd
# dump_modify     1 unwrap yes
dump            2 all custom ${dcdfreq} ${outputname}.dump id type x y z vx vy vz ix iy iz
dump_modify     2 append no

# Compute profile water and ions
compute lwater water chunk/atom bin/1d z lower $(lz/v_nlayers) units box
fix profw water  ave/chunk 1 $(v_nsteps) $(v_nsteps) lwater  density/number density/mass  temp file ${outputname}water.profiles adof 2

# Compute prodile Oxygen and hydrogen
compute loxygen oxygen chunk/atom bin/1d z lower $(lz/v_nlayers) units box
fix tempoxy oxygen ave/chunk 1 $(v_nsteps) $(v_nsteps) loxygen density/number density/mass temp adof $(v_dofo) file ${outputname}oxygen.profiles
compute lhydrogen hydrogen chunk/atom bin/1d z lower $(lz/v_nlayers) units box
fix temphyd hydrogen ave/chunk 1 $(v_nsteps) $(v_nsteps) lhydrogen density/number density/mass temp adof $(v_dofh) file ${outputname}hydrogen.profiles

run   $(v_nsteps)
write_data      ${outputname}.data pair ij

compute  1 all rdf

# Invalid fix styles
compute  1 all temp

fix  12 all what

fix nvt all nvt water temp "hello"

invalid_command

run  2 

fix # Here to show what happens with in-complete commands.