dislocation_periodic_array - 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
import random
# 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('dislocation_periodic_array')
1.2. Display calculation description and theory
[3]:
# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))
dislocation_periodic_array calculation style
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
Introduction
The dislocation_periodic_array calculation constructs an atomic system with a periodic array of dislocations configuration. A single dislocation is inserted into an otherwise perfect crystal, and the system is kept periodic in the two system box directions that are within the dislocation’s slip plane. The system is then statically relaxed with the atoms at the boundary perpendicular to the slip plane held fixed.
Version notes
2019-07-30: Notebook added.
2020-05-22: Notebook updated for iprPy version 0.10 and tested for hcp
2020-09-22: Notebook updated to reflect that calculation method has changed to now use atomman.defect.Dislocation. Setup and parameter definition cleaned up and streamlined.
2022-03-11: Notebook updated to reflect version 0.11.
Additional dependencies
Disclaimers
This calculation was designed to be general enough to properly generate a dislocation for any crystal system but has not been fully tested yet for extreme cases.
Method and Theory
System orientation considerations
Properly constructing a periodic array of dislocations atomic configuration requires careful consideration of dislocation solutions and atomic system boundaries. Solutions for straight dislocations based on elasticity often follow the convention of using a Cartesian system (\(x', y', z'\)) in which the dislocation line is oriented along the \(z'\)-axis, and the slip plane taken to be the \(y'=0\) plane. The dislocation’s Burgers vector, \(\vec{b}\), is then in the \(x'z'\)-plane, with edge component in the \(x'\)-direction and screw component in the \(z'\)-direction. When the dislocation slips, the dislocation line will move in the \(x'\)-direction.
For any such dislocation solution, there will be a shearing along the slip plane resulting in disregistry, i.e. a relative displacement between the top and bottom halves. This disregistry has limits such that it is \(0\) for \(x' \to -\infty\) and \(\vec{b}\) for \(x' \to +\infty\).
Within an atomic system, the dislocation line should be aligned with one of the system’s box vectors making the dislocation infinitely long and initially perfectly straight. The slip plane can then be defined as containing that box vector and another one. This results in the third box vector being the only one with a component parallel to the slip plane’s normal.
For LAMMPS-based simulations, the most convenient orientation to use is to align the dislocation with the \(\vec{a}\) box vector, and to define the slip plane as containing both \(\vec{a}\) and \(\vec{b}\). Given the limits that LAMMPS places on how system boxes can be defined, this results in favorable alignment of the system to the LAMMPS Cartesian system (\(x, y, z\)). The dislocation line will be along the \(x\)-axis, the slip plane normal parallel to the \(z\)-axis, and dislocation motion will be in the \(y\) direction. Thus, the LAMMPS coordinates corresponds to a rotation of the theory coordinates such that \(x'=y, y'=z, z'=x\).
Linear displacements solution
The simplest way to insert a dislocation is to cut the system in half along the slip plane and apply equal but opposite linear displacements, \(\vec{u}\), to the two halves with end conditions
\(\vec{u}(y=-\frac{Ly}{2}) = 0\)
\(\vec{u}(y=\frac{Ly}{2}) = \pm \frac{\vec{b}}{2}\)
Applying these displacements results in a disregistry along the slip plane that ranges from \(0\) to \(\vec{b}\). While the two \(y\) boundaries of the system both correspond to a perfect crystal, they are misaligned from each other by \(\frac{\vec{b}}{2}\). A coherent periodic boundary along the \(\vec{b}\) box vector can be established by adding or subtracting \(\frac{\vec{b}}{2}\) from \(\vec{b}\).
Note that with dislocations containing an edge component, a half-plane of atoms either needs to be inserted or removed to ensure boundary compatibility. Here, this is accomplished by always shifting \(\vec{b}\) to be shorter in the \(y\) direction, and removing excess atoms by identifying (near) duplicates.
Using dislocation solutions
A slightly more complicated, but ultimately more efficient, way of creating a periodic array of dislocations system is to combine the linear displacements solultion above with a more accurate linear elastic dislocation solution. The linear solution is used for the atoms at the free surfaces in the \(z\) direction, and for ensuring periodicity across the \(\vec{b}\) box vector direction. The linear elastic dislocation solution is then used for atoms in the middle of the system to construct an initial dislocation.
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. dislocation_array()
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 dislocation_array(lammps_command: str,
ucell: am.System,
potential: lmp.Potential,
C: am.ElasticConstants,
burgers: Union[list, np.ndarray],
ξ_uvw: Union[list, np.ndarray],
slip_hkl: Union[list, np.ndarray],
mpi_command: Optional[str] = None,
m: Union[list, np.ndarray] = [0,1,0],
n: Union[list, np.ndarray] = [0,0,1],
sizemults=None,
amin: float = None,
bmin: float = None,
cmin: float = None,
shift: Union[list, np.ndarray, None] = None,
shiftscale: bool = False,
shiftindex: int = None,
tol: float = 1e-8,
etol: float = 0.0,
ftol: float = 0.0,
maxiter: int = 10000,
maxeval: int = 100000,
dmax: float = uc.set_in_units(0.01, 'angstrom'),
annealtemp: float = 0.0,
annealsteps: Optional[int] = None,
randomseed: Optional[int] = None,
boundarywidth: float = 0.0,
boundaryscale: bool = False,
cutoff: Optional[float] = None,
linear: bool = False) -> dict:
"""
Creates and relaxes a dislocation monopole system.
Parameters
----------
lammps_command :str
Command for running LAMMPS.
ucell : atomman.System
The unit cell to use as the seed for generating the dislocation
monopole system.
potential : atomman.lammps.Potential
The LAMMPS implemented potential to use.
C : atomman.ElasticConstants
The elastic constants associated with the bulk crystal structure
for ucell.
burgers : array-like object
The dislocation's Burgers vector given as a Miller or
Miller-Bravais vector relative to ucell.
ξ_uvw : array-like object
The dislocation's line direction given as a Miller or
Miller-Bravais vector relative to ucell.
slip_hkl : array-like object
The dislocation's slip plane given as a Miller or Miller-Bravais
plane relative to ucell.
mpi_command : str or None, optional
The MPI command for running LAMMPS in parallel. If not given, LAMMPS
will run serially.
m : array-like object, optional
The m unit vector for the dislocation solution. m, n, and ξ
(dislocation line) should be right-hand orthogonal. Default value
is [0,1,0] (y-axis).
n : array-like object, optional
The n unit vector for the dislocation solution. m, n, and ξ
(dislocation line) should be right-hand orthogonal. Default value
is [0,0,1] (z-axis). n is normal to the dislocation slip plane.
sizemults : tuple, optional
The size multipliers to use when generating the system. Values are
limited to being positive integers. The multipliers for the two
non-periodic directions must be even. If not given, the default
multipliers will be 2 for the non-periodic directions and 1 for the
periodic direction.
amin : float, optional
A minimum thickness to use for the a box vector direction of the
final system. Default value is 0.0. For the non-periodic
directions, the resulting vector multiplier will be even. If both
amin and sizemults is given, then the larger multiplier for the two
will be used.
bmin : float, optional
A minimum thickness to use for the b box vector direction of the
final system. Default value is 0.0. For the non-periodic
directions, the resulting vector multiplier will be even. If both
bmin and sizemults is given, then the larger multiplier for the two
will be used.
cmin : float, optional
A minimum thickness to use for the c box vector direction of the
final system. Default value is 0.0. For the non-periodic
directions, the resulting vector multiplier will be even. If both
cmin and sizemults is given, then the larger multiplier for the two
will be used.
shift : float, optional
A rigid body shift to apply to the rotated cell prior to inserting
the dislocation. Should be selected such that the ideal slip plane
does not correspond to any atomic planes. Is taken as absolute if
shiftscale is False, or relative to the rotated cell's box vectors
if shiftscale is True. Cannot be given with shiftindex. If
neither shift nor shiftindex is given then shiftindex = 0 is used.
shiftindex : float, optional
The index of the identified optimum shifts based on the rotated
cell to use. Different values allow for the selection of different
atomic planes neighboring the slip plane. Note that shiftindex
values only apply shifts normal to the slip plane; best shifts for
non-planar dislocations (like bcc screw) may also need a shift in
the slip plane. Cannot be given with shiftindex. If neither shift
nor shiftindex is given then shiftindex = 0 is used.
shiftscale : bool, optional
If False (default), a given shift value will be taken as absolute
Cartesian. If True, a given shift will be taken relative to the
rotated cell's box vectors.
tol : float
A cutoff tolerance used with obtaining the dislocation solution.
Only needs to be changed if there are issues with obtaining a
solution.
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.
annealtemp : float, optional
The temperature to perform a dynamic relaxation at. Default is 0.0,
which will skip the dynamic relaxation.
annealsteps : int, optional
The number of time steps to run the dynamic relaxation for. Default
is None, which will run for 10000 steps if annealtemp is not 0.0.
randomseed : int or None, optional
Random number seed used by LAMMPS in creating velocities and with
the Langevin thermostat. Default is None which will select a
random int between 1 and 900000000.
boundarywidth : float, optional
The width of the boundary region to apply. Default value is 0.0,
i.e. no boundary region. All atoms in the boundary region will
have their atype values changed.
boundaryscale : bool, optional
If False (Default), the boundarywidth will be taken as absolute.
If True, the boundarywidth will be taken relative to the magnitude
of the unit cell's a box vector.
linear : bool, optional
If True, then only linear displacements will be used and not the
dislocation solution. Using only linear displacements is useful
for screw dislocations and dislocations with large stacking fault
distances. If False (default) then the dislocation solution will
be used for the middle displacements and linear displacements only
in the boundary region.
cutoff : float, optional
Cutoff distance to use for identifying duplicate atoms to remove.
For dislocations with an edge component, applying the displacements
creates an extra half-plane of atoms that will have (nearly)
identical positions with other atoms after altering the boundary
conditions. Default value is 0.5 Angstrom.
Returns
-------
dict
Dictionary of results consisting of keys:
- **'dumpfile_base'** (*str*) - The filename of the LAMMPS dump file
for the relaxed base system.
- **'symbols_base'** (*list of str*) - The list of element-model
symbols for the Potential that correspond to the base system's
atypes.
- **'dumpfile_disl'** (*str*) - The filename of the LAMMPS dump file
for the relaxed dislocation monopole system.
- **'symbols_disl'** (*list of str*) - The list of element-model
symbols for the Potential that correspond to the dislocation
monopole system's atypes.
- **'dislocation'** (*atomman.defect.Dislocation*) - The Dislocation
object used to generate the monopole system.
- **'E_total_disl'** (*float*) - The total potential energy of the
dislocation monopole system.
"""
# Construct dislocation configuration generator
dislocation = am.defect.Dislocation(ucell, C, burgers, ξ_uvw, slip_hkl,
m=m, n=n, shift=shift, shiftindex=shiftindex,
shiftscale=shiftscale, tol=tol)
# Generate the base and dislocation systems
base_system, disl_system = dislocation.periodicarray(sizemults=sizemults,
amin=amin, bmin=bmin, cmin=cmin,
shift=shift,
shiftindex=shiftindex,
shiftscale=shiftscale,
boundarywidth=boundarywidth,
boundaryscale=boundaryscale,
linear=linear,
cutoff=cutoff,
return_base_system=True)
# Initialize results dict
results_dict = {}
# Save initial perfect system
base_system.dump('atom_dump', f='base.dump')
results_dict['dumpfile_base'] = 'base.dump'
results_dict['symbols_base'] = base_system.symbols
# Save dislocation generator
results_dict['dislocation'] = dislocation
# Relax system
relaxed = disl_relax(lammps_command, disl_system, potential,
mpi_command=mpi_command, annealtemp=annealtemp,
annealsteps=annealsteps, randomseed=randomseed,
etol=etol, ftol=ftol, maxiter=maxiter,
maxeval=maxeval, dmax=dmax)
# Save relaxed dislocation system with original box vects
system_disl = am.load('atom_dump', relaxed['dumpfile'], symbols=disl_system.symbols)
system_disl.box_set(vects=disl_system.box.vects, origin=disl_system.box.origin)
system_disl.dump('atom_dump', f='disl.dump')
results_dict['dumpfile_disl'] = 'disl.dump'
results_dict['symbols_disl'] = system_disl.symbols
results_dict['E_total_disl'] = relaxed['E_total']
# Cleanup files
Path('0.dump').unlink()
Path(relaxed['dumpfile']).unlink()
for dumpjsonfile in Path('.').glob('*.dump.json'):
dumpjsonfile.unlink()
return results_dict
2.2. disl_relax()
[5]:
def disl_relax(lammps_command: str,
system: am.System,
potential: lmp.Potential,
mpi_command: Optional[str] = None,
annealtemp: float = 0.0,
annealsteps: Optional[int] = None,
randomseed: Optional[int] = None,
etol: float = 0.0,
ftol: float = 1e-6,
maxiter: int = 10000,
maxeval: int = 100000,
dmax: float = uc.set_in_units(0.01, 'angstrom')) -> dict:
"""
Sets up and runs the disl_relax.in LAMMPS script for relaxing a
dislocation monopole 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.
annealtemp : float, optional
The temperature to perform a dynamic relaxation at. Default is 0.0,
which will skip the dynamic relaxation.
annealsteps : int, optional
The number of time steps to run the dynamic relaxation for. Default
is None, which will run for 10000 steps if annealtemp is not 0.0.
randomseed : int or None, optional
Random number seed used by LAMMPS in creating velocities and with
the Langevin thermostat. Default is None which will select a
random int between 1 and 900000000.
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:
- **'logfile'** (*str*) - The name of the LAMMPS log file.
- **'dumpfile'** (*str*) - The name of the LAMMPS dump file
for the relaxed system.
- **'E_total'** (*float*) - The total potential energy for the
relaxed 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='system.dat',
potential=potential)
lammps_variables['atomman_system_pair_info'] = system_info
lammps_variables['anneal_info'] = anneal_info(annealtemp, annealsteps,
randomseed, potential.units)
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
lammps_variables['group_move'] = ' '.join(np.array(range(1, system.natypes // 2 + 1), dtype=str))
# Set dump_modify format based on dump_modify_version
if lammps_date < datetime.date(2016, 8, 3):
lammps_variables['dump_modify_format'] = '"%d %d %.13e %.13e %.13e %.13e"'
else:
lammps_variables['dump_modify_format'] = 'float %.13e'
# Write lammps input script
lammps_script = 'disl_relax.in'
template = read_calc_file('iprPy.calculation.dislocation_periodic_array',
'disl_relax.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)
thermo = output.simulations[-1]['thermo']
# Extract output values
results = {}
results['logfile'] = 'log.lammps'
results['dumpfile'] = '%i.dump' % thermo.Step.values[-1]
results['E_total'] = uc.set_in_units(thermo.PotEng.values[-1],
lammps_units['energy'])
return results
2.3. anneal_info()
[6]:
def anneal_info(temperature: float = 0.0,
runsteps: Optional[int] = None,
randomseed: Optional[int] = None,
units: str = 'metal') -> str:
"""
Generates LAMMPS commands for thermo anneal.
Parameters
----------
temperature : float, optional
The temperature to relax at (default is 0.0).
randomseed : int or None, optional
Random number seed used by LAMMPS in creating velocities and with
the Langevin thermostat. (Default is None which will select a
random int between 1 and 900000000.)
units : str, optional
The LAMMPS units style to use (default is 'metal').
Returns
-------
str
The generated LAMMPS input lines for performing a dynamic relax.
Will be '' if temperature==0.0.
"""
# Return nothing if temperature is 0.0 (don't do thermo anneal)
if temperature == 0.0:
return ''
# Generate velocity, fix nvt, and run LAMMPS command lines
else:
if randomseed is None:
randomseed = random.randint(1, 900000000)
if runsteps is None:
runsteps = 10000
start_temp = 2 * temperature
tdamp = 100 * lmp.style.timestep(units)
timestep = lmp.style.timestep(units)
info = '\n'.join([
'velocity move create %f %i mom yes rot yes dist gaussian' % (start_temp, randomseed),
'fix nvt all nvt temp %f %f %f' % (temperature, temperature,
tdamp),
'timestep %f' % (timestep),
'thermo %i' % (runsteps),
'run %i' % (runsteps),
])
return info
2.4. disl_relax.template file
[7]:
with open('disl_relax.template', 'w') as f:
f.write("""#LAMMPS input script for relaxing a dislocation
<atomman_system_pair_info>
group move type <group_move>
group hold subtract all move
compute peatom all pe/atom
dump first all custom <maxeval> *.dump id type x y z c_peatom
dump_modify first format <dump_modify_format>
thermo_style custom step pe
fix nomove hold setforce 0.0 0.0 0.0
<anneal_info>
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. Elastic constants
C_dict is a dictionary containing the unique elastic constants for the potential and crystal structure defined above.
C is an atomman.ElasticConstants object built from C_dict.
[11]:
C_dict = {}
C_dict['C11'] = uc.set_in_units(247.86, 'GPa')
C_dict['C12'] = uc.set_in_units(147.83, 'GPa')
C_dict['C44'] = uc.set_in_units(124.84, 'GPa')
# -------------- Derived parameters -------------- #
# Build ElasticConstants object from C_dict terms
C = am.ElasticConstants(**C_dict)
3.5. Defect parameters
burgers is the crystallographic Miller Burgers vector for the dislocation.
ξ_uvw is the Miller [uvw] line vector direction for the dislocation. The angle between burgers and ξ_uvw determines the dislocation’s character
slip_hkl is the Miller (hkl) slip plane for the dislocation.
m is the Cartesian vector of the final system that the dislocation solution’s m vector (in-plane, perpendicular to ξ) should align with. Limited to being parallel to one of the three Cartesian axes.
n is the Cartesian vector of the final system that the dislocation solution’s n vector (slip plane normal) should align with. Limited to being parallel to one of the three Cartesian axes.
shift is a rigid body shift to apply to the atoms in the system. This controls how the atomic positions align with the ideal position of the dislocation core, which is at coordinates (0,0) for the two Cartesian axes aligned with m and n.
shiftscale allows for shift to be defined relative to the cell created by rotating ucell to coincide with the dislocation solution orientation. This is useful as it allows for shift values to be defined relative to the defect type and crystal prototype rather than on a per-crystal basis.
shiftindex alternate to specifying shift values, the shiftindex allows for one of the identified suggested shift values to be used that will position the slip plane halfway between two planes of atoms. Note that shiftindex values only shift atoms in the slip plane normal direction and may not be the ideal positions for some dislocation cores.
[12]:
# fcc a/2 <110>{111} dislocations
burgers = 0.5 * np.array([ 1., -1., 0.])
slip_hkl = np.array([1, 1, 1])
# Line direction determines dislocation character
ξ_uvw = np.array([ 1, 1,-2]) # 90 degree edge
#ξ_uvw = np.array([ 1, 0,-1]) # 60 degree mixed
#ξ_uvw = np.array([ 1,-2, 1]) # 30 degree mixed
#ξ_uvw = np.array([ 1,-1, 0]) # 0 degree screw
# Best choice for m + n as it works for non-cubic systems
m = [0,1,0]
n = [0,0,1]
# Specify shift or shiftindex
shift = None
shiftscale = True
shiftindex = 0
3.6. Calculation-specific parameters
boundarywidth sets the minimum width of the fixed-atom boundary region.
boundaryscale flag indicating if boundarywidth is absolute (False) or scaled relative to ucell’s a lattice parameter (True).
duplicatecutoff Cutoff distance to use for identifying duplicate atoms to remove. For dislocations with an edge component, applying the displacements creates an extra half-plane of atoms that will have (nearly) identical positions with other atoms after altering the boundary conditions.
onlylinear boolean flag. If False (default) the atoms in the middle of the system will be displaced according to an elasticity solution for a dislocation core and boundary atoms displaced by a linear gradient. If True, all atoms are displaced by the linear gradient which may be necessary for pure screw dislocations.
annealtemperature is the temperature at which to relax (anneal) the dislocation system. If annealtemperature is 0.0, then only a static relaxation will be performed. Default value is 0.0.
annealsteps is the number of nvt iteration steps to perform at the given temperature. Default value is 0 if annealtemperature is zero, 10000 otherwise.
randomseed allows for the random seed used in generating initial atomic velocities for a dynamic relaxation to be specified. This is an integer between 1 and 900000000. Default value is None, which will randomly pick a number in that range.
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]:
boundarywidth = 3
boundaryscale = True
duplicatecutoff = uc.set_in_units(0.5, 'angstrom')
onlylinear = False
annealtemperature = 50.0
annealsteps = 10000
randomseed = None
energytolerance = 0.0
forcetolerance = uc.set_in_units(1e-6, 'eV/angstrom')
maxiterations = 10000
maxevaluations = 100000
maxatommotion = uc.set_in_units(0.01, 'angstrom')
3.7. System modifications
sizemults list of three sets of two integers specifying how many times the ucell vectors of \(a\), \(b\) and \(c\) are replicated in positive and negative directions when creating system.
amin, bmin, cmin specifies minimum widths in the three directions. The corresponding sizemults will be increased to ensure these widths if needed.
[14]:
sizemults = [1, 40, 40]
amin = uc.set_in_units(0.0, 'angstrom')
bmin = uc.set_in_units(0.0, 'angstrom')
cmin = uc.set_in_units(0.0, '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.
[15]:
results_dict = dislocation_array(lammps_command, ucell, potential, C,
burgers, ξ_uvw, slip_hkl,
mpi_command = mpi_command,
m = m,
n = n,
shift = shift,
shiftscale = shiftscale,
shiftindex = shiftindex,
sizemults = sizemults,
amin = amin,
bmin = bmin,
cmin = cmin,
etol = energytolerance,
ftol = forcetolerance,
maxiter = maxiterations,
maxeval = maxevaluations,
dmax = maxatommotion,
annealtemp = annealtemperature,
annealsteps = annealsteps,
randomseed = randomseed,
boundarywidth = boundarywidth,
boundaryscale = boundaryscale,
cutoff = duplicatecutoff,
linear = onlylinear)
print(results_dict.keys())
dict_keys(['dumpfile_base', 'symbols_base', 'dislocation', 'dumpfile_disl', 'symbols_disl', 'E_total_disl'])
4.2. Report results
Values returned in the results_dict:
‘dumpfile_base’ (str) - The filename of the LAMMPS dump file for the relaxed base system.
‘symbols_base’ (list of str) - The list of element-model symbols for the Potential that correspond to the base system’s atypes.
‘dumpfile_disl’ (str) - The filename of the LAMMPS dump file for the relaxed dislocation monopole system.
‘symbols_disl’ (list of str) - The list of element-model symbols for the Potential that correspond to the dislocation monopole system’s atypes.
‘dislocation’ (atomman.defect.Dislocation) - The Dislocation object used to generate the monopole system.
‘E_total_disl’ (float) - The total potential energy of the dislocation monopole system.
[16]:
length_unit = 'angstrom'
energy_unit = 'eV'
pressure_unit = 'GPa'
energy_per_length_unit = energy_unit+'/'+length_unit
[17]:
print('pre-ln factor (alpha) =', uc.get_in_units(results_dict['dislocation'].dislsol.preln, energy_per_length_unit), energy_per_length_unit)
print('K_tensor =')
print(uc.get_in_units(results_dict['dislocation'].dislsol.K_tensor, pressure_unit), pressure_unit)
pre-ln factor (alpha) = 0.38115202169620926 eV/angstrom
K_tensor =
[[ 81.07544348 0. -12.76027037]
[ 0. 123.86918865 0. ]
[-12.76027037 0. 127.31323026]] GPa
[18]:
print('Perfect system is saved as', results_dict['dumpfile_base'])
print('with symbols', results_dict['symbols_base'])
print('Defect system is saved as', results_dict['dumpfile_disl'])
print('with symbols', results_dict['symbols_disl'])
Perfect system is saved as base.dump
with symbols ('Ni',)
Defect system is saved as disl.dump
with symbols ('Ni', 'Ni')