{
 "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\n",
      "---   ---------------   ---------   ---------   ----\n",
      "  1   -7.921682118187                   -0.69    5.9\n",
      "  2   -7.926161830555       -2.35       -1.22    1.0\n",
      "  3   -7.926837794412       -3.17       -2.37    1.8\n",
      "  4   -7.926861520074       -4.62       -3.04    2.9\n",
      "  5   -7.926861644079       -6.91       -3.41    1.9\n",
      "  6   -7.926861669404       -7.60       -3.78    1.6\n",
      "  7   -7.926861680431       -7.96       -4.33    1.2\n",
      "  8   -7.926861681752       -8.88       -4.81    2.0\n",
      "  9   -7.926861681859       -9.97       -5.22    1.5\n",
      " 10   -7.926861681871      -10.91       -5.90    1.6\n",
      " 11   -7.926861681873      -11.86       -7.07    2.1\n",
      " 12   -7.926861681873      -13.47       -7.42    3.5\n",
      " 13   -7.926861681873      -15.05       -8.29    1.4\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": [],
   "cell_type": "code",
   "source": [
    "using AtomsIO\n",
    "\n",
    "# Read a file using [AtomsIO](https://github.com/mfherbst/AtomsIO.jl),"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "which directly yields an AbstractSystem."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            log10(ΔE)   log10(Δρ)   Diag\n",
      "---   ---------------   ---------   ---------   ----\n",
      "  1   -7.921699031998                   -0.69    5.6\n",
      "  2   -7.926167595907       -2.35       -1.22    1.0\n",
      "  3   -7.926837017823       -3.17       -2.37    1.9\n",
      "  4   -7.926861520391       -4.61       -3.03    2.6\n",
      "  5   -7.926861642562       -6.91       -3.39    1.9\n",
      "  6   -7.926861669757       -7.57       -3.79    1.9\n",
      "  7   -7.926861680378       -7.97       -4.33    1.2\n",
      "  8   -7.926861681766       -8.86       -4.86    2.0\n",
      "  9   -7.926861681862      -10.02       -5.29    1.8\n",
      " 10   -7.926861681871      -11.01       -5.86    2.0\n",
      " 11   -7.926861681873      -11.92       -6.99    1.8\n",
      " 12   -7.926861681873      -13.48       -7.67    3.2\n",
      " 13   -7.926861681873      -15.05       -8.20    2.9\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "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": 5
  },
  {
   "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\n",
      "---   ---------------   ---------   ---------   ----\n",
      "  1   -7.921684822302                   -0.69    5.6\n",
      "  2   -7.926171465626       -2.35       -1.22    1.0\n",
      "  3   -7.926840765557       -3.17       -2.37    1.8\n",
      "  4   -7.926864901222       -4.62       -3.04    2.6\n",
      "  5   -7.926865055580       -6.81       -3.41    2.2\n",
      "  6   -7.926865081823       -7.58       -3.81    1.9\n",
      "  7   -7.926865091371       -8.02       -4.29    1.2\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": 6
  },
  {
   "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": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "second_system = atomic_system(model)"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "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": 8
    }
   ],
   "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": 8
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.4"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.4",
   "language": "julia"
  }
 },
 "nbformat": 4
}