# surface_energy_static - Methodology and code

Python imports

[1]:

# Standard library imports
from pathlib import Path
import shutil
import datetime
from copy import deepcopy
from math import floor
from typing import Optional, Tuple, 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

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

[2]:

# Load the calculation being demoed


### 1.2. Display calculation description and theory

[3]:

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


#### surface_energy_static calculation style

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

##### Introduction

The surface_energy_static calculation style evaluates the formation energy for a free surface by slicing an atomic system along a specific plane.

###### Version notes

• 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: Calculation updated to use atomman.defect.FreeSurface class. Setup and parameter definition streamlined.

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

###### Disclaimers
• NIST disclaimers

• Other atomic configurations at the free surface for certain planar cuts may have lower energies. The atomic relaxation will find a local minimum, which may not be the global minimum. Additionally, the material cut is planar perfect and therefore does not explore the effects of atomic roughness.

#### Method and Theory

First, an initial system is generated. This is accomplished by

1. Starting with a unit cell system.

2. Generating a transformed system by rotating the unit cell such that the new system’s box vectors correspond to crystallographic directions, and filled in with atoms to remain a perfect bulk cell when the three boundaries are periodic.

3. All atoms are shifted by a fractional amount of the box vectors if needed.

4. A supercell system is constructed by combining multiple replicas of the transformed system.

Two LAMMPS simulations are then performed that apply an energy/force minimization on the system, and the total energy of the system after relaxing is measured, $$E_{total}$$. In the first simulation, all of the box’s directions are kept periodic (ppp), while in the second simulation two are periodic and one is non-periodic (ppm). This effectively slices the system along the boundary plane creating two free surfaces, each with surface area

$A = \left| \vec{a_1} \times \vec{a_2} \right|,$

where $$\vec{a_1}$$ and $$\vec{a_2}$$ are the two lattice vectors corresponding to the periodic in-plane directions.

The formation energy of the free surface, $$E_{f}^{surf}$$, is computed in units of energy over area as

$E_{f}^{surf} = \frac{E_{total}^{surf} - E_{total}^{0}} {2 A}.$

The calculation method allows for the specification of which of the three box dimensions the cut is made along. If not specified, the default behavior is to make the $$\vec{c}$$ vector direction non-periodic. This choice is due to the limitations of how LAMMPS defines triclinic boxes. $$\vec{c}$$ vector is the only box vector that is allowed to have a component in the Cartesian z direction. Because of this, the other two box vectors are normal to the z-axis and therefore will be in the cut plane.

