point_defect_static calculation style

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

Description updated: 2019-07-26

Introduction

The point_defect_static calculation style computes the formation energy of a point defect by comparing the energies of a system before and after a point defect is inserted. The resulting defect system is analyzed using a few different metrics to help characterize if the defect reconfigures to a different structure upon relaxation.

Version notes

Additional dependencies

Disclaimers

  • NIST disclaimers

  • The computed values of the point defect formation energies and elastic dipole tensors are sensitive to the size of the system. Larger systems minimize the interaction between the defects, and the affect that the defects have on the system’s pressure. Infinite system formation energies can be estimated by measuring the formation energy for multiple system sizes, and extrapolating to 1/natoms = 0.

  • Because only a static relaxation is performed, the final configuration might not be the true stable configuration. Additionally, the stable configuration may not correspond to any of the standard configurations characterized by the point-defect records in the iprPy/library. Running multiple configurations increases the likelihood of finding the true stable state, but it does not guarantee it. Alternatively, a dynamic simulation or a genetic algorithm may be more thorough.

  • The metrics used to identify reconfigurations are not guaranteed to work for all crystals and defects. Most notably, the metrics assume that the defect’s position coincides with a high symmetry site in the lattice.

  • The current version assumes that the initial defect-free base system is an elemental crystal structure. The formation energy expression needs to be updated to handle multi-component crystals.

Method and Theory

The method starts with a bulk initial system, and relaxes the atomic positions with a LAMMPS simulation that performs an energy/force minimization. The cohesive energy, \(E_{coh}\), is taken by dividing the system’s total energy by the number of atoms in the system.

A corresponding defect system is then constructed using the atomman.defect.point() function. The defect system is relaxed using the same energy/force minimization as was done with the bulk system. The formation energy of the defect, \(E_{f}^{ptd}\), is obtained as

\[E_{f}^{ptd} = E_{total}^{ptd} - E_{coh} * N^{ptd},\]

where \(E_{total}^{ptd}\) is the total potential energy of the relaxed defect system, and \(N^{ptd}\) is the number of atoms in the defect system.

The elastic dipole tensor, \(P_{ij}\), is also estimated for the point defect. \(P_{ij}\) is a symmetric second rank tensor that characterizes the elastic nature of the defect. Here, \(P_{ij}\) is estimated using [1, 2]

\[P_{ij} = -V \langle \sigma_{ij} \rangle,\]

where \(V\) is the system cell volume and \(\langle \sigma_{ij} \rangle\) is the residual stress on the system due to the defect, which is computed using the measured cell stresses (pressures) of the defect-free system, \(\sigma_{ij}^{0}\), and the defect-containing system, \(\sigma_{ij}^{ptd}\)

\[\langle \sigma_{ij} \rangle = \sigma_{ij}^{ptd} - \sigma_{ij}^{0}.\]

The atomman.defect.point() method allows for the generation of four types of point defects:

  • vacancy, where an atom at a specified location is deleted.

  • interstitial, where an extra atom is inserted at a specified location (that does not correspond to an existing atom).

  • substitutional, where the atomic type of an atom at a specified location is changed.

  • dumbbell interstitial, where an atom at a specified location is replaced by a pair of atoms equidistant from the original atom’s position.

Point defect complexes (clusters, balanced ionic defects, etc.) can also be constructed piecewise from these basic types.

The final defect-containing system is analyzed using a few simple metrics to determine whether or not the point defect configuration has relaxed to a lower energy configuration:

  • centrosummation adds up the vector positions of atoms relative to the defect’s position for all atoms within a specified cutoff. In most simple crystals, the defect positions are associated with high symmetry lattice sites in which the centrosummation about that point in the defect-free system will be zero. If the defect only hydrostatically displaces neighbor atoms, then the centrosummation should also be zero for the defect system. This is computed for all defect types.

\[\vec{cs} = \sum_i^N{\left( \vec{r}_i - \vec{r}_{ptd} \right)}\]
  • position_shift is the change in position of an interstitial or substitutional atom following relaxation of the system. A non-zero value indicates that the defect atom has moved from its initially ideal position.

\[\Delta \vec{r} = \vec{r}_{ptd} - \vec{r}_{ptd}^{0}\]
  • db_vect_shift compares the unit vector associated with the pair of atoms in a dumbbell interstitial before and after relaxation. A non-zero value indicates that the dumbbell has rotated from its ideal configuration.

