# REAX potential for PPAN system
# .....

units		real

atom_style	charge
atom_modify	map yes
read_data	datafile 

pair_style	reaxff/kk lmp_control
pair_coeff	* * ffield C N H

#Neighbor list
neighbor	2.5 bin
neigh_modify	delay 0 every 5 check no

#Emin
# Variable for minimization
variable        etol equal 1.0e-10
variable        ftol equal 1.0e-10
variable        maxiter equal 50000
variable        maxeval equal 10000
variable        dmax equal 1.0e-1

# Energy Minimization and equilibration
fix             1 all qeq/reaxff/kk 1 0.0 10.0 1.0e-6 reaxff maxiter 1000


minimize        ${etol} ${ftol} ${maxiter} ${maxeval}
min_style       cg

dump            1 all atom 1 dump.reax1
dump_modify     1 scale no
run 		0
undump		1

#Initial temp
velocity all create 300 8732958
reset_timestep	0
restart         1000 refile1 refile2

#Property computes
compute e all pair reaxff
variable eb equal c_e[1]
variable ea equal c_e[2]
variable elp equal c_e[3]
variable emol equal c_e[4]
variable ev equal c_e[5]
variable epen equal c_e[6]
variable ecoa equal c_e[7]
variable ehb equal c_e[8]
variable et equal c_e[9]
variable eco equal c_e[10]
variable ew equal c_e[11]
variable ep equal c_e[12]
variable efi equal c_e[13]
variable eqeq equal c_e[14]
variable stp equal "step"
variable Pe equal "pe"

variable Ke equal "ke"
variable Te equal "etotal"
variable Temp equal "temp"
variable P equal "press"
variable Px equal "pxx"
variable Py equal "pyy"
variable Pz equal "pzz"
variable Lx equal "lx"
variable Ly equal "ly"
variable Lz equal "lz"
variable Dens equal "density"

#NVT MD Nose-Hoover thermostat
fix             5 all reaxff/bonds/kk 1000 bonds.out
fix             6 all reaxff/species/kk 1 10 1000 mol.out position 1000 mol.pos

fix             7a all print 1000 "${stp} ${Pe}  ${eb}  ${ea}  ${elp}  ${emol}  ${ev}  ${epen}  ${ecoa}  ${ehb}  ${et}  ${eco}  ${ew}  ${ep}  ${efi}  ${eqeq}" file energy_comp
fix             7b all print 1000 "${stp} ${Pe}  ${Ke}  ${Te}  ${Temp}  ${P}  ${Px}  ${Py}  ${Pz}  ${Lx}  ${Ly}  ${Lz}  ${Dens}" file Thermo


thermo 		100
thermo_style	custom step time pe ke etotal vol temp press pxx pyy pzz lx ly lz density

dump            1 all atom 1000 dump.reax_md
dump_modify	1 scale no

# ——————————————————————————————————————
#  Pyrolysis workflow for one target temperature
# ——————————————————————————————————————

# 1) Settings
variable    dt        equal "0.25e-3"        # timestep in ps (0.25 fs)
timestep    0.25

# target temperature (K) and ramp rate (K/ps)
variable    Ttarget   equal "2800"
variable    ramprate  equal "50"             # K/ps

# compute number of steps for heating and cooling
variable    ramp_time_ps  equal "(v_Ttarget - 300.0)/v_ramprate"
variable    ramp_steps    equal "v_ramp_time_ps/v_dt"

# hold at target T for 2 ns (2000 ps)
variable    hold_ps      equal "2000"
variable    hold_steps   equal "v_hold_ps/v_dt"

# cool back to 300 K at the same rate
variable    cool_time_ps equal "(v_Ttarget - 300.0)/v_ramprate"
variable    cool_steps   equal "v_cool_time_ps/v_dt"

# final NVT equilibration at 300 K for 500 ps
variable    equil_ps     equal "400"
variable    equil_steps  equal "v_equil_ps/v_dt"

variable    image        equal "v_ramp_steps + v_hold_steps + v_cool_steps + v_equil_steps"

# 3) NVT heat ramp from 300 → ${Ttarget}
fix             heat all nvt temp 300.0 ${Ttarget} 100.0
run             ${ramp_steps}
write_restart   refile.heat
write_data      data.heat
unfix	        heat

# 4) Hold at high T in NVT
fix             hold all nvt temp ${Ttarget} ${Ttarget} 100
run             ${hold_steps}
write_restart   refile.hold
write_data      data.hold
unfix           hold

# 5) NPT cooling back to 300 K, 1 atm
fix             cool all npt temp ${Ttarget} 300.0 100.0 iso 1.0 1.0 1000.0
run             ${cool_steps}
write_restart   refile.cool
write_data      data.cool
unfix           cool

# 6) Freeze volume and get stress-free snapshots
fix             eq all nvt temp 300.0 300.0 100.0
run             ${equil_steps}


run             80000 start 0 stop ${image}
write_restart   refile.*
write_data      data.*

run             80000 start 0 stop ${image}
write_restart   refile.*
write_data      data.*

run             80000 start 0 stop ${image}
write_restart   refile.*
write_data      data.*

run             80000 start 0 stop ${image}
write_restart   refile.*
write_data      data.*

run             80000 start 0 stop ${image}
write_restart   refile.*
write_data      data.*


# — Repeat blocks 3–7 for each Ttarget = 2200, 2500, 2800 —