fix qeq/comb command

Dear:
I am using lammps to simulate ZrO2 materials, and using the fix qeq/comb command in the input file for relexation.
But when I set the Nevery=50 or 100, the tempeture can not reach 300K and vibrate a lot; And when I set the Nevery=500 or 1000, the tempeture can reach 300K, but the tempeture will vibrate every 500 or 1000 timesteps. Is it normal?
I read the inference "# Tao Liang, Tzu-Ray Shan, Yu-Ting Cheng, Bryce D. Devine, Mark Noordhoek, Yangzhong Li, Zhize Lu, Simon R. Phillpot,and Susan B. Sinnott, “Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials”, Materials Science and Engineering: Reports, 74 (2013) 255-279. " It said that “Treating the charge of each ion as a particle in the dunamic charge equilibration approach results in the free energy associated with charges not being conserved, which is analogous to the non-conserved conditon in the energy minimization of ions”. Will it confict when I use the NVE command?

looking forward to your reply
Thank you!

units metal
dimension 3
boundary p p p
atom_style charge

lattice custom 5.147 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 basis 0.5 0.5 0 basis 0.5 0 0.5 basis 0 0.5 0.5 &
basis 0.25 0.25 0.25 basis 0.75 0.25 0.25 basis 0.25 0.75 0.25 basis 0.75 0.75 0.25 &
basis 0.25 0.25 0.75 basis 0.75 0.25 0.75 basis 0.25 0.75 0.75 basis 0.75 0.75 0.75
region box block 0 5 0 5 0 10
create_box 2 box
create_atoms 2 box basis 1 1 basis 2 1 basis 3 1 basis 4 1&
basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2 basis 10 2 basis 11 2 basis 12 2

mass 1 91.224
mass 2 15.999

group zr type 1
compute ch1 zr property/atom q
compute q1 zr reduce ave c_ch1

group o type 2
compute ch2 o property/atom q
compute q2 o reduce ave c_ch2

pair_style comb3 polar_off

#pair_style comb3 polar_off
pair_coeff * * ffield.comb3 Zr O

thermo 100
thermo_style custom step pe ke c_q1 c_q2 temp

min_style sd
minimize 1.0e-8 1.0e-8 100000 100000
velocity all create 500.0 50000 rot yes dist gaussian
fix 1 all nvt temp 300.0 300.0 0.2 drag 0.2

fix zr02 all qeq/comb 1000 0.0001 file fq.out

compute dong all ke/atom
compute shi all pe/atom

reset_timestep 0
timestep 0.001

dump 1 all custom 1000 re.xyz id type x y z c_dong c_shi
run 50000

unfix 1
write_restart re.rs

Increase qeq frequency to every 1 should help with the equilibration. You might also need a stronger thermostat (smaller tdamp).

Ray