Berendsen barostat and NaN problem for reax/c system

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

Hi Oleg,

fix 4 all press/berendsen aniso 0.0 0.0 1000.0

Step Temp KinEng PotEng Pxx Pyy Pzz TotEng Volume
    1000 406.67685 209.71517 -17215.957 36217.572 27084.292
15370.268 -17006.241 1889.3166
    1001 406.1878 209.46298 -17215.698 28932.519 21402.734
12020.98 -17006.235 2377.8346
    1002 409.32288 211.07968 -15388.542 -nan -nan
-nan -15177.462 -nan
    1003 413.54222 213.25551 -706.71283 -nan -nan
-nan -493.45732 -nan
    1004 413.54222 213.25551 -706.71283 -nan -nan
-nan -493.45732 -nan

I know the problem already appeared in this list, but could not find the

All you need is a larger pdamp value. Please see the doc page on why
berendsen barostat pdamp is related to bulk modulus and why you would
need a larger pdamp for your system.


Ray, thanks for your reply!

I'll try that. As long as needed Pdamp is a function of bulk modulus, it should not depend on system size (except the periodic boundary conditions influence), right?


21.09.2012, 19:01, "Ray Shan" <[email protected]>:


You are welcome. Right, it does not. Default bulk modulus is set to
10, but PETN has a much larger value, hence the increased pdamp. You
can either increase pdamp or modulus - results are the same.


Ok, thx again.


22.09.2012, 03:12, "Ray Shan" <[email protected]>: