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.00013457378763965757

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.017612223704441047

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457378763965757
Displaced dipole:  0.017612223704441047
Polarizability :   1.7746797492080704

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.0011111926746701255 - 0.0031029296713632453im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.29039608310873166 - 0.106753494607299im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.03860478290403846 - 0.014000627389795173im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.07630346910565387 + 0.026995785539201143im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620694094929146e-7 -2.5020146129106006e-7 … -4.88064323154221e-8 -2.502015103633533e-7; -6.293615522880875e-7 -4.793035281477448e-7 … -1.2077840355881162e-7 -4.793035617640851e-7; … ; 1.3626255707144997e-7 1.1356225666911447e-7 … 4.711986313943294e-8 1.1356224583976462e-7; 1.3305467940277348e-7 1.363263598292854e-7 … 5.2006821995565663e-8 1.363263210042436e-7;;; -2.502014430887211e-7 -1.7404900770236488e-7 … -3.556493412995826e-8 -1.7404904396548402e-7; -4.793035145507272e-7 -3.679730836633563e-7 … -1.1202890398891496e-7 -3.6797310904121925e-7; … ; 1.1356225701674932e-7 1.1211094424460717e-7 … 1.0980677993408265e-7 1.1211093424601567e-7; 1.3632637051537242e-7 1.2550864659317588e-7 … 5.0949059660424804e-8 1.2550861723768724e-7;;; -4.880640291007727e-8 -3.5564914875582676e-8 … -1.1412139593440352e-8 -3.556492419007575e-8; -1.207783869265392e-7 -1.1202888877687138e-7 … -9.786220868856838e-8 -1.1202890044826197e-7; … ; 4.711987120838806e-8 1.098067967728301e-7 … 2.493822324985087e-7 1.0980678973176259e-7; 5.200684537537811e-8 5.094907959143934e-8 … 4.910625547326893e-8 5.094907107203143e-8;;; … ;;; 6.475765281690107e-8 4.110335950410908e-8 … -6.681697621694953e-9 4.110337017625309e-8; 1.68975600815876e-7 6.653287491971814e-8 … -1.3294450106045551e-7 6.653287702105374e-8; … ; 3.2677611022344195e-8 1.6017361784452986e-7 … 3.847034997123091e-7 1.6017361344374432e-7; -9.141988780633268e-8 -3.58303409533708e-8 … 7.037279728474913e-8 -3.583033597231567e-8;;; -4.880644236797853e-8 -3.556494529765918e-8 … -1.141215016178204e-8 -3.556495363397724e-8; -1.2077841503483883e-7 -1.1202891648268882e-7 … -9.7862231114047e-8 -1.1202892583545715e-7; … ; 4.7119863285352735e-8 1.0980677637374522e-7 … 2.4938219170433874e-7 1.0980677106735768e-7; 5.200681933304681e-8 5.0949056425871874e-8 … 4.9106237385809276e-8 5.0949048028108633e-8;;; -2.502015295554556e-7 -1.7404907138156636e-7 … -3.5564949760355584e-8 -1.7404910752418527e-7; -4.793035780498058e-7 -3.6797313114597645e-7 … -1.1202892043226281e-7 -3.679731575820468e-7; … ; 1.1356224525315363e-7 1.1211093319682511e-7 … 1.0980677226857566e-7 1.1211092431455717e-7; 1.3632631013238613e-7 1.2550860240487384e-7 … 5.094904927353027e-8 1.255085733433301e-7;;;;], δHψ = Matrix{ComplexF64}[[-0.0006463666736284985 + 0.0018049348991481258im 0.03178891453926747 - 0.0898387134280544im 1.199759929156 + 0.5561255046341826im 0.3472459588641802 - 1.0467462486716639im; 0.1860764158144633 + 0.06855039100750034im -1.4993557886465734 - 0.63135364993726im 0.2657794263991552 + 0.12263848904192272im 0.07658491739004082 - 0.2313178340397604im; … ; 0.046037455575326104 + 0.0168376314147263im 0.019518818415273614 + 0.006851495007219983im -0.017480375451211238 - 0.011086703145477377im -0.00575077024396883 + 0.025263337716800605im; -0.06822853169489974 - 0.024144029925419076im -0.02410604051719184 - 0.00859061476076375im 0.01774771409779361 - 0.027734140809244138im -0.03913308324451685 + 0.10030425905973007im]], δVind = [0.0019654319554949745 0.0020810004560778647 … 0.0031536657563999045 0.0020810008970083165; 0.003889130655301682 0.004779268969379111 … 0.010276671946165047 0.0047792694753278735; … ; -0.006150846257496305 -0.008269602165450838 … -0.01730296796187828 -0.008269601428677577; 0.00023868271596364497 -0.00039869711785493886 … -0.003664693441982729 -0.00039869664344201433;;; 0.00208100032606176 0.0022232277411621602 … 0.0034201199432946673 0.002223228236725695; 0.00477926882992651 0.005734952702176801 … 0.01090359796458065 0.005734953263244861; … ; -0.008269602538145876 -0.009601007886016111 … -0.012623110932421178 -0.009601007134073624; -0.00039869728435196813 -0.0010730132109081454 … -0.0042306801437026935 -0.001073012683062241;;; 0.0031536637510105 0.0034201179675041244 … 0.004832504854704787 0.0034201188034873294; 0.010276669547808436 0.010903595855871986 … 0.012272140939192157 0.010903596577701107; … ; -0.0173029711855794 -0.012623112950266564 … -0.009620265129936915 -0.01262311207692296; -0.0036646956288909696 -0.004230682081109542 … -0.005491108829215544 -0.004230681232186096;;; … ;;; -0.0023365619741520945 -0.0028154646385163116 … 0.005693323647270163 -0.002815465385666356; -0.012913094277302152 -0.01784118181790994 … 0.01149834259957508 -0.01784118307926784; … ; -0.02517446036239741 -0.012810883941633211 … -0.008395872296536601 -0.012810882979514712; 0.007597976186239108 0.01022540988875629 … -0.0054497301292915825 0.010225408908389171;;; 0.003153666184155816 0.00342012064524438 … 0.004832508721297502 0.003420121493729097; 0.010276672264673202 0.010903598577896923 … 0.01227214385697618 0.010903599386002065; … ; -0.017302966755188068 -0.012623109835403414 … -0.009620262754913202 -0.012623108869765443; -0.00366469277036186 -0.004230679104471266 … -0.005491105646353533 -0.0042306782692065;;; 0.0020810010292288228 0.0022232285284156904 … 0.0034201212611709346 0.0022232290227819553; 0.004779269608157884 0.0057349535669491206 … 0.010903599232939903 0.00573495412344792; … ; -0.0082696010570734 -0.009601006377402044 … -0.012623109244889645 -0.009601005570588694; -0.00039869647229021977 -0.001073012300782233 … -0.004230678615270362 -0.0010730117715667686;;;;], δρ0 = [-3.602140249429085e-7 -2.4888427406664115e-7 … -4.8528365015723615e-8 -2.4888427137402484e-7; -5.972249775885376e-7 -4.606977816183556e-7 … -1.1905765455146896e-7 -4.606977774539236e-7; … ; 1.2735786517185243e-7 1.083633636964803e-7 … 4.640715427652082e-8 1.083633616577054e-7; 1.035200232671971e-7 1.1957957694998955e-7 … 5.069504510334723e-8 1.195795760732394e-7;;; -2.488842720998356e-7 -1.7310672376502194e-7 … -3.5356972875335224e-8 -1.731067222431641e-7; -4.6069777983944653e-7 -3.5721703400903647e-7 … -1.109143655356361e-7 -3.572170313544011e-7; … ; 1.0836336144273146e-7 1.0846493692269454e-7 … 1.0862233026279272e-7 1.0846493487060563e-7; 1.19579575927553e-7 1.1611143939177404e-7 … 5.0208287206508576e-8 1.1611143809637693e-7;;; -4.852836628385258e-8 -3.535697378242842e-8 … -1.1341310605879969e-8 -3.535697440051157e-8; -1.1905766071484822e-7 -1.1091436862302264e-7 … -9.770666446428932e-8 -1.1091437233041971e-7; … ; 4.640715465902799e-8 1.086223318987673e-7 … 2.484523986757228e-7 1.0862233142377471e-7; 5.0695045598506317e-8 5.02082888485902e-8 … 4.926856582996679e-8 5.020828820842628e-8;;; … ;;; 6.434571286384687e-8 4.0838629877025236e-8 … -6.63792818192439e-9 4.0838629459386403e-8; 1.697503579474113e-7 6.683085807318296e-8 … -1.3340480133792148e-7 6.683085581414494e-8; … ; 3.2709969123129004e-8 1.6023520544020314e-7 … 3.8390139794498895e-7 1.602352095819333e-7; -9.271366412094116e-8 -3.6336211230512995e-8 … 7.126307793964095e-8 -3.633621057641195e-8;;; -4.852836779647458e-8 -3.535697589923046e-8 … -1.1341311669390637e-8 -3.5356975546383194e-8; -1.1905766053878395e-7 -1.1091437063802163e-7 … -9.770666556553571e-8 -1.109143720542753e-7; … ; 4.640715621439279e-8 1.0862233281569719e-7 … 2.484524029878225e-7 1.086223340947122e-7; 5.0695047793499e-8 5.020829131353275e-8 … 4.9268569118920925e-8 5.0208290818409634e-8;;; -2.488842743166122e-7 -1.7310672595853e-7 … -3.535697437295394e-8 -1.7310672423850637e-7; -4.6069778178103906e-7 -3.5721703590878896e-7 … -1.1091436913354377e-7 -3.572170342354844e-7; … ; 1.0836336368775014e-7 1.0846493862240915e-7 … 1.0862233327012396e-7 1.0846493768616807e-7; 1.1957957690090974e-7 1.1611144063490217e-7 … 5.0208289637771464e-8 1.1611143971883052e-7;;;;], δeigenvalues = [[0.0004950056397550781, 0.09679864450676383, 0.034030367873601974, 0.03920273601905468]], δoccupation = [[0.0, 0.0, 0.0, 0.0]], δεF = 0.0, ε_adj = 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.15475965283188264 -0.15461620270727358 … -0.1541875310068298 -0.15461620270718293; -0.15461620270710752 -0.15447229890639894 … -0.15404229468374772 -0.15447229890631037; … ; -0.1541875310071786 -0.15404229468426642 … -0.15360833735506507 -0.15404229468416747; -0.1546162027073548 -0.1544722989066498 … -0.15404229468399214 -0.15447229890655562;;; -0.15461620270701196 -0.1544722989063048 … -0.15404229468364708 -0.15447229890621259; -0.15447229890613817 -0.15432794309744519 … -0.15389660283620007 -0.15432794309735515; … ; -0.15404229468399558 -0.15389660283671866 … -0.15346128701853637 -0.15389660283661852; -0.1544722989063847 -0.15432794309769507 … -0.1538966028364436 -0.15432794309759956;;; -0.1541875310064986 -0.154042294683582 … -0.15360833735438215 -0.15404229468348782; -0.15404229468341413 -0.1538966028361315 … -0.15346128701795686 -0.15389660283603956; … ; -0.15360833735473253 -0.15346128701847853 … -0.15302191314298297 -0.1534612870183767; -0.1540422946836612 -0.15389660283638198 … -0.15346128701820114 -0.15389660283628456;;; … ;;; -0.15347880586341398 -0.1533313919733149 … -0.1528909249361858 -0.15333139197322457; -0.15333139197314274 -0.1531835057510774 … -0.15274162878778688 -0.15318350575098996; … ; -0.1528909249365488 -0.15274162878832584 … -0.1522955395995801 -0.15274162878822695; -0.15333139197340365 -0.15318350575134246 … -0.15274162878804434 -0.15318350575124834;;; -0.1541875310073772 -0.15404229468446123 … -0.15360833735527482 -0.15404229468437097; -0.15404229468429195 -0.15389660283700973 … -0.15346128701884879 -0.15389660283692205; … ; -0.15360833735563123 -0.15346128701937814 … -0.15302191314389527 -0.15346128701927972; -0.15404229468454683 -0.15389660283726844 … -0.1534612870191005 -0.15389660283717455;;; -0.1546162027074531 -0.1544722989067465 … -0.15404229468409583 -0.1544722989066562; -0.15447229890657938 -0.1543279430978867 … -0.15389660283664858 -0.1543279430977987; … ; -0.1540422946844474 -0.153896602837171 … -0.15346128701899522 -0.15389660283707254; -0.1544722989068295 -0.15432794309814055 … -0.1538966028368958 -0.15432794309804673]), 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.015163897562580887 -0.012550129044052543 … -0.00496686796877001 -0.012550128974748312; -0.012550129011953796 -0.010362594969968865 … -0.004135158744168625 -0.010362594917469834; … ; -0.004966868143943603 -0.0041351589124744275 … -0.002069780011858048 -0.004135158857808239; -0.012550129038539703 -0.010362595006626946 … -0.0041351587381326296 -0.010362594930296797;;; -0.012550128982369308 -0.010362594933432924 … -0.004135158672445105 -0.01036259488104108; -0.01036259494275343 -0.008622038866836688 … -0.0038975111078590023 -0.008622038830856008; … ; -0.004135158853492344 -0.003897511233609957 … -0.0034841067060037215 -0.0038975111822722296; -0.01036259493156662 -0.008622038871193432 … -0.003897511041953715 -0.008622038800216803;;; -0.004966868030918408 -0.004135158731777353 … -0.0020697798777563714 -0.004135158769684172; -0.004135158872317737 -0.003897511170956782 … -0.0034841067813678084 -0.003897511253736336; … ; -0.0020697800223237047 -0.003484106742399214 … -0.005620579712934885 -0.003484106730264807; -0.004135158769203376 -0.0038975111258236434 … -0.0034841065871324835 -0.0038975110965847196;;; … ;;; -0.006547897898260345 -0.004876340077297015 … -0.001486644036589344 -0.00487634005740074; -0.004876340052251273 -0.0026779887966291447 … -0.004123054211753725 -0.0026779887395214437; … ; -0.001486644062692845 -0.004123054124172286 … -0.007170895106344393 -0.004123054189100226; -0.004876340060438369 -0.002677988810293806 … -0.00412305414566025 -0.00267798877545244;;; -0.004966868132067499 -0.00413515888361163 … -0.002069780005653637 -0.00413515886192678; -0.00413515887396385 -0.0038975112195656974 … -0.0034841068083934857 -0.0038975112522931; … ; -0.00206978007363855 -0.003484106769012344 … -0.005620579783786201 -0.0034841067922110245; -0.004135158895283281 -0.0038975112620980657 … -0.003484106738962373 -0.0038975112327709737;;; -0.01255012905269016 -0.01036259501037039 … -0.004135158782619783 -0.01036259495724745; -0.010362594975740762 -0.008622038899466201 … -0.003897511187544803 -0.00862203887717869; … ; -0.004135158912573341 -0.003897511276172566 … -0.0034841067704246406 -0.0038975112508675015; -0.010362595011324405 -0.008622038950243304 … -0.003897511165431859 -0.008622038890565339])], 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.009532892325608923 -0.006919175177518232 … 0.0006640013870785177 -0.006919175108123324; -0.006919175145253285 -0.00473167846867998 … 0.0014956972380634386 -0.004731678416092381; … ; 0.000664001211556204 0.0014956970692388567 … 0.003561046585640702 0.0014956971240039935; -0.006919175172086411 -0.004731678505588888 … 0.0014956972438550462 -0.004731678429164565;;; -0.006919175115573401 -0.0047316784320500585 … 0.0014956973098874619 -0.004731678379565982; -0.004731678441203781 -0.0029911474616241863 … 0.0017333331994732644 -0.002991147425553467; … ; 0.0014956971284918064 0.001733333073203669 … 0.0021467157259862854 0.0017333331246415385; -0.0047316784302634685 -0.002991147466230786 … 0.0017333332651351075 -0.00299114739515865;;; 0.0006640013252612997 0.0014956972506203004 … 0.003561046720425305 0.0014956972128076559; 0.0014956971102478653 0.0017333331364440966 … 0.002146715651201735 0.001733333053756469; … ; 0.003561046575507613 0.0021467156896486077 … 1.0229563037512934e-5 0.00214671570188485; 0.0014956972131152293 0.001733333181326741 … 0.0021467158451928385 0.001733333210663087;;; … ;;; -0.0009170698361046434 0.0007544785291964335 … 0.004144158319871035 0.0007544785491830526; 0.0007544785544143992 0.002952824655866686 … 0.0015077453896018157 0.0029528247130618172; … ; 0.004144158293404575 0.001507745476644241 … -0.0015401021833924885 0.0015077454118151944; 0.0007544785459664277 0.0029528246419370146 … 0.0015077454554379135 0.0029528246768724996;;; 0.0006640012232336061 0.0014956970979067546 … 0.0035610465916353644 0.0014956971196818658; 0.0014956971077239268 0.0017333330869569391 … 0.002146715623284132 0.0017333330543172166; … ; 0.003561046523294042 0.002146715662135864 … 1.022949127387103e-5 0.002146715639035605; 0.0014956970861496714 0.0017333330441658612 … 0.002146715692463585 0.0017333330735868503;;; -0.0069191751863354 -0.0047316785094292275 … 0.0014956971992640313 -0.004731678456215971; -0.004731678474632316 -0.00299114749469518 … 0.0017333331193389893 -0.002991147472319683; … ; 0.0014956970689589495 0.0017333330301887275 … 0.0021467156611064834 0.001733333055592241; -0.004731678510466065 -0.0029911475457261357 … 0.0017333331412047698 -0.002991147485954356]), 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 = [-4.824169479849485e-7 + 1.5368085919870506e-6im 0.0017564158769558522 + 0.006952728848040345im … 0.001165160220493915 + 0.004809860873498852im -0.0017226001709848237 - 0.0070128301887597356im; 0.019961827052924958 + 0.007177586943151946im -0.018401287382080217 - 0.010305334445661472im … 0.009422622270193942 + 0.0003498980009959313im -0.016363690759827132 - 0.00212087121172097im; … ; 0.022750645565684693 + 0.008208226714078354im -0.019258470019909615 - 0.005674717535081686im … 0.012857850463070149 + 0.006087003141037451im -0.019958517341484794 - 0.008484775572164403im; -0.019976969301976524 - 0.007125092132105068im 0.016349155071419107 + 0.002139270316532979im … -0.011018913461724684 - 0.006970093823046864im 0.018388187164767616 + 0.010323481804627787im;;; 0.006540095638496514 + 0.0007847837077175467im -0.004863417415982892 - 0.0045284911577260474im … 0.000517505755031607 - 0.003152886023171298im -0.0028242093323979905 + 0.003656127008447475im; -0.02120580239208683 - 0.006692917630061603im 0.01772235621392028 + 0.008142773056921264im … -0.009036516296583034 - 0.0007255426454155663im 0.016410048649386853 + 0.0028726799599554046im; … ; -0.018295951281465578 - 0.006913653756985749im 0.015459432973467685 + 0.004825092853720521im … -0.01058835129587346 - 0.005177968740355785im 0.01600894152374279 + 0.007031103263216926im; 0.013545982996346578 + 0.0057516240349961im -0.011460651170607024 - 0.0023078667506040633im … 0.008158607230539426 + 0.0054321377982699025im -0.012773624617411227 - 0.007577943982421962im;;; -0.0044688554086634155 - 0.0005074975632946909im 0.0034243853451129515 + 0.0017396701426453608im … -0.001093530834748098 + 0.0014316050286340832im 0.002724169634870108 - 0.0010706361535964878im; 0.013291825343258852 + 0.004003346899182881im -0.011169160035094956 - 0.00442601658956682im … 0.006395151592921592 + 0.0007565944324509793im -0.010619701894264222 - 0.002220057025505733im; … ; 0.010998289998458451 + 0.004330247496099591im -0.009464112565862076 - 0.0031464622244859854im … 0.006927474331674857 + 0.0034949182046177907im -0.009779997354057605 - 0.004414563520439728im; -0.007149705097248035 - 0.003316726018716298im 0.006025845699236759 + 0.0017314761499684509im … -0.004558519324036804 - 0.0031660571989436323im 0.006575378347552724 + 0.003937614663807459im;;; … ;;; -0.0019516171131450202 - 0.0002969388269695505im 0.0015080563677296933 - 0.00029617960208470873im … -0.0012793441457121639 - 0.000946359592144963im 0.001783098829180488 + 0.0008082949767704493im; -0.0031462863864344043 - 0.0014604657954247605im 0.002943818887776914 + 0.0017725448675877326im … -0.0019154504019923122 - 0.00024377256581240177im 0.0027051957378269927 + 0.0008140956142676688im; … ; -0.007840451943407874 - 0.002598192804151394im 0.006979113574084315 + 0.00200869732253791im … -0.005465728029531162 - 0.002296396811898034im 0.007147231035935311 + 0.002683152799914222im; 0.006437888438891276 + 0.0019717183549360128im -0.005534849900115308 - 0.001263741049412178im … 0.004197210242477411 + 0.0019452028384316859im -0.005773629603863557 - 0.002222223238442415im;;; 0.004437922485666754 + 0.0006277793093184854im -0.002691123075450067 + 0.000960716452968882im … 0.0018443719578511218 + 0.0018055162460006411im -0.003391008821536271 - 0.001849639465448011im; 0.007175944007760612 + 0.0032242607581799785im -0.0066043310663319935 - 0.0038447312201926186im … 0.003932816170650861 + 0.0004431868104552675im -0.006055262291885962 - 0.001638808499982999im; … ; 0.013936260439646408 + 0.004704381234066954im -0.011926386935614812 - 0.0034599979436849802im … 0.008534316972287851 + 0.003699243177899711im -0.012242405088306405 - 0.004727913908962366im; -0.013265083156316245 - 0.004096024012733894im 0.010590112568662585 + 0.002312916927336167im … -0.007020825847468109 - 0.0034793266356420664im 0.011139664035331988 + 0.0045188734096625206im;;; -0.00652200087482779 - 0.0008806693953109662im 0.0027897325009257077 - 0.0035535062138681173im … -0.002108781648497906 - 0.003487085723752209im 0.004828422495019684 + 0.0046311852645997045im; -0.013552043890735518 - 0.005717748175881151im 0.012794941187746708 + 0.007515112637816804im … -0.0069748639536489606 - 0.00046308212003228936im 0.011483101269114471 + 0.0022451684264418164im; … ; -0.0209225788523877 - 0.007248240253244685im 0.01752112260675448 + 0.005087747367169803im … -0.011773427376852113 - 0.005328631718946675im 0.018070749011601755 + 0.0072935214551225706im; 0.021199156433283097 + 0.006727300168396124im -0.016387498761141997 - 0.0029356078677978963im … 0.01022022880211379 + 0.005694373128236625im -0.01770039837666159 - 0.00820554321641068im],)])]), [8.519873157033032e-7 4.6272647367521155e-7 … 2.4177775097539868e-8 4.627264654496907e-7; 4.627264698654967e-7 2.501447392518402e-7 … 1.3582991898794509e-8 2.5014473519137064e-7; … ; 2.417777778537623e-8 1.3582993635138632e-8 … 1.565920974304921e-9 1.358299307116744e-8; 4.6272647302091804e-7 2.5014474208710425e-7 … 1.3582991836526018e-8 2.5014473618345994e-7;;; 4.627264663541919e-7 2.501447364260218e-7 … 1.3582991158853651e-8 2.5014473237384784e-7; 2.5014473714688914e-7 1.3895606471834618e-7 … 1.1279984726839361e-8 1.3895606286843244e-7; … ; 1.358299302664486e-8 1.1279985868608492e-8 … 7.937751232522723e-9 1.127998540248192e-8; 2.501447362816675e-7 1.3895606494235368e-7 … 1.1279984128446896e-8 1.3895606129314725e-7;;; 2.4177776051136793e-8 1.3582991770961501e-8 … 1.5659206594024116e-9 1.3582992162030212e-8; 1.3582993220858086e-8 1.1279985299745712e-8 … 7.937751770071272e-9 1.127998605134934e-8; … ; 1.5659209988794797e-9 7.93775149211953e-9 … 3.57185518593873e-8 7.937751405570221e-9; 1.358299215707091e-8 1.1279984889951502e-8 … 7.937750384650454e-9 1.1279984624473837e-8;;; … ;;; 5.791842362196305e-8 2.2816027711467047e-8 … 5.620368484032885e-10 2.2816027418134717e-8; 2.281602734221682e-8 3.4892793024767214e-9 … 1.3458508660269695e-8 3.4892790705837897e-9; … ; 5.620368788902079e-10 1.3458507762430367e-8 … 7.728314346703486e-8 1.3458508428036842e-8; 2.2816027462919943e-8 3.4892793579647757e-9 … 1.345850798271106e-8 3.489279216486279e-9;;; 2.417777760315079e-8 1.3582993337370914e-8 … 1.5659209597346721e-9 1.3582993113656022e-8; 1.3582993237839147e-8 1.1279985741095527e-8 … 7.937751962833753e-9 1.1279986038244846e-8; … ; 1.5659211193806848e-9 7.937751681942647e-9 … 3.5718553282257514e-8 7.93775184741162e-9; 1.358299345778433e-8 1.1279986127271863e-8 … 7.937751467606214e-9 1.1279985860994078e-8;;; 4.6272647470039183e-7 2.5014474237663874e-7 … 1.3582992295481755e-8 2.501447382679134e-7; 2.501447396982573e-7 1.389560663959684e-7 … 1.1279985450356858e-8 1.389560652500828e-7; … ; 1.3582993636160005e-8 1.127998625505933e-8 … 7.937751692016928e-9 1.1279986025300429e-8; 2.5014474245041584e-7 1.389560690066336e-7 … 1.1279985249577619e-8 1.3895606593833676e-7;;;;], Matrix{ComplexF64}[[-0.0878304788189526 + 0.24526063088722008im 0.32192864140358807 - 0.909803036575945im -2.5493755956908098e-6 - 5.74956199677103e-6im 1.5591002342765895e-6 + 6.943772534223126e-6im; 0.05959112721220831 - 0.16640416449798762im 0.014185451497743834 - 0.040089532029642855im 0.17599981808186332 - 0.3789753985977211im -0.3312436672174695 - 0.10981875246134075im; … ; -0.01893254636933372 + 0.052867846457396235im -0.005130688240808128 + 0.014499846467005884im 0.001657063771596655 - 0.00745143074138698im -0.0007351420816532955 + 0.00022579581104574992im; 0.03394618528501495 - 0.09479241077594971im 0.010164295482406471 - 0.028725330148737217im 0.0002443999287395372 + 0.02093809997112832im -0.013640082790263243 - 0.007117192016330768im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239441275813093, [[-0.5541104248769135, -0.010678400639348329, 0.14495661553954473, 0.1449567918399179]], 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.620694094929146e-7, -6.293615522880875e-7, -2.2951203779506708e-7, 4.99402569030401e-7, 7.923311526967417e-7, 5.851071902275466e-7, 6.516971181729082e-7, 1.0189309314882702e-6, 3.6385931063167383e-7, -2.2010931004801185e-6  …  2.7810178550638047e-6, 1.2016953037930177e-6, 2.445628273422015e-8, 1.0078537009706923e-7, 5.535306956651433e-7, 4.534172167111358e-7, 4.631915210356978e-9, -1.0437052375536711e-7, 1.1211092431455717e-7, 1.255085733433301e-7], resid_history = [0.2493920943846231, 0.0037665483817179166, 0.0002852540577154798, 4.681303066966953e-6, 5.464529781309187e-8, 6.547191373133949e-11, 4.3509640605866655e-8, 1.2152859246844163e-9], converged = true, n_iter = 8, residual_norm = 1.2152859246844163e-9, maxiter = 100, tol = 1.0e-8, s = 0.9355762919943543, residual = [8.310238311016665e-7, 5.631588970613943e-7, 1.1990073719801048e-7, -1.7679840494895444e-7, -2.0919070387719782e-7, -1.2282758683129306e-7, -1.13653356246359e-7, -1.5395324453578462e-7, -4.966590606953656e-8, 2.7752156952970247e-7  …  7.534991587098935e-7, 3.3506583654900253e-7, 7.297626939805432e-9, 3.263089141655259e-8, 1.9325839867030094e-7, 1.7326244892269251e-7, 2.0684166907466365e-9, -6.139910788558812e-8, 1.0651570543942866e-7, 3.5625408160095644e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.44734234638729e-7, -4.2669341595551263e-8, 1.2152859246844163e-9, -0.00020444591269588118, 4.6855202770845754e-6, -5.530222440718846e-8, 6.547191373133949e-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}}([[-6.83153775527254e-7, -2.1619588544158727e-6, -7.277596129915345e-7, 1.2032314411965527e-6, 1.494144842107618e-6, 1.060205549435474e-6, 1.261819337404885e-6, 1.8913492552778758e-6, 5.4133701534794e-7, -2.5873887622088975e-6  …  3.248393952127769e-6, 1.2857468374903273e-6, 3.0041552670505763e-8, 1.4525802596387116e-7, 7.881165941098168e-7, 5.422895198340295e-7, 4.854730758449192e-9, -1.268939490897842e-7, 1.9288802131981047e-7, 3.202356432591413e-7], [1.2057821183438345e-6, 2.68602606264832e-6, 8.795495920425312e-7, -1.467129932329863e-6, -1.8331118554344074e-6, -1.2804016728004955e-6, -1.4972387027439838e-6, -2.2431205395310217e-6, -6.521057663574671e-7, 3.15339429441545e-6  …  -3.142160298980095e-6, -1.2078369751499662e-6, -2.7983793538888718e-8, -1.3602190709757214e-7, -7.293622803012574e-7, -4.757256591801528e-7, -3.856759558870073e-9, 9.41780258208082e-8, -1.3606730485824633e-7, -1.2230952896090128e-7]]), H = [1.1342289298637815 0.011692357039065257 … 0.0 0.0; 0.12716389992392407 1.0265922736359114 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1413351491931722 0.12599917461665366 … 0.0 0.0; 0.0 1.019295431774588 … 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.9257125182388806
Interacting polarizability:     1.7736548495209685

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 = 1.00e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.09e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 3.61e-12
└ * number of operations = 11
Non-interacting polarizability: 1.925712518238774
Interacting polarizability:     1.773654844007169

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