point_defect_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, Union, Tuple
# 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('point_defect_static')
1.2. Display calculation description and theory
[3]:
# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))
point_defect_static calculation style
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
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
2020-12-30 Version 0.10+ update
2022-03-11: Notebook updated to reflect version 0.11.
Additional dependencies
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
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]
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}\)
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.
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.
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.
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.
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. calc()
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 calc(lammps_command: str,
system: am.System,
potential: lmp.Potential,
point_kwargs: Union[list, dict],
cutoff: float,
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'),
tol: float = uc.set_in_units(1e-5, 'angstrom')) -> dict:
"""
Adds one or more point defects to a system and evaluates the defect
formation energy. Evaluates a relaxed system containing a point defect
to determine if the defect structure has transformed to a different
configuration.
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).
cutoff : float
Cutoff distance to use in identifying neighbor atoms.
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 simulation 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).
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:
- **'E_pot'** (*float*) - The per-atom potential 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.
- **'pij_tensor'** (*numpy.ndarray of float*) - The elastic dipole
tensor associated with the defect.
- **'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.
- **'has_reconfigured'** (*bool*) - Flag indicating if the structure
has been identified as relaxing to a different defect configuration.
- **'centrosummation'** (*numpy.ndarray of float*) - The centrosummation
parameter used for evaluating if the configuration has relaxed.
- **'position_shift'** (*numpy.ndarray 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.ndarray of float*) - The db_vect_shift
parameter used for evaluating if the configuration has relaxed.
Only given for dumbbell-style defects.
"""
# Run ptd_energy to refine values
results_dict = pointdefect(lammps_command,
system,
potential,
point_kwargs,
mpi_command = mpi_command,
etol = etol,
ftol = ftol,
maxiter = maxiter,
maxeval = maxeval,
dmax = dmax)
# Run check_ptd_config
results_dict2 = check_ptd_config(results_dict['system_ptd'],
point_kwargs,
cutoff, tol)
results_dict.update(results_dict2)
return results_dict
2.2. pointdefect
[5]:
def pointdefect(lammps_command: str,
system: am.System,
potential: lmp.Potential,
point_kwargs: Union[list, dict],
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')) -> dict:
"""
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 simulation 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_pot'** (*float*) - The per-atom potential 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.
- **'pij_tensor'** (*numpy.ndarray of float*) - The elastic dipole
tensor associated with the defect.
- **'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.
"""
# 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',
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'] = 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
lammps_script = 'min.in'
template = read_calc_file('iprPy.calculation.point_defect_static',
'min.template')
with open(lammps_script, 'w') as f:
f.write(filltemplate(template, lammps_variables, '<', '>'))
# Run lammps to relax perfect.dat
output = lmp.run(lammps_command, script_name=lammps_script,
mpi_command=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_pot = 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',
potential=potential)
lammps_variables['atomman_system_pair_info'] = system_info
# Write lammps input script
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 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_pot * system_ptd.natoms
# Compute strain tensor
pij = -(pressure_base - pressure_ptd) * system_base.box.volume
# Cleanup files
for fname in Path.cwd().glob('atom.*'):
fname.unlink()
for dumpjsonfile in Path.cwd().glob('*.dump.json'):
dumpjsonfile.unlink()
# Return results
results_dict = {}
results_dict['E_pot'] = E_pot
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
2.3. check_ptd_config()
[6]:
def check_ptd_config(system: am.System,
point_kwargs: Union[list, dict],
cutoff: float,
tol: float = uc.set_in_units(1e-5, 'angstrom')) -> dict:
"""
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.ndarray of float*) - The centrosummation
parameter used for evaluating if the configuration has relaxed.
- **'position_shift'** (*numpy.ndarray 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.ndarray 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}
2.4. min.template file
[7]:
with open('min.template', 'w') as f:
f.write("""# LAMMPS input script that performs a simple energy minimization
box tilt large
<atomman_system_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. 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.
[8]:
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).
[9]:
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. 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.
[10]:
# Create ucell by loading prototype record
ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc')
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
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.
[11]:
# Add a vacancy by deleting the atom at [0,0,0] - works for any crystal structure
point_kwargs = [{
'ptd_type': 'v',
'pos': np.array([0.0, 0.0, 0.0]),
'scale': True,
}]
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.
system is an atomman.System to perform the scan on (required).
[12]:
sizemults = [12, 12, 12]
# Generate system by supersizing ucell
system = ucell.supersize(*sizemults)
print('# of atoms in system =', system.natoms)
# of atoms in system = 6912
3.6 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.
[13]:
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.
[14]:
results_dict = calc(lammps_command, system, potential, point_kwargs,
cutoff = ucell.r0() * 1.1,
mpi_command = mpi_command,
etol = energytolerance,
ftol = forcetolerance,
maxiter = maxiterations,
maxeval = maxevaluations,
dmax = maxatommotion)
print(results_dict.keys())
dict_keys(['E_pot', 'E_ptd_f', 'E_total_base', 'E_total_ptd', 'pij_tensor', 'system_base', 'system_ptd', 'dumpfile_base', 'dumpfile_ptd', 'has_reconfigured', 'centrosummation'])
4.2. Report results
Values returned in the results_dict:
‘E_pot’ (float) - The per-atom potential 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.
‘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.
[15]:
energy_unit = 'eV'
print('E_pot =', uc.get_in_units(results_dict['E_pot'], energy_unit), energy_unit)
print('E_ptd_f =', uc.get_in_units(results_dict['E_ptd_f'], energy_unit), energy_unit)
E_pot = -4.449999998376013 eV
E_ptd_f = 1.6000474976244732 eV
[16]:
# Show info for the base systems
print(results_dict['system_base'].natoms, 'atoms in base system')
print('E_total_base =', uc.get_in_units(results_dict['E_total_base'], energy_unit), energy_unit)
print('Dumped to', results_dict['dumpfile_base'])
6912 atoms in base system
E_total_base = -30758.399988775 eV
Dumped to perfect.dump
[17]:
# Show info for the ptd systems
print(results_dict['system_ptd'].natoms, 'atoms in ptd system')
print('E_total_ptd =', uc.get_in_units(results_dict['E_total_ptd'], energy_unit), energy_unit)
print('Dumped to', results_dict['dumpfile_ptd'])
6911 atoms in ptd system
E_total_ptd = -30752.349941279 eV
Dumped to defect.dump
[18]:
# Show the computed elastic dipole tensor
pij = deepcopy(results_dict['pij_tensor'])
print(f'pij ({energy_unit})=')
print(uc.get_in_units(pij, energy_unit))
print()
# Filter out near-zero values
pij[np.isclose(pij, 0.0)] = 0.0
print(f'pij refined ({energy_unit})=')
print(uc.get_in_units(pij, energy_unit))
pij (eV)=
[[-2.76326923e+00 4.59715560e-13 3.33427702e-13]
[ 4.59715560e-13 -2.76326923e+00 -6.90836219e-14]
[ 3.33427702e-13 -6.90836219e-14 -2.76326923e+00]]
pij refined (eV)=
[[-2.76326923 0. 0. ]
[ 0. -2.76326923 0. ]
[ 0. 0. -2.76326923]]
[19]:
length_unit = 'angstrom'
# Show characterization anaysis
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 = [5.89576226e-14 5.77315973e-14 5.77315973e-14] angstrom