# LAMMPS Input Script for Water Droplet on Anatase Surface

# Initialize simulation
units           real
boundary        p p m
atom_style      full
bond_style      harmonic
angle_style     harmonic
special_bonds   lj/coul 0.0 0.0 0.5

# Read anatase data file
read_data       anatase.data extra/atom/types 2 extra/bond/types 1 extra/angle/types 1 extra/bond/per/atom 4 extra/angle/per/atom 4 extra/special/per/atom 4

# Adjust box size in z-direction
change_box      all z final 0.0 150.0

# Define masses for water atoms (SPC/E model)
mass            3 15.9994  # Oxygen
mass            4 1.008    # Hydrogen

# Load water molecule template
molecule        water spce.mol

# Define spherical region for water placement
region droplet sphere 50.0 50.0 100.0 20.0 units box  # Center at (50, 50, 100), radius = 20 Å

# Create water molecules in the liquid region
create_atoms 0 random 512 23462134 droplet mol water 85342

# Define groups
group           H2O type 3 4
group           slab type 1 2

# Bond coefficients for SPC/E water
bond_coeff      1 450.0 1.0   # Bond type 1: Oxygen-Hydrogen bond

# Angle coefficients for SPC/E water
angle_coeff     1 55.0 109.47  # Angle type 1: H-O-H angle

# Pair coefficients
pair_style      lj/charmm/coul/charmm 10 12
pair_coeff      1 1 0.0 0.0   # Ti-Ti (anatase)
pair_coeff      1 2 0.0 0.0   # Ti-O (anatase)
pair_coeff      2 2 0.0 0.0   # O-O (anatase)
pair_coeff      3 3 0.1553 3.166   # O-O (water-water, SPC/E model)
pair_coeff      4 4 0.0 0.0       # H-H (water-water, no LJ interaction)
pair_coeff      1 3 1.50 8.0     # Ti-O (anatase-water)
pair_coeff      2 3 1.50 8.2    # O-O (anatase-water)
pair_coeff      1 4 2.0 7.2       # Ti-H (anatase-water)
pair_coeff      2 4 2.0 7.2       # O-H (anatase-water)

# Enable geometric mixing rule for cross-interactions
pair_modify     mix geometric

# Neighbor list settings (optimize for performance)
neighbor        10.0 bin
neigh_modify    every 1 delay 0 check yes

# Initialize velocities
velocity        H2O create 100.0 465135

# Output settings
reset_timestep  0
thermo          1000
thermo_style    custom step temp press vol density

# SHAKE constraints for water molecules
fix myshk H2O shake 1.0e-6 100 10000 b 1 a 1

# Minimize the system to remove high energies
minimize 1.0e-6 1.0e-8 10000 100000

# Thermostat for the entire system
fix mynvt all nvt temp 100.0 300.0 200.0

# Fix the slab atoms in place (prevent movement)
fix freeze_slab slab setforce 0.0 0.0 0.0

# Dump trajectory
dump            mydmp all atom 1000 dump.lammpstrj

# Set timestep
timestep 0.25

# Equilibrate the system
run 5000

# Run simulation
run 10000

# Write final data file
write_data      data.anatase_water pair ij