# input script for SWNT simulation in LAMMPS
#-------10--------20--------30--------40--------50--------60--------70--------80

# Variables
variable input index water-model-418-molecules.txt

# INITIALIZATION
units		    real
boundary            s s s
atom_style          full

processors 3 3 2

# ATOM DEFINITION

read_data           ${input}

# FORCE FIELD

pair_style          reax/c NULL
pair_coeff          * * ffield.reax H O

# SETTINGS

neighbor	    1.0 bin
neigh_modify	    delay 5

# ------------------------------ RELAXATION  -----------------------------------

# ------------------------------------------------------------------------------
variable peng equal pe
variable keng equal ke
variable temp equal temp
variable Pr equal press
variable istep equal step

group hydrogen type 1
group oxygen type 2

# ------- calculate mean square distance ---------------------------------------
compute 1 oxygen msd 
compute 2 hydrogen msd
variable msd_o equal c_1[4]
variable msd_h equal c_2[4]
fix 1 oxygen print 100 "${istep},${msd_o}" file msd-o.txt screen no
fix 2 hydrogen print 100 "${istep},${msd_h}" file msd-h.txt screen no

# ------- calculate RDF --------------------------------------------------------
compute 10 all rdf 50 1 1
compute 11 all rdf 50 1 2
compute 12 all rdf 50 2 2

compute 21 all coord/atom 1.2 1 2
compute 22 all reduce sum c_21[1]
compute 23 all reduce sum c_21[2]

variable NH equal count(hydrogen)
variable NO equal count(oxygen)

variable Hcoord equal c_22/v_NH
variable Ocoord equal c_23/v_NO


fix 10 all ave/time 10 10 50000 c_10 file water-rdf11.rdf mode vector
fix 11 all ave/time 10 10 50000 c_11 file water-rdf12.rdf mode vector
fix 12 all ave/time 10 10 50000 c_12 file water-rdf22.rdf mode vector

fix 13 all print 1000 "${istep},${Hcoord},${Ocoord}" file coordination.txt screen no

# ------- calculate mass density -----------------------------------------------
fix 20 oxygen ave/spatial 10 10 50000 z lower 0.5 density/mass file density-o.profile units box
fix 21 hydrogen ave/spatial 10 10 50000 z lower 0.5 density/mass file density-h.profile units box

# ------------------------------------------------------------------------------
thermo_style custom step cpu temp pe vol pzz

# output for VMD
#dump 		1 all custom 5000 dump.water id type xs ys zs
dump 		2 all custom 5000 dump.water-coord id type x y z

fix 30 all print 1000 "${istep}, ${peng}, ${keng}, ${temp},${Pr}" file water-etp-data.txt screen no

thermo 100
# ------------------------------------------------------------------------------

timestep 0.1

fix 40 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

fix 41 all wall/reflect xlo 0.0 xhi 25.0 ylo 0.0 yhi 25.0 zlo 0.0 zhi 20.0 units box

fix 42 all reax/c/species 2 5 1000 species.out cutoff 1 2 0.35

fix 60 all nve
run 100000
unfix 60

velocity all create 300.0 12345 mom yes rot no
fix 60 all nvt temp 300.0 300.0 100.0
run 100000


print "Water relaxation is complete!"