I am trying to use LAMMPS to simulate a displacement cascade in a uranium dioxide lattice. My difficulty lies in implementing a particular potential for LAMMPS to use. My potential consists of a ZBL term for short-range interactions and a coulombic + Buckingham term for intermediate-range interactions. I have splined these two regions together to create a smooth potential for each of the three interacting atom pairs (U-U, O-O, O-U). For each atom pair, I have an analytic expression for this potential.

For input into LAMMPS, I thought that the best way to implement this potential was to use pair_style hybrid/overlay to superpose (1) pair_style coul/long and (2) a table from my above splined potential from which I have subtracted the coulombic term. My motivation was to recreate my original spline and to preserve long-range coulombic interactions. For each atom pair, I input a separate table. I used 11,000 points for these input tables and I used the spline option for interpolation. One noteworthy issue is that the repulsive ZBL potential leads to a very high gradient for values in my input tables that are close to 0. Because of this, I used 10,000 points for this core ‘high gradient’ region.

I implemented this in LAMMPS and I used the pair_write command to output a table containing the interpolated values of the potential created by the hybrid/overlay described above. The graph of this LAMMPS-generated potential has huge osciallations in the region corresponding roughly to the ZBL core.

To trouble-shoot, I examined each part of the hybrid/overlay potential separately. I discovered 2 Problems.

First, I used pair_write to examine LAMMPS’ interpolations of coul/long. That is, I defined all interactions to be only via coul/long (without inputting my tables) and then I output the interpolated potentials using pair_write. I discovered that the potentials for all three atom pairs (U-U, O-O, O-U) were identical. The tables for each were the same despite my specifying different charges for uranium and oxygen. I see no other possible explanation for this.

Then, I defined all interactions to be strictly according to my input tables (without any coulombic interaction). All 3 tables output by pair_write contain huge oscillations near the core region and do not look anything like a splined interpolation. However, the output tables are reasonably close to my input tables at longer distances. I played with the number and spacing of my table points but this did not change anything. When I changed the interpolating option from spline to linear, it removed the oscillations at short ranges; however, instead, at short ranges, the potential does not decrease exponentially but rather increases in the approximate ZBL region before sharply decreasing to more realistic values at greater distances.

I would like to know

- why LAMMPS treats the Coulombic potential of all three atom pairs the same;
- how to submit a table that LAMMPS can interpolate correctly. Of course, if there is a better way to approach my problem, then I would be happy to hear this.

Thank you for your attention.