v0.2.7/examples/ase.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Creating slabs with ASE\n",
    "\n",
    "ASE is short for the\n",
    "[atomistic simulation environment](https://wiki.fysik.dtu.dk/ase/index.html),\n",
    "a Python package to simplify the process of setting up, running and\n",
    "analysing results from atomistic simulations across different programs.\n",
    "Extremely powerful in this respect are the routines this code provides\n",
    "for setting up complicated systems (including surface-adsorption scenarios,\n",
    "defects, nanotubes, etc).\n",
    "See also the [ASE installation instructions](https://wiki.fysik.dtu.dk/ase/install.html).\n",
    "\n",
    "This example shows how to use ASE to setup a particular gallium arsenide surface\n",
    "and run the resulting calculation in DFTK.\n",
    "If you are less interested in having access to the full playground of options in DFTK,\n",
    "but more interested in performing analysis in ASE itself,\n",
    "have a look at [asedftk](https://github.com/mfherbst/asedftk).\n",
    "This package provides an ASE-compatible calculator class based on DFTK,\n",
    "such that one may write the usual Python scripts against ASE,\n",
    "but the calculations are still run in DFTK.\n",
    "\n",
    "The particular example we consider the (1, 1, 0) GaAs surface separated by vacuum\n",
    "with the setup slightly adapted from [^RCW2001].\n",
    "\n",
    "[^RCW2001]:\n",
    "    D. Raczkowski, A. Canning, and L. W. Wang\n",
    "    *Thomas-Fermi charge mixing for obtaining self-consistency in density functional calculations*\n",
    "    Phys. Rev. B **64**, 121101(R)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Parameters of the calculation. Since this surface is far from easy to converge,\n",
    "we made the problem simpler by choosing a smaller `Ecut` and smaller values\n",
    "for `n_GaAs` and `n_vacuum`.\n",
    "More interesting settings are `Ecut = 15` and `n_GaAs = n_vacuum = 20`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "miller = (1, 1, 0)   # Surface Miller indices\n",
    "n_GaAs = 2           # Number of GaAs layers\n",
    "n_vacuum = 4         # Number of vacuum layers\n",
    "Ecut = 5             # Hartree\n",
    "kgrid = (4, 4, 1);   # Monkhorst-Pack mesh"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "Use ASE to build the structure:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using PyCall\n",
    "\n",
    "ase_build = pyimport(\"ase.build\")\n",
    "a = 5.6537  # GaAs lattice parameter in Ångström (because ASE uses Å as length unit)\n",
    "gaas = ase_build.bulk(\"GaAs\", \"zincblende\", a=a)\n",
    "surface = ase_build.surface(gaas, miller, n_GaAs, 0, periodic=true);"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "Get the amount of vacuum in Ångström we need to add"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "d_vacuum = maximum(maximum, surface.cell) / n_GaAs * n_vacuum\n",
    "surface = ase_build.surface(gaas, miller, n_GaAs, d_vacuum, periodic=true);"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Write an image of the surface and embed it as a nice illustration:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "pyimport(\"ase.io\").write(\"surface.png\", surface * (3, 3, 1),\n",
    "                         rotation=\"-90x, 30y, -75z\")"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "<img src=\"surface.png\" width=500 height=500 />"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Use the `load_atoms` and `load_lattice` functions\n",
    "to convert to DFTK datastructures.\n",
    "These two functions not only support importing ASE atoms into DFTK,\n",
    "but a few more third-party datastructures as well.\n",
    "Typically the imported `atoms` use a bare Coulomb potential,\n",
    "such that appropriate pseudopotentials need to be attached in a post-step:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using DFTK\n",
    "\n",
    "atoms = load_atoms(surface)\n",
    "atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional=\"pbe\")) => position\n",
    "         for (el, position) in atoms]\n",
    "lattice = load_lattice(surface);"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "We model this surface with (quite large a) temperature of 0.01 Hartree\n",
    "to ease convergence. Try lowering the SCF convergence tolerance (`tol`\n",
    "or the `temperature` to see the full challenge of this system."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Free energy       Eₙ-Eₙ₋₁     ρout-ρin   Diag\n",
      "---   ---------------   ---------   --------   ----\n",
      "  1   -16.60904549584         NaN   2.60e-01    4.8 \n",
      "  2   -16.72296863664   -1.14e-01   8.41e-02    1.0 \n",
      "  3   -16.66425534939    5.87e-02   5.14e-02    2.3 \n",
      "  4   -16.73220156006   -6.79e-02   1.93e-02    2.0 \n",
      "  5   -16.73631937580   -4.12e-03   5.13e-03    1.8 \n",
      "  6   -16.73635765509   -3.83e-05   3.00e-03    3.0 \n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "model = model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe],\n",
    "                  temperature=0.001, smearing=DFTK.Smearing.Gaussian())\n",
    "basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\n",
    "\n",
    "scfres = self_consistent_field(basis, tol=1e-4, mixing=KerkerMixing());"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Energy breakdown:\n    Kinetic             5.8405521 \n    AtomicLocal         -105.5505203\n    AtomicNonlocal      2.3482321 \n    Ewald               35.5044300\n    PspCorrection       0.2016043 \n    Hartree             49.5126002\n    Xc                  -4.5931249\n    Entropy             -0.0001312\n\n    total               -16.736357655094\n"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "scfres.energies"
   ],
   "metadata": {},
   "execution_count": 7
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.4"
  },
  "kernelspec": {
   "name": "julia-1.5",
   "display_name": "Julia 1.5.4",
   "language": "julia"
  }
 },
 "nbformat": 4
}