<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fe-C diffusion and elastodiffusivity\n",
    "Taking data from <NAME>, <NAME>, <NAME>, <NAME> and <NAME>, Acta mater. **59** (2011) p. 6963 [doi:10.1016/j.actamat.2011.07.048](http://dx.doi.org/10.1016/j.actamat.2011.07.048)\n",
    "\n",
    "Fe in the body-centered cubic phase, $a_0 = 0.28553\\text{ nm}$; C sit at octahedral sites, where the transition states between octahedral sites are represented by tetrahedral sites. The data is obtained from an EAM potential, where $C_{11} = 243\\text{ GPa}$, $C_{12}=145\\text{ GPa}$, and $C_{44} = 116\\text{ GPa}$. The tetrahedral transition state is 0.816 eV above the octahedral site, and the attempt frequency is taken as 10 THz ($10^{13}\\text{ Hz}$).\n",
    "\n",
    "The dipole tensors can be separated into *parallel* and *perpendicular* components; the parallel direction points towards the closest Fe atoms for the C, while the perpendicular components lie in the interstitial plane. For the octahedral, the parallel component is 8.03 eV, and the perpendicular is 3.40 eV; for the tetrahedral transition state, the parallel component is 4.87 eV, and the perpendicular is 6.66 eV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../')\n",
    "import numpy as np\n",
    "import onsager.crystal as crystal\n",
    "import onsager.OnsagerCalc as onsager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create BCC lattice (lattice constant in nm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Lattice:\n",
      "  a1 = [-0.142765  0.142765  0.142765]\n",
      "  a2 = [ 0.142765 -0.142765  0.142765]\n",
      "  a3 = [ 0.142765  0.142765 -0.142765]\n",
      "#Basis:\n",
      "  (Fe) 0.0 = [ 0.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "a0 = 0.28553\n",
    "Fe = crystal.Crystal.BCC(a0, \"Fe\")\n",
    "print(Fe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add carbon interstitial sites at octahedral sites in the lattice. This code (1) gets the set of symmetric Wyckoff positions corresponding to the single site $[00\\frac12]$ (first translated into unit cell coordinates), and then adds that new basis to our Fe crystal to generate a *new* crystal structure, that we name \"FeC\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Lattice:\n",
      "  a1 = [-0.142765  0.142765  0.142765]\n",
      "  a2 = [ 0.142765 -0.142765  0.142765]\n",
      "  a3 = [ 0.142765  0.142765 -0.142765]\n",
      "#Basis:\n",
      "  (Fe) 0.0 = [ 0.  0.  0.]\n",
      "  (C) 1.0 = [ 0.   0.5  0.5]\n",
      "  (C) 1.1 = [ 0.5  0.5  0. ]\n",
      "  (C) 1.2 = [ 0.5  0.   0.5]\n"
     ]
    }
   ],
   "source": [
    "uoct = np.dot(Fe.invlatt, np.array([0, 0, 0.5*a0]))\n",
    "FeC = Fe.addbasis(Fe.Wyckoffpos(uoct), [\"C\"])\n",
    "print(FeC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we construct a *diffuser* based on our interstitial. We need to create a `sitelist` (which will be the Wyckoff positions) and a `jumpnetwork` for the transitions between the sites. There are tags that correspond to the unique states and transitions in the diffuser. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diffuser for atom 1 (C)\n",
      "#Lattice:\n",
      "  a1 = [-0.142765  0.142765  0.142765]\n",
      "  a2 = [ 0.142765 -0.142765  0.142765]\n",
      "  a3 = [ 0.142765  0.142765 -0.142765]\n",
      "#Basis:\n",
      "  (Fe) 0.0 = [ 0.  0.  0.]\n",
      "  (C) 1.0 = [ 0.   0.5  0.5]\n",
      "  (C) 1.1 = [ 0.5  0.5  0. ]\n",
      "  (C) 1.2 = [ 0.5  0.   0.5]\n",
      "states:\n",
      "i:+0.000,+0.500,+0.500\n",
      "transitions:\n",
      "i:+0.000,+0.500,+0.500^i:+0.500,+1.000,+0.500\n",
      "\n"
     ]
    }
   ],
   "source": [
    "chem = 1  # 1 is the index corresponding to our C atom in the crystal\n",
    "sitelist = FeC.sitelist(chem)\n",
    "jumpnetwork = FeC.jumpnetwork(chem, 0.6*a0)  # 0.6*a0 is the cutoff distance for finding jumps\n",
    "FeCdiffuser = onsager.Interstitial(FeC, chem, sitelist, jumpnetwork)\n",
    "print(FeCdiffuser)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we assemble our data: the energies, prefactors, and dipoles for the C atom in Fe, matched to the *representative* states: these are the first states in the lists, which are also identified by the tags above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}<commit_msg>Fe-C notebook.<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fe-C diffusion and elastodiffusivity\n",
    "Taking data from <NAME>, <NAME>, <NAME>, <NAME> and C. Domain, Acta mater. **59** (2011) p. 6963 [doi:10.1016/j.actamat.2011.07.048](http://dx.doi.org/10.1016/j.actamat.2011.07.048)\n",
    "\n",
    "Fe in the body-centered cubic phase, $a_0 = 0.28553\\text{ nm}$; C sit at octahedral sites, where the transition states between octahedral sites are represented by tetrahedral sites. The data is obtained from an EAM potential, where $C_{11} = 243\\text{ GPa}$, $C_{12}=145\\text{ GPa}$, and $C_{44} = 116\\text{ GPa}$. The tetrahedral transition state is 0.816 eV above the octahedral site, and the attempt frequency is taken as 10 THz ($10^{13}\\text{ Hz}$).\n",
    "\n",
    "The dipole tensors can be separated into *parallel* and *perpendicular* components; the parallel direction points towards the closest Fe atoms for the C, while the perpendicular components lie in the interstitial plane. For the octahedral, the parallel component is 8.03 eV, and the perpendicular is 3.40 eV; for the tetrahedral transition state, the parallel component is 4.87 eV, and the perpendicular is 6.66 eV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../')\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "%matplotlib inline\n",
    "import onsager.crystal as crystal\n",
    "import onsager.OnsagerCalc as onsager\n",
    "from scipy.constants import physical_constants\n",
    "kB = physical_constants['Boltzmann constant in eV/K'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create BCC lattice (lattice constant in nm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Lattice:\n",
      "  a1 = [-0.142765  0.142765  0.142765]\n",
      "  a2 = [ 0.142765 -0.142765  0.142765]\n",
      "  a3 = [ 0.142765  0.142765 -0.142765]\n",
      "#Basis:\n",
      "  (Fe) 0.0 = [ 0.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "a0 = 0.28553\n",
    "Fe = crystal.Crystal.BCC(a0, \"Fe\")\n",
    "print(Fe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Elastic constants converted from GPa ($10^9$ J/m$^3$) to eV/(atomic volume)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S11 = 0.1023 S12 = -0.0382 S44 = 0.1187\n"
     ]
    }
   ],
   "source": [
    "stressconv = 1e9*1e-27*Fe.volume/physical_constants['electron volt'][0]\n",
    "c11, c12, c44 = 243*stressconv, 145*stressconv, 116*stressconv\n",
    "s11, s12, s44 = (c11+c12)/((c11-c12)*(c11+2*c12)), -c12/((c11-c12)*(c11+2*c12)), 1/c44\n",
    "print('S11 = {:.4f} S12 = {:.4f} S44 = {:.4f}'.format(s11, s12, s44))\n",
    "stensor = np.zeros((3,3,3,3))\n",
    "for a in range(3):\n",
    "    for b in range(3):\n",
    "        for c in range(3):\n",
    "            for d in range(3):\n",
    "                if a==b and b==c and c==d: stensor[a,b,c,d] = s11\n",
    "                elif a==b and c==d: stensor[a,b,c,d] = s12\n",
    "                elif (a==d and b==c) or (a==c and b==d): stensor[a,b,c,d] = s44/4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add carbon interstitial sites at octahedral sites in the lattice. This code (1) gets the set of symmetric Wyckoff positions corresponding to the single site $[00\\frac12]$ (first translated into unit cell coordinates), and then adds that new basis to our Fe crystal to generate a *new* crystal structure, that we name \"FeC\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Lattice:\n",
      "  a1 = [-0.142765  0.142765  0.142765]\n",
      "  a2 = [ 0.142765 -0.142765  0.142765]\n",
      "  a3 = [ 0.142765  0.142765 -0.142765]\n",
      "#Basis:\n",
      "  (Fe) 0.0 = [ 0.  0.  0.]\n",
      "  (C) 1.0 = [ 0.5  0.   0.5]\n",
      "  (C) 1.1 = [ 0.   0.5  0.5]\n",
      "  (C) 1.2 = [ 0.5  0.5  0. ]\n"
     ]
    }
   ],
   "source": [
    "uoct = np.dot(Fe.invlatt, np.array([0, 0, 0.5*a0]))\n",
    "FeC = Fe.addbasis(Fe.Wyckoffpos(uoct), [\"C\"])\n",
    "print(FeC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we construct a *diffuser* based on our interstitial. We need to create a `sitelist` (which will be the Wyckoff positions) and a `jumpnetwork` for the transitions between the sites. There are tags that correspond to the unique states and transitions in the diffuser. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diffuser for atom 1 (C)\n",
      "#Lattice:\n",
      "  a1 = [-0.142765  0.142765  0.142765]\n",
      "  a2 = [ 0.142765 -0.142765  0.142765]\n",
      "  a3 = [ 0.142765  0.142765 -0.142765]\n",
      "#Basis:\n",
      "  (Fe) 0.0 = [ 0.  0.  0.]\n",
      "  (C) 1.0 = [ 0.5  0.   0.5]\n",
      "  (C) 1.1 = [ 0.   0.5  0.5]\n",
      "  (C) 1.2 = [ 0.5  0.5  0. ]\n",
      "states:\n",
      "i:+0.500,+0.000,+0.500\n",
      "transitions:\n",
      "i:+0.000,+0.500,+0.500^i:+0.500,+0.500,+1.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "chem = 1  # 1 is the index corresponding to our C atom in the crystal\n",
    "sitelist = FeC.sitelist(chem)\n",
    "jumpnetwork = FeC.jumpnetwork(chem, 0.6*a0)  # 0.6*a0 is the cutoff distance for finding jumps\n",
    "FeCdiffuser = onsager.Interstitial(FeC, chem, sitelist, jumpnetwork)\n",
    "print(FeCdiffuser)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we assemble our data: the energies, prefactors, and dipoles for the C atom in Fe, matched to the *representative* states: these are the first states in the lists, which are also identified by the tags above.\n",
    "\n",
    "*A note about units:* If $\\nu_0$ is in THz, and $a_0$ is in nm, then $a_0^2\\nu_0 = 10^{-2}\\text{ cm}^2/\\text{s}$. Thus, we multiply by `Dconv` = $10^{-2}$ so that our diffusivity is output in cm<sup>2</sup>/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eneT: [ 0.816]\n",
      "dipoleT: [[ 6.66  0.    0.  ]\n",
      " [ 0.    4.87  0.  ]\n",
      " [ 0.    0.    6.66]]\n",
      "ene: [ 0.]\n",
      "pre: [ 1.]\n",
      "dipole: [[ 3.4   0.    0.  ]\n",
      " [ 0.    8.03  0.  ]\n",
      " [ 0.    0.    3.4 ]]\n",
      "preT: [ 0.1]\n"
     ]
    }
   ],
   "source": [
    "Dconv = 1e-2\n",
    "vu0 = 10*Dconv\n",
    "Etrans = 0.816\n",
    "dipoledict = {'Poctpara': 8.03, 'Poctperp': 3.40, \n",
    "              'Ptetpara': 4.87, 'Ptetperp': 6.66}\n",
    "FeCthermodict = {'pre': np.ones(len(sitelist)), 'ene': np.zeros(len(sitelist)),\n",
    "                 'preT': vu0*np.ones(len(jumpnetwork)), \n",
    "                 'eneT': Etrans*np.ones(len(jumpnetwork))}\n",
    "# now to construct the site and transition dipole tensors; we use a \"direction\"--either\n",
    "# the site position or the jump direction--to determine the parallel and perpendicular\n",
    "# directions.\n",
    "for dipname, Pname, direction in zip(('dipole', 'dipoleT'), ('Poct', 'Ptet'), \n",
    "                                     (np.dot(FeC.lattice, FeC.basis[chem][sitelist[0][0]]),\n",
    "                                      jumpnetwork[0][0][1])):\n",
    "    # identify the non-zero index in our direction:\n",
    "    paraindex = [n for n in range(3) if not np.isclose(direction[n], 0)][0]\n",
    "    Ppara, Pperp = dipoledict[Pname + 'para'], dipoledict[Pname + 'perp']\n",
    "    FeCthermodict[dipname] = np.diag([Ppara if i==paraindex else Pperp \n",
    "                                      for i in range(3)])\n",
    "for k,v in FeCthermodict.items():\n",
    "    print('{}: {}'.format(k, v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We look at the diffusivity $D$, the elastodiffusivity $d$, and the activation volume tensor $V$ over a range of temperatures from 300K to 1200K."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "<KEY>",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11027b518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Trange = np.linspace(300, 1200, 91)\n",
    "Tlabels = Trange[0::30]\n",
    "Dlist = []\n",
    "for T in Trange:\n",
    "    beta = 1./(kB*T)\n",
    "    D = FeCdiffuser.diffusivity(FeCthermodict['pre'], beta*FeCthermodict['ene'], \n",
    "                                FeCthermodict['preT'], beta*FeCthermodict['eneT'])\n",
    "    Dlist.append(D[0,0])\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(1./Trange, np.array(Dlist))\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_ylabel('$D$ [cm$^2$/s]')\n",
    "ax1.set_xlabel('$T^{-1}$ [K$^{-1}$]')\n",
    "ax2 = ax1.twiny()\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax2.set_xticks([1./t for t in Tlabels])\n",
    "ax2.set_xticklabels([\"{:.0f}K\".format(t) for t in Tlabels])\n",
    "ax2.set_xlabel('$T$ [K]')\n",
    "ax2.grid(False)\n",
    "ax2.tick_params(axis='x', top='on', direction='in', length=6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Activation volume. We plot the isotropic value (change in diffusivity with respect to pressure), but also the $V_{xxxx}$, $V_{xxyy}$, and $V_{xyxy}$ terms. Interestingly, the $V_{xxxx}$ term is negative---which indicates that diffusivity along the $[100]$ direction *increases* with compressive stress in the $[100]$ direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "<KEY>",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1102b20b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Trange = np.linspace(300, 1200, 91)\n",
    "Tlabels = Trange[0::30]\n",
    "Dlist = []\n",
    "Vlist = []\n",
    "for T in Trange:\n",
    "    beta = 1./(kB*T)\n",
    "    D, dD = FeCdiffuser.elastodiffusion(FeCthermodict['pre'], beta*FeCthermodict['ene'], \n",
    "                                        [beta*FeCthermodict['dipole']],\n",
    "                                        FeCthermodict['preT'], beta*FeCthermodict['eneT'], \n",
    "                                        [beta*FeCthermodict['dipoleT']])\n",
    "    Dlist.append(D[0,0])\n",
    "    Vtensor = (kB*T/(D[0,0]))*np.tensordot(dD, stensor, axes=((2,3),(0,1)))\n",
    "    Vlist.append([np.trace(np.trace(Vtensor))/9, \n",
    "                  Vtensor[0,0,0,0], Vtensor[0,0,1,1], Vtensor[0,1,0,1]])\n",
    "V = np.array(Vlist)\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(1./Trange, V[:,0], 'k', 1./Trange, V[:,1], 'b', \n",
    "         1./Trange, V[:,2], 'r', 1./Trange, V[:,3], 'g')\n",
    "ax1.set_yscale('linear')\n",
    "ax1.set_ylabel('$V$ [atomic volume]')\n",
    "ax1.set_xlabel('$T^{-1}$ [K$^{-1}$]')\n",
    "ax2 = ax1.twiny()\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax2.set_xticks([1./t for t in Tlabels])\n",
    "ax2.set_xticklabels([\"{:.0f}K\".format(t) for t in Tlabels])\n",
    "ax2.set_xlabel('$T$ [K]')\n",
    "ax2.grid(False)\n",
    "ax2.tick_params(axis='x', top='on', direction='in', length=6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "<KEY>IZE<KEY>",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11061a400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Trange, (V[:,1]-V[:,2])/V[:,3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
{'fraction_non_alphanumeric': 0.20576153360016852, 'fraction_numerical': 0.06140720455024226, 'mean_word_length': 2.3656504785537043, 'pattern_counts': {'":': 204, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4397288', 'n_tokens_mistral': 8269, 'n_tokens_neox': 7490, 'n_words': 1705}