Greetings, everyone.
My task is to perform PKA simulations in PuxLa(1-x)PO4 crystals. Unfortunately I have been unable to find almost any articles on how it is common to run such simulations in ionic solids (much less code examples in the supplement).
Using Morse-Coulomb potentials it was possible to run a simulation, but I doubt its physicality, so I decided to switch to bucl-coul + zbl.
I was helped to make Buckingham-Coulomb potentials for the set La, Pu, P, O. I managed to optimise the geometry of my supercell in the npt ensemble (p=0, T=300K) in a periodic simulation box. After performing the equilibration, I prescribed write_restart and further worked with it.
The essence of the problem: even setting a relatively low velocity to the Pu atom ( corresponding to10 keV) , the system goes into uncontrolled increase of kinetic energy and temperature. Judging by the animation a few oxygen atoms are knocked out and start frantically flinging around the cell, and Pu itself flies a much shorter distance than it did with the Morse potential. The whole thing ends with a Lost Atoms error.
I varied the time step, set a short inner radius for ZBL, did a frequent recalculation of the list of neighbours, tried to make a Langevin bath for the boundaries, but nothing helped.
Does your simulation account for the short-range problem with Buckingham potentials? That is – the exp(-r) repulsion is finite at short distances while the r^{-6} attraction is not, so there is some cutoff within which two particles experience unbounded attraction.
This would be an especially big problem for bombardment simulations where particles can instantaneously have very high energies. I can’t see a reliable way to prevent particles from falling into the unphysical short-range attraction and then crashing the simulation. I hope you’ve accounted for that.
I had hoped to correct for particle sticking due to the Buckingham potential by using the ZBL nuclear repulsion potential. In my code this is done by hybrid/overlaying the buck/coul/long and zbl potentials. I understand correctly that with this way of setting the potentials we do not need to spline the functions U(r)?
With hybrid/overlay you are just adding two potentials, with a (spline) table you can define a piecewise potential where different interactions are used for different distances. So you are comparing apples and oranges.
You’re absolutely right, Dr Kohlmeyer!
I think I realised the error with ZBL and Buckingham. Because their sum at small r will be the sum of a large positive and a large negative, I definitely need to spline and use the table potential. Thanks for helping me get to this point!
Another issue is that the tabulated potentials I obtained contain a significant number (about 20%) of points where there is a mismatch between U and -grad(U). I generated the potentials with the ATSIM library for Python, the code should be reliable. Do you think it is necessary to solve this problem?
Yes, the junction between the ZBL part and the repulsive part is an important feature. When it is poorly done it can lead to local minima in the potential function that are completely unphysical. The common practice when developping ML potentials is to substract the repulsive part from the data and fit the model (and hence the junction) on that.
In your case this is not an ML potential so I would manually get some points from the repulsive part and manually perform the junction.
There is some python code in tools/tabulate folder and some examples for generating tables. The file pair_hybrid_tabulate.py is an example for a piecewise potential function (repulsive Morse and attractive LJ) with a smooth transition using a switching function.
The file pair_zbladd_tabulate.py shows a case for adding ZBL to an EAM model with pair style hybrid/overlay, where you have to subtract out the pairwise additive component of the EAM potential similar what @tomasfbouvier was referring to. For your purposes, you could use pair_hybrid_tabulate.py as a template and then copy over the ZBL function from the pair _zbladd_tabulate.py and use the buckingham function for the attractive part. To add coulomb interactions, I would use pair style hybrid/overlay so you don’t have to include coulomb in the table file and can easily switch between different Coulomb variants.