# LAMMPS input file for MPD-TMA crosslinking simulation
# Based on the provided molecular structures


# --- INITIALIZATION ---
units           real
boundary        p p p
atom_style      full
bond_style       harmonic
angle_style      harmonic
dihedral_style   opls
improper_style   cvff
pair_style       lj/cut/coul/long 12.0
kspace_style    pppm 1.0e-4
special_bonds   lj/coul 0.0 0.0 0.5


# --- READ DATA FILE ---
# Assuming you have prepared a data file with MPD and TMA molecules
read_data       MPD_TMA_crosslinked_FULL.data extra/bond/per/atom 1 

# Set up groups for TMC and MPD molecules

group CC1   type 1 # 1  CC1
group CM1   type 2 # 2  CM1
group CM2   type 3 # 3  CM2
group CM3   type 4 # 4  CM3
group CT1   type 5 # 5  CT1
group CT2   type 6 # 6  CT2
group HA    type 7 # 7  HA
group HM    type 8 # 8  HM
group HO    type 9 # 9  HO
group HT    type 10 # 10  HT
group NA1   type 11 # 11  NA1
group OC1   type 12 # 12  OC1
group OH    type 13 # 13  OH

# Charges 

set type 1 charge  0.6797 # CC1
set type 2 charge  0.1716 # CM1
set type 3 charge  -0.230 # CM2
set type 4 charge  -0.1565 # CM3
set type 5 charge  -0.1956 # CT1
set type 6 charge  0.014 # CT2
set type 7 charge  -0.23105  # HA
set type 8 charge  0.10745 # HM
set type 9 charge  0.4315 # H0
set type 10 charge 0.1816 # HT
set type 11 charge -0.4621 # NA1  
set type 12 charge -0.5761 # # OC1  
set type 13 charge -0.5351 # OH 


# --- MOLECULE TEMPLATES ---
# Load molecule templates (make sure these files have correct type numbering)
molecule pre_rxn pre_reacted_mpd_tma.data
molecule post_rxn post_reacted_mpd_tma.data


# Set up fix bond/react
fix rxn all bond/react  stabilization yes statted_grp 0.03 react rxn1 all 1 0 4.5 pre_rxn post_rxn mpd_tma_map.txt prob 1.0 12345 


# --- SETTINGS ---
neighbor        2.0 bin
neigh_modify    every 1 delay 0 check yes

# --- OUTPUT SETTINGS ---
thermo_style    custom step temp press pe ke etotal density vol
thermo          1000

# Dump configurations
dump 1 all atom 1000 dump.lammpstrj

# --- SIMULATION ---
timestep        1.0  # 1 fs

# Initial minimization
minimize 1.0e-4 1.0e-6 1000 10000

# Short NVT equilibration
velocity all create 300.0 12345
fix nvt all nvt temp 300.0 300.0 100.0
run 10000
unfix nvt

# Production run with crosslinking
fix npt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run 100000