diatom_scan calculation style

Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.

Introduction

The diatom_scan calculation style evaluates the interaction energy between two atoms at varying distances. This provides a measure of the isolated pair interaction of two atoms providing insights into the strengths of the attraction/repulsion and the effective range of interatomic spacings. This scan also gives insight into the computational smoothness of the potential’s functional form.

Version notes

  • 2019-07-30: Notebook added.

  • 2020-05-22: Version 0.10 update - potentials now loaded from database.

  • 2020-09-22: Setup and parameter definition streamlined. Method and theory expanded.

Additional dependencies

Disclaimers

  • NIST disclaimers

  • No 3+ body interactions are explored with this calculation as only two atoms are used.

Method and Theory

Two atoms are placed in an otherwise empty system. The total energy of the system is evaluated for different interatomic spacings. This provides a means of evaluating the pair interaction component of an interatomic potential, which is useful for a variety of reasons

  • The diatom_scan is a simple calculation that can be used to fingerprint a given interaction. This can be used to help determine if two different implementations produce the same resulting potential when direct comparisons of the potential parameters is not feasible.

  • For a potential to be suitable for radiation studies, the extreme close-range interaction energies must be prohibitively repulsive while not being so large that the resulting force on the atoms will eject them from the system during integration. The diatom_scan results provide a means of evaluating the close-range interactions.

  • The smoothness of the potential is also reflected in the diatom_scan energy results. Numerical derivatives of the measured points can determine the order of smoothness as well as the approximate r values where discontinuities occur.

  • Evaluating large separation values provides a means of identifying the energy of the isolated atoms, given that the separation exceeds the potential’s cutoff. The isolated_atom calculation is an alternative method for obtaining this.

Demonstration

1. Setup

1.1. Library imports

Import libraries needed by the Notebook. The external libraries used are:

[1]:
# Standard library imports
from pathlib import Path
import os
import datetime
from math import floor

# http://www.numpy.org/
import numpy as np

# https://github.com/usnistgov/atomman
import atomman as am
import atomman.lammps as lmp
import atomman.unitconvert as uc

# https://github.com/usnistgov/iprPy
import iprPy

print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)
Notebook last executed on 2020-09-22 using iprPy version 0.10.2

Import additional libraries for plotting. The external libraries used are:

[2]:
import bokeh
print('Bokeh version =', bokeh.__version__)
from bokeh.plotting import figure, output_file, show
from bokeh.embed import components
from bokeh.resources import Resources
from bokeh.io import output_notebook
output_notebook()
Bokeh version = 1.3.4
Loading BokehJS ...

1.2. Default calculation setup

[3]:
# Specify calculation style
calc_style = 'diatom_scan'

# If workingdir is already set, then do nothing (already in correct folder)
try:
    workingdir = workingdir

# Change to workingdir if not already there
except:
    workingdir = Path('calculationfiles', calc_style)
    if not workingdir.is_dir():
        workingdir.mkdir(parents=True)
    os.chdir(workingdir)

# Initialize connection to library
library = iprPy.Library(load=['lammps_potentials'])

2. Assign values for the calculation’s run parameters

2.1. Specify system-specific paths

  • lammps_command is the LAMMPS command to use (required).

  • mpi_command MPI command for running LAMMPS in parallel. A value of None will run simulations serially.

[4]:
lammps_command = 'lmp_serial'
mpi_command = None

2.2. Load interatomic potential

  • potential_name gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.

  • potential is an atomman.lammps.Potential object (required).

[5]:
potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'

# Retrieve potential and parameter file(s)
potential = library.get_lammps_potential(id=potential_name, getfiles=True)

2.3. Specify calculation-specific run parameters

  • symbols is the element or pair of element model symbols to use for the diatom.

  • rmin is the minimum r spacing to use.

  • rmax is the minimum r spacing to use.

  • rsteps is the number of r spacing steps to evaluate.

[6]:
symbols = 'Ni'
rmin = uc.set_in_units(0.02, 'angstrom')
rmax = uc.set_in_units(6.0, 'angstrom')
rsteps = 300

3. Define calculation function(s) and generate template LAMMPS script(s)

3.1. run0.template

[7]:
with open('run0.template', 'w') as f:
    f.write("""#LAMMPS input script that evaluates a system's energy without relaxing

<atomman_system_pair_info>

thermo_style custom step pe
thermo_modify format float %.13e

run 0""")

3.2. diatom()

