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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457367761609862
Displaced dipole:  0.01761222112394483
Polarizability :   1.7746794801560928

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.00195728605107362 + 0.002651782971097589im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.24738002831708594 + 0.1858206586639194im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.03294127643410866 + 0.024519770647517352im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.06530420562198984 - 0.04781581993024235im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620694476268427e-7 -2.5020154900844625e-7 … -4.880642016935065e-8 -2.5020148258798945e-7; -6.293616412697945e-7 -4.793036478135164e-7 … -1.2077843017626117e-7 -4.793035965165797e-7; … ; 1.3626259648695757e-7 1.1356228674836029e-7 … 4.711991507427953e-8 1.1356229883616078e-7; 1.3305471121170025e-7 1.3632634644896892e-7 … 5.2006864025091114e-8 1.3632639447427066e-7;;; -2.502015561421357e-7 -1.7404913541173031e-7 … -3.556494231692436e-8 -1.7404908719281752e-7; -4.793036512629118e-7 -3.679732295861531e-7 … -1.12028947290053e-7 -3.679731914295112e-7; … ; 1.1356228484398611e-7 1.1211096695711244e-7 … 1.0980683957687331e-7 1.1211097835404355e-7; 1.3632634037451929e-7 1.2550859860955793e-7 … 5.094908971806775e-8 1.255086341536478e-7;;; -4.880645354645205e-8 -3.556496582281385e-8 … -1.1412153526811193e-8 -3.556495507769662e-8; -1.2077845703771582e-7 -1.1202897052856693e-7 … -9.786225917697826e-8 -1.120289596512165e-7; … ; 4.7119904819190415e-8 1.098068196482711e-7 … 2.4938229363611526e-7 1.0980683040976105e-7; 5.200683882916597e-8 5.094906792804069e-8 … 4.9106276913364135e-8 5.094907875996984e-8;;; … ;;; 6.47575665765145e-8 4.1103315609114924e-8 … -6.681694261636078e-9 4.1103300708081306e-8; 1.689755348453664e-7 6.653283677227324e-8 … -1.329444843306276e-7 6.653282833587191e-8; … ; 3.2677655324542714e-8 1.6017369843633156e-7 … 3.847037227848879e-7 1.601737157792116e-7; -9.141993779916525e-8 -3.5830347738569894e-8 … 7.037288306932101e-8 -3.583035000531557e-8;;; -4.880641836303221e-8 -3.556493955110257e-8 … -1.1412145611886794e-8 -3.556492923061405e-8; -1.2077843342390482e-7 -1.1202894821048753e-7 … -9.7862240924867e-8 -1.1202893832701919e-7; … ; 4.711991602547002e-8 1.0980684160719647e-7 … 2.49382342498249e-7 1.0980685486602279e-7; 5.200686712841558e-8 5.094909350230309e-8 … 4.910630013053276e-8 5.094910532260593e-8;;; -2.502014764180604e-7 -1.7404907720346913e-7 … -3.556492949117433e-8 -1.740490290587072e-7; -4.79303594539937e-7 -3.679731873051395e-7 … -1.1202893559636979e-7 -3.6797314937407126e-7; … ; 1.1356230090629165e-7 1.1211098163616247e-7 … 1.0980685384140446e-7 1.1211099399856929e-7; 1.3632640027032272e-7 1.255086431137602e-7 … 5.094910377273045e-8 1.2550867924840783e-7;;;;], δHψ = Matrix{ComplexF64}[[0.001138528536623848 - 0.0015425086569378543im -0.08874287760464745 - 0.034730758129501126im -1.2703252998847634 - 0.3640976525196995im 0.5412140400465947 - 0.09042852989891581im; -0.1584710362995449 - 0.1192078841084443im -0.6803966742812625 + 1.4777467407398093im -0.28136308318076664 - 0.08010377137791289im 0.12041195497463124 - 0.019919128864645547im; … ; -0.039242832662536586 - 0.029376065686115407im 0.007490671089821814 - 0.019282467095989218im 0.012431914826243319 + 0.0076964408547678265im -0.0030919895539639995 - 0.004824554637539802im; 0.05839180992423619 + 0.04276053391555883im -0.009379959846263572 + 0.023810074685733824im -0.08264719138175712 + 0.001606580426756953im 0.07724973613343353 - 0.05156560895030607im]], δVind = [0.0019654319959856244 0.0020810009789342778 … 0.003153663849418723 0.0020810004283650135; 0.003889130820308781 0.004779269643463492 … 0.010276669238547085 0.004779269051810523; … ; -0.006150846319012278 -0.00826960095908111 … -0.017302965479208753 -0.008269601982973414; 0.00023868256945617106 -0.0003986967212768963 … -0.0036646951402109387 -0.0003986973436541394;;; 0.0020810010474015847 0.002223229080866506 … 0.0034201190469710483 0.0022232284668848934; 0.004779269729252726 0.005734954223068618 … 0.010903596300217392 0.005734953568528019; … ; -0.008269600846854108 -0.009601004955338163 … -0.012623110032262578 -0.009601005983895987; -0.0003986966473262482 -0.0010730118631592602 … -0.00423068033162547 -0.0010730125510664493;;; 0.003153665993560313 0.0034201212851389285 … 0.004832504459739289 0.0034201202293670688; 0.01027667150512699 0.010903598489293299 … 0.012272140430258847 0.01090359746130777; … ; -0.017302962242955806 -0.012623107687718212 … -0.009620264164982179 -0.012623108704287744; -0.003664692751700666 -0.004230678024770746 … -0.0054911072504238235 -0.004230679073699631;;; … ;;; -0.0023365590791093608 -0.0028154620865670364 … 0.005693310745346065 -0.002815461186645234; -0.012913091811040328 -0.017841181056974995 … 0.011498336154800731 -0.017841180025730298; … ; -0.02517445860903001 -0.012810886475336118 … -0.008395876571787506 -0.01281088729930795; 0.0075979804767501305 0.010225415484558086 … -0.005449735675478325 0.01022541711540136;;; 0.0031536636089552336 0.003420118641241562 … 0.004832500422308024 0.003420117564750005; 0.010276668799285538 0.010903595730591329 … 0.012272137483542478 0.010903594659806644; … ; -0.017302965880189033 -0.01262311047943491 … -0.009620266240950127 -0.012623111401974817; -0.003664695343938558 -0.004230680694551307 … -0.005491109999452818 -0.004230681721924484;;; 0.002081000363496741 0.002223228314939308 … 0.0034201177023271585 0.0022232276970050764; 0.004779268965085429 0.00573495337213193 … 0.010903594924543682 0.0057349527122804765; … ; -0.008269602077857353 -0.00960100621216316 … -0.012623111282467623 -0.009601007221141766; -0.00039869741115174627 -0.0010730127101989822 … -0.004230681630268983 -0.0010730134014964386;;;;], δρ0 = [-3.6021406375094134e-7 -2.488843035357549e-7 … -4.852837920817725e-8 -2.488843023483143e-7; -5.972250813563619e-7 -4.606978685944977e-7 … -1.1905770272286563e-7 -4.606978650709564e-7; … ; 1.2735791006158537e-7 1.0836340742014296e-7 … 4.6407200654287335e-8 1.0836340862398319e-7; 1.0352006960339984e-7 1.1957961652976905e-7 … 5.069506954024128e-8 1.1957961710784686e-7;;; -2.488843028645733e-7 -1.7310674856378503e-7 … -3.5356986797950014e-8 -1.7310674784995145e-7; -4.6069786649238457e-7 -3.5721711068534147e-7 … -1.109144158142867e-7 -3.57217107793727e-7; … ; 1.0836340687074906e-7 1.0846498454484817e-7 … 1.0862238872969798e-7 1.0846498607941539e-7; 1.195796162489926e-7 1.1611147497374669e-7 … 5.0208313942203735e-8 1.1611147543228161e-7;;; -4.8528379690294875e-8 -3.535698712916402e-8 … -1.1341321738028706e-8 -3.535698725565203e-8; -1.1905770504612753e-7 -1.1091441802370184e-7 … -9.770671039977546e-8 -1.1091441703140528e-7; … ; 4.640719862658408e-8 1.0862238621156013e-7 … 2.4845247986049434e-7 1.0862238877242985e-7; 5.069506883647224e-8 5.0208313143745355e-8 … 4.926859519513753e-8 5.020831411301439e-8;;; … ;;; 6.434571113654503e-8 4.083862634682572e-8 … -6.637936889753314e-9 4.083862530692551e-8; 1.697503474664701e-7 6.683083655788062e-8 … -1.3340484599634131e-7 6.683083354690994e-8; … ; 3.2710006143105656e-8 1.6023525981496226e-7 … 3.83901489169535e-7 1.6023526655968103e-7; -9.27136639920165e-8 -3.6336203215107415e-8 … 7.126310673353483e-8 -3.633620003814984e-8;;; -4.852838088132825e-8 -3.5356987967017146e-8 … -1.1341322913558565e-8 -3.5356988527220676e-8; -1.1905770845944246e-7 -1.109144202307149e-7 … -9.77067125909998e-8 -1.1091442026881165e-7; … ; 4.6407200718732686e-8 1.0862238777077037e-7 … 2.4845248562077853e-7 1.0862239277427516e-7; 5.0695069995988704e-8 5.0208314116546356e-8 … 4.926859800688014e-8 5.0208316057224795e-8;;; -2.4888430397048623e-7 -1.7310674930664436e-7 … -3.535698750742407e-8 -1.7310674877810726e-7; -4.6069786855947033e-7 -3.572171119180554e-7 … -1.1091441642012201e-7 -3.5721710940666477e-7; … ; 1.0836340932450815e-7 1.084649869131559e-7 … 1.0862239272972173e-7 1.0846498934341545e-7; 1.1957961706874562e-7 1.1611147577158987e-7 … 5.020831569116524e-8 1.1611147671910414e-7;;;;], δeigenvalues = [[0.0004950055338191922, 0.09679860443212389, 0.03371867864345182, 0.04923855983314967]], δ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.15475965281880386 -0.15461620269405899 … -0.15418753099406765 -0.15461620269425752; -0.15461620269412948 -0.15447229889328462 … -0.15404229467108085 -0.15447229889348355; … ; -0.15418753099393076 -0.15404229467087163 … -0.153608337342129 -0.15404229467107347; -0.15461620269418813 -0.15447229889334305 … -0.1540422946711403 -0.1544722988935424;;; -0.15461620269389997 -0.15447229889305555 … -0.1540422946708458 -0.154472298893252; -0.15447229889312472 -0.15432794308429443 … -0.1538966028234926 -0.1543279430844911; … ; -0.1540422946707116 -0.15389660282328743 … -0.15346128700556086 -0.15389660282348727; -0.15447229889318395 -0.15432794308435324 … -0.15389660282355277 -0.15432794308455064;;; -0.1541875309933587 -0.15404229467030117 … -0.153608337341544 -0.1540422946704966; -0.15404229467036948 -0.1538966028229466 … -0.15346128700521092 -0.1538966028231421; … ; -0.1536083373414118 -0.15346128700500875 … -0.15302191312996716 -0.1534612870052079; -0.15404229467042946 -0.15389660282300618 … -0.1534612870052723 -0.15389660282320275;;; … ;;; -0.1534788058504663 -0.15333139196021964 … -0.15289092492354958 -0.1533313919604307; -0.15333139196029572 -0.15318350573808404 … -0.15274162877524936 -0.15318350573829537; … ; -0.15289092492340142 -0.15274162877502312 … -0.15229553958675177 -0.15274162877523795; -0.15333139196035558 -0.15318350573814365 … -0.15274162877531033 -0.15318350573835576;;; -0.1541875309943894 -0.15404229467133002 … -0.15360833734259774 -0.15404229467153555; -0.15404229467140382 -0.15389660282397907 … -0.1534612870062687 -0.15389660282418496; … ; -0.15360833734245385 -0.15346128700604894 … -0.15302191313103244 -0.1534612870062582; -0.1540422946714627 -0.15389660282403764 … -0.15346128700632847 -0.15389660282424422;;; -0.15461620269441814 -0.1544722988935729 … -0.15404229467137578 -0.15447229889377437; -0.15447229889364486 -0.1543279430848137 … -0.15389660282402468 -0.15432794308501555; … ; -0.15404229467123567 -0.1538966028238106 … -0.15346128700609665 -0.1538966028240155; -0.15447229889370337 -0.15432794308487194 … -0.1538966028240839 -0.15432794308507428]), 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.015163898132798385 -0.012550129622197433 … -0.004966868767957299 -0.0125501295924145; -0.012550129617083957 -0.010362595618274543 … -0.004135159606481831 -0.010362595562910444; … ; -0.0049668687838788144 -0.0041351596252847105 … -0.0020697812084557003 -0.004135159651371174; -0.01255012959034653 -0.01036259556305592 … -0.004135159637539253 -0.010362595564977768;;; -0.012550129600283275 -0.01036259558439341 … -0.0041351596099610454 -0.010362595560341543; -0.01036259558743112 -0.00862203957971965 … -0.0038975120485400406 -0.00862203953084159; … ; -0.004135159612676137 -0.003897512025171048 … -0.0034841076793300507 -0.0038975120586785704; -0.01036259555064658 -0.008622039509405622 … -0.003897512044878932 -0.008622039518347982;;; -0.004966868793564934 -0.004135159634780826 … -0.002069781149758923 -0.004135159641456311; -0.004135159662351134 -0.0038975121034493677 … -0.0034841076766671903 -0.0038975120782127823; … ; -0.0020697811531808297 -0.0034841076315254155 … -0.005620580540315256 -0.0034841076835974228; -0.004135159615408137 -0.0038975120134858623 … -0.003484107664213848 -0.00389751206061098;;; … ;;; -0.006547897672237972 -0.004876339737471494 … -0.001486645274462339 -0.004876339650030223; -0.004876339661932651 -0.002677988136233526 … -0.004123054917225304 -0.002677988056230875; … ; -0.0014866450651616079 -0.004123054772844177 … -0.007170895752543574 -0.004123054879152715; -0.004876339736729695 -0.0026779882529624453 … -0.004123054929793989 -0.002677988100026918;;; -0.004966868855966326 -0.004135159686169711 … -0.002069781282813695 -0.004135159727213109; -0.004135159729749829 -0.0038975121467937366 … -0.0034841077219171747 -0.0038975121438668907; … ; -0.0020697812086724458 -0.003484107662265259 … -0.005620580615463342 -0.003484107761115462; -0.004135159674534213 -0.0038975120623745987 … -0.0034841077845420296 -0.003897512151512522;;; -0.012550129626237026 -0.010362595606944806 … -0.004135159664926644 -0.010362595593048512; -0.010362595608722582 -0.008622039593953827 … -0.0038975120610397003 -0.00862203955035682; … ; -0.004135159669796735 -0.0038975120782249826 … -0.003484107758552345 -0.00389751213016669; -0.01036259558500577 -0.008622039542932915 … -0.0038975121259907756 -0.00862203956479342])], 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.009532892882747633 -0.006919175742448526 … 0.0006640006006533817 -0.006919175712864101; -0.006919175737405411 -0.004731679103871343 … 0.0014956963884171 -0.004731679048706168; … ; 0.0006640005848688397 0.001495696369823359 … 0.003561045401979119 0.0014956963435350567; -0.006919175710726576 -0.0047316790487111175 … 0.001495696357300253 -0.004731679050832306;;; -0.006919175720375384 -0.004731679069761309 … 0.001495696385172809 -0.0047316790459058684; -0.00473167907286802 -0.002991148161356391 … 0.0017333322714997 -0.0029911481126750055; … ; 0.0014956963825919992 0.0017333322950738063 … 0.0021467147656354657 0.0017333322613664439; -0.004731679036142682 -0.0029911480911011482 … 0.001733332275100718 -0.0029911481002409066;;; 0.0006640005757546788 0.001495696360897648 … 0.0035610454612609007 0.0014956963540267362; 0.0014956963332591168 0.0017333322171364084 … 0.002146714768648296 0.0017333322421774836; … ; 0.0035610454579712214 0.002146714813992187 … 1.0228748672952773e-5 0.0021467147617210336; 0.0014956963801422175 0.0017333323070403226 … 0.0021467147810402985 0.00173333225971864;;; … ;;; -0.0009170695971346 0.0007544788821172294 … 0.004144157094634266 0.0007544789693474472; 0.0007544789575800494 0.0029528253292556614 … 0.0015077446966677635 0.0029528254090469813; … ; 0.004144157304083184 0.0015077448412750705 … -0.0015401028167633477 0.001507744734751704; 0.0007544788827231648 0.0029528252124671785 … 0.001507744684038182 0.002952825365190598;;; 0.0006640005123225842 0.0014956963084798922 … 0.003561045327152388 0.0014956962672309634; 0.0014956962648260824 0.0017333321727595599 … 0.002146714722340519 0.0017333321754805149; … ; 0.003561045401437522 0.002146714782212148 … 1.0228672459552599e-5 0.0021467146831526955; 0.0014956963199828852 0.0017333322571201335 … 0.002146714659655962 0.0017333321677756256;;; -0.006919175746847304 -0.004731679092830041 … 0.0014956963296772173 -0.0047316790791351974; -0.004731679094679621 -0.0029911481761098185 … 0.0017333322584679937 -0.002991148132714651; … ; 0.0014956963249472925 0.001733332241496707 … 0.0021467146858773502 0.0017333321893501078; -0.00473167907102129 -0.002991148125147138 … 0.0017333321934577435 -0.002991148147209981]), 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 = [3.2601433080683068e-6 - 1.6250971353469305e-6im 0.004336446854911009 + 0.003650056518692218im … 0.003005585306381865 + 0.002413351960088193im -0.004385143562398747 - 0.0035717586173295273im; -0.02028550211610194 - 0.015010726581115756im 0.014193679749653448 + 0.010232319728720646im … -0.011450824149650493 - 0.008580909516868677im 0.019306349808501138 + 0.014464262598201105im; … ; -0.02163361408809148 - 0.016033351942549158im 0.0191371319979318 + 0.014266501866050517im … -0.010388238254106306 - 0.007637107842420614im 0.017381214987918708 + 0.012813966206482223im; 0.020314753179653144 + 0.014956841957439083im -0.019294049639212718 - 0.014486268990039217im … 0.007318487276001192 + 0.005273554692935903im -0.014181550318411151 - 0.010254470034372389im;;; 0.008320316432255273 + 0.0005270982942748932im -0.007416921906359134 - 0.0024516789303550123im … -0.00039359548565716135 - 0.0015456540893078726im -0.0023037428074174954 + 0.0017803198427811004im; 0.011860408791329352 + 0.012060835633266254im -0.008910722980382882 - 0.008529783515056347im … 0.008089429980123544 + 0.006945842161718814im -0.01220362045693785 - 0.0112544103357176im; … ; 0.019937889655863412 + 0.013639188576618104im -0.017452402329025042 - 0.012109016727463921im … 0.00980781315343589 + 0.006707094506672324im -0.01607400943863707 - 0.010968921292926934im; -0.021627978231410466 - 0.012656736206365973im 0.01846870023668315 + 0.011679538576305843im … -0.007626010310311692 - 0.004654226574329196im 0.015175957654479802 + 0.00895500801425849im;;; -0.005660793229398414 - 0.00038176540250598633im 0.004757696175419782 + 0.0010016797163437006im … -0.0008702473754894236 + 0.0006616470647138322im 0.0030012755704193848 - 0.00045109846341136775im; -0.005499543554816579 - 0.006660365278092792im 0.004275249754057631 + 0.005019942969572291im … -0.004247139960538945 - 0.0042218784040434516im 0.00565399229128515 + 0.006160202494579364im; … ; -0.013248011852017277 - 0.008560191150002167im 0.01184343011984922 + 0.007763653497102048im … -0.007301015306823722 - 0.004694714520690796im 0.01105109387966699 + 0.007108321268137731im; 0.013270092120118149 + 0.00719347233724252im -0.01144018387449481 - 0.006580000436274164im … 0.005700476356542309 + 0.0031207355936520367im -0.010061657473631557 - 0.005439740754669291im;;; … ;;; -0.002536431843916442 - 0.00011330348041363724im 0.0017995877489677516 - 0.00019446223456682653im … -0.0019088779032763778 - 0.0004472287881732205im 0.0024898100960797772 + 0.00037635857703137253im; 0.00649535646090766 + 0.003303091737794342im -0.005394140000671205 - 0.0026648322500150957im … 0.0044749974216174125 + 0.0024406922565168255im -0.005993151997451304 - 0.003160255534067415im; … ; 0.004694527048381235 + 0.004468995573464836im -0.004457762997430332 - 0.004208782143742798im … 0.002911290229971169 + 0.0027976065431636254im -0.004036157771486814 - 0.0038601711349316586im; -0.0022065643223167478 - 0.0031200228905409392im 0.002295465655232502 + 0.0030110369601975724im … -0.001060939094945811 - 0.0016454842902810865im 0.0016964873504495322 + 0.002515658356807265im;;; 0.005716679551095139 + 0.0002959552456314243im -0.0030557763379691783 + 0.0005357413994815014im … 0.0028842082645021486 + 0.0008854777378329686im -0.004811763130075265 - 0.0009165751132231951im; -0.013313166125749835 - 0.0071258902708875235im 0.010107757667451962 + 0.005367679287571657im … -0.007393915623291834 - 0.004409529886602006im 0.011486341524220602 + 0.006507810071555475im; … ; -0.009493467694197472 - 0.008336606272642964im 0.008696722518360078 + 0.0075760451921850075im … -0.005247483153032331 - 0.004572138729650271im 0.007904282408249709 + 0.006920758750558307im; 0.005456828809196226 + 0.006728189656922495im -0.0056078924790339485 - 0.006232248639423472im … 0.0025537726162484726 + 0.002932970588758309im -0.004229538148525579 - 0.005092285424293266im;;; -0.0083613063672129 - 0.0004660112194113772im 0.0023595434757255835 - 0.0018687722105624992im … -0.003749822193173027 - 0.0017456418908461716im 0.007471900805633314 + 0.0023624699142783473im; 0.02163706004265429 + 0.01264271742658375im -0.015205906815758449 - 0.00890499347206549im … 0.010724109858691691 + 0.007102839174511525im -0.018498531518676958 - 0.011629229555854572im; … ; 0.016581156434826792 + 0.013439539426229666im -0.014817454639329007 - 0.011952131494353im … 0.008293283368639225 + 0.006616739448238632im -0.013439055617302832 - 0.01081197478788431im; -0.011851587318311058 - 0.012074976684679207im 0.01217354415013272 + 0.011304412338796865im … -0.0049913216460943455 - 0.004497200570704639im 0.008881153320468042 + 0.008580214005049875im],)])]), [8.519874193223702e-7 4.627265422936703e-7 … 2.4177787360144746e-8 4.627265387588062e-7; 4.6272654168675045e-7 2.501447893940999e-7 … 1.3583000794934454e-8 2.5014478511203703e-7; … ; 2.4177787604441065e-8 1.3583000988916374e-8 … 1.5659237841961261e-9 1.3583001258038465e-8; 4.627265385133621e-7 2.5014478512329844e-7 … 1.3583001115340209e-8 2.501447852719345e-7;;; 4.6272653969273516e-7 2.501447867736048e-7 … 1.358300083082652e-8 2.501447849133537e-7; 2.5014478700856143e-7 1.3895610137059506e-7 … 1.127999326784873e-8 1.3895609885758433e-7; … ; 1.3583000858836126e-8 1.1279993055669082e-8 … 7.937758174946051e-9 1.1279993359901236e-8; 2.5014478416350137e-7 1.389560977554703e-7 … 1.127999323461e-8 1.389560982152318e-7;;; 2.4177787753063464e-8 1.3583001086879561e-8 … 1.5659236463623236e-9 1.358300115574762e-8; 1.3583001371315264e-8 1.1279993766401467e-8 … 7.93775815595188e-9 1.1279993537267192e-8; … ; 1.5659236543966443e-9 7.937757833973989e-9 … 3.571856847528345e-8 7.937758205383886e-9; 1.3583000887022168e-8 1.1279992949572574e-8 … 7.93775806712781e-9 1.1279993377448225e-8;;; … ;;; 5.791841728341387e-8 2.2816022701428147e-8 … 5.620382941556824e-10 2.281602141228347e-8; 2.28160215877644e-8 3.4892766208654915e-9 … 1.3458515892397697e-8 3.4892762960052626e-9; … ; 5.620380497053307e-10 1.345851441227906e-8 … 7.728316559143408e-8 1.3458515502098507e-8; 2.2816022690489477e-8 3.4892770948564963e-9 … 1.3458516021244765e-8 3.489276473845523e-9;;; 2.417778871053638e-8 1.3583001617038364e-8 … 1.5659239588046452e-9 1.3583002040466101e-8; 1.3583002066636994e-8 1.1279994159952839e-8 … 7.937758478709713e-9 1.1279994133376943e-8; … ; 1.5659237847034978e-9 7.937758053230965e-9 … 3.5718569984447925e-8 7.937758758295735e-9; 1.3583001497001456e-8 1.1279993393463412e-8 … 7.937758925391767e-9 1.1279994202796263e-8;;; 4.627265427731125e-7 2.501447885178155e-7 … 1.3583001397882349e-8 2.501447874430196e-7; 2.5014478865531995e-7 1.3895610210244185e-7 … 1.1279993381342144e-8 1.389560998609393e-7; … ; 1.3583001448127938e-8 1.1279993537375056e-8 … 7.937758740012876e-9 1.1279994008987599e-8; 2.501447868209684e-7 1.3895609947923743e-7 … 1.1279993971070608e-8 1.3895610060318337e-7;;;;], Matrix{ComplexF64}[[0.15470708162487842 - 0.20960126261578432im -0.8987057202168979 - 0.3517209683113156im -2.138974509674092e-6 - 6.591397196956554e-6im 5.290203795685547e-6 + 1.6868422352052334e-5im; -0.10496549131024252 + 0.14221003503321855im -0.039600504601833196 - 0.0154982056509518im -0.1153682856274968 + 0.40125929652981224im -0.02868797066413442 - 0.17022997592563438im; … ; 0.03334832087646421 - 0.045181190726102875im 0.014322990722777665 + 0.005605501405594566im -0.002003481048618291 + 0.011013791662659428im -0.0029675442083513526 - 0.006846707132158903im; -0.05979376772546839 + 0.08101018444787216im -0.028374940692966318 - 0.011104927364204242im 0.004965623155972688 - 0.0395775268525448im 0.014892920171546743 + 0.02887157985363175im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239440689057155, [[-0.5541104211005515, -0.010678392680591564, 0.1449566935239377, 0.14495692537598134]], 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.620694476268427e-7, -6.293616412697945e-7, -2.2951210140690395e-7, 4.994025778737152e-7, 7.923312259204114e-7, 5.851072722589693e-7, 6.51697140081709e-7, 1.0189308631736487e-6, 3.638591650266266e-7, -2.2010932851887605e-6  …  2.781018125552819e-6, 1.2016955065958249e-6, 2.4456402087597785e-8, 1.0078545121141891e-7, 5.535307624709201e-7, 4.534172550350118e-7, 4.631924509514797e-9, -1.0437050458808499e-7, 1.1211099399856929e-7, 1.2550867924840783e-7], resid_history = [0.24939209661766476, 0.003766548902380488, 0.0002852742886401705, 4.680651223415575e-6, 6.479628520275775e-8, 5.094622121201125e-11, 4.5301985417821874e-8, 1.2214220697582768e-9], converged = true, n_iter = 8, residual_norm = 1.2214220697582768e-9, maxiter = 100, tol = 1.0e-8, s = 0.9356355703822955, residual = [7.656305343744348e-7, 5.431502784253346e-7, 1.2053806804487726e-7, -1.8326300993757801e-7, -2.228461164673338e-7, -1.35683827178738e-7, -1.3120013289550736e-7, -1.836370168054008e-7, -5.980484117665099e-8, 3.356086923177352e-7  …  5.864270690807978e-7, 2.636898922038434e-7, 5.7883607281522355e-9, 2.609709213825244e-8, 1.5676610165488332e-7, 1.4338015550087262e-7, 1.7474552315036511e-9, -5.282583663101538e-8, 9.336231541758966e-8, 3.2006148853361145e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.3200277562404963e-7, -4.469998652369094e-8, 1.2214220697582768e-9, -0.00020446019030966012, 4.6837096000863675e-6, -6.543468399374197e-8, 5.094622121201125e-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}}([[-2.78103100092434e-7, -1.912981889686593e-6, -6.707562002422784e-7, 1.092474031380886e-6, 1.3248246673002302e-6, 9.463888256286226e-7, 1.1480141420678934e-6, 1.711980544239268e-6, 4.6888139791503856e-7, -2.109405591312813e-6  …  2.7624862543270688e-6, 1.0726322263499884e-6, 2.6050468430269863e-8, 1.308407539299873e-7, 7.123262636245113e-7, 4.76317220980111e-7, 4.167836740107501e-9, -1.1418702976781742e-7, 1.8649597159666555e-7, 3.334951792812279e-7], [6.399424694309241e-7, 1.925814559923886e-6, 6.486362898516487e-7, -1.066640414914207e-6, -1.3039011198848719e-6, -9.12794279227142e-7, -1.0826607179349303e-6, -1.6137039370118954e-6, -4.5157978320668417e-7, 2.0731922766471927e-6  …  -1.956134716888577e-6, -7.355037381947421e-7, -1.8263054888958325e-8, -9.439744184560102e-8, -5.105292288254089e-7, -3.20294538064385e-7, -2.4809678606885434e-9, 6.346814547405438e-8, -9.785072990860817e-8, -7.675282599558382e-8]]), H = [1.1132012525111072 0.010982291325309462 … 0.0 0.0; 0.13767150670182135 1.0218043265456114 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1216819836075853 0.13631211331129076 … 0.0 0.0; 0.0 1.0130991101089153 … 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.9257125331111717
Interacting polarizability:     1.7736548627547701

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.28e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 1.34e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.32e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 2.15e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125331110558
Interacting polarizability:     1.7736548557857854

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