## 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. surface_energy_static()

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 surface_energy_static(lammps_command: str,
ucell: am.System,
potential: lmp.Potential,
hkl: Union[list, np.ndarray],
mpi_command: Optional[str] = None,
sizemults: Union[list, tuple, None] = None,
minwidth: float = None,
even: bool = False,
conventional_setting: str = 'p',
cutboxvector: str = 'c',
atomshift: Union[list, np.ndarray, None] = None,
shiftindex: Optional[int] = None,
etol: float = 0.0,
ftol: float = 0.0,
maxiter: int = 10000,
maxeval: int = 100000,
dmax: float = uc.set_in_units(0.01, 'angstrom')) -> dict:
"""
Evaluates surface formation energies by slicing along one periodic
boundary of a bulk system.

Parameters
----------
lammps_command :str
Command for running LAMMPS.
ucell : atomman.System
The crystal unit cell to use as the basis of the stacking fault
configurations.
potential : atomman.lammps.Potential
The LAMMPS implemented potential to use.
hkl : array-like object
The Miller(-Bravais) crystal fault plane relative to ucell.
mpi_command : str, optional
The MPI command for running LAMMPS in parallel.  If not given, LAMMPS
will run serially.
sizemults : list or tuple, optional
The three System.supersize multipliers [a_mult, b_mult, c_mult] to use on the
rotated cell to build the final system. Note that the cutboxvector sizemult
must be an integer and not a tuple.  Default value is [1, 1, 1].
minwidth : float, optional
If given, the sizemult along the cutboxvector will be selected such that the
width of the resulting final system in that direction will be at least this
value. If both sizemults and minwidth are given, then the larger of the two
in the cutboxvector direction will be used.
even : bool, optional
A True value means that the sizemult for cutboxvector will be made an even
number by adding 1 if it is odd.  Default value is False.
conventional_setting : str, optional
Allows for rotations of a primitive unit cell to be determined from
(hkl) indices specified relative to a conventional unit cell.  Allowed
settings: 'p' for primitive (no conversion), 'f' for face-centered,
'i' for body-centered, and 'a', 'b', or 'c' for side-centered.  Default
behavior is to perform no conversion, i.e. take (hkl) relative to the
given ucell.
cutboxvector : str, optional
Indicates which of the three system box vectors, 'a', 'b', or 'c', to
cut with a non-periodic boundary (default is 'c').
atomshift : array-like object, optional
A Cartesian vector shift to apply to all atoms.  Can be used to shift
atoms perpendicular to the fault plane to allow different termination
planes to be cut.  Cannot be given with shiftindex.
shiftindex : int, optional
Allows for selection of different termination planes based on the
preferred shift values determined by the underlying fault generation.
Cannot be given with atomshift. If neither atomshift nor shiftindex
given, then shiftindex will be set to 0.
etol : float, optional
The energy tolerance for the structure minimization. This value is
unitless. (Default is 0.0).
ftol : float, optional
The force tolerance for the structure minimization. This value is in
units of force. (Default is 0.0).
maxiter : int, optional
The maximum number of minimization iterations to use (default is
10000).
maxeval : int, optional
The maximum number of minimization evaluations to use (default is
100000).
dmax : float, optional
The maximum distance in length units that any atom is allowed to relax
in any direction during a single minimization iteration (default is
0.01 Angstroms).

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

- **'dumpfile_base'** (*str*) - The filename of the LAMMPS dump file
of the relaxed bulk system.
- **'dumpfile_surf'** (*str*) - The filename of the LAMMPS dump file
of the relaxed system containing the free surfaces.
- **'E_total_base'** (*float*) - The total potential energy of the
relaxed bulk system.
- **'E_total_surf'** (*float*) - The total potential energy of the
relaxed system containing the free surfaces.
- **'A_surf'** (*float*) - The area of the free surface.
- **'E_pot'** (*float*) - The per-atom potential energy of the relaxed bulk
system.
- **'E_surf_f'** (*float*) - The computed surface formation energy.

Raises
------
ValueError
For invalid cutboxvectors
"""
# Construct free surface configuration generator
surf_gen = am.defect.FreeSurface(hkl, ucell, cutboxvector=cutboxvector,
conventional_setting=conventional_setting)

# Check shift parameters
if shiftindex is not None:
assert atomshift is None, 'shiftindex and atomshift cannot both be given'
atomshift = surf_gen.shifts[shiftindex]
elif atomshift is None:
atomshift = surf_gen.shifts[0]

# Generate the free surface configuration
system = surf_gen.surface(shift=atomshift, minwidth=minwidth,
sizemults=sizemults, even=even)
A_surf= surf_gen.surfacearea

# Identify lammps_date version
lammps_date = lmp.checkversion(lammps_command)['date']

# Evaluate system with free surface
surf_results = relax_system(lammps_command, system, potential,
mpi_command=mpi_command, etol=etol, ftol=ftol,
maxiter=maxiter, maxeval=maxeval, dmax=dmax,
lammps_date=lammps_date)

# Extract results from system with free surface
dumpfile_surf = 'surface.dump'
shutil.move(surf_results['finaldumpfile'], dumpfile_surf)
shutil.move('log.lammps', 'surface-log.lammps')
E_total_surf = surf_results['potentialenergy']

# Evaluate perfect system (all pbc removes cut)
system.pbc = [True, True, True]
perf_results = relax_system(lammps_command, system, potential,
mpi_command=mpi_command, etol=etol, ftol=ftol,
maxiter=maxiter, maxeval=maxeval, dmax=dmax,
lammps_date=lammps_date)

