Simulation Crashes Randomly with 'Non-numeric Pressure' ERROR

[Antoine]

Hello everyone,
I hope you’re doing well. I’m encountering a problem where my simulation randomly stops with the following error:
ERROR: Non-numeric pressure

Interestingly, I’m able to continue the simulation using the last restart file, and I’ve checked—there’s no atom overlap issue. The error occurs unpredictably: sometimes after 1 ns, sometimes earlier or later.

I came across this related topic:

Very rare non-numeric pressure due to atom overlaps?

Since you are using hybrid/overlay, you could just not have the lj/cut pair style be applied to those pairs. Yes, because you need to divide by this (squared) distance. No they are included, if they are listed under pair_coeff. Yes, with hybrid/overlay you can just not define pair_coeff for lj/cut for those pairs. With regular hybrid, you can use pair_coeff none and with plain lj/cut you can set the cutoff to 0 or you could use neigh_modify exclude.

I initially used a null coefficient line, but even after omitting it, the error persists.

I’m using LAMMPS version: 29 Aug 2024 - Update 1.
Here is my ff after modification

pair_style hybrid/overlay lj/cut 12.0 coul/long 12.0
pair_modify tail yes

# pair_coeff * * lj/cut 0.0 0.0

pair_coeff 2 4 lj/cut 5.95E-04 3.2607
pair_coeff 2 3 lj/cut 5.60E-04 3.2696
pair_coeff 2 6 lj/cut 5.95E-04 3.2607
pair_coeff 2 5 lj/cut 5.30E-04 3.2338
pair_coeff 4 4 lj/cut 6.18E-02 3.4478
pair_coeff 3 4 lj/cut 4.54E-02 3.626
pair_coeff 4 6 lj/cut 6.18E-02 3.4478
pair_coeff 1 4 lj/cut 7.06E-04 5.4523
pair_coeff 4 9 lj/cut 1.04E-03 4.4456
pair_coeff 4 5 lj/cut 8.72E-01 2.8954
pair_coeff 3 3 lj/cut 1.24E+00 2.7351
pair_coeff 3 6 lj/cut 4.54E-02 3.626
pair_coeff 1 3 lj/cut 8.69E-04 4.3654
pair_coeff 3 9 lj/cut 1.46E-03 4.3654
pair_coeff 3 5 lj/cut 5.26E-03 4.2407
pair_coeff 6 6 lj/cut 6.18E-02 3.4478
pair_coeff 1 6 lj/cut 8.69E-04 4.3654
pair_coeff 6 9 lj/cut 1.04E-03 4.4456
pair_coeff 5 6 lj/cut 8.72E-01 2.8954
pair_coeff 1 5 lj/cut 8.76E-04 4.3654
pair_coeff 5 9 lj/cut 6.04E-04 4.4723
pair_coeff 5 5 lj/cut 1.54E-01 3.1627

bond_style harmonic
bond_coeff 1 554.1349 1.0000
bond_coeff 2 554.1349 1.0000

angle_style harmonic
angle_coeff 1 45.7696 109.47

pair_coeff * * coul/long
kspace_style ewald 1e-5
special_bonds coul 0 0 0 lj 1 1 1

for input file :

processors * * *

timer off
units real
boundary p p p
atom_style full

timestep 1.0

group water type 5 7

fix wshake water shake 0.0001 20 0 b 2 a 1

fix NPT all npt temp 300.0 300.0 100.0 x 1.0 1.0 1000.0 y 1.0 1.0 1000.0 z 1.0 1.0 1000.0 couple none

Thanks in advance.

[akohlmey]

@Antoine

pair_style hybrid/overlay lj/cut 12.0 coul/long 12.0

Why use pair style hybrid/overlay and not lj/cut/coul/long?

@Antoine

special_bonds coul 0 0 0 lj 1 1 1

Why this setting?

@Antoine

kspace_style ewald 1e-5

Why Ewald and not PPPM?

You have atom types without LJ parameters. Have you checked that those atoms cannot come close to another due to their charges?

[Antoine]

Thank you it work :)) , i think the problem was the special_bonds coul 0 0 0 lj 1 1 1 command

[akohlmey]

@Antoine

i think the problem was the special_bonds coul 0 0 0 lj 1 1 1 command

That doesn’t answer my question(s).

Since I provided help with your problem, please also have the courtesy to answer my questions.

[Antoine]

i’m pretty sure i took it directly from here (the layout but not the same parameters)
also when changing ewald to ppm , the simulation become Way much faster

https://doi.org/10.1063/1.5042755

SupplMat.pdf (223.3 KB)

[akohlmey]

The main issue that I see is that the authors of the paper used wolf summation for coulomb instead of a long-range solver that you are using. Those can be quite different.

It is rather strange that they use the special_bonds command that they do, but that is part of the force field, obviously.