So I’ve looked at your earlier data sets and tried a few things.
I had to refresh myself on the details of how pair table works.
Here is what I learned.
The basic issue is that you are trying to tabulate a gigantic range
of energy/force values that vary rapidly in r. One of your tables
goes up to 10^60 at r = 1.0e-6, if I recall. That is hard #1. It’s
compounded by the fact you are doing it at tiny separations. When
LAMMPS converts small r to r^2 (for efficiency), you have even
less coverage of that range with tabulation points. That is hard #2.
I think we could implement a couple options in LAMMPS to help
with this case. Maybe allow for internal spacings in r, rather than r^2.
Or allow preservation of the non-uniform spacings that you have setup
to cover small r better. Both of those would come at a computational cost.
Or we could have an option to use Hermite splines instead of cubic.
I think those will be less oscillatory when eng/force is varying so
rapidly. And wiggles in one spline will not bleed into the others
which is what cubic splines do.
However, here are 3 more simple/concrete things to recognize and try:
a) I believe your tabulated energy and force values are nearly
all inconsistent. I.e. F != -dE/dr. So even if you ran with
a well-splined potential you would have very bad energy
conservation. You should look carefully at how you are generating
b) I see no reason to include tabulated values at tiny separations (in eV)
that are larger than 10^6 or 10^9. Eng = 10^9 would only
be approached if you shot a GeV particle into a head-on collision?
It seems like you would have to use timesteps 1/10^5 smaller
(i.e. sqrt(10^9)) than normal to accurately integrate that interaction? So I can’t
imagine you will do efficient dynamics that sample those
tiny separations of a potential. Excluding them makes
the numerics of the spline fitting a lot easier.
c) For example, attached are several scripts and plots.
table.py - build a table for a LJ potential:
run with python table.py -m LJ -l LJ100 -n 101 -r R 0.3 2.5 > table.txt
to get a tabulation with just 100 points that goes all the way up to ~10^7
in energy and 10^8 in force
table.txt - the tabulation file produced
in.table - simple LAMMPS run that reads the table
and splines it at 1000 pts and uses pair_write to write out the potential
plot.txt = the LAMMPS output potential
plot 1,2,3.pdf = plots of tabulated vs LAMMPS in 3 different R regimes
the results are quite good everywhere, no wiggles
there is essentially no visible difference between tabulated and pair
given that there are 10x more points in pair and it is thus
replacing a straight dotted line with 10 pts on a curve
If I pushed the inner separation to 0.1, then Eng/force ->10^14
and I was unsuccessful avoiding all wiggles, though they were
all in the R regime where force/energy values were humongous.
It’s possible using more points would have helped, but again
I say avoid the issue if at all possible.
Hope these ideas help,
plot.txt (29.1 KB)
in.table (294 Bytes)
plot1.pdf (12.6 KB)
plot2.pdf (13 KB)
plot3.pdf (13 KB)
table.py (3.27 KB)
table.txt (3.91 KB)