E_vs_r_scan calculation style

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


The E_vs_r_scan calculation style calculation creates a plot of the cohesive energy vs interatomic spacing, \(r\), for a given atomic system. The system size is uniformly scaled (\(b/a\) and \(c/a\) ratios held fixed) and the energy is calculated at a number of sizes without relaxing the system. All box sizes corresponding to energy minima are identified.

This calculation was created as a quick method for scanning the phase space of a crystal structure with a given potential in order to identify starting guesses for further structure refinement calculations.

Version notes

  • 2018-07-09: Notebook added.

  • 2019-07-30: Description updated and small changes due to iprPy version.

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

  • 2020-09-22: Setup and parameter definitions streamlined.

Additional dependencies


  • NIST disclaimers

  • The minima identified by this calculation do not guarantee that the associated crystal structure will be stable as no relaxation is performed by this calculation. Upon relaxation, the atomic positions and box dimensions may transform the system to a different structure.

  • It is possible that the calculation may miss an existing minima for a crystal structure if it is outside the range of \(r\) values scanned, or has \(b/a\), \(c/a\) values far from the ideal.

Method and Theory

An initial system (and corresponding unit cell system) is supplied. The \(r/a\) ratio is identified from the unit cell. The system is then uniformly scaled to all \(r_i\) values in the range to be explored and the energy for each is evaluated using LAMMPS and “run 0” command, i.e. no relaxations are performed.

In identifying energy minima along the curve, only the explored values are used without interpolation. In this way, the possible energy minima structures are identified for \(r_i\) where \(E(r_i) < E(r_{i-1})\) and \(E(r_i) < E(r_{i+1})\).


1. Setup

1.1. Library imports

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

# Standard library imports
from pathlib import Path
import os
import datetime
from copy import deepcopy
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:

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
Bokeh version = 1.3.4
Loading BokehJS ...

1.2. Default calculation setup

# Specify calculation style
calc_style = 'E_vs_r_scan'

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

# Change to workingdir if not already there
    workingdir = Path('calculationfiles', calc_style)
    if not workingdir.is_dir():

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

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

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. Load initial unit cell system

  • ucell is an atomman.System representing a fundamental unit cell of the system (required). Here, this is generated using the load parameters and symbols.

# Create ucell by loading prototype record
ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.52)

avect =  [ 3.520,  0.000,  0.000]
bvect =  [ 0.000,  3.520,  0.000]
cvect =  [ 0.000,  0.000,  3.520]
origin = [ 0.000,  0.000,  0.000]
natoms = 4
natypes = 1
symbols = ('Ni',)
pbc = [ True  True  True]
per-atom properties = ['atype', 'pos']
     id   atype  pos[0]  pos[1]  pos[2]
      0       1   0.000   0.000   0.000
      1       1   0.000   1.760   1.760
      2       1   1.760   0.000   1.760
      3       1   1.760   1.760   0.000

2.4. Modify system

  • sizemults list of three integers specifying how many times the ucell vectors of \(a\), \(b\) and \(c\) are replicated in creating system.

  • system is an atomman.System to perform the scan on (required).

sizemults = [3, 3, 3]

# Generate system by supersizing ucell
system = ucell.supersize(*sizemults)
print('# of atoms in system =', system.natoms)
# of atoms in system = 108

2.5. Specify calculation-specific run parameters

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

rmin = uc.set_in_units(2.0, 'angstrom')
rmax = uc.set_in_units(6.0, 'angstrom')
rsteps = 200

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

3.1. run0.template

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


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.2. e_vs_r()