# Extract results from perfect system
dumpfile_base = 'perfect.dump'
shutil.move(perf_results['finaldumpfile'], dumpfile_base)
shutil.move('log.lammps', 'perfect-log.lammps')
E_total_base = perf_results['potentialenergy']

# Compute the free surface formation energy
E_surf_f = (E_total_surf - E_total_base) / (2 * A_surf)

# Save values to results dictionary
results_dict = {}

results_dict['dumpfile_base'] = dumpfile_base
results_dict['dumpfile_surf'] = dumpfile_surf
results_dict['E_total_base'] = E_total_base
results_dict['E_total_surf'] = E_total_surf
results_dict['A_surf'] = A_surf
results_dict['E_pot'] = E_total_base / system.natoms
results_dict['E_surf_f'] = E_surf_f

return results_dict


### 2.2. relax_system

[5]:

def relax_system(lammps_command: str,
system: am.System,
potential: lmp.Potential,
mpi_command: Optional[str] = None,
etol: float = 0.0,
ftol: float = 0.0,
maxiter: int = 10000,
maxeval: int = 100000,
dmax: float = uc.set_in_units(0.01, 'angstrom'),
lammps_date: Optional[datetime.date] = None) -> dict:
"""
Sets up and runs the min.in LAMMPS script for performing an energy/force
minimization to relax a system.

Parameters
----------
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.
etol : float, optional
The energy tolerance for the structure minimization. This value is
unitless. (Default is 0.0).
ftol : float, optional
The force tolerance for the structure minimization. This value is in
units of force. (Default is 0.0).
maxiter : int, optional
The maximum number of minimization iterations to use (default is
10000).
maxeval : int, optional
The maximum number of minimization evaluations to use (default is
100000).
dmax : float, optional
The maximum distance in length units that any atom is allowed to relax
in any direction during a single minimization iteration (default is
0.01 Angstroms).
lammps_date : datetime.date or None, optional
The date version of the LAMMPS executable.  If None, will be identified
from the lammps_command (default is None).

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

- **'logfile'** (*str*) - The name of the LAMMPS log file.
- **'initialdatafile'** (*str*) - The name of the LAMMPS data file
used to import an inital configuration.
- **'initialdumpfile'** (*str*) - The name of the LAMMPS dump file
corresponding to the inital configuration.
- **'finaldumpfile'** (*str*) - The name of the LAMMPS dump file
corresponding to the relaxed configuration.
- **'potentialenergy'** (*float*) - The total potential energy of
the relaxed system.
"""

# Ensure all atoms are within the system's box
system.wrap()

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

#Get lammps version date
if lammps_date is None:
lammps_date = lmp.checkversion(lammps_command)['date']

# Define lammps variables
lammps_variables = {}
system_info = system.dump('atom_data', f='system.dat',
potential=potential)
lammps_variables['atomman_system_pair_info'] = system_info
lammps_variables['etol'] = etol
lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])
lammps_variables['maxiter'] = maxiter
lammps_variables['maxeval'] = maxeval
lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])

# Set dump_modify format based on dump_modify_version
if lammps_date < datetime.date(2016, 8, 3):
lammps_variables['dump_modify_format'] = '"%i %i %.13e %.13e %.13e %.13e"'
else:
lammps_variables['dump_modify_format'] = 'float %.13e'

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

# Run LAMMPS
output = lmp.run(lammps_command, script_name=lammps_script,
mpi_command=mpi_command)

# Extract output values
thermo = output.simulations[-1]['thermo']
results = {}
results['logfile'] = 'log.lammps'
results['initialdatafile'] = 'system.dat'
results['initialdumpfile'] = 'atom.0'
results['finaldumpfile'] = 'atom.%i' % thermo.Step.values[-1]
results['potentialenergy'] = uc.set_in_units(thermo.PotEng.values[-1],
lammps_units['energy'])

return results


### 2.3. min.template file

[6]:

with open('min.template', 'w') as f:
f.write("""#LAMMPS input script that performs an energy minimization

box tilt large

<atomman_system_pair_info>

thermo_style custom step lx ly lz pxx pyy pzz pe
thermo_modify format float %.13e

compute peatom all pe/atom

min_modify dmax <dmax>

dump dumpit all custom <maxeval> atom.* id type x y z c_peatom
dump_modify dumpit format <dump_modify_format>

minimize <etol> <ftol> <maxiter> <maxeval>""")


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