\[\Delta \vec{db} = \frac{\vec{r}_{db1} - \vec{r}_{db2}}{|\vec{r}_{db1} - \vec{r}_{db2}|} - \frac{\vec{r}_{db1}^0 - \vec{r}_{db2}^0}{|\vec{r}_{db1}^0 - \vec{r}_{db2}^0|}\]

If any of the metrics have values not close to (0,0,0), then there was likely an atomic configuration relaxation.

The final defect system and the associated perfect base system are retained in the calculation’s archive. The calculation’s record reports the base system’s cohesive energy, the point defect’s formation energy, and the values of any of the reconfiguration metrics used.

Demonstration

1. Setup

1.1. Library imports

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

[1]:
# Standard library imports
from pathlib import Path
import os
import sys
import uuid
import shutil
import datetime
from copy import deepcopy

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

# https://github.com/usnistgov/DataModelDict
from DataModelDict import DataModelDict as DM

# 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 2019-07-29 using iprPy version 0.9.0

1.2. Default calculation setup

[2]:
# Specify calculation style
calc_style = 'point_defect_static'

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

2. Assign values for the calculation’s run parameters

2.1. Specify system-specific paths

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

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

[3]:
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_file gives the path to the potential_LAMMPS reference record to use. Here, this parameter is automatically generated using potential_name and librarydir.

  • potential_dir gives the path for the folder containing the artifacts associated with the potential (i.e. eam.alloy file). Here, this parameter is automatically generated using potential_name and librarydir.

  • potential is an atomman.lammps.Potential object (required). Here, this parameter is automatically generated from potential_file and potential_dir.

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

# Define potential_file and potential_dir using librarydir and potential_name
potential_file = Path(iprPy.libdir, 'potential_LAMMPS', f'{potential_name}.json')
potential_dir = Path(iprPy.libdir, 'potential_LAMMPS', potential_name)

# Initialize Potential object using potential_file and potential_dir.
potential = lmp.Potential(potential_file, potential_dir)
print('Successfully loaded potential', potential)
Successfully loaded potential 1999--Mishin-Y--Ni--LAMMPS--ipr1

2.3. Load initial unit cell system

  • prototype_name gives the name of the crystal_prototype reference record in the iprPy library to load.

  • symbols is a list of the potential’s elemental model symbols to associate with the unique atom types of the loaded system.

  • box_parameters is a list of the a, b, c lattice constants to assign to the loaded file.

  • load_file gives the path to the atomic configuration file to load for the ucell system. Here, this is generated automatically using prototype_name and librarydir.

  • load_style specifies the format of load_file. Here, this is automatically set for crystal_prototype records.

  • load_options specifies any other keyword options for properly loading the load_file. Here, this is automatically set for crystal_prototype records.

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

[5]:
prototype_name = 'A1--Cu--fcc'
symbols = ['Ni']
box_parameters = uc.set_in_units([3.52, 3.52, 3.52], 'angstrom')

# Define load_file using librarydir and prototype_name
load_file = Path(iprPy.libdir, 'crystal_prototype', f'{prototype_name}.json')

# Define load_style and load_options for crystal_prototype records
load_style = 'system_model'
load_options = {}

# Create ucell by loading prototype record
ucell = am.load(load_style, load_file, symbols=symbols, **load_options)

# Rescale ucell using box_parameters
ucell.box_set(a=box_parameters[0], b=box_parameters[1], c=box_parameters[2], scale=True)

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

2.4. Specify the defect parameters

  • pointdefect_name gives the name of a point-defect reference record in the iprPy library containing point defect input parameters.

  • pointdefect_file gives the path to a point_defect reference containing point defect input parameters. Here, this is built automatically using pointdefect_name and librarydir.

  • point_kwargs (required) is a dictionary or list of dictonaries containing parameters for generating the defect. Here, values are extracted from pointdefect_file. Allowed keywords are:

    • ptd_type indicates which defect type to generate: ‘v’ for vacancy, ‘i’ for interstitial, ‘s’ for substitutional, or ‘db’ for dumbbell.

    • atype is the atom type to assign to the defect atom (‘i’, ‘s’, ‘db’ ptd_types).

    • pos specifies the position for adding the defect atom (all ptd_types).

    • ptd_id specifies the id of an atom in the initial system where the defect is to be added. Alternative to using pos (‘v’, ‘s’, ‘db’ ptd_types).

    • db_vect gives the vector associated with the dumbbell interstitial to generate (‘db’ ptd_type).

    • scale indicates if pos and db_vect are in absolute (False) or box-relative (True) coordinates. Default is False.

    • atol is the absolute tolerance for position-based searching. Default is 1e-3 angstroms.