def e_vs_r(lammps_command, system, potential,
           mpi_command=None, ucell=None,
           rmin=uc.set_in_units(2.0, 'angstrom'),
           rmax=uc.set_in_units(6.0, 'angstrom'), rsteps=200):
    Performs a cohesive energy scan over a range of interatomic spaces, r.

    lammps_command :str
        Command for running LAMMPS.
    system : atomman.System
        The system to perform the calculation on.
    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.
    ucell : atomman.System, optional
        The fundamental unit cell correspodning to system.  This is used to
        convert system dimensions to cell dimensions. If not given, ucell will
        be taken as system.
    rmin : float, optional
        The minimum r spacing to use (default value is 2.0 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 200).

        Dictionary of results consisting of keys:

        - **'r_values'** (*numpy.array of float*) - All interatomic spacings,
          r, explored.
        - **'a_values'** (*numpy.array of float*) - All unit cell a lattice
          constants corresponding to the values explored.
        - **'Ecoh_values'** (*numpy.array of float*) - The computed cohesive
          energies for each r value.
        - **'min_cell'** (*list of atomman.System*) - Systems corresponding to
          the minima identified in the Ecoh_values.
    # Build filedict if function was called from iprPy
        assert __name__ == pkg_name
        calc = iprPy.load_calculation(calculation_style)
        filedict = calc.filedict
        filedict = {}

    # Make system a deepcopy of itself (protect original from changes)
    system = deepcopy(system)

    # Set ucell = system if ucell not given
    if ucell is None:
        ucell = system

    # Calculate the r/a ratio for the unit cell
    r_a = r_a_ratio(ucell)

    # Get ratios of lx, ly, and lz of system relative to a of ucell
    lx_a = system.box.a / ucell.box.a
    ly_a = system.box.b / ucell.box.a
    lz_a = system.box.c / ucell.box.a
    alpha = system.box.alpha
    beta =  system.box.beta
    gamma = system.box.gamma

    # Build lists of values
    r_values = np.linspace(rmin, rmax, rsteps)
    a_values = r_values / r_a
    Ecoh_values = np.empty(rsteps)

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

        # Rescale system's box
        a = a_values[i]
        system.box_set(a = a * lx_a,
                       b = a * ly_a,
                       c = a * lz_a,
                       alpha=alpha, beta=beta, gamma=gamma, scale=True)

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

        # Define lammps variables
        lammps_variables = {}
        system_info = system.dump('atom_data', f='atom.dat',
        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
            output = lmp.run(lammps_command, lammps_script, mpi_command)
            Ecoh_values[i] = np.nan
            thermo = output.simulations[0]['thermo']

            if output.lammps_date < datetime.date(2016, 8, 1):
                Ecoh_values[i] = uc.set_in_units(thermo.peatom.values[-1],
                Ecoh_values[i] = uc.set_in_units(thermo.v_peatom.values[-1],

        # Rename log.lammps
            shutil.move('log.lammps', 'run0-'+str(i)+'-log.lammps')

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

    # Find unit cell systems at the energy minimums
    min_cells = []
    for i in range(1, rsteps - 1):
        if (Ecoh_values[i] < Ecoh_values[i-1]
            and Ecoh_values[i] < Ecoh_values[i+1]):
            a = a_values[i]
            cell = deepcopy(ucell)
            cell.box_set(a = a,
                         b = a * ucell.box.b / ucell.box.a,
                         c = a * ucell.box.c / ucell.box.a,
                         alpha=alpha, beta=beta, gamma=gamma, scale=True)

    # Collect results
    results_dict = {}
    results_dict['r_values'] = r_values
    results_dict['a_values'] = a_values
    results_dict['Ecoh_values'] = Ecoh_values
    results_dict['min_cell'] = min_cells

    return results_dict

3.3. r_a_ratio()

def r_a_ratio(ucell):
    Calculates the r/a ratio by identifying the shortest interatomic spacing, r,
    for a unit cell.

    ucell : atomman.System
        The unit cell system to evaluate.

        The shortest interatomic spacing, r, divided by the unit cell's a
        lattice parameter.
    r_a = ucell.box.a
    for i in range(ucell.natoms):
        for j in range(i):
            dmag = np.linalg.norm(ucell.dvect(i,j))
            if dmag < r_a:
                r_a = dmag
    return r_a / ucell.box.a

4. Run calculation function(s)

results_dict = e_vs_r(lammps_command, system, potential,
                      mpi_command = mpi_command,
                      ucell = ucell,
                      rmin = rmin,
                      rmax = rmax,
                      rsteps = rsteps)
dict_keys(['r_values', 'a_values', 'Ecoh_values', 'min_cell'])

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.

length_unit = 'angstrom'
energy_unit = 'eV'

5.2. Plot Ecoh vs r

Ecoh = uc.get_in_units(results_dict['Ecoh_values'], energy_unit)
r = uc.get_in_units(results_dict['r_values'], length_unit)

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

plot = figure(title = f'Cohesive Energy vs. Interatomic Spacing',
              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, Ecoh, line_width = 2)


5.3. List minimum energy configurations

for mincell in results_dict['min_cell']:
    print('Possible minimum near:')
    print('a =', uc.get_in_units(mincell.box.a, length_unit), length_unit)
    print('b =', uc.get_in_units(mincell.box.b, length_unit), length_unit)
    print('c =', uc.get_in_units(mincell.box.c, length_unit), length_unit)
Possible minimum near:
a = 3.510660803076929 angstrom
b = 3.510660803076929 angstrom
c = 3.510660803076929 angstrom

Possible minimum near:
a = 7.376651646951118 angstrom
b = 7.376651646951118 angstrom
c = 7.376651646951118 angstrom