Sentaurus Materials Workbench Technical Descriptions

Operating with defect lists

Introduction

The defect list commands define the types of defect and set up the environment to calculate the energies and entropies of defects. Each of these commands (see Defect list types) indicates the defect type.

Defect lists usage

Once a particular type of defect has been chosen by the use of the particular defect list types, the usage is very similar.

  • The defect list finds all the possible defects in the specified reference MaterialSpecifications and finds out which ones are equivalent for the given reference. Then, it creates a list on non equivalent defects.

  • By default, all non equivalent defects are included. A different approach can be followed by filtering the defects. The available filters are listed for each particular defect, being ‘filterByWyckoff’, ‘filterByLatticeSpecies’, and ‘filterByDistance’.

    • ‘filterByWyckoff’ is intended to let users pick up particular defects in the defect list.

    • ‘filterByLatticeSpecies’ is intended to select inequivalent defects by element name

    • ‘filterByDistance’ is used to pick up defects at a particular distance from an origin.

  • Once the defects have been chosen and optionally filtered, the energy and entropy calculations are run by calling the ‘update’ command on them.

  • The entropy calculations are optional, and can be disabled by defining and assumed formation entropy in MaterialSpecifications. By default they are performed.

  • Once the defects have been updated they can be saved with a regular nlsave command. These saved files can be open and visualized.

  • A simple, text-based, visualization of the computed values can be obtained with the displayDefectCharacterization command.

  • Further preprocessing on the values computed to generate a Sentaurus Process file can be done using the writeSentaurusParameters command.

Defect list types

The defect list types are:

Directory names

All the defect list automatically create directories with the logs and results of the ChargedPointDefect calculations that are run for each particular defect. Such directories names have the following conventions.

  1. Defect name, in particular cluster, interstitial, pair, split, substitutional and vacancy

  2. Defect description, usually including the defect Wyckoff type. See each particular defect for the description.

  3. Defect site, different for each defect, even if they are equivalent.

  4. Material specifications hash, that uniquely identify the particular material specification options used as reference.

For instance, the following shows directory names created by SMW:

pair/CV/C^Si_0+Si_0/001+000/8510ca94d3b960045791f06a3741606f
vacancy/Si_0/004/8510ca94d3b960045791f06a3741606f
substitutional/C^Si_0/004/8510ca94d3b960045791f06a3741606f
interstitial/C_0/004/8510ca94d3b960045791f06a3741606f
cluster/AsV4/Si_0+Si_0+Si_0+Si_0+As^Si_0/004+002+015+008/8510ca94d3b960045791f06a3741606f

Usage example

This example script shows how to perform a simple study of vacancies and interstitials in bulk silicion, displaying the results on screen and extracting a Sentaurus Process KMC file.

from SMW import *

# Setup the reference bulk unit cell configuration.
lattice = FaceCenteredCubic(5.4306*Angstrom)

elements = [Silicon, Silicon]

fractional_coordinates = [[ 0.  ,  0.  ,  0.  ],
                          [ 0.25,  0.25,  0.25]]

bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates,
)

# Setup the calculator.
basis_set = getattr(LDABasis, 'Silicon_SingleZeta')
calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=LDA.PZ,
    numerical_accuracy_parameters=NumericalAccuracyParameters(
        k_point_sampling=MonkhorstPackGrid(1, 1, 1))
)

# Create a MaterialSpecifications out out it
# An easier alternative is to use an already specified MaterialsSpecifications.
# For this example, we trade off accuracy for speed and create a small cell
# with 1x1x1 k-points and LDA instead of calling a more accurate specification
# from the database.
reference = MaterialSpecifications(
    pristine_configuration=bulk_configuration,
    formation_energy_calculator=calculator,
    supercell_repetitions=(1, 1, 1),
    sentaurus_material_name="Silicon",
)

# Create a list with all inequivalent vacancies and investigate charges 0 and 1
vacancies = VacancyList(reference, [0, 1])
nlprint(vacancies)

# Create a list with all non equivalent interstitials, charges [-1, 0]
interstitials = InterstitialList(reference, Silicon, [-1, 0])
nlprint(interstitials)

# Update, i.e., perform calculations
vacancies.update()
interstitials.update()

# Save
nlsave('defects.hdf5', vacancies)
nlsave('defects.hdf5', interstitials)

# Show results on screen
displayDefectCharacterization([vacancies, interstitials], [])

# Generate a simple Sentaurus Process KMC file with them.
writeSentaurusParameters([vacancies, interstitials], [], SentaurusProcessKMC, 'KMC.pdb')

defect_lists.py

It is important to note that this example is only intended to provide an overview of the creation of defect lists, and its use. In particular, it does not extensively study the defects (by adding other charges, and more defects, e.g., split interstials and substitutionals), and it uses a fast calculator.

