Pair distance > table outer cutoff for Hertzian-cut potential

Respected LAMMPS users and developer,
I am using lammps-15Jun2023 to simulate a hertzian-cut potential with table potential, just like harmonic-cut potential. I have to guarantee

  1. a very high accuracy to achieve mechanical equilibrium (so I increase the number to 3000000 and and have proven it with pair_write where the deviation from the input emerges only approaching for extremely close 0)
  2. a zero force when there is no contact for neighbours. However, I always meet the error just like “Pair distance > table outer cutoff: ijtype 1 1 dist 1”. As I understand, I have set the outer boundary to be the sum of radius of neighbors but fail to avoid this error (even though I add a new line of “100000 1.01 0 0” in the last).
    I have tried my best, but still failed. Thank you for any generous help!!!

My code is like

pair_style table spline 3000000
pair_coeff 1 1 R11.table Hertzian_R11 1
pair_coeff 1 2 R12.table Hertzian_R12 1.2
pair_coeff 2 2 R22.table Hertzian_R22 1.4

The input of R11.table is like:

Hertzian_R11
N 100000

1 1e-6 38811.861989976234 98009.85149988899
2 0.0002 38785.07070938662 97969.25300356193

99999 0.99999999 8.057611055120241e-35 6.048042357555986e-19
100000 1 0 0

Why don’t you just extend the tables beyond the selected cutoff? or choose the cutoff just a tiny bit shorter so the ambiguity at the cutoff has no impact.

If you need very high precision with a custom potential, you may want to have a look at pair style lepton. There you can enter a custom expression and it is evaluated analytically (including the derivative for the force, which is computed symbolically).

1 Like

No, both do not work. I tried to use

pair_coeff 1 1 R11.table Hertzian_R11 0.999

or

N 100001

100001 1.01 0 0 (even for 100001 5 0 0)

but with the same error. “Pair distance > table outer cutoff: ijtype 1 1 dist 1”
As for the lepton, it is very convenient to use, but still with new error and lower efficiency than 2.5~5 times of using table potential. The code is like:

####### used the low-precision table potential to pre-equilibrate with the code shown before (this part runs well)
minimize 0.0 1e-11 100000 10000000000
write_data bidispersed_01.data
######## using lepton to optimize the mechanical equilibrium
pair_style lepton/sphere 3
pair_coeff 1 1 "step(r0-r)*0.4*sqrt(riTj/r0)*(r0-r)^2.5*2e5; r0=0.5+0.5; riTj=0.5*0.5"
pair_coeff 1 2 "step(r0-r)*0.4*sqrt(riTj/r0)*(r0-r)^2.5*2e5; r0=0.5+0.7; riTj=0.5*0.7"
pair_coeff 2 2 "step(r0-r)*0.4*sqrt(riTj/r0)*(r0-r)^2.5*2e5; r0=0.7+0.7; riTj=0.7*0.7"
minimize 0.0 1e-11 100000 10000000000
write_data bidispersed_02.data

,but with nan in pressure and Fnorm and in the dump.atom file once I use lepton:

Step Temp TotEng Press Fnorm Fmax
1493 0 nan -nan -nan 0.81983669
2000 3.2276003e-07 0.00072780886 -nan 546.94971 18.995159
3000 6.2803188e-07 0.00056606878 -nan 464.19947 17.950523
4000 -nan -nan -nan 394.56802 14.247259
5000 1.2008342e-07 nan -nan -nan 13.905959

74000 3.075035e-27 7.7645077e-27 -nan 1.0049589e-11 2.2633653e-13
74375 -nan -nan -nan 9.9879596e-12 1.625703e-13
Loop time of 236.439 on 96 procs for 73029 steps with 88776 atoms
80.5% CPU use with 96 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = force tolerance
Energy initial, next-to-last, final =
nan nan 4.62486220433769e-27
Force two-norm initial, final = -nan 9.9879596e-12
Force max component initial, final = 113.94563 1.625703e-13
Final line search alpha, max atom move = 0 0
Iterations, force evaluations = 73029 79239
Generated 0 of 10 mixed pair_coeff terms from geometric mixing rule
ERROR on proc 0: Non-numeric atom coords - simulation unstable (…/domain_omp.cpp:58)
Last command: write_data bidispersed_0.9996_0.data

