Combining table potential with hybrid buck/coul/long/cs + zbl

I am modeling a Li-Al-O system with a core shell model for Oxygen and buck/coul potential for all the 3 species. I have added zbl for all the core particles without adding zbl interactions for the shell of oxygen. My potential description looks like this:
units metal
dimension 3
boundary p p p
atom_style full
#atom_modify map hash

fix csinfo all property/atom i_CSID
read_data smallgrain_stoich.lmp fix csinfo NULL CS-Info

group cores type 1 2 3 # Li core Al core O core
group shells type 4 # O shell

pair_style hybrid/overlay buck/coul/long/cs 10 12 zbl 0.5 2
pair_coeff * * buck/coul/long/cs 0.0 1.000 0.00
pair_coeff 1 4 buck/coul/long/cs 632.1018 0.2906 0.00 #12 #Li core O shell
pair_coeff 2 4 buck/coul/long/cs 1109.92381 0.31540 0.00 #10 #Al core O shell
pair_coeff 4 4 buck/coul/long/cs 12420.5 0.2215 29.07 #25 #O shell O shell

bond_style harmonic
bond_coeff 1 31.0 0.0
kspace_style ewald 2.0e-3

pair_coeff 1 2 zbl 3 13

pair_coeff 1 3 zbl 3 8

pair_coeff 2 3 zbl 13 8

pair_coeff 1 1 zbl 3 3
pair_coeff 2 2 zbl 13 13
pair_coeff 3 3 zbl 8 8

Now, using this potential I was trying to model grain boundaries and simulate cascade displacements caused by a PKA but my system is experiencing 10^19 K temperatures and ends in oxygen core shell separation error. To fix the situation, I am trying to modify the potential by changing the oxygen shell-shell (pair_coeff 4 4 ) interaction into a splined potential using the table style. I want to do the same for Li core-O shell and Al core- O shell interactions but keep all other interactions as they are in the present form. Can you kindly suggest me if this is advisable? If so, could you guide me on how to implement the table potential along with buck/coul/long/cs + ZBL?

What makes you believe that this would be a meaningful correction? Which implies that you have already identified the cause for the unexpected behavior. However, your post implies that this is not the case.

You have a very complex setup. Can this run non-PKA simulations without issues and with conserving energy sufficiently well? Why use such ultra-low accuracy for kspace convergence?

Are you sure the problem is not caused by using too large a timestep?

Yes this potential runs non pka simulations fine in a single crystal. As soon as I introduce grain boundaries, and run a normal equilibration, the bond atoms separate in oxygen core-shell. I was thinking of freezing the oxygen core-shell to avoid the separation. The time step is already extremely small and limited by adaptive Dt variable timestep function as well with D_max=0.03 Angstrom.

This sounds like there is a fundamental problem with your force field parameters/settings or something related, that only doesn’t show for a perfect crystal due to symmetry.
The fact that this happens with a very small timestep, makes this even more suspicious.

You need to try different minimal subsystems to systematically test all aspects of your force field and whether they behave as they should. There is no quick-do-this-not-that kind of fix for such a situation and guessing or trying random things to suppress the unwanted behavior will only move you away from science and toward doing computer animation. You need to fully understand what is causing the issue before you can try to figure out how to address the problem. This is going to be a laborious process and there is very little that can be done from remote to help.