# Combining potentials¶

It is often useful to describe the different regions of a pair-potential using different functional forms. This section of the user-guide will show ways in which different potentials can be combined.

First, examples will be given demonstrating the use multi-range potentials described earlier.

Methods for linking potential-forms using splines will then be considered.

## Multi-range potential examples¶

This section shows how potential-forms can be combined using the multi-range syntax introduced here: Multi-range potentials.

### Example: truncating a potential¶

The starting structure for a simulation may contain overlapping atoms (for example if generated from random coordinates). In these cases it is useful to move atoms apart before starting the simulation proper.

In principle this should be possible by using an energy minimisation run, however some potential forms are badly behaved at very small separations (for instance the Buckingham catastrophe) leading to unexpected results.

One way of removing atom overlap is to temporarily substitute the potentials with some that are guaranteed to work, even when two atoms are completely on top of each other. An example of a potential-form that could be used for this purpose is now described (this is available natively in LAMMPS as pair_style soft):

$V_{ij}(r_{ij}) = A \left[ 1+ \cos \left( \frac{\pi r_{ij}}{r_c} \right) \right]$

Let’s start by plotting this function with sensible parameter values of $$A$$ = 10 and $$r_c$$ = 1.6:

To make this useful, the potential is truncated at $$r_c$$ (i.e. $$r_{ij} < r_c$$). This means that the potential is repulsive below $$r_c$$ and doesn’t act above:

Now, having been suitably truncated, the potential will gently push atoms apart in an MD or energy minimisation run. This can be implemented for potable as follows; first define the cosine function:

[Potential-Form]
soft(r, A, rc) = A * (1+cos((pi*r)/rc))


The truncation at $$r_c$$ can then be applied using a multi-range [Pair] deinition.

[Pair]
Si-O : soft 10.0 1.6 >1.6 as.zero
O-O  : soft  5.0 2.4 >2.4 as.zero


Here the Si-O interaction have $$A$$ = 10.0 and an $$r_c$$ value of 1.6Å (the position of the Si-O peak in a silica glass radial distribution function) after which the Zero (zero) potential-form takes over. Similarly, the O-O pair has a weaker repulsion ($$A=5.0$$) and a cutoff-radius, $$r_c$$, of 2.4Å (again this is about the position of the O-O peak in the RDF for silica glass).

In this way the cosine function has been truncated. Specifying an Si-Si interaction will be left as an exercise for the reader.

The complete input file for this example can be downloaded here soft_a.aspot and is:

[Tabulation]
target : GULP
cutoff : 10.0
dr : 0.01

[Pair]
Si-O : soft 10.0 1.6 >1.6 as.zero
O-O  : soft  5.0 2.4 >2.4 as.zero

[Potential-Form]
soft(r, A, rc) = A * (1+cos((pi*r)/rc))


### Example: truncating using if() statement¶

This example shows an alternative to using the multi-range syntax to implement the soft cosine potential described above. The exprtk package used to evaluate expressions in the [Potential-Form] section, support an if ... then ... else construct through its if() function. This has the form:

if (CONDITION, FORMULA_IF_TRUE, FORMULA_IF_FALSE)


This means the $$r > r_c$$ condition could have been included in [Potential-Form] like this:

[Potential-Form]
cos_form(r, A, rc) = A * (1+cos((pi*r)/rc))
soft(r, A, rc) = if(r>rc, 0, cos_form(r, A, rc))


This simplifies the definition of the [Pair] section by avoiding the repetition of the $$r_c$$ parameter when defining the multi-range version of the potential:

[Pair]
Si-O : soft 10.0 1.6
O-O  : soft  5.0 2.4


The potable input using this version of the potential-form can be downloaded here: soft_if.aspot.

### Example: parametrising a model using published spline coefficients¶

The UO2 potential model by Morelon describes U-O and O-O interactions [Morelon2003]. The former employs the Born-Mayer (bornmayer) potential-form while the latter uses the combination of a 3rd and 5th order polynomial spline to link an $$A \exp\left( - \frac{r_{ij}}{\rho}\right)$$ term to a $$-\frac{C}{r_{ij}^6}$$ term. The spline-coefficients for the Morelon model are available in a paper by Potashnikov et al. [Potashnikov2011]. They give the O-O interaction as:

(8)$\begin{split}V(r_{ij}) = \begin{cases} 11272.6 \exp\left(\frac{0.1363}{r_{ij}}\right) &: r_{ij} < 1.2 \\ 479.955 -1372.53 r_{ij} +1562.22 r_{ij}^2 -881.969 r_{ij}^3 +246.435 r_{ij}^4 -27.2447 r_{ij}^5 &: 1.2 \leq r_{ij} < 2.1 \\ 42.8917 -55.4965 r_{ij} +23.0774 r_{ij}^2 -31.13140 r_{ij}^3 &: 2.1 \leq r_{ij} < 2.6 \\ -\frac{134}{r_{ij}^6} &: r_{ij} \geq 2.6 \\ \end{cases}\end{split}$