[7]:

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

[8]:

potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'

# Retrieve potential and parameter file(s) using atomman


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

[9]:

# Create ucell by loading prototype record

print(ucell)

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


### 3.4. Defect parameters

• hkl gives the Miller (hkl) or Miller-Bravais (hkil) plane to create the free surface on.

• cutboxvector specifies which of the three box vectors (‘a’, ‘b’, or ‘c’) is to be made non-periodic to create the free surface.

• shiftindex can be used for complex crystals to specify different termination planes.

[10]:

hkl = [1,0,0]
cutboxvector = 'c'
shiftindex = 0


### 3.5. System modifications

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

• minwidth specifies a minimum width that the system should be along the cutboxvector direction. The given sizemult in that direction will be increased if needed to ensure that the system is at least this wide.

[11]:

sizemults = [5, 5, 10]
minwidth = uc.set_in_units(0.0, 'angstrom')


### 3.5. Calculation-specific parameters

• energytolerance is the energy tolerance to use during the minimizations. This is unitless.

• forcetolerance is the force tolerance to use during the minimizations. This is in energy/length units.

• maxiterations is the maximum number of minimization iterations to use.

• maxevaluations is the maximum number of minimization evaluations to use.

• maxatommotion is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units.

[12]:

energytolerance = 1e-8
forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')
maxiterations = 10000
maxevaluations = 100000
maxatommotion = uc.set_in_units(0.01, 'angstrom')


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

[13]:

results_dict = surface_energy_static(lammps_command, ucell, potential, hkl,
mpi_command = mpi_command,
sizemults = sizemults,
minwidth = minwidth,
cutboxvector = cutboxvector,
shiftindex = shiftindex,
etol = energytolerance,
ftol = forcetolerance,
maxiter = maxiterations,
maxeval = maxevaluations,
dmax = maxatommotion)
print(results_dict.keys())

dict_keys(['dumpfile_base', 'dumpfile_surf', 'E_total_base', 'E_total_surf', 'A_surf', 'E_pot', 'E_surf_f'])


### 4.2. Report results

Values returned in the results_dict:

• ‘dumpfile_base’ (str) - The filename of the LAMMPS dump file of the relaxed bulk system.

• ‘dumpfile_surf’ (str) - The filename of the LAMMPS dump file of the relaxed system containing the free surfaces.

• ‘E_total_base’ (float) - The total potential energy of the relaxed bulk system.

• ‘E_total_surf’ (float) - The total potential energy of the relaxed system containing the free surfaces.

• ‘A_surf’ (float) - The area of the free surface.

• ‘E_pot’ (float) - The per-atom potential energy of the relaxed bulk system.

• ‘E_surf_f’ (float) - The computed surface formation energy.

[14]:

energy_unit = 'eV'
area_unit = 'nm^2'
energy_area_unit = 'mJ/m^2'

print('E_pot =  ', uc.get_in_units(results_dict['E_pot'], energy_unit), energy_unit)
print('A_surface =', uc.get_in_units(results_dict['A_surf'], area_unit), area_unit)
print('E_surface_f =', uc.get_in_units(results_dict['E_surf_f'], energy_area_unit), energy_area_unit)

E_pot =   -4.4499999983449 eV
A_surface = 3.0975990100882553 nm^2
E_surface_f = 1877.9483735246017 mJ/m^2

[15]:

# Show info for the base systems
print('E_total_base =', uc.get_in_units(results_dict['E_total_base'], energy_unit), energy_unit)
print('Dumped to', results_dict['dumpfile_base'])

E_total_base = -4449.9999983449 eV
Dumped to perfect.dump

[16]:

# Show info for the base systems
print('E_total_surf =', uc.get_in_units(results_dict['E_total_surf'], energy_unit), energy_unit)
print('Dumped to', results_dict['dumpfile_surf'])

E_total_surf = -4377.3846462123 eV
Dumped to surface.dump