Reaxff c/h/o energy conservation and temperature drift in NVE

Dear LAMMPS community,

This question is regarding the Reaxff c/h/o combustion force field (Chenoweth, Van duin, and Goddard, 2008). I am simulating Graphene Oxide (GO), and using the NVE ensemble with NEMD to calculate thermal conductivity. There are issues with energy conservation in the NVE ensemble. For example, one system is GO with 50% oxygen (hydroxyl) concentration. The system is first equilibrated using NPT at 300K for 4ns, and a timestep of 0.5fs. Then it is transitioned from NVT to NVE for 25ps each. When the calculation is done in NVE, and a heat flux is applied (using the M-P algorithm), the temperature drifts from 300K - ~370K in the span of 300ps, as shown in the first attached file. Most likely the equilibration and timestep used are sufficient. This issue is less pronounced when the temperature gradient is removed.

I have tested the same methodology (NEMD) with Silicon using the SW potential, and the temperature does not drift at all, as shown by the second attached file.

Please let me know if you can be of any help,