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:
DefectClusterList for defects containing an agglomeration of any or several of VacancyList, InterstitialList, SplitInterstitialList and SubstitutionalList.
InterstitialList for interstial defects occupying empty spaces between lattice sites.
DefectPairList for substitutional defects paired with vacancy or self interstitial defects.
SplitInterstitialList for split interstitial defects, i.e., interstitial defects that make extra room by displacing a lattice atom out their site
SubstitutionalList for substitutional defects that replace an element in a lattice site.
VacancyList for vacancy defects.
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.
Defect name, in particular cluster, interstitial, pair, split, substitutional and vacancy
Defect description, usually including the defect Wyckoff type. See each particular defect for the description.
Defect site, different for each defect, even if they are equivalent.
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')
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
}
Transition Path calculation¶
For transition path calculations, needed to compute migration barriers, see the TransitionPathList page.
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.
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
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.
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.
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
Define reference (MaterialSpecifications) imported from database or defined by user.
Calculate either heat of formation or specific element-rich condition in AtomicChemicalPotentialList
Update calculated chemical potential in reference (MaterialSpecifications) system.
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],[])
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())
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:
Define reference (MaterialSpecifications) imported from database or defined by user.
Calculate defects
Calculate GTP chemical potential using grandThermodynamicPotential
Update calculated chemical potential in reference (MaterialSpecifications) system.
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)
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.
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
Creation of a suitable amorphous material using amorphize.
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.
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()
orSentaurusMaterialsWorkbench.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.Computing migration paths between the selected defects using TransitionPathList.
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')
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)
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)
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)
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)
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.