Dear LAMMPS developers and users!

I repeatedly encounter a problem using Berendsen press-conserving procedure with PETN described by reax/c forcefield. Here is an input file and log for a small test system, the one we are interested in (32x4x4 unit cells) shows the same behaviour. Basically, pressures and volume become NaN right after fix is applied.

Input file:

units real

atom_style charge

variable n_x equal 3

variable n_y equal 1

variable n_z equal 1

variable temperature equal 300

variable temperature_x2 equal ${temperature}*2

log log.petn.reaxc.T${temperature}K.{n_x}x{n_y}x${n_z}.txt

read_data data.reax

replicate 3 1 1

velocity all create ${temperature_x2} 9999

pair_style reax/c NULL

pair_coeff * * ffield.reax 1 2 3 4

timestep 0.1

fix 1 all nve

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

fix 3 all reax/c/bonds 100 1 100 bonds.petn.reaxc.T${temperature}K.{n_x}x{n_y}x${n_z}.txt

dump d1 all cfg 1000 cfg.petn.T${temperature}K.{n_x}x{n_y}x${n_z}.*.cfg id type xs ys zs

dump_modify d1 element C H O N

dump d2 all custom 1000 dump.petn.T${temperature}K.{n_x}x{n_y}x${n_z}.*.dump id type x y z

dump_modify d2 element C H O N

thermo 100

thermo_style custom step temp ke pe pxx pyy pzz etotal vol

run 1000

fix 4 all press/berendsen aniso 0.0 0.0 1000.0

thermo 1

run 11000

unfix 4

thermo 100

compute my_ke all ke/atom

compute my_press all stress/atom

restart 10000 restart_petn_T${temperature}K.{n_x}x{n_y}x${n_z}.*

run 10000