nlprint output

nlprint detects only 1 inequivalent position for the vacancy, and 2 inequivalent positions for the interstitial (the tetrahedral and hexagonal positions).

+------------------------------------------------------------------------------+
| VacancyList defect summary                                                   |
+------------------------------------------------------------------------------+
| 1 inequivalent defect positions:                                             |
|   * Wyckoff 0 (2 defects) is Si_0, site 001:                                 |
|     - Si vacancy of type 0, site 1 at (0.2500, 0.2500, 0.2500)               |
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| InterstitialList defect summary                                              |
+------------------------------------------------------------------------------+
| 2 inequivalent defect positions:                                             |
|   * Wyckoff 0 (4 defects) is Si_0, site 004:                                 |
|     - Si interstitial of type 0 at (0.6250, 0.6250, 0.6250)                  |
|   * Wyckoff 1 (2 defects) is Si_1, site 002:                                 |
|     - Si interstitial of type 1 at (0.5000, 0.5000, 0.5000)                  |
+------------------------------------------------------------------------------+

Calculation update

The calculation starts with an update of the vacancy and interstitial lists. For example, the vacancy list calculations are reported as:

+------------------------------------------------------------------------------+
| VacancyList                                                                  |
+------------------------------------------------------------------------------+
| 1 defect(s) will be calculated.                                              |
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| - Reloading 0 tasks for vacancy/Si_0/001.                                    |
| - Tasks to be done (14):                                                     |
|   * : Pristine phonons calculation. Repetitions: 1 1 1                       |
|   * : Update pristine configuration. k-points: 1 1 1                         |
|   * : Update defect configuration for phonon calculations. Charge: 0         |
|       Repetitions: 1 1 1 Relax coordinates: True                             |
|   * : Update defect. Charge: 0 Repetitions: 1 1 1 Relax coordinates: True    |
|   * : Update pristine configuration. k-points: 1 1 1                         |
|   * : Total energy of pristine configuration. k-points: 1 1 1                |
|   * : Valence band maximum. k-points: 1 1 1                                  |
|   * : Update defect. Charge: 0 Repetitions: 1 1 1 Relax coordinates: True    |
|   * : Total energy of defect. Charge: 0 Repetitions: 1 1 1                   |
|   * : Band shift correction. Charge: 0 Repetitions: 1 1 1                    |
|   * : Update defect. Charge: 1 Repetitions: 1 1 1 Relax coordinates: True    |
|   * : Total energy of defect. Charge: 1 Repetitions: 1 1 1                   |
|   * : Band shift correction. Charge: 1 Repetitions: 1 1 1                    |
|   * : Periodic charge correction. Charge: 1 Repetitions: 1 1 1               |
| Folder: vacancy/Si_0/001/69ac938a14462d7dcd9be9da8da30356/                   |
+------------------------------------------------------------------------------+
| VacancyList. Simulations: 0/1. Tasks: 0/14. Cores per item 1.                |
+------------------------------------------------------------------------------+
| Name                                               Tasks      Status Cores   |
| vacancy/Si_0/001                                    0/14     Running     1   |
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| VacancyList. Simulations: 1/1. Tasks: 14/14. Cores per item 1.               |
+------------------------------------------------------------------------------+
| Name                                               Tasks      Status Cores   |
| vacancy/Si_0/001                                   14/14   Converged     1   |
+------------------------------------------------------------------------------+

Complete information about the type and number of tasks is given, together where the folder with the final results are located. The message about convergence status informs us about the correct finishing of the simulation, plus the final status for the simulation: all results are converged.

Calculation results

Once all the calculation results are available, the displayDefectCharacterization function provides an overview of the available results.

# Show results on screen
displayDefectCharacterization([vacancies, interstitials], [])

producing:

Legend                                                                 Si
                                              ---------------------------

Total energies (eV)                           ---------------------------
Si (2)                                                         -337.019
vacancy/Si_0/001                                               -170.119
vacancy/Si_0/001.p1                                            -170.216
interstitial/Si_0/004.n1                                       -517.566
interstitial/Si_0/004                                          -514.321
interstitial/Si_1/002.n1                                       -517.566
interstitial/Si_1/002                                          -514.321

Given chemical potentials (eV)                ---------------------------
Si                                                             -168.510

Element chemical potentials (eV)                                 NOTYPE

Formation energies (eV)                       ---------------------------
vacancy/Si_0/001.p1                                              -3.656
vacancy/Si_0/001                                                 -1.692

interstitial/Si_1/002                                            -8.735
interstitial/Si_0/004                                            -8.735
interstitial/Si_1/002.n1                                         -8.419
interstitial/Si_0/004.n1                                         -3.633

Fermi-level dependent formation energies (eV) ---------------------------
Eg (Extracted)                                                    2.443
Eg (Specified)                                                     None

vacancy/Si_0/001                            q  Ef@Ev  Ef@Ei  Ef@Ec  Etrap
                                              ---------------------------
                                         ( 0) -1.692 -1.692 -1.692      -
                                         ( 1) -3.656 -2.435 -1.213  0.479

                                       Trans.       EF-Evbm            Ef
                                              ---------------------------
                                      ( 1/ 1)         0.000        -3.656
                                      ( 1/ 0)         1.964        -1.692
                                      ( 0/ 0)         2.443        -1.692

interstitial/Si_0/004                       q  Ef@Ev  Ef@Ei  Ef@Ec  Etrap
                                              ---------------------------
                                         (-1) -3.633 -4.854 -6.076      -
                                         ( 0) -8.735 -8.735 -8.735 -2.659

                                       Trans.       EF-Evbm            Ef
                                              ---------------------------
                                      ( 0/ 0)         0.000        -8.735

interstitial/Si_1/002                       q  Ef@Ev  Ef@Ei  Ef@Ec  Etrap
                                              ---------------------------
                                         (-1) -8.419 -9.641 -10.862      -
                                         ( 0) -8.735 -8.735 -8.735  2.128

                                       Trans.       EF-Evbm            Ef
                                              ---------------------------
                                      ( 0/ 0)         0.000        -8.735
                                      ( 0/-1)         0.315        -8.735
                                      (-1/-1)         2.443       -10.862

Sentaurus Process Parameters

The writeSentaurusParameters function creates the following file:

# Generated by Sentaurus Materials Workbench
pdbSet Silicon KMC Int Ef   -8.73457
pdbSet Silicon KMC Int D0FS 9.53815
pdbSet Silicon KMC Vac Ef   -1.69186
pdbSet Silicon KMC Vac D0FS 0.131049
pdbSetDoubleArray Silicon KMC Int e0 {
                  IM 0.315462
}
pdbSetDoubleArray Silicon KMC Vac e0 {
                  VP 1.96434
}

KMC.pdb

Transition Path calculation

For transition path calculations, needed to compute migration barriers, see the TransitionPathList page.

introbar

Chemical Potentials in Compound Materials

The total energy of a compound material is always smaller than the sum of the chemical potentials of the elemental solids. Otherwise, the compound material cannot exist in a stable state. For example, the total energy of the 64-atom TiN supercell, which is composed of 32 Ti and 32 N atoms, is smaller than 32 times the sum of the Ti chemical potential in the Ti elemental solid and the N chemical potential in the \(N_2\) gas phase.

\[E_{TiN} = \mu_{Ti} + \mu_N + dH\]

where \(\mu_{Ti}\) and \(\mu_N\) are elemental chemical potential of Ti and N. \(E_{TiN}\) is the total energy of one TiN pair which is calculated by (total energy of 64-atom TiN)/32 in reference system. In defect formation energy calculation, it is needed to reflect \(dH\) (enthalpy of formation) to elemental chemical potential so that the following equation can be satisfied

\[E_{TiN} = \mu_{Ti}' + \mu_{N}'\]

Sentaurus Materials Workbench provides different methods to determine the correct chemical potentials of elements in a compound material:

  • Heat of formation

  • Specific element-rich condition

  • Grand thermodynamic potential

Heat of Formation

The heat of formation method assumes the equal enthalpy change for each element atom from the elemental solid to the compound material.

