API reference

Spinney: A Python package for first-principles Point Defect calculations.

Spinney is a collection of Python modules aimed for the analysis and postprocessing of first-principles calculations of point defects in solids.

Spinney can assists with the major tasks necessary for the characterization of point defects in solids. The classes and functions that it implements can be divided into the following groups, which are related to the several steps necessary for processing the ab-initio calculations:

• General high-level interface for point-defect calculations

  spinney.structures.pointdefect

  Contains the PointDefect class which offers a convenient interface to calculate the properties of a point defect, such as its formation energy including finite-size corrections.

• Determination of the possible values of equilibrium chemical potentials

  spinney.thermodynamics.chempots

  Contains the Range class which is able to determine the possible equilibrium values of the chemical potentials given a set of competing phases.

  It also contains classes describing the chemical potentials of common gas phases, such as \(\mathrm{O}_2\), \(\mathrm{H}_2\), \(\mathrm{N}_2\), \(\mathrm{Cl}_2\) and \(\mathrm{F}_2\), as a function of temperature and pressure employing experimental data and empirical formulas.

• Correction schemes for electrostatic finite-size effects in supercells

  • spinney.defects.kumagai

    Implements the correction scheme proposed by Kumagai and Oba in Phys. Rev. B 89, 195205 (2014)

  • spinney.defects.fnv

    Implements the correction scheme proposed by Freysoldt, Neugebauer and Van de Walle in Phys. Rev. Lett. 102, 016402 (2009)

• Calculation of equilibrium defect properties

  • spinney.defects.diagrams

    Contains the spinney.defects.diagrams.Diagram class which allows to plot the defect formation energies as a function of the Fermi level and calculate charge transition levels.

  • spinney.defects.concentration

    Contains the spinney.defects.concentration.EquilibriumConcentration class that allows to calculate equilibrium properties, such as defects and carriers concentrations and the Fermi level position.

• General-purpose tools

  • spinney.tools.formulas

    Contains some helper functions useful for dealing with chemical formulas.

  • spinney.tools.reactions

    Contains helper functions useful for calculating reaction energies.

• Support for first-principles codes

  The Spinney package currently offers interfaces for these computer codes:

  • VASP: spinney.io.vasp

  • WIEN2k: spinney.io.wien2k

General high-level interface for point-defect calculations

• spinney.structures.pointdefect

• spinney.structures.defectivesystem

The pointdefect module

Module implementing a PointDefect class to store and process information concerning a system with point defects.

class spinney.structures.pointdefect.DummyAseCalculator(atoms)

Dummy calculator to be used when the calculations are done with a code not supported by ase. Use this dummy calculator only to pass custom Atoms objects to the PointDefect class. Any other use of ought to be avoided.

Parameters

atoms (ase.Atoms) – the Atoms instance to be used with a PointDefect instance

class spinney.structures.pointdefect.PointDefect(ase_atoms)

PointDefect class.

It represents a 3D-periodic system containing one or more point defects.

Parameters

ase_atoms (instance of ase.Atoms) –

the Atoms object created from a completed point defect calculation.

Note:

it is assumed that all the Atoms object have units of Angstrom for lengths and eV for energies. This is the default in ase.

my_name

an alphanumeric label for the instance

Type

string, optional

defect_position

the scaled positions of the defect in the supercell

Type

1D numpy array

defect_charge

the charge state of the defect

Type

float

pristine_system

the ase.Atoms describing the pristine system

Type

ase.Atoms instance

parent_compounds

the keys are: pristine, chemical_formula1, … the value associated to pristine is the ase.Atoms instance representing the pristine system. chemical_formula1 etc. are the chemical formulas of the other compounds involved in the creation of the point defect. These are added using set_parent_elements(). Their values are the ase.Atoms instances representing these compounds.

Excluding the pristine system, these data are not necessary, but can be used to test the chemical potential values to be used in the calculation of the defect formation energy.

For example, for an Oxygen vacancy in MgO, the keys of parent_compounds would be: ‘pristine’, ‘O2’, ‘Mg2’, which means the ase Atoms object for ‘Mg2’ contains 2 Mg atoms; i.e. it is the primitive HCP cell.

Type

dict

parent_elements

each element is an instance of an ase.Atoms object representing the compound describing the standard state of the elements involved in the defective system.

For example, for an Oxygen vacancy in MgO, these would be the ase.Atoms instances for ‘O2’ and ‘Mg2’

Type

dict, optional

Examples

To initialize a PointDefect instance, it is only necessary to have an initialized ase.Atoms object with attached a calculator that supports the Atoms.get_total_energy() method.

Suppose that the output, obtained by a first-principle code, describing a defective system is saved in the file output.fmt. And that format fmt can be read by ase.io.read(). Then the following snippet can be used to initialize a PointDefect object representing the defective system.

>>> import ase.io
>>> defect = ase.io.read('output.fmt')
>>> pdf = PointDefect(defect)

add_correction_scheme_data(**kwargs)

Add the extra data needed in order to calculate the defect formation energy with the choosen scheme

Parameters

kwargs (dict) –

• For the correction schemes ‘ko’ and ‘fnv’:

  keys: potential_pristine, potential_defective

  the values are numpy arrays.

  • for ‘ko’:

    the arrays have to contain the value of the electrostatic potential at the ionic sites. The order has to be the same of the one used in the ase.Atoms objects employed in the initialization of the PointDefect instance and employed in set_pristine_system()

  • for ‘fnv’:

    the arrays have to contain the electrostatic potential on a 3D grid. The file has to match the supercell used to initialize the PointDefect instance.

• For ‘fnv’:

  axis the unit cell axis along which the electrostatic potential will be averaged.

  defect_density (optional) a 3D array with the charge density that can be used to model the defect-induced one. This will be used to fit the model charge to the defect-induced charge density.

  x_comb (optional) a float between 0 and 1. Weight of the exponential function with respect to the Gaussian function in modeling the defect-induced charge density. Default x_comb = 0: the charge density is a pure Gaussian.

  gamma (optional) a float. The parameter of the exponential function. Default value is 1.

  beta (optional) a float. The parameter of the Gaussian function. Default value is 1.

  shift_tol (optional) a float representing the tolerance to be used in order to locate the defect position along axis. Default value: \(1e-5 \times \mathrm{length cell parameter of axis}\)

  e_tol (optional) a float, break condition for the iterative calculation of the correction energy. Value in Hartree. Default: 1e-6 Ha.

• For ko:

  distance_tol (optional) rounding tolerance for comparing distances, in units of Angstrom. Default value is 5e-2 Angstroms.

  e_tol (optional) a float, break condition for the iterative calculation of the correction energy. Value in eV. Default: 1e-6 eV.

calculate_finite_size_correction(verbose=False)

Calculate the energy correction for finite-size effects employing the method chosen with set_finite_size_correction_scheme.

Parameters

verbose (bool) – If True, several details are returned as a dictionary. If False, only the correction energy is returned. This has to be added to the energy of the defective system.

get_defect_formation_energy(include_corr=False)

Returns the formation energy of the defective system.

Parameters

include_corr (bool) – If True, the formation energy is already corrected for finite-size effects.

Returns

energy – the defect formation energy

Return type

float

set_Eg(value)

Parameters

value (float) – The system band gap value.

set_chemical_potential_ranges(ranges)

For each element involved in the creation of the defective system, specify the minimum and maximum value that the chemical potential can have.

If set_parent_elements() was used, it will be checked that the input elements are the same.

Parameters

ranges (dictionary of 2-ples) –

