Pair Potential Tabulation

Pair potentials are tabulated using PairTabulation objects, a tabulation class is provided for each target (note an alternative procedural interface is also provided atsim.potentials.writePotentials()):

Using Tabulation objects

The constructor of PairTabulation classes have the following basic signature:

PairTabulation(self, potentials, cutoff, nr)

Where:

Therefore to create a LAMMPS_PairTabulation with a 10 Å cutoff with 5000 rows from a list of potentials stored in the variable potentials the following would be used:

tabulation = LAMMPS_PairTabulation(potentials, 10, 5000)

The pair potential model can then be written to a file-like object using the write() method:

with open("tabulation.lmptab", "w") as outfile:
    tabulation.write(outfile)

See also

Potential Objects

Potential objects should implement the following interface:

class PotentialInterface
speciesA

(str): Attribute giving first species in pair being described by pair-potential

speciesB

(str): Attribute giving second species in pair described by pair-potential

energy(self, r)

Calculate energy between atoms for given separation.

Parameters:r (float) – Separation between atoms of speciesA and speciesB
Returns:Energy in eV for given separation.
Return type:float
force(self, r)

Calculate force (-dU/dr) for interaction at a given separation.

Parameters:r (float) – Separation
Returns:-dU/dr at r in eV per Angstrom.
Return type:float

In most cases the atsim.potentials.Potential class provided in atsim.potentials can be used. This wraps a python callable that returns potential energy as a function of separation to provide the values returned by the energy() method. The forces calculated by the force() method are obtained by taking the numerical derivative of the wrapped function.

Example: Instantiating atsim.potentials.Potential Objects

The following example shows how a Born-Mayer potential function can be described and used to create a Potential object for the interaction between Gd and O. The Born-Mayer potential is given by:

\[U_{\text{Gd-O}}(r_{ij}) = A \exp\left( \frac{- r_{ij}}{{\rho}} \right)\]

Where \(U_{\text{Gd-O}}(r_{ij})\) is the potential energy between atoms \(i\) and \(j\) of types Gd and O, separated by \(r_{ij}\). The parameters \(A\) and \(\rho\) will be taken as 1000.0 and 0.212.

The Gd-O potential function can be defined as:

import math
from atsim.potentials import Potential

def bornMayer_Gd_O(rij):
    energy = 1000.0 * math.exp(-rij/0.212)
    return energy

This is then passed to Potential’s constructor along with the species names:

pot = Potential('Gd', 'O', bornMayer_Gd_O)

The energy and force at a separation of 1Å can then be obtained by calling the energy() and force() methods:

>>> pot.energy(1.0)
8.942132960434881
>>> pot.force(1.0)
42.17987245936639

Predefined Potential Forms

In the previous example, a function named bornMayer_Gd_O() was defined for a single pair-interaction, with the potential parameters hard-coded within the function. Explicitly defining a function for each interaction quickly becomes tedious for anything but the smallest parameter sets. In order to make the creation of functions using standard potential forms easier, a set of function factories are provided within the atsim.potentials.potentialsforms module.

Using the potentialsforms module, the function:

import math

def bornMayer_Gd_O(rij):
    energy = 1000.0 * math.exp(-rij/0.212)
    return energy

can be rewritten as:

from atsim.potentials import potentialforms
bornMayer_Gd_O = potentialsforms.bornmayer(1000.0, 0.212)

See API reference for list of available potential forms: atsim.potentials.potentialforms

Combining Potential Forms

Pair interactions are often described using a combination of standard potential forms. This was seen for the Basak potentials used within the Python API Getting Started example, where the oxygen-uranium pair potential was the combination of a Buckingham and Morse potential forms. This combination was made using the plus() function. This returns a callable which, when invoked, returns the sum of the values returned by the callables originally passed to plus().

The combination functions listed below will return a wrapped function that correctly evaluate the first and second derivatives of the combined callables. That is, when the callables provide .deriv() and .deriv2() methods, these will, where possible be used in the evaluation. In this way accurate analytical derivatives can be combined and will appear in the resulting tabulation. If any callable does not implement these methods, the system will revert to using numerical evaluation of derivatives.

Spline Interpolation

The SplinePotential class can be used to smoothly interpolate between two different potential forms within the same potential curve: one potential function acts below a given cutoff (referred to as the detachment point) and the other potential function takes over at larger separations (acting above a second cutoff called the attachment point). An exponential interpolating spline acts between the detachment and attachment points to provide a smooth transition between the two potential curves.

See also

  • Splining with potable - description of how to do splining with potable rather than using the Python API.

The SplinePotential class aims to automatically determine spline coefficients such that the resultant, interpolated, potential curve is continuous in its first and second derivatives. The analytical form of the interpolating spline is (where \(r_{ij}\) is interatomic separation and \(B_{0..5}\) are the spline coefficients calculated by the SplinePotential class):

\[U(r_{ij}) = \exp \left( B_0 + B_1 r_{ij} + B_2 r_{ij}^2 + B_3 r_{ij}^3 + B_4 r_{ij}^4 + B_5 r_{ij}^5 \right)\]

