#================================ Simulation Setting ================================

units           real
boundary        p p p   #s means shrink-wrapped
atom_style      charge

read_restart    restart2.reaxff

pair_style      reax/c NULL checkqeq yes safezone 1.5 mincap 100
pair_coeff      * * NM.ff C H N O
fix             fqeq all qeq/reax 1 0.0 10.0 1e-6 reax/c

compute         reax all pair reax/c
variable        eb      equal c_reax[1]         # bond energy
variable        ea      equal c_reax[2]         # atom energy
variable        elp     equal c_reax[3]         # lone-pair energy
variable        emol    equal c_reax[4]         # molecule energy (always 0.0)
variable        ev      equal c_reax[5]         # valence angle energy
variable        epen    equal c_reax[6]         # double-bond valence angle penalty
variable        ecoa    equal c_reax[7]         # valence angle conjugation energy
variable        ehb     equal c_reax[8]         # hydrogen bond energy
variable        et      equal c_reax[9]         # torsion energy
variable        eco     equal c_reax[10]        # conjugation energy
variable        ew      equal c_reax[11]        # van der Waals energy
variable        ep      equal c_reax[12]        # Coulomb energy
variable        efi     equal c_reax[13]        # electric field energy (always 0.0)
variable        eqeq    equal c_reax[14]        # charge equilibration energy

fix             species all reax/c/species 1 1 1000 out.species3 element C H N O

fix             fnveB31 all nve

thermo          1000
thermo_style    custom step cpuremain temp epair etotal enthalpy press &
                             v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa &
                             v_ehb v_et v_eco v_ew v_ep v_efi v_eqeq
							 
rerun data.d1064_2500K.lammpstrj dump x y z q vx vy vz format native

unfix           fnveB31
unfix			species
unfix			fqeq

print			"all done"