Thermal stability of Polymer Compounds using ReaxFF

Dear All,

I am a beginner with Reactive MD and am working on the thermal stability of polymer compounds composed of polyester and organic additive. The idea is to study the reaction pathways and evolution of reactive species during heating using ReaxFF.

I employ the force field parameters from https://doi.org/10.1021/acs.jpca.7b05257 for PHCO, and I want to validate if it is suitable for my system. The plan is to do a test for a small system, i.e. a polymer chain with an additive molecule. First, I do non-reactive MD with NVT ensemble equilibrated at 300K. Then, I plan to conduct reactive MD simulations heating from 300K to 1000K/1500K/2000K. To validate the parameters, I will compare the cohesive energies from ReaxFF and from DFT using CASTEP, like plotting E/atom against the number of atoms. I am not sure if it is enough for validation. Could anyone help me out?

Besides, I am wondering how can I get the evolution information of species. Should I dump molecule ID?

FYI, my input script has been attached.

#package gpu 4 gpuID 0 3
#units 		real

#atom_style	charge
#read_data	data.mix_0.5

#pair_style      reax/c NULL 
#pair_coeff	* * ffield.reax.PHCO H C O H P 
#newton		on

#fix 		1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

#minimize	1.0e-4 1.0e-6 100 10000
#write_restart	restart.minimizer

#clear

package gpu 4 gpuID 0 3
units 		real

atom_style	charge
read_restart	restart.minimizer

pair_style      reax/c NULL 
pair_coeff	* * ffield.reax.PHCO H C O H P 
newton		on

neighbor	2 bin
neigh_modify	every 10 delay 0 check no

restart		1000000 restart
reset_timestep	0

fix		1 all nve
#fix 		1 all nvt temp 300.0 1000.0 100.0 
#fix		1 all npt temp 1000.0 1000.0 100.0 iso 1.0 1.0 1000.0
fix 		2 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

thermo 		10
thermo_style	custom step temp cpu pe etotal evdwl ecoul vol density

timestep	0.1
dump		1 all atom 10000 dump_mix-0.5pct.lammpstrj
run		10000000

undump		1 

Thank you very much for reading my questions!

Please note that this is more a question about the science you want to do rather than how the software you are using works, and thus this is not a best place/category to ask this kind of question.

You should probably look up in the publication of the specific set of ReaxFF parameters that you are planning to use, what the authors of that publication used to validate their results. And probably also in the publication(s) describing ReaxFF in general. Another place to look would be studies on similar materials.

Molecule IDs in LAMMPS have no structural meaning. They are just a label for atoms so they can be grouped. It is an input parameter and does not change unless you explicitly change them. The reaxff pair style will not make any changes.
To monitor the formation and changes of compounds you can use fix reaxff/species command — LAMMPS documentation and fix reaxff/bonds command — LAMMPS documentation
Please note that these two commands (and the pair style and the QEq fix) were recently renamed when the entire ReaxFF code underwent a significant refactoring. With older versions of LAMMPS those have different names (the new code is backward compatible, i.e. old input files a supposed to continue to work, but the (online) documentation only shows the new names).
In the older LAMMPS versions you have to replace “reaxff” with “reax/c”.

Hello Axel,

Thanks for the hints! I’ll try the two fix commands you suggested. :grinning: