Problem of reax/c with NVE and NVT

Dear All,
      I am trying to simulate the graphene with reax/c implemented in
Lammps. But I have encountered some difficulties. Two steps of
simulations are wanted to be carried: NVT, then NVT. But I find that
the temperature is rising unreasonably and the size of system is
expanding unlimitedly. The input and output files are in the
attachment. Anyone can give some suggestions?

Thanks in advance!
                                                                Wang Jian

The input script is as follows

units real
dimension 3
newton on
boundary s p s
atom_style charge
neighbor 0.3 nsq
neigh_modify check yes
read_data graph_reaxc.xyz
velocity all create 500 1212127641 mom yes rot yes dist
gaussian units box

pair_style reax/c NULL checkqeq no
pair_coeff * * ffield.reax.cho 1

fix NVT all nvt temp 500 500 1
thermo_style custom step temp etotal vol
thermo_modify lost warn
thermo 10000

# Run
timestep 0.5
dump data1 all custom 10000 outputnvt.xyz type x y z
run 100000

unfix NVT
fix nve all nve
dump data2 all custom 10000 outputnve.xyz type x y z
run 100000

The structure data is as follows.

2360 atoms
     1 atom types

     -11.541138 13.053984 xlo xhi
    -125.433333 125.906667 ylo yhi
      -1.675000 1.675000 zlo zhi

Masses

         1 12

Atoms

     1 1 0.000000 2.459512 -124.960000 0.000000
     2 1 0.000000 4.919024 -124.960000 0.000000
     3 1 0.000000 7.378536 -124.960000 0.000000
     4 1 0.000000 9.838049 -124.960000 0.000000
     5 1 0.000000 12.297561 -124.960000 0.000000
     6 1 0.000000 -9.838049 -124.960000 0.000000
     7 1 0.000000 -7.378536 -124.960000 0.000000
     8 1 0.000000 -4.919024 -124.960000 0.000000

Part of log file is

LAMMPS (17 Oct 2012)
……
memory usage per processor = 1.90524 Mbytes
Step Temp TotEng Volume
       0 500 -422493.8 3.9355046
   10000 496.96668 -447001.31 34366.804
   20000 511.68336 -446813.17 52925.553
   30000 489.01261 -447035.91 52192.758
   40000 494.57279 -447055.84 49979.958
   50000 504.45062 -446935.76 52020.052
   60000 516.99174 -446829.53 55296.869
   70000 511.90547 -446981.81 60563.274
   80000 493.29511 -447077.87 65691.772
   90000 497.178 -446963.27 69837.58
  100000 494.79137 -447110.48 67269.375
Loop time of 2077.01 on 4 procs for 100000 steps with 2360 atoms

Pair time (\) = 1987\.14 \(95\.6728\) Neigh time \() = 12.6195 (0.607578)
Comm time (\) = 58\.2898 \(2\.80642\) Outpt time \() = 0.036369 (0.00175102)
Other time (%) = 18.9314 (0.911474)

Nlocal: 590 ave 595 max 584 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Nghost: 196.25 ave 200 max 190 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Neighs: 32088.8 ave 32362 max 31777 min
Histogram: 1 0 0 0 0 2 0 0 0 1

Total # of neighbors = 128355
Ave neighs/atom = 54.3877
Neighbor list builds = 6838
Dangerous builds = 53

unfix NVT

fix nve all nve
dump data2 all custom 10000 outputnve.xyz type x y z

run 100000
Memory usage per processor = 1.90524 Mbytes
Step Temp TotEng Volume
  100000 494.79137 -447110.48 67647.011
  110000 1792.6493 -427406.65 63336.634
  120000 19927.294 -13282.889 29380787
  130000 20029.52 82340.806 7.0560577e+08
  140000 20103.372 89890.917 2.5201139e+09
  150000 20106.79 92315.111 5.4685969e+09
  160000 20113.849 93225.42 9.5450081e+09
  170000 20136.624 93850.547 1.4749347e+10
  180000 20113.54 94312.333 2.1081615e+1

in_script.zip (360 KB)

three comments:

a) you have dangerous neighborlist builds. that is never a good sign.
    your neighborlist skin of 0.3 is very small for a system with real units.
    why did you change the default? there is no reason to?

b) have you tested using a smaller time step?

c) please keep in mind that reax is not an ab initio calculation, so
  you cannot just use any parameter set for any kind of molecule.
  more specifically, the one you are using was - to the best of my
  knowledge - not parameterized for graphene (rather for small
  molecules in the gas phase at high temperature), and thus is
  likely to produce bogus results.

axel.