The dump.file writes:

ITEM: TIMESTEP
1500
ITEM: NUMBER OF ATOMS
88776
ITEM: BOX BOUNDS pp pp pp
0.0000000000000000e+00 3.5000000000000000e+02
0.0000000000000000e+00 3.5000000000000000e+02
-1.0000000000000001e-01 1.0000000000000001e-01
ITEM: ATOMS id type x y vx vy fx fy radius
1 1 175 175 0 0 0 0 0.5
2 1 -nan -nan 0 0 0 0 0.5
3 1 34.46526057685814 83.16471202871583 -0.005441 -0.003598 -1.281175656 -0.8701474842 0.5
4 1 -nan -nan 0 0 0 0 0.5

I will try to optimize the configuration first to avoid “nan” , but can do nothing about the efficiency. I checked the file with pair_write, which indeed ouputs the high-accuracy interaction between all types (remarkable and bravo!) except that the last line to be 0 or nan in the last line.

1000 1 nan -nan

I record these for future potential readers.

Best regards to Akohlmey,
Yang

It must work. There must be something wrong in your tables, that is not evident from remote. Can you prepare a minimal runnable input deck (ideally with much fewer table point) that reproduces the error you see and upload it here? or upload it elsewhere and provide the link?
I am travelling for a bit, but will try to look into it as soon as I find time.

It is not meant to be efficient, just a way to get something running and alternate option to generate tables with pair_write. If you get NaN already for the energy, there is something not quite right with your expression.

1 Like

R11.table (63.8 KB)
R12.table (63.6 KB)
R15.table (63.6 KB)
R22.table (63.6 KB)
R25.table (63.6 KB)

bidispersed_table.lj (4.4 KB)
table_error.out (265.2 KB)

This error sometimes do not happen but with higher probability for long-time running or performing energy minimization for many times. I have tried to make the original code as short as possible, and I hope this code is easy to understand.

Really? I am using Hertzian-cut potential what LAMMPS uses but with the expression. The expression could output the exact energy and force consistent with the theoretical formulation when I use pair_write to check.
bidispersed_Lepton.lj (5.1 KB)
lepton_error.out (40.2 KB)

Have fun in traveling!

After some digging through the source code, I have found that currently your table must go beyond the cutoff used. This usually doesn’t get triggered easily since people use a moderate number of internal tabulation points. Technically, this behavior would not be needed for linear interpolation because the interpolation does not need to access data beyond the current table entry, but the check in the code for linear and spline (which does access the next set of table data) is the same.

Taking your R11.table file and editing the header to:

Hertzian_S_0.9994_R11
N    1003

1 1.0000010000287557e-06 39927.94331129907 99879.88608988757

and then changing the end of the table to:

1001 0.9994 0 0
1002 0.9994001 0 0
1003 0.9994002 0 0

should avoid the “outer cutoff” error. Similar edits would apply to the other tables.

As for the Lepton expression. The problem is here (r0-r)^2.5. This turns to NaN for any r > r0. It needs to be replaced with abs(r0-r)^2.5 and the calculation can run without the NaN errors.

1 Like

Thank you so much for your help so far! Did you run it successfully?I have done what you suggested but still failed whether I use spline or linear interpolation.
By the way, the table I gave you is to use non-uniform separation to create the table to put more focus on the slight contact between neighbors. Before asking for help, I tried to create the table with 1e6 lines of uniform separation and by this way I got less chances of error but with unwanted attractive interaction within the range of cutoff and poor accuracy.

Hertzian_S_0.9994_R11
N    100002 R 1e-6 0.9994
...
100001 0.9994 0 0
100002 0.9994001 0 0

Yes, you are right. I succeed in running it this time but I need to run ~4e8 steps to achieve mechanical equilibrium of force on every disk < 1e-7 , which will take me too much time for me to use lepton to run.