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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457372412583617
Displaced dipole:  0.017612221445940414
Polarizability :   1.7746795170066247

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.002088861841218813 + 0.0025494268928271955im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.2376681465705623 + 0.19809100978728117im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.03165937336585814 + 0.026153984909799713im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.06280301293584839 - 0.051056571279886054im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620694686950157e-7 -2.502015208988996e-7 … -4.880643992424451e-8 -2.502015398988057e-7; -6.293617692318201e-7 -4.793037033520796e-7 … -1.207784562166784e-7 -4.793037158831241e-7; … ; 1.3626262875054775e-7 1.1356231520023552e-7 … 4.711990421620507e-8 1.1356230944181477e-7; 1.3305480779000657e-7 1.3632644471987943e-7 … 5.2006856910589215e-8 1.3632642773018982e-7;;; -2.502015572598869e-7 -1.7404910438117997e-7 … -3.556495518839889e-8 -1.7404911901146973e-7; -4.793037276732396e-7 -3.679732566714281e-7 … -1.1202896973989965e-7 -3.6797326847082504e-7; … ; 1.1356230919493716e-7 1.1211099668919404e-7 … 1.0980684149112248e-7 1.1211099103470271e-7; 1.3632641681223045e-7 1.2550867637383276e-7 … 5.0949086116347416e-8 1.2550866320198907e-7;;; -4.880644812486106e-8 -3.556495480226647e-8 … -1.1412156389901314e-8 -3.556495796316286e-8; -1.2077846104602703e-7 -1.1202896696084244e-7 … -9.786227736999894e-8 -1.1202897190705945e-7; … ; 4.711990777182588e-8 1.0980684700301303e-7 … 2.493823167462002e-7 1.0980684238643505e-7; 5.2006853379851846e-8 5.09490888737335e-8 … 4.910627530454809e-8 5.0949084125491116e-8;;; … ;;; 6.475759024596672e-8 4.110332136772168e-8 … -6.6816955356925395e-9 4.110332639204923e-8; 1.6897558799016164e-7 6.653285844855848e-8 … -1.329445070828391e-7 6.653286543724049e-8; … ; 3.267764945955621e-8 1.6017372123924915e-7 … 3.847037312390639e-7 1.601737205931883e-7; -9.141995301873431e-8 -3.583035811278965e-8 … 7.037286932436873e-8 -3.5830356554906337e-8;;; -4.88064231669429e-8 -3.5564936887013896e-8 … -1.1412148792288464e-8 -3.556493836631615e-8; -1.207784434719159e-7 -1.1202895511599032e-7 … -9.786225743372043e-8 -1.1202895125282321e-7; … ; 4.7119910554733986e-8 1.0980685635613559e-7 … 2.493823448622633e-7 1.0980685413584667e-7; 5.200687093889147e-8 5.094910466505507e-8 … 4.910628807769353e-8 5.094910006929876e-8;;; -2.502015032888123e-7 -1.7404906583173603e-7 … -3.556494474269358e-8 -1.7404907823546348e-7; -4.793036907966964e-7 -3.679732329899083e-7 … -1.1202895877113916e-7 -3.6797323818827226e-7; … ; 1.135623164477427e-7 1.1211100251788523e-7 … 1.0980684678090679e-7 1.1211099788916805e-7; 1.3632645573289474e-7 1.255087049982756e-7 … 5.094909304565663e-8 1.255086920752197e-7;;;;], δHψ = Matrix{ComplexF64}[[0.0012150645084975482 - 0.0014829693617182723im 0.0733870956194369 + 0.06079358903551314im 0.6730429065325509 + 0.5082183807329im -1.9360142749231395 + 0.2764663087207388im; -0.15224095080731995 - 0.12706802120644567im 1.1092530667394596 - 1.1900558223809856im 0.14897343077664588 + 0.11233354971732014im -0.4260749275859421 + 0.061771719186566834im; … ; -0.03770732542599567 - 0.03132268766907227im -0.013156583865887753 + 0.01596345067790439im -0.016326667643758223 - 0.013526989366169199im 0.023720381950772362 - 0.004721845442535805im; 0.056155064074073266 + 0.04565824676263983im 0.01636958934347574 - 0.019670705234312746im -0.046386524165142495 - 0.08521970341790969im -0.057947187947431734 - 0.02337511877411427im]], δVind = [0.0019654321487573034 0.0020810007859916797 … 0.0031536654995419004 0.0020810009684168484; 0.003889132015711743 0.0047792705521961165 … 0.01027667296080088 0.004779270761561474; … ; -0.006150849655515805 -0.008269605323623377 … -0.01730297038667809 -0.008269605046595935; 0.00023868173728396449 -0.00039869802653323076 … -0.0036646953038651756 -0.0003986978324303219;;; 0.0020810011054421193 0.002223228739867735 … 0.00342012065404958 0.002223228946015855; 0.004779270923130625 0.005734955092105082 … 0.010903599832980941 0.005734955311162152; … ; -0.008269604707053324 -0.00960100948587641 … -0.012623112499341723 -0.009601009207210989; -0.0003986976745804457 -0.0010730134238156328 … -0.004230680623847601 -0.0010730132045397103;;; 0.003153665974292688 0.0034201205638017 … 0.004832506652879163 0.003420120925239251; 0.010276673553929272 0.010903599854315226 … 0.0122721437362461 0.010903600134952793; … ; -0.01730296887083415 -0.012623112375370121 … -0.009620265738248469 -0.012623112072012332; -0.003664694707160558 -0.004230680644947642 … -0.005491107649008199 -0.004230680335244655;;; … ;;; -0.0023365596607211633 -0.0028154621804148076 … 0.005693315820944901 -0.0028154624432615814; -0.012913093376423435 -0.017841182025692102 … 0.011498340255840917 -0.017841181921022672; … ; -0.025174465850774253 -0.01281088941912854 … -0.008395877013929826 -0.012810888983929965; 0.007597980823400806 0.010225416729062842 … -0.005449734920570078 0.010225416301869782;;; 0.0031536643661921587 0.0034201187612226676 … 0.0048325041462748105 0.0034201191483825225; 0.01027667169877433 0.010903597775340524 … 0.012272141993920325 0.010903598384180337; … ; -0.017302971916146416 -0.012623114467973954 … -0.009620267194091292 -0.012623114077866019; -0.003664696489668799 -0.004230682474908237 … -0.005491109672178055 -0.004230682200542361;;; 0.0020810006485092137 0.0022232282264382313 … 0.0034201197863470306 0.002223228429609753; 0.004779270391734888 0.005734954473455172 … 0.010903598981982959 0.0057349547281248065; … ; -0.008269605629197768 -0.009601010460611934 … -0.012623113547575547 -0.009601010161114115; -0.00039869818315723554 -0.0010730139870035623 … -0.00423068160952241 -0.0010730137825993593;;;;], δρ0 = [-3.60214076111616e-7 -2.4888431161820706e-7 … -4.852837861044341e-8 -2.4888431122990196e-7; -5.972250839357239e-7 -4.606978673500606e-7 … -1.1905769645547547e-7 -4.6069786579015884e-7; … ; 1.273578932908109e-7 1.0836339337117459e-7 … 4.640718367570945e-8 1.0836339102801017e-7; 1.0352005285875781e-7 1.1957960462262e-7 … 5.0695060461557826e-8 1.195796021189614e-7;;; -2.488843118135314e-7 -1.7310675410434935e-7 … -3.535698683182483e-8 -1.7310675468828722e-7; -4.6069786539157843e-7 -3.572171057114221e-7 … -1.1091441155965132e-7 -3.5721710709506405e-7; … ; 1.0836339349133749e-7 1.0846497446036506e-7 … 1.0862237538185194e-7 1.0846497189574461e-7; 1.1957960336693306e-7 1.1611146609790443e-7 … 5.02083069176198e-8 1.1611146361572875e-7;;; -4.852837917285732e-8 -3.535698680672051e-8 … -1.1341321079286344e-8 -3.53569867658655e-8; -1.1905769561492313e-7 -1.1091440913160291e-7 … -9.770670908001655e-8 -1.1091441115599172e-7; … ; 4.6407189059507596e-8 1.0862238046865342e-7 … 2.484524697268952e-7 1.0862237836971479e-7; 5.069506257509457e-8 5.020830925212638e-8 … 4.926858991978212e-8 5.020830749022824e-8;;; … ;;; 6.434571860134366e-8 4.08386311014897e-8 … -6.6379338322169495e-9 4.0838631930193475e-8; 1.697503711512274e-7 6.683085050492121e-8 … -1.3340483358195422e-7 6.683085587056792e-8; … ; 3.270999552133314e-8 1.6023525292301595e-7 … 3.839014785561481e-7 1.602352554651916e-7; -9.271366977568574e-8 -3.63362065372278e-8 … 7.126310034650237e-8 -3.633620661030391e-8;;; -4.852837807880021e-8 -3.535698679192054e-8 … -1.1341319519775628e-8 -3.5356985005736e-8; -1.1905769583806238e-7 -1.1091441350905034e-7 … -9.770670267995784e-8 -1.1091440671243379e-7; … ; 4.6407185866883266e-8 1.086223763230492e-7 … 2.484524690357025e-7 1.086223765626443e-7; 5.0695062290034806e-8 5.020830879611193e-8 … 4.926858907927275e-8 5.020830718254934e-8;;; -2.4888431050655027e-7 -1.7310675430794624e-7 … -3.535698533312359e-8 -1.731067524852777e-7; -4.60697866865417e-7 -3.57217110613036e-7 … -1.1091440871445004e-7 -3.5721710536117517e-7; … ; 1.0836339189427779e-7 1.0846497224203956e-7 … 1.08622373842347e-7 1.0846497064544823e-7; 1.1957960366048903e-7 1.1611146597131382e-7 … 5.0208305688217965e-8 1.1611146374483468e-7;;;;], δeigenvalues = [[0.0004950053638938527, 0.09679861677934747, 0.04483541595504245, 0.010279956695096073]], δ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.154759652815687 -0.15461620269098098 … -0.15418753099084556 -0.15461620269109222; -0.154616202690968 -0.15447229889016306 … -0.15404229466781302 -0.15447229889027284; … ; -0.15418753099090257 -0.1540422946678836 … -0.1536083373390002 -0.15404229466799815; -0.15461620269111565 -0.1544722988903095 … -0.15404229466796468 -0.1544722988904222;;; -0.1546162026913279 -0.15447229889052322 … -0.15404229466817876 -0.15447229889063457; -0.1544722988905097 -0.15432794308171996 … -0.15389660282078124 -0.15432794308182993; … ; -0.15404229466823702 -0.15389660282085355 … -0.1534612870029957 -0.15389660282096873; -0.15447229889065858 -0.15432794308186767 … -0.15389660282093462 -0.15432794308198075;;; -0.15418753099132754 -0.15404229466831104 … -0.15360833733942864 -0.15404229466842337; -0.15404229466829686 -0.1538966028209156 … -0.15346128700305245 -0.1538966028210267; … ; -0.15360833733948837 -0.15346128700312686 … -0.15302191312796273 -0.15346128700324338; -0.15404229466844813 -0.1538966028210658 … -0.1534612870032085 -0.15389660282117998;;; … ;;; -0.1534788058457561 -0.15333139195555057 … -0.15289092491870854 -0.1533313919556632; -0.15333139195553872 -0.15318350573336845 … -0.15274162877035927 -0.15318350573348027; … ; -0.15289092491876358 -0.1527416287704261 … -0.15229553958198305 -0.15274162877054268; -0.15333139195568551 -0.1531835057335143 … -0.15274162877051 -0.15318350573362857;;; -0.15418753099019544 -0.15404229466717587 … -0.15360833733828153 -0.15404229466728742; -0.1540422946671637 -0.1538966028197792 … -0.1534612870019045 -0.15389660281988982; … ; -0.15360833733833693 -0.15346128700197245 … -0.15302191312679525 -0.15346128700208755; -0.15404229466730998 -0.1538966028199247 … -0.15346128700205486 -0.15389660282003767;;; -0.15461620269075815 -0.15447229888995181 … -0.15404229466760164 -0.15447229889006311; -0.1544722988899392 -0.15432794308114778 … -0.1538966028202036 -0.15432794308125783; … ; -0.15404229466765773 -0.15389660282027273 … -0.15346128700240835 -0.15389660282038722; -0.154472298890086 -0.1543279430812935 … -0.1538966028203543 -0.15432794308140613]), 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.015163898024801514 -0.012550129515990018 … -0.004966868512025692 -0.012550129455777145; -0.012550129494408578 -0.010362595475213846 … -0.00413515941100376 -0.01036259543942143; … ; -0.00496686846306934 -0.004135159343871024 … -0.0020697807437770085 -0.0041351592950367826; -0.012550129455613991 -0.010362595451353485 … -0.004135159309537703 -0.010362595384851889;;; -0.012550129480733203 -0.01036259546544116 … -0.00413515941333424 -0.010362595429600708; -0.01036259543269452 -0.008622039392815821 … -0.003897511893439366 -0.00862203939831015; … ; -0.004135159349905174 -0.003897511839819125 … -0.0034841074464097115 -0.0038975117894947744; -0.010362595428499606 -0.008622039422332192 … -0.003897511819997015 -0.008622039364000689;;; -0.004966868555160242 -0.004135159442943678 … -0.0020697809548601933 -0.004135159412176236; -0.004135159390613987 -0.0038975118455438287 … -0.003484107587457542 -0.0038975118841904017; … ; -0.002069780891266152 -0.0034841075408017378 … -0.0056205804502080245 -0.0034841075046228727; -0.004135159406073239 -0.0038975119108768516 … -0.0034841075153439482 -0.00389751184229216;;; … ;;; -0.006547898000628943 -0.00487633997043913 … -0.0014866447832769538 -0.004876340033514658; -0.0048763400606764 -0.0026779884707579336 … -0.004123054641193299 -0.0026779886106808115; … ; -0.0014866447548405677 -0.004123054683067031 … -0.007170895690584468 -0.004123054728586211; -0.004876339989583334 -0.002677988441031256 … -0.004123054762861765 -0.002677988447174822;;; -0.004966868541927505 -0.004135159479648216 … -0.002069780808373517 -0.004135159344075694; -0.004135159415420007 -0.0038975119590057487 … -0.0034841074630518054 -0.00389751180741661; … ; -0.002069780800703546 -0.0034841074585017936 … -0.005620580438596069 -0.0034841074668415933; -0.004135159375463112 -0.0038975118816070415 … -0.003484107460594397 -0.0038975118070472903;;; -0.012550129486662662 -0.010362595496527422 … -0.004135159332924599 -0.010362595407347306; -0.010362595473681634 -0.008622039478837065 … -0.0038975118411239463 -0.008622039395071059; … ; -0.0041351593117825765 -0.0038975117910890174 … -0.0034841074141113496 -0.0038975117603541817; -0.010362595414212705 -0.008622039411965739 … -0.0038975117483020854 -0.008622039337619994])], 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.009532892771633893 -0.0069191756331631014 … 0.0006640008598070776 -0.006919175573061445; -0.006919175611568561 -0.004731678957689087 … 0.0014956965871630017 -0.004731678922006445; … ; 0.0006640009087065029 0.001495696654225072 … 0.0035610458697866135 0.0014956967029447661; -0.00691917557292155 -0.004731678933975136 … 0.00149569668847743 -0.004731678867586228;;; -0.006919175598253231 -0.004731678948276724 … 0.0014956965844666473 -0.0047316789125475996; -0.0047316789155164 -0.002991147971878093 … 0.0017333324293117335 -0.0029911479774823897; … ; 0.0014956966478375383 0.0017333324828596158 … 0.0021467150011209753 0.0017333325330687807; -0.00473167891147034 -0.002991148001542151 … 0.0017333325026007906 -0.0029911479433237242;;; 0.0006640008161905245 0.001495696554724927 … 0.0035610456582749936 0.0014956965853800417; 0.0014956966070688842 0.0017333324770729339 … 0.0021467148600164123 0.0017333324383152554; … ; 0.003561045721809333 0.0021467149065977483 … 1.0228840784608749e-5 0.0021467149426600955; 0.0014956965914584477 0.0017333324115897256 … 0.0021467149319740197 0.001733332480060231;;; … ;;; -0.0009170699208153662 0.0007544786538186644 … 0.004144157590660695 0.000754478590630504; 0.0007544785635933008 0.0029528249994468427 … 0.0015077449775898564 0.0029528248594121375; … ; 0.00414415761904207 0.0015077449356492335 … -0.0015401027500355283 0.0015077448900134805; 0.0007544786345395956 0.0029528250290277204 … 0.001507744855770733 0.0029528250227698846;;; 0.0006640008305553559 0.0014956965191555358 … 0.0035610458059087803 0.001495696654616508; 0.0014956965833960133 0.0017333323647474104 … 0.0021467149855700918 0.0017333325162259433; … ; 0.003561045813523351 0.0021467149900521025 … 1.0228853564019091e-5 0.0021467149815972003; 0.0014956966232066915 0.001733332442000623 … 0.0021467149878772037 0.001733332516447409;;; -0.006919175603612951 -0.004731678978791582 … 0.0014956966654534106 -0.004731678889722738; -0.004731678955933027 -0.002991148057327128 … 0.00173333248220483 -0.0029911479736711728; … ; 0.0014956966865393943 0.0017333325321705366 … 0.002146715034006645 0.0017333325627908805; -0.0047316788966108685 -0.0029911479906015188 … 0.0017333325748760337 -0.002991147916368406]), 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.9503006165922754e-6 - 2.1791069514107254e-6im -0.0013978902108993459 + 0.0003803271251316403im … -0.0008931144926521528 + 0.0001927550441858896im 0.0013318019794348662 - 0.00031550851840468543im; -0.011047799452102203 - 0.00908940676727455im 0.011962181814579944 + 0.008924517547002362im … -0.0067154424233821005 - 0.006130492972372052im 0.010360149892251869 + 0.009328666295423428im; … ; -0.015068290550402645 - 0.012413632758853642im 0.013340255901582842 + 0.011289423828355537im … -0.00934341024460536 - 0.007384579177916196im 0.013891436582517904 + 0.011151893621127378im; 0.011132225092441191 + 0.009082129733043603im -0.010362049834184456 - 0.009364782558895019im … 0.007943149452646923 + 0.0058778810267100445im -0.011963920209160643 - 0.008959534675119004im;;; 0.0008718402083373537 + 0.002431324664852242im 0.0003074809928018304 - 0.0016520608380871462im … 0.00078950046797984 + 0.00036548995263132884im -0.0012941118012783476 - 0.001247475749213462im; 0.01062848950865577 + 0.007717822750542854im -0.010172619729812265 - 0.00711386655990714im … 0.005629892334556453 + 0.004802984993265039im -0.009140169746814026 - 0.007372700899276449im; … ; 0.013799827868117618 + 0.011705894148800404im -0.011835562917831081 - 0.010253333254711408im … 0.00824395286815189 + 0.006686559967634851im -0.012268180848158674 - 0.010145451648737862im; -0.01169512235049847 - 0.01057055083841071im 0.009796108829266103 + 0.009245463388861813im … -0.006820618354632536 - 0.0053969627685758435im 0.010828384537045505 + 0.00898607514636411im;;; -0.0005788528608098003 - 0.0016787136838219492im 0.00011494728444778803 + 0.0012351295468039565im … -0.0004839529513788961 - 0.0005023281310356856im 0.000665974422679478 + 0.0010979334951953437im; -0.006930276764620161 - 0.004844104711101232im 0.006150849972150098 + 0.0041693698623594885im … -0.0036144076505730013 - 0.0028663722477255227im 0.005718124705094261 + 0.004277067789248072im; … ; -0.00924079532776061 - 0.00800390444675669im 0.0079506441639016 + 0.007017417578016312im … -0.005790773735318841 - 0.00482158600403396im 0.008199207577604404 + 0.006955571472090255im; 0.007728069486366248 + 0.007164507864967034im -0.006299766628111294 - 0.00603070048306332im … 0.004421146834628885 + 0.003716902370444009im -0.0067322517962997395 - 0.005922939660926429im;;; … ;;; -0.00030714594451062443 - 0.0006971697707206011im 0.0003711501328807617 + 0.0005566904023229902im … -3.060117432871485e-5 - 0.0004161924847552037im 0.00015421258846731246 + 0.0006108909193859792im; 0.0037794129823498935 + 0.0034658115551399657im -0.0034414328700834544 - 0.003030230317382612im … 0.002283436166561394 + 0.0022155704081755446im -0.0032532280960934495 - 0.0030771845594460383im; … ; 0.004785545574537361 + 0.0036895069726016003im -0.004256054348589154 - 0.0033554877940100893im … 0.0033808361811868306 + 0.0025211491234888323im -0.004388243123553794 - 0.0033224270596915162im; -0.0032559170373671403 - 0.0022978204581289067im 0.002796534711632169 + 0.0020801357770573754im … -0.0022260194500653852 - 0.0014786748214712772im 0.002984573704493748 + 0.00203309713694467im;;; 0.0006688562183673384 + 0.0016205608736002538im -0.0007433535384806217 - 0.00103056297645089im … -7.057295007926786e-5 + 0.0005868614371483642im -0.0001920733600401098 - 0.0011679259470151406im; -0.0077886394596886795 - 0.007110092188946384im 0.0067922895114219785 + 0.0058607537693083755im … -0.003960955520044678 - 0.003779052907736096im 0.006359408807594802 + 0.0059685042615931im; … ; -0.00882768122677857 - 0.006915192679537566im 0.007604342719156524 + 0.00610476166530187im … -0.005564686040970292 - 0.004225642516545091im 0.007852928665479117 + 0.006042799329624598im; 0.0068702893512509475 + 0.004898969205721561im -0.005658829593094135 - 0.004339488590120317im … 0.004074790208683333 + 0.0028039170206200677im -0.006091374971037572 - 0.0042315634504046275im;;; -0.0009549993133242625 - 0.0024057829073936407im 0.0013802536607929392 + 0.0011825283297642968im … 0.00041997303832033846 - 0.0006080178666043796im -0.00022167615148746544 + 0.0015873167795906126im; 0.011701522467102085 + 0.010552867939901319im -0.010864566979677429 - 0.00893923121844589im … 0.005920090556648319 + 0.0055670262731927975im -0.009831949393242372 - 0.009198127599172642im; … ; 0.013430664417327383 + 0.010732787291604174im -0.01154573193476197 - 0.00948924287829208im … 0.008077313679490418 + 0.00624719257096539im -0.011978403686310656 - 0.009381363928266986im; -0.010623092452177522 - 0.007736209965670639im 0.009104811055076614 + 0.0074204047174662845im … -0.006530643886359027 - 0.0046326932849921im 0.010137228862195018 + 0.007160865465834807im],)])]), [8.519873996973358e-7 4.6272652968819587e-7 … 2.4177783433167316e-8 4.627265225416913e-7; 4.6272652712674955e-7 2.501447783292669e-7 … 1.3582998778263973e-8 2.5014477556095344e-7; … ; 2.417778268199161e-8 1.3582998085680648e-8 … 1.5659226930204116e-9 1.3582997581878348e-8; 4.627265225223294e-7 2.501447764838254e-7 … 1.3582997731478526e-8 2.5014477134035725e-7;;; 4.627265255036582e-7 2.5014477757341426e-7 … 1.3582998802304672e-8 2.5014477480139777e-7; 2.5014477504068225e-7 1.389560917611049e-7 … 1.1279991859595977e-8 1.3895609204359555e-7; … ; 1.3582998147934459e-8 1.1279991372745466e-8 … 7.937756513601557e-9 1.1279990915819832e-8; 2.501447747162291e-7 1.3895609327866146e-7 … 1.127999119277084e-8 1.3895609027959533e-7;;; 2.417778409501953e-8 1.3582999107777612e-8 … 1.5659231886945167e-9 1.3582998790359744e-8; 1.3582998567909543e-8 1.1279991424723035e-8 … 7.937757519650091e-9 1.1279991775618347e-8; … ; 1.5659230393588988e-9 7.937757186869085e-9 … 3.571856666570456e-8 7.937756928818592e-9; 1.3582998727397212e-8 1.1279992017919642e-8 … 7.937757005288399e-9 1.1279991395200474e-8;;; … ;;; 5.791842649277127e-8 2.2816026136063262e-8 … 5.620377204835329e-10 2.281602706598743e-8; 2.2816027466432823e-8 3.4892779792397646e-9 … 1.3458513062657509e-8 3.4892785474103396e-9; … ; 5.620376872713997e-10 1.3458513491926517e-8 … 7.728316347008805e-8 1.3458513958565174e-8; 2.281602641830422e-8 3.4892778585296726e-9 … 1.3458514309940804e-8 3.4892778834786694e-9;;; 2.4177783891976423e-8 1.3582999486442903e-8 … 1.5659228447078364e-9 1.3582998087791565e-8; 1.358299882382483e-8 1.1279992454912754e-8 … 7.937756632303057e-9 1.1279991078543734e-8; … ; 1.5659228266984756e-9 7.937756599848043e-9 … 3.571856643250897e-8 7.937756659335338e-9; 1.3582998411603856e-8 1.1279991752165879e-8 … 7.937756614776036e-9 1.1279991075190315e-8;;; 4.627265262073977e-7 2.501447799777355e-7 … 1.3582997972751485e-8 2.5014477308022846e-7; 2.5014477821076743e-7 1.3895609618381385e-7 … 1.127999138459489e-8 1.3895609187705656e-7; … ; 1.3582997754636857e-8 1.1279990930293676e-8 … 7.937756283226866e-9 1.1279990651234983e-8; 2.501447736112236e-7 1.3895609274568668e-7 … 1.1279990541806141e-8 1.3895608892326125e-7;;;;], Matrix{ComplexF64}[[0.1651070590295301 - 0.20151086706091445im 0.743196480285787 + 0.6156611119534352im -7.354121875824737e-7 + 7.564293933030227e-7im -1.0047132585633454e-6 - 1.1319256940017412e-6im; -0.11202165648710333 + 0.13672087224223328im 0.03274817165095799 + 0.027128459370732275im 0.1606190142909395 - 0.21275638315817016im 0.0874017946866887 + 0.6145451158240974im; … ; 0.03559011727983978 - 0.04343724268499973im -0.01184458229014232 - 0.009812006476117432im -0.0015541565971410242 - 0.0014690918161941616im 0.00024135730265498602 + 0.013320013609357863im; -0.06381332408125452 + 0.07788327365289538im 0.023465034860092673 + 0.019438344734690343im 0.017257250393887008 - 0.003276322037582867im 0.003284638486574856 - 0.04117557677434142im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944094016106, [[-0.5541104206551859, -0.010678398148035262, 0.14495657618124813, 0.14495712769243724]], 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.620694686950157e-7, -6.293617692318201e-7, -2.2951214330451243e-7, 4.994027123492682e-7, 7.92331422081445e-7, 5.851074027127658e-7, 6.51697268726488e-7, 1.0189310831089882e-6, 3.6385921526907054e-7, -2.2010939094118293e-6  …  2.7810188536755433e-6, 1.2016957943921937e-6, 2.4456415384407077e-8, 1.0078550790177606e-7, 5.535309207564831e-7, 4.5341735375320544e-7, 4.63190461452002e-9, -1.0437054902818173e-7, 1.1211099788916805e-7, 1.255086920752197e-7], resid_history = [0.2493920927407076, 0.003766507751907085, 0.00028503912759627795, 4.68326204639584e-6, 1.948232402485184e-8, 6.197895500119976e-11, 3.896557905762414e-8, 1.2598706915451453e-9], converged = true, n_iter = 8, residual_norm = 1.2598706915451453e-9, maxiter = 100, tol = 1.0e-8, s = 0.9346690219387168, residual = [7.861562708345523e-7, 4.894016305009542e-7, 8.991724421022918e-8, -1.0445286425820919e-7, -8.642680485111505e-8, -2.7658781156109254e-8, -2.1852626834734856e-9, 3.0238758107315e-8, 1.9346258859372784e-8, -1.5317699246944987e-7  …  1.2420799618100502e-6, 5.399698465119955e-7, 1.1312742804108643e-8, 4.859741221633277e-8, 2.784273323955718e-7, 2.407253171335016e-7, 2.7205581794232555e-9, -7.48909801562142e-8, 1.1949542656052912e-7, 3.6836240409199306e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.6855491875924217e-7, -3.8169574311034866e-8, 1.2598706915451453e-9, -0.0002042900104891065, 4.689377840510274e-6, -1.9724445143053414e-8, 6.197895500119976e-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}}([[-3.823596419423574e-7, -1.2722813047107187e-6, -5.546203581737052e-7, 1.2948066491804728e-6, 2.1527865023142546e-6, 1.7003695664883136e-6, 2.0542467489160656e-6, 3.412246466285862e-6, 1.2484848509939678e-6, -7.633947627303487e-6  …  9.165967325121665e-6, 3.770889094594622e-6, 7.454072487231991e-8, 3.040884986407088e-7, 1.6361158874959848e-6, 1.2728860966693514e-6, 1.2074981593387672e-8, -2.5480304893968646e-7, 2.7120538173765413e-7, 3.670496240414561e-7], [1.1995058451576826e-6, 2.627884537637502e-6, 1.0769345789351302e-6, -2.4536593463678842e-6, -4.007689666582871e-6, -3.1191871335707443e-6, -3.7326303282635606e-6, -6.172199901639044e-6, -2.251565467263272e-6, 1.3707876995868905e-5  …  -1.5532449262847762e-5, -6.364832568462534e-6, -1.2507631374415942e-7, -5.079520422047099e-7, -2.723504215789746e-6, -2.103281649110121e-6, -1.9599935232966412e-8, 3.9878505934001635e-7, -3.9241980100771735e-7, -3.71994599499617e-7]]), H = [1.1806020698246427 0.004473343067801553 … 0.0 0.0; 0.10619291223087683 1.0243066460741375 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.185368373916862 0.09621932418532403 … 0.0 0.0; 0.0 1.0203206818435446 … 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.9257125358891891
Interacting polarizability:     1.7736548658359303

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 = 5.27e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 6.87e-11
[ Info: GMRES linsolve in iteration 2; step 2: normres = 9.93e-12
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 2.19e-13
└ * number of operations = 12
Non-interacting polarizability: 1.9257125358891625
Interacting polarizability:     1.773654868892374

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