For each element involved in the creation of the point defect, an element of ranges[element] is the 2-ple:

(minimum value atomic chemical potential, maximum value atomic chemical potential)

With ‘atomic chemical potential’ it is intended that the values are referred to a single atom, and not to the formula unit of the corresponding element used in parent_elements

set_chemical_potential_values(chem_pots, force=False)

Set the chemical potential values, for each species involved in forming the defective system, to be used in the calculation of the defect formation energy.

If these are not set, an exception will be raised. If chemical_potential_ranges is not None, it will be checked if the given chemical potentials are within these ranges; if not, an exception will be raised. In any case, the given chemical potential values will be tested againts the parent elements chemical potentials, if these are given. Use force=True to bypass these checks.

Parameters

• chem_pots (dict of floats) – each element is the chemical potential of one element (value give per atom) The keys must be the same used in self.parent_elements, if the parent elements were set.

• force (bool) – If True, the inserted values will be used for the calculations, even tough they are not physical valid values

set_defect_charge(charge)

Parameters

charge (float) – The formal charge assigned to the point defect. It is assumed that only one localized charge is present in the supercell.

set_defect_position(position)

Set the fractional coordinates of the point defects in the system.

Parameters

position (1D array) – the scaled position of the point defect with respect to self.get_cell()

set_dielectric_tensor(value)

Parameters

value (2D array or float) – The value of the system dielectric tensor (or constant).

set_fermi_level_value_from_vbm(value)

Set the value of the Fermi level with respect to the system valence band maximum.

This value will be used to calculate the defect formation energy of charged defects.

Parameters

value (float) –

set_finite_size_correction_scheme(scheme)

Set the correction scheme for finite-size effects in point defect calculations.

Parameters

scheme (string) –

The correction scheme to use. Possible values:

• ’ko’ : Kumagai and Oba, PRB 89, 195205 (2014)

• ’fnv’: Freysoldt, Neugebauer, and Van de Walle, Phys. Rev. Lett. 102, 016402 (2009)

set_parent_elements(elements)

Set the ase.Atoms objects representing the compounds, in their reference state, for the elements involved in the formation of the point defects.

Parameters

elements (list or tuple) – each element of the sequence is an Atoms object. For example, if the defective system consists of the C vacancy- N_C complex in diamond, then: elements = (ase.Atoms for C (e.g. C diamond), ase.Atoms for N2)

set_pristine_system(ase_pristine)

Parameters

ase_pristine (ase.Atoms) – The ase.Atoms object representing the pristine system.

set_vbm(value)

Parameters

value (float) – The system valence band maximum, which determines the minimum value of the electron chemical potential.

The defectivesystem module

class spinney.structures.defectivesystem.DefectiveSystem(data_path, calculator)

Container class describing a system with point defects.

Parameters

• data_path (string) –

  path to the folder containing the results of the point defect calculations. It is expected a directory tree like this:

  "data_path"
   ├── data_defects
   │   ├── "defect_name"
   │   │   ├── "charge_state"
   │   │   │   └── "files"
   │   │   ├── "charge_state"
   │   │   │   └── "files"
   │   │   ├── "charge_state"
   │   │   │   └── "files"
   │   │   └── "charge_state"
   │   │       └── "files"
   │   ├── "defect_name"
   ...
   └── pristine
       └── "files"
  

  • defect_name is a string describing the point defect.

  • charge_state must be the charge state of the considered defect.

  • files are the data needed for calculating defect formation energies. These depends on the calculator in use:

  • VASP: at least there should be the OUTCAR file.

  • WIEN2k: at least there should be the case.struct and case.scf files.

  For all calculators a file named position.txt containing the fractional coordinates of the defective site must be present.

• calculator (string) – the code used to calculate the data

vbm

the valence band maximum of the host material

Type

float

dielectric_tensor

the dielectric tensor of the host material

Type

float or 2D array

chemical_potentials

a dictionary whose keys are the parent elements forming the defective system and whose values are the chemical potentials to be used in the calculation of the defect formation energy

Type

dict

correction_scheme

specifies the correction scheme for finite-size-effects to be used

Type

str

data

collects the calculated formation energies for each processed point defect

Type

pandas DataFrame object

point_defects

a list of PointDefect objects corresponding to the processed point defects

Type

list

gap_range

a tuple containing the valence band maximum and conduction band minimum of the pristine system

Type

tuple

extended_gap_range

a tuple containing the valence band maximum and conduction band minimum for an extetnded band gap of the pristine system. Used to initialize a Spinney Diagram object

Type

tuple

diagram

an object representing the defect formation energies as a function of the electron chemical potential

Type

Spinney Diagram object

calculate_energies(verbose=True)

Calculate defect formation energies :param verbose: if True, information about the process is printed :type verbose: bool

property defects_degeneracy_numbers

A dictionary o Represents the degeneracy for each type of defect in each of its charge states. The order has to match that of charge_states[defect_type].

write_formation_energies(out_file)

Write the defect formation energies to a file in a format used by Spynney.

Parameters

out_file (str) – name of the file used to write the energies

Determination of the possible values of equilibrium chemical potentials

• spinney.thermodynamics.chempots

The chempots module

Module with general tools and classes for handling chemical potential-related quantities.

The Range class allows to calculate valid chemical potential values given competing phases.

class spinney.thermodynamics.chempots.ChlorineChemPot(energy_units='eV', pressure_units='Pa')

Class for modelling the chemical potential of the \(\mathrm{Cl}_2\) gas molecule

Parameters

• energy_units (string) – the units to be used for the energy

• pressure_units (string) – the units to be used for the pressure

class spinney.thermodynamics.chempots.FluorineChemPot(energy_units='eV', pressure_units='Pa')

Class for modelling the chemical potential of the \(\mathrm{F}_2\) gas molecule

Parameters

• energy_units (string) – the units to be used for the energy

• pressure_units (string) – the units to be used for the pressure

class spinney.thermodynamics.chempots.HydrogenChemPot(energy_units='eV', pressure_units='Pa')

Class for modelling the chemical potential of the \(\mathrm{H}_2\) gas molecule

Parameters

• energy_units (string) – the units to be used for the energy

• pressure_units (string) – the units to be used for the pressure

class spinney.thermodynamics.chempots.IdealGasChemPot(energy_units='eV', pressure_units='Pa')

Class for modelling the chemical potential of an ideal gas molecule

Parameters

• energy_units (string) – the units to be used for the energy

• pressure_units (string) – the units to be used for the pressure

G_diff_Shomate_Eq(T)

Calculate the standard Gibbs free energy at temperature T with respect to the standard one at 0 K using Shomate equation.

Parameters

T (float) – the temperature

Returns

gibbs_energy – the standard Gibbs free energy at T in the units of self.energy_units

Return type

float

get_ideal_gas_chemical_potential_Shomate(mu_standard, partial_pressure, T)

Given the standard chemical potential at 0K, e.g. calculated with DFT, returns the value at given temperature obtained using the Shomate equation and the ideal gas formulas.

\[\mu(T, p) = h(0, p^\circ) + [g(T, p^\circ) - g(0, p^\circ)] + k_BT\ln\left(\frac{p}{p^\circ}\right)\]

Parameters

• mu_standard (float) – standard-state chemical potential of the molecule at 0K

• partial_pressure (array or float) – the partial pressure

• T (array or float) – the gas temperature

Returns

chem_pot – The chemical potentials as a 2D numpy array of shape (len(T), len(partial_pressure)) if both T and partial_pressure are arrays; otherwise a 1D numpy array of length len(partial_pressure)