The SplinePotential has a number of applications, for example:

  • certain potential forms can become attractive in an unphysical manner at small separations (an example is the so-called Buckingham catastrophe); SplinePotential can be used to combine an appropriate repulsive potential at short separations whilst still using the other form for equilibrium and larger separations.
  • similarly different potential forms may be better able to express certain separations than others. For instance the zbl() potential is often used to describe the high energy interactions found in radiation damage cascades but must be combined with another potential to describe equilibrium properties.

The atsim.potentials.spline.Buck4_SplinePotential can also be used to connect two potential functions. The splined region of this potentialform is described via an instance of atsim.potentials.spline.Buck4_Spline.

Both atsim.potentials.spline.Buck4_SplinePotential and atsim.potentials.spline.SplinePotential inherit from atsim.potentials.spline.Custom_SplinePotential. This provides the useful property splineCoefficients which can be used to access the coefficients used to describe the polynomial connecting the two potential functions. These are often quoted in journal papers as they allow the same spline to be reproduced exactly by readers.

Example: Splining ZBL Potential on to Buckingham Potential

As mentioned above, for certain parameterisations, popular potential forms can exhibit unphysical behaviour for some interatomic separations.

See also

A popular model for the description of silicate and phosphate systems is that due to van Beest, Kramer and van Santen (the BKS potential set) [1]. In the current example, the Si-O interaction from this model will be considered. This uses the Buckingham potential form with the following parameters:

  • A = 18003.7572 eV
  • \(\rho\) = 0.205204 Å
  • C = 133.5381 eV \(Å^6\)
  • Charges:
    • Si = 2.4 e
    • O = -1.2 e

The following plot shows the combined coulomb and short-range contributions for this interaction plotted as a function of separation. The large C term necessary to describe the equilibrium properties of silicates means that as \(r_{ij}\) gets smaller, the \(\frac{C}{r_{ij}^6}\) overwhelms the repulsive Born-Mayer component of the Buckingham potential meaning that it turns over. This creates only a relatively shallow minimum arround the equilibrium Si-O separation. Within simulations containing high velocities (e.g. high temperatures or collision cascades) atoms could easily enter the very negative, attractive portion of the potential at low \(r_{ij}\) - effectively allowing atoms to collapse onto each other. In order to overcome this deficiency a ZBL potential will be splined onto the Si-O interaction within this example.

../_images/bks_plot.png

Fig. 8 Plot of BKS Si-O potential showing the short-range (bks_buck) component, electrostatic (bks_coul) and the effective Si-O interaction (bks_buck + bks_coul). This shows that this potential turns over at small separations making it unsuitable for use where high energies may be experienced such as high-temperature or radiation damage cascade simulations.

The first step to using SplinePotential is to choose appropriate detachment and attachment points. This is perhaps best done plotting the two potential functions to be splined. The potentials module contains the convenience functions atsim.potentials.plot() and atsim.potentials.plotToFile() to make this task easier. The following piece of code first defines the ZBL and Buckingham potentials before plotting them into the files zbl.dat and bks_buck.dat. These files each contain two, space delimited, columns giving \(r_{ij}\) and energy, and may be easily plotted in Excel or GNU Plot.

from atsim.potentials import potentialforms
import atsim.potentials

zbl = potentialforms.zbl(14, 8)
bks_buck = potentialforms.buck(18003.7572, 1.0/4.87318, 133.5381)

atsim.potentials.plot( 'bks_buck.dat', 0.1, 10.0, bks_buck, 5000)
atsim.potentials.plot( 'zbl.dat', 0.1, 10.0, zbl, 5000)

Plotting these files show that detachmentX and attachmentX values of 0.8 and 1.4 may be appropriate. The zbl and bks_buck functions can then be splined between these points as follows:

spline = atsim.potentials.SplinePotential(zbl, bks_buck, 0.8, 1.4)

Plot data can then be created for the combined functions with the interpolating spline:

atsim.potentials.plot( 'spline.dat', 0.1, 10.0, spline, 5000)

Plotting the splined Si-O potential together with the original buck and zbl functions allows the smooth transition between the two functions to be observed, as shown in the following function:

../_images/bks_zbl_plot.png

Fig. 9 Plot of BKS Si-O interaction showing the short-range (buck) and ZBL functions plotted with the curve generated by SplinePotential (spline). This joins them with a an interpolating spline acting between the detachment point at \(r_{ij} = 0.8\) and re-attachment point at \(r_{ij} = 1.4\) shown by dashed lines.

Finally, the potential can be tabulated in a format suitable for LAMMPS:

bks_SiO = atsim.potentials.Potential('Si', 'O', spline)
tabulation = atsim.potentials.pair_tabulation.LAMMPS_PairTabulation(
    [bks_SiO],
    10.0, 5000)
with open('bks_SiO.lmptab', 'w') as outfile:
    tabulation.write(outfile)
[1]Van Beest, B. W. H., Kramer, G. J., & van Santen, R. A. (1990). Force fields for silicas and aluminophosphates based on ab initio calculations. Physical Review Letters , 64 (16), 1955–1958. http://dx.doi.org/doi:10.1103/PhysRevLett.64.1955