Polarizability by linear response

We compute the polarizability of a Helium atom. The polarizability is defined as the change in dipole moment

\[μ = ∫ 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,

using DFTK
using LinearAlgebra
using PseudoPotentialData

a = 10.
lattice = a * I(3)  # cube of ``a`` bohrs
pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.largecore.gth")
# Helium at the center of the box
atoms     = [ElementPsp(:He, pseudopotentials)]
positions = [[1/2, 1/2, 1/2]]


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 = [(r[1] - a/2) for r in r_vectors_cart(basis)]
    sum(rr .* ρ) * basis.dvol
end;

Using finite differences

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

model  = model_DFT(lattice, atoms, positions;
                   functionals=LDA(), symmetries=false)
basis  = PlaneWaveBasis(model; Ecut, kgrid)
scfres = self_consistent_field(basis; tol)
μref   = dipole(basis, scfres.ρ)

-0.0001345744216435374

Then in a small uniform field:

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

0.01761222116765141

polarizability = (με - μref) / ε

println("Reference dipole:  $μref")
println("Displaced dipole:  $με")
println("Polarizability :   $polarizability")

Reference dipole:  -0.0001345744216435374
Displaced dipole:  0.01761222116765141
Polarizability :   1.7746795589294948

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 (also known as density-functional perturbation theory) to compute this analytically; we refer to standard textbooks for the formalism. In the following, $χ_0$ is the independent-particle polarizability, and $K$ the Hartree-exchange-correlation kernel. We denote with $δV_{\rm ext}$ an external perturbing potential (like in this case the uniform electric field).

# `δVext` is the potential from a uniform field interacting with the dielectric dipole
# of the density.
δVext = [-(r[1] - a/2) for r in r_vectors_cart(basis)]
δVext = cat(δVext; dims=4)

54×54×54×1 Array{Float64, 4}:
[:, :, 1, 1] =
  5.0       5.0       5.0       5.0      …   5.0       5.0       5.0
  4.81481   4.81481   4.81481   4.81481      4.81481   4.81481   4.81481
  4.62963   4.62963   4.62963   4.62963      4.62963   4.62963   4.62963
  4.44444   4.44444   4.44444   4.44444      4.44444   4.44444   4.44444
  4.25926   4.25926   4.25926   4.25926      4.25926   4.25926   4.25926
  4.07407   4.07407   4.07407   4.07407  …   4.07407   4.07407   4.07407
  3.88889   3.88889   3.88889   3.88889      3.88889   3.88889   3.88889
  3.7037    3.7037    3.7037    3.7037       3.7037    3.7037    3.7037
  3.51852   3.51852   3.51852   3.51852      3.51852   3.51852   3.51852
  3.33333   3.33333   3.33333   3.33333      3.33333   3.33333   3.33333
  ⋮                                      ⋱                      
 -3.33333  -3.33333  -3.33333  -3.33333  …  -3.33333  -3.33333  -3.33333
 -3.51852  -3.51852  -3.51852  -3.51852     -3.51852  -3.51852  -3.51852
 -3.7037   -3.7037   -3.7037   -3.7037      -3.7037   -3.7037   -3.7037
 -3.88889  -3.88889  -3.88889  -3.88889     -3.88889  -3.88889  -3.88889
 -4.07407  -4.07407  -4.07407  -4.07407     -4.07407  -4.07407  -4.07407
 -4.25926  -4.25926  -4.25926  -4.25926  …  -4.25926  -4.25926  -4.25926
 -4.44444  -4.44444  -4.44444  -4.44444     -4.44444  -4.44444  -4.44444
 -4.62963  -4.62963  -4.62963  -4.62963     -4.62963  -4.62963  -4.62963
 -4.81481  -4.81481  -4.81481  -4.81481     -4.81481  -4.81481  -4.81481

[:, :, 2, 1] =
  5.0       5.0       5.0       5.0      …   5.0       5.0       5.0
  4.81481   4.81481   4.81481   4.81481      4.81481   4.81481   4.81481
  4.62963   4.62963   4.62963   4.62963      4.62963   4.62963   4.62963
  4.44444   4.44444   4.44444   4.44444      4.44444   4.44444   4.44444
  4.25926   4.25926   4.25926   4.25926      4.25926   4.25926   4.25926
  4.07407   4.07407   4.07407   4.07407  …   4.07407   4.07407   4.07407
  3.88889   3.88889   3.88889   3.88889      3.88889   3.88889   3.88889
  3.7037    3.7037    3.7037    3.7037       3.7037    3.7037    3.7037
  3.51852   3.51852   3.51852   3.51852      3.51852   3.51852   3.51852
  3.33333   3.33333   3.33333   3.33333      3.33333   3.33333   3.33333
  ⋮                                      ⋱                      
 -3.33333  -3.33333  -3.33333  -3.33333  …  -3.33333  -3.33333  -3.33333
 -3.51852  -3.51852  -3.51852  -3.51852     -3.51852  -3.51852  -3.51852
 -3.7037   -3.7037   -3.7037   -3.7037      -3.7037   -3.7037   -3.7037
 -3.88889  -3.88889  -3.88889  -3.88889     -3.88889  -3.88889  -3.88889
 -4.07407  -4.07407  -4.07407  -4.07407     -4.07407  -4.07407  -4.07407
 -4.25926  -4.25926  -4.25926  -4.25926  …  -4.25926  -4.25926  -4.25926
 -4.44444  -4.44444  -4.44444  -4.44444     -4.44444  -4.44444  -4.44444
 -4.62963  -4.62963  -4.62963  -4.62963     -4.62963  -4.62963  -4.62963
 -4.81481  -4.81481  -4.81481  -4.81481     -4.81481  -4.81481  -4.81481

[:, :, 3, 1] =
  5.0       5.0       5.0       5.0      …   5.0       5.0       5.0
  4.81481   4.81481   4.81481   4.81481      4.81481   4.81481   4.81481
  4.62963   4.62963   4.62963   4.62963      4.62963   4.62963   4.62963
  4.44444   4.44444   4.44444   4.44444      4.44444   4.44444   4.44444
  4.25926   4.25926   4.25926   4.25926      4.25926   4.25926   4.25926
  4.07407   4.07407   4.07407   4.07407  …   4.07407   4.07407   4.07407
  3.88889   3.88889   3.88889   3.88889      3.88889   3.88889   3.88889
  3.7037    3.7037    3.7037    3.7037       3.7037    3.7037    3.7037
  3.51852   3.51852   3.51852   3.51852      3.51852   3.51852   3.51852
  3.33333   3.33333   3.33333   3.33333      3.33333   3.33333   3.33333
  ⋮                                      ⋱                      
 -3.33333  -3.33333  -3.33333  -3.33333  …  -3.33333  -3.33333  -3.33333
 -3.51852  -3.51852  -3.51852  -3.51852     -3.51852  -3.51852  -3.51852
 -3.7037   -3.7037   -3.7037   -3.7037      -3.7037   -3.7037   -3.7037
 -3.88889  -3.88889  -3.88889  -3.88889     -3.88889  -3.88889  -3.88889
 -4.07407  -4.07407  -4.07407  -4.07407     -4.07407  -4.07407  -4.07407
 -4.25926  -4.25926  -4.25926  -4.25926  …  -4.25926  -4.25926  -4.25926
 -4.44444  -4.44444  -4.44444  -4.44444     -4.44444  -4.44444  -4.44444
 -4.62963  -4.62963  -4.62963  -4.62963     -4.62963  -4.62963  -4.62963
 -4.81481  -4.81481  -4.81481  -4.81481     -4.81481  -4.81481  -4.81481

;;; … 

[:, :, 52, 1] =
  5.0       5.0       5.0       5.0      …   5.0       5.0       5.0
  4.81481   4.81481   4.81481   4.81481      4.81481   4.81481   4.81481
  4.62963   4.62963   4.62963   4.62963      4.62963   4.62963   4.62963
  4.44444   4.44444   4.44444   4.44444      4.44444   4.44444   4.44444
  4.25926   4.25926   4.25926   4.25926      4.25926   4.25926   4.25926
  4.07407   4.07407   4.07407   4.07407  …   4.07407   4.07407   4.07407
  3.88889   3.88889   3.88889   3.88889      3.88889   3.88889   3.88889
  3.7037    3.7037    3.7037    3.7037       3.7037    3.7037    3.7037
  3.51852   3.51852   3.51852   3.51852      3.51852   3.51852   3.51852
  3.33333   3.33333   3.33333   3.33333      3.33333   3.33333   3.33333
  ⋮                                      ⋱                      
 -3.33333  -3.33333  -3.33333  -3.33333  …  -3.33333  -3.33333  -3.33333
 -3.51852  -3.51852  -3.51852  -3.51852     -3.51852  -3.51852  -3.51852
 -3.7037   -3.7037   -3.7037   -3.7037      -3.7037   -3.7037   -3.7037
 -3.88889  -3.88889  -3.88889  -3.88889     -3.88889  -3.88889  -3.88889
 -4.07407  -4.07407  -4.07407  -4.07407     -4.07407  -4.07407  -4.07407
 -4.25926  -4.25926  -4.25926  -4.25926  …  -4.25926  -4.25926  -4.25926
 -4.44444  -4.44444  -4.44444  -4.44444     -4.44444  -4.44444  -4.44444
 -4.62963  -4.62963  -4.62963  -4.62963     -4.62963  -4.62963  -4.62963
 -4.81481  -4.81481  -4.81481  -4.81481     -4.81481  -4.81481  -4.81481

[:, :, 53, 1] =
  5.0       5.0       5.0       5.0      …   5.0       5.0       5.0
  4.81481   4.81481   4.81481   4.81481      4.81481   4.81481   4.81481
  4.62963   4.62963   4.62963   4.62963      4.62963   4.62963   4.62963
  4.44444   4.44444   4.44444   4.44444      4.44444   4.44444   4.44444
  4.25926   4.25926   4.25926   4.25926      4.25926   4.25926   4.25926
  4.07407   4.07407   4.07407   4.07407  …   4.07407   4.07407   4.07407
  3.88889   3.88889   3.88889   3.88889      3.88889   3.88889   3.88889
  3.7037    3.7037    3.7037    3.7037       3.7037    3.7037    3.7037
  3.51852   3.51852   3.51852   3.51852      3.51852   3.51852   3.51852
  3.33333   3.33333   3.33333   3.33333      3.33333   3.33333   3.33333
  ⋮                                      ⋱                      
 -3.33333  -3.33333  -3.33333  -3.33333  …  -3.33333  -3.33333  -3.33333
 -3.51852  -3.51852  -3.51852  -3.51852     -3.51852  -3.51852  -3.51852
 -3.7037   -3.7037   -3.7037   -3.7037      -3.7037   -3.7037   -3.7037
 -3.88889  -3.88889  -3.88889  -3.88889     -3.88889  -3.88889  -3.88889
 -4.07407  -4.07407  -4.07407  -4.07407     -4.07407  -4.07407  -4.07407
 -4.25926  -4.25926  -4.25926  -4.25926  …  -4.25926  -4.25926  -4.25926
 -4.44444  -4.44444  -4.44444  -4.44444     -4.44444  -4.44444  -4.44444
 -4.62963  -4.62963  -4.62963  -4.62963     -4.62963  -4.62963  -4.62963
 -4.81481  -4.81481  -4.81481  -4.81481     -4.81481  -4.81481  -4.81481

[:, :, 54, 1] =
  5.0       5.0       5.0       5.0      …   5.0       5.0       5.0
  4.81481   4.81481   4.81481   4.81481      4.81481   4.81481   4.81481
  4.62963   4.62963   4.62963   4.62963      4.62963   4.62963   4.62963
  4.44444   4.44444   4.44444   4.44444      4.44444   4.44444   4.44444
  4.25926   4.25926   4.25926   4.25926      4.25926   4.25926   4.25926
  4.07407   4.07407   4.07407   4.07407  …   4.07407   4.07407   4.07407
  3.88889   3.88889   3.88889   3.88889      3.88889   3.88889   3.88889
  3.7037    3.7037    3.7037    3.7037       3.7037    3.7037    3.7037
  3.51852   3.51852   3.51852   3.51852      3.51852   3.51852   3.51852
  3.33333   3.33333   3.33333   3.33333      3.33333   3.33333   3.33333
  ⋮                                      ⋱                      
 -3.33333  -3.33333  -3.33333  -3.33333  …  -3.33333  -3.33333  -3.33333
 -3.51852  -3.51852  -3.51852  -3.51852     -3.51852  -3.51852  -3.51852
 -3.7037   -3.7037   -3.7037   -3.7037      -3.7037   -3.7037   -3.7037
 -3.88889  -3.88889  -3.88889  -3.88889     -3.88889  -3.88889  -3.88889
 -4.07407  -4.07407  -4.07407  -4.07407     -4.07407  -4.07407  -4.07407
 -4.25926  -4.25926  -4.25926  -4.25926  …  -4.25926  -4.25926  -4.25926
 -4.44444  -4.44444  -4.44444  -4.44444     -4.44444  -4.44444  -4.44444
 -4.62963  -4.62963  -4.62963  -4.62963     -4.62963  -4.62963  -4.62963
 -4.81481  -4.81481  -4.81481  -4.81481     -4.81481  -4.81481  -4.81481

Then:

\[δρ = χ_0 δV = χ_0 (δV_{\rm ext} + K δρ),\]

which implies

\[δρ = (1-χ_0 K)^{-1} χ_0 δV_{\rm ext}.\]

From this we identify the polarizability operator to be $χ = (1-χ_0 K)^{-1} χ_0$. Numerically, we apply $χ$ to $δV = -x$ by solving a linear equation (the Dyson equation) iteratively.

First we do this using the DFTK.solve_ΩplusK_split function provided by DFTK, which uses an adaptive Krylov subspace algorithm ^[HS2025]:

# Multiply δVext times the Bloch waves, then solve the Dyson equation:
δVψ = DFTK.multiply_ψ_by_blochwave(scfres.basis, scfres.ψ, δVext)
res = DFTK.solve_ΩplusK_split(scfres, -δVψ; verbose=true)

(δψ = Matrix{ComplexF64}[[0.0018322683424264032 + 0.002739654816142262im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.2586156838226855 - 0.16983557541115998im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.03422644949034864 - 0.022691341767800424im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.06710562756663412 + 0.045252918514107676im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206941494755575e-7 -2.502015100364878e-7 … -4.880642547335821e-8 -2.502014760275074e-7; -6.293616494036539e-7 -4.793036461964749e-7 … -1.207784505577762e-7 -4.793036192506421e-7; … ; 1.3626262382449146e-7 1.1356231791582643e-7 … 4.711994058732591e-8 1.1356232246040363e-7; 1.3305476665616598e-7 1.3632640121567523e-7 … 5.200687546243365e-8 1.3632642570350452e-7;;; -2.5020152043774494e-7 -1.7404909878582718e-7 … -3.5564945840766545e-8 -1.7404907438358128e-7; -4.793036534369e-7 -3.6797322633884953e-7 … -1.1202896796630411e-7 -3.6797320767418493e-7; … ; 1.1356231703647192e-7 1.1211100534774695e-7 … 1.0980687688982592e-7 1.1211101055523718e-7; 1.363263937456586e-7 1.2550864970698067e-7 … 5.0949105280469304e-8 1.2550866822336934e-7;;; -4.880644494406551e-8 -3.556495842840299e-8 … -1.1412157440064507e-8 -3.55649528376492e-8; -1.2077846731171146e-7 -1.1202898035037433e-7 … -9.786228063606738e-8 -1.1202897521179736e-7; … ; 4.7119934821255035e-8 1.098068637399198e-7 … 2.4938234811612797e-7 1.098068691897136e-7; 5.200685935547242e-8 5.0949091741093704e-8 … 4.910629537743492e-8 5.094909723947654e-8;;; … ;;; 6.475756811562236e-8 4.110331136415303e-8 … -6.68169822472093e-9 4.110330481209636e-8; 1.6897553478168972e-7 6.653282534195096e-8 … -1.329445099809433e-7 6.653282421845394e-8; … ; 3.267767554360692e-8 1.6017373249386057e-7 … 3.847037456755863e-7 1.6017373729944244e-7; -9.141993908997658e-8 -3.583034407915674e-8 … 7.037288704491634e-8 -3.583034715162134e-8;;; -4.880641945180553e-8 -3.5564939856647484e-8 … -1.1412151157112735e-8 -3.5564933652918363e-8; -1.207784474772575e-7 -1.1202896494329537e-7 … -9.786226608117877e-8 -1.1202895624601524e-7; … ; 4.711993383864382e-8 1.0980687152909891e-7 … 2.493823698061366e-7 1.0980687501002867e-7; 5.200687585034334e-8 5.094910712842332e-8 … 4.9106307763931384e-8 5.094911179690007e-8;;; -2.5020146484067366e-7 -1.7404905873677742e-7 … -3.5564936446588474e-8 -1.7404903350222938e-7; -4.793036116625335e-7 -3.6797319746088506e-7 … -1.1202895836040225e-7 -3.679731754101621e-7; … ; 1.1356232106663377e-7 1.1211100933313244e-7 … 1.0980687862822358e-7 1.1211101307213645e-7; 1.3632643233012773e-7 1.2550867841766587e-7 … 5.0949112505618316e-8 1.2550869621233512e-7;;;;], δHψ = Matrix{ComplexF64}[[-0.0010658077747930882 - 0.0015936227759341592im 0.05680841936906145 - 0.07651355517484931im 0.3285984347352903 - 1.0029011841249287im -0.11290388875255636 + 1.223789022704744im; -0.16583005730248968 + 0.10873812984187507im -1.2473127107153574 - 1.0444549588840324im 0.07237202727986157 - 0.22194640184102601im -0.024447652954613503 + 0.27112702885964063im; … ; -0.04092980656338081 + 0.026975990763637808im 0.016639908082481017 + 0.012289795102406305im -0.011709368897524625 + 0.005533182414484905im -0.00024169411419535476 - 0.01994027022893233im; 0.06000817839666706 - 0.040460910267267315im -0.02051272010728462 - 0.015301315863239694im -0.04117239539669053 - 0.12855445493128942im -0.015292899410304754 - 0.02115504582929931im]], δVind = [0.0019654317757027027 0.002081000590163893 … 0.0031536637771642424 0.0020810003089786182; 0.0038891307238482645 0.004779269357888311 … 0.01027666902201051 0.004779269061077429; … ; -0.006150847087537447 -0.008269601785407216 … -0.01730296287449031 -0.008269602366363642; 0.00023868214478133345 -0.00039869733827023164 … -0.003664695171734414 -0.000398697656210091;;; 0.0020810006817958827 0.0022232285116161165 … 0.003420118683420018 0.0022232281956799265; 0.004779269470115776 0.0057349537310421314 … 0.010903595690099934 0.005734953394077389; … ; -0.008269601588888356 -0.009601005706784058 … -0.01262310939850391 -0.009601006269369085; -0.0003986972400368636 -0.001073012669242569 … -0.004230680472820717 -0.0010730130224955818;;; 0.0031536650890714155 0.0034201199651426547 … 0.004832503079041746 0.0034201194256240964; 0.01027667038744626 0.010903596944093367 … 0.012272139241390772 0.010903596412070357; … ; -0.017302960864815538 -0.012623108167338207 … -0.00962026437011413 -0.012623108716294535; -0.0036646937756823657 -0.004230679223553521 … -0.005491107438645727 -0.004230679776199325;;; … ;;; -0.0023365591009449683 -0.002815461940384876 … 0.005693311020742097 -0.0028154614659701895; -0.01291309221244622 -0.017841182242332405 … 0.0114983364156872 -0.017841181445122468; … ; -0.025174453782722643 -0.012810886299816662 … -0.008395876255866255 -0.012810886830733429; 0.0075979807616899766 0.010225416726363835 … -0.005449735102368184 0.010225417350528336;;; 0.003153663472088633 0.003420118154104357 … 0.004832500446677541 0.003420117644804927; 0.010276668593529636 0.010903594991916739 … 0.012272137275608077 0.010903594604677137; … ; -0.017302964560265734 -0.012623110404538547 … -0.009620265915934232 -0.012623111022727562; -0.0036646956217347764 -0.004230681118249593 … -0.005491109432441571 -0.004230681685938001;;; 0.0020810002179709767 0.0022232279894244702 … 0.0034201177975493537 0.0022232276780614134; 0.004779268945385431 0.005734953132920266 … 0.010903594796346377 0.005734952818320164; … ; -0.008269602632843723 -0.009601006776716447 … -0.01262311059730754 -0.009601007387462099; -0.00039869775644811546 -0.0010730132461143531 … -0.004230681415228084 -0.0010730136002580811;;;;], δρ0 = [-3.6021403110694096e-7 -2.4888428028886556e-7 … -4.852837723284552e-8 -2.4888427978547155e-7; -5.972250799946921e-7 -4.606978717902057e-7 … -1.1905771619856915e-7 -4.6069786934454276e-7; … ; 1.273579329397026e-7 1.0836343256598744e-7 … 4.6407226784202206e-8 1.0836343154593099e-7; 1.035201155738034e-7 1.1957965320574523e-7 … 5.06950848965663e-8 1.1957965321126262e-7;;; -2.488842791933971e-7 -1.7310673213337796e-7 … -3.535698690584798e-8 -1.731067321580656e-7; -4.606978706411531e-7 -3.5721711759742175e-7 … -1.1091443212392905e-7 -3.57217116969696e-7; … ; 1.0836343357938724e-7 1.084650156915527e-7 … 1.0862242641141044e-7 1.0846501579662058e-7; 1.1957965411526324e-7 1.1611150634517688e-7 … 5.0208331298638934e-8 1.1611150669904393e-7;;; -4.8528376439939334e-8 -3.535698633124097e-8 … -1.1341325154235902e-8 -3.5356986327751053e-8; -1.190577178269587e-7 -1.1091443258556207e-7 … -9.770672994079223e-8 -1.109144325014961e-7; … ; 4.640722611163481e-8 1.0862242324367063e-7 … 2.484525327480597e-7 1.0862242444364282e-7; 5.069508394509679e-8 5.020832974520521e-8 … 4.9268613891654095e-8 5.020833015228259e-8;;; … ;;; 6.434570049507968e-8 4.0838617837319113e-8 … -6.637938583149967e-9 4.0838618409404994e-8; 1.6975033835648028e-7 6.683082288213937e-8 … -1.3340485968815428e-7 6.683082455854992e-8; … ; 3.271002713028142e-8 1.6023529500167108e-7 … 3.8390153332682707e-7 1.6023529434166055e-7; -9.27136714911937e-8 -3.633620060682858e-8 … 7.126312080237842e-8 -3.6336200901986036e-8;;; -4.852837657470272e-8 -3.5356986867997396e-8 … -1.1341325288679459e-8 -3.535698626019815e-8; -1.1905771709663468e-7 -1.1091443446941418e-7 … -9.770672981005825e-8 -1.1091443090152654e-7; … ; 4.640721883597459e-8 1.0862241671008711e-7 … 2.4845252412145857e-7 1.0862241597660523e-7; 5.0695081396904936e-8 5.020832781875643e-8 … 4.9268611882989245e-8 5.020832742323461e-8;;; -2.4888428037405844e-7 -1.731067335541244e-7 … -3.5356987034661895e-8 -1.731067327628604e-7; -4.6069787032684073e-7 -3.5721711934163697e-7 … -1.1091443118904834e-7 -3.5721711552633096e-7; … ; 1.083634283500115e-7 1.0846501118564406e-7 … 1.0862242102801896e-7 1.0846500988441313e-7; 1.1957965140507574e-7 1.1611150438327438e-7 … 5.020832990938703e-8 1.161115041077018e-7;;;;], δeigenvalues = [[0.0004950057641342383, 0.09679856885873216, 0.04049242819693278, 0.036904408382939485]], δoccupation = [[0.0, 0.0, 0.0, 0.0]], δεF = 0.0, ε = DFTK.DielectricAdjoint{Array{Float64, 4}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}, Float64, Vector{Vector{Float64}}, StaticArraysCore.SVector{3, Float64}}(Hamiltonian(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), HamiltonianBlock[DFTK.DftHamiltonianBlock(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Any[DFTK.FourierMultiplication{Float64, Vector{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.0, 0.19739208802178715, 0.7895683520871486, 1.7765287921960844, 3.1582734083485944, 4.934802200544679, 7.106115168784338, 9.67221231306757, 12.633093633394378, 15.988759129764759  …  20.133992978222288, 16.38354330580833, 13.027877809437953, 10.066996489111146, 7.5008993448279115, 5.329586376588253, 3.5530575843921683, 2.1713129682396586, 1.184352528130723, 0.5921762640653614]), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.1603906580688546 0.1602471565738079 … 0.15981840036267833 0.16024715657380792; 0.16024715657380803 0.16010321540768782 … 0.15967315066597979 0.16010321540768782; … ; 0.15981840036267841 0.1596731506659797 … 0.15923916395256382 0.1596731506659797; 0.16024715657380809 0.16010321540768785 … 0.1596731506659798 0.16010321540768785;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.1599588345026577 … 0.15952744714353234 0.1599588345026577; … ; 0.15967315066597973 0.15952744714353229 … 0.15909210945052638 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772;;; 0.1598184003626783 0.15967315066597965 … 0.15923916395256382 0.15967315066597965; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256385 0.15909210945052635 … 0.15865272241895537 0.15909210945052635; 0.1596731506659798 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; … ;;; 0.15910963392556968 0.15896221057980836 … 0.1585217272926462 0.15896221057980836; 0.15896221057980842 0.15881431920357322 … 0.15837242838914242 0.15881431920357322; … ; 0.15852172729264621 0.15837242838914237 … 0.157926332522532 0.15837242838914237; 0.15896221057980844 0.15881431920357328 … 0.1583724283891425 0.15881431920357328;;; 0.1598184003626783 0.15967315066597962 … 0.15923916395256382 0.15967315066597962; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256382 0.15909210945052635 … 0.15865272241895534 0.15909210945052635; 0.15967315066597979 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.15995883450265772 … 0.15952744714353237 0.15995883450265772; … ; 0.1596731506659797 0.15952744714353229 … 0.15909210945052635 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772]), DFTK.NonlocalOperator{Float64, Matrix{ComplexF64}, Matrix{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Matrix{ComplexF64}(undef, 7809, 0), Matrix{Float64}(undef, 0, 0)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [-0.15475965282627968 -0.1546162027014435 … -0.15418753100175617 -0.15461620270183546; -0.1546162027010446 -0.15447229890010616 … -0.15404229467820021 -0.15447229890049863; … ; -0.15418753100254745 -0.15404229467939787 … -0.1536083373509718 -0.15404229467979613; -0.15461620270223392 -0.15447229890129838 … -0.15404229467940575 -0.15447229890169228;;; -0.15461620270200463 -0.15447229890107023 … -0.15404229467917038 -0.15447229890146066; -0.15447229890066963 -0.15432794309174705 … -0.15389660283124915 -0.15432794309213788; … ; -0.1540422946799647 -0.15389660283245157 … -0.15346128701504805 -0.15389660283284842; -0.15447229890186062 -0.15432794309294093 … -0.1538966028324566 -0.1543279430933333;;; -0.1541875310020756 -0.15404229467892794 … -0.15360833735048726 -0.15404229467931943; -0.1540422946785234 -0.15389660283100828 … -0.15346128701358352 -0.1538966028314001; … ; -0.15360833735128923 -0.15346128701479772 … -0.15302191314008712 -0.15346128701519576; -0.15404229467972325 -0.1538966028322111 … -0.15346128701480027 -0.15389660283260462;;; … ;;; -0.1534788058560306 -0.1533313919656818 … -0.1528909249293029 -0.15333139196609247; -0.15333139196527382 -0.15318350574295767 … -0.15274162878040792 -0.15318350574336884; … ; -0.15289092493011225 -0.15274162878163333 … -0.15229553959366934 -0.15274162878205141; -0.1533313919664999 -0.15318350574418696 … -0.15274162878165223 -0.15318350574460005;;; -0.15418753100058977 -0.15404229467743213 … -0.15360833734899318 -0.15404229467783415; -0.15404229467702946 -0.15389660282950426 … -0.15346128701208125 -0.1538966028299068; … ; -0.1536083373497921 -0.1534612870132905 … -0.153021913138582 -0.15346128701369943; -0.1540422946782364 -0.15389660283071419 … -0.15346128701330525 -0.1538966028311184;;; -0.1546162027012574 -0.15447229890031774 … -0.15404229467841876 -0.15447229890071368; -0.1544722988999182 -0.15432794309099038 … -0.1538966028304935 -0.15432794309138678; … ; -0.15404229467921168 -0.15389660283169337 … -0.1534612870142909 -0.15389660283209575; -0.1544722989011129 -0.15432794309218806 … -0.1538966028317045 -0.1543279430925859]), DFTK.RealSpaceMultiplication{Float64, SubArray{Float64, 3, Array{Float64, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [-0.015163898301249128 -0.012550129805831669 … -0.004966869076253849 -0.012550129789814846; -0.012550129807199842 -0.010362595827786758 … -0.0041351599677066486 -0.01036259578855573; … ; -0.004966869045167167 -0.004135160007451542 … -0.0020697818641760616 -0.00413515998314251; -0.012550129782871799 -0.010362595795265155 … -0.004135160044623287 -0.010362595793539363;;; -0.01255012979618061 -0.010362595799304673 … -0.004135160007365085 -0.010362595798784974; -0.010362595814085239 -0.008622039820763917 … -0.003897512473948602 -0.008622039806199579; … ; -0.0041351600267413415 -0.0038975125346421857 … -0.003484108297133604 -0.003897512536344991; -0.010362595798660803 -0.008622039804935634 … -0.003897512550055001 -0.008622039817008549;;; -0.004966869046886575 -0.004135159974125191 … -0.002069781764239018 -0.004135159972406378; -0.004135160010332145 -0.0038975124903150985 … -0.0034841081745642048 -0.003897512485240673; … ; -0.002069781842370665 -0.0034841082324482357 … -0.005620581061040689 -0.0034841082551748817; -0.004135159981082659 -0.0038975124715628488 … -0.0034841082285475074 -0.0038975124888355785;;; … ;;; -0.006547897700514364 -0.004876339614894927 … -0.0014866456613682234 -0.0048763396396385664; -0.004876339579092328 -0.0026779878284728573 … -0.004123055249892015 -0.0026779878684136046; … ; -0.001486645624624521 -0.0041230552003549545 … -0.00717089603794548 -0.004123055192849855; -0.004876339667187573 -0.002677987948705742 … -0.004123055195134985 -0.0026779879674540797;;; -0.004966869012995416 -0.004135159978630933 … -0.0020697817549874408 -0.004135159927130977; -0.0041351599888414204 -0.0038975125273119082 … -0.003484108165392602 -0.0038975124436009864; … ; -0.002069781642239319 -0.0034841081121328114 … -0.0056205809649689294 -0.0034841080984977764; -0.004135159888244439 -0.003897512408020897 … -0.003484108168630322 -0.0038975123910042282;;; -0.012550129780313612 -0.010362595805515597 … -0.0041351599928579836 -0.010362595773301248; -0.010362595801500103 -0.008622039839843049 … -0.0038975124496824077 -0.008622039777262096; … ; -0.00413515991335745 -0.003897512442466336 … -0.003484108196696028 -0.0038975124134406043; -0.010362595756007464 -0.008622039777694702 … -0.003897512503264418 -0.008622039765953288])], DFTK.FourierMultiplication{Float64, Vector{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.0, 0.19739208802178715, 0.7895683520871486, 1.7765287921960844, 3.1582734083485944, 4.934802200544679, 7.106115168784338, 9.67221231306757, 12.633093633394378, 15.988759129764759  …  20.133992978222288, 16.38354330580833, 13.027877809437953, 10.066996489111146, 7.5008993448279115, 5.329586376588253, 3.5530575843921683, 2.1713129682396586, 1.184352528130723, 0.5921762640653614]), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [-0.009532893058674202 -0.006919175933467272 … 0.0006640002846683146 -0.006919175917842385; -0.006919175934436403 -0.004731679320205101 … 0.001495696020072922 -0.004731679281366537; … ; 0.000664000314963796 0.0014956959791302925 … 0.00356104473741597 0.001495696003041059; -0.006919175911297631 -0.004731679288875683 … 0.0014956959419507752 -0.004731679287543799;;; -0.006919175924377375 -0.004731679292687245 … 0.0014956959794441785 -0.004731679292557956; -0.004731679307067048 -0.002991148409853279 … 0.0017333318383345934 -0.0029911483956797675; … ; 0.001495695959273696 0.001733331776438527 … 0.002146714138344724 0.0017333317743388724; -0.004731679292833576 -0.0029911483952188463 … 0.001733331761020827 -0.0029911484076841424;;; 0.0006640003137161222 0.001495696012926517 … 0.0035610448378375426 0.001495696014253838; 0.0014956959771241837 0.0017333318222089874 … 0.0021467142623786736 0.0017333318268915876; … ; 0.0035610447589039538 0.0021467142032803917 … 1.0228217827559088e-5 0.002146714180155703; 0.0014956960051739077 0.0017333318397584215 … 0.0021467142071786775 0.0017333318220921733;;; … ;;; -0.0009170696309752906 0.0007544789992316378 … 0.004144156701975067 0.0007544789740773269; 0.0007544790354422713 0.0029528256321427007 … 0.0015077443588424845 0.0029528255917907822; … ; 0.004144156737909445 0.0015077444071540805 … -0.0015401031090828306 0.0015077444142410981; 0.0007544789461209796 0.002952825510680577 … 0.0015077444123552872 0.002952825491519153;;; 0.0006640003490931204 0.0014956960099165514 … 0.0035610448485832026 0.001495696061014495; 0.0014956960001088525 0.0017333317867161968 … 0.002146714273052547 0.0017333318700245796; … ; 0.0035610449605324077 0.002146714325103027 … 1.022831540442061e-5 0.002146714338329139; 0.001495696099498938 0.0017333319047972869 … 0.002146714268590889 0.0017333319214097513;;; -0.006919175907763142 -0.004731679298145688 … 0.0014956959947029008 -0.004731679266327244; -0.004731679293730485 -0.002991148428175711 … 0.0017333318633564608 -0.0029911483659911635; … ; 0.0014956960734105688 0.001733331869372576 … 0.002146714239539416 0.0017333318979959348; -0.0047316792494325015 -0.0029911483672250445 … 0.0017333318085635025 -0.0029911483558814796]), DFTK.NonlocalOperator{Float64, Matrix{ComplexF64}, Matrix{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Matrix{ComplexF64}(undef, 7809, 0), Matrix{Float64}(undef, 0, 0)), nothing, @NamedTuple{ψ_reals::Array{ComplexF64, 3}}[(ψ_reals = [-1.792692006397377e-6 - 1.879532065461179e-6im -0.002154995869600181 - 0.008159182992142588im … -0.0015172495788999365 - 0.0056214230507541756im 0.0022032693521072864 + 0.008211304357733946im; -0.018902859458598476 + 0.01260986936314339im 0.01734916209429515 - 0.005964264095899431im … -0.008233186020221588 + 0.01006931509052314im 0.014794324120422928 - 0.015555443924626901im; … ; -0.02098737947912027 + 0.013969046840300387im 0.017505875067525614 - 0.013584121839729273im … -0.011819870732101986 + 0.00568096421376741im 0.018383342893404944 - 0.010292785144521564im; 0.01886398707111428 - 0.012648117713826595im -0.014804339724046366 + 0.015538271943687995im … 0.010321354353083064 - 0.002339148954774984im -0.017359158920407264 + 0.005947557497362703im;;; 0.00545418221166115 - 0.00018245089472699282im -0.0019336873465758263 + 0.004882932412904674im … 0.0021414508601889882 + 0.0038280303126769343im -0.004488391923904938 - 0.004707313298633637im; 0.012890538642494485 - 0.010635156433053267im -0.012116278683501673 + 0.005787255569993712im … 0.006146564498497724 - 0.008069640370904881im -0.010471141866381338 + 0.011961205448850765im; … ; 0.019035952533628567 - 0.01198287146837983im -0.01579536829642167 + 0.011488664605111568im … 0.010760809178041552 - 0.005016737746543401im -0.016484176631544976 + 0.008905108678250753im; -0.01926247318508363 + 0.01086757475824681im 0.01459814488592166 - 0.01207805742568557im … -0.009432520930171415 + 0.00234905937070618im 0.01624343512628981 - 0.005904343621906451im;;; -0.003733168551375142 + 0.00010265774341876739im 0.0021321726205590463 - 0.0017048090482220746im … -0.001737180163597194 - 0.001829169661050486im 0.003009543116255037 + 0.001586155576833117im; -0.006712245923662295 + 0.006137338755796289im 0.006212702557012756 - 0.003877590060804658im … -0.003491338422018701 + 0.004847229597456153im 0.005524112634174556 - 0.006460819799757898im; … ; -0.01255006837445687 + 0.007571155048925159im 0.010682403909971747 - 0.007323300379816244im … -0.007762001850540735 + 0.003537050587103942im 0.011078387391253636 - 0.005838243188250956im; 0.01184208127380264 - 0.00627085350456307im -0.009366231338039481 + 0.00654113106698403im … 0.006408128966417358 - 0.0017730597298363147im -0.010054971851898568 + 0.003957633029010718im;;; … ;;; -0.001609780841671514 + 0.00010609232861639027im 0.0015282840523817992 + 0.0005525846303940959im … -0.0006552363768785833 + 0.0009687813222265778im 0.0011836183692828387 - 0.0007410858326961068im; 0.005686973315956705 - 0.002990021746619121im -0.005143820657287287 + 0.0020845847240438756im … 0.0033409842448935663 - 0.0027160785796837635im -0.004844652268781002 + 0.0032073126514171935im; … ; 0.005204663678462503 - 0.004007345154386582im -0.004601629807293648 + 0.004008814748055794im … 0.0037437545360556605 - 0.002132664417029576im -0.004812127483788089 + 0.003218785059265791im; -0.0029742445790000177 + 0.0028022239882999614im 0.0025151230197022305 - 0.0030372960092633476im … -0.0021774604362605407 + 0.0009745108838637457im 0.002814203106587497 - 0.0019146872166343402im;;; 0.003671273716598434 - 0.000199741975583911im -0.0029512812202858382 - 0.0014961657738222904im … 0.0007046750310346513 - 0.0019308050833490992im -0.002073713513186282 + 0.0017956746982491522im; -0.011795522387527423 + 0.006346454756186203im 0.010006780852082967 - 0.004033892558409724im … -0.005537951228596696 + 0.004932552082255842im 0.0093178973626252 - 0.006617927083327im; … ; -0.010108148014552471 + 0.007469701300561211im 0.008635944717519547 - 0.007238274847668404im … -0.006426503342121579 + 0.0034814408246005187im 0.009031742440006554 - 0.0057531852262641785im; 0.006758655690616188 - 0.006061755655937394im -0.0055724592006675695 + 0.006384625693709723im … 0.004361571832280418 - 0.0016879704193410445im -0.006261064776888767 + 0.0038008377778841838im;;; -0.0054055718625030825 + 0.00025753161693109234im 0.004428077210848355 + 0.00462333382373311im … -4.125830273191287e-5 + 0.003919054573759201im 0.0018729948654577435 - 0.004968276544969449im; 0.01925199284117718 - 0.01089552536010947im -0.016211437285149184 + 0.005955621730656448im … 0.007860121389942838 - 0.008141127117404168im -0.014566043907279447 + 0.01213039464685807im; … ; 0.01685303828254397 - 0.01189202350578851im -0.014081824671231234 + 0.011417357929107063im … 0.00977590489909395 - 0.004975629270527988im -0.014770568342055646 + 0.008833737862196902im; -0.012900775894723448 + 0.010607224512477292im 0.010503080441476688 - 0.011909670252751386im … -0.007719100701966054 + 0.002277625016210978im 0.012148381135951535 - 0.005735180584483325im],)])]), [8.519874499329698e-7 4.6272656408871753e-7 … 2.4177792090598293e-8 4.6272656218771556e-7; 4.6272656425110153e-7 2.5014480559851237e-7 … 1.3583004521541882e-8 2.501448025642503e-7; … ; 2.417779161360795e-8 1.3583004931574953e-8 … 1.5659253239800133e-9 1.3583004680791738e-8; 4.6272656136366167e-7 2.50144803083175e-7 … 1.3583005315065153e-8 2.5014480294970185e-7;;; 4.6272656294325196e-7 2.5014480339560705e-7 … 1.3583004930686571e-8 2.5014480335541846e-7; 2.5014480453879263e-7 1.3895611376367914e-7 … 1.127999713038874e-8 1.389561130148697e-7; … ; 1.3583005130580821e-8 1.1279997681462201e-8 … 7.937762581546914e-9 1.1279997696922973e-8; 2.5014480334581517e-7 1.389561129498771e-7 … 1.1279997821408526e-8 1.3895611357060613e-7;;; 2.417779163998531e-8 1.3583004587760597e-8 … 1.565925089304647e-9 1.3583004570028788e-8; 1.3583004961292512e-8 1.1279997278990888e-8 … 7.937761707296433e-9 1.1279997232916905e-8; … ; 1.565925272775267e-9 7.937762120164256e-9 … 3.5718578932778643e-8 7.937762282267064e-9; 1.358300465953672e-8 1.1279997108728456e-8 … 7.937762092341024e-9 1.1279997265558854e-8;;; … ;;; 5.7918418076397776e-8 2.2816020894285736e-8 … 5.620387460356709e-10 2.2816021259077584e-8; 2.2816020366451477e-8 3.4892753711687533e-9 … 1.3458519302724154e-8 3.4892755333507134e-9; … ; 5.620387031213888e-10 1.345851879489773e-8 … 7.7283175362942e-8 1.345851871795824e-8; 2.2816021665234687e-8 3.489275859388662e-9 … 1.3458518741386932e-8 3.4892759355186397e-9;;; 2.417779111996755e-8 1.3583004634246533e-8 … 1.5659250675808827e-9 1.3583004102938762e-8; 1.3583004739580935e-8 1.1279997614907959e-8 … 7.937761641876316e-9 1.1279996854847681e-8; … ; 1.5659248028206145e-9 7.937761261991596e-9 … 3.571857700341443e-8 7.93776116473895e-9; 1.3583003701761916e-8 1.127999653179405e-8 … 7.937761664971727e-9 1.1279996377290022e-8;;; 4.6272656106004927e-7 2.5014480387598744e-7 … 1.3583004781020673e-8 2.5014480138441885e-7; 2.501448035654079e-7 1.3895611474461984e-7 … 1.1279996910064291e-8 1.3895611152707316e-7; … ; 1.358300396084344e-8 1.1279996844543454e-8 … 7.93776186515397e-9 1.1279996581002551e-8; 2.5014480004685e-7 1.3895611154931015e-7 … 1.1279997396564884e-8 1.389561109456402e-7;;;;], Matrix{ComplexF64}[[-0.14482554869327116 - 0.2165468324577476im 0.5753032203209336 - 0.7748585778204558im 8.00407275583946e-7 - 1.1632923045198e-6im 4.600753831678916e-6 - 3.3543222016757152e-6im; 0.09826107969472447 + 0.14692245783957528im 0.025350121595468043 - 0.03414330957614817im -0.3170187742788721 - 0.10402316797192403im 0.386458074164333 + 0.03613659683248368im; … ; -0.031218279280165958 - 0.046678362673813305im -0.009168810975199013 + 0.012349195319959911im -0.011896782138976852 + 8.46949640585939e-5im 0.005558712266365617 + 0.0013425572226068im; 0.05597458852520891 + 0.08369462394572431im 0.018164116273407492 - 0.024464700851738158im 0.0490201133023294 - 0.006040542485377242im -0.010141034682531658 - 0.005504496458331454im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944070382873, [[-0.5541104229972931, -0.01067839107928148, 0.14495657035137074, 0.14495680204030933]], 1.0e-6, DFTK.BandtolBalanced{Float64, Vector{Vector{Float64}}}([[0.6300127939257145]], 1.1102230246251565e-16, Inf, 1.0e-6), 100, [0.0, 0.0, 0.0]), info_gmres = (x = [-3.6206941494755575e-7, -6.293616494036539e-7, -2.295121252107262e-7, 4.994025659838541e-7, 7.923312381075025e-7, 5.851072996783246e-7, 6.516971590544323e-7, 1.0189308596776472e-6, 3.638591402552996e-7, -2.201093323270999e-6  …  2.7810181883655596e-6, 1.2016955532822936e-6, 2.4456439114193986e-8, 1.0078548421755235e-7, 5.535307822198788e-7, 4.534172539953229e-7, 4.631917018570319e-9, -1.0437049982185616e-7, 1.1211101307213645e-7, 1.2550869621233512e-7], resid_history = [0.24939209548847494, 0.003766548609279519, 0.0002852632754436215, 4.6810178600517805e-6, 3.902845589101416e-8, 6.830208687304678e-11, 4.185823954209537e-8, 1.069527265273756e-9], converged = true, n_iter = 8, residual_norm = 1.069527265273756e-9, maxiter = 100, tol = 1.0e-8, s = 0.9355928859642837, residual = [9.176268738389954e-7, 6.329017582349217e-7, 1.367425487728674e-7, -2.0345122393956156e-7, -2.425991773445646e-7, -1.4428936346802728e-7, -1.3582501600386124e-7, -1.8616107402348812e-7, -6.009040676378617e-8, 3.354731049668492e-7  …  7.776989548804456e-7, 3.472612490145683e-7, 7.578802985101463e-9, 3.3952931261507924e-8, 2.0197109399864353e-7, 1.8227917697613426e-7, 2.189983849915268e-9, -6.53109192501668e-8, 1.1379906248933432e-7, 3.8326229419649837e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.17911354011892e-7, -4.0860921687205376e-8, 1.069527265273756e-9, -0.00020445246025578174, 4.68465890156713e-6, -3.9533495181045924e-8, 6.830208687304678e-11, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], V = KrylovKit.OrthonormalBasis{Vector{Float64}}([[-5.09508638031772e-7, -2.199076022947106e-6, -7.579824370251051e-7, 1.252257615187114e-6, 1.5497019519185292e-6, 1.1063338981373843e-6, 1.3261445730816041e-6, 1.9805837193497777e-6, 5.564100168622662e-7, -2.596395995102225e-6  …  3.177327875015803e-6, 1.2431535996859378e-6, 2.9531659294534158e-8, 1.4512365548913963e-7, 7.846258284882747e-7, 5.271471521792388e-7, 4.617201709521459e-9, -1.231849290424347e-7, 1.9334175520713754e-7, 3.140958206659681e-7], [8.777057802900597e-7, 2.190047940359829e-6, 7.282558085955668e-7, -1.2137675502336053e-6, -1.5102289979496046e-6, -1.0567859564466488e-6, -1.2404867602322453e-6, -1.8514786872994914e-6, -5.295264631380461e-7, 2.504111396838794e-6  …  -2.242360709510154e-6, -8.482932649153535e-7, -2.0409688505104453e-8, -1.0269221537548119e-7, -5.502395514263288e-7, -3.45170492555366e-7, -2.6452266305915045e-9, 6.46995696081743e-8, -9.378480896834223e-8, -3.5198115024770256e-8]]), H = [1.1090909807573242 0.015019477769608297 … 0.0 0.0; 0.1328374021589868 1.032854445952726 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1170177165155402 0.13774147577411044 … 0.0 0.0; 0.0 1.0240731668167375 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]))

From the result of solve_ΩplusK_split we can easily compute the polarisabilities:

println("Non-interacting polarizability: $(dipole(basis, res.δρ0))")
println("Interacting polarizability:     $(dipole(basis, res.δρ))")

Non-interacting polarizability: 1.9257125252248755
Interacting polarizability:     1.7736548556866976

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

Manual solution of the Dyson equations

To see what goes on under the hood, we also show how to manually solve the Dyson equation using KrylovKit:

using KrylovKit

# Apply ``(1- χ_0 K)``
function dielectric_operator(δρ)
    δV = apply_kernel(basis, δρ; scfres.ρ)
    χ0δV = apply_χ0(scfres, δV).δρ
    δρ - χ0δV
end

# Apply ``χ_0`` once to get non-interacting dipole
δρ_nointeract = apply_χ0(scfres, δ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 Main conflicts with an existing identifier.
[ Info: GMRES linsolve starts with norm of residual = 4.19e+00
[ Info: GMRES linsolve in iteration 1; step 1: normres = 2.49e-01
[ Info: GMRES linsolve in iteration 1; step 2: normres = 3.77e-03
[ Info: GMRES linsolve in iteration 1; step 3: normres = 2.85e-04
[ Info: GMRES linsolve in iteration 1; step 4: normres = 4.69e-06
[ Info: GMRES linsolve in iteration 1; step 5: normres = 1.09e-08
[ Info: GMRES linsolve in iteration 1; step 6: normres = 6.27e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 8.69e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 6.99e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 3.10e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125252247633
Interacting polarizability:     1.7736548493841404

We obtain the identical result to above.

HS2025M. F. Herbst and B. Sun. Efficient Krylov methods for linear response in plane-wave electronic structure calculations. arXiv 2505.02319