Python API Getting Started

The following example gives a complete python script showing how the potential API can be used to tabulate potentials for DL_POLY .

See also

  • Also see Quick-Start which allows you to achieve the same using the potable tool without requiring programming experience.

The following example (basak_tabulate.py) shows how the UO2 potential model of Basak [1] can be tabulated:

  • the U + O interaction within this model combines Buckingham and Morse potential forms. Although DL_POLY natively supports both potential forms they cannot be combined with the code itself. By creating a TABLE file the Basak model can be described to DL_POLY.
  • when executed from the command line this script will write tabulated potentials into a file named TABLE.
#! /usr/bin/env python

from atsim.potentials import Potential, plus, potentialforms
from atsim.potentials.pair_tabulation import DLPoly_PairTabulation


def makePotentialObjects():
    # O-O Interaction:
    # Buckingham
    # A = 1633.00510, rho = 0.327022, C = 3.948790
    f_OO = potentialforms.buck(1633.00510, 0.327022, 3.948790)

    # U-U Interaction:
    # Buckingham
    # A = 294.640000, rho = 0.327022, C = 0.0
    f_UU = potentialforms.buck(294.640000, 0.327022, 0.0)

    # O-U Interaction
    # Buckingham + Morse.
    # Buckingham:
    # A = 693.648700, rho = 693.648700, C = 0.0
    # Morse:
    # D0 = 0.577190, alpha = 1.6500, r0 = 2.36900
    buck_OU = potentialforms.buck(693.648700, 0.327022, 0.0)
    morse_OU = potentialforms.morse(1.6500, 2.36900, 0.577190)

    # Compose the buckingham and morse functions into a single function
    # using the atsim.potentials.plus() function
    f_OU = plus(buck_OU, morse_OU)

    # Construct list of Potential objects
    potential_objects = [
        Potential('O', 'O', f_OO),
        Potential('U', 'U', f_UU),
        Potential('O', 'U', f_OU)
    ]
    return potential_objects


def main():
    potential_objects = makePotentialObjects()
    # Tabulate into file called TABLE
    # using short-range cutoff of 6.5 Angs with grid
    # increment of 1e-3 Angs (6500 grid points)

    tabulation = DLPoly_PairTabulation(potential_objects,
                                       6.5, 6500)

    with open('TABLE', 'w') as outfile:
        tabulation.write(outfile)


if __name__ == '__main__':
    main()

Tabulating the Potentials

Defining the Potentials

The first step to tabulating pair potentials is to define Potential objects (see Potential Objects). Normally this involves creating a python function for the desired pair interaction before passing this to the atsim.potentials.Potential() constructor to provide labels for the species pair pertinent to the interaction.

  • The functions f_OO and f_UU use the Buckingham form and are created using buck() function factory (see Predefined Potential Forms for more on the pre-defined forms provided):

    def makePotentialObjects():
        # O-O Interaction:
        # Buckingham
        # A = 1633.00510, rho = 0.327022, C = 3.948790
        f_OO = potentialforms.buck(1633.00510, 0.327022, 3.948790)
    
        # U-U Interaction:
        # Buckingham
        # A = 294.640000, rho = 0.327022, C = 0.0
        f_UU = potentialforms.buck(294.640000, 0.327022, 0.0)
    
  • The O-U interaction is a little more tricky to define as Buckingham and Morse potentials need to be combined. Pre-canned implementations of both of these are provided in atsim.potentials.potentialforms as buck() and morse(). Two functions are created, one for each component of the O-U interaction and stored in the buck_OU and morse_OU variables:

        # O-U Interaction
        # Buckingham + Morse.
        # Buckingham:
        # A = 693.648700, rho = 693.648700, C = 0.0
        # Morse:
        # D0 = 0.577190, alpha = 1.6500, r0 = 2.36900
        buck_OU = potentialforms.buck(693.648700, 0.327022, 0.0)
        morse_OU = potentialforms.morse(1.6500, 2.36900, 0.577190)
    
  • These are then composed into the desired function, f_OU, using the plus() function (see Combining Potential Forms):

        f_OU = plus(buck_OU, morse_OU)
    

Make TABLE File

The table file is written from the main() function of basak_tabulate.py

def main():
    potential_objects = makePotentialObjects()
    # Tabulate into file called TABLE
    # using short-range cutoff of 6.5 Angs with grid
    # increment of 1e-3 Angs (6500 grid points)

    tabulation = DLPoly_PairTabulation(potential_objects,
                                       6.5, 6500)

    with open('TABLE', 'w') as outfile:
        tabulation.write(outfile)
  • First the makePotentialObjects() function is called, returning a list of Potential objects that are stored in the potential_objects variable.

  • An instance of DLPoly_PairTabulation is created by passing this list of potentials a cut-off value of 6.5Å and specifying 6500 rows (i.e. a grid increment of 0.001 Å) to its constructor:

        tabulation = DLPoly_PairTabulation(potential_objects,
                                           6.5, 6500)
    
  • The write() method of the Tabulation object is then called with the file object into which the tabulation is written:

        with open('TABLE', 'w') as outfile:
            tabulation.write(outfile)
    
  • Now run the basak_tabulate.py file (making sure you have installed atsim.potentials first):

    python basak_tabulate.py
    
  • This will create a DL_POLY TABLE file in the working directory.

