Issues with ReaxFF simulation: "Non-numeric pressure" error and Van der Waals warning with Ca

LAMMPS LAMMPS General Discussion
1

Dear all,

I’m a beginner working with ReaxFF reactive force fields. Recently, I came across a 2022 paper by van Duin that provides ReaxFF parameters highly relevant to my system of interest—a hydrous carbonate mixture of (Ca, Na₂, Mg)CO₃ + H₂O. I would like to use the force field from this work to study my system. The original paper and force field files are attached at the end of this post.

However, I’ve encountered a few issues I haven’t been able to resolve, and I’d appreciate any advice from the community:

Problem 1: Simulation becomes unstable under NPT
After running for a while, the simulation crashes with the following error:

ERROR: Non-numeric pressure - simulation unstable (src/fix_nh.cpp:1069)

Following suggestions from forum threads, I have tried:

  • Decreasing the timestep
  • Lowering the temperature and pressure
  • Enlarging the simulation box
  • Adding more steps of energy minimization
  • Performing an initial NVT equilibration before switching to NPT

Unfortunately, the simulation still fails with the same error after some steps under the NPT ensemble.

Problem 2: Van der Waals warning for Ca element
At the beginning of the simulation, I consistently get the following warning:

WARNING: Van der Waals parameters for element CA indicate inner wall+shielding, but earlier atoms indicate a different van der Waals method. This may cause division-by-zero errors. Keeping van der Waals setting for earlier atoms. (src/REAXFF/reaxff_ffield.cpp:251)

I have seen similar warnings mentioned in the forum former posts, but I couldn’t find a definitive solution. I now suspect that this mismatch in van der Waals settings might be causing the instability and ultimately the “Non-numeric pressure” error. If that’s the case, is there a known workaround or correction I should apply?

Simulation details:

  • LAMMPS version: 29 Sep 2021
  • System composition: Na₂CO₃ + CaCO₃ + MgCO₃ + H₂O
  • Initial structure: Built with Packmol
  • Force field: From the 2022 van Duin paper (attached)
  • Input script: Based on the examples/reaxff/H2O input file from the official LAMMPS distribution, with additional energy minimization steps and customized output settings. The full input script and log files are attached.

#----------------- Init Section -----------------
boundary p p p
units real
atom_style full

#----------------- Atom Definition Section -----------------

read_data “in.data”

#----------------- Settings Section -----------------

variable x index 1
variable y index 1
variable z index 1

replicate $x $y $z

pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * reax_ff.carbonate C Ca H Mg Na O
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes

#----------------- Run Section -----------------

timestep 0.25

velocity all create 1 372748 rot yes dist gaussian loop local
thermo 100
thermo_style custom step temp press density vol
thermo_modify flush yes

dump mydump all dcd 100 dump.dcd

fix 1 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 400

min_style cg
minimize 1e-4 1e-6 1000 10000

fix fxnpt all npt temp 1 300 100.0 iso 1000 1000 1000
run 5000
unfix fxnpt

min_style cg
minimize 1e-4 1e-6 1000 10000

fix fxnpt all npt temp 300 500 100.0 iso 1000 1000 1000
run 5000
unfix fxnpt

fix fxnpt all npt temp 500 500 100.0 iso 1000 1000 1000
run 400000

minimize 1e-4 1e-6 1000 10000

write_data steady.data
write_restart steady.restart

Any help or insight would be greatly appreciated. If there’s anything I should clarify or additional files I should provide, please let me know.

Best regards,

Dasgupta et al - 2022 - Development and application of ReaxFF methodology for understanding the chemical dynamics of metal c.pdf (5.6 MB)
in.data (26.8 KB)
reax_ff.carbonate (12.4 KB)
SI.pdf (143.3 KB)
slurm-99561.out (107.8 KB)

2

Warning about inner shielding is because some parameters with nonzero shielding parameters were pasted into a force field without shielding. Reaxff turns off shielding with zero parameters to avoid division by zero.

For the rest see (most important are reaxff Kokkos instead of legacy serial version and check for polarization catastrophe )

3

Thank you very much for your previous suggestions. I followed your advice and used the ReaxFF with KOKKOS acceleration. Initially, the simulation ran smoothly, but after around 69,000 steps, I encountered the following error again:

ERROR: Non-numeric pressure - simulation unstable (src/fix_nh.cpp:1069)

I will attach the relevant input files once more for your reference.

In addition, I checked for the possibility of a polarization catastrophe as you suggested. Specifically, I examined the condition η > 7.2 × γ for QEq. In my force field, the Ca atom has an etaEEM value of 6.3136 and a gammaEEM value of 0.9605, which results in:

6.3136 < 7.2 × 0.9605 = 6.9156

This seems to violate the stability criterion you mentioned. Could this be a potential cause of the simulation instability and the non-numeric pressure error?

I truly appreciate your help and would be grateful for any further advice you could provide.
log.lammps (90.0 KB)
system.in (1.2 KB)

4

Yes. Your qeq charge equilibration matrix is singular so in the fix qeq iterative solver it divides by 0 if it was in exact arithmetic but since this calculation is being done in approximate floating point precision (ie. double) it divides by a very small number close to 0 which gives a very big number close to infinity then BOOM game over at next iteration.

I’m surprised the original article didn’t have same problem. Maybe your system/geometry is more sensitive to polarization compared to theirs.