if partial_pressure is an array and T a float; otherwise a 2D array of shape (len(T), 1) if partial_pressure is a float and T an array; finally if both T and partial_pressure are float, the result is a float

Return type

numpy array/float

class spinney.thermodynamics.chempots.NitrogenChemPot(energy_units='eV', pressure_units='Pa')

Class for modelling the chemical potential of the \(\mathrm{N}_2\) gas molecule

Parameters

• energy_units (string) – the units to be used for the energy

• pressure_units (string) – the units to be used for the pressure

class spinney.thermodynamics.chempots.OxygenChemPot(energy_units='eV', pressure_units='Pa')

Class for modelling the chemical potential of the \(\mathrm{O}_2\) gas molecule

Parameters

• energy_units (string) – the units to be used for the energy

• pressure_units (string) – the units to be used for the pressure

class spinney.thermodynamics.chempots.Range(coeff_equalities, const_equalities, coeff_inequalities, const_inequalities, bounds)

Class for finding the allowed ranges for the chemical potentials of given elements given the competing phases.

Parameters

• coeff_equalities (2D tuples) – the coefficients of the linear equalities. For each equality there should be an array with the corresponding coefficients. Each element in the tuple must contain the same number of elements. If a variable appears at least once in the constraint equations, it must be set to 0 in all the other equations.

• const_equalities (1D tuple) –

  the constant values of the linear equalities. const_equalities[i] has to be the constant value of the linear equation with coefficients coeff_equalities[i].

  Example:

  We have the linear constraints given by these equations:

  \begin{align*} ax + by + cz &= d \\ a'x + b'y + c'z &= d' \\ a''x + b''y = d'' \end{align*}

  then:

  coeff_equalities = ((a,b,c), (a’, b’, c’), (a’’, b’’, 0))

  const_equalities = (d, d’, d’’)

  Warning

  If a variable is present anywhere in either the equality or inequality constraints, then it must be explicitly given even if its value is zero (the coefficient of z in the last equation of the previous example).

  coeff_inequalities and const_inequalities are analogously defined, but for the inequality conditions.

• bounds (tuple of tuples) – the bounds of the independent variables. If not present, they will be set to (None, None) for each independent variable, which will be treated as -inf and +inf

number_of_variables

number of chemical potentials

Type

int

variables_extrema

shape = (number elements, 2) for each element, returns the minimum and maximum possible values for its chemical potential

Type

2D numpy array

variables_extrema_2d

shape = (2, 2) for each element, returns the minimum and maximum possible values for its chemical potential. Calculated after an intersection with a plane has been performed.

Type

2D numpy array

mu_labels

the symbols labeling the elements for which we caculate the chemical potentials

Type

tuple

Notes

ALL INEQUALITIES SHOULD BE IN THE FORM:

coeff_inequalities * x <= const_inequalities

Examples

Suppose that we want to find the chemical potential extrema for Ti and O in rutile TiO2.

We have the following constraints:

\begin{align*} \mu_\mathrm{Ti} + 2\mu_\mathrm{O} &= \mu_\mathrm{TiO_2}\\ \mu_\mathrm{Ti} + \mu_\mathrm{O} &\leq \mu_\mathrm{TiO}\\ 2\mu_\mathrm{Ti} + 3\mu_\mathrm{O} &\leq \mu_\mathrm{Ti_2O_3}\\ \mu_\mathrm{Ti} &\leq \mu_\mathrm{Ti}(HCP)\\ \mu_\mathrm{O} &\leq \frac{1}{2} \mu_\mathrm{O_2} \end{align*}

So, we will have to type:

>>> coeff_equalities = ((1,2))
>>> const_equalities = (mu_TiO2, )
>>> coeff_inequalities = ((1, 1), (2, 3))
>>> const_inequalities = (mu_TiO, mu_Ti2O3)
>>> bounds = ((None, mu_Ti(HCP)), (None, 1/2 mu_O2))

check_value_variables(variables, eq_tol=1e-06)

Checks whether a set of variables is within the feasible region.

Parameters

• variables (array) – the length is given by no_variables: the number of elements in the system

• eq_tol (float) – tolerance on equality conditions

Returns

result – a booleand and a numpy array of shape self.no_variables times the number of constraints

• First element:

  True, if variables is within the feasible region, False, otherwise.

• Second element:

  Array of shape: no_variables times the number of constraints where the number of contraints is the number of equalities plus inequalities plus bounds conditions.

  This array is a 2D array whose rows represent the variables and whose columns represents the various constraints, the order is:

  equalities, inequalities and bound down, bound up

  in the same order given by the user at initialization. The values of such matrices are bool.

  Example:

  Suppose we have two variables, one equality and two inequality and 2 bounds constraints. If both variables satisfy the equality, first inequality and bounds. But variable 2 does not satisfy the second inequality; the funcion will return a np.ndarray result of shape (2,5):

  result = [[True, True, True, True, True],
            [True, True, False, True, True]]
  

Return type

tuple

plot_feasible_region_on_plane(plane_axes, plane_values=None, save_plot=False, tol=1e-06, **kwargs)

Plot the feasibile region intersection with the plane specified by plane_axes.

Parameters

• plane_axes (1D array of length 2) – the array specifies the indices of the variables which will form the axes of the intersection plane.

• plane_values (1D array of length self.number_of_variables - 2) – Constant value specifying the plane.

• save_plot (bool) – if True, the plot will be saved

• tol (float) – a numerical tolerance for determining the points giving the feasible region

• kwargs (dict, optional) – optional key:value pairs for plotting the diagram

Returns

Return type

the figure instance of the plot.

Example

For example, if we are considering a ternary system A, B, C, and we want to intersect the feasible region with the plane B = b, using the axes given by variables A and B, then plane_axes = [0, 2], plane_values = [b]. This specifies the plane \(B = b\).

By default plane_values is filled with zeros.

set_chemical_potential_labels(labels)

Set the name of the independent variables.

Parameters

labels (list) – the names of the independent variables, one should use the same order used in self.coeff_*

set_color_map(cmap)

Set the color map to be used in the plot

Parameters

cmap (matplotlib.cm class) – one of the color maps offered by matplotlib. By default Set1 is used, which is fine when less than 10 compounds are considered

set_compound_dict(compounds_dict)

Set a dictionary which contains the compound names relative to the equality, inequality and bound conditions.

Valid names are: formula + _ + description.

description is any string that described the compound

Note

Note that the validity of the name format is not checked internally and should be done by the user.

compounds_dict : dict

{'equality' : [compoundA1, ...],
 'inequality' : [compound1, ...],
 'bound' : [compoundA, ...]}

the order of the elements in the dictionary values has to be the same of the order used in the relative const_equalities, const_inequalities and bounds lists

spinney.thermodynamics.chempots.find_partial_pressure_given_mu(mu, mu_standard, T, pa=1e-12, pb=1000000000.0, energy_units='eV', pressure_units='Pa')

Finds the partial pressure required to change the chemical potential of the ideal gas from mu_standard to mu at temperature T.

Parameters

• mu (float) – the desired value of the chemical potential

• mu_standard (float) – the standard state chemical potential of the ideal gas

• T (float) – the temperature of interest

• pa (float) – the lower bound for the partial pressure

• pb (float) – the upperbound for the partial pressure

• energy_units (string) – the unit of energy used for the value of mu and mu_standard

• pressure_units (string) – the unit of pressure

Returns

p – the partial pressure giving mu

Return type

float

spinney.thermodynamics.chempots.get_chem_pots_conditions(file, order, equality_compounds, delimiter=None)