Using the TABLE File in DL_POLY

A set of DL_POLY files are provided allowing a simple NPT molecular dynamics equilibration simulation to be run against the TABLE file created in the previous step using writePotentials. Copy the files linked from the following table into the same directory as the TABLE file:

File Description
CONFIG 4×4×4 UO2:sub:2 super-cell containing 768 atoms.
CONTROL Defines 300K equilibration run under NPT ensemble lasting 10ps.
FIELD File defining potentials and charges.
  • The FIELD file contains the directives relevant to the TABLE file:

    UO2.cif. Supercell: 4 x 4 x 4
    units eV
    molecules 1
    UO2.cif. Supercell: 4 x 4 x 4
    nummols 1
    atoms 768
                 O     15.999400     -1.200000     512    0
                 U    238.028910      2.400000     256    0
    finish
    vdw 3
    O O tab
    U U tab
    O U tab
    CLOSE
    
  • The following lines define the atom multiplicity and charges (O=-1.2e and U=2.4e):

    nummols 1
    atoms 768
                 O     15.999400     -1.200000     512    0
                 U    238.028910      2.400000     256    0
    finish
    
  • The vdw section states that the O-O, U-U and O-U interactions should be read from the TABLE file:

    vdw 3
    O O tab
    U U tab
    O U tab
    CLOSE
    
  • Once all the files are in the same directory, the simulation can be started by invoking DL_POLY:

    DLPOLY.Z
    

Quick-Start: LAMMPS

Once the potential model has been defined as a series of Potential creating tabulations for different codes in different formats is fairly simple. The script described in this example is given in basak_tabulate_lammps.py. This contains the same potential definition as the previous example, however the main() function has been modified to create a table suitable for LAMMPS :

def main():
    potential_objects = makePotentialObjects()
    # Tabulate into file called Basak.lmptab
    # using short-range cutoff of 6.5 Angs with grid
    # increment of 1e-3 Angs (6500 grid points)
    tabulation = LAMMPS_PairTabulation(potential_objects, 6.5, 6500) # <-- The tabulation class has been changed

    with open('Basak.lmptab', 'w') as outfile:  # <-- Filename changed from 'TABLE'
        tabulation.write(outfile)

Only the two highlighted lines have been changed:

  1. the first changes the tabulation class to LAMMPS_PairTabulation. This describes the same interface as the the previous DLPoly_PairTabulation class meaning it is a drop in replacement.
  2. the second changes the output filename to Basak.lmptab

Running the file creates the Basak.lmptab file:

python basak_tabulate_lammps.py

Using Basak.lmptab in LAMMPS

LAMMPS input files are provided for use with the table file:

Copy these files into the same directory as Basak.lmptab, the simulation can then be run using:

lammps -in equilibrate.lmpin -log equilibrate.lmpout -echo both

The section of equilibrate.lmpin which defines the potential model and makes use of the table file is as follows:

variable O equal 1
variable U equal 2

set type $O charge -1.2
set type $U charge 2.4

kspace_style pppm 1.0e-6

pair_style hybrid/overlay coul/long ${SR_CUTOFF} table linear 6500 pppm
pair_coeff * * coul/long
pair_coeff $O $O table Basak.lmptab O-O
pair_coeff $O $U table Basak.lmptab O-U
pair_coeff $U $U table Basak.lmptab U-U

Notes:

  1. As LAMMPS uses ID numbers to define species the variable commands associate:
    • index 1 with variable $O
    • index 2 with $U to aid readability.
  1. The set type SPECIES_ID charge lines define the charges of oxygen and uranium.
  1. Uses the hybrid/overlay pair_style to combine the coul/long and table styles.

    pair_style hybrid/overlay coul/long ${SR_CUTOFF} table linear 6500 pppm
    
    • The coul/long style is used to calculate electrostatic interactions using the pppm kspace_style defined previously.

    • table linear 6500 pppm:

      • linear interpolation of table values should be used
      • all 6500 rows of the table are employed
      • corrections appropriate to the pppm kspace_style will be applied.
  2. Means that electrostatic interactions should be calculated between all pairs of ions.

    pair_coeff * * coul/long
    
  3. Each pair_coeff reads an interaction from the Basak.lmptab file.

    pair_coeff $O $O table Basak.lmptab O-O
    pair_coeff $O $U table Basak.lmptab O-U
    pair_coeff $U $U table Basak.lmptab U-U
    
  • The general form is:
    • pair_coeff SPECIES_ID_1 SPECIES_ID_2 table TABLE_FILENAME TABLE_KEYWORD
    • Here the SPECIES_IDs use the $O and $U variables defines earlier.
    • TABLE_KEYWORD - the table file contains multiple blocks, each defining a single interaction.
    • The TABLE_KEYWORD is the title of the block. The writePotentials() function creates labels of the form LABEL_A-LABEL_B albeit with the species sorted into alphabetical order. This label format is described in greater detail here.
[1]Basak, C. (2003). Classical molecular dynamics simulation of UO2 to predict thermophysical properties. Journal of Alloys and Compounds, 360 (1-2), 210–216. http://dx.doi.org/doi:10.1016/S0925-8388(03)00350-5