{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# AtomsBase integration"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "[AtomsBase.jl](https://github.com/JuliaMolSim/AtomsBase.jl) is a common interface\n",
    "for representing atomic structures in Julia. DFTK directly supports using such\n",
    "structures to run a calculation as is demonstrated here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using DFTK"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Feeding an AtomsBase AbstractSystem to DFTK\n",
    "In this example we construct a silicon system using the `ase.build.bulk` routine\n",
    "from the [atomistic simulation environment](https://wiki.fysik.dtu.dk/ase/index.html)\n",
    "(ASE), which is exposed by [ASEconvert](https://github.com/mfherbst/ASEconvert.jl)\n",
    "as an AtomsBase `AbstractSystem`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "FlexibleSystem(Si₂, periodic = TTT):\n    bounding_box      : [       0    2.715    2.715;\n                            2.715        0    2.715;\n                            2.715    2.715        0]u\"Å\"\n\n    Atom(Si, [       0,        0,        0]u\"Å\")\n    Atom(Si, [  1.3575,   1.3575,   1.3575]u\"Å\")\n"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "# Construct bulk system and convert to an AbstractSystem\n",
    "using ASEconvert\n",
    "system_ase = ase.build.bulk(\"Si\")\n",
    "system = pyconvert(AbstractSystem, system_ase)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "To use an AbstractSystem in DFTK, we attach pseudopotentials, construct a DFT model,\n",
    "discretise and solve:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime\n",
      "---   ---------------   ---------   ---------   ----   ------\n",
      "  1   -7.921671305758                   -0.69    5.8         \n",
      "  2   -7.926164107424       -2.35       -1.22    1.0    250ms\n",
      "  3   -7.926837685413       -3.17       -2.37    1.6    325ms\n",
      "  4   -7.926861531114       -4.62       -3.04    2.9    298ms\n",
      "  5   -7.926861649114       -6.93       -3.43    1.6    250ms\n",
      "  6   -7.926861671688       -7.65       -3.84    1.8    288ms\n",
      "  7   -7.926861680394       -8.06       -4.32    1.2    240ms\n",
      "  8   -7.926861681757       -8.87       -4.87    2.0    258ms\n",
      "  9   -7.926861681862       -9.98       -5.28    2.0    265ms\n",
      " 10   -7.926861681871      -11.04       -5.79    1.8    270ms\n",
      " 11   -7.926861681873      -11.76       -6.43    1.8    258ms\n",
      " 12   -7.926861681873      -13.07       -7.40    2.0    268ms\n",
      " 13   -7.926861681873      -14.75       -8.18    3.1    311ms\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "system = attach_psp(system; Si=\"hgh/lda/si-q4\")\n",
    "\n",
    "model  = model_LDA(system; temperature=1e-3)\n",
    "basis  = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])\n",
    "scfres = self_consistent_field(basis, tol=1e-8);"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "If we did not want to use ASE we could of course use any other package\n",
    "which yields an AbstractSystem object. This includes:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Reading a system using AtomsIO"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime\n",
      "---   ---------------   ---------   ---------   ----   ------\n",
      "  1   -7.921719603968                   -0.69    5.6         \n",
      "  2   -7.926168395527       -2.35       -1.22    1.0    251ms\n",
      "  3   -7.926836444110       -3.18       -2.37    1.9    343ms\n",
      "  4   -7.926861501842       -4.60       -3.01    2.5    290ms\n",
      "  5   -7.926861632765       -6.88       -3.33    1.8    256ms\n",
      "  6   -7.926861666827       -7.47       -3.73    1.8    253ms\n",
      "  7   -7.926861680698       -7.86       -4.40    1.2    278ms\n",
      "  8   -7.926861681828       -8.95       -5.03    2.1    267ms\n",
      "  9   -7.926861681851      -10.64       -5.13    1.9    262ms\n",
      " 10   -7.926861681871      -10.70       -5.87    1.1    237ms\n",
      " 11   -7.926861681873      -11.87       -6.87    2.0    291ms\n",
      " 12   -7.926861681873      -13.28       -7.34    3.0    305ms\n",
      " 13   -7.926861681873   +  -14.75       -8.31    2.2    270ms\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "using AtomsIO\n",
    "\n",
    "# Read a file using [AtomsIO](https://github.com/mfherbst/AtomsIO.jl),\n",
    "# which directly yields an AbstractSystem.\n",
    "system = load_system(\"Si.extxyz\")\n",
    "\n",
    "# Now run the LDA calculation:\n",
    "system = attach_psp(system; Si=\"hgh/lda/si-q4\")\n",
    "model  = model_LDA(system; temperature=1e-3)\n",
    "basis  = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])\n",
    "scfres = self_consistent_field(basis, tol=1e-8);"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "The same could be achieved using [ExtXYZ](https://github.com/libAtoms/ExtXYZ.jl)\n",
    "by `system = Atoms(read_frame(\"Si.extxyz\"))`,\n",
    "since the `ExtXYZ.Atoms` object is directly AtomsBase-compatible."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Directly setting up a system in AtomsBase"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime\n",
      "---   ---------------   ---------   ---------   ----   ------\n",
      "  1   -7.921699595531                   -0.69    5.8         \n",
      "  2   -7.926169244709       -2.35       -1.22    1.0    256ms\n",
      "  3   -7.926839506380       -3.17       -2.37    1.6    288ms\n",
      "  4   -7.926864878291       -4.60       -3.00    2.8    344ms\n",
      "  5   -7.926865032982       -6.81       -3.29    1.9    284ms\n",
      "  6   -7.926865076868       -7.36       -3.70    1.6    296ms\n",
      "  7   -7.926865091899       -7.82       -4.42    1.4    253ms\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "using AtomsBase\n",
    "using Unitful\n",
    "using UnitfulAtomic\n",
    "\n",
    "# Construct a system in the AtomsBase world\n",
    "a = 10.26u\"bohr\"  # Silicon lattice constant\n",
    "lattice = a / 2 * [[0, 1, 1.],  # Lattice as vector of vectors\n",
    "                   [1, 0, 1.],\n",
    "                   [1, 1, 0.]]\n",
    "atoms  = [:Si => ones(3)/8, :Si => -ones(3)/8]\n",
    "system = periodic_system(atoms, lattice; fractional=true)\n",
    "\n",
    "# Now run the LDA calculation:\n",
    "system = attach_psp(system; Si=\"hgh/lda/si-q4\")\n",
    "model  = model_LDA(system; temperature=1e-3)\n",
    "basis  = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])\n",
    "scfres = self_consistent_field(basis, tol=1e-4);"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Obtaining an AbstractSystem from DFTK data"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "At any point we can also get back the DFTK model as an\n",
    "AtomsBase-compatible `AbstractSystem`:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "FlexibleSystem(Si₂, periodic = TTT):\n    bounding_box      : [       0     5.13     5.13;\n                             5.13        0     5.13;\n                             5.13     5.13        0]u\"a₀\"\n\n    Atom(Si, [  1.2825,   1.2825,   1.2825]u\"a₀\")\n    Atom(Si, [ -1.2825,  -1.2825,  -1.2825]u\"a₀\")\n"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "second_system = atomic_system(model)"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "Similarly DFTK offers a method to the `atomic_system` and `periodic_system` functions\n",
    "(from AtomsBase), which enable a seamless conversion of the usual data structures for\n",
    "setting up DFTK calculations into an `AbstractSystem`:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "FlexibleSystem(Si₂, periodic = TTT):\n    bounding_box      : [       0  5.13155  5.13155;\n                          5.13155        0  5.13155;\n                          5.13155  5.13155        0]u\"a₀\"\n\n    Atom(Si, [ 1.28289,  1.28289,  1.28289]u\"a₀\")\n    Atom(Si, [-1.28289, -1.28289, -1.28289]u\"a₀\")\n"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "lattice = 5.431u\"Å\" / 2 * [[0 1 1.];\n",
    "                           [1 0 1.];\n",
    "                           [1 1 0.]];\n",
    "Si = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\n",
    "atoms     = [Si, Si]\n",
    "positions = [ones(3)/8, -ones(3)/8]\n",
    "\n",
    "third_system = atomic_system(lattice, atoms, positions)"
   ],
   "metadata": {},
   "execution_count": 7
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.5",
   "language": "julia"
  }
 },
 "nbformat": 4
}