Quick-Start

The following provides a quick example of how a pair potential model can be defined and then tabulated for different simulation codes. For more advanced features such as EAM potentials, defining custom potential forms, splining or using the python interface, please see User Guide.

In this example we will define Basak’s [Basak2003] potential model for UO2. This has been chosen as it presents an issue that is often solved by using tabulated potentials. The U-O interaction combines Buckingham and Morse potential forms. Although most simulation codes natively provide these as analytical forms, some don’t then allow them to be used in combination for the same pair-interaction. As a result, it becomes necessary to combine them externally and feed them into the code in tabulated form. This tutorial will show how this can be achieved using atsim.potentials.

  1. Installation (see Installation for more detail):

    pip install atsim.potentials
    
  2. Write input file
    • In the Basak model [Basak2003] the O-O and U-U interactions are defined entirely using the Buckingham potential form. Whilst the O-U pair is the sum of a Buckingham and Morse potential (see below for full details of the potential model). Both forms are provided by atsim.potentials and are specified as:
    (1)\[\begin{split}V_\text{Buck}^\text{atsim}(r_{ij}) & = A_{ij} \exp\left( -\frac{r_{ij}}{\rho_{ij}} \right) - \frac{C_{ij}}{r_{ij}^6} \\ V_\text{Morse}^\text{atsim}(r_{ij}) & = D_{ij} \left[ \exp \{ - 2 \gamma_{ij} (r_{ij} - r_{ij}^* \} - 2 \exp \{ - \gamma_{ij} (r_{ij} - r_{ij}^* \} \right] \\\end{split}\]
    • The Buckingham potential is parametrised using values for \(A_ij\), \(\rho_{ij}\) and \(C_{ij}\), specific to each species pair.
    • The Morse potential takes \(D_{ij}\), \(\gamma_{ij}\)and \(r_{ij}^*\).
    • Parameters for the Basak model are given in Table 1.
    Table 1 Basak parameters for use with standard form of Buckingham and Morse potentials
    Parameters O-O U-U O-U
    \(A_{ij}\)/eV 1633.010243 294.640906 693.650934
    \(\rho_{ij}\) 0.327022 0.327022 0.327022
    \(C_{ij}\)/\(\text{eV}Å^6\) 3.948787 0.0 0.0
    \(D_{ij}\)/eV NA NA 0.577190
    \(\gamma_{ij}\)/\(Å^{-1}\) NA NA 1.650000
    \(r_{ij}^*\) NA NA 2.369000
    • The potable tool is used to generate table files.

    • It accepts input in a straightforward format, reminiscent of .ini configuration files.

    • The potential parameters from Table 1 have been transferred into a file suitable for potable: basak.aspot:

      [Tabulation]
      target :  DL_POLY
      cutoff : 6.5
      nr : 652
      
      [Pair]
      O-O = as.buck 1633.010242995040 0.327022 3.948787
      U-U = as.buck 294.640906285709 0.327022 0.0
      O-U = sum(as.buck 693.650933805978 0.327022 0.0, 
      		  as.morse 1.65 2.369 0.577189831995)
      
  3. Generate tabulated potential
    • Download the basak.aspot file and then generate a tabulation in DL_POLY format (see below for information on other formats) by running this command:

      potable basak.aspot TABLE
      
    • This will generate a tabulation in the file named TABLE.

And that’s it! The rest of this page now gives details on what you’ve just done. Read on for more.

What are tabulated potentials?

Pair potentials are functions that relate potential energy to the separation of two interacting particles: \(U_{ij}(r_{ij})\). These are typically defined as equations, with a particular analytical form, that can be tailored to the chemistry of a given pair of species through a set of parameters (e.g. \(A_ij\), \(\rho_{ij}\) and \(C_ij\) in the case of the Buckingham form shown above).

Using these analytical potentials, an entire potential model, comprising of many interactions, can be described in a very compact form. Simulation codes come with large libraries of analytical potential forms, even so, it is impossible for all codes to support all potential forms. As a way of providing flexibility, and as an alternative to requiring users to edit and recompile simulation codes whenever they want to use an unusual form, most support table files.

Tabulated potential approximate a \(U_{ij}(r_{ij})\) functions as a series of [separation, potential-energy] points that are stored in a file, readable by the code. For values not stored in the file, the simulation code performs interpolation to obtain energies and forces for in-between values.

Fig. 1 shows the process followed in producing one of these tables. The desired analytical potential is defined and energy and forces (the first derivative of energy with respect to separation) are sampled at regular intervals. These samples are then written to a file in the format required by the simulation code.

The aim of atsim.potentials is to make this process simple, flexible and transferable across simulation codes.

../_images/tabulated_potential.svg

Fig. 1 Tabulated potentials are obtained by sampling a mathematical formula at regular separations.

Potential model

Before moving on to the tabulation procedure, it is useful to understand what it is we’re trying to tabulate.

Basak’s [Basak2003] model employs the following potential form:

(2)\[V(r_{ij}) = V_\mathrm{Coul}(r_{ij}) + V_\mathrm{Buck}(r_{ij}) + V_\mathrm{Morse}(r_{ij})\]

Where \(V(r_{ij})\) is the potential energy of two atoms (\(i\) and \(j\)) separated by \(r_{ij}\). The first term in eqn. (2) defines the long range electrostatic interaction between the two atoms (with charges \(q_i\) and \(q_j\)):

(3)\[V_\mathrm{Coul}(r_{ij}) = \frac{q_i q_j}{4\pi \epsilon_0 r_{ij}}\]

Where \(\epsilon_0\) is the permittivity of free space.

In a periodic system like UO2, it is necessary to employ various mathematical tricks to get the Coulomb sum to converge quickly and reliably (see for instance Ewald, cell-multipole or particle-mesh methods). As a result the Coulomb part of a potential model isn’t normally included in a tabulation file and therefore it won’t appear further in this example.

This leaves the short-range components of \(V(r_{ij})\) from eqn. (2) for us to bother about, namely \(V_\mathrm{Buck}(r_{ij})\) and \(V_\mathrm{Morse}(r_{ij})\). In Basak’s [Basak2003] paper these are defined as follows:

(4)\[\begin{split}V_\mathrm{Buck}(r_{ij}) & = f_0 b_{ij} \exp\left( \frac{a_{ij} - r_{ij}}{b_{ij}}\right) - \frac{C_{ij}}{r_{ij}^6} \\ V_\mathrm{Morse}(r_{ij}) & = f_0 d_{ij} \left[ \exp \{ - 2 \gamma_{ij} (r_{ij} - r_{ij}^*) \} - 2 \exp \{ - \gamma_{ij} (r_{ij} - r_{ij}^*) \} \right] \\\end{split}\]

The \(f_0\), \(a_{ij}\), \(b_{ij}\), \(C_{ij}\), \(D_{ij}\), \(\gamma_{ij}\) and \(r_{ij}^*\) specifying each pair interaction are given in the following table.

Table 2 Potential parameters for Basak model taken from original paper [Basak2003] (units have been converted for certain values).
Parameters O-O U-U O-U
\(a_{ij}\)/\(Å\) 3.82 3.26 3.54
\(b_{ij}\)/\(Å\) 0.327022 0.327022 0.327022
\(C_{ij}\)/\(\text{eV} Å^6\) 3.948787 0.0 0.0
\(d_{ij}\)/\(r_{ij}^*\) NA NA 13.6765
\(\gamma_{ij}\)/\(r_{ij}^*\) NA NA 1.65
\(r_{ij}^*\) NA NA 2.369
\(f_0\)/\(\text{eV} Å^{-1}\) 0.042203 0.042203 0.042203

The atsim.potentials package comes pre-loaded with a large number of potential-forms, so that you don’t have to constantly redefine functional forms (see List of Potential Forms). In this tutorial you will use the pre-defined versions of the Buckingham and Morse potentials.

If you compare the two definitions of the Buckingham and Morse potentials given in eqns. (1) and (4) you will see there are some differences. The definitions of the Morse potential are almost identical, allowing \(\gamma_{ij}\) and \(r_{ij}^*\) to be used directly with the atsim Morse form. \(D_{ij}\), is however obtained as the product of \(f_0\) and \(d_{ij}\)

(5)\[D_ij = f_0 d_{ij}\]

The Buckingham potential used by Basak has more significant differences to the atsim standard. The \(C_{ij}\) parameter can be used directly in both versions, however we need to manipulate the function from eqn. (4) to show how \(A_{ij}\) and \(\rho_{ij}\) may be obtained from the original Basak parameter set (Table 2).

Let’s do this now by rearranging the Buckingham potential from eqn. (4):

(6)\[\begin{split}V_\mathrm{Buck}(r_{ij}) & = f_0 b_{ij} \exp\left( \frac{a_{ij} - r_{ij}}{b_{ij}}\right) - \frac{C_{ij}}{r_{ij}^6} \\ & = f_0 b_{ij} \exp\left( \frac{a_{ij}}{b_{ij}} - \frac{r_{ij}}{b_{ij}} \right) - \frac{C_{ij}}{r_{ij}^6} \\ & = f_0 b_{ij} \frac{\exp\left( \frac{a_{ij}}{b_{ij}} \right)}{\exp \left( \frac{r_{ij}}{b_{ij}} \right)} - \frac{C_{ij}}{r_{ij}^6} \\ & = f_0 b_{ij} \exp \left( \frac{ a_{ij} }{ b_{ij} } \right) \exp \left( -\frac{r_{ij}}{b_{ij}} \right) - \frac{C_{ij}}{r_{ij}^6} \\\end{split}\]

By comparing coefficients between this equation, eqn. (6) and the atsim form, eqn. (1), it becomes obvious that the Basak parameters can be brought into the standard form using these relationships:

(7)\[\begin{split}A_{ij} & = f_0 b_{ij} \exp \left ( \frac{a_{ij}}{b_{ij}} \right ) \\ \rho_{ij} & = b_{ij} \\\end{split}\]

Using these relationships together with \(D_{ij}\) from eqn. (5) the table of potential parameters given above was obtained (Table 1).

Writing the potential definition

Using the parameters from Table 1, the model was described in the basak.aspot file:

[Tabulation]
target :  DL_POLY
cutoff : 6.5
nr : 652

[Pair]
O-O = as.buck 1633.010242995040 0.327022 3.948787
U-U = as.buck 294.640906285709 0.327022 0.0
O-U = sum(as.buck 693.650933805978 0.327022 0.0, 
		  as.morse 1.65 2.369 0.577189831995)

This file contains two configuration blocks:

  • [Tabulation]: this defines the table’s output format

    • target : DL_POLY means a DL_POLY TABLE file will be produced.
    • cutoff : 6.5 gives the maximum separation to include in the tabulation (\(r_{ij} = 6.5\) Å).
    • nr : 652 the tabulation should contain 652 rows.
  • [Pair]: this section defines the O-O, U-U and O-U pair interactions:

    • The basic form of each line is SPECIES_A-SPECIES_B = POTENTIAL_FORM PARAMETER_1 PARAMETER_2 ... PARAMETER_N
      • Where SPECIES_A and SPECIES_B define each pair of interacting species.
      • POTENTIAL_FORM is a label identifying the functional form of the interaction. Here the as.buck and as.morse types are used. The as. prefix indicates these are standard forms provided by atsim.potentials (see List of Potential Forms for a complete list).

The Buckingham (buck) potential takes three parameters, separated by spaces:

as.buck A ⍴ C

And the Morse (morse) parameters are:

as.morse ɣ r* D

By comparing the input file and Table 1, it should be apparent how the pair-interactions have been parametrised.

The O-U interaction makes use of a potential modifier to combine the as.buck and as.morse forms. This is achieved using the sum() modifier which adds up all the contributions of the comma separated list of potentials defined inside the brackets.

Generating the tabulation

Save the basak.aspot to your drive. Then generate the TABLE file by executing the potable command:

potable basak.aspot TABLE

This will interpret the potential model described above and write it in DL_POLY format to a file named TABLE.

For an example of how this TABLE file can be used in a DL_POLY simulation see Using the TABLE File in DL_POLY.

Specifying other tabulation targets

Once the potential model has been defined, creating tabulations for different codes in different formats is simple.

LAMMPS

To target LAMMPS it is just a case of changing the target option in the [Tabulation] section of the basak.aspot file to LAMMPS e.g.

[Tabulation]
target : LAMMPS
cutoff : 6.5
nr : 652

A LAMMPS table file named Basak.lmptab would then be generated by re-running potable:

potable basak.aspot Basak.lmptab

An example of how to use this file in LAMMPS is given here: Using Basak.lmptab in LAMMPS.

GULP

In the previous section a new tabulation target was specified by editing the basak.aspot file. For temporary changes, the potable allows configuration options to be overridden from the command line. Let’s do this now to make a table file suitable for GULP. This is achieved with the –override-item option, like this:

potable --override-item Tabulation:target=GULP basak.aspot potentials.lib

The potentials.lib file can be used with the following GULP input file to run an energy minimisation: basak.gin