energy_check - 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('energy_check')

1.2. Display calculation description and theory

[3]:

# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))

energy_check calculation style

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

Idea suggested by Udo v. Toussaint (Max-Planck-Institute f. Plasmaphysics)

Introduction

The energy_check calculation style provides a quick check if the energy of an atomic configuration matches with an expected one.

Version notes

Additional dependencies

Disclaimers

NIST disclaimers
Small variations in the energy are to be expected due to numerical precisions.

Method and Theory

The calculation performs a quick run 0 (no relaxation) energy calculation on a given atomic configuration using a given potential and compares the computed potential energy versus an expected energy value.

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. energy_check()

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 energy_check(lammps_command: str,
                 system: am.System,
                 potential: lmp.Potential,
                 mpi_command: Optional[str] = None) -> dict:
    """
    Performs a quick run 0 calculation to evaluate the potential energy of a
    configuration.

    Parameters
    ----------
    lammps_command :str
        Command for running LAMMPS.
    system : atomman.System
        The atomic configuration to evaluate.
    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:
        - **'E_pot'** (*float*) - The per-atom potential energy of the system.
    """

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

    # Define lammps variables
    lammps_variables = {}
    system_info = system.dump('atom_data', f='init.dat',
                              potential=potential)
    lammps_variables['atomman_system_pair_info'] = system_info

    # Fill in lammps input script
    template = read_calc_file('iprPy.calculation.energy_check', 'run0.template')
    script = filltemplate(template, lammps_variables, '<', '>')

    # Run LAMMPS
    output = lmp.run(lammps_command, script=script,
                     mpi_command=mpi_command, logfile=None)

    # Extract output values
    thermo = output.simulations[-1]['thermo']
    results = {}
    results['E_pot'] = uc.set_in_units(thermo.v_peatom.values[-1],
                                       lammps_units['energy'])

    return results

2.2. run0.template file

[5]:

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

box tilt large

<atomman_system_pair_info>

variable peatom equal pe/atoms

thermo_style custom step lx ly lz pxx pyy pzz pe v_peatom
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)

3.3. System

system is an atomman.System representing a fundamental unit cell of the system (required). Here, this is loaded as the ucell from a relaxed_crystal record.
expected_potential_energy is the expected per-atom potential energy for the system. Not needed for the calculation itself, but used here to compare with the computed value. This is taken from the relaxed_crystal record.

[8]:

# Fetch a relaxed crystal record from the database
potdb = am.library.Database(local=False)
crystal = potdb.get_relaxed_crystal(potential_LAMMPS_id=potential.id, family='A1--Cu--fcc', standing='good')

# Set ucell from the crystal record
system = crystal.ucell

# Set the expected potential energy from the crystal record
expected_potential_energy = crystal.potential_energy

Multiple matching record retrieved from remote
#  family               symbols  alat    Ecoh    method  standing
 1 A1--Cu--fcc          Ni        3.5200 -4.4500 dynamic good
 2 A1--Cu--fcc          Ni        7.3760  0.0119 dynamic good
Please select one:1

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.

[9]:

results_dict = energy_check(lammps_command, system, potential, mpi_command=mpi_command)
print(results_dict.keys())

dict_keys(['E_pot'])

4.2. Report results

Values returned in the results_dict:

‘E_pot’ is the computed average potential energy across all atoms.

[10]:

energy_unit = 'eV'
print('Measured potential energy:', uc.get_in_units(results_dict['E_pot'], energy_unit), energy_unit)
if expected_potential_energy is not None:
    print('Expected potential energy:', uc.get_in_units(expected_potential_energy, energy_unit), energy_unit)

Measured potential energy: -4.4499999983489 eV
Expected potential energy: -4.44999999835575 eV