The following figure plots the O-O interaction to show how its constituent parts relate to each other.

Fig. 4 The O-O interaction of the Morelon model is made up of several potential-forms acting over different ranges. The upper plot shows the un-trimmed versions of the individual functions. The lower plot shows how these combine to form the overall interaction.

Note: in a simulation the O-O interaction would also include a repulsive electrostatic term which isn’t plotted here.

Using this information we can write the Morelon model as potable input (and can be downloaded as morelon.aspot) :

[Tabulation]
target : LAMMPS
cutoff : 10.0
nr : 1001

[Pair]
O-U : as.bornmayer 566.498 0.42056
O-O : as.bornmayer 11272.6 0.1363
>1.2
as.polynomial 479.955 -1372.53 1562.22 -881.969 246.435 -27.2447
>2.1
as.polynomial 42.8917 -55.4965 23.0774 -3.13140
>2.6
as.buck 0.0 1.0 134.0


There are a few features here that are worth noting:

• The four distinct regions are defined as multi-range potentials.

• The third and fifth order polynomials are both defined using the as.polynomial form:

• The order of the polynomial is implicitly defined by the number of parameters.
• The $$-\frac{C}{r_{ij}^6}$$ term is implemented using the Buckingham (buck) form. In order to only use the attractive part of the potential an A parameter of 0.0 is used as its first argument and a non-zero rho is then specified.

## Splining¶

Splines are curves that are used to smoothly connect different functional forms across different ranges of a potential-form. Typically splining curves are polynomials. Care must be taken when determining the coefficients of the spline function as it is important that potentials are smooth and do not have steps in their derivatives. For this reason spline parameters are determined so that at its end-points, it gives the same potential-energy and derivatives (1st and 2nd) as the functions it joins.

In the previous example (see Example: parametrising a model using published spline coefficients), spline coefficients were provided but their calculation can be quite involved. As a result potable provides services to make this easier.

Splining is performed in the [Pair] section of the input through the spline() potential-modifier.

spline() takes a single argument which has the form of a multi-range potential definition. This has three distinct parts, representing the three ranges of the splined potential:

1. The region described by the starting potential.
2. The interpolation region given by the spline.
3. The final part defined by the end potential.

In the description that follows, the separation at which the starting potential ends will be referred to as the detachment-point. Similarly, the end of the splined region is the attachment-point.

Bearing in mind its similarity to the multi-range syntax spline() definitions have this basic form:

spline(POT_A_DEFN >R_DETACH SPLINE_DEFN >R_ATTACH POT_B_DEFN)


Where:

• POT_A_DEFN - this is the potential definition for the starting-potential. This is typically a potential-form label followed by potential parameters (see [Pair] section).
• R_DETACH - detachment-point (note: as this is a multi-range potential definition the inclusive form >=R_DETACH is also valid).
• SPLINE_DEFN - the type of spline to be used along with any parameters it takes. This will be described properly in the next paragraph.
• R_ATTACH - attachment point (again >=R_ATTACH is also accepted).
• POT_B_DEFN - the potential definition of the end-potential.

SPLINE_DEFN has the same format as a potential-definition: an identifying label followed by a list of parameters (note: some spline-types do not take any parameters):

SPLINE_LABEL PARAM_1 PARAM_2 ... PARAM_N


Currently SPLINE_LABEL can be exp_spline or buck4_spline. These spline-types will now be described.

### Exponential Spline exp_spline¶

The spline used when exp_spline is specified is given as:

(9)$V(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)$

Where $$B_{0..5}$$ are the spline coefficients calculated automatically during splining.

The exponential spline does not take any parameters.

The functional form given in eqn. (9) is also available as a potential-form. Although the user is then responsible for calculating their own spline-coefficients. See: Exponential Spline (exp_spline).

#### Example: splining to the zbl potential form using exp_spline¶

This example will show how to use spline() with the exp_spline type to add the repulsive Ziegler-Biersack-Littmark (zbl) to a Buckingham (buck) potential that is unphysical at small separations.

For certain parameterisations, popular potential forms can exhibit unphysical behaviour for some interatomic separations. 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) [VanBeest1990]. 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 plot in Fig. 5 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 around 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.

Fig. 5 Plot showing Si-O pair-potential and its electrostatic and short-range components. A small separations it becomes attractive which is unphysical.

The first step to using exp_spline is to choose appropriate detachment and attachment points. This is perhaps best done by plotting the two potential functions to be splined. The as.buck and as.zbl curves in Fig. 6a, show that detachment and attachment at separations of 0.8 and 1.4 may be appropriate.

