Input and output for Onsager transport calculation

The Onsager calculators currently include two computational approaches to determining transport coefficients: an “interstitial” calculation, and a “vacancy-mediated” calculator. Below we describe the

0. Assumptions used in transport model that are necessary to understand the data to be input, and the limitations of the results;

1. Crystal class setup needed to initiate a calculation;

2. Interstitial calculator setup needed for an single mobile species calculation, or,

3. Vacancy-mediated calculator setup needed for a vacancy-mediated substitutional solute calculation;

4. the creation of VASP-style input files to be run to generate input data;

5. proper Formatting of input data to be compatible with the calculators; and

6. Interpretation of output which includes how to convert output into transport coefficients.

This follows the overall structure of a transport coefficient calculation. Broadly speaking, these are the steps necessary to compute transport coefficients:

1. Identify the crystal to be considered; this requires mapping whatever defects are to be considered mobile onto appropriate Wyckoff sites in the crystal, even if those exact sites are not occupied by true atoms.

2. Generate lists of symmetry unrelated “defect states” and “defect state transitions,” along with the appropriate “calculator object.”

3. Construct input files for total energy calculations to be run outside of the Onsager codebase; extract appropriate energy and frequency information from those runs.

4. Input the data in a format that the calculator can understand, and transform those energies and frequencies into rates at a given temperature assuming Arrhenius behavior.

5. Transform the output into physically relevant quantities (Onsager coefficients, solute diffusivities, mobilities, or drag ratios) with appropriate units.

Assumptions used in Onsager

The Onsager code computes transport of defects on an infinite crystalline lattice. Currently, the code requires that the particular defects can be mapped onto Wyckoff positions in a crystal. This does not require that the defect be an atom occupying various Wyckoff positions (though that obviously is captured), but merely that the defect have the symmetry and transitions that can be equivalently described by an “object” that occupies Wyckoff positions. Simple examples include vacancies, substitutional solutes, simple interstitial atoms, as well as more complex cases such as split vacancy defects (e.g.: a V-O\(_\text{i}\)-V split double vacancy with oxygen interstitial in a closed-packed crystal; the entire defect complex can be mapped on to the Wyckoff position of the oxygen interstitial). In order to calculate diffusion, a few assumptions are made:

• defects are dilute: we never consider more than one defect at a time in an “infinite” periodic crystal; the vacancy-mediated diffuser uses one vacancy and one solute.

• defects diffuse via a Markovian process: defect states are well-defined, and the transition time from state-to-state is much longer than the equilibration time in a state, so that the evolution of the system is described by the Master equation with time-independent rates.

• defects do not alter the underlying symmetry of the crystal: while the defect itself can have a lower symmetry (according to its Wyckoff position), the presence of a defect does not lead to a global phase transformation to a different crystal; moreover, the crystal maintains translational invariance so that the energy of the system with defect(s) is unchanged under translations.

All of these assumptions are usually good: the dilute limit is valid without strong interactions (such as site blocking), Markovian processes are valid as long as barriers are a few times \(k_\text{B}T\), and we are not currently aware of any (simple) defects that induce phase transformations.

Furthermore, relaxation around a defect (or defect cluster) is allowed, but the assumption is that all of the atomic positions can be easily mapped back to “perfect” crystal lattice sites. This is an “off-lattice” model. In some cases, it can be possible to incorporate “new” states, especially metastable states, that are only accessible by a defect.

Finally, the code requires that all diffusion happens on a single sublattice. This sublattice is defined by a single chemical species; it can include multiple Wyckoff positions. But the current algorithms assume that transitions do not result in the creation of antisite defects (where a chemical species is on an “incorrect” sublattice).

Crystal class setup

The assumption of translational invariance of our defects is captured by the use of a Crystal object. Following the standard definition of a crystal, we need to specify (a) three lattice vectors, and (b) at least one basis position, corresponding to at least one site. The crystal needs to contain at least the Wyckoff positions on a single sublattice corresponding to the diffusing defects. It can be useful for it to contain more atoms that act as “spectator” atoms: they do not participate in diffusion, but define both the underlying symmetry of the crystal, and if atomic-scale calculations will be used to compute configuration and transition-state energies, are necessary to define the energy landscape of diffusion.

The lattice vectors of the underlying crystal set the units of length in the transport coefficients. Hence, if the vectors are entered in units of nm, this corresponds to a factor of \(10^{-18}\text{ m}^2\) in the transport coefficients. This should also be considered when including factors of volume per atom as well.

• The lattice vectors are given by three vectors, \(\mathbf{a}_1\), \(\mathbf{a}_2\), \(\mathbf{a}_3\) in Cartesian coordinates. In python, these are input when creating a Crystal either as a list of three numpy vectors, or as a square numpy matrix. Note: if you enter the three vectors as a matrix, remember that it assumes the vectors are column vectors. That is, if amat is the matrix, then amat[:,0] is \(\mathbf{a}_1\), amat[:,1] is \(\mathbf{a}_2\), and amat[:,2] is \(\mathbf{a}_3\). This may not be what you’re expecting. The main recommendation is to enter the lattice vectors as a list (or tuple) of three numpy vectors.

• The atomic basis is given by a list of lists of numpy vectors of positions in unit cell coordinates. For a given basis, then basis[0] is a list of all positions for the first chemical element in the crystal, basis[1] is the second chemical element, and so on. If you only have a single chemical element, you may enter a list of ``numpy`` vectors.

• An optional spin degree of freedom can be included. This is a list of objects, with one for each chemical element. These can be either scalar or vectors, with the assumption that they transform as those objects under group operations. If not included, the spins are all assumed to be equal to 0. Inclusion of these additional degrees of freedom (currently) only impacts the reduction of the unit cell, and the construction of the space group operations.

• We also take in, strictly for bookkeeping purposes, a list of names for the chemical elements. This is an optional input, but recommended for readability.

Once initialized, two main internal operations take place:

1. The unit cell is reduced and optimized. Reduction is a process where we try to find the smallest unit cell representation for the Crystal. This means that the four-atom “simple cubic” unit cell of face-centered cubic can be input, and the code will reduce it to the standard single-atom primitive cell. The reduction algorithm can end up with “unusual” choices of lattice vectors, so we also optimize the lattice vectors so that they are as close to orthogonal as possible, and ordered from smallest to largest. The atomic basis may be shifted uniformly so that if an inversion operation is present, then the inversion center is the origin. Neither choice changes the representation of the crystal; however, the reduction operation can be skipped by including the option noreduce=True.

2. Full symmetry analysis is performed, including: automated construction of space group generator operators, partitioning of basis sites into symmetry related Wyckoff positions, and determination of point group operations for every basis site. All of these operations are automated, and make no reference to crystallographic tables. The algorithm cannot identify which space group it has generated, nor which Wyckoff positions are present. The algorithm respects both chemistry and spin; this also makes spin a useful manipulation tool to artificially lower symmetry for testing purposes as needed.

Note: Crystals can also be constructed by manipulating existing Crystal objects. A useful case is for the interstitial diffuser: when working “interactively,” it is often easier to first make the underlying “spectator” crystal, and then have that Crystal construct the set of Wyckoff positions for a single site in the crystal, and then add that to the basis. Crystal objects are intended to be read-only, so these manipulations result in the creation of a new Crystal object.

A few quick examples:

import numpy as np
import sys
sys.path.extend(['.', '..'])
from onsager import crystal

Face-centered cubic crystal, vacancy-diffusion

Face-centered cubic crystals could be created either by entering the primitive basis:

a0 = 1.
FCCcrys = crystal.Crystal([a0*np.array([0,0.5,0.5]),
                           a0*np.array([0.5,0,0.5]),
                           a0*np.array([0.5,0.5,0])],
                          [np.array([0.,0.,0.])], chemistry=['fcc'])
print(FCCcrys)

#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]

or by entering the simple cubic unit cell with four atoms:

FCCcrys2 = crystal.Crystal(a0*np.eye(3),
                           [np.array([0.,0.,0.]), np.array([0,0.5,0.5]),
                            np.array([0.5,0,0.5]), np.array([0.5,0.5,0])],
                           chemistry=['fcc'])
print(FCCcrys2)

#Lattice:
  a1 = [ 0.5  0.   0.5]
  a2 = [ 0.   0.5  0.5]
  a3 = [-0.5  0.   0.5]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]

The effect of noreduce can be seen by regenerating the FCC crystal using the simple cubic unit cell:

FCCcrys3 = crystal.Crystal(a0*np.eye(3),
                           [np.array([0.,0.,0.]), np.array([0,0.5,0.5]),
                            np.array([0.5,0,0.5]), np.array([0.5,0.5,0])],
                           chemistry=['fcc'], noreduce=True)
print(FCCcrys3)

#Lattice:
  a1 = [ 1.  0.  0.]
  a2 = [ 0.  1.  0.]
  a3 = [ 0.  0.  1.]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]
  (fcc) 0.1 = [ 0.   0.5  0.5]
  (fcc) 0.2 = [ 0.5  0.   0.5]
  (fcc) 0.3 = [ 0.5  0.5  0. ]

Rocksalt crystal, vacancy-diffusion

Two chemical species, with interpenetrating FCC lattices. In MgO, we would allow for V\(_\text{O}\) (oxygen vacancies) to diffuse, with Mg as a “spectator species”:

MgO = crystal.Crystal([a0*np.array([0,0.5,0.5]),
                       a0*np.array([0.5,0,0.5]),
                       a0*np.array([0.5,0.5,0])],
                      [[np.array([0.,0.,0.])], [np.array([0.5,0.5,0.5])]],
                      chemistry=['Mg', 'O'])
print(MgO)

#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (Mg) 0.0 = [ 0.  0.  0.]
  (O) 1.0 = [ 0.5  0.5  0.5]

Face-centered cubic crystal, interstitial diffusion

Interstitials in FCC crystals usually diffuse through a network of octahedral and tetrahedral sites. We can use the Wyckoffpos(u) function in a crystal to generate a list of equivalent sites corresponding to the interstitial positions, and the addbasis() function to create a new crystal with these interstitial sites.

octbasis = FCCcrys.Wyckoffpos(np.array([0.5, 0.5, 0.5]))
tetbasis = FCCcrys.Wyckoffpos(np.array([0.25, 0.25, 0.25]))
FCCcrysint = FCCcrys.addbasis(octbasis + tetbasis, ['int'])
print(octbasis)
print(tetbasis)
print(FCCcrysint)

[array([ 0.5,  0.5,  0.5])]
[array([ 0.75,  0.75,  0.75]), array([ 0.25,  0.25,  0.25])]
#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]
  (int) 1.0 = [ 0.5  0.5  0.5]
  (int) 1.1 = [ 0.75  0.75  0.75]
  (int) 1.2 = [ 0.25  0.25  0.25]

Interstitial calculator setup

The Interstitial calculator is designed for systems where we have a single defect species that diffuses throughout the crystal. This includes single vacancy diffusion, and interstitial solute diffusivity. As for any diffusion calculator, we need to define the configurations that the defect will sample, and the transition states of the defect. In the case of a single defect species,

• configurations are simply the Wyckoff positions of the particular sublattice (specified by a chemistry index);

• transition states are pairs of configurations with a displacement vector that connects the initial to the final system.

We use the sitelist(chemistry) function to construct a list of lists of indices for a given chemistry; the lists of indices are all symmetrically equivalent crystal basis indices, and each list is symmetrically inequivalent: this is a space group partitioning into equivalent Wyckoff positions.

