isolated_atom - Methodology and code

Python imports

[1]:
# Standard library imports
from pathlib import Path
import datetime
from typing import Optional, Union

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

# https://ipython.org/
from IPython.display import display, Markdown

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

# https://github.com/usnistgov/iprPy
import iprPy
from iprPy.tools import read_calc_file

print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)
Notebook last executed on 2023-07-31 using iprPy version 0.11.6

1. Load calculation and view description

1.1. Load the calculation

[2]:
# Load the calculation being demoed
calculation = iprPy.load_calculation('isolated_atom')

1.2. Display calculation description and theory

[3]:
# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))

isolated_atom calculation style

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

Introduction

The isolated_atom calculation style evaluates the base energies of all atomic models associated with an interatomic potential. For some potentials, the isolated energy values are necessary to properly compute the cohesive energy of crystal structures. This also provides a simple test whether a potential implementation is compatible with a version of LAMMPS.

Version notes
  • 2020-09-22: Notebook first added.

  • 2022-03-11: Notebook updated to reflect version 0.11.

Additional dependencies
Disclaimers
  • NIST disclaimers

  • Some potentials have two cutoffs with atomic energies outside the first being the “isolated” energy while outside the second have zero energy. The first isolated energy values for those potentials can be found using the diatom_scan calculation instead.

Method and Theory

The calculation loops over all symbol models of the potential and creates a system with a single particle inside a system with non-periodic boundary conditions. The potential energy of each unique isolated atom is evaluated without relaxation/integration.

The cohesive energy, \(E_{coh}\), of a crystal structure is given as the per-atom potential energy of the crystal structure at equilibrium \(E_{crystal}/N\) relative to the potential energy of the same atoms infinitely far apart, \(E_i^{\infty}\)

\[E_{coh} = \frac{E_{crystal} - \sum{N_i E_{i}^{\infty}}}{N},\]

Where the \(N_i\) values are the number of each species \(i\) and \(\sum{N_i} = N\).

For most potentials, \(E_i^{\infty}=0\) meaning that the measured potential energy directly corresponds to the cohesive energy. However, this is not the case for all potentials as some have offsets either due to model artifacts or because it allowed for a better fitted model.

2. Define calculation functions and generate files

This section defines the calculation functions and associated resource files exactly as they exist inside the iprPy package. This allows for the code used to be directly visible and modifiable by anyone looking to see how it works.

2.1. isolated_atom()

This is the primary function for the calculation. The version of this function built in iprPy can be accessed by calling the calc() method of an object of the associated calculation class.

[4]:
def isolated_atom(lammps_command: str,
                  potential: am.lammps.Potential,
                  mpi_command: Optional[str] = None) -> dict:
    """
    Evaluates the isolated atom energy for each elemental model of a potential.

    Parameters
    ----------
    lammps_command :str
        Command for running LAMMPS.
    potential : atomman.lammps.Potential
        The LAMMPS implemented potential to use.
    mpi_command : str, optional
        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS
        will run serially.

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

        - **'energy'** (*dict*) - The computed potential energies for each
          symbol.
    """
    # Initialize dictionary
    energydict = {}

    # Initialize single atom system
    box = am.Box.cubic(a=1)
    atoms = am.Atoms(atype=1, pos=[[0.5, 0.5, 0.5]])
    system = am.System(atoms=atoms, box=box, pbc=[False, False, False])

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

    # Define lammps variables
    lammps_variables = {}

    # Loop over symbols
    for symbol in potential.symbols:
        system.symbols = symbol

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

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

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

        # Run lammps and extract data
        output = lmp.run(lammps_command, script_name=lammps_script,
                         mpi_command=mpi_command)
        energy = output.simulations[0]['thermo'].PotEng.values[-1]
        energydict[symbol] = uc.set_in_units(energy, lammps_units['energy'])

    # Collect results
    results_dict = {}
    results_dict['energy'] = energydict

    return results_dict

2.2. run0.template file

[5]:
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. Specify input parameters

3.1. 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.

[6]:
lammps_command = 'lmp'
mpi_command = None

# Optional: check that LAMMPS works and show its version
print(f'LAMMPS version = {am.lammps.checkversion(lammps_command)["version"]}')
LAMMPS version = 15 Sep 2022

3.2. 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).

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

# Retrieve potential and parameter file(s) using atomman
potential = am.load_lammps_potential(id=potential_name, getfiles=True)

4. Run calculation and view results

4.1. Run calculation

All primary calculation method functions take a series of inputs and return a dictionary of outputs.

[8]:
results_dict = isolated_atom(lammps_command, potential, mpi_command=mpi_command)
print(results_dict.keys())
dict_keys(['energy'])

4.2. Report results

Values returned in the results_dict: - ‘energy’ is a dictionary containing the computed isolated atom energy values for each of the potential’s symbol models.

[9]:
energy_unit = 'eV'
for symbol, energy in results_dict['energy'].items():
    print(symbol, uc.get_in_units(energy, energy_unit), energy_unit)
Ni -3.0970029925997e-11 eV