v0.2.5/examples/polarizability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Polarizability by linear response\n",
    "\n",
    "We compute the polarizability of a Helium atom. The polarizability\n",
    "is defined as the change in dipole moment\n",
    "$$\n",
    "\\mu = \\int r ρ(r) dr\n",
    "$$\n",
    "with respect to a small uniform electric field ``E = -x``.\n",
    "\n",
    "We compute this in two ways: first by finite differences (applying a\n",
    "finite electric field), then by linear response. Note that DFTK is\n",
    "not really adapted to isolated atoms because it uses periodic\n",
    "boundary conditions. Nevertheless we can simply embed the Helium\n",
    "atom in a large enough box (although this is computationally wasteful).\n",
    "\n",
    "As in other tests, this is not fully converged, convergence\n",
    "parameters were simply selected for fast execution on CI,"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using DFTK\n",
    "using LinearAlgebra\n",
    "\n",
    "a = 10.\n",
    "lattice = a * I(3)  # cube of ``a`` bohrs\n",
    "He = ElementPsp(:He, psp=load_psp(\"hgh/lda/He-q2\"))\n",
    "atoms = [He => [[1/2; 1/2; 1/2]]]  # Helium at the center of the box\n",
    "\n",
    "kgrid = [1, 1, 1]  # no kpoint sampling for an isolated system\n",
    "Ecut = 30\n",
    "tol = 1e-8\n",
    "\n",
    "# dipole moment of a given density (assuming the current geometry)\n",
    "function dipole(ρ)\n",
    "    basis = ρ.basis\n",
    "    rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]\n",
    "    d = sum(rr .* ρ.real) * basis.model.unit_cell_volume / prod(basis.fft_size)\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Polarizability by finite differences\n",
    "We first compute the polarizability by finite differences.\n",
    "First compute the dipole moment at rest:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            Eₙ-Eₙ₋₁     ρout-ρin   Diag\n",
      "---   ---------------   ---------   --------   ----\n",
      "  1   -2.769906494251         NaN   2.99e-01    9.0 \n",
      "  2   -2.771190105730   -1.28e-03   4.76e-02    1.0 \n",
      "  3   -2.771212590881   -2.25e-05   3.56e-03    1.0 \n",
      "  4   -2.771212899540   -3.09e-07   7.41e-04    1.0 \n",
      "  5   -2.771212981178   -8.16e-08   7.84e-05    2.0 \n",
      "  6   -2.771212981532   -3.54e-10   4.07e-05    1.0 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-0.00013578956544657336"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "model = model_LDA(lattice, atoms; symmetries=false)\n",
    "basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\n",
    "res = self_consistent_field(basis, tol=tol)\n",
    "μref = dipole(res.ρ)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then in a small uniform field:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     Energy            Eₙ-Eₙ₋₁     ρout-ρin   Diag\n",
      "---   ---------------   ---------   --------   ----\n",
      "  1   -2.770005185494         NaN   2.99e-01    9.0 \n",
      "  2   -2.771277911455   -1.27e-03   4.75e-02    1.0 \n",
      "  3   -2.771299883896   -2.20e-05   3.67e-03    1.0 \n",
      "  4   -2.771300293005   -4.09e-07   7.09e-04    1.0 \n",
      "  5   -2.771300362204   -6.92e-08   3.61e-05    2.0 \n",
      "  6   -2.771300362318   -1.14e-10   2.00e-05    1.0 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "0.01762258255369235"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "ε = .01\n",
    "model_ε = model_LDA(lattice, atoms; extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],\n",
    "                    symmetries=false)\n",
    "basis_ε = PlaneWaveBasis(model_ε, Ecut; kgrid=kgrid)\n",
    "res_ε = self_consistent_field(basis_ε, tol=tol)\n",
    "με = dipole(res_ε.ρ)"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reference dipole:  -0.00013578956544657336\n",
      "Displaced dipole:  0.01762258255369235\n",
      "Polarizability :   1.775837211913892\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "polarizability = (με - μref) / ε\n",
    "\n",
    "println(\"Reference dipole:  $μref\")\n",
    "println(\"Displaced dipole:  $με\")\n",
    "println(\"Polarizability :   $polarizability\")"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "The result on more converged grids is very close to published results.\n",
    "For example [DOI 10.1039/C8CP03569E](https://doi.org/10.1039/C8CP03569E)\n",
    "quotes **1.65** with LSDA and **1.38** with CCSD(T)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Polarizability by linear response\n",
    "Now we use linear response to compute this analytically; we refer to standard\n",
    "textbooks for the formalism. In the following, ``\\chi_0`` is the\n",
    "independent-particle polarizability, and ``K`` the\n",
    "Hartree-exchange-correlation kernel. We denote with ``dV_{\\rm ext}`` an external\n",
    "perturbing potential (like in this case the uniform electric field). Then:\n",
    "$$\n",
    "d\\rho = \\chi_0 dV = \\chi_0 (dV_{\\rm ext} + K d\\rho),\n",
    "$$\n",
    "which implies\n",
    "$$\n",
    "d\\rho = (1-\\chi_0 K)^-1 \\chi_0 dV_{\\rm ext}.\n",
    "$$\n",
    "From this we identify the polarizability operator to be ``\\chi = (1-\\chi_0 K)^-1 \\chi_0``.\n",
    "Numerically, we apply ``\\chi`` to ``dV = -x`` by solving a linear equation\n",
    "(the Dyson equation) iteratively."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: using KrylovKit.basis in module ##274 conflicts with an existing identifier.\n",
      "┌ Info: GMRES linsolve in iter 1; step 1: normres = 2.494245705328e-01\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:55\n",
      "┌ Info: GMRES linsolve in iter 1; step 2: normres = 3.770663921762e-03\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; step 3: normres = 2.985072405936e-04\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; step 4: normres = 5.448461702985e-06\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; step 5: normres = 4.738323704907e-07\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; step 6: normres = 1.877275072317e-09\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; step 7: normres = 1.262034512755e-11\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; step 8: normres = 3.944597968025e-13\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 1; finised at step 8: normres = 3.944597968025e-13\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:96\n",
      "┌ Info: GMRES linsolve in iter 2; step 1: normres = 8.873941545557e-07\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:55\n",
      "┌ Info: GMRES linsolve in iter 2; step 2: normres = 9.692220031648e-08\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 2; step 3: normres = 1.953654345662e-09\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 2; step 4: normres = 2.780611375355e-11\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 2; step 5: normres = 6.378716681182e-13\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 2; finised at step 5: normres = 6.378716681182e-13\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:96\n",
      "┌ Info: GMRES linsolve in iter 3; step 1: normres = 4.533952806588e-11\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:55\n",
      "┌ Info: GMRES linsolve in iter 3; step 2: normres = 4.940722609283e-12\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 3; step 3: normres = 4.053619397579e-14\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89\n",
      "┌ Info: GMRES linsolve in iter 3; finised at step 3: normres = 4.053619397579e-14\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:96\n",
      "┌ Info: GMRES linsolve converged at iteration 3, step 3:\n",
      "│ *  norm of residual = 4.4209930061973646e-14\n",
      "│ *  number of operations = 18\n",
      "└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:127\n",
      "Non-interacting polarizability: 1.9293430757915033\n",
      "Interacting polarizability:     1.7766945841079098\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "using KrylovKit\n",
    "\n",
    "# KrylovKit cannot deal with the density as a 3D array, so we need `vec` to vectorize\n",
    "# and `devec` to \"unvectorize\"\n",
    "devec(arr) = from_real(basis, reshape(arr, size(res.ρ.real)))\n",
    "\n",
    "# Apply (1- χ0 K) to a vectorized dρ\n",
    "function dielectric_operator(dρ)\n",
    "    dρ = devec(dρ)\n",
    "    dv = apply_kernel(basis, dρ; ρ=res.ρ)\n",
    "    χ0dv = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dv...)[1]\n",
    "    vec((dρ - χ0dv).real)\n",
    "end\n",
    "\n",
    "# dVext is the potential from a uniform field interacting with the dielectric dipole\n",
    "# of the density.\n",
    "dVext = from_real(basis, [-a * (r[1] - 1/2) for r in r_vectors(basis)])\n",
    "\n",
    "# Apply χ0 once to get non-interacting dipole\n",
    "dρ_nointeract = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dVext)[1]\n",
    "\n",
    "# Solve Dyson equation to get interacting dipole\n",
    "dρ = devec(linsolve(dielectric_operator, vec(dρ_nointeract.real), verbosity=3)[1])\n",
    "\n",
    "println(\"Non-interacting polarizability: $(dipole(dρ_nointeract))\")\n",
    "println(\"Interacting polarizability:     $(dipole(dρ))\")"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "As expected, the interacting polarizability matches the finite difference\n",
    "result. The non-interacting polarizability is higher."
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.3"
  },
  "kernelspec": {
   "name": "julia-1.5",
   "display_name": "Julia 1.5.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}