The transition states are stored as a jumpnetwork, which is a list of lists of tuples of transitions: (initial index, final index, deltax) where the indices are self-explanatory, and deltax is a Cartesian vector corresponding to the translation from the initial state to the final state. The transitions in each list is equivalent by symmetry, and the separate lists are symmetrically inequivalent. Note also that reverse transitions are included: (final index, initial index, -deltax). While the jumpnetwork can be constructed “by hand,” it is recommended to use the jumpnetwork() function inside of a crystal to automate the generation, and then remove “spurious” transitions that are identified.

The algorithm in jumpnetwork() is rather simple: a transition is included if

• the distance between the initial and final state is less than a cutoff distance, and

• the line segment between the initial and final state does not come within a minimum distance of other defect states, and

• the line segment between the initial and final state does not come within a minimum distance of any atomic site in the crystal.

The first criterion identifies “close” jumps, while the second criterion eliminates “long” transitions between states when an intermediate configuration may be possible (i.e., \(\text{A}\to\text{B}\) when \(\text{A}\to\text{C}\to\text{B}\) would be more likely as the state C is “close” to the line connecting A to B), and the final criterion elimates transitions that takes the defect too close to a “spectator” atom in the crystal.

The interstitial diffuser also identifies unique tags for all configurations and transition states. The interstitial tags for configurations are strings with i: followed by unit cell coordinates of site to three decimal digits. The interstitial tags for transition states are strings with i: followed by the unit cell coordinates of the initial state, a ^, and the unit cell coordinates of the final state. When one pretty-prints the interstitial diffuser object, the symmetry unique tags are printed. Note that all of the symmetry equivalent tags are stored in the object, and can be used to identify configurations and transition states, and this is the preferred method for indexing, rather than relying on the particular index into the corresponding lists. The interstitial diffuser calculator contains dictionaries that can be used to convert from tags to indices and vice versa.

Finally, YAML interfaces to output the sitelist and jumpnetwork for an interstitial diffuser are includes; combined with the YAML output of the Crystal, this allows for a YAML serialized representation of the diffusion object.

from onsager import OnsagerCalc

Face-centered cubic crystal, vacancy-diffusion

We identify the vacancy sites with the crystal sites in the lattice.

chem = 0
FCCsitelist = FCCcrys.sitelist(chem)
print(FCCsitelist)

[[0]]

chem = 0
FCCjumpnetwork = FCCcrys.jumpnetwork(chem, cutoff=a0*0.78)
for n, jn in enumerate(FCCjumpnetwork):
    print('Jump type {}'.format(n))
    for (i,j), dx in jn:
        print('  {} -> {} dx= {}'.format(i,j,dx))

Jump type 0
  0 -> 0 dx= [ 0.  -0.5 -0.5]
  0 -> 0 dx= [-0.   0.5  0.5]
  0 -> 0 dx= [ 0.5  0.   0.5]
  0 -> 0 dx= [-0.5 -0.  -0.5]
  0 -> 0 dx= [ 0.5  0.  -0.5]
  0 -> 0 dx= [-0.5 -0.   0.5]
  0 -> 0 dx= [ 0.  -0.5  0.5]
  0 -> 0 dx= [-0.   0.5 -0.5]
  0 -> 0 dx= [-0.5  0.5  0. ]
  0 -> 0 dx= [ 0.5 -0.5 -0. ]
  0 -> 0 dx= [-0.5 -0.5  0. ]
  0 -> 0 dx= [ 0.5  0.5 -0. ]

chem = 0
FCCvacancydiffuser = OnsagerCalc.Interstitial(FCCcrys, chem, FCCsitelist, FCCjumpnetwork)
print(FCCvacancydiffuser)

Diffuser for atom 0 (fcc)
#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]
states:
i:+0.000,+0.000,+0.000
transitions:
i:+0.000,+0.000,+0.000^i:-1.000,+0.000,+0.000

Rocksalt crystal, vacancy-diffusion

Two chemical species, with interpenetrating FCC lattices. In MgO, we would allow for V\(_\text{O}\) (oxygen vacancies) to diffuse, with Mg as a “spectator species”.

chem = 1
MgOsitelist = MgO.sitelist(chem)
print(MgOsitelist)

[[0]]

chem = 1
MgOjumpnetwork = MgO.jumpnetwork(chem, cutoff=a0*0.78)
for n, jn in enumerate(MgOjumpnetwork):
    print('Jump type {}'.format(n))
    for (i,j), dx in jn:
        print('  {} -> {} dx= {}'.format(i,j,dx))

Jump type 0
  0 -> 0 dx= [ 0.   0.5 -0.5]
  0 -> 0 dx= [-0.  -0.5  0.5]
  0 -> 0 dx= [-0.5  0.   0.5]
  0 -> 0 dx= [ 0.5 -0.  -0.5]
  0 -> 0 dx= [ 0.  -0.5 -0.5]
  0 -> 0 dx= [-0.   0.5  0.5]
  0 -> 0 dx= [-0.5 -0.5  0. ]
  0 -> 0 dx= [ 0.5  0.5 -0. ]
  0 -> 0 dx= [-0.5  0.  -0.5]
  0 -> 0 dx= [ 0.5 -0.   0.5]
  0 -> 0 dx= [-0.5  0.5  0. ]
  0 -> 0 dx= [ 0.5 -0.5 -0. ]

chem = 1
MgOdiffuser = OnsagerCalc.Interstitial(MgO, chem, MgOsitelist, MgOjumpnetwork)
print(MgOdiffuser)

Diffuser for atom 1 (O)
#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (Mg) 0.0 = [ 0.  0.  0.]
  (O) 1.0 = [ 0.5  0.5  0.5]
states:
i:+0.500,+0.500,+0.500
transitions:
i:+0.500,+0.500,+0.500^i:+0.500,-0.500,+1.500

Face-centered cubic crystal, interstitial diffusion

Interstitials in FCC crystals usually diffuse through a network of octahedral and tetrahedral sites. Nominally, diffusion should occur through an octahedral-tetrahedral jumps, but we can extend the cutoff distance to find additinoal jumps between tetrahedrals.

chem = 1
FCCintsitelist = FCCcrysint.sitelist(chem)
print(FCCintsitelist)

[[0], [1, 2]]

chem = 1
FCCintjumpnetwork = FCCcrysint.jumpnetwork(chem, cutoff=a0*0.51)
for n, jn in enumerate(FCCintjumpnetwork):
    print('Jump type {}'.format(n))
    for (i,j), dx in jn:
        print('  {} -> {} dx= {}'.format(i,j,dx))

Jump type 0
  0 -> 2 dx= [ 0.25 -0.25  0.25]
  2 -> 0 dx= [-0.25  0.25 -0.25]
  0 -> 1 dx= [-0.25  0.25 -0.25]
  1 -> 0 dx= [ 0.25 -0.25  0.25]
  0 -> 2 dx= [-0.25  0.25  0.25]
  2 -> 0 dx= [ 0.25 -0.25 -0.25]
  0 -> 2 dx= [-0.25 -0.25 -0.25]
  2 -> 0 dx= [ 0.25  0.25  0.25]
  0 -> 2 dx= [ 0.25  0.25 -0.25]
  2 -> 0 dx= [-0.25 -0.25  0.25]
  0 -> 1 dx= [ 0.25 -0.25 -0.25]
  1 -> 0 dx= [-0.25  0.25  0.25]
  0 -> 1 dx= [-0.25 -0.25  0.25]
  1 -> 0 dx= [ 0.25  0.25 -0.25]
  0 -> 1 dx= [ 0.25  0.25  0.25]
  1 -> 0 dx= [-0.25 -0.25 -0.25]
Jump type 1
  2 -> 1 dx= [ 0.   0.   0.5]
  1 -> 2 dx= [-0.  -0.  -0.5]
  2 -> 1 dx= [ 0.  -0.5  0. ]
  1 -> 2 dx= [-0.   0.5 -0. ]
  2 -> 1 dx= [-0.5  0.   0. ]
  1 -> 2 dx= [ 0.5 -0.  -0. ]
  2 -> 1 dx= [ 0.   0.5  0. ]
  1 -> 2 dx= [-0.  -0.5 -0. ]
  2 -> 1 dx= [ 0.5  0.   0. ]
  1 -> 2 dx= [-0.5 -0.  -0. ]
  2 -> 1 dx= [ 0.   0.  -0.5]
  1 -> 2 dx= [-0.  -0.   0.5]

chem = 1
FCCintdiffuser = OnsagerCalc.Interstitial(FCCcrysint, chem,
                                          FCCintsitelist, FCCintjumpnetwork)
print(FCCintdiffuser)

Diffuser for atom 1 (int)
#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]
  (int) 1.0 = [ 0.5  0.5  0.5]
  (int) 1.1 = [ 0.75  0.75  0.75]
  (int) 1.2 = [ 0.25  0.25  0.25]
states:
i:+0.500,+0.500,+0.500
i:+0.750,+0.750,+0.750
transitions:
i:+0.500,+0.500,+0.500^i:+0.250,+1.250,+0.250
i:+0.250,+0.250,+0.250^i:+0.750,+0.750,-0.250

The YAML representation is intended to combine both the structural information necessary to construct the (1) crystal, (2) chemistry index of the diffusing defect, (3) sitelist, and (4) jumpnetwork; and the energies, prefactors, and elastic dipoles (derivative of energy with respect to strain) for the symmetry representatives of configurations and jumps. This will become input for the diffuser when computing transport coefficients as a function of temperature, as well as derivatives with respect to strain (elastodiffusion tensor, activation volume tensor).

print(FCCintdiffuser.crys.simpleYAML() +
      'chem: {}\n'.format(FCCintdiffuser.chem) +
      FCCintdiffuser.sitelistYAML(FCCintsitelist) +
      FCCintdiffuser.jumpnetworkYAML(FCCintjumpnetwork))

basis:
- - !numpy.ndarray [0.0, 0.0, 0.0]
- - !numpy.ndarray [0.5, 0.5, 0.5]
  - !numpy.ndarray [0.75, 0.75, 0.75]
  - !numpy.ndarray [0.25, 0.25, 0.25]
chemistry: [fcc, int]
lattice: !numpy.ndarray
- [0.0, 0.5, 0.5]
- [0.5, 0.0, 0.5]
- [0.5, 0.5, 0.0]
lattice_constant: 1.0
spins: null
threshold: 1.0e-08
chem: 1
Dipole:
- !numpy.ndarray
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
- !numpy.ndarray
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
Energy: [0, 0]
Prefactor: [1, 1]
sitelist:
- [0]
- [1, 2]
DipoleT:
- !numpy.ndarray
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
- !numpy.ndarray
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
  - [0.0, 0.0, 0.0]