[6]:
pointdefect_name = 'A1--Cu--fcc--vacancy'
#pointdefect_name = 'A1--Cu--fcc--1nn-divacancy'
#pointdefect_name = 'A1--Cu--fcc--2nn-divacancy'
#pointdefect_name = 'A1--Cu--fcc--100-dumbbell'
#pointdefect_name = 'A1--Cu--fcc--110-dumbbell'
#pointdefect_name = 'A1--Cu--fcc--111-dumbbell'
#pointdefect_name = 'A1--Cu--fcc--octahedral-interstitial'
#pointdefect_name = 'A1--Cu--fcc--tetrahedral-interstitial'
#pointdefect_name = 'A1--Cu--fcc--crowdion-interstitial'

# Define pointdefect_file using librarydir and pointdefect_name
pointdefect_file = Path(iprPy.libdir, 'point_defect', f'{pointdefect_name}.json')

# Parse pointdefect_file using iprPy.input.interpret()
defectinputs = {'ucell':ucell, 'pointdefect_file':pointdefect_file}
iprPy.input.subset('pointdefect').interpret(defectinputs)

# Extract point_kwargs
point_kwargs = defectinputs['point_kwargs']
print('point_kwargs =')
point_kwargs
point_kwargs =
[6]:
[{'ptd_type': 'v', 'pos': array([0., 0., 0.]), 'scale': False}]

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

[7]:
sizemults = [10, 10, 10]

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

2.6. Specify calculation-specific run 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.

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

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

3.1. min.template