\[\mu_{Ti}' = \mu_{Ti} + \frac{dH}{2}\]
\[\mu_{N}' = \mu_N + \frac{dH}{2}\]

Specific Element-Rich Condition

The rich condition assumes that the compound forms in a specific reservoir, that is, an elemental solid. In that case, the chemical potentials of the elements in the compound can be set to those from the reservoir systems. For example, the Ti chemical potential for TiN is set to the chemical potential from the Ti elemental solid, which is called the Ti-rich condition.

\[\mu_{Ti}' = \mu_{Ti}\]
(26)\[\mu'_N = \mu_N + dH = E_{TiN} - \mu_{Ti}\]

Heat or formation and Specific Element-Rich condition can be selected by AtomicChemicalPotentialList which takes MaterialSpecifications and AtomicChemicalPotential of each element as arguments.

If the AtomicChemicalPotential of all elements is defined, SMW will calculate chemical potentials for Heat or formation condition while, if one element is defined, SMW assumes Specific Element-Rich condition for the element.

The script flow is

  1. Define reference (MaterialSpecifications) imported from database or defined by user.

  2. Calculate either heat of formation or specific element-rich condition in AtomicChemicalPotentialList

  3. Update calculated chemical potential in reference (MaterialSpecifications) system.

  4. Calculate defects with updated reference (MaterialSpecifications).

Example script for Heat of Formation in TiN

from SMW import *
filename="Reference/TiN.hdf5"
bulk_configuration=nlread(filename, BulkConfiguration)[-1]
calculator=bulk_configuration.calculator()

# Define the reference system- TiN
TiN = MaterialSpecifications(
    pristine_configuration=bulk_configuration,
    supercell_repetitions=(1, 1, 1),
    formation_energy_calculator=calculator(),
    relaxation_calculator=calculator(),
    atomic_chemical_potentials=None,
    bulk_modulus= 288.0 * GPa,
    assumed_formation_entropy= 0 * boltzmann_constant,
    sentaurus_material_name='TitaniumNitride',
    optimize_geometry_parameters=OptimizeGeometryParameters(max_forces=0.01*eV/Ang)

)

################################################################################################
# Chemical potential set-up: Heat of Formation
################################################################################################

Heat_of_formation   = AtomicChemicalPotentialList(TiN,
                                   [AtomicChemicalPotential(Titanium),
                                    AtomicChemicalPotential(Nitrogen)])

# update atomic chemical potential based on the information in AtomicChemicalPotential class
Heat_of_formation.update()
nlprint(Heat_of_formation)

# update chemical potentials in reference system (MaterialSpecifications) with Heat of formation option
TiN_HOF = TiN(atomic_chemical_potentials=Heat_of_formation.calculatedAtomicChemicalPotentials())

################################################################################################
# Defect definition
################################################################################################
# Ti and N vacancies in TiN
vacancy = VacancyList(TiN_HOF, processes_per_defect=14)
vacancy.update()
nlsave('defects.hdf5', vacancy)

displayDefectCharacterization([vacancy],[])

heat_of_formation.py

In the example script,

Heat_of_formation.update()
nlprint(Heat_of_formation)

will calculate and print new chemical potentials for Ti and N under Heat of Formation condition. The output will be as follows:

+------------------------------------------------------------------------------+
| AtomicChemicalPotentialList summary                                          |
+------------------------------------------------------------------------------+
| 2 elements provided:                                                         |
| Calculated total pristine energy: -14501.81 eV                               |
+------------------------------------------------------------------------------+
| Elemental energy, entropy: mu_el, s_el. Chemical potential, entropy: mu, s   |
| Ti: mu_el = -179.24 eV, s_el= None, mu = -181.07 eV, s = None.               |
| N: mu_el = -270.27 eV, s_el= None, mu = -272.11 eV, s = None.                |
+------------------------------------------------------------------------------+
| Type: Heat of formation                                                      |
+------------------------------------------------------------------------------+

Example script for Ti rich-condition in TiN

The script flow is the same as one for heat of formation but, now only AtomicChemicalPotential(Titanium) is defined in AtomicChemicalPotentialList.

SMW calculates the chemical potential of Ti from Ti elemental solid and the chemical potential of N is calculated by equation (26)

Ti_rich   = AtomicChemicalPotentialList(TiN,[AtomicChemicalPotential(Titanium)])
Ti_rich.update()
nlprint(Ti_rich)

# update reference material with updated chemical potential - Ti rich condition
TiN_rich = TiN(atomic_chemical_potentials=Ti_rich.calculatedAtomicChemicalPotentials())

Ti_rich.py

which output will be as follows:

+------------------------------------------------------------------------------+
| AtomicChemicalPotentialList summary                                          |
+------------------------------------------------------------------------------+
| 1 elements provided:                                                         |
| Calculated total pristine energy: -14501.81 eV                               |
+------------------------------------------------------------------------------+
| Elemental energy, entropy: mu_el, s_el. Chemical potential, entropy: mu, s   |
| Ti: mu_el = -179.24 eV, s_el= None, mu = -179.24 eV, s = None.               |
| N: mu_el = None, s_el= None, mu = -273.95 eV, s = None.                      |
+------------------------------------------------------------------------------+
| Type: Rich Ti                                                                |
+------------------------------------------------------------------------------+

Grand Thermodynamic Potential (GTP)

Sentaurus Materials Workbench uses a self-consistent iterative method to calculate the chemical potentials with (m-1) stoichiometric balance equations.

GTP requires defect energies for vacancies, anti-site substitutionals and interstitials for all elements in the reference material.

The script flow is:

  1. Define reference (MaterialSpecifications) imported from database or defined by user.

  2. Calculate defects

  3. Calculate GTP chemical potential using grandThermodynamicPotential

  4. Update calculated chemical potential in reference (MaterialSpecifications) system.

  5. Calculate defects with updated reference (MaterialSpecifications).

All tasks in defect calculations have been already done in the step 2. Thus, in step 5, SMW will only proceed the update in formation energy calculation using GTP chemical potentials.

Example script for TiN defect calculation using GTP:

from SMW import * 

filename="Reference/TiN.hdf5"

bulk_configuration=nlread(filename, BulkConfiguration)[-1]
calculator=bulk_configuration.calculator()


# Define the reference material - TiN
TiN = MaterialSpecifications(
    pristine_configuration=bulk_configuration,
    supercell_repetitions=(1, 1, 1),
    formation_energy_calculator=calculator(),
    relaxation_calculator=calculator(),
    atomic_chemical_potentials=None,
    bulk_modulus= 288.0 * GPa,
    assumed_formation_entropy= 0 * boltzmann_constant,
    sentaurus_material_name='TitaniumNitride',
    optimize_geometry_parameters=OptimizeGeometryParameters(max_forces=0.01*eV/Ang)

)


# Ti and N vacancies in TiN
vacancy = VacancyList(TiN)
nlsave('defects.hdf5', vacancy)

# Anti-sites. Ti substitutional on N site
Ti_N = SubstitutionalList(TiN, Titanium)
nlsave('defects.hdf5', Ti_N)

# Anti-sites. N substitutional on Ti site
N_Ti  = SubstitutionalList(TiN, Nitrogen)
nlsave('defects.hdf5', N_Ti)

# Ti interstitials
Ti_interstitial = InterstitialList(TiN, Titanium)
nlsave('defects.hdf5', Ti_interstitial)

# N interstitials
N_interstitial = InterstitialList(TiN, Nitrogen)
nlsave('defects.hdf5', N_interstitial)

# print defects info
if processIsMaster():
    nlprint(vacancy)
    nlprint(Ti_N)
    nlprint(N_Ti)
    nlprint(Ti_interstitial)
    nlprint(N_interstitial)

# update and run all defects
defect_lists=[vacancy, Ti_N, N_Ti, Ti_interstitial, N_interstitial]

updateAllDefectsAndTransitions(defect_lists=defect_lists, processes_per_defect=14)

##########################################################################################
# Update results with GTP
##########################################################################################

# calculate atomic chemical potentials using GTP
GTP=grandThermodynamicPotential(defect_lists)

# print GTP chemical potentials
if processIsMaster():
    for x in GTP:
        print(x.element().name(), x.internalEnergy(), x.entropy())

# Update GTP potentials in reference system
TiN_GTP=TiN(atomic_chemical_potentials=GTP)

# Re-calculate energies based on GTP

# Ti and N vacancies in TiN  
vacancy = VacancyList(TiN_GTP)
nlsave('defects.hdf5', vacancy)

# Anti-sites. Ti substitutional on N site
Ti_N = SubstitutionalList(TiN_GTP, Titanium)
nlsave('defects.hdf5', Ti_N)

# Anti-sites. N substitutional on Ti site
N_Ti  = SubstitutionalList(TiN_GTP, Nitrogen)
nlsave('defects.hdf5', N_Ti)

# Ti interstitials
Ti_interstitial = InterstitialList(TiN_GTP, Titanium)
nlsave('defects.hdf5', Ti_interstitial)

# N interstitials
N_interstitial = InterstitialList(TiN_GTP, Nitrogen)
nlsave('defects.hdf5', N_interstitial)

# update and run all defects 
defect_lists=[vacancy, Ti_N, N_Ti, Ti_interstitial, N_interstitial]

updateAllDefectsAndTransitions(defect_lists=defect_lists, processes_per_defect=14)

GTP.py

That produces:

Nitrogen -271.794577432 eV None
Titanium -181.387067179 eV None

Summarized chemical potentials in the TiN example

Chemical Potential

\(\mu_{Ti}\) (eV)

\(\mu_N\) (eV)

Elemental

-179.24

-270.27

Heat of formation

-181.07

-272.11

Ti rich-condition

-179.24

-273.95

N rich-condition

-182.91

-270.27

GTP

-181.39

-271.79

Note

The tutorial is generated to explain different chemical potentials options and how to set them up in the script. Elemental chemical potentials in this example are calculated from elemental structures from database which are not relaxed by pseudopotential and basis sets used for the TiN reference. Thus, elemental chemical potentials could be change once they are fully relaxed.

introbar

Diffusivity calculation in amorphous materials

Diffusivity calculations in amorphous materials can be computed using different procedures. The procedure described here can be defined as the “frozen amorphous approximation”. It assumes that, given a particular amorphous configuration, one can detect the positions where the impurities will be as interstitials in the structures, and simulate the migration paths between them. The assumption is that the amorphous will not change significatively during the migration of the impurities, so that the available paths can be determined and computed in advanced in a way similar to a crystalline structure without need to run dynamic simulations.

Alternative approaches using molecular dynamics can be used, although they will tend to be more computationally intensive. In this case one will proceed in a way similar to the computation of diffusivity in a liquid material, as explained in Diffusion in Liquids from Molecular Dynamics Simulations.

Procedure

Calculation of diffusivity of impurities (i.e., non native atoms) in amorphous takes the following steps

  1. Creation of a suitable amorphous material using amorphize.

  2. Introducing defects at different positions in the material. An approach is to add the impurity as an “interstitial” in the configuration previously created using InterstitialList. An alternative approach would be to create them as “substitutional” defects using SubstitutionalList. Substitutional and interstitial concepts might not be well defined for amorphous materials, but the goal of this approach is to have a distribution of probable impurity positions. Also, the “frozen amorphous approximation” treats the amorphous material as a crystal with no symmetries.

  3. Optionally filtering the created defects. Since the amorphous structure contains no symmetries all defects will be considered non equivalent. This might create a big number of possible configurations, involving very comutationally expensive calculations, especially when considering that all the transitions between neighboring configurations are also to be computed. A way to reduce the number of configurations is to look for configurations similar enough and approximate them as equivalent. This is done with the filterByEquivalentEnvironments filter. (See SentaurusMaterialsWorkbench.DefectList.InterstitialList.InterstitialList.filterByEquivalentEnvironments() or SentaurusMaterialsWorkbench.DefectList.SubstitutionalList.SubstitutionalList.filterByEquivalentEnvironments()). The function provides full control on how many final inequivalent defects are desired, assigning the rest of them to the most similar inequivalent ones.

  4. Computing migration paths between the selected defects using TransitionPathList.

  5. Extracting diffusivity with the monteCarloDiffusivityCalculation function by setting the amorphous_material option to ‘True’.

Example: Carbon diffusivity in amorphous silicon

Note

The following example is used for illustration purposes only and contains several assumptions (use of force fiels, very small system size, use of 1 only realization of the amorphous configuration, avoiding phonon calculations, etc.), despite the “frozen amorphous approximation” one and the assumption of “interstitial” positions for the impurities there, that prevents it from providing a realistic measurement of C diffusion in \(\alpha\)-Si.

Creation of amorphous material

Once a suitable crystalline reference for silicon has been defined, it can be amorphized with the command:

method = MolecularDynamicsMeltQuench(amorphous_quality=0.75)
amorphous = amorphize(reference, method, 'amorphous-Si.hdf5', 'amorphous')

create-amorphous.py

Producing something similar to:

+------------------------------------------------------------------------------+
| Density: 2.32958070454 g/cm**3                                               |
| Amorphization assessment: 0.811735950565 [0-1]                               |
| Coordination quality after geometry relaxation: 0.796875 [0-1]               |
| Probably a valid, but defected, amorphous.                                   |
+------------------------------------------------------------------------------+

Relaxation of the amorphous material

The amorphous creation fixes the amorphous density to the requested one. An extra step ensuring that the system relaxes to the correct density is required.

# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------
bulk_configuration = nlread("amorphous-Si-FF.hdf5", BulkConfiguration)[-1]

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------

constraints = [BravaisLatticeConstraint()]

bulk_configuration = OptimizeGeometry(
    bulk_configuration,
    max_forces=0.05*eV/Ang,
    max_stress=0.1*GPa,
    max_steps=200,
    max_step_length=0.2*Ang,
    constraints=constraints,
    trajectory_filename="trajectory.hdf5",
    optimizer_method=LBFGS(),
)
nlsave('relaxed.hdf5', bulk_configuration)
nlprint(bulk_configuration)

relax-amorphous.py

The optimization slighly decreases the total volume after 18 steps:

+------------------------------------------------------------------------------+
|       Step Step Length      Volume  Max. Force  Stress Error        Enthalpy |
|                  (Ang)    (Ang**3)    (eV/Ang)         (GPa)            (eV) |
+------------------------------------------------------------------------------+
| OPT      0  0.0000e+00    1281.249  4.0362e-02       2.87314     -274.484506 |
| OPT      1  2.0000e-01    1265.188  6.7261e-01       0.73503     -274.686882 |
| OPT      2  2.1925e-02    1265.290  8.3871e-01       1.37784     -274.697253 |
| OPT      3  7.7712e-03    1265.095  4.6161e-01       1.26136     -274.724506 |
| OPT      4  1.5944e-02    1264.650  6.8857e-01       0.99890     -274.743216 |
| OPT      5  1.5221e-02    1264.377  4.0008e-01       1.02232     -274.762084 |
| OPT      6  1.0119e-01    1262.518  4.8816e-01       1.05634     -274.829821 |
| OPT      7  1.0991e-01    1260.244  9.7923e-01       0.56901     -274.864684 |
| OPT      8  4.0849e-02    1259.524  6.5189e-01       0.26991     -274.900963 |
| OPT      9  2.6241e-02    1259.162  2.4004e-01       0.12602     -274.917765 |
| OPT     10  1.9388e-02    1259.009  3.3161e-01       0.25610     -274.924518 |
| OPT     11  9.5720e-03    1259.092  1.8791e-01       0.27491     -274.928307 |
| OPT     12  1.1892e-02    1259.691  4.3707e-01       0.22745     -274.930769 |
| OPT     13  7.7098e-03    1260.053  1.2431e-01       0.20002     -274.933578 |
| OPT     14  3.6644e-03    1260.482  1.1018e-01       0.14468     -274.935018 |
| OPT     15  1.2022e-02    1261.625  9.6900e-02       0.06681     -274.937680 |
| OPT     16  9.3826e-03    1262.513  6.1651e-02       0.04835     -274.938993 |
| OPT     17  6.7729e-03    1263.010  6.7651e-02       0.07240     -274.939766 |
| OPT     18  1.9191e-03    1262.978  4.8074e-02       0.08901     -274.940150 |
+------------------------------------------------------------------------------+
| Geometry optimization converged in 18 steps.                                 |
+------------------------------------------------------------------------------+

Impurity inclusion and filtering

Impurities can be created by the simple use of an InterstitialList and the selection of a fixed number of them to reduce the number of inequivalent positions. In this example a very limited number of only 3 different positions is assumed. An existing material configuration is modified to adapt to this simple example.

from SMW import *
from SentaurusMaterialsWorkbench.MaterialSpecifications.DefectMaterialCollection.Silicon_216_DFT_SCAN import Silicon_216_DFT_SCAN

amorphous = nlread("relaxed.hdf5", BulkConfiguration)[-1]

specs = Silicon_216_DFT_SCAN(
    formation_energy_calculator=TremoloXCalculator(parameters=Tersoff_SiC_2005()),
    pristine_configuration=amorphous,
    supercell_repetitions=(1, 1, 1),
    assumed_formation_entropy=5*boltzmann_constant,
    assumed_transition_prefactor=1e13/Second,
)

ints = InterstitialList(
    material_specifications=specs,
    interstitial_element=Carbon,
    )

filtered = ints.filterByEquivalentEnvironments(3)
nlprint(filtered)
filtered.update()

nlsave("aSi-interstitials.hdf5", filtered)

C-aSi-interstitials.py

As expected 3 interstitials are found:

+------------------------------------------------------------------------------+
| InterstitialList defect summary                                              |
+------------------------------------------------------------------------------+
| 3 inequivalent defect positions:                                             |
|   * Wyckoff 40 (46 defects) is C_40, site 040:                               |
|     - C interstitial of type 40 at (0.3521, 0.6386, 0.1991)                  |
|   * Wyckoff 106 (71 defects) is C_106, site 106:                             |
|     - C interstitial of type 106 at (0.6514, 0.5564, 0.5461)                 |
|   * Wyckoff 185 (72 defects) is C_185, site 185:                             |
|     - C interstitial of type 185 at (0.9471, 0.7745, 0.9669)                 |
+------------------------------------------------------------------------------+

and all of them are converged after updating:

+------------------------------------------------------------------------------+
| InterstitialList. Simulations: 3/3. Tasks: 15/15. Cores per item 1.          |
+------------------------------------------------------------------------------+
| Name                                               Tasks      Status Cores   |
| interstitial/C_40/040                                5/5   Converged     1   |
| interstitial/C_106/106                               5/5   Converged     1   |
| interstitial/C_185/185                               5/5   Converged     1   |
+------------------------------------------------------------------------------+

Transition path computation

The transition paths between all the inequivalent impurity configurations (3 in this example) is computed next by a simple call to TransitionPathList. Some extra options are used to improve the convergence ratio of simulations

from SMW import *

ints = nlread("aSi-interstitials.hdf5", InterstitialList)[-1]
nlprint(ints)

nebs = TransitionPathList(ints, optimize_transition_path_parameters=
    OptimizeTransitionPathParameters(number_of_images=3, number_of_fixed_atoms=0))
nebs.update()
nlprint(nebs)

nlsave("C-aSi-nebs.hdf5", nebs)

C-aSi-NEBs.py

Sentaurus Materials Workbench finds 6 extra defects (plus the three already calculated) that are required as states for the NEBs, and updates them all:

+------------------------------------------------------------------------------+
| TransitionPathDefectsList. Simulations: 3/9. Tasks: 15/45. Cores per item 1. |
+------------------------------------------------------------------------------+
| Name                                               Tasks      Status Cores   |
| interstitial/C_4/004                                 0/5     Running     1   |
| interstitial/C_24/024                                0/5      Queued     1   |
| interstitial/C_38/038                                0/5      Queued     1   |
| interstitial/C_40/040                                5/5      Queued     1   |
| interstitial/C_80/080                                0/5      Queued     1   |
| interstitial/C_98/098                                0/5      Queued     1   |
| interstitial/C_99/099                                0/5      Queued     1   |
| interstitial/C_106/106                               5/5      Queued     1   |
| interstitial/C_185/185                               5/5      Queued     1   |
+------------------------------------------------------------------------------+

After that, 6 different paths are run. Only 4 pahts converge:

+------------------------------------------------------------------------------+
| TransitionPathList. Simulations: 6/6. Tasks: 12/12. Cores per item 1.        |
+------------------------------------------------------------------------------+
| Name                                               Tasks      Status Cores   |
| interstitial/C_40/040 -> interstitial/C_80/080       2/2 Unconverged     1   |
| interstitial/C_40/040 -> interstitial/C_38/038       2/2   Converged     1   |
| interstitial/C_40/040 -> interstitial/C_24/024       2/2 Unconverged     1   |
| interstitial/C_106/106 -> interstitial/C_99/099      2/2   Converged     1   |
| interstitial/C_106/106 -> interstitial/C_98/098      2/2   Converged     1   |
| interstitial/C_185/185 -> interstitial/C_4/004       2/2   Converged     1   |
+------------------------------------------------------------------------------+

With the following migration energies printed by nlprint:

+------------------------------------------------------------------------------+
 | TransitionPathList paths summary: 6 paths                                    |
 +------------------------------------------------------------------------------+
 |   Path: #0, Wyckoffs 40 -> 80.                                               |
 |     Starts: interstitial/C_40/040                                            |
 |     Ends:   interstitial/C_80/080                                            |
 |          Converged?: no . Em (F) =  0.27 eV, Em (R) =  0.15 eV.              |
 |   Path: #1, Wyckoffs 40 -> 38.                                               |
 |     Starts: interstitial/C_40/040                                            |
 |     Ends:   interstitial/C_38/038                                            |
 |          Converged?: yes. Em (F) =  0.27 eV, Em (R) =  0.00 eV.              |
 |   Path: #2, Wyckoffs 40 -> 24.                                               |
 |     Starts: interstitial/C_40/040                                            |
 |     Ends:   interstitial/C_24/024                                            |
 |          Converged?: no . Em (F) =  0.67 eV, Em (R) =  2.96 eV.              |
 |   Path: #3, Wyckoffs 106 -> 99.                                              |
 |     Starts: interstitial/C_106/106                                           |
 |     Ends:   interstitial/C_99/099                                            |
 |          Converged?: yes. Em (F) =  2.89 eV, Em (R) =  0.00 eV.              |
 |   Path: #4, Wyckoffs 106 -> 98.                                              |
 |     Starts: interstitial/C_106/106                                           |
 |     Ends:   interstitial/C_98/098                                            |
 |          Converged?: yes. Em (F) =  0.83 eV, Em (R) =  2.44 eV.              |
 |   Path: #5, Wyckoffs 185 -> 4.                                               |
 |     Starts: interstitial/C_185/185                                           |
 |     Ends:   interstitial/C_4/004                                             |
 |          Converged?: yes. Em (F) =  0.34 eV, Em (R) =  0.00 eV.              |
 +------------------------------------------------------------------------------+

It can be seen that there is a very big variation in migration energies between the different paths, stressing the importance of performing Monte Carlo sampling to properly collapse all the values into a diffusivity measurement.

Diffusivity computation

The final diffusion in the amorphous sample is provided by the physical collapsing of all the possible paths (all of them with different starting and end formation energies, and different transition barriers), by means of Monte Carlo sampling. This is done with the monteCarloDiffusivityCalculation function.

from SMW import *

defects = nlread("aSi-interstitials.hdf5", InterstitialList)[-1]
nebs = nlread("C-aSi-nebs.hdf5", TransitionPathList)[-1]

diffusivity = monteCarloDiffusivityCalculation(
    [defects],
    [nebs],
    number_of_samples=200,
    number_of_jumps=1000,
    amorphous_material=True,
    discard_faulty=False,
)

print("Diffusivity is", diffusivity)

C-aSi-AKMC.py

That produces the following result:

Diffusivity is (PhysicalQuantity(1.67255734304,cm**2/s), PhysicalQuantity(2.32500019595,eV))

Giving us the diffusivity in silicon for this application: 2.32 eV.

Note

This small example is very sensitive to variation, so that you can expect different results than the ones shown here.