units           metal                   # eV, Ang, ps
variable        kcal_mol equal 43.36e-3 # eV
atom_style      full
newton on
boundary        p p p

bond_style      harmonic
angle_style     harmonic
dihedral_style  none
improper_style  none

read_data       water.dump

velocity        all set 0.0 0.0 0.0
group           water type 1 2

# TIP4P/Ice
bond_coeff      1 44.0 0.9572  # 1000 kcal/mol ~= 44 eV
angle_coeff     1 44.0 104.52
set             type 1 charge -1.1794
set             type 2 charge  0.5897

thermo_style      custom step temp epair etotal ndanger
thermo            10

neighbor          4.0 bin
neigh_modify      every 1 delay 0 check yes
timestep          1e-3  # ps


##### MINIMIZING BOND LENGTH AND ANGLE TO EQUILIBRIUM ONES ######
pair_style        lj/cut 10.0
pair_coeff        *   *    0.00000             0.0000

min_style cg
minimize 1e-15 1e-15 100000 10000000

write_data        01_fired.dump
#################################################################


##### INTRODUCING TIP4P INTERACTION #############################
pair_style        lj/cut/tip4p/cut 1 2 1 1 0.1577 10.0 8.5
pair_coeff        *   *    0.00000             0.0000
pair_coeff        1   1  $(0.21084*v_kcal_mol) 3.1668

fix               1 all viscous 0.01
fix               3 all shake 1e-4 20 100 b 1 a 1 t 1 2
fix               4 all nve

run 10000
write_data        02_tip4p.dump
#################################################################