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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457444136012807
Displaced dipole:  0.017612219624110977
Polarizability :   1.7746794065471103

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.0028516356318156913 - 0.0016526027371530426im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.1528820610966118 - 0.26898563243703677im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.02044718924448773 - 0.03561262818985947im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.040851376122264776 + 0.06987243342195046im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206941472042736e-7 -2.5020146073002197e-7 … -4.880643392088655e-8 -2.5020151788210615e-7; -6.293615903504907e-7 -4.79303555213192e-7 … -1.2077841232876656e-7 -4.793035940261842e-7; … ; 1.3626256358845857e-7 1.1356226414079805e-7 … 4.711985836604992e-8 1.1356224843867984e-7; 1.3305470809367533e-7 1.3632638517169482e-7 … 5.200682477524527e-8 1.3632633940399915e-7;;; -2.5020149958843844e-7 -1.740490471171463e-7 … -3.556494490870236e-8 -1.7404908812188095e-7; -4.793035826586342e-7 -3.679731365456474e-7 … -1.1202892180294085e-7 -3.679731632892226e-7; … ; 1.1356225280766019e-7 1.1211094494030481e-7 … 1.0980677627842596e-7 1.1211093027138596e-7; 1.3632635415873497e-7 1.2550863828845395e-7 … 5.094905590808936e-8 1.255086036834377e-7;;; -4.880642599532495e-8 -3.556493332964046e-8 … -1.1412146188072954e-8 -3.556494210826313e-8; -1.2077840841768785e-7 -1.1202891512973813e-7 … -9.78622264789637e-8 -1.1202892015027998e-7; … ; 4.711986401821743e-8 1.0980679325625428e-7 … 2.4938222149679455e-7 1.0980678118968965e-7; 5.200683268310747e-8 5.094907086951679e-8 … 4.910624750689193e-8 5.0949059868594996e-8;;; … ;;; 6.475760105475906e-8 4.11033268434383e-8 … -6.681690500415645e-9 4.1103339983835915e-8; 1.6897558907451448e-7 6.653287416872969e-8 … -1.3294447536700713e-7 6.653288007687687e-8; … ; 3.2677605668600745e-8 1.6017365165333996e-7 … 3.847035800886122e-7 1.6017363742275295e-7; -9.141994025054822e-8 -3.583036458956688e-8 … 7.037282416512763e-8 -3.583036048021275e-8;;; -4.880641665047642e-8 -3.556492593820806e-8 … -1.1412144599390483e-8 -3.5564935635285315e-8; -1.207783998850517e-7 -1.120289058169418e-7 … -9.786222359292357e-8 -1.1202891425561904e-7; … ; 4.711986961030621e-8 1.0980680121528551e-7 … 2.4938223407972304e-7 1.098067888795373e-7; 5.200684099240775e-8 5.094907809253409e-8 … 4.91062540446657e-8 5.094906763539903e-8;;; -2.5020147888485436e-7 -1.740490312565895e-7 … -3.556494215640561e-8 -1.740490736590418e-7; -4.793035660033461e-7 -3.679731224265302e-7 … -1.1202892007534789e-7 -3.6797315238397537e-7; … ; 1.1356225953404351e-7 1.1211095077379629e-7 … 1.0980678073642996e-7 1.1211093629306972e-7; 1.3632637034858118e-7 1.2550865047562e-7 … 5.094906028519016e-8 1.2550861616481646e-7;;;;], δHtotψ = Matrix{ComplexF64}[[-0.0016587603976778442 + 0.0009612983054348767im -0.0019933177016583196 + 0.09527618028681722im 0.3875367220588255 + 0.5651231713069881im 1.7334917394615634 - 0.21438275360367243im; 0.09786750284621809 + 0.17246895626101233im 1.6217387478922407 + 0.12899052373670605im 0.08569053809617234 + 0.1251371778205066im 0.38423705282747456 - 0.045942880188426884im; … ; 0.024292642994850948 + 0.04257723143239283im -0.02068283398279997 - 0.00038074459410666444im -0.0032475783694832384 - 0.020585418936767165im -0.02453953830567146 + 0.004655602690361586im; -0.03652491040875863 - 0.06248197859647901im 0.025584192934696486 + 0.0005925958655514004im 0.06999926233560132 - 0.12228562586257002im 0.013627514004018047 + 0.020019985090577088im]], δVind = [0.0019654319269149815 0.0020810003900544613 … 0.0031536657786317567 0.002081000886864597; 0.003889130812944778 0.004779269107319041 … 0.010276672189206304 0.004779269648606469; … ; -0.006150847065807485 -0.00826960309336812 … -0.01730296974831866 -0.008269602261117042; 0.0002386825013748913 -0.0003986973995692496 … -0.003664693703412945 -0.00039869683188676155;;; 0.002081000729797381 0.0022232281569839695 … 0.0034201207944744073 0.0022232287130460435; 0.00477926948392742 0.00573495337592686 … 0.01090359902195648 0.005734953985482636; … ; -0.008269602559058188 -0.009601007854795935 … -0.012623110640106775 -0.009601007006838868; -0.0003986970186244511 -0.0010730129521913446 … -0.004230679444685211 -0.0010730123238006895;;; 0.0031536652142718866 0.0034201194886540904 … 0.0048325073395505 0.0034201204723804724; 0.01027667156649905 0.010903597619579927 … 0.01227214315481981 0.010903598694568114; … ; -0.017302970268373036 -0.012623111822318272 … -0.00962026400728231 -0.012623110900315199; -0.003664694325115219 -0.004230680716526147 … -0.005491107061038242 -0.004230679756539574;;; … ;;; -0.0023365598292188045 -0.00281546202134837 … 0.005693319659475317 -0.0028154628834782744; -0.012913091707773083 -0.017841177791503272 … 0.011498340582945316 -0.01784117889094065; … ; -0.02517447033729039 -0.012810887506163148 … -0.008395874564988113 -0.012810886614145898; 0.007597978664779531 0.010225413267026895 … -0.0054497330683267 0.010225411889663385;;; 0.0031536645680555336 0.0034201187718088796 … 0.0048325061864584445 0.0034201197493369936; 0.01027667089059418 0.010903596968561819 … 0.012272142284938479 0.01090359795002193; … ; -0.017302970889068132 -0.012623112500884007 … -0.009620264609681644 -0.012623111596847359; -0.0036646949871163292 -0.004230681423451751 … -0.005491107794861014 -0.004230680435560694;;; 0.002081000548229469 0.0022232279498234212 … 0.0034201204242018547 0.0022232285106829276; 0.004779269277796766 0.005734953149163762 … 0.010903598613010027 0.005734953756329191; … ; -0.008269602805839603 -0.009601008127687506 … -0.012623110965931453 -0.009601007271960092; -0.00039869721356898926 -0.001073013173639259 … -0.004230679765393959 -0.0010730125393129727;;;;], δρ0 = [-3.6021403032035583e-7 -2.488842773369866e-7 … -4.852836498348541e-8 -2.488842753754335e-7; -5.972249957658856e-7 -4.6069779731414863e-7 … -1.1905765911873405e-7 -4.606977929191414e-7; … ; 1.2735786366468338e-7 1.0836336468885849e-7 … 4.6407148171347275e-8 1.0836335905603595e-7; 1.035200323278221e-7 1.1957958593291017e-7 … 5.0695046331985034e-8 1.1957958328489327e-7;;; -2.488842780042195e-7 -1.731067287802703e-7 … -3.5356973841179835e-8 -1.7310672668425746e-7; -4.606977964603059e-7 -3.572170508193088e-7 … -1.1091437180586896e-7 -3.57217045811554e-7; … ; 1.0836335991621532e-7 1.0846493953117856e-7 … 1.0862233004679124e-7 1.0846493404162787e-7; 1.1957958264718628e-7 1.1611144663110167e-7 … 5.0208289924927015e-8 1.1611144415437078e-7;;; -4.852836655091287e-8 -3.5356975765144646e-8 … -1.1341310978992756e-8 -3.535697459815051e-8; -1.1905766210986774e-7 -1.1091437721627033e-7 … -9.770666627347213e-8 -1.1091437329505699e-7; … ; 4.640715142001069e-8 1.0862233669789622e-7 … 2.484524092330976e-7 1.0862233230949989e-7; 5.069504707944167e-8 5.0208292491600256e-8 … 4.926856999855595e-8 5.020829062834491e-8;;; … ;;; 6.434571411830215e-8 4.0838631215157775e-8 … -6.637925917881335e-9 4.0838631391182365e-8; 1.6975037594800423e-7 6.683086929639635e-8 … -1.3340479681045862e-7 6.683087012492519e-8; … ; 3.2709957612005626e-8 1.6023520864371073e-7 … 3.8390141030212704e-7 1.6023520426651418e-7; -9.27136804031912e-8 -3.633622050427653e-8 … 7.126307879696843e-8 -3.6336221450241564e-8;;; -4.852836503578407e-8 -3.535697432040218e-8 … -1.1341311381008688e-8 -3.535697395806738e-8; -1.1905765959643562e-7 -1.1091437351365568e-7 … -9.770666783114765e-8 -1.1091437274209812e-7; … ; 4.640715500318974e-8 1.0862234023825056e-7 … 2.4845241202924497e-7 1.0862233545885196e-7; 5.0695049604121174e-8 5.02082944472662e-8 … 4.9268571987397386e-8 5.020829306012116e-8;;; -2.4888427416584395e-7 -1.731067255718438e-7 … -3.535697397709692e-8 -1.731067245115723e-7; -4.6069779275611976e-7 -3.572170466909078e-7 … -1.1091437268781383e-7 -3.572170442196973e-7; … ; 1.0836336366303588e-7 1.0846494270615952e-7 … 1.0862233219585611e-7 1.0846493731890096e-7; 1.1957958669431026e-7 1.161114497442934e-7 … 5.020829161167215e-8 1.1611144743160318e-7;;;;], δeigenvalues = [[0.0004950062517053606, 0.09679860386726949, 0.047294193144899424, 0.01305271012460123]], δ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.15475965281846835 -0.15461620269434528 … -0.15418753099246665 -0.15461620269328913; -0.15461620269321474 -0.15447229889299605 … -0.15404229466888672 -0.15447229889193; … ; -0.1541875309947692 -0.1540422946723409 … -0.15360833734172466 -0.1540422946712846; -0.15461620269443807 -0.15447229889421665 … -0.1540422946701381 -0.15447229889316375;;; -0.154616202694051 -0.15447229889383293 … -0.154042294669737 -0.1544722988927708; -0.15447229889269742 -0.1543279430844975 … -0.15389660282179124 -0.1543279430834257; … ; -0.1540422946720493 -0.15389660282526013 … -0.15346128700565906 -0.15389660282419743; -0.1544722988939247 -0.1543279430857222 … -0.1538966028230462 -0.15432794308466313;;; -0.1541875309939695 -0.15404229467154804 … -0.1536083373408963 -0.15404229467047306; -0.15404229467040068 -0.15389660282361722 … -0.15346128700396805 -0.1538966028225327; … ; -0.1536083373432324 -0.1534612870074729 … -0.15302191313054483 -0.15346128700639716; -0.15404229467163885 -0.15389660282485307 … -0.15346128700523412 -0.15389660282378104;;; … ;;; -0.15347880584856657 -0.15333139195895984 … -0.15289092492034687 -0.15333139195788112; -0.15333139195780224 -0.15318350573623402 … -0.15274162877143335 -0.15318350573514436; … ; -0.15289092492270487 -0.1527416287749714 … -0.1522955395847656 -0.15274162877389452; -0.15333139195905768 -0.1531835057374861 … -0.15274162877272054 -0.15318350573641162;;; -0.1541875309930171 -0.15404229467058755 … -0.1536083373399423 -0.15404229466952274; -0.1540422946694457 -0.15389660282265458 … -0.1534612870030104 -0.15389660282157944; … ; -0.15360833734226792 -0.15346128700649925 … -0.15302191312958058 -0.15346128700543588; -0.15404229467068326 -0.15389660282388892 … -0.15346128700427794 -0.1538966028228281;;; -0.15461620269357162 -0.15447229889334968 … -0.15404229466925692 -0.15447229889229247; -0.15447229889221697 -0.1543279430840132 … -0.15389660282130938 -0.15432794308294598; … ; -0.15404229467156397 -0.15389660282477022 … -0.1534612870051738 -0.1538966028237137; -0.15447229889344366 -0.154327943085237 … -0.153896602822565 -0.15432794308418343]), 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.015163897733979144 -0.01255012919116992 … -0.004966868022397866 -0.01255012915087498; -0.012550129167807556 -0.010362595118052709 … -0.004135158840695981 -0.010362595088471024; … ; -0.004966868042291504 -0.00413515887498641 … -0.0020697798254417934 -0.004135158754981709; -0.012550129161889807 -0.010362595119845015 … -0.00413515874070373 -0.010362595055976085;;; -0.012550129158250618 -0.010362595116096336 … -0.004135158775908343 -0.010362595062784804; -0.010362595101657955 -0.008622039052195836 … -0.0038975112380040706 -0.00862203900357188; … ; -0.004135158776854409 -0.0038975112475901556 … -0.0034841066816820422 -0.0038975111351723358; -0.010362595067086713 -0.008622039020697891 … -0.0038975111334502823 -0.008622038949156727;;; -0.004966868056570746 -0.004135158877645602 … -0.0020697799413963935 -0.004135158803961475; -0.00413515888049202 -0.003897511335759501 … -0.003484106812405241 -0.003897511259285699; … ; -0.0020697799173172017 -0.0034841068079063163 … -0.005620579829871721 -0.003484106726813399; -0.004135158795828976 -0.0038975112546766454 … -0.0034841067399836332 -0.003897511173556436;;; … ;;; -0.006547898295396584 -0.004876340406828347 … -0.0014866437747180243 -0.004876340430259056; -0.004876340403056667 -0.002677989092229679 … -0.004123054136925888 -0.0026779891194706248; … ; -0.00148664369451171 -0.004123054090243931 … -0.007170895126613964 -0.0041230540295997585; -0.004876340451639418 -0.002677989158108266 … -0.004123054057508013 -0.0026779892101119057;;; -0.004966868083641133 -0.004135158878133971 … -0.0020697800372662 -0.004135158857309361; -0.004135158855846666 -0.0038975112809380734 … -0.003484106867436813 -0.003897511273317598; … ; -0.002069780003005398 -0.003484106853907128 … -0.005620579834797283 -0.003484106766953377; -0.004135158850118937 -0.0038975112830560155 … -0.0034841067829269193 -0.003897511230235842;;; -0.01255012917780277 -0.010362595120939833 … -0.004135158833500095 -0.010362595095016112; -0.010362595094922121 -0.008622039026284157 … -0.0038975112700735766 -0.008622039012972597; … ; -0.004135158847289053 -0.003897511301662114 … -0.003484106713169598 -0.003897511193420811; -0.01036259510439122 -0.008622039045173464 … -0.003897511185368279 -0.008622038992275774])], 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.009532892483592882 -0.006919175311707312 … 0.0006640013478138116 -0.0069191752703561885; -0.006919175287214269 -0.004731678603360934 … 0.0014956971563970843 -0.004731678572713213; … ; 0.0006640013256177099 0.0014956971186523907 … 0.003561046785397369 0.0014956972397133853; -0.006919175282519792 -0.004731678606373819 … 0.0014956972551379748 -0.004731678541451982;;; -0.006919175278493745 -0.004731678602241614 … 0.0014956972203342914 -0.004731678547867931; -0.004731678586667558 -0.002991147634035636 … 0.0017333330837370313 -0.002991147584339899; … ; 0.0014956972170760193 0.001733333070681999 … 0.002146715763185275 0.0017333331841625243; -0.004731678553323557 -0.002991147603762378 … 0.0017333331870359346 -0.002991147531162145;;; 0.0006640013121380507 0.0014956971167860086 … 0.003561046670271134 0.0014956971915451086; 0.0014956971150870335 0.0017333329841556448 … 0.0021467156341531124 0.0017333330617139682; … ; 0.00356104669201425 0.0021467156351471185 … 1.0229458538810708e-5 0.0021467157173157866; 0.0014956971985119845 0.0017333330640026557 … 0.002146715705308705 0.0017333331461948966;;; … ;;; -0.0009170702183934811 0.0007544782140201742 … 0.004144158597581296 0.0007544781916681859; 0.0007544782189495116 0.0029528243751095282 … 0.001507745480783182 0.0029528243489582388; … ; 0.004144158675429636 0.0015077455239270246 … -0.0015401021888475763 0.001507745585648086; 0.0007544781691113476 0.002952824307978915 … 0.0015077455589139477 0.0029528242570497493;;; 0.0006640012860200679 0.0014956971172580937 … 0.0035610465753553146 0.001495697139147518; 0.001495697140687374 0.0017333330399397192 … 0.002146715580079191 0.0017333330486353346; … ; 0.0035610466072905048 0.0021467155901199724 … 1.022945457747771e-5 0.002146715678137095; 0.0014956971451775929 0.001733333036587431 … 0.0021467156633215986 0.0017333330904684228;;; -0.00691917529756653 -0.004731678606601859 … 0.001495697163222627 -0.0047316785796208995; -0.004731678579451275 -0.00299114760763965 … 0.0017333330521494177 -0.00299114759326086; … ; 0.0014956971471266813 0.0017333330170999544 … 0.0021467157321829423 0.0017333331263977734; -0.0047316785901470315 -0.0029911476277527556 … 0.0017333331355991362 -0.002991147573801492]), 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.9441204607791507e-6 + 1.4440071613458788e-6im -0.004097846938932913 + 0.006980966785927305im … -0.0028839746899847105 + 0.00480351300021486im 0.004191176703209178 - 0.007026197131220047im; 0.009363583430952234 + 0.016268545397133828im -0.006112744148666178 - 0.018828333362709157im … 0.007178899300744676 + 0.005665824899624914im -0.010969491473819062 - 0.010624422724051516im; … ; 0.011252204708052351 + 0.019480726573283962im -0.010670439930499888 - 0.01564588635974328im … 0.005398472288400378 + 0.01258642686473995im -0.00900313194792414 - 0.01846250250568989im; -0.00944389396057367 - 0.016183089335906266im 0.010951268275727382 + 0.010610439058779838im … -0.0032125666849193673 - 0.012261581791767982im 0.006093742917848069 + 0.01881752202002874im;;; 0.0018281506787459919 - 0.0024752364988986433im 0.0013224475729399596 - 0.002637896178545142im … 0.002355232266557734 - 0.0037971430126079657im -0.003533050370131636 + 0.005567371843994083im; -0.009588890954574266 - 0.013288399546669773im 0.0064495981252257195 + 0.014311086539039772im … -0.006160204325549631 - 0.0047105786968887915im 0.009575601885926583 + 0.009030303803658064im; … ; -0.009474192768519907 - 0.017546041141587226im 0.008836826320273801 + 0.01397630018574454im … -0.004572324492539699 - 0.011147180635766769im 0.00752839132405582 + 0.016186897154267993im; 0.0074733483168885366 + 0.01614922983484473im -0.008162438792594299 - 0.01089722792741474im … 0.0026092020229893445 + 0.01041787556990222im -0.005036045769548914 - 0.01617983557543556im;;; -0.001274684650366688 + 0.0017046634706220957im 5.839727911459856e-5 + 0.000230972991191849im … -0.001381366597232076 + 0.002149257801547149im 0.0017245164360907932 - 0.0025843545721814043im; 0.0060761530728026345 + 0.007792959612565943im -0.004400504188556893 - 0.007789939163435079im … 0.003967217441779574 + 0.0030367730535269923im -0.005708519813737615 - 0.005580332128145769im; … ; 0.0059264385012730274 + 0.011562514469802578im -0.005554070983114421 - 0.009432011986273334im … 0.0030919471092148185 + 0.007839599957750694im -0.004802077844424608 - 0.010702412068864992im; -0.0043146481295027106 - 0.010132024283194086im 0.004368536536600524 + 0.007337034412484041im … -0.001643829045931414 - 0.006642213230284293im 0.0030605066796049385 + 0.009547181915893614im;;; … ;;; -0.0004978459862840586 + 0.000706058463939724im 0.000743074193627654 - 0.0011473755401200337im … 0.00018888279176393642 - 0.0003774822158001161im 8.749018138747957e-5 - 4.0570436617669685e-5im; -0.002039489974037777 - 0.004830074184530957im 0.0015359193549157723 + 0.004739729130491371im … -0.001742281194684708 - 0.00236471509494492im 0.0021048886728265593 + 0.0037794371298241923im; … ; -0.0037294563994898482 - 0.0055797080088336im 0.0035608206701431442 + 0.004704398028886175im … -0.0022341033564586043 - 0.004348333150344447im 0.0031605375973710846 + 0.0053803144163187796im; 0.0028720823098060387 + 0.003640318440494189im -0.002812405396076082 - 0.0027592669974508275im … 0.001395669076681687 + 0.003046732230164273im -0.002243553459574361 - 0.003719766102390563im;;; 0.0011665675885924379 - 0.0016119461808429404im -0.001622689755111201 + 0.002511818343277894im … -0.000573173197938644 + 0.001052619237949496im 4.543616831139737e-5 - 0.0003050745458566558im; 0.004397241226769686 + 0.010073877763661223im -0.003145314031667012 - 0.009493801891262551im … 0.003290033873706334 + 0.00395610418951628im -0.004454523356942948 - 0.007283080873435965im; … ; 0.00673679916484176 + 0.010465490519339538im -0.0062333064237512815 - 0.008512188987166163im … 0.003534446414576972 + 0.00723939452508791im -0.005480776449922147 - 0.009783035634361385im; -0.005996568054553939 - 0.007852399004240194im 0.005625791785966106 + 0.0056334996630908645im … -0.002321682608940884 - 0.005722825708277186im 0.004316441411661613 + 0.007844826809166724im;;; -0.00174400884297568 + 0.0023817377364958843im 0.0034240473058821724 - 0.0054893335752118485im … 0.00163205760892989 - 0.002816966755934138im -0.001434663551350308 + 0.002717998678435328im; -0.007490277812642783 - 0.016140734011925957im 0.0050944804032909095 + 0.01614971420346121im … -0.005592693677141374 - 0.005480162384707355im 0.00822228138327515 + 0.01086736736262044im; … ; -0.01020034294962413 - 0.016565532337241277im 0.009406680040548277 + 0.01320608279092933im … -0.004898946778745361 - 0.010704537031696736im 0.008097600684059902 + 0.01541725865894997im; 0.009576982545616052 + 0.013300150451793277im -0.009520528760676498 - 0.009059653440145955im … 0.003177267488517116 + 0.009648262800296731im -0.006392204169156965 - 0.014343854885107827im],)])]), [8.519873468495262e-7 4.627264911361701e-7 … 2.4177775920395706e-8 4.627264863536821e-7; 4.627264883633578e-7 2.501447507051685e-7 … 1.3582992894629652e-8 2.501447484172106e-7; … ; 2.4177776225644933e-8 1.3582993248388833e-8 … 1.5659205365565983e-9 1.358299201034782e-8; 4.627264876609957e-7 2.5014475084379407e-7 … 1.358299186305046e-8 2.5014474590393704e-7;;; 4.6272648722906643e-7 2.50144750553844e-7 … 1.3582992226240802e-8 2.501447464305553e-7; 2.501447494371414e-7 1.3895607424842393e-7 … 1.1279985908507688e-8 1.3895607174846728e-7; … ; 1.3582992235998825e-8 1.1279985995541359e-8 … 7.937751059042036e-9 1.1279984974835688e-8; 2.5014474676327923e-7 1.3895607262899466e-7 … 1.1279984959197704e-8 1.38956068950769e-7;;; 2.4177776444738707e-8 1.3582993275822135e-8 … 1.5659208088442074e-9 1.3582992515653697e-8; 1.3582993305189714e-8 1.12799867960843e-8 … 7.937751991448994e-9 1.127998610173384e-8; … ; 1.565920752300321e-9 7.937751959361263e-9 … 3.571855420777192e-8 7.937751380950057e-9; 1.3582992431754108e-8 1.1279986059887309e-8 … 7.937751474889013e-9 1.1279985323345455e-8;;; … ;;; 5.7918434759211243e-8 2.2816032569737264e-8 … 5.620365425548152e-10 2.2816032915179733e-8; 2.28160325141349e-8 3.489280502798011e-9 … 1.3458507893172109e-8 3.4892806134113232e-9; … ; 5.620364488785719e-10 1.3458507414612766e-8 … 7.728314416101908e-8 1.3458506792917074e-8; 2.2816033230390647e-8 3.489280770305511e-9 … 1.3458507079018661e-8 3.4892809814724266e-9;;; 2.4177776860106852e-8 1.3582993280863055e-8 … 1.5659210339699216e-9 1.3582993066019757e-8; 1.3582993050931771e-8 1.1279986298327556e-8 … 7.937752383970365e-9 1.1279986229136304e-8; … ; 1.5659209535158212e-9 7.937752287469051e-9 … 3.571855430669016e-8 7.937751667254077e-9; 1.3582992991843272e-8 1.1279986317560918e-8 … 7.937751781188594e-9 1.1279985837972836e-8;;; 4.627264895496682e-7 2.501447509284634e-7 … 1.3582992820391121e-8 2.5014474892342797e-7; 2.501447489161597e-7 1.3895607291620474e-7 … 1.1279986199683821e-8 1.389560722317967e-7; … ; 1.358299296264518e-8 1.1279986486496614e-8 … 7.937751283633085e-9 1.1279985503708578e-8; 2.5014474964853534e-7 1.389560738873779e-7 … 1.1279985430596325e-8 1.389560711676955e-7;;;;], Matrix{ComplexF64}[[-0.22539792064467384 + 0.13062445309755502im -0.020186406793019328 + 0.964868966617586im -3.2118272463328258e-6 - 2.1839026393936597e-6im -4.717161547272393e-6 - 1.428924313531833e-6im; 0.15292773427670492 - 0.08862593576997532im -0.0008894943600977204 + 0.042515956786397545im 0.17838526933599552 - 0.1226434531910793im -0.06971432246349486 - 0.5468710067515388im; … ; -0.048586283847210074 + 0.02815712202659541im 0.0003217176053218031 - 0.01537745325223806im -0.002906158185619162 - 0.005758644125289843im -0.00010941701259026418 - 0.010088774023065903im; 0.08711553970925859 - 0.05048591262984641im -0.0006373474475292417 + 0.030463929125907528im 0.025780818965715872 + 0.0254051514793816im -0.002831604407551748 + 0.026405388998122437im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239441093319806, [[-0.5541104230270715, -0.010678398839324664, 0.14495661918835592, 0.144957568496499]], 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.6206941472042736e-7, -6.293615903504907e-7, -2.2951204892843918e-7, 4.994026351440958e-7, 7.923312898961975e-7, 5.851073400770467e-7, 6.516972249739351e-7, 1.0189309802481383e-6, 3.6385926752019413e-7, -2.2010932985229395e-6  …  2.781018064528627e-6, 1.2016954017724932e-6, 2.4456287079177784e-8, 1.0078534889580351e-7, 5.535306647174871e-7, 4.534171564015417e-7, 4.631842275701132e-9, -1.0437056418286488e-7, 1.1211093629306972e-7, 1.2550861616481646e-7], resid_history = [0.24939209475943971, 0.003766548014548187, 0.0002852313016831659, 4.681767248765144e-6, 3.07864769702868e-8, 5.026268575858442e-11, 4.2682040731651864e-8, 1.1998835765507805e-9], converged = true, n_iter = 8, residual_norm = 1.1998835765507805e-9, maxiter = 100, tol = 1.0e-8, s = 0.9354576109319984, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.4021004224864117e-7, -4.155150258192966e-8, 1.1998835765507805e-9, -0.00020442967070176765, 4.686583249431763e-6, -3.110648750248981e-8, 5.026268575858442e-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}}([[-9.977897561276084e-7, -2.4263216516193884e-6, -7.999888308159799e-7, 1.341097689820894e-6, 1.6974800845109415e-6, 1.2008666141666502e-6, 1.4080234522448565e-6, 2.107648399922971e-6, 6.168652089861627e-7, -3.0443721399618867e-6  …  3.2557402672452284e-6, 1.2865232545078923e-6, 2.983983527766485e-8, 1.4277584988036175e-7, 7.628850208256957e-7, 5.100404453330553e-7, 4.350149112926988e-9, -1.0845025588392483e-7, 1.545075466967414e-7, 1.6758835125700227e-7], [1.6745396818773424e-6, 3.2006466528895283e-6, 1.0315443391106483e-6, -1.7424420604305589e-6, -2.21388729329466e-6, -1.543632763950271e-6, -1.782643955946e-6, -2.667817215197563e-6, -7.905967182431872e-7, 3.927943093280577e-6  …  -3.336280160871001e-6, -1.2822317960075876e-6, -2.9604208708718937e-8, -1.4274839174688924e-7, -7.534331266651586e-7, -4.739412340271882e-7, -3.5609685136287716e-9, 7.931074078818808e-8, -9.896621708652669e-8, 6.127227420458639e-8]]), H = [1.1439287031517402 0.013417653285527261 … 0.0 0.0; 0.1264013648486794 1.0341259344759215 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1508910386869953 0.12691363765597094 … 0.0 0.0; 0.0 1.0268021394762776 … 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.9257125304498006
Interacting polarizability:     1.7736548607888876

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, δρ))")

[ 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 = 6.22e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.02e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 4.17e-12
└ * number of operations = 11
Non-interacting polarizability: 1.925712530449709
Interacting polarizability:     1.7736548569644124

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