# Polarizability by linear response¶

We compute the polarizability of a Helium atom. The polarizability is defined as the change in dipole moment $$\mu = \int r ρ(r) dr$$ with respect to a small uniform electric field E = -x.

We compute this in two ways: first by finite differences (applying a finite electric field), then by linear response. Note that DFTK is not really adapted to isolated atoms because it uses periodic boundary conditions. Nevertheless we can simply embed the Helium atom in a large enough box (although this is computationally wasteful).

As in other tests, this is not fully converged, convergence parameters were simply selected for fast execution on CI,

In [1]:
using DFTK
using LinearAlgebra

a = 10.
lattice = a * I(3)  # cube of a bohrs
He = ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))
atoms = [He => [[1/2; 1/2; 1/2]]]  # Helium at the center of the box

kgrid = [1, 1, 1]  # no k-point sampling for an isolated system
Ecut = 30
tol = 1e-8

# dipole moment of a given density (assuming the current geometry)
function dipole(basis, ρ)
rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
sum(rr .* ρ) * basis.dvol
end;


## Using finite differences¶

We first compute the polarizability by finite differences. First compute the dipole moment at rest:

In [2]:
model = model_LDA(lattice, atoms; symmetries=false)
basis = PlaneWaveBasis(model; Ecut, kgrid)
res   = self_consistent_field(basis, tol=tol)
μref  = dipole(basis, res.ρ)

n     Energy            Eₙ-Eₙ₋₁     ρout-ρin   α      Diag
---   ---------------   ---------   --------   ----   ----
1   -2.769623389309         NaN   2.95e-01   0.80    8.0
2   -2.771174788933   -1.55e-03   5.19e-02   0.80    1.0
3   -2.771212876710   -3.81e-05   1.07e-03   0.80    2.0
4   -2.771212980625   -1.04e-07   1.78e-04   0.80    2.0
5   -2.771212981724   -1.10e-09   3.48e-05   0.80    2.0

Out[2]:
-0.0001361390138976282

Then in a small uniform field:

In [3]:
ε = .01
model_ε = model_LDA(lattice, atoms; extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
symmetries=false)
basis_ε = PlaneWaveBasis(model_ε; Ecut, kgrid)
res_ε = self_consistent_field(basis_ε, tol=tol)
με = dipole(basis_ε, res_ε.ρ)

n     Energy            Eₙ-Eₙ₋₁     ρout-ρin   α      Diag
---   ---------------   ---------   --------   ----   ----
1   -2.770002807670         NaN   2.99e-01   0.80    9.0
2   -2.771276930867   -1.27e-03   4.75e-02   0.80    1.0
3   -2.771299473233   -2.25e-05   4.39e-03   0.80    2.0
4   -2.771300361930   -8.89e-07   9.54e-05   0.80    2.0
5   -2.771300362350   -4.20e-10   1.03e-05   0.80    2.0

Out[3]:
0.01761830703783743
In [4]:
polarizability = (με - μref) / ε

println("Reference dipole:  $μref") println("Displaced dipole:$με")
println("Polarizability :   $polarizability")  Reference dipole: -0.0001361390138976282 Displaced dipole: 0.01761830703783743 Polarizability : 1.7754446051735058  The result on more converged grids is very close to published results. For example DOI 10.1039/C8CP03569E quotes 1.65 with LSDA and 1.38 with CCSD(T). ## Using linear response¶ Now we use linear response to compute this analytically; we refer to standard textbooks for the formalism. In the following, \chi_0 is the independent-particle polarizability, and K the Hartree-exchange-correlation kernel. We denote with \delta V_{\rm ext} an external perturbing potential (like in this case the uniform electric field). Then: $$\delta\rho = \chi_0 \delta V = \chi_0 (\delta V_{\rm ext} + K \delta\rho),$$ which implies $$\delta\rho = (1-\chi_0 K)^-1 \chi_0 \delta V_{\rm ext}.$$ From this we identify the polarizability operator to be \chi = (1-\chi_0 K)^{-1} \chi_0. Numerically, we apply \chi to \delta V = -x by solving a linear equation (the Dyson equation) iteratively. In [5]: using KrylovKit # Apply (1- χ0 K) function dielectric_operator(δρ) δV = apply_kernel(basis, δρ; ρ=res.ρ) χ0δV = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, δV) δρ - χ0δV end # δVext is the potential from a uniform field interacting with the dielectric dipole # of the density. δVext = [-a * (r[1] - 1/2) for r in r_vectors(basis)] δVext = cat(δVext; dims=4) # Apply χ0 once to get non-interacting dipole δρ_nointeract = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, δVext) # Solve Dyson equation to get interacting dipole δρ = linsolve(dielectric_operator, δρ_nointeract, verbosity=3)[1] println("Non-interacting polarizability:$(dipole(basis, δρ_nointeract))")
println("Interacting polarizability:     \$(dipole(basis, δρ))")

WARNING: using KrylovKit.basis in module ##299 conflicts with an existing identifier.
┌ Info: GMRES linsolve in iter 1; step 1: normres = 2.491879897268e-01
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:55
┌ Info: GMRES linsolve in iter 1; step 2: normres = 3.716159829337e-03
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; step 3: normres = 1.189406142900e-04
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; step 4: normres = 6.519644366479e-06
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; step 5: normres = 6.321696599733e-07
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; step 6: normres = 1.868577004875e-09
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; step 7: normres = 1.415566786213e-11
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; step 8: normres = 1.216373501087e-12
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 1; finished at step 8: normres = 1.216373501087e-12
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:96
┌ Info: GMRES linsolve in iter 2; step 1: normres = 9.623760450371e-07
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:55
┌ Info: GMRES linsolve in iter 2; step 2: normres = 3.872342023505e-08
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 2; step 3: normres = 2.242307819033e-09
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 2; step 4: normres = 1.947660226669e-11
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 2; step 5: normres = 6.561759879075e-13
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 2; finished at step 5: normres = 6.561759879075e-13
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:96
┌ Info: GMRES linsolve in iter 3; step 1: normres = 1.869676299332e-11
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:55
┌ Info: GMRES linsolve in iter 3; step 2: normres = 1.075829261256e-12
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:89
┌ Info: GMRES linsolve in iter 3; finished at step 2: normres = 1.075829261256e-12
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:96
┌ Info: GMRES linsolve converged at iteration 3, step 2:
│ *  norm of residual = 1.071364497536787e-12
│ *  number of operations = 17
└ @ KrylovKit /home/runner/.julia/packages/KrylovKit/YPiz7/src/linsolve/gmres.jl:127
Non-interacting polarizability: 1.9277533317790276
Interacting polarizability:     1.775160610781069


As expected, the interacting polarizability matches the finite difference result. The non-interacting polarizability is higher.