[8]:
def diatom(lammps_command, potential, symbols,
           mpi_command=None,
           rmin=uc.set_in_units(0.02, 'angstrom'),
           rmax=uc.set_in_units(6.0, 'angstrom'), rsteps=300):
    """
    Performs a diatom energy scan over a range of interatomic spaces, r.

    Parameters
    ----------
    lammps_command :str
        Command for running LAMMPS.
    potential : atomman.lammps.Potential
        The LAMMPS implemented potential to use.
    symbols : list
        The potential symbols associated with the two atoms in the diatom.
    mpi_command : str, optional
        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS
        will run serially.
    rmin : float, optional
        The minimum r spacing to use (default value is 0.02 angstroms).
    rmax : float, optional
        The maximum r spacing to use (default value is 6.0 angstroms).
    rsteps : int, optional
        The number of r spacing steps to evaluate (default value is 300).

    Returns
    -------
    dict
        Dictionary of results consisting of keys:

        - **'r_values'** (*numpy.array of float*) - All interatomic spacings,
          r, explored.
        - **'energy_values'** (*numpy.array of float*) - The computed potential
          energies for each r value.
    """
    # Build filedict if function was called from iprPy
    try:
        assert __name__ == pkg_name
        calc = iprPy.load_calculation(calculation_style)
        filedict = calc.filedict
    except:
        filedict = {}

    # Build lists of values
    r_values = np.linspace(rmin, rmax, rsteps)
    energy_values = np.empty(rsteps)

    # Define atype based on symbols
    symbols = iprPy.tools.aslist(symbols)
    if len(symbols) == 1:
        atype = [1, 1]
    elif len(symbols) == 2:
        atype = [1, 2]
    else:
        raise ValueError('symbols must have one or two values')

    # Initialize system (will shift second atom's position later...)
    box = am.Box.cubic(a = rmax + 1)
    atoms = am.Atoms(atype=atype, pos=[[0.1, 0.1, 0.1], [0.1, 0.1, 0.1]])
    system = am.System(atoms=atoms, box=box, pbc=[False, False, False], symbols=symbols)

    # Add charges if required
    if potential.atom_style == 'charge':
        system.atoms.prop_atype('charge', potential.charges(system.symbols))

    # Get lammps units
    lammps_units = lmp.style.unit(potential.units)

    # Define lammps variables
    lammps_variables = {}

    # Loop over values
    for i in range(rsteps):

        # Shift second atom's x position
        system.atoms.pos[1] = np.array([0.1 + r_values[i], 0.1, 0.1])

        # Save configuration
        system_info = system.dump('atom_data', f='diatom.dat',
                                  potential=potential,
                                  return_pair_info=True)
        lammps_variables['atomman_system_pair_info'] = system_info

        # Write lammps input script
        template_file = 'run0.template'
        lammps_script = 'run0.in'
        template = iprPy.tools.read_calc_file(template_file, filedict)
        with open(lammps_script, 'w') as f:
            f.write(iprPy.tools.filltemplate(template, lammps_variables,
                                             '<', '>'))

        # Run lammps and extract data
        try:
            output = lmp.run(lammps_command, lammps_script, mpi_command)
        except:
            energy_values[i] = np.nan
        else:
            energy = output.simulations[0]['thermo'].PotEng.values[-1]
            energy_values[i] = uc.set_in_units(energy, lammps_units['energy'])

    if len(energy_values[np.isfinite(energy_values)]) == 0:
        raise ValueError('All LAMMPS runs failed. Potential likely invalid or incompatible.')

    # Collect results
    results_dict = {}
    results_dict['r_values'] = r_values
    results_dict['energy_values'] = energy_values

    return results_dict

4. Run calculation function(s)

[9]:
results_dict = diatom(lammps_command, potential, symbols,
                      mpi_command = mpi_command,
                      rmin = rmin,
                      rmax = rmax,
                      rsteps = rsteps)
[10]:
results_dict.keys()
[10]:
dict_keys(['r_values', 'energy_values'])

5. Report results

5.1. Define units for outputting values

  • length_unit is the unit of length to display values in.

  • energy_unit is the unit of energy to display values in.

[11]:
length_unit = 'angstrom'
energy_unit = 'eV'

5.2. Plot E vs r

[12]:
energy = uc.get_in_units(results_dict['energy_values'], energy_unit)
r = uc.get_in_units(results_dict['r_values'], length_unit)

Emin = floor(energy.min())
if Emin < -10:
    Emin = -10

plot = figure(title = f'Diatom energy scan for {potential_name}',
              plot_width = 800,
              plot_height = 600,
              x_range = [uc.get_in_units(rmin, 'angstrom'), uc.get_in_units(rmax, 'angstrom')],
              y_range = [Emin, 0],
              x_axis_label=f'r ({length_unit})',
              y_axis_label=f'Cohesive Energy ({energy_unit}/atom)')

plot.line(r, energy, line_width = 2, legend = symbols)
plot.legend.location = "bottom_right"

show(plot)