# Initialization

units real
dimension 3
atom_style full
boundary p p p
read_data packmol_waterbox.dat

# create groups
group hy type 1      # H_water
group ox type 2      # O_water

# set charges
set group hy charge 0.4238   # SPC/E model partial charge
set group ox charge -0.8476  # SPC/E model partial charge

# potential parameters
pair_style lj/cut/coul/long 8.5
pair_coeff 1 1*2 0.0000 0.0000 # LJ epsilon, sigma of OH, HH = 0.0
pair_coeff 2 2 0.1553 3.1660 # SPC/E model LJ epsion and sigma for O-O interaction in real units

bond_style harmonic
bond_coeff 1 7.7 1.000       # SPC/E model equilibrium bond length

angle_style harmonic
angle_coeff 1 1.8 109.470    # SPC/E model equilibrium bending angle

kspace_style pppm 1.0e-5


neighbor 2.0 bin
neigh_modify delay 0 every 1 delay 10 check yes


variable  T equal 298.0           # equilibirum temperature in Kelvin
variable  P equal 1.0             # pressure in atmosphere
variable  dt equal 1.0            # time step in femto second
variable  tdamp equal 10.0*${dt}
variable  pdamp equal 100.0*${dt}

thermo          100
thermo_style    custom step pe etotal press temp
thermo_modify   norm no

fix 1 all box/relax iso 1.0 vmax 0.001
minimize 1.0e-5 1.0e-5 10000 10000
min_modify line quadratic
write_restart restart.min.*

reset_timestep 0
timestep ${dt}

velocity     all create $T 4928459 dist gaussian
fix          SHAKE all shake 1e-5 20 0 b 1 a 1          # it is important to use shake algorithm for water molecules. Shake algorithm makes water molecules behave as rigid bodies.
fix          NPT all npt temp $T $T ${tdamp} iso $P $P ${pdamp} # try to use npt after minimization to bring the density of water molecules to its equilibrium value at the specified temperature
dump         XYZ all xyz 500 equil_npt.*.xyz
run 20000
write_restart restart.equil_npt.*
unfix NPT
unfix SHAKE
undump XYZ