EnergyT: [0, 0]
PrefactorT: [1, 1]
jumpnetwork:
- - !!python/tuple
    - !!python/tuple [0, 2]
    - !numpy.ndarray [0.25, -0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [2, 0]
    - !numpy.ndarray [-0.25, 0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [0, 1]
    - !numpy.ndarray [-0.25, 0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [1, 0]
    - !numpy.ndarray [0.25, -0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [0, 2]
    - !numpy.ndarray [-0.25, 0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [2, 0]
    - !numpy.ndarray [0.25, -0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [0, 2]
    - !numpy.ndarray [-0.25, -0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [2, 0]
    - !numpy.ndarray [0.25, 0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [0, 2]
    - !numpy.ndarray [0.25, 0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [2, 0]
    - !numpy.ndarray [-0.25, -0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [0, 1]
    - !numpy.ndarray [0.25, -0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [1, 0]
    - !numpy.ndarray [-0.25, 0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [0, 1]
    - !numpy.ndarray [-0.25, -0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [1, 0]
    - !numpy.ndarray [0.25, 0.25, -0.25]
  - !!python/tuple
    - !!python/tuple [0, 1]
    - !numpy.ndarray [0.25, 0.25, 0.25]
  - !!python/tuple
    - !!python/tuple [1, 0]
    - !numpy.ndarray [-0.25, -0.25, -0.25]
- - !!python/tuple
    - !!python/tuple [2, 1]
    - !numpy.ndarray [0.0, 0.0, 0.5]
  - !!python/tuple
    - !!python/tuple [1, 2]
    - !numpy.ndarray [-0.0, -0.0, -0.5]
  - !!python/tuple
    - !!python/tuple [2, 1]
    - !numpy.ndarray [0.0, -0.5, 0.0]
  - !!python/tuple
    - !!python/tuple [1, 2]
    - !numpy.ndarray [-0.0, 0.5, -0.0]
  - !!python/tuple
    - !!python/tuple [2, 1]
    - !numpy.ndarray [-0.5, 0.0, 0.0]
  - !!python/tuple
    - !!python/tuple [1, 2]
    - !numpy.ndarray [0.5, -0.0, -0.0]
  - !!python/tuple
    - !!python/tuple [2, 1]
    - !numpy.ndarray [0.0, 0.5, 0.0]
  - !!python/tuple
    - !!python/tuple [1, 2]
    - !numpy.ndarray [-0.0, -0.5, -0.0]
  - !!python/tuple
    - !!python/tuple [2, 1]
    - !numpy.ndarray [0.5, 0.0, 0.0]
  - !!python/tuple
    - !!python/tuple [1, 2]
    - !numpy.ndarray [-0.5, -0.0, -0.0]
  - !!python/tuple
    - !!python/tuple [2, 1]
    - !numpy.ndarray [0.0, 0.0, -0.5]
  - !!python/tuple
    - !!python/tuple [1, 2]
    - !numpy.ndarray [-0.0, -0.0, 0.5]

Vacancy-mediated calculator setup

For the vacancy mediated diffuser, the configurations and transition states are more complicated. First, we have three types of configurations:

1. Vacancy, sufficiently far away from the solute to have zero interaction energy.

2. Solute, sufficiently far away from the vacancy to have zero interaction energy.

3. Vacancy-solute complexes.

The vacancies and solutes are assumed to be able to occupy the same sites in the crystal, and that neither the vacancy or solute lowers the underlying symmetry of the site. This is a rephrasing of our previous assumption that the symmetry of the defect can be mapped onto the symmetry of the crystal Wyckoff position. There are cases where this is not true: that is, some solutes, when substituted into a crystal, will relax in a way that breaks symmetry. While mathematically this can be treated, we do not currently have an implementation that supports this.

The complexes are only considered out to a finite distance; this is called the “thermodynamic range.” It is defined in terms of “shells,” which is the number of “jumps” from the solute in order to reach the vacancy. We include one more shell out, called the “kinetic range,” which are complexes that include transitions to complexes in the thermodynamic range.

When we consider transition states, we have three types of transition states:

1. Vacancy transitions, sufficiently far away from the solute to have zero interaction energy.

2. Vacancy-solute complex transitions, where only the vacancy changes position (both between complexes in the thermodynamic range, and between the kinetic and thermodynamic range).

3. Vacancy-solute complex transitions, where the vacancy and solute exchange place.

These are called, in the “five-frequency framework”, omega-0, omega-1, and omega-2 jumps, respectively. The five-frequency model technically identifies omega-1 jumps as only between complexes in the thermodynamic range, while the two “additional” jump types, omega-3 and omega-4, connect complexes in the kinetic range to the thermodynamic range. Operationally, we combine omega-1, -3, and -4 into a single set.

To make a diffuser, we need to

1. Identify the sitelist of the vacancies (and hence, solutes),

2. Identify the jumpnetwork of the vacancies

3. Determine the thermodynamic range

then, the diffuser automatically constructs the complexes out to the thermodynamic range, and the full jumpnetworks.

The vacancy-mediated diffuser also identifies unique tags for all configurations and transition states. The tags for configurations are strings with

• v: followed by unit cell coordinates of site to three decimal digits for the vacancy;

• s: followed by unit cell coordinates of site to three decimal digits for the solute;

• s:...-v:... for a solute-vacancy complex.

The transition states are strings with

• omega0: + (initial vacancy configuration) + ^ + (final vacancy configuration);

• omega1: + (initial solute-vacancy configuration) + ^ + (final vacancy configuration);

• omega2: + (initial solute-vacancy configuration) + ^ + (final solute-vacancy configuration).

When one pretty-prints the vacancy-mediated diffuser object, the symmetry unique tags are printed. Note that all of the symmetry equivalent tags are stored in the object, and can be used to identify configurations and transition states, and this is the preferred method for indexing, rather than relying on the particular index into the corresponding lists. The vacancy-mediated diffuser calculator contains dictionaries that can be used to convert from tags to indices and vice versa.

Face-centered cubic crystal, vacancy mediated-diffusion

We construct the Onsager equivalent of the classic five-frequency model. We can use the sitelist and jumpnetwork that we already constructed for the vacancy by itself. Note that the omega-1 list contains four jumps: one that is the normally identified “omega-1”, and three others that correspond to vacancy “escapes” from the first neighbor complex: to the second, third, and fourth neighbors. In the classic five-frequency model, these rates are all forced to be equal.

chem = 0
fivefreqdiffuser = OnsagerCalc.VacancyMediated(FCCcrys, chem,
                                               FCCsitelist, FCCjumpnetwork, 1)
print(fivefreqdiffuser)

Diffuser for atom 0 (fcc), Nthermo=1
#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]
vacancy configurations:
v:+0.000,+0.000,+0.000
solute configurations:
s:+0.000,+0.000,+0.000
solute-vacancy configurations:
s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000
omega0 jumps:
omega0:v:+0.000,+0.000,+0.000^v:-1.000,+0.000,+0.000
omega1 jumps:
omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+0.000,-1.000,+0.000
omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-2.000,+1.000,+0.000
omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,-1.000^v:-1.000,+1.000,-1.000
omega1:s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000^v:-2.000,+0.000,+0.000
omega2 jumps:
omega2:s:+0.000,+0.000,+0.000-v:+1.000,+0.000,+0.000^s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000

An HDF5 representation of the diffusion calculator can be stored for efficient reconstruction of the object, as well as passing between machines. The HDF5 representation includes everything: the underlying Crystal, the sitelist and jumpnetworks, all of the precalculation and analysis needed for diffusion. This greatly speeds up the construction of the calculator.

import h5py
# replace '/dev/null' with your file of choice, and remove backing_store=False
# to read and write to an HDF5 file.
f = h5py.File('/dev/null', 'w', driver='core', backing_store=False)
fivefreqdiffuser.addhdf5(f)  # adds the diffuser to the HDF5 file

# how to read in (after opening `f` as an HDF5 file)
fivefreqcopy = OnsagerCalc.VacancyMediated.loadhdf5(f)  # creates a new diffuser from HDF5
f.close()  # close up the HDF5 file
print(fivefreqcopy)

Diffuser for atom 0 (fcc), Nthermo=1
#Lattice:
  a1 = [ 0.   0.5  0.5]
  a2 = [ 0.5  0.   0.5]
  a3 = [ 0.5  0.5  0. ]
#Basis:
  (fcc) 0.0 = [ 0.  0.  0.]
vacancy configurations:
v:+0.000,+0.000,+0.000
solute configurations:
s:+0.000,+0.000,+0.000
solute-vacancy configurations:
s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000
omega0 jumps:
omega0:v:+0.000,+0.000,+0.000^v:-1.000,+0.000,+0.000
omega1 jumps:
omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+0.000,-1.000,+0.000
omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-2.000,+1.000,+0.000
omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,-1.000^v:-1.000,+1.000,-1.000
omega1:s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000^v:-2.000,+0.000,+0.000
omega2 jumps:
omega2:s:+0.000,+0.000,+0.000-v:+1.000,+0.000,+0.000^s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000

VASP-style input files

At this stage, we have the diffusion “calculator” necessary to compute diffusion, but we need to determine appropriate atomic-scale data to act as input into our calculators. There are two primary steps: (1) constructing appropriate “supercells” containing defect configurations and transition states to be computed, and (2) extracting the appropriate information from those calculations to use in the diffuser. This section deals with the former; the next section will deal with the latter.

The tags are the most straightforward way to identify structures as they are computed, and hence they serve as the mechanism for communicating data into the calculators. To make supercells with defects, we take advantage of the supercell module in Onsager; both calculators contain a makesupercell() function that returns dictionaries of supercells, tags, and appropriate information. Currently, to transform these into usable input files, the automator module can convert such dictionaries into tarballs with an appropriate directory structure, files containing information about appropriate tags for the different configurations, a Makefile that converts CONTCAR output into appropriate POS input for the nudged-elastic band calculation.

Both makesupercell() commands require an input supercell definition, which is a \(3\times3\) integer matrix of column vectors; if N is such a matrix, then the supercell vectors are the columns of A = np.dot(a, N), so that \(\mathbf A_1\) has components N[:,0] in direct coordinates.

from onsager import automator
import tarfile

Face-centered cubic crystal, interstitial diffusion

We will need to construct (and relax) appropriate intersitial sites, and the transition states between them.

help(FCCintdiffuser.makesupercells)

Help on method makesupercells in module onsager.OnsagerCalc:

makesupercells(super_n) method of onsager.OnsagerCalc.Interstitial instance
    Take in a supercell matrix, then generate all of the supercells needed to compute
    site energies and transitions (corresponding to the representatives).

    :param super_n: 3x3 integer matrix to define our supercell
    :return superdict: dictionary of ``states``, ``transitions``, ``transmapping``,
        and ``indices`` that correspond to dictionaries with tags.

        * superdict['states'][i] = supercell of site;
        * superdict['transitions'][n] = (supercell initial, supercell final);
        * superdict['transmapping'][n] = ((site tag, groupop, mapping), (site tag, groupop, mapping))
        * superdict['indices'][tag] = index of tag, where tag is either a state or transition tag.

N = np.array([[-2,2,2],[2,-2,2],[2,2,-2]])  # 32 atom FCC supercell
print(np.dot(FCCcrys.lattice, N))
FCCintsupercells = FCCintdiffuser.makesupercells(N)

[[ 2.  0.  0.]
 [ 0.  2.  0.]
 [ 0.  0.  2.]]

help(automator.supercelltar)

Help on function supercelltar in module onsager.automator:

supercelltar(tar, superdict, filemode=436, directmode=509, timestamp=None, INCARrelax='...', INCARNEB='...', KPOINTS='...', basedir='', statename='relax.', transitionname='neb.', IDformat='{:02d}', JSONdict='tags.json', YAMLdef='supercell.yaml')
    Takes in a tarfile (needs to be open for writing) and a supercelldict (from a
    diffuser) and creates the full directory structure inside the tarfile. Best used in
    a form like

    ::

        with tarfile.open('supercells.tar.gz', mode='w:gz') as tar:
            automator.supercelltar(tar, supercelldict)

    :param tar: tarfile open for writing; may contain other files in advance.
    :param superdict: dictionary of ``states``, ``transitions``, ``transmapping``, ``indices`` that
        correspond to dictionaries with tags; the final tag ``reference`` is the basesupercell
        for calculations without defects.

        * superdict['states'][i] = supercell of state;
        * superdict['transitions'][n] = (supercell initial, supercell final);
        * superdict['transmapping'][n] = ((site tag, groupop, mapping), (site tag, groupop, mapping))
        * superdict['indices'][tag] = (type, index) of tag, where tag is either a state or transition tag; or...
        * superdict['indices'][tag] = index of tag, where tag is either a state or transition tag.
        * superdict['reference'] = (optional) supercell reference, no defects

    :param filemode: mode to use for files (default: 664)
    :param directmode: mode to use for directories (default: 775)
    :param timestamp: UNIX time for files; if None, use current time (default)
    :param INCARrelax: contents of INCAR file to use for relaxation; must contain {system} to be replaced
        by tag value (default: automator.INCARrelax)
    :param INCARNEB: contents of INCAR file to use for NEB; must contain {system} to be replaced
        by tag value (default: automator.INCARNEB)
    :param KPOINTS: contents of KPOINTS file (default: gamma-point only calculation);
        if None or empty, no KPOINTS file at all
    :param basedir: prepended to all files/directories (default: '')
    :param statename: prepended to all state names, before 2 digit number (default: relax.)
    :param transitionname: prepended to all transition names, before 2 digit number  (default: neb.)
    :param IDformat: format for integer tags (default: {:02d})
    :param JSONdict: name of JSON file storing the tags corresponding to each directory (default: tags.json)
    :param YAMLdef: YAML file containing full definition of supercells, relationship, etc. (default: supercell.yaml);
        set to None to not output. **may want to change this to None for the future**

with tarfile.open('io-test-int.tar.gz', mode='w:gz') as tar:
    automator.supercelltar(tar, FCCintsupercells)

tar = tarfile.open('io-test-int.tar.gz', mode='r:gz')

tar.list()

?rw-rw-r-- 0/0        244 2017-07-11 16:14:30 INCAR.relax
?rw-rw-r-- 0/0        305 2017-07-11 16:14:30 INCAR.NEB
?rw-rw-r-- 0/0         31 2017-07-11 16:14:30 KPOINTS
?rwxrwxr-x 0/0       1344 2017-07-11 16:14:30 trans.pl
?rwxrwxr-x 0/0       6283 2017-07-11 16:14:30 nebmake.pl
?rw-rw-r-- 0/0      25975 2017-07-11 16:14:30 Vasp.pm
?rwxrwxr-x 0/0          0 2017-07-11 16:14:30 relax.01/
?rw-rw-r-- 0/0       2267 2017-07-11 16:14:30 relax.01/POSCAR
?rw-rw-r-- 0/0        258 2017-07-11 16:14:30 relax.01/INCAR
?rw-rw-r-- 0/0         35 2017-07-11 16:14:30 relax.01/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 relax.01/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 relax.01/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:30 relax.00/
?rw-rw-r-- 0/0       2267 2017-07-11 16:14:30 relax.00/POSCAR
?rw-rw-r-- 0/0        258 2017-07-11 16:14:30 relax.00/INCAR
?rw-rw-r-- 0/0         35 2017-07-11 16:14:30 relax.00/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 relax.00/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 relax.00/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:30 neb.01/
?rw-rw-r-- 0/0       2298 2017-07-11 16:14:30 neb.01/POS.init
?rw-rw-r-- 0/0       2296 2017-07-11 16:14:30 neb.01/POS.final
?rw-rw-r-- 0/0        342 2017-07-11 16:14:30 neb.01/INCAR
?rw-rw-r-- 0/0         58 2017-07-11 16:14:30 neb.01/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 neb.01/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 neb.01/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:30 neb.00/
?rw-rw-r-- 0/0       2298 2017-07-11 16:14:30 neb.00/POS.init
?rw-rw-r-- 0/0       2296 2017-07-11 16:14:30 neb.00/POS.final
?rw-rw-r-- 0/0        342 2017-07-11 16:14:30 neb.00/INCAR
?rw-rw-r-- 0/0         58 2017-07-11 16:14:30 neb.00/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 neb.00/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:30 neb.00/POTCAR -> ../POTCAR
?rw-rw-r-- 0/0        191 2017-07-11 16:14:30 neb.00/trans.init
?rw-rw-r-- 0/0        191 2017-07-11 16:14:30 neb.00/trans.final
?rw-rw-r-- 0/0        191 2017-07-11 16:14:30 neb.01/trans.init
?rw-rw-r-- 0/0        191 2017-07-11 16:14:30 neb.01/trans.final
?rw-rw-r-- 0/0       1170 2017-07-11 16:14:30 Makefile
?rw-rw-r-- 0/0         14 2017-07-11 16:14:30 relax.00/NEBlist
?rw-rw-r-- 0/0          7 2017-07-11 16:14:30 relax.01/NEBlist
?rw-rw-r-- 0/0        213 2017-07-11 16:14:30 tags.json
?rw-rw-r-- 0/0    1243794 2017-07-11 16:14:30 supercell.yaml

Contents of the Makefile:

with tar.extractfile('Makefile') as f:
    print(f.read().decode('ascii'))

# Makefile to construct NEB input from relaxation output
# we set this so that the makefile doesn't use builtin implicit rules
MAKEFLAGS = -rk

makeneb := "./nebmake.pl"
transform := "./trans.pl"

Nimages ?= 1

.PHONY: help

target := $(foreach neb, $(wildcard neb.*), $(neb)/01/POSCAR)
target: $(target)

help:
        @echo "# Creates input POSCAR for NEB runs, once relaxation runs are complete"
        @echo "# Uses CONTCAR in relaxation directories to create initial run geometry"
        @echo "# environment variable: Nimages (default: $(Nimages))"
        @echo "# target files:"
        @echo $(target) | sed 's/ / /g'
        @echo "# default target: all"

neb.%: neb.%/01/POSCAR neb.%/POSCAR.init neb.%/POSCAR.final

neb.%/01/POSCAR: neb.%/POSCAR.init neb.%/POSCAR.final
        @$(makeneb) $^ $(Nimages)

neb.%/POSCAR.init:
        @$(transform) $^ > $@

neb.%/POSCAR.final:
        @$(transform) $^ > $@

###############################################################
# structure of NEB runs:
neb.00/POSCAR.init: neb.00/trans.init relax.00/CONTCAR
neb.00/POSCAR.final: neb.00/trans.final relax.00/CONTCAR
neb.01/POSCAR.init: neb.01/trans.init relax.01/CONTCAR
neb.01/POSCAR.final: neb.01/trans.final relax.00/CONTCAR

Contents of the tags.json file:

with tar.extractfile('tags.json') as f:
    print(f.read().decode('ascii'))

{
    "neb.00": "i:+0.250,+0.250,+0.250^i:+0.750,+0.750,-0.250",
    "neb.01": "i:+0.500,+0.500,+0.500^i:+0.250,+1.250,+0.250",
    "relax.00": "i:+0.750,+0.750,+0.750",
    "relax.01": "i:+0.500,+0.500,+0.500"
}

Contents of one POSCAR file for relaxation of a configuration:

with tar.extractfile('relax.00/POSCAR') as f:
    print(f.read().decode('ascii'))

i:+0.750,+0.750,+0.750 fcc(32),int_i(1)
1.0
   2.0000000000000000    0.0000000000000000    0.0000000000000000
   0.0000000000000000    2.0000000000000000    0.0000000000000000
   0.0000000000000000    0.0000000000000000    2.0000000000000000
32 1
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.2500000000000000  0.2500000000000000  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.0000000000000000
  0.7500000000000000  0.7500000000000000  0.0000000000000000
  0.2500000000000000  0.0000000000000000  0.2500000000000000
  0.5000000000000000  0.2500000000000000  0.2500000000000000
  0.7500000000000000  0.5000000000000000  0.2500000000000000
  0.0000000000000000  0.7500000000000000  0.2500000000000000
  0.5000000000000000  0.0000000000000000  0.5000000000000000
  0.7500000000000000  0.2500000000000000  0.5000000000000000
  0.0000000000000000  0.5000000000000000  0.5000000000000000
  0.2500000000000000  0.7500000000000000  0.5000000000000000
  0.7500000000000000  0.0000000000000000  0.7500000000000000
  0.0000000000000000  0.2500000000000000  0.7500000000000000
  0.2500000000000000  0.5000000000000000  0.7500000000000000
  0.5000000000000000  0.7500000000000000  0.7500000000000000
  0.0000000000000000  0.2500000000000000  0.2500000000000000
  0.2500000000000000  0.5000000000000000  0.2500000000000000
  0.5000000000000000  0.7500000000000000  0.2500000000000000
  0.7500000000000000  0.0000000000000000  0.2500000000000000
  0.2500000000000000  0.2500000000000000  0.5000000000000000
  0.5000000000000000  0.5000000000000000  0.5000000000000000
  0.7500000000000000  0.7500000000000000  0.5000000000000000
  0.0000000000000000  0.0000000000000000  0.5000000000000000
  0.5000000000000000  0.2500000000000000  0.7500000000000000
  0.7500000000000000  0.5000000000000000  0.7500000000000000
  0.0000000000000000  0.7500000000000000  0.7500000000000000
  0.2500000000000000  0.0000000000000000  0.7500000000000000
  0.7500000000000000  0.2500000000000000  0.0000000000000000
  0.0000000000000000  0.5000000000000000  0.0000000000000000
  0.2500000000000000  0.7500000000000000  0.0000000000000000
  0.5000000000000000  0.0000000000000000  0.0000000000000000
  0.3750000000000000  0.3750000000000000  0.3750000000000000

tar.close()

Face-centered cubic crystal, vacancy mediated-diffusion

We will need to construct (and relax) appropriate vacancy, solute, and solute-vacancy complexes, and the transition states between them. The commands are nearly identical to the interstitial diffuser; the primary difference is the larger number of configurations and files.

help(fivefreqdiffuser.makesupercells)

Help on method makesupercells in module onsager.OnsagerCalc:

makesupercells(super_n) method of onsager.OnsagerCalc.VacancyMediated instance
    Take in a supercell matrix, then generate all of the supercells needed to compute
    site energies and transitions (corresponding to the representatives).

    Note: the states are lone vacancy, lone solute, solute-vacancy complexes in
    our thermodynamic range. Note that there will be escape states are endpoints of
    some omega1 jumps. They are not relaxed, and have no pre-existing tag. They will
    only be output as a single endpoint of an NEB run; there may be symmetry equivalent
    duplicates, as we construct these supercells on an as needed basis.

    We've got a few classes of warnings (from most egregious to least) that can issued
    if the supercell is too small; the analysis will continue despite any warnings:

    1. Thermodynamic shell states map to different states in supercell
    2. Thermodynamic shell states are not unique in supercell (multiplicity)
    3. Kinetic shell states map to different states in supercell
    4. Kinetic shell states are not unique in supercell (multiplicity)

    The lowest level can still be run reliably but runs the risk of errors in escape transition
    barriers. Extreme caution should be used if any of the other warnings are raised.

    :param super_n: 3x3 integer matrix to define our supercell
    :return superdict: dictionary of ``states``, ``transitions``, ``transmapping``,
        ``indices`` that correspond to dictionaries with tags; the final tag
        ``reference`` is the basesupercell for calculations without defects.

        * superdict['states'][i] = supercell of state;
        * superdict['transitions'][n] = (supercell initial, supercell final);
        * superdict['transmapping'][n] = ((site tag, groupop, mapping), (site tag, groupop, mapping))
        * superdict['indices'][tag] = (type, index) of tag, where tag is either a state or transition tag.
        * superdict['reference'] = supercell reference, no defects

N = np.array([[-3,3,3],[3,-3,3],[3,3,-3]])  # 108 atom FCC supercell
print(np.dot(FCCcrys.lattice, N))
fivefreqsupercells = fivefreqdiffuser.makesupercells(N)

[[ 3.  0.  0.]
 [ 0.  3.  0.]
 [ 0.  0.  3.]]

with tarfile.open('io-test-fivefreq.tar.gz', mode='w:gz') as tar:
    automator.supercelltar(tar, fivefreqsupercells)

tar = tarfile.open('io-test-fivefreq.tar.gz', mode='r:gz')

tar.list()

?rw-rw-r-- 0/0        244 2017-07-11 16:14:53 INCAR.relax
?rw-rw-r-- 0/0        305 2017-07-11 16:14:53 INCAR.NEB
?rw-rw-r-- 0/0         31 2017-07-11 16:14:53 KPOINTS
?rwxrwxr-x 0/0       1344 2017-07-11 16:14:53 trans.pl
?rwxrwxr-x 0/0       6283 2017-07-11 16:14:53 nebmake.pl
?rw-rw-r-- 0/0      25975 2017-07-11 16:14:53 Vasp.pm
?rw-rw-r-- 0/0       6844 2017-07-11 16:14:53 POSCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 relax.00/
?rw-rw-r-- 0/0       6784 2017-07-11 16:14:53 relax.00/POSCAR
?rw-rw-r-- 0/0        258 2017-07-11 16:14:53 relax.00/INCAR
?rw-rw-r-- 0/0         35 2017-07-11 16:14:53 relax.00/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 relax.00/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 relax.00/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 relax.02/
?rw-rw-r-- 0/0       6845 2017-07-11 16:14:53 relax.02/POSCAR
?rw-rw-r-- 0/0        258 2017-07-11 16:14:53 relax.02/INCAR
?rw-rw-r-- 0/0         35 2017-07-11 16:14:53 relax.02/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 relax.02/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 relax.02/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 relax.01/
?rw-rw-r-- 0/0       6807 2017-07-11 16:14:53 relax.01/POSCAR
?rw-rw-r-- 0/0        281 2017-07-11 16:14:53 relax.01/INCAR
?rw-rw-r-- 0/0         58 2017-07-11 16:14:53 relax.01/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 relax.01/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 relax.01/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 neb.00/
?rw-rw-r-- 0/0       6822 2017-07-11 16:14:53 neb.00/POS.init
?rw-rw-r-- 0/0       6820 2017-07-11 16:14:53 neb.00/POS.final
?rw-rw-r-- 0/0        349 2017-07-11 16:14:53 neb.00/INCAR
?rw-rw-r-- 0/0         65 2017-07-11 16:14:53 neb.00/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.00/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.00/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 neb.02/
?rw-rw-r-- 0/0       6845 2017-07-11 16:14:53 neb.02/POS.init
?rw-rw-r-- 0/0       6843 2017-07-11 16:14:53 neb.02/POS.final
?rw-rw-r-- 0/0        372 2017-07-11 16:14:53 neb.02/INCAR
?rw-rw-r-- 0/0         88 2017-07-11 16:14:53 neb.02/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.02/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.02/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 neb.04/
?rw-rw-r-- 0/0       6845 2017-07-11 16:14:53 neb.04/POS.init
?rw-rw-r-- 0/0       6843 2017-07-11 16:14:53 neb.04/POSCAR.final
?rw-rw-r-- 0/0        372 2017-07-11 16:14:53 neb.04/INCAR
?rw-rw-r-- 0/0         88 2017-07-11 16:14:53 neb.04/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.04/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.04/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 neb.01/
?rw-rw-r-- 0/0       6845 2017-07-11 16:14:53 neb.01/POS.init
?rw-rw-r-- 0/0       6843 2017-07-11 16:14:53 neb.01/POSCAR.final
?rw-rw-r-- 0/0        372 2017-07-11 16:14:53 neb.01/INCAR
?rw-rw-r-- 0/0         88 2017-07-11 16:14:53 neb.01/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.01/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.01/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 neb.03/
?rw-rw-r-- 0/0       6845 2017-07-11 16:14:53 neb.03/POS.init
?rw-rw-r-- 0/0       6843 2017-07-11 16:14:53 neb.03/POSCAR.final
?rw-rw-r-- 0/0        372 2017-07-11 16:14:53 neb.03/INCAR
?rw-rw-r-- 0/0         88 2017-07-11 16:14:53 neb.03/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.03/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.03/POTCAR -> ../POTCAR
?rwxrwxr-x 0/0          0 2017-07-11 16:14:53 neb.05/
?rw-rw-r-- 0/0       6868 2017-07-11 16:14:53 neb.05/POS.init
?rw-rw-r-- 0/0       6866 2017-07-11 16:14:53 neb.05/POS.final
?rw-rw-r-- 0/0        395 2017-07-11 16:14:53 neb.05/INCAR
?rw-rw-r-- 0/0        111 2017-07-11 16:14:53 neb.05/incar.sed
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.05/KPOINTS -> ../KPOINTS
?rw-rw-r-- 0/0          0 2017-07-11 16:14:53 neb.05/POTCAR -> ../POTCAR
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.00/trans.init
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.00/trans.final
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.01/trans.init
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.02/trans.init
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.02/trans.final
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.03/trans.init
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.04/trans.init
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.05/trans.init
?rw-rw-r-- 0/0        420 2017-07-11 16:14:53 neb.05/trans.final
?rw-rw-r-- 0/0       1447 2017-07-11 16:14:53 Makefile
?rw-rw-r-- 0/0          7 2017-07-11 16:14:53 relax.00/NEBlist
?rw-rw-r-- 0/0         35 2017-07-11 16:14:53 relax.01/NEBlist
?rw-rw-r-- 0/0        710 2017-07-11 16:14:53 tags.json
?rw-rw-r-- 0/0    3458541 2017-07-11 16:14:53 supercell.yaml

Contents of Makefile:

with tar.extractfile('Makefile') as f:
    print(f.read().decode('ascii'))

# Makefile to construct NEB input from relaxation output
# we set this so that the makefile doesn't use builtin implicit rules
MAKEFLAGS = -rk

makeneb := "./nebmake.pl"
transform := "./trans.pl"

Nimages ?= 1

.PHONY: help

target := $(foreach neb, $(wildcard neb.*), $(neb)/01/POSCAR)
target: $(target)

help:
        @echo "# Creates input POSCAR for NEB runs, once relaxation runs are complete"
        @echo "# Uses CONTCAR in relaxation directories to create initial run geometry"
        @echo "# environment variable: Nimages (default: $(Nimages))"
        @echo "# target files:"
        @echo $(target) | sed 's/ / /g'
        @echo "# default target: all"

neb.%: neb.%/01/POSCAR neb.%/POSCAR.init neb.%/POSCAR.final

neb.%/01/POSCAR: neb.%/POSCAR.init neb.%/POSCAR.final
        @$(makeneb) $^ $(Nimages)

neb.%/POSCAR.init:
        @$(transform) $^ > $@

neb.%/POSCAR.final:
        @$(transform) $^ > $@

###############################################################
# structure of NEB runs:
neb.00/POSCAR.init: neb.00/trans.init relax.00/CONTCAR
neb.00/POSCAR.final: neb.00/trans.final relax.00/CONTCAR
neb.01/POSCAR.init: neb.01/trans.init relax.01/CONTCAR
neb.02/POSCAR.init: neb.02/trans.init relax.01/CONTCAR
neb.02/POSCAR.final: neb.02/trans.final relax.01/CONTCAR
neb.03/POSCAR.init: neb.03/trans.init relax.01/CONTCAR
neb.04/POSCAR.init: neb.04/trans.init relax.01/CONTCAR
neb.05/POSCAR.init: neb.05/trans.init relax.01/CONTCAR
neb.05/POSCAR.final: neb.05/trans.final relax.01/CONTCAR

Contents of the tags.json file:

with tar.extractfile('tags.json') as f:
    print(f.read().decode('ascii'))

{
    "neb.00": "omega0:v:+0.000,+0.000,+0.000^v:-1.000,+0.000,+0.000",
    "neb.01": "omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,-1.000^v:-1.000,+1.000,-1.000",
    "neb.02": "omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+0.000,-1.000,+0.000",
    "neb.03": "omega1:s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000^v:-2.000,+0.000,+0.000",
    "neb.04": "omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-2.000,+1.000,+0.000",
    "neb.05": "omega2:s:+0.000,+0.000,+0.000-v:+1.000,+0.000,+0.000^s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000",
    "relax.00": "v:+0.000,+0.000,+0.000",
    "relax.01": "s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000",
    "relax.02": "s:+0.000,+0.000,+0.000"
}

Contents of one POSCAR file for relaxation of a configuration:

with tar.extractfile('relax.01/POSCAR') as f:
    print(f.read().decode('ascii'))

s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000 fcc(106),solute(1)
1.0
   3.0000000000000000    0.0000000000000000    0.0000000000000000
   0.0000000000000000    3.0000000000000000    0.0000000000000000
   0.0000000000000000    0.0000000000000000    3.0000000000000000
106 1
Direct
  0.1666666666666667  0.1666666666666667  0.0000000000000000
  0.3333333333333333  0.3333333333333333  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.0000000000000000
  0.6666666666666666  0.6666666666666666  0.0000000000000000
  0.8333333333333333  0.8333333333333333  0.0000000000000000
  0.1666666666666667  0.0000000000000000  0.1666666666666667
  0.3333333333333333  0.1666666666666667  0.1666666666666667
  0.5000000000000000  0.3333333333333333  0.1666666666666667
  0.6666666666666666  0.5000000000000000  0.1666666666666667
  0.8333333333333333  0.6666666666666666  0.1666666666666667
  0.0000000000000000  0.8333333333333333  0.1666666666666667
  0.3333333333333333  0.0000000000000000  0.3333333333333333
  0.5000000000000000  0.1666666666666667  0.3333333333333333
  0.6666666666666666  0.3333333333333333  0.3333333333333333
  0.8333333333333333  0.5000000000000000  0.3333333333333333
  0.0000000000000000  0.6666666666666666  0.3333333333333333
  0.1666666666666667  0.8333333333333333  0.3333333333333333
  0.5000000000000000  0.0000000000000000  0.5000000000000000
  0.6666666666666666  0.1666666666666667  0.5000000000000000
  0.8333333333333333  0.3333333333333333  0.5000000000000000
  0.0000000000000000  0.5000000000000000  0.5000000000000000
  0.1666666666666667  0.6666666666666666  0.5000000000000000
  0.3333333333333333  0.8333333333333333  0.5000000000000000
  0.6666666666666666  0.0000000000000000  0.6666666666666666
  0.8333333333333333  0.1666666666666667  0.6666666666666666
  0.0000000000000000  0.3333333333333333  0.6666666666666666
  0.1666666666666667  0.5000000000000000  0.6666666666666666
  0.3333333333333333  0.6666666666666666  0.6666666666666666
  0.5000000000000000  0.8333333333333333  0.6666666666666666
  0.8333333333333333  0.0000000000000000  0.8333333333333333
  0.0000000000000000  0.1666666666666667  0.8333333333333333
  0.1666666666666667  0.3333333333333333  0.8333333333333333
  0.3333333333333333  0.5000000000000000  0.8333333333333333
  0.5000000000000000  0.6666666666666666  0.8333333333333333
  0.6666666666666666  0.8333333333333333  0.8333333333333333
  0.0000000000000000  0.1666666666666667  0.1666666666666667
  0.1666666666666667  0.3333333333333333  0.1666666666666667
  0.3333333333333333  0.5000000000000000  0.1666666666666667
  0.5000000000000000  0.6666666666666666  0.1666666666666667
  0.6666666666666666  0.8333333333333333  0.1666666666666667
  0.8333333333333333  0.0000000000000000  0.1666666666666667
  0.1666666666666667  0.1666666666666667  0.3333333333333333
  0.3333333333333333  0.3333333333333333  0.3333333333333333
  0.5000000000000000  0.5000000000000000  0.3333333333333333
  0.6666666666666666  0.6666666666666666  0.3333333333333333
  0.8333333333333333  0.8333333333333333  0.3333333333333333
  0.0000000000000000  0.0000000000000000  0.3333333333333333
  0.3333333333333333  0.1666666666666667  0.5000000000000000
  0.5000000000000000  0.3333333333333333  0.5000000000000000
  0.6666666666666666  0.5000000000000000  0.5000000000000000
  0.8333333333333333  0.6666666666666666  0.5000000000000000
  0.0000000000000000  0.8333333333333333  0.5000000000000000
  0.1666666666666667  0.0000000000000000  0.5000000000000000
  0.5000000000000000  0.1666666666666667  0.6666666666666666
  0.6666666666666666  0.3333333333333333  0.6666666666666666
  0.8333333333333333  0.5000000000000000  0.6666666666666666
  0.0000000000000000  0.6666666666666666  0.6666666666666666
  0.1666666666666667  0.8333333333333333  0.6666666666666666
  0.3333333333333333  0.0000000000000000  0.6666666666666666
  0.6666666666666666  0.1666666666666667  0.8333333333333333
  0.8333333333333333  0.3333333333333333  0.8333333333333333
  0.0000000000000000  0.5000000000000000  0.8333333333333333
  0.1666666666666667  0.6666666666666666  0.8333333333333333
  0.3333333333333333  0.8333333333333333  0.8333333333333333
  0.5000000000000000  0.0000000000000000  0.8333333333333333
  0.0000000000000000  0.3333333333333333  0.0000000000000000
  0.1666666666666667  0.5000000000000000  0.0000000000000000
  0.3333333333333333  0.6666666666666666  0.0000000000000000
  0.5000000000000000  0.8333333333333333  0.0000000000000000
  0.6666666666666666  0.0000000000000000  0.0000000000000000
  0.0000000000000000  0.3333333333333333  0.3333333333333333
  0.1666666666666667  0.5000000000000000  0.3333333333333333
  0.3333333333333333  0.6666666666666666  0.3333333333333333
  0.5000000000000000  0.8333333333333333  0.3333333333333333
  0.6666666666666666  0.0000000000000000  0.3333333333333333
  0.8333333333333333  0.1666666666666667  0.3333333333333333
  0.1666666666666667  0.3333333333333333  0.5000000000000000
  0.3333333333333333  0.5000000000000000  0.5000000000000000
  0.5000000000000000  0.6666666666666666  0.5000000000000000
  0.6666666666666666  0.8333333333333333  0.5000000000000000
  0.8333333333333333  0.0000000000000000  0.5000000000000000
  0.0000000000000000  0.1666666666666667  0.5000000000000000
  0.3333333333333333  0.3333333333333333  0.6666666666666666
  0.5000000000000000  0.5000000000000000  0.6666666666666666
  0.6666666666666666  0.6666666666666666  0.6666666666666666
  0.8333333333333333  0.8333333333333333  0.6666666666666666
  0.0000000000000000  0.0000000000000000  0.6666666666666666
  0.1666666666666667  0.1666666666666667  0.6666666666666666
  0.5000000000000000  0.3333333333333333  0.8333333333333333
  0.6666666666666666  0.5000000000000000  0.8333333333333333
  0.8333333333333333  0.6666666666666666  0.8333333333333333
  0.0000000000000000  0.8333333333333333  0.8333333333333333
  0.1666666666666667  0.0000000000000000  0.8333333333333333
  0.3333333333333333  0.1666666666666667  0.8333333333333333
  0.6666666666666666  0.3333333333333333  0.0000000000000000
  0.8333333333333333  0.5000000000000000  0.0000000000000000
  0.0000000000000000  0.6666666666666666  0.0000000000000000
  0.1666666666666667  0.8333333333333333  0.0000000000000000
  0.3333333333333333  0.0000000000000000  0.0000000000000000
  0.5000000000000000  0.1666666666666667  0.0000000000000000
  0.8333333333333333  0.3333333333333333  0.1666666666666667
  0.0000000000000000  0.5000000000000000  0.1666666666666667
  0.1666666666666667  0.6666666666666666  0.1666666666666667
  0.3333333333333333  0.8333333333333333  0.1666666666666667
  0.5000000000000000  0.0000000000000000  0.1666666666666667
  0.6666666666666666  0.1666666666666667  0.1666666666666667
  0.0000000000000000  0.0000000000000000  0.0000000000000000

tar.close()

Formatting of input data

Once the atomic-scale data from an appropriate total energy calculation is finished, the data needs to be input into formats that the appropriate diffusion calculator can understand. There are some common definitions between the two, but some differences as well.

In all cases, we work with the assumption that our states are thermally occupied, and our rates are Arrhenius. That means that the (relative) probability of any state can be written as

\[\rho = Z^{-1}\rho^0 \exp(-E/k_\text{B}T)\]

for the partition function \(Z\), a site entropic term \(\rho^0 = \exp(S/k_\text{B})\), and energy \(E\). The transition rate from state A to state B is given by

\[\lambda(\text{A}\to\text{B}) = \frac{\nu^\text{T}_{\text{A}-\text{B}}}{\rho^0_A} \exp(-(E^\text{T}_{\text{A}-\text{B}} - E_\text{A})/k_\text{B}T)\]

where \(E^\text{T}_{\text{A}-\text{B}}\) is the energy of the transition state between A and B, and \(\nu^\text{T}_{\text{A}-\text{B}}\) is the prefactor for the transition state.

If we assume harmonic transition state theory, then we can write the site entropic term \(\rho^0\) as

\[\rho^0 = \frac{\prod \nu^{\text{perfect-supercell}}}{\prod \nu^{\text{defect-supercell}}}\]

where \(\nu\) are the vibrational eigenvalues of the corresponding supercells, and the prefactor for the transition state is

\[\nu^\text{T} = \frac{\prod \nu^{\text{perfect-supercell}}}{\prod_{\nu^2>0} \nu^{\text{transition state}}}\]

where we take the product over the real vibrational frequencies in the transition state (there should be one imaginary mode). From a practical point of view, the perfect-supercell cancels out; we will often set \(\rho^0\) to 1 for a single state (so that the other \(\rho^0\) are relative probabilities), and then \(\nu^\text{T}\) becomes more similar to the attempt frequency for the particular jumps. The definitions above map most simply onto a “hopping atom” approximation for the jump rates: the \(3\times3\) force-constant matrix is computed for the atom that is moving in the transition, and its eigenvalues are used to determine the modes \(\nu\).

Note the units: \(\rho^0\) is unitless, while \(\nu^\text{T}\) has units of inverse time; this means that the inverse time unit in the computed transport coefficients will come from \(\nu^\text{T}\) values. If they are entered in THz, that contributes \(10^{12}\text{ s}^{-1}\).

Because we normalize our probabilities, our energies and transition state energies are relative to each other. In all of our calculations, we will multiply energies by \(\beta=(k_\text{B}T)^{-1}\) to get a unitless values as inputs for our diffusion calculators. This means that the diffusers do not have direct information about temperature; explicit temperature factors that appear in the Onsager coefficients must be included by hand from the output transport coefficients. It also means that the calculators do not have a “unit” of energy; rather, \(k_\text{B}T\) and the energies must be in the same units.

Face-centered cubic crystal, interstitial diffusion

We need to compute prefactors and energies for our interstitial diffuser. We can also include information about elastic dipoles (derivatives of energy with respect to strain) in order to compute derivatives of diffusivity with respect to strain (elastodiffusion).

help(FCCintdiffuser.diffusivity)

Help on method diffusivity in module onsager.OnsagerCalc:

diffusivity(pre, betaene, preT, betaeneT, CalcDeriv=False) method of onsager.OnsagerCalc.Interstitial instance
    Computes the diffusivity for our element given prefactors and energies/kB T.
    Also returns the negative derivative of diffusivity with respect to beta (used to compute
    the activation barrier tensor) if CalcDeriv = True
    The input list order corresponds to the sitelist and jumpnetwork

    :param pre: list of prefactors for unique sites
    :param betaene: list of site energies divided by kB T
    :param preT: list of prefactors for transition states
    :param betaeneT: list of transition state energies divided by kB T
    :return D[3,3]: diffusivity as a 3x3 tensor
    :return DE[3,3]: diffusivity times activation barrier (if CalcDeriv == True)

The ordering in the lists pre, beteene, preT and betaeneT corresponds to the sitelist and jumpnetwork lists. The tags can be used to determine the proper indices. The most straightforward way to store this in python is a dictionary, where the key is the tag, and the value is a list of [prefactor, energy]. The advantage of this is that it can be easily transformed to and from JSON for simple serialization.

To see a full list of all tags in the dictionary, the tags member of a diffuser gives a dictionary of all tags, ordered to match the structure of sitelist and jumpnetwork.

FCCintdiffuser.tags

{'states': [['i:+0.500,+0.500,+0.500'],
  ['i:+0.750,+0.750,+0.750', 'i:+0.250,+0.250,+0.250']],
 'transitions': [['i:+0.500,+0.500,+0.500^i:+0.250,+1.250,+0.250',
   'i:+0.250,+0.250,+0.250^i:+0.500,-0.500,+0.500',
   'i:+0.500,+0.500,+0.500^i:+0.750,-0.250,+0.750',
   'i:+0.750,+0.750,+0.750^i:+0.500,+1.500,+0.500',
   'i:+0.500,+0.500,+0.500^i:+1.250,+0.250,+0.250',
   'i:+0.250,+0.250,+0.250^i:-0.500,+0.500,+0.500',
   'i:+0.500,+0.500,+0.500^i:+0.250,+0.250,+0.250',
   'i:+0.250,+0.250,+0.250^i:+0.500,+0.500,+0.500',
   'i:+0.500,+0.500,+0.500^i:+0.250,+0.250,+1.250',
   'i:+0.250,+0.250,+0.250^i:+0.500,+0.500,-0.500',
   'i:+0.500,+0.500,+0.500^i:-0.250,+0.750,+0.750',
   'i:+0.750,+0.750,+0.750^i:+1.500,+0.500,+0.500',
   'i:+0.500,+0.500,+0.500^i:+0.750,+0.750,-0.250',
   'i:+0.750,+0.750,+0.750^i:+0.500,+0.500,+1.500',
   'i:+0.500,+0.500,+0.500^i:+0.750,+0.750,+0.750',
   'i:+0.750,+0.750,+0.750^i:+0.500,+0.500,+0.500'],
  ['i:+0.250,+0.250,+0.250^i:+0.750,+0.750,-0.250',
   'i:+0.750,+0.750,+0.750^i:+0.250,+0.250,+1.250',
   'i:+0.250,+0.250,+0.250^i:-0.250,+0.750,-0.250',
   'i:+0.750,+0.750,+0.750^i:+1.250,+0.250,+1.250',
   'i:+0.250,+0.250,+0.250^i:+0.750,-0.250,-0.250',
   'i:+0.750,+0.750,+0.750^i:+0.250,+1.250,+1.250',
   'i:+0.250,+0.250,+0.250^i:+0.750,-0.250,+0.750',
   'i:+0.750,+0.750,+0.750^i:+0.250,+1.250,+0.250',
   'i:+0.250,+0.250,+0.250^i:-0.250,+0.750,+0.750',
   'i:+0.750,+0.750,+0.750^i:+1.250,+0.250,+0.250',
   'i:+0.250,+0.250,+0.250^i:-0.250,-0.250,+0.750',
   'i:+0.750,+0.750,+0.750^i:+1.250,+1.250,+0.250']]}

In this example, the energy of the octahedral site is 0, with a base prefactor of 1. The tetrahedral site has an energy of 0.5 (eV) above, with a higher relative vibrational degeneracy of 2. The transition state energy from octahedral to tetrahedral is 1.0 (eV) with a prefactor of 10 (THz); and the transition state energy from tetrahedral to tetrahedral is 2.0 (eV) with a prefactor of 50 (THz).

FCCintdata = {
    'i:+0.500,+0.500,+0.500': [1., 0.],
    'i:+0.750,+0.750,+0.750': [2., 0.5],
    'i:+0.500,+0.500,+0.500^i:+0.750,+0.750,-0.250': [10., 1.0],
    'i:+0.750,+0.750,+0.750^i:+1.250,+1.250,+0.250': [50., 2.0]
}

# Conversion from dictionary to lists for a given kBT
# We go through the tags in order, and find one in our data set.
kBT = 0.25  # eV; a rather high temperature
pre = [FCCintdata[t][0] for taglist in FCCintdiffuser.tags['states']
       for t in taglist if t in FCCintdata]
betaene = [FCCintdata[t][1]/kBT for taglist in FCCintdiffuser.tags['states']
           for t in taglist if t in FCCintdata]
preT = [FCCintdata[t][0] for taglist in FCCintdiffuser.tags['transitions']
        for t in taglist if t in FCCintdata]
betaeneT = [FCCintdata[t][1]/kBT for taglist in FCCintdiffuser.tags['transitions']
            for t in taglist if t in FCCintdata]
print(pre,betaene,preT,betaeneT,sep='\n')

[1.0, 2.0]
[0.0, 2.0]
[10.0, 50.0]
[4.0, 8.0]

DFCCint, dDFCCint = FCCintdiffuser.diffusivity(pre, betaene, preT, betaeneT, CalcDeriv=True)
print(DFCCint, dDFCCint, sep='\n')

[[  6.48557013e-02   1.30104261e-18  -1.30104261e-18]
 [  1.30104261e-18   6.48557013e-02   1.30104261e-18]
 [ -1.30104261e-18   1.30104261e-18   6.48557013e-02]]
[[  2.35630632e-01   3.46944695e-18  -3.46944695e-18]
 [  3.46944695e-18   2.35630632e-01   3.46944695e-18]
 [ -3.46944695e-18   3.46944695e-18   2.35630632e-01]]

The interpretation of this output will be described below.

Face-centered cubic crystal, vacancy mediated-diffusion

We will need to compute prefactors and energies for our vacancy, solute, and solute-vacancy complexes, and the transition states between them. The difference compared with the interstitial case is that complex prefactors and energies are excess quantities. That means for a complex, its \(\rho^0\) is the product of \(\rho^0\) for the solute state, the vacancy state, and the excess; the energy \(E\) is the sum of the energy of the solute state, the vacancy state, and the excess. However for the transition states, the prefactors and energies are “absolute”.

help(fivefreqdiffuser.Lij)

Help on method Lij in module onsager.OnsagerCalc:

Lij(bFV, bFS, bFSV, bFT0, bFT1, bFT2, large_om2=100000000.0) method of onsager.OnsagerCalc.VacancyMediated instance
    Calculates the transport coefficients: L0vv, Lss, Lsv, L1vv from the scaled free energies.
    The Green function entries are calculated from the omega0 info. As this is the most
    time-consuming part of the calculation, we cache these values with a dictionary
    and hash function.

    :param bFV[NWyckoff]: beta*eneV - ln(preV) (relative to minimum value)
    :param bFS[NWyckoff]: beta*eneS - ln(preS) (relative to minimum value)
    :param bFSV[Nthermo]: beta*eneSV - ln(preSV) (excess)
    :param bFT0[Nomega0]: beta*eneT0 - ln(preT0) (relative to minimum value of bFV)
    :param bFT1[Nomega1]: beta*eneT1 - ln(preT1) (relative to minimum value of bFV + bFS)
    :param bFT2[Nomega2]: beta*eneT2 - ln(preT2) (relative to minimum value of bFV + bFS)
    :param large_om2: threshold for changing treatment of omega2 contributions (default: 10^8)
    :return Lvv[3, 3]: vacancy-vacancy; needs to be multiplied by cv/kBT
    :return Lss[3, 3]: solute-solute; needs to be multiplied by cv*cs/kBT
    :return Lsv[3, 3]: solute-vacancy; needs to be multiplied by cv*cs/kBT
    :return Lvv1[3, 3]: vacancy-vacancy correction due to solute; needs to be multiplied by cv*cs/kBT

The vacancy-mediated diffuser expects combined \(\beta F := (E - TS)/k_\text{B}T\), so that our probabilities and rates are proportional to \(\exp(-\beta F)\). This is complicated to directly construct, so we have the intermediate function preene2betafree(), which is best used by feeding a dictionary of arrays:

help(fivefreqdiffuser.preene2betafree)

Help on function preene2betafree in module onsager.OnsagerCalc:

preene2betafree(kT, preV, eneV, preS, eneS, preSV, eneSV, preT0, eneT0, preT1, eneT1, preT2, eneT2, **ignoredextraarguments)
    Read in a series of prefactors (exp(S/k_B)) and energies, and return
    beta F for energies and transition state energies. Used to provide scaled values
    to Lij().
    Can specify all of the entries using a dictionary; e.g., ``preene2betafree(kT, **data_dict)``
    and then send that output as input to Lij: ``Lij(*preene2betafree(kT, **data_dict))``
    (we ignore extra arguments so that a dictionary including additional entries can be passed)

    :param kT: temperature times Boltzmann's constant kB
    :param preV: prefactor for vacancy formation (prod of inverse vibrational frequencies)
    :param eneV: vacancy formation energy
    :param preS: prefactor for solute formation (prod of inverse vibrational frequencies)
    :param eneS: solute formation energy
    :param preSV: excess prefactor for solute-vacancy binding
    :param eneSV: solute-vacancy binding energy
    :param preT0: prefactor for vacancy transition state
    :param eneT0: energy for vacancy transition state (relative to eneV)
    :param preT1: prefactor for vacancy swing transition state
    :param eneT1: energy for vacancy swing transition state (relative to eneV + eneS + eneSV)
    :param preT2: prefactor for vacancy exchange transition state
    :param eneT2: energy for vacancy exchange transition state (relative to eneV + eneS + eneSV)
    :return bFV: beta*eneV - ln(preV) (relative to minimum value)
    :return bFS: beta*eneS - ln(preS) (relative to minimum value)
    :return bFSV: beta*eneSV - ln(preSV) (excess)
    :return bFT0: beta*eneT0 - ln(preT0) (relative to minimum value of bFV)
    :return bFT1: beta*eneT1 - ln(preT1) (relative to minimum value of bFV + bFS)
    :return bFT2: beta*eneT2 - ln(preT2) (relative to minimum value of bFV + bFS)

Even this is a bit complicated; so we use an additional function that maps the tags into the appropriate lists, tags2preene():

help(fivefreqdiffuser.tags2preene)

Help on method tags2preene in module onsager.OnsagerCalc:

tags2preene(usertagdict, VERBOSE=False) method of onsager.OnsagerCalc.VacancyMediated instance
    Generates energies and prefactors based on a dictionary of tags.

    :param usertagdict: dictionary where the keys are tags, and the values are tuples: (pre, ene)
    :param VERBOSE: (optional) if True, also return a dictionary of missing tags, duplicate tags, and bad tags
    :return thermodict: dictionary of ene's and pre's corresponding to usertagdict
    :return missingdict: dictionary with keys corresponding to tag types, and the values are
      lists of lists of symmetry equivalent tags that are missing
    :return duplicatelist: list of lists of tags in usertagdict that are (symmetry) duplicates
    :return badtaglist: list of all tags in usertagdict that aren't found in our dictionary

In this example, we have a vacancy-solute binding energy of -0.25 (eV), a vacancy jump barrier of 1.0 (eV) with a prefactor of 10 (THz), an “omega-1” activation barrier of 0.75 (eV) which is a transition state energy of 0.75-0.25 = 0.5, an omega-2 activation barrier of 0.5 (eV) which is a transition state energy of 0.5-0.25 = 0.25, and all of the “omega-3/-4” escape jumps with a transition state energy of 1-0.25/2 = 0.875 (eV).

fivefreqdata = {
    'v:+0.000,+0.000,+0.000': [1., 0.],
    's:+0.000,+0.000,+0.000': [1., 0.],
    's:+0.000,+0.000,+0.000-v:+0.000,-1.000,+0.000': [1., -0.25],
    'omega0:v:+0.000,+0.000,+0.000^v:+0.000,+1.000,+0.000': [10., 1.],
    'omega1:s:+0.000,+0.000,+0.000-v:+1.000,+0.000,-1.000^v:+1.000,+1.000,-1.000': [10., 0.5],
    'omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+1.000,+0.000,+0.000': [20., 0.875],
    'omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-1.000,+2.000,+0.000': [20., 0.875],
    'omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,+0.000^v:+0.000,+2.000,+0.000': [20., 0.875],
    'omega2:s:+0.000,+0.000,+0.000-v:+0.000,-1.000,+0.000^s:+0.000,+0.000,+0.000-v:+0.000,+1.000,+0.000':
    [10., 0.25]
}

# Conversion from dictionary to lists for a given kBT
# note that we can nest the mapping functions.
kBT = 0.25  # eV; a rather high temperature
fivefreqpreene = fivefreqdiffuser.tags2preene(fivefreqdata)
fivefreqbetaF = fivefreqdiffuser.preene2betafree(kBT, **fivefreqpreene)
L0vv, Lss, Lsv, L1vv = fivefreqdiffuser.Lij(*fivefreqbetaF)
print(L0vv, Lss, Lsv, L1vv, sep='\n')

[[ 0.18315639 -0.         -0.        ]
 [-0.          0.18315639 -0.        ]
 [-0.         -0.          0.18315639]]
[[ 1.34561077  0.          0.        ]
 [ 0.          1.34561077  0.        ]
 [ 0.          0.          1.34561077]]
[[-0.40965028  0.          0.        ]
 [ 0.         -0.40965028  0.        ]
 [ 0.          0.         -0.40965028]]
[[ 6.24182702  0.          0.        ]
 [ 0.          6.24182702  0.        ]
 [ 0.          0.          6.24182702]]

The interpretation of this output will be described below.

Interpretation of output

The final step is to take the output from the diffuser calculator, and convert this into physical quantities: solute diffusivity, elastodiffusivity, Onsager coefficients, drag ratios, and so on.

There are two underlying definitions that we use to define our transport coefficients:

\[\mathbf j = -\underline D \nabla c\]

defines the solute diffusivity as the tensorial transport coefficient that relates defect concentration gradients to defect fluxes, and

\[\mathbf j^\text{s} = -\underline L^\text{ss}\nabla\mu^\text{s} - \underline L^\text{sv}\nabla\mu^\text{v}\]

\[\mathbf j^\text{v} = -\underline L^\text{vv}\nabla\mu^\text{v} - \underline L^\text{sv}\nabla\mu^\text{s}\]

defines the Onsager coefficients as the tensorial transport coefficients that relate solute and vacancy chemical potential gradients to solute and vacancy fluxes. We use these equation to also define the units of our transport coefficients. Fluxes are in units of (number)/area/time, so with concentration in (number)/volume, diffusivity has units of area/time. If the chemical potential is written in units of energy, the Onsager coefficients have units of (number)/length/energy/time. If the chemical potentials will instead have units of energy/volume, then the corresponding Onsager coefficients have units of area/energy/time.

Below are more specific details about the different calculators and the output available.

Interstitial diffusivity

The interstitial diffuser outputs a diffusivity tensor that has the units of squared length based on the lengths in the corresponding Crystal, and inverse time units corresponding to the rates that are given as input: the ratio of transition state prefactors to configuration prefactors. In a crystalline system, it is typical to specify the lattice vectors in either nm (\(10^{-9}\text{ m}\)) or Å (\(10^{-10}\text{ m}\)), and the prefactors of rates are often THz (\(10^{12}\text{ s}\)), while diffusivity is often reported in either \(\text{m}^2/\text{s}\) or \(\text{cm}^2/\text{s}\). The conversion factors are

\[1\text{ nm}^2\cdot\text{THz} = 10^{-6}\text{ m}^2/\text{s} = 10^{-2}\text{ cm}^2/\text{s}\]

\[1\text{ A}^2\cdot\text{THz} = 10^{-8}\text{ m}^2/\text{s} = 10^{-4}\text{ cm}^2/\text{s}\]

It it worth noting that this model of diffusion assumes that the “interstitial” form of the defect is its ground state configuration (or at least one of the configurations used in the derivation of the diffusivity is a ground state configuration). This is generally the case for the diffusion of a vacancy, or light interstitial elements; however, the are materials where a solute has a lower energy as a substitutional defect, but can occupy an interstitial site and diffuse from there. This requires knowledge of the relative occupancy of the two states. Using Kroger-Vink notation, let [B] be the total solute concentration, and \([\text{B}_\text{A}]\) and \([\text{B}_\text{i}]\) the substitutional and interstitial concentrations, then

\[D_\text{B} = \left\{[\text{B}_\text{i}]D_\text{int} + [\text{B}_\text{A}]D_\text{sub}\right\}/[\text{B}]\]

for interstitial diffusivity \(D_\text{int}\) and substitutional diffusivity \(D_\text{sub}\). The relative occupancies may be determined by global thermal equilibrium or local thermal equilibrium. The latter is more complex, and relies on knowledge of local defect processes and conditions, and is not discussed further here. For global thermal equilibrium, if we know the energy of the ground state substitutional defect \(E_\text{sub}\) and the lowest energy configuration used by the diffuser \(E_\text{int}\), then

\[[\text{B}_\text{i}]/[\text{B}] = (1 + \exp((E_\text{int}-E_\text{sub})/k_\text{B}T)^{-1} \approx \exp(-(E_\text{int}-E_\text{sub})/k_\text{B}T)\]

and

\[[\text{B}_\text{A}]/[\text{B}] = (1 + \exp(-(E_\text{int}-E_\text{sub})/k_\text{B}T)^{-1} \approx 1\]

where the approximations are valid when \(E_\text{int}-E_\text{sub}\gg k_\text{B}T\).

Derivatives of diffusivity: activation barrier tensor

At any given temperature, the temperature dependence of the diffusivity can be taken as an Arrhenius form,

\[\underline D = \underline D_0 \exp(-\beta \underline E^\text{act})\]

for inverse temperature \(\beta = (k_\text{B}T)^{-1}\), and the activation barrier, \(\underline E^\text{act}\) can also display anisotropy. Note that in this expression, the exponential is taken on a per-component basis, not as a true tensor exponential. We can compute \(Q\) by taking the per-component logarithmic derivative with respect to inverse temperature,

\[\underline E^\text{act} = -\underline D^{-1/2}\frac{d\underline D}{d\beta}\underline D^{-1/2}\]

The diffusivity() function with CalcDeriv=True returns a second tensorial quantity, dD which when multiplied by \(k_\text{B}T\), gives \(d\underline D/d\beta\). Hence, to compute the activation barrier tensor, we evaluate:

print(np.dot(np.linalg.inv(DFCCint), kBT*dDFCCint))

[[  9.08288044e-01  -4.84706356e-18   4.84706356e-18]
 [ -4.84706356e-18   9.08288044e-01  -4.84706356e-18]
 [  4.84706356e-18  -4.84706356e-18   9.08288044e-01]]

In this case, as the matrices are isotropic, we can use \(\underline D^{-1}\) rather than \(\underline D^{-1/2}\) which must be computed via diagonalization.

This tensor has the same energy units as the variable kBT.

Given the barriers for diffusion, one might have expected that \(\underline E^\text{act}\) would be 1, as that is the transition state energy to go from octahedral to tetrahedral. However, the activation barrier is approximately the rate-limiting transition state energy minus the average configuration energy. Since we’ve chosen a large temperature, the tetrahedral sites have non-negligible occupation, which raises the average energy. As the temperature decreases, the activation energy will approach 1.

Derivatives of diffusivity: elastodiffusion and activation volume tensor

The derivative with respect to strain is the fourth-rank elastodiffusivity tensor \(\underline d\), where

\[d_{abcd} = \frac{dD_{ab}}{d\varepsilon_{cd}}\]

This is returned by the elastodiffusion function, which requires the elastic dipole tensors be included in the function call as well. The elastic dipoles have the same units of energies, and so are input as \(\beta\underline P\), which is unitless. The returned tensor has the same units as the diffusivity.

The activation volume tensor (logarithmic derivative of diffusivity with respect to stress) can be computed from the elastodiffusivity tensor if the compliance tensor \(\underline S\) is known; then,

\[V^\text{act}_{abcd} = k_\text{B}T \sum_{ijkl=1}^3 (\underline D^{-1/2})_{ai} d_{ijkl} (\underline D^{-1/2})_{bj} S_{klcd}\]

The units of this quantity are given by the units of \(k_\text{B}T\) (energy) multiplied by the units of \(\underline S\) (inverse pressure). Typically, \(k_\text{B}T\) will be known in eV and \(\underline S\) in GPa\(^{-1}\), so the conversion factor

\[1\text{ eV}\cdot\text{GPa}^{-1} = 1.6022\times10^{-19}\text{ J}\cdot10^{-9}\text{ m}^3/\text{J} = 0.16022\text{ nm}^3 = 160.22\text{ A}^3\]

can be useful.

Vacancy-mediated diffusivity

The interstitial diffuser outputs a diffusivity tensor that has the units of squared length based on the lengths in the corresponding Crystal, and inverse time units corresponding to the rates that are given as input: the ratio of transition state prefactors to configuration prefactors. In a crystalline system, it is typical to specify the lattice vectors in either nm (\(10^{-9}\text{ m}\)) or Å (\(10^{-10}\text{ m}\)), and the prefactors of rates are often THz (\(10^{12}\text{ s}\)). The quantities L0vv, Lss, Lsv, and L1vv output by the Lij function all have the units of area/time, so the the conversion factors below are often useful:

\[1\text{ nm}^2\cdot\text{THz} = 10^{-6}\text{ m}^2/\text{s} = 10^{-2}\text{ cm}^2/\text{s}\]

\[1\text{ A}^2\cdot\text{THz} = 10^{-8}\text{ m}^2/\text{s} = 10^{-4}\text{ cm}^2/\text{s}\]

To convert the four quantities into \(\underline L^\text{vv}\), \(\underline L^\text{ss}\), and \(\underline L^\text{sv}\), some additional information is required.

First, in the dilute limit, \(\underline L^\text{ss}\) and \(\underline L^\text{sv}\) are proportional to \((k_\text{B}T)^{-1}c^\text{v}c^\text{s}\); none of these quantities are known to the diffuser, and the two concentrations are essentially independent variables that must be supplied. The concentrations in these cases are fractional concentrations, not per volume. Finally, if the Onsager coefficients are for chemical potential specified as energies (not energies per volume), the quantities need to be divided by the volume per atom, and the final quantity has the appropriate units. Hence,

• \(\underline L^\text{ss}\) = Lss*(solute concentration)*(vacancy concentration)/(volume)/kBT

• \(\underline L^\text{sv}\) = Lsv*(solute concentration)*(vacancy concentration)/(volume)/kBT

where the concentration quantities are fractional.

The vacancy \(\underline L^\text{vv}\) is more complicated, as it has a leading order term that is independent of solute, and a first order correction that is linear in the solute concentration. Hence,

• \(\underline L^\text{vv}\) = (L0vv + L1vv*(solute concentration))*(vacancy concentration)/(volume)/kBT

Drag ratio

The drag ratio is the unitless (tensorial) quantity \(\underline L^\text{sv}(\underline L^\text{ss})^{-1}\). Because of the identical prefactors in front of both terms in the dilute limit, this is given by

print(np.dot(Lsv, np.linalg.inv(Lss)))

[[-0.30443445  0.          0.        ]
 [ 0.         -0.30443445  0.        ]
 [ 0.          0.         -0.30443445]]

The vacancy wind factor \(G=\underline L^\text{As}(\underline L^\text{ss})^{-1}\) is related to the drag ratio by simple transformations.

Solute diffusivity in the dilute limit

The solute diffusivity can also be computed for the dilute limit as well. The general relation between \(\underline D^\text{s}\) and the Onsager transport coefficients is

\[\underline D^\text{s} = k_\text{B}T\Omega\left\{(c^\text{s})^{-1}\underline L^\text{ss} - (1-c^\text{s}-c^\text{v})^{-1}\underline L^\text{As}\right\}\left(1+\frac{d\ln\gamma^\text{s}}{d\ln c^\text{s}}\right)\]

where \(\Omega\) is the volume per atom and \(\gamma^\text{s}\) is the solute activity:

\[\mu^\text{s} = \mu^\text{s}_0 + k_\text{B}T\ln\left(\gamma^\text{s}c^\text{s}/c^\text{s}_0\right)\]

In the dilute limit, \(\gamma^\text{s}\to 1\), and thus

\[\underline D^\text{s} = k_\text{B}T\Omega(c^\text{s})^{-1}\underline L^\text{ss}\]

Conveniently, this cancels most of the “missing” prefactors we put in to compute the Onsager coefficient; hence,

• \(\underline D^\text{s}\) = Lss*(vacancy concentration)

where the concentration quantities are fractional. In the case of global thermal equilibrium, the vacancy concentration is the equilibrium concentration \(\exp(-(E^\text{v}_\text{form} - TS^\text{v}_\text{form})/k_\text{B}T)\).

A similar argument holds for the vacancy diffusivity in the dilute limit

• \(\underline D^\text{v}\) = L0vv + (solute concentration)*L1vv

The off-diagonal diffusivity terms are more complex as (1) they are non-symmetric (\(\underline D^\text{sv} \ne \underline D^\text{vs}\)), and (2) the vacancy-dependency of the solute activity and the solute-dependence of the vacancy activity needs to be known to properly include thermodynamic factors.