[9]:
with open('min.template', 'w') as f:
    f.write("""# LAMMPS input script that performs a simple energy minimization

box tilt large

<atomman_system_info>

<atomman_pair_info>

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

compute peatom all pe/atom

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

min_modify dmax <dmax>

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

3.2. pointdefect()

[10]:
def pointdefect(lammps_command, system, potential, point_kwargs,
                mpi_command=None, etol=0.0, ftol=0.0, maxiter=10000,
                maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom')):
    """
    Adds one or more point defects to a system and evaluates the defect
    formation energy.

    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.
    point_kwargs : dict or list of dict
        One or more dictionaries containing the keyword arguments for
        the atomman.defect.point() function to generate specific point
        defect configuration(s).
    mpi_command : str, optional
        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS
        will run serially.
    sim_directory : str, optional
        The path to the directory to perform the simuation in.  If not
        given, will use the current working directory.
    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:

        - **'E_coh'** (*float*) - The cohesive energy of the bulk system.
        - **'E_ptd_f'** (*float*) - The point.defect formation energy.
        - **'E_total_base'** (*float*) - The total potential energy of the
          relaxed bulk system.
        - **'E_total_ptd'** (*float*) - The total potential energy of the
          relaxed defect system.
        - **'system_base'** (*atomman.System*) - The relaxed bulk system.
        - **'system_ptd'** (*atomman.System*) - The relaxed defect system.
        - **'dumpfile_base'** (*str*) - The filename of the LAMMPS dump file
          for the relaxed bulk system.
        - **'dumpfile_ptd'** (*str*) - The filename of the LAMMPS dump file
          for the relaxed defect system.
    """
    try:
        # Get script's location if __file__ exists
        script_dir = Path(__file__).parent
    except:
        # Use cwd otherwise
        script_dir = Path.cwd()

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

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

    # Define lammps variables
    lammps_variables = {}
    system_info = system.dump('atom_data', f='perfect.dat',
                              units=potential.units,
                              atom_style=potential.atom_style)
    lammps_variables['atomman_system_info'] = system_info
    lammps_variables['atomman_pair_info'] = potential.pair_info(system.symbols)
    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'] = dmax

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

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

    # Run lammps to relax perfect.dat
    output = lmp.run(lammps_command, lammps_script, mpi_command)

    # Extract LAMMPS thermo data.
    thermo = output.simulations[0]['thermo']
    E_total_base = uc.set_in_units(thermo.PotEng.values[-1],
                                   lammps_units['energy'])
    E_coh = E_total_base / system.natoms

    pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
    pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
    pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
    pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])
    pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
    pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
    pressure_base = np.array([[pxx, pxy, pxz], [pxy, pyy, pyz], [pxz, pyz, pzz]])

    # Rename log file
    shutil.move('log.lammps', 'min-perfect-log.lammps')

    # Load relaxed system from dump file and copy old box vectors because
    # dump files crop the values.
    last_dump_file = 'atom.' + str(thermo.Step.values[-1])
    system_base = am.load('atom_dump', last_dump_file, symbols=system.symbols)
    system_base.box_set(vects=system.box.vects)
    system_base.dump('atom_dump', f='perfect.dump')

    # Add defect(s)
    system_ptd = deepcopy(system_base)
    if not isinstance(point_kwargs, (list, tuple)):
        point_kwargs = [point_kwargs]
    for pkwargs in point_kwargs:
        system_ptd = am.defect.point(system_ptd, **pkwargs)

    # Update lammps variables
    system_info = system_ptd.dump('atom_data', f='defect.dat',
                                  units = potential.units,
                                  atom_style = potential.atom_style)
    lammps_variables['atomman_system_info'] = system_info

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

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

    # Extract lammps thermo data
    thermo = output.simulations[0]['thermo']
    E_total_ptd = uc.set_in_units(thermo.PotEng.values[-1],
                                  lammps_units['energy'])
    pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
    pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
    pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
    pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])
    pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
    pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
    pressure_ptd = np.array([[pxx, pxy, pxz], [pxy, pyy, pyz], [pxz, pyz, pzz]])

    # Rename log file
    shutil.move('log.lammps', 'min-defect-log.lammps')

    # Load relaxed system from dump file and copy old vects as
    # the dump files crop the values
    last_dump_file = 'atom.'+str(thermo.Step.values[-1])
    system_ptd = am.load('atom_dump', last_dump_file, symbols=system_ptd.symbols)
    system_ptd.box_set(vects=system.box.vects)
    system_ptd.dump('atom_dump', f='defect.dump')

    # Compute defect formation energy
    E_ptd_f = E_total_ptd - E_coh * system_ptd.natoms

    # Compute strain tensor
    pij = -(pressure_base - pressure_ptd) * system_base.box.volume

    # Cleanup files
    for fname in script_dir.glob('atom.*'):
        fname.unlink()
    for dumpjsonfile in script_dir.glob('*.dump.json'):
        dumpjsonfile.unlink()

    # Return results
    results_dict = {}
    results_dict['E_coh'] = E_coh
    results_dict['E_ptd_f'] = E_ptd_f
    results_dict['E_total_base'] = E_total_base
    results_dict['E_total_ptd'] = E_total_ptd
    results_dict['pij_tensor'] = pij
    results_dict['system_base'] = system_base
    results_dict['system_ptd'] = system_ptd
    results_dict['dumpfile_base'] = 'perfect.dump'
    results_dict['dumpfile_ptd'] = 'defect.dump'

    return results_dict


3.3. check_ptd_config()

[11]:
def check_ptd_config(system, point_kwargs, cutoff,
                     tol=uc.set_in_units(1e-5, 'angstrom')):
    """
    Evaluates a relaxed system containing a point defect to determine if the
    defect structure has transformed to a different configuration.

    Parameters
    ----------
    system : atomman.System
        The relaxed defect system.
    point_kwargs : dict or list of dict
        One or more dictionaries containing the keyword arguments for
        the atomman.defect.point() function to generate specific point
        defect configuration(s).
    cutoff : float
        Cutoff distance to use in identifying neighbor atoms.
    tol : float, optional
        Absolute tolerance to use for identifying if a defect has
        reconfigured (default is 1e-5 Angstoms).

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

        - **'has_reconfigured'** (*bool*) - Flag indicating if the structure
          has been identified as relaxing to a different defect configuration.
        - **'centrosummation'** (*numpy.array of float*) - The centrosummation
          parameter used for evaluating if the configuration has relaxed.
        - **'position_shift'** (*numpy.array of float*) - The position_shift
          parameter used for evaluating if the configuration has relaxed.
          Only given for interstitial and substitutional-style defects.
        - **'db_vect_shift'** (*numpy.array of float*) - The db_vect_shift
          parameter used for evaluating if the configuration has relaxed.
          Only given for dumbbell-style defects.
    """

    # Check if point_kwargs is a list
    if not isinstance(point_kwargs, (list, tuple)):
        pos = point_kwargs['pos']

    # If it is a list of 1, use that set
    elif len(point_kwargs) == 1:
        point_kwargs = point_kwargs[0]
        pos = point_kwargs['pos']

    # If it is a list of two (divacancy), use the first and average position
    elif len(point_kwargs) == 2:
        pos = (np.array(point_kwargs[0]['pos'])
               + np.array(point_kwargs[1]['pos'])) / 2
        point_kwargs = point_kwargs[0]

    # More than two not supported by this function
    else:
        raise ValueError('Invalid point defect parameters')

    # Initially set has_reconfigured to False
    has_reconfigured = False

    # Calculate distance of all atoms from defect position
    pos_vects = system.dvect(system.atoms.pos, pos)
    pos_mags = np.linalg.norm(pos_vects, axis=1)

    # Calculate centrosummation by summing up the positions of the close atoms
    centrosummation = np.sum(pos_vects[pos_mags < cutoff], axis=0)

    if not np.allclose(centrosummation, np.zeros(3), atol=tol):
        has_reconfigured = True

    # Calculate shift of defect atom's position if interstitial or substitutional
    if point_kwargs['ptd_type'] == 'i' or point_kwargs['ptd_type'] == 's':
        position_shift = system.dvect(system.natoms-1, pos)

        if not np.allclose(position_shift, np.zeros(3), atol=tol):
            has_reconfigured = True

        return {'has_reconfigured': has_reconfigured,
                'centrosummation': centrosummation,
                'position_shift': position_shift}

    # Investigate if dumbbell vector has shifted direction
    elif point_kwargs['ptd_type'] == 'db':
        db_vect = point_kwargs['db_vect'] / np.linalg.norm(point_kwargs['db_vect'])
        new_db_vect = system.dvect(-2, -1)
        new_db_vect = new_db_vect / np.linalg.norm(new_db_vect)
        db_vect_shift = db_vect - new_db_vect

        if not np.allclose(db_vect_shift, np.zeros(3), atol=tol):
            has_reconfigured = True

        return {'has_reconfigured': has_reconfigured,
                'centrosummation': centrosummation,
                'db_vect_shift': db_vect_shift}

    else:
        return {'has_reconfigured': has_reconfigured,
                'centrosummation': centrosummation}


4. Run calculation function(s)

4.1. Generate point defect system and evaluate the energy

[12]:
results_dict = pointdefect(lammps_command, system, potential, point_kwargs,
                           mpi_command = mpi_command,
                           etol = energytolerance,
                           ftol = forcetolerance,
                           maxiter = maxiterations,
                           maxeval = maxevaluations,
                           dmax = maxatommotion)
[13]:
results_dict.keys()
[13]:
dict_keys(['E_coh', 'E_ptd_f', 'E_total_base', 'E_total_ptd', 'pij_tensor', 'system_base', 'system_ptd', 'dumpfile_base', 'dumpfile_ptd'])

4.2. Characterize if the defect has reconfigured

[14]:
cutoff = 1.05*ucell.box.a
results_dict.update(check_ptd_config(results_dict['system_ptd'],
                                     point_kwargs,
                                     cutoff))
[15]:
results_dict.keys()
[15]:
dict_keys(['E_coh', 'E_ptd_f', 'E_total_base', 'E_total_ptd', 'pij_tensor', 'system_base', 'system_ptd', 'dumpfile_base', 'dumpfile_ptd', 'has_reconfigured', 'centrosummation'])

5. Report results

5.1. Define units for outputting values

  • length_unit is the unit of length to display relaxed lattice constants in.

  • energy_unit is the unit of energy to display cohesive energies in.

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

5.2. Print \(E_{coh}\) and \(E_{ptd}^f\)

[17]:
print('Ecoh =  ', uc.get_in_units(results_dict['E_coh'], energy_unit), energy_unit)
print('Eptd_f =', uc.get_in_units(results_dict['E_ptd_f'], energy_unit), energy_unit)
Ecoh =   -4.449999998346 eV
Eptd_f = 1.5999212636525046 eV

5.3. Check configuration parameters

[18]:
print('Has the system (likely) reconfigured?', results_dict['has_reconfigured'])
if 'centrosummation' in results_dict:
    print('centrosummation =', uc.get_in_units(results_dict['centrosummation'], length_unit), length_unit)
if 'position_shift' in results_dict:
    print('position_shift = ', uc.get_in_units(results_dict['position_shift'], length_unit), length_unit)
if 'db_vect_shift' in results_dict:
    print('db_vect_shift =  ', uc.get_in_units(results_dict['db_vect_shift'], length_unit), length_unit)
Has the system (likely) reconfigured? False
centrosummation = [-1.89367102e-12 -1.89381844e-12 -1.89359639e-12] angstrom