###################################
# LAMMPS input file for TIP4P/2005f water model
###################################

## Define variables in units real
variable eps_OO equal 0.1852
variable sigma_OO equal 3.1644
variable eps_HH equal 0.000000
variable sigma_HH equal 0.000000
variable eps_OH equal 0.000000
variable sigma_OH equal 0.000000

variable Temp equal 293
variable press equal 1

variable seed1 equal 12345 # Seed for random generation

variable dt_i equal 0.2 # Timestep initializing
variable Nti equal 2*300 # Time $(ns) for NPT and NVT initialization
variable Ni equal v_Nti/v_dt_i # Nr of timesteps for NPT and NVT initialization
variable Nf equal 1000 # 0 # Dump every N steps

## Global parameters
units real
dimension 3
boundary p p p
atom_style full
bond_style morse
special_bonds lj/coul 0.0 0.0 0.5
angle_style harmonic

## Create simulation box
pair_style   lj/cut/tip4p/long 1 2 1 1 0.1546 8.5
read_data    3000_TIP4Pf_medium_density_space_aroundedges.data

## Set up potential
pair_coeff 1 1 ${eps_OO} ${sigma_OO} # Pair style for O-O $(O type 1)
pair_coeff 2 2 ${eps_HH} ${sigma_HH} # Pair style for H-H $(H type 2)
pair_modify mix geometric tail yes # Geometric mixing of pair parameters
pair_coeff 1 2 ${eps_OH} ${sigma_OH} # Pair style for O-H 
kspace_style pppm/tip4p 1.0e-5

## Bond and angle potentials
bond_coeff   1  103.38934  2.287  0.9419
angle_coeff  1  43.95435   107.4

## Create dumps
thermo ${Nf}
thermo_style custom step time temp press etotal ke pe vol density

## Minimize
minimize 1.0e-4 1.0e-6 100 1000
reset_timestep 0

## Set how neighbour lists are updated
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

## Initial velocities
velocity all create ${Temp} ${seed1}
timestep ${dt_i}

# Fix NVT
fix fnvt all nvt temp ${Temp} ${Temp} $(100*dt)

# Compute temp after subracting COM velocity and angular velocity
compute rototemp all temp/rotate
fix temproto all ave/time 100 1 100 c_rototemp mode scalar file temp_roto.dat

run ${Ni}

quit