{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Modelling a gallium arsenide surface\n",
    "\n",
    "This example shows how to use the\n",
    "[atomistic simulation environment](https://wiki.fysik.dtu.dk/ase/index.html),\n",
    "or ASE for short,\n",
    "to set up a particular gallium arsenide surface\n",
    "and run the resulting calculation in DFTK.\n",
    "The particular example we consider the (1, 1, 0) GaAs surface separated by vacuum."
   ],
   "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=\"https://docs.dftk.org/stable/surface.png\" width=500 height=500 />"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Use the `load_atoms`, `load_positions` 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",
    "positions = load_positions(surface)\n",
    "lattice   = load_lattice(surface)\n",
    "atoms = map(load_atoms(surface)) do el\n",
    "    if el.symbol == :Ga\n",
    "        ElementPsp(:Ga, psp=load_psp(\"hgh/pbe/ga-q3.hgh\"))\n",
    "    elseif el.symbol == :As\n",
    "        ElementPsp(:As, psp=load_psp(\"hgh/pbe/as-q5.hgh\"))\n",
    "    else\n",
    "        error(\"Unsupported element: $el\")\n",
    "    end\n",
    "end;"
   ],
   "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` or try `mixing=KerkerMixing()`\n",
    "to see the full challenge of this system."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            log10(ΔE)   log10(Δρ)   Diag\n",
      "---   ---------------   ---------   ---------   ----\n",
      "  1   -16.58751649151                   -0.58    4.7\n",
      "  2   -16.72549999485       -0.86       -1.01    1.1\n",
      "  3   -16.73067512899       -2.29       -1.58    3.7\n",
      "  4   -16.73128215007       -3.22       -2.16    2.9\n",
      "  5   -16.73133383662       -4.29       -2.61    3.0\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "model = model_DFT(lattice, atoms, positions, [:gga_x_pbe, :gga_c_pbe],\n",
    "                  temperature=0.001, smearing=DFTK.Smearing.Gaussian())\n",
    "basis = PlaneWaveBasis(model; Ecut, kgrid)\n",
    "\n",
    "scfres = self_consistent_field(basis, tol=1e-4, mixing=LdosMixing());"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Energy breakdown (in Ha):\n    Kinetic             5.8605554 \n    AtomicLocal         -105.6249290\n    AtomicNonlocal      2.3500153 \n    Ewald               35.5044300\n    PspCorrection       0.2016043 \n    Hartree             49.5751872\n    Xc                  -4.5981931\n    Entropy             -0.0000040\n\n    total               -16.731333836616"
     },
     "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.8.1"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.1",
   "language": "julia"
  }
 },
 "nbformat": 4
}