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 1replicate $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 yesdump 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 10000fix fxnpt all npt temp 1 300 100.0 iso 1000 1000 1000
run 5000
unfix fxnptmin_style cg
minimize 1e-4 1e-6 1000 10000fix fxnpt all npt temp 300 500 100.0 iso 1000 1000 1000
run 5000
unfix fxnptfix fxnpt all npt temp 500 500 100.0 iso 1000 1000 1000
run 400000minimize 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)