Pyrolysis of polyacrylonitrile

I am modelling the pyrolysis of polyacrylonitrile and ladder polyacrylonitrile. In the process of pyrolysis, gases usually are evolved like H2, N2, HCN and NH3. However, in my pyrolysis process, only H2 seems to be formed. I wonder if there is a problem with my in.reaxc file that seems to prevent LAMMPS from detecting the other gases significantly.

I am attaching the control file and ffield file too

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 —

ffield.txt (8.5 KB)

control.txt (4.6 KB)

LAMMPS does not care what kind of components you have, it just executes the selected functionality with the given parameters. Whether those parameters are suitable or the choice of simulation conditions is not a LAMMPS problem but a science problem. It is impossible to know from remote and just looking at your input what the thermodynamic and (reaction) kinetic properties of your system are and how high the probability for certain compounds to appear in the given time is. Please recall that atomic scale MD simulations run on tiny systems with very short simulation times compared to experiments, so if higher barriers need to be crossed for certain reactions, you may not see them without running for longer times or using some enhanced sampling method due to the low probability of the reactions to happen.