#LAMMPS INPUT FOR TIP4P WATER
echo screen
dimension 3
units real
atom_style full

boundary p p p
variable Nrun equal 12500000
variable Ndump equal ${Nrun}/1000
variable Ntermo equal ${Nrun}/1000
variable Nf        equal 100
variable StepRescale equal 2000
variable Nr equal 10
########################
#variables for density, etc.
variable watMoleMass equal 18.0153 # /(g/mol)
variable nAvog equal 6.0221415e23 # Avogadro's number
variable watMoleculeMass equal (${watMoleMass}/${nAvog}) # /(g/molecule)
variable A3_in_cm3 equal 1e-24 # Angstrom^3 in cm^3
variable nAtoms equal atoms
variable nMolecules equal v_nAtoms/3
########################

variable create_temp string no
variable temp_rescale string no
variable nvt string no
variable donadio string no
variable npt string no
variable npzAT string yes
variable nhchains string no

#READ SECTION
#read_data water_test.lll
read_restart restart.wat
#read_data water2.dat

#timestep of 2 fs
timestep 2.0


#TEMPERATURE AND PRESSURE
include 'temperature'

#PAIR POTENTIAL PARAMETERS
include 'interactions'

special_bonds  lj/coul 0.0 0.0 0.0


#VELOCITY CREATION
if  "${create_temp}==yes" then "velocity all create ${Temp} 1234" &
"velocity all zero linear"

#NEIGHBOR SECTION
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes one 50000

#MAINTAIN RIGIDITY OF WATER MOLECULES
fix constrain all shake 1.0e-4 100 0 b 1 a 1

#INTEGRATE SECTION
if "${temp_rescale}==yes" then "fix verlet all nve" &
"fix trescale all temp/rescale ${StepRescale} ${Temp} ${Temp} 0.0001 1.0"

if "${npt}==yes" then "fix verlet all nve" &
"fix trescale all temp/csvr ${Temp} ${Temp} ${StepRescale} 123321" &
"fix prescale all press/berendsen iso ${Press} ${Press} ${StepRescale}"

if "${nvt}==yes" then "fix termostat all nvt temp ${Temp} ${Temp} ${StepRescale}"

if "${donadio}==yes" then "fix verlet all nve" &
"fix trescale all temp/csvr ${Temp} ${Temp} ${StepRescale} 123321" 

if "${npzAT}==yes" then "fix integrate all npt temp ${Temp} ${Temp} ${StepRescale} z ${Press} ${Press} ${StepRescale}"

if "${nhchains}==yes" then "fix integrate all npt temp ${Temp} ${Temp} ${StepRescale} iso ${Press} ${Press} ${StepRescale}"
fix removeMomentum all momentum 1 linear 1 1 1 


compute T all temp
fix TempAve all ave/time ${Nf} ${Nr} ${Ntermo} c_T

variable P equal press
fix PressAve all ave/time ${Nf} ${Nr} ${Ntermo} v_P

compute PE all pe pair kspace
variable PE_Mol equal c_PE/v_nMolecules
fix PEAve_Mol all ave/time ${Nf} ${Nr} ${Ntermo} v_PE_Mol

variable Dens equal v_nMolecules*${watMoleculeMass}/(vol*${A3_in_cm3})
fix DensAve all ave/time ${Nf} ${Nr} ${Ntermo} v_Dens

thermo ${Ntermo}
thermo_style custom step temp f_TempAve press f_PressAve f_PEAve_Mol f_DensAve 
thermo_modify flush yes



dump trj all custom ${Ndump} wat.lammpstrj id mol type x y z #xu yu ix iy
dump_modify trj sort id 


run ${Nrun}

write_restart restart.wat
#write_data water2.dat nocoeff