Reads file and returns the data needed to initialize a Range class as well the dictionary describing the various compounds.

Parameters

• file (str) –

  the location of the file with the required information The format of this file must be:

  Label_compound(delimiter)energy

  Label_compound is a string, whose different entries must be separated by an underscore.

  delimiter is a string specifying the delimiter used to separate Label_compound from energy

  energy is a float representing the (formation) energy of the compound

• order (array) – ordered list of the chemical symbols to be used in specifying the equality and inequality conditions

• equality_compounds (array of strings) – the compounds in file corresponding to the equality conditions

Returns

result – The first 5 elements are needed to initialize a Range instance.

The last element is the argument for the method Range.set_compound_dict()

Return type

tuple of 6 elements

spinney.thermodynamics.chempots.get_conditions_from_file(file, order, delimiter=None)

Reads and parses a file and returns a triplet with the elements, coefficients, and the constants. These values can be used for the inequality or equality conditions in Range.

Lines starting with ‘#’ will be skipped.

Parameters

• file (str) – the file that has to be read. The format has to be: Compound_formula_unit(delimiter)Energy_of_the_compound_per_fu

• Note – All inequalities are indented as “<=”: pay attention to prepare the input file in this way

• order (array) – the name of the elements in the order that has to be used for the coefficients. Eg, for the Mn-O system, order can be [‘Mn’, ‘O’] or [‘O’, ‘Mn’]

• delimiter (str) – the symbol used to separate compound names from energy values.

Returns

compound_dict – the keys are the compound identifiers as read from the file. The elements of the tuple are:

• the formula unit of the compound

• the stoichiometric coefficients for the elements in the compound

• the (formation) energy per atom of the compound

Return type

dictionary of 3-ple:

spinney.thermodynamics.chempots.ideal_gas_chemical_potential(mu_standard, partial_pressure, T, energy_units='eV', pressure_units='Pa')

Temperature and pressure-dependent chemical potential of an ideal gas

\[\mu = \mu^\circ + k_B T \ln \left(\frac{p}{p^\circ}\right)\]

Parameters

• mu_standard (float) – the gas molecule chemical potential at standard pressure

• partial_pressure (array or float) – the gas partial pressure with respect to the standard pressure

• T (array or float) – the temperature range of interest

• energy_units (string) – the units in which energy is expressed

• pressure_unit (string) – the units in which pressure is expressed

Returns

chem_pots – The chemical potentials as a 2D numpy array of shape (len(T), len(partial_pressure)) if both T and partial_pressure are arrays; otherwise a 1D numpy array of length len(partial_pressure) if partial_pressure is an array and T a float; otherwise a 2D array of shape (len(T), 1) if partial_pressure is a float and T an array; finally if both T and partial_pressure are float, the result is a float

Return type

array/float

spinney.thermodynamics.chempots.parse_mu_from_string(string, delimiter=None)

Reads and parses a string and returns a dictionary with the elements, coefficients, and the constants of the equality/inequality conditions.

Parameters

• string (str) –

  the string that has to be parsed. The format has to be:

  Formula_unit(delimiter)Energy_of_the_compound_per_fu # comments

  Example:

  SrTiO3 -12.3

  will return:

  result = {'Sr':1/5, 'Ti':1/5, 'O':3/5, 'energy':-12.3}

• delimiter (str) – the symbol used to separate compound names from energy values

Returns

result – as in the example above

Return type

dict

Correction schemes for electrostatic finite-size effects in supercells

• spinney.defects.kumagai

• spinney.defects.fnv

The kumagai module

Implementation of the correction scheme for charged point defects of Kumagai and Oba:

25. Kumagai and F. Oba, PRB 89, 195205 (2014)

class spinney.defects.kumagai.KumagaiCorr(cell, positions_defective, positions_pristine, defect_position, defect_formal_charge, dielectric_constant, dft_core_potential_def, dft_core_potential_prist, direct_cutoff=10, reciprocal_cutoff=1, alpha=None, length_units='Angstrom', energy_units='eV', tol_en=1e-06, min_steps=2, tol_dist=0.01)

Implementation of the Kumagai correction scheme of: PRB 89, 195205 (2014)

Parameters

• cell (2D numpy array) – cell parameters defective supercell

• positions_defective/positions_pristine (2D numpy array) – fractional coordinates of the atoms in the DEFECTIVE and PRISTINE supercells of same size and shape.

• defect_position (1D numpy array) – fractional coordinates of the point defect

• defect_formal_charge (float) – charge of the point defect

• dielectric_constant (2D numpy array) – the dielectric tensor

• dft_core_potential_def/dft_core_potential_prist (1D numpy array) – potential at the ionic sites calculated by first-principles for DEFECTIVE and PRISTINE systems, respectively. Same atomic ordering as in the corresponding positions_* arrays

• direct_cutoff/reciprocal_cutoff/alpha – see the spinney.defects.madelung.Ewald class

• length_units/energy_units (strings) – the unit of length and energy used for cell and dft_core_potential_*

• tol_dist (float) – rounding tolerance for comparing distances, in units of length_units

atomic_mapping

element [0] has the indexes of the pristine system atoms with corresponding atoms in the defective system

element [1] has the indexes of the defective system mapped by map[0]. Where map is returned by map_atoms_pristine_defect()

Type

2-ple of 1D arrays

atoms_in_sampling_region

the first element contains the indices of the atoms of the pristine system in the sampling region, the second element contains the corresponding atoms of the defective system

Type

2-ple

difference_potential_vs_distance

element [0] is an array containing the distances from the point defects where the electrostatic core potential has been sampled

element [1] is the difference of the defective and pristine core potentials at the point whose distance from the defect is in element [0]

Type

2-ple of 1D arrays

difference_potential_vs_distance_sampling_region

as above, but only the sampling region is considered

Type

2-ple of 1D arrays

alignment_potential_vs_distance_sampling_region

element [0] is an array containing the distances from the defect of the sites within the sampling region.

element [1] is the corresponding alignment potential: site potential defect - site potential pristine - ewald potential

Type

2-ple of 1D arrays

ewald_potential_vs_distance_sampling_region

the ewald potential calculated at the sites in the sampling region. Structure analogous to alignment_potential_vs_distance_sampling_region

Type

2-ple of 1D arrays

grouped_atom_by_distance

the first element of the tuple reports the distances from the defect The second element is another 2-ple; each element of this tuple contains a list of lists. Each inner list contains the atomic index of the atoms that are at a given distance (within tol_dist) from the defective site for the pristine and defective supercell.

Type

2-ple

property difference_potential_vs_distance

Pot_def - Pot_pristine

spinney.defects.kumagai.kumagai_sampling_region(cell, atom_coordinates, defect_position)

Given the scaled coordinates of atoms in the supercell, returns the indices of the atoms belonging to the sampling region as defined in: Y. Kumagai and F. Oba, PRB 89, 195205 (2014)

Parameters

• cell (2D numpy array) – each row represents the Cartesian coordinates of the primitive vectors

• atom_coordinates (2D numpy array) – each row represents the fractional coordinates of the atoms in the supercell

• defect_position (1D numpa array) – the scaled coordinates of the defect atom in the supercell. It will be used to center the supercell

Returns

in_region, sphere_radius – in_region is a list with the indexes of the atoms in the sampling region.

sphere_radius is a float and it is the radius of the sphere inscribed in the Wigner-Seitz cell.

Return type

tuple

spinney.defects.kumagai.map_atoms_pristine_defect(cell, scaled_pos_prist, scaled_pos_defect, dist_tol=0.01)