Fig. 6 The result of splining using the exp_spline potential.

We now have enough information to spline the two potentials together. This gives is defined in the [Pair] section:

[Pair]
Si-O : spline(
as.zbl 14 8
>=0.8
exp_spline
>=1.4
as.buck 18003.7572 0.205204 133.5381 )


The result of splining can be seen in Fig. 6b (here it plotted with the Coulomb interaction included).

The complete potable input can be downloaded here: exp_spline.aspot.

### Buckingham-4 Spline buck4_spline¶

When used with buck4_spline the overall potential has the form:

(10)$\begin{split}V(r_{ij}) = \begin{cases} V_\text{PotA}(r_{ij}) & : 0 \leq r_{ij} \leq r_\text{detach}\\ a_0 + a_1 r_{ij} +a_2 r_{ij}^2+a_3 r_{ij}^3+a_4 r_{ij}^4+a_5 r_{ij}^5 & : r_\text{detach} < r_{ij} < r_\text{min}\\ b_0 +b_1 r_{ij}+b_2 r_{ij}^2+b_3 r_{ij}^3 & : r_\text{min} \leq r_{ij} < r_\text{attach}\\ V_\text{PotB}(r_{ij}) & : r_{ij} \geq r_\text{attach}\\ \end{cases}\end{split}$

The spline form can be seen to be the combination of a fifth order and cubic polynomial. As for other spline types the coefficients ($$a_{0 \ldots 5}$$ and $$b_{0 \ldots 3}$$) are chosen so that first and second derivatives (with respect to separation) are continuous through the attachment and detachment points. In addition, the coefficients are determined to give a stationary point where the two polynomials meet at $$r_\text{min}$$.

The SPLINE_DEFN part of the spline() definition for buck4_spline is:

buck4_spline $$r_\text{min}$$

#### Relationship to as.buck4 potential-form¶

As its name suggests, the buck4_spline is closely related to the Buckingham-4 (buck4) potential-form. When used as a potential-form, the $$V_\text{PotA}(r_{ij})$$ and $$V_\text{PotB}(r_{ij})$$ terms from eqn. (10) are pre-defined as the first and last terms of the Buckingham potential:

$\begin{split}V_\text{PotA}(r_{ij}) & = A_{ij} \exp\left(- \frac{r_{ij}}{\rho} \right) \\ V_\text{PotB}(r_{ij}) & = - \frac{C}{r_{ij}} \\\end{split}$

The potential-form is therefore shorthand for the following:

spline(as.buck A RHO 0.0 >R_DETACH buck4_spline R_MIN >R_ATTACH as.buck 0 1.0 C)


And would be specified as:

as.buck4 A RHO C R_DETACH R_MIN R_ATTACH


Where: A, RHO and C are the potential parameters and R_DETACH, R_MIN and R_ATTACH define the splined region.

In most cases users will prefer the as.buck4 form unless they have a specific reason to change the initial and final potential functions.

#### Example: redefining the Morelon model using buck4_spline¶

In this example the Morelon potential model will be re-implemented first using the spline() modifier then with the as.buck4.

This model was discussed in a previous example: Example: parametrising a model using published spline coefficients where O-O and O-U interactions were defined like this:

[Pair]
O-U : as.bornmayer 566.498 0.42056
O-O : as.bornmayer 11272.6 0.1363
>1.2
as.polynomial 479.955 -1372.53 1562.22 -881.969 246.435 -27.2447
>2.1
as.polynomial 42.8917 -55.4965 23.0774 -3.13140
>2.6
as.buck 0.0 1.0 134.0


From this it can be seen that the O-U interaction is relatively simple whilst the O-O interaction has four distinct parts (also defined in eqn. (8)), and when examined more closely it can be seen it matches the four range Buckingham form.

The polynomial terms, with their pre-supplied coefficients, can be replaced by using the spline() modifier with the buck4_spline spline type. This is now shown and can be downloaded as: morelon_buck4_spline.aspot

 1 2 3 4 5 6 7 8 [Pair] O-U : as.bornmayer 566.498 0.42056 O-O : spline( as.bornmayer 11272.6 0.1363 >1.2 buck4_spline 2.1 >2.6 as.buck 0.0 1.0 134.0 ) 

To change the starting potential and end potential lines 4 and 8 would be edited. To edit the detach and attach points then lines 5 and 7 would be altered. To change the position of the stationary point in the splined region, the value of 2.1 at the end of line 6 would be edited.

This input file could also be re-written to use the buck4 potential-form, like this:

[Pair]
O-U : as.bornmayer 566.498 0.42056
O-O : as.buck4 11272.6 0.1363 134.0 1.2 2.1 2.6


This version of the model can be downloaded here: morelon_buck4.aspot