Silicon, Aluminum, Iron, and Graphene are characterized with a popular Fortran implementation of the self-consistent density functional theory -- Quantum ESPRESSO. An introduction to DFT and using QE is given, and results for Silicon are compared with those of previous calculations using the empirical pseudopotential and tight-binding methods.

Self-Consistent Density Functional Theory Calculations with Quantum ESPRESSO.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<h1>Self-Consistent Density Functional Theory Calculations with Quantum ESPRESSO</h1>\n",
    "<br> William M. Medlar <br>\n",
    "*Department of Physics, University of North Texas, Denton, Texas*\n",
    "<br><br> Silicon, Aluminum, Iron, and Graphene are characterized with a popular Fortran implementation of the self-consistent density functional <br> theory -- Quantum ESPRESSO. An introduction to DFT and using QE is given, and results for Silicon are compared with those of previous <br> calculations using the empirical pseudopotential and tight-binding methods.<br>\n",
    "    \n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I. Introduction\n",
    "\n",
    "The Quantum ESPRESSO software package is a huge, state-of-the-art, open-source library of high-performance code for use with quantum chemistry calculations. In this paper I will introduce the plane wave self-consistent field (pw.x) program to calculate the electronic ground state of Silicon, then expand to metals with Aluminum, spin-polarized systems with Iron, and finish off with a more complicated example of Graphene."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. Density Functional Theory\n",
    "\n",
    "The electronic ground state of a solid determines the bonding characteristics of its nuclei$^1$. Materials with a completely full valence shell (\"closed-shell systems\") are in very stable configurations and tend to form closely packed solids -- these are, relatively, easy to study. On the other hand, open-shell systems do not have their electrons on a tight leash, and thus are much harder to computationally model.\n",
    "\n",
    "The electronic ground state here is a many body problem, in which each electron and nuclei experiences and affects every other particle in the system, and thus are often approximated using methods such as pseudopotentials. Density functional theory, on the other hand, is an exact theory of a many body system.\n",
    "\n",
    "DFT applies to any system under an external potential $V_{ext}(\\vec{r})$ with a Hamiltonian:\n",
    "\n",
    "$$\n",
    "\\hat{H} = - \\frac{\\hbar^2}{2 m_e} \\sum\\limits_i \\nabla^2 + \\sum\\limits_i V_{ext}(\\vec{r_i}) + \\frac{1}{2} \\sum\\limits_{i \\neq j} \\frac{e^2}{\\| \\vec{r_i} - \\vec{r_j} \\|}\n",
    "$$\n",
    "\n",
    "DFT finds its roots in the Hohenberg and Kohn theorems:\n",
    "\n",
    "1. If two systems of electrons have the same ground state density $n(\\vec{r})$, the difference between their potentials must be constant.\n",
    "2. The functional of the ground state energy of the system $E[n]$ can be defined in terms of the density $n(\\vec{r})$ for any external potential $V_{ext}(\\vec{r_i})$.\n",
    "\n",
    "See Clark$^2$ for an easy to follow proof of these two theorems.\n",
    "\n",
    "The trouble here lies in finding this functional, denoted $F_{HK}[n]$. Luckily Kohn and Sham$^3$ devised an ingenious method that allows us to solve a completely different problem, yet retain the exactness of DFT. Their assumption was that the original interacting system of electrons was equal to a fictitious non-interacting system with the complicated exchange and correlation parameters wrapped up in a neat approximate functional of the density -- the local density approximation.\n",
    "\n",
    "Now we've got a full set of equations:\n",
    "\n",
    "1. The density of the fictitious system: $n(\\vec{r}) = \\sum\\limits_{s} \\sum\\limits_{i=1, N} \\|\\Psi_i^s(\\vec{r})\\|^2$\n",
    "2. The non-interacting kinetic energy: $T_s = \\frac{1}{2} \\sum\\limits_{s} \\sum\\limits_{i=1, N} \\|\\nabla \\Psi_i^s(\\vec{r})\\|^2$\n",
    "3. The classic electronic Coulomb (Hartree) energy: $E_{\\text{Hartree}}[n] = \\frac{1}{2} \\iint d^3 rd^3 r' \\frac{ n(\\vec{r}) \\: n(\\vec{r'}) }{\\| \\vec{r} - \\vec{r'} \\|}$\n",
    "4. And finally the full HK functional: $E_{KS}[n] = T_s[n] + \\int d^3r \\: V_{ext} \\: (\\vec{r}) n(\\vec{r}) + E_{\\text{Hartree}}[n] + E_{||} + E_{XC}[n]$\n",
    "\n",
    "To solve this beastly functional we can view it as a problem of minimization -- take a guess as to the density, vary the wavefunctions to end up with an eigenvalue problem involving the Kohn-Sham Hamiltonian, and solve for the eigenvalues using standard numerical routines. From here we can derive the new electron density and run the calculation again until it is self-consistent.\n",
    "\n",
    "$$\n",
    "H_{KS}^{\\sigma}(\\vec{r}) = -\\frac{1}{2} \\nabla^2 + V_{KS}^{\\sigma}(\\vec{r}) \\:\\:\\: \\text{where}\n",
    "\\\\ V_{KS}^{\\sigma}(\\vec{r}) = V_{ext}(\\vec{r}) + \\frac{\\delta E_{\\text{Hartree}}}{\\delta n (\\vec{r}, \\vec{\\sigma})} + \\frac{\\delta E_{XC}}{\\delta n (\\vec{r}, \\vec{\\sigma})} = V_{ext}(\\vec{r}) + V_{\\text{Hartree}}(\\vec{r}) + V_{XC}^{\\sigma}(\\vec{r})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Implementation\n",
    "\n",
    "The hardest part, converting our variational equations to computable code, has already been done for us. It's simply a matter of properly manipulating the input files for Quantum ESPRESSO programs from here on out. To begin, we'll start with something simple, like..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Silicon\n",
    "\n",
    "Our first question is \"how is the crystal structure defined?\". Silicon is a homogenous face-centered cubic lattice with two atoms per irreducable unit cell. Under the ```&system``` section, this corresponds to:\n",
    "- ```ibrav = 2``` to signify fcc\n",
    "- ```celldm(1) = 10.2``` only one lattice (dimensional) parameter (in au) is necessary to define our Bravais lattice\n",
    "- ```nat = 2``` for two atoms per unit cell\n",
    "- ```ntyp = 1``` as our Silicon lattice is homogenous (one atomic species)\n",
    "- ```ecutwfc = 12``` we'll get to this parameter in a moment\n",
    "\n",
    "The ```ATOMIC_SPECIES``` section contains our atomic species (Si), its mass in amu (28.086), and the filename of our pseudopotential (Si.pz-vbc.UPF). The section below that, ```ATOMIC_POSITIONS alat```, contains the locations of the atoms in the unit cell in terms of our lattice constant. For Silicon that's:\n",
    "```\n",
    " Si 0.00 0.00 0.00 \n",
    " Si 0.25 0.25 0.25 \n",
    "```\n",
    "\n",
    "Finally we reach the first parameter we'll vary: our k-point grid. Currently it should be set to \n",
    "```\n",
    "K_POINTS automatic\n",
    " 4 4 4 1 1 1\n",
    "```\n",
    "\n",
    "The ```automatic 4 4 4``` corresponds to an automatically generated (via a Monkhorst-Pack grid) finite number of sample points (in this case 4) in the first Brillouin zone to measure. We want to ensure that we have a sufficiently dense grid of k-points to account for the atom's periodicity, however every k-point we specify will be used to solve the Kohn-Sham Hamiltonian, increasing computation time, so we must find a happy medium.\n",
    "\n",
    "Run the plane wave self-consistent field program with\n",
    "\n",
    "```\n",
    "pw.x -in <input file> > <output file>\n",
    "```  \n",
    "\n",
    "taking note of the line  \n",
    "\n",
    "```!    total energy    =    -15.79449593 Ry```.  \n",
    "\n",
    "This is our final ground state energy after iterating until the density calculation became self-consistent. From here we can vary our number of k-points from 2, 4, 6, ... until we acheive a sufficient level of convergence. Fig. 1. contains a plot of the number of k-points versus the total energy -- it is very clear where the k-points converge and additional points no longer have much effect on accuracy.\n",
    "\n",
    "K-points aren't the only parameter we should test for convergence, we can also modify our kinetic energy cutoff ```ecutwfc```. This parameter is responsible for determining the size of our plane wave basis, with a value given in Rydbergs. Perform the same convergence study by incrementing ```ecutwfc``` from 12, 16, 20, ... -- this will probably take significantly longer than the k-point study as the size of the plane wave basis is responsible for the size of the Hamiltonian matrix. Fig. 2. contains a plot of ```ecutwfc``` values agains the total energy; again it is easy to see where we should set our parameters.\n",
    "\n",
    "The last parameter we can adjust is our lattice parameter, ```celldm(1)```, although this is not a converging value. Adjust the k-point grid and plane wave basis to our converged parameters and plot incremental values of ```celldm(1)``` from 9.8, 9.9, ... , 10.7 au. Fig. 3. contains a plot of the lattice parameter versus the total energy. Our true lattice parameter will lie where the lattice is most stable, i.e., has the lowest energy configuration. You can use the ```ev.x``` program with the calculated energies to extract an accurate value for the lattice parameter. Set this as your ```celldm(1)``` before continuing.\n",
    "\n",
    "To determine the band structure of Silicon, we can use the same ```pw.x``` program with a few adjusted parameters:\n",
    "- under ```&control``` set ```calculation = 'bands' ``` to perform a band structure calculation\n",
    "- under ```&system``` specify the number of bands to calculate with ```nbnd = 8```\n",
    "- set ```KPOINTS``` to ```tpiba```  to specify a k-point grid in units of $\\frac{2 \\pi}{a}$with a k-point grid along the high symmetry path $L \\rightarrow \\Gamma \\rightarrow X \\rightarrow K \\rightarrow \\Gamma$\n",
    "\n",
    "A plottable file can be generated using the ```bands.x``` program with the ```&bands > filband``` parameter set to the output file name and ```lsym = .false.``` Use the interactive program ```plotband.x``` with the output file. Fig. 4. shows the calculated band structure of Silicon with converged parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aluminum\n",
    "\n",
    "Aluminum has a simpler lattice than Silicon with only one atom per unit cell, but don't let that fool you, we have to accout for the fact that Aluminum is a metal. Valence bands and a sparse k-grid will not suffice. Again we should seek k-point convergence; results on Fig. 5.\n",
    "\n",
    "We're met with a few new variables: ```occupations```, ```smearing```, and ```degauss```. ```smearing``` can take on values of \"gauss\", \"marzari-vanderbilt\", or \"methfessel-paxton\". These options are various algorithms to smear out the Fermi level of the metal, with the ```degauss``` variable as a measure of the smearing strength, regardless of the algorithm used. Degauss is commonly set somewhere between 0.1 and 0.5, but we'll test out a few and measure convergence. See Fig. 6. for the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iron\n",
    "\n",
    "Iron is doubly unique compared to our previous materials. Not only is it magnetic, but it also requires a special kind of pseudopotential as its localized 3D atomic states are very hard. This special pseudopotential is called \"Ultrasoft PP\". With Ultrasoft we'll have to manually set our ```ecutrho```, as it should be much larger than the default value we've been using this whole time -- 4x ```ecutwfc```. We must also manually set the number of bands, ```nbnd = 8```. As always, we test for k-convergence. See Fig. 7. for results.\n",
    "\n",
    "In the non-magnetic case the number of k-points is doubled, one set for the spin up and one for the spin down states.. You can check this by setting ```verbosity``` to \"high\" to see all the output pw.x has to offer.\n",
    "\n",
    "Take note of a new variable, ```starting_magnetization``` and a new line in the return value of pw.x, \"total\" and \"absolute magnetization\". Since we only have one magnetic atom per unit cell, the only possible structure here is ferromagnetic. To reach antiferromagnetic, you'll have to expand your unit cell to include a *second* Iron atom with an equal and opposite starting magnetization. This corresponds to an input file containing:\n",
    "```\n",
    "&system\n",
    " [...]\n",
    " ibrav = 1,\n",
    " nat = 2,\n",
    " starting_magnetization(1) =  0.6,\n",
    " starting_magnetization(2) = -0.6\n",
    "ATOMIC_POSITIONS crystal\n",
    " Fe 0.0 0.0 0.0,\n",
    " Fe 0.5 0.5 0.5\n",
    "```\n",
    "Note that we still have one type of atom (```ntyp = 1```), but we now have two atoms sharing a unit cell, both with equal and opposite magnetization.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV. Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Silicon\n",
    "\n",
    "K-point convergence at 0.005% was acheived with ```nk = 6```, while the plane wave basis acheived 0.005% convergence at ```ecutwfc = 32```. Fig. 1. contains a plot of nk versus total energy until the point of convergence; Fig. 2. contains a plot of ecutwfc in Ry versus total energy until the point of convergence. The calculated lattice parameter was found to be 10.2087 a.u., a slight divergence from the standard value of 10.26 a.u. With only a handful of data points, it's fair to say this may not be wholly accurate. The lattice parameter, celldm(1) is plotted against total energy in Fig. 3. The band structure is given in Fig. 4. \n",
    "\n",
    "\n",
    "<img src=\"http://i.imgur.com/8XU3tSB.png\">\n",
    "$$\\scriptsize \\text{Fig. 1. Convergence of k-vector grid for Silicon.}$$\n",
    "\n",
    "<img src=\"http://i.imgur.com/bcMNhLT.png\">\n",
    "$$\\scriptsize \\text{Fig. 2. Convergence of plane waves for Silicon.}$$\n",
    "\n",
    "<img src=\"http://i.imgur.com/liNF2pL.png\">\n",
    "$$\\scriptsize \\text{Fig. 3. Lattice constant versus energy configuration of Silicon.}$$\n",
    "\n",
    "<img src=\"http://i.imgur.com/4VURRKg.png\">\n",
    "$$\\scriptsize \\text{Fig. 4. plotband.x-derived band structure of Silicon with nk = 6 and ecutwfc = 32.}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aluminum\n",
    "\n",
    "K-point convergence to 0.005% was found to be nk = 10 for all three smearing methods, see Fig. 5. Large (> 0.5 Ry) values of ```degauss``` cause divergence among the energy, as the Fermi level becomes too smeared. Methfessel-Paxton diverged the least with high degauss, and should be the go-to algorithm for any calculations that require large amounts of smearing without loss of generality. For standard metals, however, I'd suggest sticking to a degauss value $\\leq$ 0.5.\n",
    "\n",
    "<img src=\"http://i.imgur.com/FmtL0DK.png\">\n",
    "$$\\scriptsize \\text{Fig. 5. Convergence of k-vector grid for Aluminum; Gauss (blue), M-V (orange), M-P (green)}$$\n",
    "\n",
    "<img src=\"http://i.imgur.com/3ApqdIy.png\">\n",
    "\n",
    "$$\\scriptsize \\text{Fig. 6. Measured range of degauss parameter for Aluminum, note that the energy diverges for larger values} \\\\ \\scriptsize \\text{Gauss (blue), M-V (orange), M-P (green)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iron\n",
    "\n",
    "Iron converged within 0.005% surprisingly quickly for k-points; as a metal I expected it to require much more than nk = 8. Fig. 7. is deceiving, be sure to look closely at the y values -- a change in K hardly affects the energy of Iron!\n",
    "\n",
    "<img src=\"http://i.imgur.com/rPmBUMi.png\">\n",
    "$$\\scriptsize \\text{Fig. 7. Convergence of k-vector grid for Iron.}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## V. References\n",
    "\n",
    "1. M.B. Nardelli, \"An introduction to Density Functional Theory and Simulations,\" http://sites.psu.edu/qe2014/wp-content/uploads/sites/9818/2014/06/Intro-to-DFT.pdf, downloaded November, 2015\n",
    "2. S.J. Clark, \"The Hohenberg-Kohn Theorems,\" http://cmt.dur.ac.uk/sjc/thesis_ppr/node12.html, downloaded November, 2015\n",
    "3. W. Kohn, L.J. Sham, \"Self-Consistent Equations Including Exchange and Correlation Effects,\" Phys. Rev. 140, A1133 (1965)"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

Note that this is not the full set of exercises, and will be updated to include them as I complete them.