Return a mapping between the indexes of the atoms in the pristine supercell and those of the atoms in the defective supercell. Both supercells have to be the same.

Parameters

• cell (2D numpy array) – the supercell

• scaled_pos_prist (2D numpy array) – fractional atomic positions of the pristine system

• scaled_pos_defect (2D numpy array) –

  fractional positions of the defective system

  Warning

  It is assumed that both pristine and defective system have the same cell: cell

• dist_tol (float/array of 3 elements) – tolerance value for which 2 distances are considered equal. If float, the value will be used to create an array with 3 elements. Each elements represent the tolerance along the corresponding cell parameter. Default value: 1% of the corresponding cell parameter.

Returns

map – map[0] has the indexes of the pristine system with correspective atoms in the defective system

map[1]`has the indexes of the defective system atoms mapped by :data:`map[0]

Return type

2-ple of 1D arrays

The fnv module

Standalone implementation of the correction scheme for charged defects proposed by Freysoldt, Neugebauer and Van de Walle.

Phys. Rev. Lett. 102, 016402 (2009)

Internally, the code use atomic units of Hartree and Bohr, as implemented in Freysoldt et al. papers.

The input values have to be inserted in units of eV and Angstrom and the final values will also be converted to these units.

class spinney.defects.fnv.FChargeDistribution(x=0, gamma=1, beta=1)

Gaussian + exponential tail radial charge distribution employed by Freysoldt et. al., Phys. Status Solidi B 248, 5 (2011)

\begin{align*} q(r) &= q \left(x N_1 e^{-r/\gamma} + (1-x) N_2 e^{-r^2/\beta^2} \right)\\ N_1 &= \frac{1}{8 \pi \gamma^3}\\ N_2 &= \frac{1}{(\sqrt \pi \beta)^3} \end{align*}

x and gamma may be found from looking at the defect wave function. The value of beta is not so important as long as the defect remains localized. A Gaussian function usually performs well enough.

Parameters

• x (float) –

• gamma (float) – units of Angstrom

• beta (float) – units of Angstrom

Notes

For localized states, x = 0; for delocalized ones, usually x is around 0.54-0.6 in semiconductors.

Direct()

The normalized charge distribution in real space

Fourier_transform()

The Fourier trasform of Direct()

Fourier_transform_gto0()

The second derivative of the charge distribution with respect to g. It approximates the distribution for \(g \rightarrow 0\):

q(g-->0) ~ q.Fourier_transform(0) + 1/2 q.Fourier_transform_gto0()

The term is used for getting the potential alignment \(V_0\) of equation 17.

classmethod get_model_D()

Model function to be used for fitting

class spinney.defects.fnv.FCorrection(supercell, q, charge_distribution, pristine_pot, defective_pot, dielectric_constant, def_position, axis_average, cutoff=500, tolerance=1e-06, shift_tol=1e-05)

Implementation of Freysoldt et al. correction.

Note

the code assumes that the units are Angstrom for lengths and eV for energies!

Parameters

• supercell (2D array) – each rows represents the cartesian coordinates of the vectors describing the supercell

• q (int) – charge state of the defect

• pristine_pot (3D numpy array) –

  this array must contain the total electrostatic potential calculated on 3D grid for the pristine system. The shape of the array is: (Na, Nb, Nc)

  where Nx is the number of grid points along the cell parameter x

• defective_pot (3D numpy array) – as pristine_pot but for the defective system

• charge_distribution (FChargeDistribution instance) – the model charge distribution used for the defect-induced charge density

• dielectric_constant (float) – dielectric constant of the bulk system

• def_position (array) – defect position in the defective supercell in fractional coordinates with respect to supercell

• axis_average (int) –

  axis which will be used to perform the plane-average of the electrostatic potential:

  • 0 for a

  • 1 for b

  • 2 for c

• cutoff (float) – cutoff for reciprocal space vectors, IN eV

• tolerance (float) – tolerance threshold used in calculating the energy terms

• shift_tol (float) – tolerance in findind the defect position on the 3D grid

generate_reciprocal_space_grid(cutoff)

Generates a regular mesh of reciprocal space points given a cutoff energy.

Parameters

cutoff (float) – maximum kinetic energy, in Hartree, of the k-point

Returns

grid – the reciprocal lattice vectors in Cartesian coordinates in units of Bohr

Return type

list

get_E_lat(all_contributions=False)

Gets the \(E^{lat}[q]\) part of the correction scheme. From equation (8) in Freysoldt et al. PSS B 248, 5 (2011).

Parameters

all_contributions (bool) – if True, the function returns the various contributions: E_lat, E_lat1, E_lat2; otherwise just the total E_lat

Returns

energy_values – The values IN eV.

Return type

float/tuple

get_correction_energy()

Returns the correction energy for finite-size effects to add to the calculated DFT energy.

Returns

E_corr – The correction term

Return type

float

get_potential_alignment()

Get the potential alignment term. I.e. the \(q \Delta V\) term in equation (13) of PSS B 248, 5 (2011)

class spinney.defects.fnv.FFitChargeDensity(cell, charge_model, density_proj_def, defect_position, axis_average, tol=1e-05)

Fits the model charge density to the DFT defective charge density averaged along an axis.

Parameters

• cell (2D numpy array) – each row reoresents the Cartesian coordinates of one cell parameter of the crystal unit cell

• charge_model (class) – model of the charge density

• density_proj_def (3D numpy array) –

  the defect-induced charge density on a 3D grid. In case of spin-polarized calculations, it is the net charge density. The shape of the array is: (Na, Nb, Nc)

  where Nx is the number of grid points along the cell parameter x

• defect_position (1D array) – the position of the defect with respect to the cell parameters (i.e. the fractional coordinates)

• axis_average (int) – the axis along which perform the plane-average

• tol (float) – numerical spatial tolerance

fit_model()

Fit the model charge density

class spinney.defects.fnv.FPlotterPot(coordinates, av_locpot_diff, av_lr_pot, av_sr_pot, pot_alignment, axis_name)

Helper class for plotting the potentials.

Parameters

• coordinates (array) – coordinates to plot the potentials along an axis, can be fractional or Cartesian

• av_locpot_diff (array) – difference DFT potential between defective and pristine system plane-averaged along an axis

• av_lr_pot (array) – the analytical long-range potential, plane-averaged over an axis

• av_sr_pot (array) – the short-range potential

• pot_alignment (float) – the final aligment of the potentials far from the defect

• axis_name (string) – the name of the axis along which the potentials were calculated

plot(title='')

Plot the potentials

Parameters

title (string, optional) – the plot title

spinney.defects.fnv.plane_average_potential(locpot, axis)

Makes the plane-average of the electrostatic potential along the chosen axis.

Parameters

• locpot (3D array) – the local electrostatic potential on a 3D mesh

• axis (int) – 0 for a 1 for b 2 for c where a,b,c are the cell parameters along which the 3D grid is defined

Returns

1D numpy array

Return type

the plane-averaged potential along the chosen axis

Calculation of equilibrium defect properties

• spinney.defects.diagrams

• spinney.defects.concentration

The diagrams module

Implements the creation of the diagrams reporting the formation energy of point defects as a function of the electron chemical potential.

class spinney.defects.diagrams.Diagram(defects_dictionary, gap_range, extended_gap_range=None, electron_mu=None)

A diagram represent the formation energies of various point defects, in various charge states, as a function of the electron chemical potential, whose value ranges from the valence band maximum to the conduction band minimum of the host material.

This class is basically a composition of PointDefectLines instances with some tools for plotting the final diagram.

Parameters

• defects_dictionary (dict of dicts) –

  for each point defect named “defect”, the value is a dictionary in the form: {charge_state1 : formation_energy_at_VBM, …}

  Example:

  For vacancy and N center in diamond:

  >>> defects_dictionary = {
          '$V$' : {-1 : value1, 0 : value2, 1 : value3},
          r'$N_C$' : {0 : value4}
          }
  

  The value of VBM must be consistent with gap_range[0]

• gap_range (array of 2 elements) – the band gap range. The first value must be consistent with the valence band maximum (VBM) used to calculate the defect formation energy in defect_dictionary

• extended_gap_range (array of 2 elements) – an extended range for the band gap. This value can be used for example when gap_range is the underestimated DFT band gap and we want to plot the defect formation energy lines on a wider gap. In this case, the difference between the extended region and the original region will be plotted in gray.

• electrom_mu (float) – pinned value of the electron chemical potential wrt the valence band maximum (gap_range[0])

labels

{defect_name : defect_label, …} where defect_name is one of the top keys if defects_dictionary

Type

dict

defects

{defect_name : Line instance representing that defect type, …}

Type

dict

plot(**kwargs)

Plot the diagram.

Parameters

kwargs (dict) – optional key:values pair for plotting the diagram

write_transition_levels(file_name)

Writes the charge transition levels of the system on a txt file.

Parameters

file_name (string) – the name of the file where the transition levels will be written

class spinney.defects.diagrams.Line(m, q, x0)

A simple class describing the equation of a line:

\[y = m(x-x_0) + q\]

where q is the intercept with the y = x0 axis

class spinney.defects.diagrams.PointDefectLines(defect_name, e_form_min_dict, gap_range, extended_gap_range=None)

This class holds the information to represent a line on the formation energy diagrams of point defects.

Parameters

• defect_name (string) – the name of the defect

• e_form_min_dict (dict) –

  {charge_state1 : formation_energy1_at_VBM, …}

  charge_state1 is an integer describing the defect charge state; formation_energy1_at_VBM is a floating-point number indicating the formation energy of that defect for an electron chemical potential equal to the valence band maximum (VBM) of the host material. This VBM value must be stored in gap_range

• gap_range (array of 2 elements) – (valence_band_maximum, conduction_band_minimum) representing the allowed range for the electron chemical potential. Note that the VBM to which e_form_dict refers to is gap_range[0]

• extended_gap_range (array of 2 elements) – an extended range for the band gap. This value can be used for example when gap_range is the underestimated DFT band gap and we want to evaluate the defect formation energy lines on a different gap.

lines

{charge_state1 : line1, …} for each defect charge state of the instance. line1 is a Line instance.

Type

dict

intersections

For each charge state of the defect, one has the dictionary containing the coordinates of the intersection point between the line corresponding to that charge state and the lines of all other charge states and band edges as well.

Type

dict of dicts

lines_limits

{defect_charge_state1 : [y_0, y_1], …} y_0 and y_1 are two numpy arrays with two elements each:

the values of the electron chemical potential and those of the defect formation energy corresponding to the two points between which the formation energy of defect_charge_state1 is the lowest one among all other defect charge states.

Type

defaultdict

transition_levels

Returns the transition levels among the possible charge states as a dictionary:

{charge_state1 : {charge_state2 : transition_level,
                  charge_state3 : transition_level, ...}, ...}

transition_level is the electrom chemical potential value at the transition level between charge states. For a complete list of intersection between all defect lines, use self.intersections

Type

dict

property intersections

Returns a dictionary of dictionaries. For each charge state of the defect, one has the dictionary containing the coordinates of the intersection point between the line corresponding to that charge state and all other charge states and band edges as well.

property lines_limits

The attribute _lines_limits is a defaultdict:

{defect_charge_state1 : [y_0, y_1], …} y_0 and y_1 are two numpy arrays with two elements each: the values of the electron chemical potential and those of the defect formation energy corresponding to the two points between which the formation energy of defect_charge_state1 is the lowest one among all other defect charge states.

This is a subset of self._intersections

plot_lines(ax='auto', **kwargs)

Given an existing axes, plots the defect formation energy lines

Parameters

• ax (matplotlib.axes.Axes instance) – the axes where the lines should be plotted

• kwargs (dict) – additional key:value pairs taken by matplotlib.pyplot.plot

property transition_levels

Returns the transition levels among the possible charge states as a dictionary:

{charge_state1{charge_state2transition_level,

charge_state3 : transition_level, …}, …}

transition_level is the electrom chemical potential value at the transition level between charge states. For a complete list of intersection between all defect lines, use self.intersections

spinney.defects.diagrams.extract_formation_energies_from_file(file_name)

Read a file with the format:

# arbitrary number of comment lines
defect_name    defect_charge_state   defect_formation_energy

Parameters

file_name (string) – the file with the data

Returns

diagram_dict – dictionary containing the information necessary to initialize an Diagram instance.

Return type

defaultdict(dict)

spinney.defects.diagrams.set_latex_globally()

Call this function to use latex synthax globally in matplotlip. Note that now every string used through matplotlib has to respect latex synthax.

The concentration module

Functions and classes for calculating defect and carrier concentrations and related quantities from their formation energies.

class spinney.defects.concentration.Carrier(dos, vbm, cbm, mu, T, dos_down=None, energy_units='eV')

Object describing a free carrier in a semiconductor.

Parameters

• dos (2D array) – the first column stores the 1-electron energy levels, the second column stores the corresponding DOS. The energy must be sorted by increasing values.

• vbm (float) – the energy level corresponding to the valence band maximum

• cbm (float) – the energy level corresponding to the conduction band minimum

• mu (float) – the electron chemical potential

• T (float or array) – the temperature range where the concentration will be calculated

• dos_down (2D array) –

  If left to None, then the spin-down electrons are considered to have the same dos as dos.

  Note:

  in case the spin down dos is reported using negative numbers, you need to flip the sign before using them in this function.

• units (str) – the units of energy

mu

electron chemical potential

Type

float

T

the temperature for which the Fermi Dirac distribution is calculated

Type

float or 1D numpy array

class spinney.defects.concentration.ConductionElectron(dos, vbm, cbm, mu, T, dos_down=None, energy_units='eV')

Class describing a free electron

get_conduction_electron_number()

Returns the number of electrons in the conduction band as a 1D numpy array of length of T

class spinney.defects.concentration.DefectConcentration(E_form, N, T, g=1, energy_units='eV')

Class describing defect concentrations in the dilute limit at the thermodynamic equilibrium:

\[n_{eq}^d(T) = \frac{N_d g}{e^{E_f(d)/k_B T} + g}\]

Or the approximated expression often used:

\[n_{eq}^d(T) = N_d g e^{(-E_f(d)/k_BT}\]

\(N_d\) is the number of sites available for the specific defect in the crystal.

Parameters

• E_form (float) – The formation energy of the point defect

• N (float) – The number of sites available for the defect per volume unit. This is usually expressed in number of sites per cm^3, or number of sites per cell.

• T (array or float) – The temperatures (in K) to use for calculating the defect concentration

• g (float) – the intrinsic degeneracy of the defect

• energy_units (string) – the employed units of energy

dilute_limit_concentration

the defect concentration at each temperature. If T was a scalar, a scalar is returned.

Type

1D numpy array of length len(T)

dilute_limit_approx_concentration

the defect concentration at each temperature, using the more approximative exponential formula. If T was a scalar, a scalar is returned.

Type

1D numpy array of length len(T)

class spinney.defects.concentration.EquilibriumConcentrations(charge_states, form_energy_vbm, vbm, e_gap, site_conc, dos, T_range, g=None, N_eff=0, units_energy='eV', dos_down=None)

Represents the defects and carriers concentrations for Fermi level values given by the charge neutrality condition in a homogeneous semiconductor:

\[\sum_i q_i n_d(q, T) + p - n = N_d\]

where \(q_i\) is the defect charge state, \(n_d\) its concentration, \(p\) and \(n\) are the concentration of free holes and electrons, respectively and \(N_d\) is the effective doping level.

Parameters

• charge_states (dict of arrays or None) –

  defect_type : list of charge states for all considered defects If equal to None, the pristine semiconductor is considered.

  Example:

  If we consider the Si vacancy and Si interstitials in the charge states -2, -1, 0, 1, 2; then:

  >>> charge_states = {'Si_int':[-2, -1, 0, 1, 2],
  ...                  'Vac_Si':[-2, -1, 0, 1, 2]}
  

• form_energy_vbm (dict of arrays or None) –

  defect_type : the defect formation energy calculated for an electron chemical potential equal to vbm. For any array, the order of the values must be the same as the one used in charge_states

  If equal to None, the pristine semiconductor is considered.

  Example:

  For the defects listed above, we would type:

  >>> form_energy_vbm = {'Si_int':[val_m2, val_m1, val_0,
  ...                              val_1, val_2],
  ...                    'Vac_Si':[valv_m2, valv_m1, valv_0,
  ...                              valv_1, valv_2]}
  

  where val_m2 is the formation energy of the Si interstitial in charge state -2, valv_m2 is the one of the Si vacancy in this charge state and so on.

• vbm (float) – the value of the valence band maximum

• e_gap (float) – the value of the band gap

• site_conc (dict) –

  defect_type : site_concentration for defect_type = ‘electron’ and ‘hole’, this value should be the concentration for the unit cell used to calculate dos.

  Example:

  For the defects listed above, we have:

  >>> site_conc = {'Si_int'  : conc_Si_int,
  ...           'Vac_Si'  : conc_Vac_Si,
  ...           'electron': conc_electrons,
  ...           'hole'    : conc_holes}
  

• dos (2D array) – The first column are the energies of 1-electron level, the second the DOS per cell. The values of vbm and e_gap must be consistent with the dos.

• T_range (array or float) – the temperature range

• g (dict of lists or None) –

  defect_type : [degeneracy charge state 1, degeneracy charge_state 2, …] each list represents the degeneracy for a given type of defect in each of its charge states. The order has to match that of charge_states[defect_type].

  If None, all the degeneracy factors are taken equal to 1. The structure is analogous to charge_states.

• N_eff (float) – the effective-doping concentration

• units_energy (str) – employed energy units

• dos_down (array) – the eventual dos for the spin-down electrons

defect_order

the ordered sequence with the defect names.

Type

tuple

charge_states

defect_name : sequence of charge states of the defect

Type

dict

formation_energies_vbm

defect_name : sequence of the defect formation energies, for every charge state listed in charge_states, for an electron chemical potential equal to vbm

Type

dict

formation_energies_equilibrium

defect_name : sequence of the defect formation energies, for every charge state listed in charge_states, for an electron chemical potential equal to equilibrium_fermi_level

Type

dict

site_conc

defect_name : effective concentration for that kind of defect

Type

dict

T

the temperatures considered for calculating the defect formation energies

Type

numpy 1D array

vbm

the value of the pristine system valence band maximum

Type

float

cbm

the value of the pristine system conduction band minimum

Type

float

N_eff

effective dopant concentration

Type

float

equilibrium_fermi_level

the value of the electron chemical potential at the equilibrium for a given temperature

Type

numpy 1D array of length len(self.T)

equilibrium_defect_concentrations

defect_name : {defect_charge_state : array, …} where array is a numpy 1D array of length len(self.T) holding the calculated defect concentration for that charge state at at given temperature

Type

dict

equilibrium_electron_concentrations

len(self.T) the value of the electron concentration at the equilibrium for a given temperature

Type

numpy 1D array of length

equilibrium_hole_concentrations

the value of the hole concentration at the equilibrium for a given temperature

Type

numpy 1D array of length len(self.T)

equilibrium_carrier_concentrations

len(self.T) the value of the carrier concentration at the equilibrium for a given temperature. If the value is positive, holes are the main carriers; otherwise electrons.

Type

numpy 1D array of length

find_root_algo

specifies which algorithm to use in order to find the roots of the charge neutrality condition. Possible values:

• brentq, newton, bisect

Type

string

get_equilibrium_defect_concentrations()

Return the equilibrium defect concentration as a function of the temperature.

Returns

concentrations

Return type

2D numpy array of shape (len(self._charge_states), len(self._T))

get_equilibrium_electron_concentrations()

Return the equilibrium electron concentration as a function of the temperature.

Returns

concentrations

Return type

1D array of length len(self.T)

get_equilibrium_fermi_level()

Return the equilibrium value of the Fermi level as a function

of the temperature.

Returns

equilibrium_fermi_level – the calculated Fermi level as a function of T

Return type

1D numpy array of length len(self.T)

class spinney.defects.concentration.FermiDiracDistribution(energy, mu, T, energy_units='eV')

Implementation of the Fermi-Dirac distribution:

\[\frac{1}{1 + e^{(E-\mu_e)/k_B T}}\]

Parameters

• energy (float or 1D darray) – the energies of the 1-electron levels

• mu (float) – the electron chemical potential

• T (float or 1D array) – the temperature

• energy_units (string) – the units in which energy and mu are expressed. T is always assumed to be in K.

mu

the electrom chemical potential

Type

float

kb

Boltzman constant in terms of energy_units

Type

float

values

The values of the Fermi distribution for those values of energy and temperature.

Type

numpy 2D array of shape (len(energy), len(T))

class spinney.defects.concentration.ValenceHole(dos, vbm, cbm, mu, T, dos_down=None, energy_units='eV')

Class describing a free hole

get_valence_holes_number()

Returns the number of holes in the valence band as a 1D numpy array of length of T.

spinney.defects.concentration.extract_formation_energies_from_file(file_name)

Read a file with the format:

# arbitrary number of comment lines
defect_name    defect_charge_state   defect_formation_energy

Parameters

file_name (string) – the file with the data

Returns

init_args – the tuple contains two elements: the two dictionaries to be used as the first two arguments for initializing a EquilibriumConcentrations instance

Return type

tuple

General-purpose tools

• spinney.tools.formulas

• spinney.tools.reactions

The formulas module

Module containing functions useful for manipulating chemical formulas

spinney.tools.formulas.count_elements(compound, total=False)

Returns the number of atoms in a chemical system.

Parameters

• compound (str) – the compound of interest

• total (bool) – If True, the total number of atoms in the system is also returned

Returns

el_count – Dictionary: atom:number of atoms for each atom in compound. If total, returns also the number of elements in compound

Return type

dict/tuple

spinney.tools.formulas.get_formula_unit(compound)

Gets the formula unit of a particular compound

Parameters

compound (string) – the formula of the compound.

Returns

formula_unit – the formula unit of compound

Return type

string

Notes

This function automatically reduces the coefficients to the smallest integers. It however preserves the number of symbols in the formula; e.g. C6H6 will return CH, but CH3COOH will return still CH3COOH.

spinney.tools.formulas.get_number_fu(compound, fu=None)

Returns how many formula units fu are present in compound

Parameters

• compound (string) – the formula of the compound

• fu (string) – the formula unit If None, the actual formula unit is used

Returns

no_units – The number of formula units in compound

Return type

float

spinney.tools.formulas.get_stoichiometry(compound, fractional=True)

Given compound, it returns its stoichiometry.

Parameters

• compound (string) – the compound’s formula

• fractional (bool) – if True, for each element is returned its molar fraction; otherwise, it is returned the number of atoms per formula

Returns

elements_count – A dictionary element:coefficient for each element in compound

Return type

dict

The reactions module

Helper functions for calculating reaction energies.

spinney.tools.reactions.calculate_defect_formation_energy(e_defect, e_pristine, chem_potentials, charge_state, E_corr=0)

Calculate the formation energy of a point defect.

Parameters

• e_defect (dict) – {formula : energy of the defective supercell}

• e_pristine (dict) – {formula : energy of the pristine supercell}

• chem_potentials (dict) – {element_name : chemical_potential} element_name is the element symbol, in the case of the electron, element_name = electron

• charge_state (int) – the formal charge state of the defect

• E_corr (float) – eventual corrections to the defect formation energy

Returns

formation_energy

Return type

float

spinney.tools.reactions.calculate_formation_energy_fu(compound_dict, components_dict)

Calculated the formation energy of a given compound per formula unit

Parameters

• compound_dict (dict) – formula : energy for each compound needed to calculate the defect formation energy: formula is a string representing the compound formula, energy is a number, representing the corresponding compound energy.

• components_dict (dict) – formula : energy for each compound needed to calculate the defect formation energy: formula is a string representing the compound formula, energy is a number, representing the corresponding compound energy.

Examples

In order to calculate the formation energy of \(\mathrm{SrTiO}_3\), compound_dict and components_dict are:

>>> compound_dict = {'SrTiO3':Ec}
>>> components_dict = {'Sr':Esr, 'Ti':Eti, 'O2':Eo2}

spinney.tools.reactions.calculate_reaction_energy(reaction, compound_energies)

Calculates the reaction energy of a compound

Parameters

reaction (dict of lists of tuples) –

Describes the reaction that will be calculated. The dictionary keys are ‘reactants’ and ‘products’. The value of ‘reactants’ is a list with tuples. Each tuple contains the FORMULA UNIT of the reactants and the number of moles for that reactant in the reaction. The value of ‘products’ is analogous, but for the products.

The compounds specified in reaction will be taken as the formula units in calculating the reaction.

compound_energiesdict of lists

Analogous structure as in reaction, but instead of the number of moles, the values are the calculated energies.

Note:

the order of the compounds in each list has to match that in reaction!

E.g. we calculated the energies of ‘Mn2’ (bcc Mn cell), ‘O2’ and ‘Mn32O48’ (Pbca space group) Then we would have:

>>> compound_energies = {
                          'reactants' : [('Mn2',E1), ('O2', E2)],
                          'products'  : [('Mn32O48', E3)]
                         }

Returns

energy – The energy of the reaction specified in reaction

Return type

float

Examples

Suppose one is interested in the reaction:

\[2\mathrm{Mn} + 3/2 \mathrm{O}_2 \rightarrow \mathrm{Mn}_2\mathrm{O}_3\]

And one has calculated the energy of Mn2 (bcc cell) and stored it in the variable E1, the energy of O2 is stored in E2 and that of Mn2O3 in E3. Then to calculate the reaction energy one can use:

>>> reaction = {
                'reactants' : [('Mn', 2), ('O2', 3/2)],
                'products'  : [('Mn2O3', 1)]
               }
>>> compound_energies = {
                         'reactants' : [('Mn2', E1), ('O2', E2)],
                         'products'  : [('Mn2O3', E3)]}
>>> calculate_reaction_energy(reaction, compound_energies)

spinney.tools.reactions.get_compound_energy_per_atom(energy, formula)

Returns the energy of the compound per atom (unit_energy/atom)

Parameters

• energy (float) – the energy of the compound

• formula (string) – the formula of the compound

Returns

energy – The energy per atom

Return type

float

spinney.tools.reactions.get_compound_energy_per_formula_unit(energy, formula, formula_unit=None)

Returns the calculated energy of the compound per formula unit (unit_energy/f.u.)

Parameters

• energy (float) – the energy of the compound with formula formula

• formula (string) – the compound formula

• formula_unit (string) – The formula unit of the compound

Returns

energy – The total energy per formula unit and the number of formula units in formula

Return type

float

Support for first-principles codes

• spinney.io.vasp

• spinney.io.wien2k

The io.vasp module

Helper functions for VASP postprocessing

spinney.io.vasp.extract_dos(vasprun_file, save_dos=False)

From the vasprun.xml file extrat the DOS.

Parameters

• vasprun_file (string) – path to the vasprun.xml file

• save_dos (bool) – if True, the extracted DOS are saved as a text file. The first column is the energy (in eV) and the second the DOS (states/cell)

Returns

dos – each element is a 2D numpy array. First column

Return type

2-ple, the DOS up and DOS down (if any)

spinney.io.vasp.extract_potential_at_core_vasp(file)

Read the VASP OUTCAR file and extract the values of the electrostatic potential evaluated at the ions positions.

Parameters

file (string) – path to the OUTCAR file

Returns

result – for each atom in the system, returns the electrostatic potential at the atomic sites

Return type

1D numpy array

The io.wien2k module

Helper functions for WIEN2k

spinney.io.wien2k.average_core_potential_wien2k(potential, r0, rmax)

Calculates the average core potential within a spherical shell with the smallest raius r0 and the largest rmax.

Parameters

• potential (1D numpy array) – the radial part of the potential corresponding to the l=0, m=0 angular component.

• r0 (float) – the smallest radius of the shell

• rmax (float) – the largest radius of the shell

Returns

float

Return type

the averaged potential in the spherical shell

spinney.io.wien2k.extract_potential_at_core_wien2k(struct, vcoul)

Extracts the average electrostatic potential within the atomic spheres for each atom in the system.

Parameters

• struct (string) – the path to the .struct file

• vcoul (string) – the path to the .vcoul file

Returns

result – for each atom in the system, returns the average electrostatic potential within the spherial region

Return type

1D numpy array

spinney.io.wien2k.prepare_ase_atoms_wien2k(struct_file, scf_file)

Prepare an ase.Atoms object compatible with the interface of the PointDefect class in Spinney.

Parameters

• struct_file (string) – the path to the WIEN2K .struct file

• total_energy (float) – the value of the total energy as found in the WIEN2K .scf file

Returns

ase_wien – the ase.Atoms object representing the system

Return type

ase.Atoms

spinney.io.wien2k.read_energy_wien2k(scf_file)

Returns the total electronic energy.

Parameters

scf_file (string) – path to the WIEN2k .scf file

Returns

energy – the energy of the system

Return type

float

spinney.io.wien2k.read_wien2k_radial_data(struct_file)

Reads from a .struct file some information about the atom-centered spheres related to the radial potential within the sphere.

Parameters

struct_file (str) – path to the struct file

Returns

No_atoms is the number of irreducible atoms in the system, for each of them, we store R0, RMT, NPT, MULTI

Return type

2D tuple array of shape (No_atoms, 4)

spinney.io.wien2k.read_wien2k_vcoul(vcoul_file)

Reads the radial part of the electrostatic potential, corresponding to the angular term l=0, m=0, inside the atomic sphere.

Parameters

vcoul_file (str) – the path of the .vcoul file

Returns

result – each list contains the radial potential for one atom

Return type

list of lists