{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Input and output for Onsager transport calculation\n", "\n", "The Onsager calculators currently include two computational approaches to determining transport coefficients: an \"interstitial\" calculation, and a \"vacancy-mediated\" calculator. Below we describe the\n", "\n", "0. **Assumptions used in transport model** that are necessary to understand the data to be input, and the limitations of the results;\n", "1. **Crystal class setup** needed to initiate a calculation;\n", "2. **Interstitial calculator setup** needed for an single mobile species calculation, or,\n", "3. **Vacancy-mediated calculator setup** needed for a vacancy-mediated substitutional solute calculation;\n", "4. the creation of **VASP-style input files** to be run to generate input data;\n", "5. proper **Formatting of input data** to be compatible with the calculators; and\n", "6. **Interpretation of output** which includes how to convert output into transport coefficients.\n", "\n", "This follows the overall structure of a transport coefficient calculation. Broadly speaking, these are the steps necessary to compute transport coefficients:\n", "\n", "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.\n", "2. Generate lists of symmetry unrelated \"defect states\" and \"defect state transitions,\" along with the appropriate \"calculator object.\"\n", "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.\n", "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.\n", "5. Transform the output into physically relevant quantities (Onsager coefficients, solute diffusivities, mobilities, or drag ratios) with appropriate units." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Assumptions used in Onsager\n", "\n", "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:\n", "\n", "* **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.\n", "* **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.\n", "* **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.\n", "\n", "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.\n", "\n", "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.\n", "\n", "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).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Crystal class setup\n", "\n", "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.\n", "\n", "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.\n", "\n", "* 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.\n", "* 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.**\n", "* 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.\n", "* 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.\n", "\n", "Once initialized, two main internal operations take place:\n", "\n", "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`.\n", "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.\n", "\n", "**Note**: `Crystal`s 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.\n", "\n", "A few quick examples:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import sys\n", "sys.path.extend(['.', '..'])\n", "from onsager import crystal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, vacancy-diffusion\n", "\n", "Face-centered cubic crystals could be created either by entering the primitive basis:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n" ] } ], "source": [ "a0 = 1.\n", "FCCcrys = crystal.Crystal([a0*np.array([0,0.5,0.5]), \n", " a0*np.array([0.5,0,0.5]), \n", " a0*np.array([0.5,0.5,0])], \n", " [np.array([0.,0.,0.])], chemistry=['fcc'])\n", "print(FCCcrys)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or by entering the simple cubic unit cell with four atoms:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#Lattice:\n", " a1 = [ 0.5 0. 0.5]\n", " a2 = [ 0. 0.5 0.5]\n", " a3 = [-0.5 0. 0.5]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n" ] } ], "source": [ "FCCcrys2 = crystal.Crystal(a0*np.eye(3), \n", " [np.array([0.,0.,0.]), np.array([0,0.5,0.5]), \n", " np.array([0.5,0,0.5]), np.array([0.5,0.5,0])], \n", " chemistry=['fcc'])\n", "print(FCCcrys2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The effect of `noreduce` can be seen by regenerating the FCC crystal using the simple cubic unit cell:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#Lattice:\n", " a1 = [ 1. 0. 0.]\n", " a2 = [ 0. 1. 0.]\n", " a3 = [ 0. 0. 1.]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n", " (fcc) 0.1 = [ 0. 0.5 0.5]\n", " (fcc) 0.2 = [ 0.5 0. 0.5]\n", " (fcc) 0.3 = [ 0.5 0.5 0. ]\n" ] } ], "source": [ "FCCcrys3 = crystal.Crystal(a0*np.eye(3), \n", " [np.array([0.,0.,0.]), np.array([0,0.5,0.5]), \n", " np.array([0.5,0,0.5]), np.array([0.5,0.5,0])], \n", " chemistry=['fcc'], noreduce=True)\n", "print(FCCcrys3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rocksalt crystal, vacancy-diffusion\n", "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\":" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (Mg) 0.0 = [ 0. 0. 0.]\n", " (O) 1.0 = [ 0.5 0.5 0.5]\n" ] } ], "source": [ "MgO = crystal.Crystal([a0*np.array([0,0.5,0.5]), \n", " a0*np.array([0.5,0,0.5]), \n", " a0*np.array([0.5,0.5,0])], \n", " [[np.array([0.,0.,0.])], [np.array([0.5,0.5,0.5])]], \n", " chemistry=['Mg', 'O'])\n", "print(MgO)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, interstitial diffusion\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[array([ 0.5, 0.5, 0.5])]\n", "[array([ 0.75, 0.75, 0.75]), array([ 0.25, 0.25, 0.25])]\n", "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n", " (int) 1.0 = [ 0.5 0.5 0.5]\n", " (int) 1.1 = [ 0.75 0.75 0.75]\n", " (int) 1.2 = [ 0.25 0.25 0.25]\n" ] } ], "source": [ "octbasis = FCCcrys.Wyckoffpos(np.array([0.5, 0.5, 0.5]))\n", "tetbasis = FCCcrys.Wyckoffpos(np.array([0.25, 0.25, 0.25]))\n", "FCCcrysint = FCCcrys.addbasis(octbasis + tetbasis, ['int'])\n", "print(octbasis)\n", "print(tetbasis)\n", "print(FCCcrysint)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interstitial calculator setup\n", "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,\n", "\n", "* configurations are simply the Wyckoff positions of the particular sublattice (specified by a `chemistry` index);\n", "* transition states are pairs of configurations with a displacement vector that connects the initial to the final system.\n", "\n", "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.\n", "\n", "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.\n", "\n", "The algorithm in `jumpnetwork()` is rather simple: a transition is included if\n", "\n", "* the distance between the initial and final state is less than a cutoff distance, and\n", "* the line segment between the initial and final state does not come within a minimum distance of other defect states, and\n", "* the line segment between the initial and final state does not come within a minimum distance of *any* atomic site in the crystal.\n", "\n", "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.\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from onsager import OnsagerCalc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, vacancy-diffusion\n", "\n", "We identify the vacancy sites with the crystal sites in the lattice." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0]]\n" ] } ], "source": [ "chem = 0\n", "FCCsitelist = FCCcrys.sitelist(chem)\n", "print(FCCsitelist)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jump type 0\n", " 0 -> 0 dx= [ 0. -0.5 -0.5]\n", " 0 -> 0 dx= [-0. 0.5 0.5]\n", " 0 -> 0 dx= [ 0.5 0. 0.5]\n", " 0 -> 0 dx= [-0.5 -0. -0.5]\n", " 0 -> 0 dx= [ 0.5 0. -0.5]\n", " 0 -> 0 dx= [-0.5 -0. 0.5]\n", " 0 -> 0 dx= [ 0. -0.5 0.5]\n", " 0 -> 0 dx= [-0. 0.5 -0.5]\n", " 0 -> 0 dx= [-0.5 0.5 0. ]\n", " 0 -> 0 dx= [ 0.5 -0.5 -0. ]\n", " 0 -> 0 dx= [-0.5 -0.5 0. ]\n", " 0 -> 0 dx= [ 0.5 0.5 -0. ]\n" ] } ], "source": [ "chem = 0\n", "FCCjumpnetwork = FCCcrys.jumpnetwork(chem, cutoff=a0*0.78)\n", "for n, jn in enumerate(FCCjumpnetwork):\n", " print('Jump type {}'.format(n))\n", " for (i,j), dx in jn:\n", " print(' {} -> {} dx= {}'.format(i,j,dx))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Diffuser for atom 0 (fcc)\n", "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n", "states:\n", "i:+0.000,+0.000,+0.000\n", "transitions:\n", "i:+0.000,+0.000,+0.000^i:-1.000,+0.000,+0.000\n", "\n" ] } ], "source": [ "chem = 0\n", "FCCvacancydiffuser = OnsagerCalc.Interstitial(FCCcrys, chem, FCCsitelist, FCCjumpnetwork)\n", "print(FCCvacancydiffuser)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rocksalt crystal, vacancy-diffusion\n", "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\"." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0]]\n" ] } ], "source": [ "chem = 1\n", "MgOsitelist = MgO.sitelist(chem)\n", "print(MgOsitelist)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jump type 0\n", " 0 -> 0 dx= [ 0. 0.5 -0.5]\n", " 0 -> 0 dx= [-0. -0.5 0.5]\n", " 0 -> 0 dx= [-0.5 0. 0.5]\n", " 0 -> 0 dx= [ 0.5 -0. -0.5]\n", " 0 -> 0 dx= [ 0. -0.5 -0.5]\n", " 0 -> 0 dx= [-0. 0.5 0.5]\n", " 0 -> 0 dx= [-0.5 -0.5 0. ]\n", " 0 -> 0 dx= [ 0.5 0.5 -0. ]\n", " 0 -> 0 dx= [-0.5 0. -0.5]\n", " 0 -> 0 dx= [ 0.5 -0. 0.5]\n", " 0 -> 0 dx= [-0.5 0.5 0. ]\n", " 0 -> 0 dx= [ 0.5 -0.5 -0. ]\n" ] } ], "source": [ "chem = 1\n", "MgOjumpnetwork = MgO.jumpnetwork(chem, cutoff=a0*0.78)\n", "for n, jn in enumerate(MgOjumpnetwork):\n", " print('Jump type {}'.format(n))\n", " for (i,j), dx in jn:\n", " print(' {} -> {} dx= {}'.format(i,j,dx))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Diffuser for atom 1 (O)\n", "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (Mg) 0.0 = [ 0. 0. 0.]\n", " (O) 1.0 = [ 0.5 0.5 0.5]\n", "states:\n", "i:+0.500,+0.500,+0.500\n", "transitions:\n", "i:+0.500,+0.500,+0.500^i:+0.500,-0.500,+1.500\n", "\n" ] } ], "source": [ "chem = 1\n", "MgOdiffuser = OnsagerCalc.Interstitial(MgO, chem, MgOsitelist, MgOjumpnetwork)\n", "print(MgOdiffuser)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, interstitial diffusion\n", "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." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0], [1, 2]]\n" ] } ], "source": [ "chem = 1\n", "FCCintsitelist = FCCcrysint.sitelist(chem)\n", "print(FCCintsitelist)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jump type 0\n", " 0 -> 2 dx= [ 0.25 -0.25 0.25]\n", " 2 -> 0 dx= [-0.25 0.25 -0.25]\n", " 0 -> 1 dx= [-0.25 0.25 -0.25]\n", " 1 -> 0 dx= [ 0.25 -0.25 0.25]\n", " 0 -> 2 dx= [-0.25 0.25 0.25]\n", " 2 -> 0 dx= [ 0.25 -0.25 -0.25]\n", " 0 -> 2 dx= [-0.25 -0.25 -0.25]\n", " 2 -> 0 dx= [ 0.25 0.25 0.25]\n", " 0 -> 2 dx= [ 0.25 0.25 -0.25]\n", " 2 -> 0 dx= [-0.25 -0.25 0.25]\n", " 0 -> 1 dx= [ 0.25 -0.25 -0.25]\n", " 1 -> 0 dx= [-0.25 0.25 0.25]\n", " 0 -> 1 dx= [-0.25 -0.25 0.25]\n", " 1 -> 0 dx= [ 0.25 0.25 -0.25]\n", " 0 -> 1 dx= [ 0.25 0.25 0.25]\n", " 1 -> 0 dx= [-0.25 -0.25 -0.25]\n", "Jump type 1\n", " 2 -> 1 dx= [ 0. 0. 0.5]\n", " 1 -> 2 dx= [-0. -0. -0.5]\n", " 2 -> 1 dx= [ 0. -0.5 0. ]\n", " 1 -> 2 dx= [-0. 0.5 -0. ]\n", " 2 -> 1 dx= [-0.5 0. 0. ]\n", " 1 -> 2 dx= [ 0.5 -0. -0. ]\n", " 2 -> 1 dx= [ 0. 0.5 0. ]\n", " 1 -> 2 dx= [-0. -0.5 -0. ]\n", " 2 -> 1 dx= [ 0.5 0. 0. ]\n", " 1 -> 2 dx= [-0.5 -0. -0. ]\n", " 2 -> 1 dx= [ 0. 0. -0.5]\n", " 1 -> 2 dx= [-0. -0. 0.5]\n" ] } ], "source": [ "chem = 1\n", "FCCintjumpnetwork = FCCcrysint.jumpnetwork(chem, cutoff=a0*0.51)\n", "for n, jn in enumerate(FCCintjumpnetwork):\n", " print('Jump type {}'.format(n))\n", " for (i,j), dx in jn:\n", " print(' {} -> {} dx= {}'.format(i,j,dx))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Diffuser for atom 1 (int)\n", "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n", " (int) 1.0 = [ 0.5 0.5 0.5]\n", " (int) 1.1 = [ 0.75 0.75 0.75]\n", " (int) 1.2 = [ 0.25 0.25 0.25]\n", "states:\n", "i:+0.500,+0.500,+0.500\n", "i:+0.750,+0.750,+0.750\n", "transitions:\n", "i:+0.500,+0.500,+0.500^i:+0.250,+1.250,+0.250\n", "i:+0.250,+0.250,+0.250^i:+0.750,+0.750,-0.250\n", "\n" ] } ], "source": [ "chem = 1\n", "FCCintdiffuser = OnsagerCalc.Interstitial(FCCcrysint, chem, \n", " FCCintsitelist, FCCintjumpnetwork)\n", "print(FCCintdiffuser)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "basis:\n", "- - !numpy.ndarray [0.0, 0.0, 0.0]\n", "- - !numpy.ndarray [0.5, 0.5, 0.5]\n", " - !numpy.ndarray [0.75, 0.75, 0.75]\n", " - !numpy.ndarray [0.25, 0.25, 0.25]\n", "chemistry: [fcc, int]\n", "lattice: !numpy.ndarray\n", "- [0.0, 0.5, 0.5]\n", "- [0.5, 0.0, 0.5]\n", "- [0.5, 0.5, 0.0]\n", "lattice_constant: 1.0\n", "spins: null\n", "threshold: 1.0e-08\n", "chem: 1\n", "Dipole:\n", "- !numpy.ndarray\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", "- !numpy.ndarray\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", "Energy: [0, 0]\n", "Prefactor: [1, 1]\n", "sitelist:\n", "- [0]\n", "- [1, 2]\n", "DipoleT:\n", "- !numpy.ndarray\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", "- !numpy.ndarray\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", " - [0.0, 0.0, 0.0]\n", "EnergyT: [0, 0]\n", "PrefactorT: [1, 1]\n", "jumpnetwork:\n", "- - !!python/tuple\n", " - !!python/tuple [0, 2]\n", " - !numpy.ndarray [0.25, -0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [2, 0]\n", " - !numpy.ndarray [-0.25, 0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 1]\n", " - !numpy.ndarray [-0.25, 0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [1, 0]\n", " - !numpy.ndarray [0.25, -0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 2]\n", " - !numpy.ndarray [-0.25, 0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [2, 0]\n", " - !numpy.ndarray [0.25, -0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 2]\n", " - !numpy.ndarray [-0.25, -0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [2, 0]\n", " - !numpy.ndarray [0.25, 0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 2]\n", " - !numpy.ndarray [0.25, 0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [2, 0]\n", " - !numpy.ndarray [-0.25, -0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 1]\n", " - !numpy.ndarray [0.25, -0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [1, 0]\n", " - !numpy.ndarray [-0.25, 0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 1]\n", " - !numpy.ndarray [-0.25, -0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [1, 0]\n", " - !numpy.ndarray [0.25, 0.25, -0.25]\n", " - !!python/tuple\n", " - !!python/tuple [0, 1]\n", " - !numpy.ndarray [0.25, 0.25, 0.25]\n", " - !!python/tuple\n", " - !!python/tuple [1, 0]\n", " - !numpy.ndarray [-0.25, -0.25, -0.25]\n", "- - !!python/tuple\n", " - !!python/tuple [2, 1]\n", " - !numpy.ndarray [0.0, 0.0, 0.5]\n", " - !!python/tuple\n", " - !!python/tuple [1, 2]\n", " - !numpy.ndarray [-0.0, -0.0, -0.5]\n", " - !!python/tuple\n", " - !!python/tuple [2, 1]\n", " - !numpy.ndarray [0.0, -0.5, 0.0]\n", " - !!python/tuple\n", " - !!python/tuple [1, 2]\n", " - !numpy.ndarray [-0.0, 0.5, -0.0]\n", " - !!python/tuple\n", " - !!python/tuple [2, 1]\n", " - !numpy.ndarray [-0.5, 0.0, 0.0]\n", " - !!python/tuple\n", " - !!python/tuple [1, 2]\n", " - !numpy.ndarray [0.5, -0.0, -0.0]\n", " - !!python/tuple\n", " - !!python/tuple [2, 1]\n", " - !numpy.ndarray [0.0, 0.5, 0.0]\n", " - !!python/tuple\n", " - !!python/tuple [1, 2]\n", " - !numpy.ndarray [-0.0, -0.5, -0.0]\n", " - !!python/tuple\n", " - !!python/tuple [2, 1]\n", " - !numpy.ndarray [0.5, 0.0, 0.0]\n", " - !!python/tuple\n", " - !!python/tuple [1, 2]\n", " - !numpy.ndarray [-0.5, -0.0, -0.0]\n", " - !!python/tuple\n", " - !!python/tuple [2, 1]\n", " - !numpy.ndarray [0.0, 0.0, -0.5]\n", " - !!python/tuple\n", " - !!python/tuple [1, 2]\n", " - !numpy.ndarray [-0.0, -0.0, 0.5]\n", "\n" ] } ], "source": [ "print(FCCintdiffuser.crys.simpleYAML() + \n", " 'chem: {}\\n'.format(FCCintdiffuser.chem) + \n", " FCCintdiffuser.sitelistYAML(FCCintsitelist) + \n", " FCCintdiffuser.jumpnetworkYAML(FCCintjumpnetwork))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vacancy-mediated calculator setup\n", "For the vacancy mediated diffuser, the configurations and transition states are more complicated. First, we have three types of configurations:\n", "\n", "1. Vacancy, sufficiently far away from the solute to have zero interaction energy.\n", "2. Solute, sufficiently far away from the vacancy to have zero interaction energy.\n", "3. Vacancy-solute complexes.\n", "\n", "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.\n", "\n", "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.\n", "\n", "When we consider transition states, we have three types of transition states:\n", "\n", "1. Vacancy transitions, sufficiently far away from the solute to have zero interaction energy.\n", "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).\n", "3. Vacancy-solute complex transitions, where the vacancy and solute exchange place.\n", "\n", "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.\n", "\n", "To make a diffuser, we need to\n", "\n", "1. Identify the `sitelist` of the vacancies (and hence, solutes),\n", "2. Identify the `jumpnetwork` of the vacancies\n", "3. Determine the thermodynamic range\n", "\n", "then, the diffuser automatically constructs the complexes out to the thermodynamic range, and the full jumpnetworks.\n", "\n", "The vacancy-mediated diffuser also identifies unique **tags** for all configurations and transition states. The tags for *configurations* are strings with \n", "\n", "* `v:` followed by unit cell coordinates of site to three decimal digits for the vacancy;\n", "* `s:` followed by unit cell coordinates of site to three decimal digits for the solute;\n", "* `s:...-v:...` for a solute-vacancy complex.\n", "\n", "The *transition states* are strings with\n", "\n", "* `omega0:` + (initial vacancy configuration) + `^` + (final vacancy configuration);\n", "* `omega1:` + (initial solute-vacancy configuration) + `^` + (final vacancy configuration);\n", "* `omega2:` + (initial solute-vacancy configuration) + `^` + (final solute-vacancy configuration).\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, vacancy mediated-diffusion\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Diffuser for atom 0 (fcc), Nthermo=1\n", "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n", "vacancy configurations:\n", "v:+0.000,+0.000,+0.000\n", "solute configurations:\n", "s:+0.000,+0.000,+0.000\n", "solute-vacancy configurations:\n", "s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000\n", "omega0 jumps:\n", "omega0:v:+0.000,+0.000,+0.000^v:-1.000,+0.000,+0.000\n", "omega1 jumps:\n", "omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+0.000,-1.000,+0.000\n", "omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-2.000,+1.000,+0.000\n", "omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,-1.000^v:-1.000,+1.000,-1.000\n", "omega1:s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000^v:-2.000,+0.000,+0.000\n", "omega2 jumps:\n", "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\n", "\n" ] } ], "source": [ "chem = 0\n", "fivefreqdiffuser = OnsagerCalc.VacancyMediated(FCCcrys, chem, \n", " FCCsitelist, FCCjumpnetwork, 1)\n", "print(fivefreqdiffuser)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 `jumpnetwork`s, all of the precalculation and analysis needed for diffusion. This greatly speeds up the construction of the calculator." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Diffuser for atom 0 (fcc), Nthermo=1\n", "#Lattice:\n", " a1 = [ 0. 0.5 0.5]\n", " a2 = [ 0.5 0. 0.5]\n", " a3 = [ 0.5 0.5 0. ]\n", "#Basis:\n", " (fcc) 0.0 = [ 0. 0. 0.]\n", "vacancy configurations:\n", "v:+0.000,+0.000,+0.000\n", "solute configurations:\n", "s:+0.000,+0.000,+0.000\n", "solute-vacancy configurations:\n", "s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000\n", "omega0 jumps:\n", "omega0:v:+0.000,+0.000,+0.000^v:-1.000,+0.000,+0.000\n", "omega1 jumps:\n", "omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+0.000,-1.000,+0.000\n", "omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-2.000,+1.000,+0.000\n", "omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,-1.000^v:-1.000,+1.000,-1.000\n", "omega1:s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000^v:-2.000,+0.000,+0.000\n", "omega2 jumps:\n", "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\n", "\n" ] } ], "source": [ "import h5py\n", "# replace '/dev/null' with your file of choice, and remove backing_store=False\n", "# to read and write to an HDF5 file.\n", "f = h5py.File('/dev/null', 'w', driver='core', backing_store=False)\n", "fivefreqdiffuser.addhdf5(f) # adds the diffuser to the HDF5 file\n", "\n", "# how to read in (after opening `f` as an HDF5 file)\n", "fivefreqcopy = OnsagerCalc.VacancyMediated.loadhdf5(f) # creates a new diffuser from HDF5\n", "f.close() # close up the HDF5 file\n", "print(fivefreqcopy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## VASP-style input files\n", "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.\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from onsager import automator\n", "import tarfile" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, interstitial diffusion\n", "We will need to construct (and relax) appropriate intersitial sites, and the transition states between them." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on method makesupercells in module onsager.OnsagerCalc:\n", "\n", "makesupercells(super_n) method of onsager.OnsagerCalc.Interstitial instance\n", " Take in a supercell matrix, then generate all of the supercells needed to compute\n", " site energies and transitions (corresponding to the representatives).\n", " \n", " :param super_n: 3x3 integer matrix to define our supercell\n", " :return superdict: dictionary of ``states``, ``transitions``, ``transmapping``,\n", " and ``indices`` that correspond to dictionaries with tags.\n", " \n", " * superdict['states'][i] = supercell of site;\n", " * superdict['transitions'][n] = (supercell initial, supercell final);\n", " * superdict['transmapping'][n] = ((site tag, groupop, mapping), (site tag, groupop, mapping))\n", " * superdict['indices'][tag] = index of tag, where tag is either a state or transition tag.\n", "\n" ] } ], "source": [ "help(FCCintdiffuser.makesupercells)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 2. 0. 0.]\n", " [ 0. 2. 0.]\n", " [ 0. 0. 2.]]\n" ] } ], "source": [ "N = np.array([[-2,2,2],[2,-2,2],[2,2,-2]]) # 32 atom FCC supercell\n", "print(np.dot(FCCcrys.lattice, N))\n", "FCCintsupercells = FCCintdiffuser.makesupercells(N)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function supercelltar in module onsager.automator:\n", "\n", "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')\n", " Takes in a tarfile (needs to be open for writing) and a supercelldict (from a\n", " diffuser) and creates the full directory structure inside the tarfile. Best used in\n", " a form like\n", " \n", " ::\n", " \n", " with tarfile.open('supercells.tar.gz', mode='w:gz') as tar:\n", " automator.supercelltar(tar, supercelldict)\n", " \n", " :param tar: tarfile open for writing; may contain other files in advance.\n", " :param superdict: dictionary of ``states``, ``transitions``, ``transmapping``, ``indices`` that\n", " correspond to dictionaries with tags; the final tag ``reference`` is the basesupercell\n", " for calculations without defects.\n", " \n", " * superdict['states'][i] = supercell of state;\n", " * superdict['transitions'][n] = (supercell initial, supercell final);\n", " * superdict['transmapping'][n] = ((site tag, groupop, mapping), (site tag, groupop, mapping))\n", " * superdict['indices'][tag] = (type, index) of tag, where tag is either a state or transition tag; or...\n", " * superdict['indices'][tag] = index of tag, where tag is either a state or transition tag.\n", " * superdict['reference'] = (optional) supercell reference, no defects\n", " \n", " :param filemode: mode to use for files (default: 664)\n", " :param directmode: mode to use for directories (default: 775)\n", " :param timestamp: UNIX time for files; if None, use current time (default)\n", " :param INCARrelax: contents of INCAR file to use for relaxation; must contain {system} to be replaced\n", " by tag value (default: automator.INCARrelax)\n", " :param INCARNEB: contents of INCAR file to use for NEB; must contain {system} to be replaced\n", " by tag value (default: automator.INCARNEB)\n", " :param KPOINTS: contents of KPOINTS file (default: gamma-point only calculation);\n", " if None or empty, no KPOINTS file at all\n", " :param basedir: prepended to all files/directories (default: '')\n", " :param statename: prepended to all state names, before 2 digit number (default: relax.)\n", " :param transitionname: prepended to all transition names, before 2 digit number (default: neb.)\n", " :param IDformat: format for integer tags (default: {:02d})\n", " :param JSONdict: name of JSON file storing the tags corresponding to each directory (default: tags.json)\n", " :param YAMLdef: YAML file containing full definition of supercells, relationship, etc. (default: supercell.yaml);\n", " set to None to not output. **may want to change this to None for the future**\n", "\n" ] } ], "source": [ "help(automator.supercelltar)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true }, "outputs": [], "source": [ "with tarfile.open('io-test-int.tar.gz', mode='w:gz') as tar:\n", " automator.supercelltar(tar, FCCintsupercells)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [], "source": [ "tar = tarfile.open('io-test-int.tar.gz', mode='r:gz')" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "?rw-rw-r-- 0/0 244 2017-07-11 16:14:30 INCAR.relax \n", "?rw-rw-r-- 0/0 305 2017-07-11 16:14:30 INCAR.NEB \n", "?rw-rw-r-- 0/0 31 2017-07-11 16:14:30 KPOINTS \n", "?rwxrwxr-x 0/0 1344 2017-07-11 16:14:30 trans.pl \n", "?rwxrwxr-x 0/0 6283 2017-07-11 16:14:30 nebmake.pl \n", "?rw-rw-r-- 0/0 25975 2017-07-11 16:14:30 Vasp.pm \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:30 relax.01/ \n", "?rw-rw-r-- 0/0 2267 2017-07-11 16:14:30 relax.01/POSCAR \n", "?rw-rw-r-- 0/0 258 2017-07-11 16:14:30 relax.01/INCAR \n", "?rw-rw-r-- 0/0 35 2017-07-11 16:14:30 relax.01/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 relax.01/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 relax.01/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:30 relax.00/ \n", "?rw-rw-r-- 0/0 2267 2017-07-11 16:14:30 relax.00/POSCAR \n", "?rw-rw-r-- 0/0 258 2017-07-11 16:14:30 relax.00/INCAR \n", "?rw-rw-r-- 0/0 35 2017-07-11 16:14:30 relax.00/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 relax.00/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 relax.00/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:30 neb.01/ \n", "?rw-rw-r-- 0/0 2298 2017-07-11 16:14:30 neb.01/POS.init \n", "?rw-rw-r-- 0/0 2296 2017-07-11 16:14:30 neb.01/POS.final \n", "?rw-rw-r-- 0/0 342 2017-07-11 16:14:30 neb.01/INCAR \n", "?rw-rw-r-- 0/0 58 2017-07-11 16:14:30 neb.01/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 neb.01/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 neb.01/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:30 neb.00/ \n", "?rw-rw-r-- 0/0 2298 2017-07-11 16:14:30 neb.00/POS.init \n", "?rw-rw-r-- 0/0 2296 2017-07-11 16:14:30 neb.00/POS.final \n", "?rw-rw-r-- 0/0 342 2017-07-11 16:14:30 neb.00/INCAR \n", "?rw-rw-r-- 0/0 58 2017-07-11 16:14:30 neb.00/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 neb.00/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:30 neb.00/POTCAR -> ../POTCAR \n", "?rw-rw-r-- 0/0 191 2017-07-11 16:14:30 neb.00/trans.init \n", "?rw-rw-r-- 0/0 191 2017-07-11 16:14:30 neb.00/trans.final \n", "?rw-rw-r-- 0/0 191 2017-07-11 16:14:30 neb.01/trans.init \n", "?rw-rw-r-- 0/0 191 2017-07-11 16:14:30 neb.01/trans.final \n", "?rw-rw-r-- 0/0 1170 2017-07-11 16:14:30 Makefile \n", "?rw-rw-r-- 0/0 14 2017-07-11 16:14:30 relax.00/NEBlist \n", "?rw-rw-r-- 0/0 7 2017-07-11 16:14:30 relax.01/NEBlist \n", "?rw-rw-r-- 0/0 213 2017-07-11 16:14:30 tags.json \n", "?rw-rw-r-- 0/0 1243794 2017-07-11 16:14:30 supercell.yaml \n" ] } ], "source": [ "tar.list()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents of the `Makefile`:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# Makefile to construct NEB input from relaxation output\n", "# we set this so that the makefile doesn't use builtin implicit rules\n", "MAKEFLAGS = -rk\n", "\n", "makeneb := \"./nebmake.pl\"\n", "transform := \"./trans.pl\"\n", "\n", "Nimages ?= 1\n", "\n", ".PHONY: help\n", "\n", "target := $(foreach neb, $(wildcard neb.*), $(neb)/01/POSCAR)\n", "target: $(target)\n", "\n", "help:\n", "\t@echo \"# Creates input POSCAR for NEB runs, once relaxation runs are complete\"\n", "\t@echo \"# Uses CONTCAR in relaxation directories to create initial run geometry\"\n", "\t@echo \"# environment variable: Nimages (default: $(Nimages))\"\n", "\t@echo \"# target files:\"\n", "\t@echo $(target) | sed 's/ / /g'\n", "\t@echo \"# default target: all\"\n", "\n", "neb.%: neb.%/01/POSCAR neb.%/POSCAR.init neb.%/POSCAR.final\n", "\n", "neb.%/01/POSCAR: neb.%/POSCAR.init neb.%/POSCAR.final\n", "\t@$(makeneb) $^ $(Nimages)\n", "\n", "neb.%/POSCAR.init:\n", "\t@$(transform) $^ > $@\n", "\n", "neb.%/POSCAR.final:\n", "\t@$(transform) $^ > $@\n", "\n", "###############################################################\n", "# structure of NEB runs:\n", "neb.00/POSCAR.init: neb.00/trans.init relax.00/CONTCAR\n", "neb.00/POSCAR.final: neb.00/trans.final relax.00/CONTCAR\n", "neb.01/POSCAR.init: neb.01/trans.init relax.01/CONTCAR\n", "neb.01/POSCAR.final: neb.01/trans.final relax.00/CONTCAR\n", "\n" ] } ], "source": [ "with tar.extractfile('Makefile') as f:\n", " print(f.read().decode('ascii'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents of the `tags.json` file:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " \"neb.00\": \"i:+0.250,+0.250,+0.250^i:+0.750,+0.750,-0.250\",\n", " \"neb.01\": \"i:+0.500,+0.500,+0.500^i:+0.250,+1.250,+0.250\",\n", " \"relax.00\": \"i:+0.750,+0.750,+0.750\",\n", " \"relax.01\": \"i:+0.500,+0.500,+0.500\"\n", "}\n", "\n" ] } ], "source": [ "with tar.extractfile('tags.json') as f:\n", " print(f.read().decode('ascii'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents of one `POSCAR` file for relaxation of a configuration:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i:+0.750,+0.750,+0.750 fcc(32),int_i(1)\n", "1.0\n", " 2.0000000000000000 0.0000000000000000 0.0000000000000000\n", " 0.0000000000000000 2.0000000000000000 0.0000000000000000\n", " 0.0000000000000000 0.0000000000000000 2.0000000000000000\n", "32 1\n", "Direct\n", " 0.0000000000000000 0.0000000000000000 0.0000000000000000\n", " 0.2500000000000000 0.2500000000000000 0.0000000000000000\n", " 0.5000000000000000 0.5000000000000000 0.0000000000000000\n", " 0.7500000000000000 0.7500000000000000 0.0000000000000000\n", " 0.2500000000000000 0.0000000000000000 0.2500000000000000\n", " 0.5000000000000000 0.2500000000000000 0.2500000000000000\n", " 0.7500000000000000 0.5000000000000000 0.2500000000000000\n", " 0.0000000000000000 0.7500000000000000 0.2500000000000000\n", " 0.5000000000000000 0.0000000000000000 0.5000000000000000\n", " 0.7500000000000000 0.2500000000000000 0.5000000000000000\n", " 0.0000000000000000 0.5000000000000000 0.5000000000000000\n", " 0.2500000000000000 0.7500000000000000 0.5000000000000000\n", " 0.7500000000000000 0.0000000000000000 0.7500000000000000\n", " 0.0000000000000000 0.2500000000000000 0.7500000000000000\n", " 0.2500000000000000 0.5000000000000000 0.7500000000000000\n", " 0.5000000000000000 0.7500000000000000 0.7500000000000000\n", " 0.0000000000000000 0.2500000000000000 0.2500000000000000\n", " 0.2500000000000000 0.5000000000000000 0.2500000000000000\n", " 0.5000000000000000 0.7500000000000000 0.2500000000000000\n", " 0.7500000000000000 0.0000000000000000 0.2500000000000000\n", " 0.2500000000000000 0.2500000000000000 0.5000000000000000\n", " 0.5000000000000000 0.5000000000000000 0.5000000000000000\n", " 0.7500000000000000 0.7500000000000000 0.5000000000000000\n", " 0.0000000000000000 0.0000000000000000 0.5000000000000000\n", " 0.5000000000000000 0.2500000000000000 0.7500000000000000\n", " 0.7500000000000000 0.5000000000000000 0.7500000000000000\n", " 0.0000000000000000 0.7500000000000000 0.7500000000000000\n", " 0.2500000000000000 0.0000000000000000 0.7500000000000000\n", " 0.7500000000000000 0.2500000000000000 0.0000000000000000\n", " 0.0000000000000000 0.5000000000000000 0.0000000000000000\n", " 0.2500000000000000 0.7500000000000000 0.0000000000000000\n", " 0.5000000000000000 0.0000000000000000 0.0000000000000000\n", " 0.3750000000000000 0.3750000000000000 0.3750000000000000\n", "\n" ] } ], "source": [ "with tar.extractfile('relax.00/POSCAR') as f:\n", " print(f.read().decode('ascii'))" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": true }, "outputs": [], "source": [ "tar.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, vacancy mediated-diffusion\n", "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." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on method makesupercells in module onsager.OnsagerCalc:\n", "\n", "makesupercells(super_n) method of onsager.OnsagerCalc.VacancyMediated instance\n", " Take in a supercell matrix, then generate all of the supercells needed to compute\n", " site energies and transitions (corresponding to the representatives).\n", " \n", " Note: the states are lone vacancy, lone solute, solute-vacancy complexes in\n", " our thermodynamic range. Note that there will be escape states are endpoints of\n", " some omega1 jumps. They are not relaxed, and have no pre-existing tag. They will\n", " only be output as a single endpoint of an NEB run; there may be symmetry equivalent\n", " duplicates, as we construct these supercells on an as needed basis.\n", " \n", " We've got a few classes of warnings (from most egregious to least) that can issued\n", " if the supercell is too small; the analysis will continue despite any warnings:\n", " \n", " 1. Thermodynamic shell states map to different states in supercell\n", " 2. Thermodynamic shell states are not unique in supercell (multiplicity)\n", " 3. Kinetic shell states map to different states in supercell\n", " 4. Kinetic shell states are not unique in supercell (multiplicity)\n", " \n", " The lowest level can still be run reliably but runs the risk of errors in escape transition\n", " barriers. Extreme caution should be used if any of the other warnings are raised.\n", " \n", " :param super_n: 3x3 integer matrix to define our supercell\n", " :return superdict: dictionary of ``states``, ``transitions``, ``transmapping``,\n", " ``indices`` that correspond to dictionaries with tags; the final tag\n", " ``reference`` is the basesupercell for calculations without defects.\n", " \n", " * superdict['states'][i] = supercell of state;\n", " * superdict['transitions'][n] = (supercell initial, supercell final);\n", " * superdict['transmapping'][n] = ((site tag, groupop, mapping), (site tag, groupop, mapping))\n", " * superdict['indices'][tag] = (type, index) of tag, where tag is either a state or transition tag.\n", " * superdict['reference'] = supercell reference, no defects\n", "\n" ] } ], "source": [ "help(fivefreqdiffuser.makesupercells)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 3. 0. 0.]\n", " [ 0. 3. 0.]\n", " [ 0. 0. 3.]]\n" ] } ], "source": [ "N = np.array([[-3,3,3],[3,-3,3],[3,3,-3]]) # 108 atom FCC supercell\n", "print(np.dot(FCCcrys.lattice, N))\n", "fivefreqsupercells = fivefreqdiffuser.makesupercells(N)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true }, "outputs": [], "source": [ "with tarfile.open('io-test-fivefreq.tar.gz', mode='w:gz') as tar:\n", " automator.supercelltar(tar, fivefreqsupercells)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "tar = tarfile.open('io-test-fivefreq.tar.gz', mode='r:gz')" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "?rw-rw-r-- 0/0 244 2017-07-11 16:14:53 INCAR.relax \n", "?rw-rw-r-- 0/0 305 2017-07-11 16:14:53 INCAR.NEB \n", "?rw-rw-r-- 0/0 31 2017-07-11 16:14:53 KPOINTS \n", "?rwxrwxr-x 0/0 1344 2017-07-11 16:14:53 trans.pl \n", "?rwxrwxr-x 0/0 6283 2017-07-11 16:14:53 nebmake.pl \n", "?rw-rw-r-- 0/0 25975 2017-07-11 16:14:53 Vasp.pm \n", "?rw-rw-r-- 0/0 6844 2017-07-11 16:14:53 POSCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 relax.00/ \n", "?rw-rw-r-- 0/0 6784 2017-07-11 16:14:53 relax.00/POSCAR \n", "?rw-rw-r-- 0/0 258 2017-07-11 16:14:53 relax.00/INCAR \n", "?rw-rw-r-- 0/0 35 2017-07-11 16:14:53 relax.00/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 relax.00/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 relax.00/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 relax.02/ \n", "?rw-rw-r-- 0/0 6845 2017-07-11 16:14:53 relax.02/POSCAR \n", "?rw-rw-r-- 0/0 258 2017-07-11 16:14:53 relax.02/INCAR \n", "?rw-rw-r-- 0/0 35 2017-07-11 16:14:53 relax.02/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 relax.02/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 relax.02/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 relax.01/ \n", "?rw-rw-r-- 0/0 6807 2017-07-11 16:14:53 relax.01/POSCAR \n", "?rw-rw-r-- 0/0 281 2017-07-11 16:14:53 relax.01/INCAR \n", "?rw-rw-r-- 0/0 58 2017-07-11 16:14:53 relax.01/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 relax.01/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 relax.01/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 neb.00/ \n", "?rw-rw-r-- 0/0 6822 2017-07-11 16:14:53 neb.00/POS.init \n", "?rw-rw-r-- 0/0 6820 2017-07-11 16:14:53 neb.00/POS.final \n", "?rw-rw-r-- 0/0 349 2017-07-11 16:14:53 neb.00/INCAR \n", "?rw-rw-r-- 0/0 65 2017-07-11 16:14:53 neb.00/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.00/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.00/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 neb.02/ \n", "?rw-rw-r-- 0/0 6845 2017-07-11 16:14:53 neb.02/POS.init \n", "?rw-rw-r-- 0/0 6843 2017-07-11 16:14:53 neb.02/POS.final \n", "?rw-rw-r-- 0/0 372 2017-07-11 16:14:53 neb.02/INCAR \n", "?rw-rw-r-- 0/0 88 2017-07-11 16:14:53 neb.02/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.02/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.02/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 neb.04/ \n", "?rw-rw-r-- 0/0 6845 2017-07-11 16:14:53 neb.04/POS.init \n", "?rw-rw-r-- 0/0 6843 2017-07-11 16:14:53 neb.04/POSCAR.final \n", "?rw-rw-r-- 0/0 372 2017-07-11 16:14:53 neb.04/INCAR \n", "?rw-rw-r-- 0/0 88 2017-07-11 16:14:53 neb.04/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.04/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.04/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 neb.01/ \n", "?rw-rw-r-- 0/0 6845 2017-07-11 16:14:53 neb.01/POS.init \n", "?rw-rw-r-- 0/0 6843 2017-07-11 16:14:53 neb.01/POSCAR.final \n", "?rw-rw-r-- 0/0 372 2017-07-11 16:14:53 neb.01/INCAR \n", "?rw-rw-r-- 0/0 88 2017-07-11 16:14:53 neb.01/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.01/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.01/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 neb.03/ \n", "?rw-rw-r-- 0/0 6845 2017-07-11 16:14:53 neb.03/POS.init \n", "?rw-rw-r-- 0/0 6843 2017-07-11 16:14:53 neb.03/POSCAR.final \n", "?rw-rw-r-- 0/0 372 2017-07-11 16:14:53 neb.03/INCAR \n", "?rw-rw-r-- 0/0 88 2017-07-11 16:14:53 neb.03/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.03/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.03/POTCAR -> ../POTCAR \n", "?rwxrwxr-x 0/0 0 2017-07-11 16:14:53 neb.05/ \n", "?rw-rw-r-- 0/0 6868 2017-07-11 16:14:53 neb.05/POS.init \n", "?rw-rw-r-- 0/0 6866 2017-07-11 16:14:53 neb.05/POS.final \n", "?rw-rw-r-- 0/0 395 2017-07-11 16:14:53 neb.05/INCAR \n", "?rw-rw-r-- 0/0 111 2017-07-11 16:14:53 neb.05/incar.sed \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.05/KPOINTS -> ../KPOINTS \n", "?rw-rw-r-- 0/0 0 2017-07-11 16:14:53 neb.05/POTCAR -> ../POTCAR \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.00/trans.init \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.00/trans.final \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.01/trans.init \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.02/trans.init \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.02/trans.final \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.03/trans.init \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.04/trans.init \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.05/trans.init \n", "?rw-rw-r-- 0/0 420 2017-07-11 16:14:53 neb.05/trans.final \n", "?rw-rw-r-- 0/0 1447 2017-07-11 16:14:53 Makefile \n", "?rw-rw-r-- 0/0 7 2017-07-11 16:14:53 relax.00/NEBlist \n", "?rw-rw-r-- 0/0 35 2017-07-11 16:14:53 relax.01/NEBlist \n", "?rw-rw-r-- 0/0 710 2017-07-11 16:14:53 tags.json \n", "?rw-rw-r-- 0/0 3458541 2017-07-11 16:14:53 supercell.yaml \n" ] } ], "source": [ "tar.list()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents of `Makefile`:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# Makefile to construct NEB input from relaxation output\n", "# we set this so that the makefile doesn't use builtin implicit rules\n", "MAKEFLAGS = -rk\n", "\n", "makeneb := \"./nebmake.pl\"\n", "transform := \"./trans.pl\"\n", "\n", "Nimages ?= 1\n", "\n", ".PHONY: help\n", "\n", "target := $(foreach neb, $(wildcard neb.*), $(neb)/01/POSCAR)\n", "target: $(target)\n", "\n", "help:\n", "\t@echo \"# Creates input POSCAR for NEB runs, once relaxation runs are complete\"\n", "\t@echo \"# Uses CONTCAR in relaxation directories to create initial run geometry\"\n", "\t@echo \"# environment variable: Nimages (default: $(Nimages))\"\n", "\t@echo \"# target files:\"\n", "\t@echo $(target) | sed 's/ / /g'\n", "\t@echo \"# default target: all\"\n", "\n", "neb.%: neb.%/01/POSCAR neb.%/POSCAR.init neb.%/POSCAR.final\n", "\n", "neb.%/01/POSCAR: neb.%/POSCAR.init neb.%/POSCAR.final\n", "\t@$(makeneb) $^ $(Nimages)\n", "\n", "neb.%/POSCAR.init:\n", "\t@$(transform) $^ > $@\n", "\n", "neb.%/POSCAR.final:\n", "\t@$(transform) $^ > $@\n", "\n", "###############################################################\n", "# structure of NEB runs:\n", "neb.00/POSCAR.init: neb.00/trans.init relax.00/CONTCAR\n", "neb.00/POSCAR.final: neb.00/trans.final relax.00/CONTCAR\n", "neb.01/POSCAR.init: neb.01/trans.init relax.01/CONTCAR\n", "neb.02/POSCAR.init: neb.02/trans.init relax.01/CONTCAR\n", "neb.02/POSCAR.final: neb.02/trans.final relax.01/CONTCAR\n", "neb.03/POSCAR.init: neb.03/trans.init relax.01/CONTCAR\n", "neb.04/POSCAR.init: neb.04/trans.init relax.01/CONTCAR\n", "neb.05/POSCAR.init: neb.05/trans.init relax.01/CONTCAR\n", "neb.05/POSCAR.final: neb.05/trans.final relax.01/CONTCAR\n", "\n" ] } ], "source": [ "with tar.extractfile('Makefile') as f:\n", " print(f.read().decode('ascii'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents of the `tags.json` file:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " \"neb.00\": \"omega0:v:+0.000,+0.000,+0.000^v:-1.000,+0.000,+0.000\",\n", " \"neb.01\": \"omega1:s:+0.000,+0.000,+0.000-v:+0.000,+1.000,-1.000^v:-1.000,+1.000,-1.000\",\n", " \"neb.02\": \"omega1:s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000^v:+0.000,-1.000,+0.000\",\n", " \"neb.03\": \"omega1:s:+0.000,+0.000,+0.000-v:-1.000,+0.000,+0.000^v:-2.000,+0.000,+0.000\",\n", " \"neb.04\": \"omega1:s:+0.000,+0.000,+0.000-v:-1.000,+1.000,+0.000^v:-2.000,+1.000,+0.000\",\n", " \"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\",\n", " \"relax.00\": \"v:+0.000,+0.000,+0.000\",\n", " \"relax.01\": \"s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000\",\n", " \"relax.02\": \"s:+0.000,+0.000,+0.000\"\n", "}\n", "\n" ] } ], "source": [ "with tar.extractfile('tags.json') as f:\n", " print(f.read().decode('ascii'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents of one `POSCAR` file for relaxation of a configuration:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "s:+0.000,+0.000,+0.000-v:+1.000,-1.000,+0.000 fcc(106),solute(1)\n", "1.0\n", " 3.0000000000000000 0.0000000000000000 0.0000000000000000\n", " 0.0000000000000000 3.0000000000000000 0.0000000000000000\n", " 0.0000000000000000 0.0000000000000000 3.0000000000000000\n", "106 1\n", "Direct\n", " 0.1666666666666667 0.1666666666666667 0.0000000000000000\n", " 0.3333333333333333 0.3333333333333333 0.0000000000000000\n", " 0.5000000000000000 0.5000000000000000 0.0000000000000000\n", " 0.6666666666666666 0.6666666666666666 0.0000000000000000\n", " 0.8333333333333333 0.8333333333333333 0.0000000000000000\n", " 0.1666666666666667 0.0000000000000000 0.1666666666666667\n", " 0.3333333333333333 0.1666666666666667 0.1666666666666667\n", " 0.5000000000000000 0.3333333333333333 0.1666666666666667\n", " 0.6666666666666666 0.5000000000000000 0.1666666666666667\n", " 0.8333333333333333 0.6666666666666666 0.1666666666666667\n", " 0.0000000000000000 0.8333333333333333 0.1666666666666667\n", " 0.3333333333333333 0.0000000000000000 0.3333333333333333\n", " 0.5000000000000000 0.1666666666666667 0.3333333333333333\n", " 0.6666666666666666 0.3333333333333333 0.3333333333333333\n", " 0.8333333333333333 0.5000000000000000 0.3333333333333333\n", " 0.0000000000000000 0.6666666666666666 0.3333333333333333\n", " 0.1666666666666667 0.8333333333333333 0.3333333333333333\n", " 0.5000000000000000 0.0000000000000000 0.5000000000000000\n", " 0.6666666666666666 0.1666666666666667 0.5000000000000000\n", " 0.8333333333333333 0.3333333333333333 0.5000000000000000\n", " 0.0000000000000000 0.5000000000000000 0.5000000000000000\n", " 0.1666666666666667 0.6666666666666666 0.5000000000000000\n", " 0.3333333333333333 0.8333333333333333 0.5000000000000000\n", " 0.6666666666666666 0.0000000000000000 0.6666666666666666\n", " 0.8333333333333333 0.1666666666666667 0.6666666666666666\n", " 0.0000000000000000 0.3333333333333333 0.6666666666666666\n", " 0.1666666666666667 0.5000000000000000 0.6666666666666666\n", " 0.3333333333333333 0.6666666666666666 0.6666666666666666\n", " 0.5000000000000000 0.8333333333333333 0.6666666666666666\n", " 0.8333333333333333 0.0000000000000000 0.8333333333333333\n", " 0.0000000000000000 0.1666666666666667 0.8333333333333333\n", " 0.1666666666666667 0.3333333333333333 0.8333333333333333\n", " 0.3333333333333333 0.5000000000000000 0.8333333333333333\n", " 0.5000000000000000 0.6666666666666666 0.8333333333333333\n", " 0.6666666666666666 0.8333333333333333 0.8333333333333333\n", " 0.0000000000000000 0.1666666666666667 0.1666666666666667\n", " 0.1666666666666667 0.3333333333333333 0.1666666666666667\n", " 0.3333333333333333 0.5000000000000000 0.1666666666666667\n", " 0.5000000000000000 0.6666666666666666 0.1666666666666667\n", " 0.6666666666666666 0.8333333333333333 0.1666666666666667\n", " 0.8333333333333333 0.0000000000000000 0.1666666666666667\n", " 0.1666666666666667 0.1666666666666667 0.3333333333333333\n", " 0.3333333333333333 0.3333333333333333 0.3333333333333333\n", " 0.5000000000000000 0.5000000000000000 0.3333333333333333\n", " 0.6666666666666666 0.6666666666666666 0.3333333333333333\n", " 0.8333333333333333 0.8333333333333333 0.3333333333333333\n", " 0.0000000000000000 0.0000000000000000 0.3333333333333333\n", " 0.3333333333333333 0.1666666666666667 0.5000000000000000\n", " 0.5000000000000000 0.3333333333333333 0.5000000000000000\n", " 0.6666666666666666 0.5000000000000000 0.5000000000000000\n", " 0.8333333333333333 0.6666666666666666 0.5000000000000000\n", " 0.0000000000000000 0.8333333333333333 0.5000000000000000\n", " 0.1666666666666667 0.0000000000000000 0.5000000000000000\n", " 0.5000000000000000 0.1666666666666667 0.6666666666666666\n", " 0.6666666666666666 0.3333333333333333 0.6666666666666666\n", " 0.8333333333333333 0.5000000000000000 0.6666666666666666\n", " 0.0000000000000000 0.6666666666666666 0.6666666666666666\n", " 0.1666666666666667 0.8333333333333333 0.6666666666666666\n", " 0.3333333333333333 0.0000000000000000 0.6666666666666666\n", " 0.6666666666666666 0.1666666666666667 0.8333333333333333\n", " 0.8333333333333333 0.3333333333333333 0.8333333333333333\n", " 0.0000000000000000 0.5000000000000000 0.8333333333333333\n", " 0.1666666666666667 0.6666666666666666 0.8333333333333333\n", " 0.3333333333333333 0.8333333333333333 0.8333333333333333\n", " 0.5000000000000000 0.0000000000000000 0.8333333333333333\n", " 0.0000000000000000 0.3333333333333333 0.0000000000000000\n", " 0.1666666666666667 0.5000000000000000 0.0000000000000000\n", " 0.3333333333333333 0.6666666666666666 0.0000000000000000\n", " 0.5000000000000000 0.8333333333333333 0.0000000000000000\n", " 0.6666666666666666 0.0000000000000000 0.0000000000000000\n", " 0.0000000000000000 0.3333333333333333 0.3333333333333333\n", " 0.1666666666666667 0.5000000000000000 0.3333333333333333\n", " 0.3333333333333333 0.6666666666666666 0.3333333333333333\n", " 0.5000000000000000 0.8333333333333333 0.3333333333333333\n", " 0.6666666666666666 0.0000000000000000 0.3333333333333333\n", " 0.8333333333333333 0.1666666666666667 0.3333333333333333\n", " 0.1666666666666667 0.3333333333333333 0.5000000000000000\n", " 0.3333333333333333 0.5000000000000000 0.5000000000000000\n", " 0.5000000000000000 0.6666666666666666 0.5000000000000000\n", " 0.6666666666666666 0.8333333333333333 0.5000000000000000\n", " 0.8333333333333333 0.0000000000000000 0.5000000000000000\n", " 0.0000000000000000 0.1666666666666667 0.5000000000000000\n", " 0.3333333333333333 0.3333333333333333 0.6666666666666666\n", " 0.5000000000000000 0.5000000000000000 0.6666666666666666\n", " 0.6666666666666666 0.6666666666666666 0.6666666666666666\n", " 0.8333333333333333 0.8333333333333333 0.6666666666666666\n", " 0.0000000000000000 0.0000000000000000 0.6666666666666666\n", " 0.1666666666666667 0.1666666666666667 0.6666666666666666\n", " 0.5000000000000000 0.3333333333333333 0.8333333333333333\n", " 0.6666666666666666 0.5000000000000000 0.8333333333333333\n", " 0.8333333333333333 0.6666666666666666 0.8333333333333333\n", " 0.0000000000000000 0.8333333333333333 0.8333333333333333\n", " 0.1666666666666667 0.0000000000000000 0.8333333333333333\n", " 0.3333333333333333 0.1666666666666667 0.8333333333333333\n", " 0.6666666666666666 0.3333333333333333 0.0000000000000000\n", " 0.8333333333333333 0.5000000000000000 0.0000000000000000\n", " 0.0000000000000000 0.6666666666666666 0.0000000000000000\n", " 0.1666666666666667 0.8333333333333333 0.0000000000000000\n", " 0.3333333333333333 0.0000000000000000 0.0000000000000000\n", " 0.5000000000000000 0.1666666666666667 0.0000000000000000\n", " 0.8333333333333333 0.3333333333333333 0.1666666666666667\n", " 0.0000000000000000 0.5000000000000000 0.1666666666666667\n", " 0.1666666666666667 0.6666666666666666 0.1666666666666667\n", " 0.3333333333333333 0.8333333333333333 0.1666666666666667\n", " 0.5000000000000000 0.0000000000000000 0.1666666666666667\n", " 0.6666666666666666 0.1666666666666667 0.1666666666666667\n", " 0.0000000000000000 0.0000000000000000 0.0000000000000000\n", "\n" ] } ], "source": [ "with tar.extractfile('relax.01/POSCAR') as f:\n", " print(f.read().decode('ascii'))" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": true }, "outputs": [], "source": [ "tar.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Formatting of input data\n", "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.\n", "\n", "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\n", "$$\\rho = Z^{-1}\\rho^0 \\exp(-E/k_\\text{B}T)$$\n\n", "\n", "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\n", "$$\\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)$$\n\n", "\n", "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.\n", "\n", "If we assume harmonic transition state theory, then we can write the site entropic term $\\rho^0$ as\n", "$$\\rho^0 = \\frac{\\prod \\nu^{\\text{perfect-supercell}}}{\\prod \\nu^{\\text{defect-supercell}}}$$\n\n", "\n", "where $\\nu$ are the vibrational eigenvalues of the corresponding supercells, and the prefactor for the transition state is\n", "$$\\nu^\\text{T} = \\frac{\\prod \\nu^{\\text{perfect-supercell}}}{\\prod_{\\nu^2>0} \\nu^{\\text{transition state}}}$$\n\n", "\n", "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$.\n", "\n", "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}$.\n", "\n", "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.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, interstitial diffusion\n", "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)." ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on method diffusivity in module onsager.OnsagerCalc:\n", "\n", "diffusivity(pre, betaene, preT, betaeneT, CalcDeriv=False) method of onsager.OnsagerCalc.Interstitial instance\n", " Computes the diffusivity for our element given prefactors and energies/kB T.\n", " Also returns the negative derivative of diffusivity with respect to beta (used to compute\n", " the activation barrier tensor) if CalcDeriv = True\n", " The input list order corresponds to the sitelist and jumpnetwork\n", " \n", " :param pre: list of prefactors for unique sites\n", " :param betaene: list of site energies divided by kB T\n", " :param preT: list of prefactors for transition states\n", " :param betaeneT: list of transition state energies divided by kB T\n", " :return D[3,3]: diffusivity as a 3x3 tensor\n", " :return DE[3,3]: diffusivity times activation barrier (if CalcDeriv == True)\n", "\n" ] } ], "source": [ "help(FCCintdiffuser.diffusivity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "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`." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'states': [['i:+0.500,+0.500,+0.500'],\n", " ['i:+0.750,+0.750,+0.750', 'i:+0.250,+0.250,+0.250']],\n", " 'transitions': [['i:+0.500,+0.500,+0.500^i:+0.250,+1.250,+0.250',\n", " 'i:+0.250,+0.250,+0.250^i:+0.500,-0.500,+0.500',\n", " 'i:+0.500,+0.500,+0.500^i:+0.750,-0.250,+0.750',\n", " 'i:+0.750,+0.750,+0.750^i:+0.500,+1.500,+0.500',\n", " 'i:+0.500,+0.500,+0.500^i:+1.250,+0.250,+0.250',\n", " 'i:+0.250,+0.250,+0.250^i:-0.500,+0.500,+0.500',\n", " 'i:+0.500,+0.500,+0.500^i:+0.250,+0.250,+0.250',\n", " 'i:+0.250,+0.250,+0.250^i:+0.500,+0.500,+0.500',\n", " 'i:+0.500,+0.500,+0.500^i:+0.250,+0.250,+1.250',\n", " 'i:+0.250,+0.250,+0.250^i:+0.500,+0.500,-0.500',\n", " 'i:+0.500,+0.500,+0.500^i:-0.250,+0.750,+0.750',\n", " 'i:+0.750,+0.750,+0.750^i:+1.500,+0.500,+0.500',\n", " 'i:+0.500,+0.500,+0.500^i:+0.750,+0.750,-0.250',\n", " 'i:+0.750,+0.750,+0.750^i:+0.500,+0.500,+1.500',\n", " 'i:+0.500,+0.500,+0.500^i:+0.750,+0.750,+0.750',\n", " 'i:+0.750,+0.750,+0.750^i:+0.500,+0.500,+0.500'],\n", " ['i:+0.250,+0.250,+0.250^i:+0.750,+0.750,-0.250',\n", " 'i:+0.750,+0.750,+0.750^i:+0.250,+0.250,+1.250',\n", " 'i:+0.250,+0.250,+0.250^i:-0.250,+0.750,-0.250',\n", " 'i:+0.750,+0.750,+0.750^i:+1.250,+0.250,+1.250',\n", " 'i:+0.250,+0.250,+0.250^i:+0.750,-0.250,-0.250',\n", " 'i:+0.750,+0.750,+0.750^i:+0.250,+1.250,+1.250',\n", " 'i:+0.250,+0.250,+0.250^i:+0.750,-0.250,+0.750',\n", " 'i:+0.750,+0.750,+0.750^i:+0.250,+1.250,+0.250',\n", " 'i:+0.250,+0.250,+0.250^i:-0.250,+0.750,+0.750',\n", " 'i:+0.750,+0.750,+0.750^i:+1.250,+0.250,+0.250',\n", " 'i:+0.250,+0.250,+0.250^i:-0.250,-0.250,+0.750',\n", " 'i:+0.750,+0.750,+0.750^i:+1.250,+1.250,+0.250']]}" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "FCCintdiffuser.tags" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": true }, "outputs": [], "source": [ "FCCintdata = {\n", " 'i:+0.500,+0.500,+0.500': [1., 0.],\n", " 'i:+0.750,+0.750,+0.750': [2., 0.5],\n", " 'i:+0.500,+0.500,+0.500^i:+0.750,+0.750,-0.250': [10., 1.0],\n", " 'i:+0.750,+0.750,+0.750^i:+1.250,+1.250,+0.250': [50., 2.0]\n", "}" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.0, 2.0]\n", "[0.0, 2.0]\n", "[10.0, 50.0]\n", "[4.0, 8.0]\n" ] } ], "source": [ "# Conversion from dictionary to lists for a given kBT\n", "# We go through the tags in order, and find one in our data set.\n", "kBT = 0.25 # eV; a rather high temperature\n", "pre = [FCCintdata[t][0] for taglist in FCCintdiffuser.tags['states'] \n", " for t in taglist if t in FCCintdata]\n", "betaene = [FCCintdata[t][1]/kBT for taglist in FCCintdiffuser.tags['states'] \n", " for t in taglist if t in FCCintdata]\n", "preT = [FCCintdata[t][0] for taglist in FCCintdiffuser.tags['transitions'] \n", " for t in taglist if t in FCCintdata]\n", "betaeneT = [FCCintdata[t][1]/kBT for taglist in FCCintdiffuser.tags['transitions'] \n", " for t in taglist if t in FCCintdata]\n", "print(pre,betaene,preT,betaeneT,sep='\\n')" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 6.48557013e-02 1.30104261e-18 -1.30104261e-18]\n", " [ 1.30104261e-18 6.48557013e-02 1.30104261e-18]\n", " [ -1.30104261e-18 1.30104261e-18 6.48557013e-02]]\n", "[[ 2.35630632e-01 3.46944695e-18 -3.46944695e-18]\n", " [ 3.46944695e-18 2.35630632e-01 3.46944695e-18]\n", " [ -3.46944695e-18 3.46944695e-18 2.35630632e-01]]\n" ] } ], "source": [ "DFCCint, dDFCCint = FCCintdiffuser.diffusivity(pre, betaene, preT, betaeneT, CalcDeriv=True)\n", "print(DFCCint, dDFCCint, sep='\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interpretation of this output will be described below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Face-centered cubic crystal, vacancy mediated-diffusion\n", "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\"." ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on method Lij in module onsager.OnsagerCalc:\n", "\n", "Lij(bFV, bFS, bFSV, bFT0, bFT1, bFT2, large_om2=100000000.0) method of onsager.OnsagerCalc.VacancyMediated instance\n", " Calculates the transport coefficients: L0vv, Lss, Lsv, L1vv from the scaled free energies.\n", " The Green function entries are calculated from the omega0 info. As this is the most\n", " time-consuming part of the calculation, we cache these values with a dictionary\n", " and hash function.\n", " \n", " :param bFV[NWyckoff]: beta*eneV - ln(preV) (relative to minimum value)\n", " :param bFS[NWyckoff]: beta*eneS - ln(preS) (relative to minimum value)\n", " :param bFSV[Nthermo]: beta*eneSV - ln(preSV) (excess)\n", " :param bFT0[Nomega0]: beta*eneT0 - ln(preT0) (relative to minimum value of bFV)\n", " :param bFT1[Nomega1]: beta*eneT1 - ln(preT1) (relative to minimum value of bFV + bFS)\n", " :param bFT2[Nomega2]: beta*eneT2 - ln(preT2) (relative to minimum value of bFV + bFS)\n", " :param large_om2: threshold for changing treatment of omega2 contributions (default: 10^8)\n", " :return Lvv[3, 3]: vacancy-vacancy; needs to be multiplied by cv/kBT\n", " :return Lss[3, 3]: solute-solute; needs to be multiplied by cv*cs/kBT\n", " :return Lsv[3, 3]: solute-vacancy; needs to be multiplied by cv*cs/kBT\n", " :return Lvv1[3, 3]: vacancy-vacancy correction due to solute; needs to be multiplied by cv*cs/kBT\n", "\n" ] } ], "source": [ "help(fivefreqdiffuser.Lij)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function preene2betafree in module onsager.OnsagerCalc:\n", "\n", "preene2betafree(kT, preV, eneV, preS, eneS, preSV, eneSV, preT0, eneT0, preT1, eneT1, preT2, eneT2, **ignoredextraarguments)\n", " Read in a series of prefactors (exp(S/k_B)) and energies, and return\n", " beta F for energies and transition state energies. Used to provide scaled values\n", " to Lij().\n", " Can specify all of the entries using a dictionary; e.g., ``preene2betafree(kT, **data_dict)``\n", " and then send that output as input to Lij: ``Lij(*preene2betafree(kT, **data_dict))``\n", " (we ignore extra arguments so that a dictionary including additional entries can be passed)\n", " \n", " :param kT: temperature times Boltzmann's constant kB\n", " :param preV: prefactor for vacancy formation (prod of inverse vibrational frequencies)\n", " :param eneV: vacancy formation energy\n", " :param preS: prefactor for solute formation (prod of inverse vibrational frequencies)\n", " :param eneS: solute formation energy\n", " :param preSV: excess prefactor for solute-vacancy binding\n", " :param eneSV: solute-vacancy binding energy\n", " :param preT0: prefactor for vacancy transition state\n", " :param eneT0: energy for vacancy transition state (relative to eneV)\n", " :param preT1: prefactor for vacancy swing transition state\n", " :param eneT1: energy for vacancy swing transition state (relative to eneV + eneS + eneSV)\n", " :param preT2: prefactor for vacancy exchange transition state\n", " :param eneT2: energy for vacancy exchange transition state (relative to eneV + eneS + eneSV)\n", " :return bFV: beta*eneV - ln(preV) (relative to minimum value)\n", " :return bFS: beta*eneS - ln(preS) (relative to minimum value)\n", " :return bFSV: beta*eneSV - ln(preSV) (excess)\n", " :return bFT0: beta*eneT0 - ln(preT0) (relative to minimum value of bFV)\n", " :return bFT1: beta*eneT1 - ln(preT1) (relative to minimum value of bFV + bFS)\n", " :return bFT2: beta*eneT2 - ln(preT2) (relative to minimum value of bFV + bFS)\n", "\n" ] } ], "source": [ "help(fivefreqdiffuser.preene2betafree)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even this is a bit complicated; so we use an additional function that maps the tags into the appropriate lists, `tags2preene()`:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on method tags2preene in module onsager.OnsagerCalc:\n", "\n", "tags2preene(usertagdict, VERBOSE=False) method of onsager.OnsagerCalc.VacancyMediated instance\n", " Generates energies and prefactors based on a dictionary of tags.\n", " \n", " :param usertagdict: dictionary where the keys are tags, and the values are tuples: (pre, ene)\n", " :param VERBOSE: (optional) if True, also return a dictionary of missing tags, duplicate tags, and bad tags\n", " :return thermodict: dictionary of ene's and pre's corresponding to usertagdict\n", " :return missingdict: dictionary with keys corresponding to tag types, and the values are\n", " lists of lists of symmetry equivalent tags that are missing\n", " :return duplicatelist: list of lists of tags in usertagdict that are (symmetry) duplicates\n", " :return badtaglist: list of all tags in usertagdict that aren't found in our dictionary\n", "\n" ] } ], "source": [ "help(fivefreqdiffuser.tags2preene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)." ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": true }, "outputs": [], "source": [ "fivefreqdata = {\n", " 'v:+0.000,+0.000,+0.000': [1., 0.],\n", " 's:+0.000,+0.000,+0.000': [1., 0.],\n", " 's:+0.000,+0.000,+0.000-v:+0.000,-1.000,+0.000': [1., -0.25],\n", " 'omega0:v:+0.000,+0.000,+0.000^v:+0.000,+1.000,+0.000': [10., 1.],\n", " '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],\n", " '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],\n", " '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],\n", " '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],\n", " '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':\n", " [10., 0.25]\n", "}" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0.18315639 -0. -0. ]\n", " [-0. 0.18315639 -0. ]\n", " [-0. -0. 0.18315639]]\n", "[[ 1.34561077 0. 0. ]\n", " [ 0. 1.34561077 0. ]\n", " [ 0. 0. 1.34561077]]\n", "[[-0.40965028 0. 0. ]\n", " [ 0. -0.40965028 0. ]\n", " [ 0. 0. -0.40965028]]\n", "[[ 6.24182702 0. 0. ]\n", " [ 0. 6.24182702 0. ]\n", " [ 0. 0. 6.24182702]]\n" ] } ], "source": [ "# Conversion from dictionary to lists for a given kBT\n", "# note that we can nest the mapping functions.\n", "kBT = 0.25 # eV; a rather high temperature\n", "fivefreqpreene = fivefreqdiffuser.tags2preene(fivefreqdata)\n", "fivefreqbetaF = fivefreqdiffuser.preene2betafree(kBT, **fivefreqpreene)\n", "L0vv, Lss, Lsv, L1vv = fivefreqdiffuser.Lij(*fivefreqbetaF)\n", "print(L0vv, Lss, Lsv, L1vv, sep='\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interpretation of this output will be described below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpretation of output\n", "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.\n", "\n", "There are two underlying definitions that we use to define our transport coefficients:\n", "$$\\mathbf j = -\\underline D \\nabla c$$\n\n", "\n", "defines the *solute diffusivity* as the tensorial transport coefficient that relates defect concentration gradients to defect fluxes, and\n", "$$\\mathbf j^\\text{s} = -\\underline L^\\text{ss}\\nabla\\mu^\\text{s} - \\underline L^\\text{sv}\\nabla\\mu^\\text{v}$$\n\n", "\n", "$$\\mathbf j^\\text{v} = -\\underline L^\\text{vv}\\nabla\\mu^\\text{v} - \\underline L^\\text{sv}\\nabla\\mu^\\text{s}$$\n\n", "\n", "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.\n", "\n", "Below are more specific details about the different calculators and the output available." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interstitial diffusivity\n", "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\n", "$$1\\text{ nm}^2\\cdot\\text{THz} = 10^{-6}\\text{ m}^2/\\text{s} = 10^{-2}\\text{ cm}^2/\\text{s}$$\n\n", "\n", "$$1\\text{ A}^2\\cdot\\text{THz} = 10^{-8}\\text{ m}^2/\\text{s} = 10^{-4}\\text{ cm}^2/\\text{s}$$\n\n", "\n", "\n", "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\n", "$$D_\\text{B} = \\left\\{[\\text{B}_\\text{i}]D_\\text{int} + [\\text{B}_\\text{A}]D_\\text{sub}\\right\\}/[\\text{B}]$$\n\n", "\n", "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\n", "$$[\\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)$$\n\n", "\n", "and\n", "$$[\\text{B}_\\text{A}]/[\\text{B}] = (1 + \\exp(-(E_\\text{int}-E_\\text{sub})/k_\\text{B}T)^{-1} \\approx 1$$\n\n", "\n", "where the approximations are valid when $E_\\text{int}-E_\\text{sub}\\gg k_\\text{B}T$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Derivatives of diffusivity: activation barrier tensor\n", "At any given temperature, the temperature dependence of the diffusivity can be taken as an Arrhenius form,\n", "$$\\underline D = \\underline D_0 \\exp(-\\beta \\underline E^\\text{act})$$\n\n", "\n", "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.\n", "We can compute $Q$ by taking the per-component logarithmic derivative with respect to inverse temperature,\n", "$$\\underline E^\\text{act} = -\\underline D^{-1/2}\\frac{d\\underline D}{d\\beta}\\underline D^{-1/2}$$\n\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 9.08288044e-01 -4.84706356e-18 4.84706356e-18]\n", " [ -4.84706356e-18 9.08288044e-01 -4.84706356e-18]\n", " [ 4.84706356e-18 -4.84706356e-18 9.08288044e-01]]\n" ] } ], "source": [ "print(np.dot(np.linalg.inv(DFCCint), kBT*dDFCCint))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "This tensor has the same energy units as the variable `kBT`.\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Derivatives of diffusivity: elastodiffusion and activation volume tensor\n", "The derivative with respect to strain is the fourth-rank *elastodiffusivity* tensor $\\underline d$, where\n", "$$d_{abcd} = \\frac{dD_{ab}}{d\\varepsilon_{cd}}$$\n\n", "\n", "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.\n", "\n", "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,\n", "$$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}$$\n\n", "\n", "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\n", "$$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$$\n\n", "\n", "can be useful." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Vacancy-mediated diffusivity\n", "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:\n", "$$1\\text{ nm}^2\\cdot\\text{THz} = 10^{-6}\\text{ m}^2/\\text{s} = 10^{-2}\\text{ cm}^2/\\text{s}$$\n\n", "\n", "$$1\\text{ A}^2\\cdot\\text{THz} = 10^{-8}\\text{ m}^2/\\text{s} = 10^{-4}\\text{ cm}^2/\\text{s}$$\n\n", "\n", "To convert the four quantities into $\\underline L^\\text{vv}$, $\\underline L^\\text{ss}$, and $\\underline L^\\text{sv}$, some additional information is required. \n", "\n", "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,\n", "\n", "* $\\underline L^\\text{ss}$ = `Lss*(solute concentration)*(vacancy concentration)/(volume)/kBT`\n", "* $\\underline L^\\text{sv}$ = `Lsv*(solute concentration)*(vacancy concentration)/(volume)/kBT`\n", "\n", "where the concentration quantities are *fractional*.\n", "\n", "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,\n", "\n", "* $\\underline L^\\text{vv}$ = `(L0vv + L1vv*(solute concentration))*(vacancy concentration)/(volume)/kBT`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Drag ratio\n", "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" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-0.30443445 0. 0. ]\n", " [ 0. -0.30443445 0. ]\n", " [ 0. 0. -0.30443445]]\n" ] } ], "source": [ "print(np.dot(Lsv, np.linalg.inv(Lss)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The vacancy wind factor $G=\\underline L^\\text{As}(\\underline L^\\text{ss})^{-1}$ is related to the drag ratio by simple transformations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solute diffusivity in the dilute limit\n", "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\n", "$$\\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)$$\n\n", "\n", "where $\\Omega$ is the volume per atom and $\\gamma^\\text{s}$ is the solute activity:\n", "$$\\mu^\\text{s} = \\mu^\\text{s}_0 + k_\\text{B}T\\ln\\left(\\gamma^\\text{s}c^\\text{s}/c^\\text{s}_0\\right)$$\n\n", "\n", "In the dilute limit, $\\gamma^\\text{s}\\to 1$, and thus\n", "$$\\underline D^\\text{s} = k_\\text{B}T\\Omega(c^\\text{s})^{-1}\\underline L^\\text{ss}$$\n\n", "\n", "Conveniently, this cancels most of the \"missing\" prefactors we put in to compute the Onsager coefficient; hence,\n", "\n", "* $\\underline D^\\text{s}$ = `Lss*(vacancy concentration)`\n", "\n", "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)$.\n", "\n", "A similar argument holds for the vacancy diffusivity in the dilute limit\n", "\n", "* $\\underline D^\\text{v}$ = `L0vv + (solute concentration)*L1vv`\n", "\n", "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." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }