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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457414491731474
Displaced dipole:  0.017612220616527646
Polarizability :   1.7746794761444962

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.0021605787704187377 + 0.0024889385571803166im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.23193363758982344 + 0.20477544128885541im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.030902132103408267 + 0.027044501567847753im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.061324357490302794 - 0.0528234321909613im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620693158885135e-7 -2.5020139590848575e-7 … -4.8806415691385935e-8 -2.502014379133178e-7; -6.293615279978113e-7 -4.793035111834006e-7 … -1.2077840260582565e-7 -4.793035439556418e-7; … ; 1.3626259559868728e-7 1.1356228723961418e-7 … 4.711987836194172e-8 1.1356227767768453e-7; 1.330547975791469e-7 1.3632644269591273e-7 … 5.200684630394566e-8 1.363264123939299e-7;;; -2.5020141711582726e-7 -1.7404899192742934e-7 … -3.5564928743162877e-8 -1.740490217244317e-7; -4.793035268666178e-7 -3.679730951796633e-7 … -1.1202890733620479e-7 -3.6797311751183e-7; … ; 1.1356227893054486e-7 1.1211096581694913e-7 … 1.098068047947702e-7 1.1211095789337196e-7; 1.363264254147432e-7 1.2550868485237667e-7 … 5.094907614341978e-8 1.2550866253686564e-7;;; -4.880640542775194e-8 -3.556491837591871e-8 … -1.1412140117573965e-8 -3.556492540793569e-8; -1.2077839265146126e-7 -1.1202889841281575e-7 … -9.786221253550475e-8 -1.1202890499851526e-7; … ; 4.7119875749360854e-8 1.0980680888496736e-7 … 2.493822681374104e-7 1.0980680460440102e-7; 5.2006850722443416e-8 5.094908426428791e-8 … 4.910626501560585e-8 5.094907853154658e-8;;; … ;;; 6.475758578004042e-8 4.11033187149123e-8 … -6.681690826686438e-9 4.1103327534106666e-8; 1.689755770108822e-7 6.65328676090206e-8 … -1.3294447446047435e-7 6.65328703684454e-8; … ; 3.267762251223193e-8 1.601736732854754e-7 … 3.8470362971914207e-7 1.6017366525310696e-7; -9.14199456226287e-8 -3.583036293921074e-8 … 7.037284199973563e-8 -3.583035983018824e-8;;; -4.880640524328108e-8 -3.5564918581804506e-8 … -1.1412142680227019e-8 -3.556492638071438e-8; -1.207783926968371e-7 -1.1202889883829284e-7 … -9.78622229673674e-8 -1.1202890921109942e-7; … ; 4.711988206901527e-8 1.0980681680304387e-7 … 2.4938227276520105e-7 1.0980680968320974e-7; 5.20068534991961e-8 5.094908756539842e-8 … 4.9106268948515716e-8 5.0949081347722256e-8;;; -2.502014163905777e-7 -1.7404899118171958e-7 … -3.55649312119915e-8 -1.74049022880833e-7; -4.79303527831496e-7 -3.6797309508247473e-7 … -1.1202891639799187e-7 -3.6797312237503535e-7; … ; 1.1356228325589951e-7 1.1211097035824805e-7 … 1.0980680740693616e-7 1.1211096152149985e-7; 1.3632642865339443e-7 1.255086877829717e-7 … 5.0949078979508225e-8 1.25508665260645e-7;;;;], δHψ = Matrix{ComplexF64}[[0.0012567815307922574 - 0.0014477843364482045im 0.06634585437767956 - 0.06840871219259585im 1.2680013592801451 - 0.6610762733828494im 0.1952478698410672 - 1.2620434398304856im; -0.14856256426486852 - 0.1313497360191304im -1.0996754343622444 - 1.198911617966199im 0.28021435098646674 - 0.14642459662136345im 0.04237204025358908 - 0.2786957258303082im; … ; -0.036800518499869234 - 0.03238325176988102im 0.014885781152599269 + 0.014364518325298609im -0.022809140473296078 + 0.018364519689247273im -0.0041345379494911105 + 0.010736036994054359im; 0.054832753571518594 + 0.04723807614120283im -0.018330490872100848 - 0.01785760629442294im -0.06824167371930065 + 0.07748120421163329im 0.014708597353760955 - 0.09575610688656244im]], δVind = [0.0019654314131352492 0.0020809998813872354 … 0.0031536644617290244 0.002081000236310543; 0.003889130369931871 0.004779268674614709 … 0.010276670912217518 0.004779269070971354; … ; -0.006150848440460172 -0.00826960436115424 … -0.017302970846311952 -0.008269603829698026; 0.00023868176378493192 -0.0003986981272782635 … -0.00366469531319161 -0.0003986977447045996;;; 0.002081000069352959 0.002223227493575364 … 0.003420119181289898 0.0022232278860037508; 0.004779268877312383 0.005734952779767446 … 0.010903597653146484 0.005734953225556966; … ; -0.008269604161436931 -0.009601009328325846 … -0.01262311245383653 -0.009601008755477387; -0.0003986979221838184 -0.0010730138518557366 … -0.004230681294283186 -0.001073013433254817;;; 0.0031536639043496066 0.003420118202568065 … 0.004832504698506603 0.0034201188568336654; 0.010276670321575 0.01090359657363244 … 0.01227214149498368 0.010903597261327636; … ; -0.01730297238254203 -0.012623113599420732 … -0.00962026573037823 -0.012623112872058164; -0.00366469597892778 -0.00423068229722132 … -0.005491109186258141 -0.004230681618397668;;; … ;;; -0.0023365593733509107 -0.0028154616232004452 … 0.0056933173601733155 -0.0028154622280827267; -0.01291309144760496 -0.017841177792764937 … 0.011498339661523818 -0.01784117871990215; … ; -0.02517446810581486 -0.012810887993823628 … -0.008395875476275095 -0.012810887420162224; 0.007597979617110109 0.010225414641305015 … -0.005449734031761603 0.010225413782475348;;; 0.0031536638435431022 0.003420118114998247 … 0.004832504258823963 0.0034201187764899914; 0.010276670320392079 0.010903596548212716 … 0.012272140999895086 0.010903597114121952; … ; -0.017302971681397765 -0.012623113422400413 … -0.009620265798153392 -0.012623112827437288; -0.003664696026500926 -0.0042306823194234915 … -0.005491109091313862 -0.004230681630765026;;; 0.0020810000533037776 0.002223227472735389 … 0.003420119093730523 0.0022232278709518097; 0.004779268873252015 0.005734952776380511 … 0.0109035973170009 0.005734953202547965; … ; -0.008269604112221766 -0.009601009256027996 … -0.012623112411038216 -0.009601008706085722; -0.00039869794851093705 -0.0010730138809334029 … -0.004230681218658778 -0.0010730134494906294;;;;], δρ0 = [-3.6021390983352183e-7 -2.488841873629119e-7 … -4.8528347267749276e-8 -2.488841884954324e-7; -5.972249040995352e-7 -4.60697725799409e-7 … -1.1905764795134771e-7 -4.606977285798378e-7; … ; 1.2735789314617446e-7 1.0836338769492678e-7 … 4.640716664797025e-8 1.0836338504686778e-7; 1.035201237644708e-7 1.1957965229448863e-7 … 5.069506553020681e-8 1.1957965173688382e-7;;; -2.4888418727867944e-7 -1.7310666086694828e-7 … -3.5356959407762985e-8 -1.7310666091837504e-7; -4.6069772582710236e-7 -3.572169940377681e-7 … -1.1091435739387893e-7 -3.5721699429484314e-7; … ; 1.0836338310003713e-7 1.0846495943001995e-7 … 1.0862235448235996e-7 1.0846495773457785e-7; 1.1957965050724578e-7 1.1611149614724825e-7 … 5.0208306927407236e-8 1.1611149583945114e-7;;; -4.852834584465577e-8 -3.5356959161729044e-8 … -1.1341305965582617e-8 -3.5356959400844254e-8; -1.1905764460195655e-7 -1.1091435763944548e-7 … -9.770665358215864e-8 -1.1091435804172606e-7; … ; 4.640716188924302e-8 1.086223509576827e-7 … 2.484524454781721e-7 1.0862235179680855e-7; 5.069506338550295e-8 5.0208305863623514e-8 … 4.926858390942979e-8 5.020830628276505e-8;;; … ;;; 6.434568447590618e-8 4.0838612017806e-8 … -6.6379254745567e-9 4.083861205698335e-8; 1.6975035310759418e-7 6.683085771769259e-8 … -1.334047910143522e-7 6.683085703709243e-8; … ; 3.270997521746791e-8 1.6023523557467465e-7 … 3.8390146148247425e-7 1.6023523421378697e-7; -9.271369276361022e-8 -3.633622236563803e-8 … 7.126309904352828e-8 -3.633622267798693e-8;;; -4.852834641040664e-8 -3.5356959956167214e-8 … -1.1341308682837168e-8 -3.5356960880734486e-8; -1.1905764525092291e-7 -1.1091435863276795e-7 … -9.77066643842384e-8 -1.1091436270767729e-7; … ; 4.6407167967586465e-8 1.0862235844275395e-7 … 2.484524492632412e-7 1.086223564789853e-7; 5.0695065638641596e-8 5.0208308666697794e-8 … 4.926858747646477e-8 5.0208308648614984e-8;;; -2.4888418816481555e-7 -1.7310666135940434e-7 … -3.535696211508648e-8 -1.7310666315309028e-7; -4.606977280470031e-7 -3.5721699494122123e-7 … -1.1091436659437385e-7 -3.572169998819147e-7; … ; 1.0836338705482785e-7 1.0846496365393966e-7 … 1.08622356864965e-7 1.0846496105889974e-7; 1.1957965266107565e-7 1.161114982384137e-7 … 5.02083094957909e-8 1.1611149773853944e-7;;;;], δeigenvalues = [[0.0004950052864777332, 0.09679859832809583, 0.029839772705108385, 0.03538878277871591]], δ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.15475965283571877 -0.15461620271096468 … -0.1541875310109779 -0.15461620271117232; -0.15461620271071153 -0.1544722989098573 … -0.154042294687657 -0.15447229891006548; … ; -0.1541875310114438 -0.15404229468838151 … -0.15360833735965623 -0.15404229468859207; -0.1546162027114128 -0.154472298910561 … -0.15404229468836855 -0.15447229891076955;;; -0.15461620271171891 -0.15447229891086856 … -0.15404229468867567 -0.15447229891107545; -0.1544722989106145 -0.15432794310177778 … -0.15389660284099177 -0.15432794310198533; … ; -0.15404229468914327 -0.15389660284171863 … -0.15346128702402026 -0.15389660284192824; -0.15447229891131686 -0.1543279431024825 … -0.1538966028417039 -0.15432794310269018;;; -0.15418753101207583 -0.15404229468901678 … -0.15360833736029178 -0.15404229468922415; -0.15404229468876 -0.15389660284133472 … -0.15346128702362946 -0.15389660284154266; … ; -0.1536083373607643 -0.15346128702436396 … -0.1530219131493626 -0.15346128702457393; -0.15404229468946817 -0.15389660284204526 … -0.15346128702434728 -0.15389660284225323;;; … ;;; -0.15347880586473917 -0.15333139197448187 … -0.15289092493780068 -0.15333139197469997; -0.15333139197421972 -0.15318350575199594 … -0.1527416287891468 -0.15318350575221487; … ; -0.15289092493828074 -0.15274162878989447 … -0.15229553960161735 -0.15274162879011569; -0.1533313919749485 -0.15318350575272732 … -0.15274162878988592 -0.15318350575294615;;; -0.15418753100951466 -0.15404229468644506 … -0.15360833735770638 -0.1540422946866583; -0.15404229468618763 -0.15389660283875162 … -0.15346128702103223 -0.15389660283896553; … ; -0.15360833735817875 -0.15346128702176748 … -0.15302191314675248 -0.1534612870219838; -0.15404229468690264 -0.15389660283946904 … -0.15346128702175751 -0.15389660283968315;;; -0.1546162027104302 -0.15447229890957465 … -0.154042294687375 -0.1544722989097843; -0.1544722989093203 -0.15432794310047834 … -0.1538966028396851 -0.15432794310068856; … ; -0.15404229468784256 -0.15389660284041248 … -0.1534612870227073 -0.15389660284062526; -0.15447229891002592 -0.15432794310118633 … -0.15389660284040108 -0.15432794310139697]), 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.015163897708148633 -0.012550129149350803 … -0.004966868132102471 -0.012550129150582432; -0.01255012911585163 -0.010362595039857966 … -0.004135158901069842 -0.010362595061527717; … ; -0.004966868225162719 -0.004135159043805244 … -0.0020697801864089996 -0.004135158974783336; -0.012550129184838131 -0.010362595132947256 … -0.004135158916957038 -0.010362595104793773;;; -0.01255012912731763 -0.010362595058060926 … -0.004135158814987726 -0.010362595046300042; -0.010362595040692962 -0.008622038956622533 … -0.0038975112102495133 -0.00862203895122518; … ; -0.004135158938305996 -0.003897511368371397 … -0.003484106911469252 -0.003897511324260846; -0.010362595082101916 -0.00862203901215538 … -0.003897511251313702 -0.008622038985188329;;; -0.004966868051222954 -0.004135158799339175 … -0.0020697799575493283 -0.004135158815366477; -0.004135158839716055 -0.003897511224128212 … -0.0034841067817425204 -0.003897511229866875; … ; -0.0020697800575198593 -0.0034841068514272224 … -0.00562057995373956 -0.003484106864670754; -0.004135158827180736 -0.003897511216996944 … -0.0034841068144317654 -0.0038975112318004297;;; … ;;; -0.006547898153054726 -0.00487634026687921 … -0.0014866440962214948 -0.0048763402602529475; -0.004876340287392089 -0.0026779889659460238 … -0.004123054264525698 -0.0026779889468549283; … ; -0.001486644115366051 -0.004123054288620851 … -0.007170895282633984 -0.004123054262792685; -0.0048763402424318765 -0.0026779889061183146 … -0.004123054270166181 -0.0026779889195917474;;; -0.004966868085600058 -0.004135158858519367 … -0.0020697802650790098 -0.004135158911818114; -0.00413515884044615 -0.0038975112345539047 … -0.003484107012949195 -0.0038975113192458964; … ; -0.0020697802226942964 -0.003484106987056663 … -0.005620579987634482 -0.0034841069444322653; -0.004135158911639791 -0.00389751132723308 … -0.00348410696683976 -0.0038975113279983327;;; -0.012550129152407469 -0.010362595080202971 … -0.004135158996243402 -0.010362595107308754; -0.01036259505268173 -0.008622038957130072 … -0.0038975114056390868 -0.008622039014924067; … ; -0.004135159019843465 -0.0038975114522313 … -0.003484106954582706 -0.003897511390042716; -0.010362595124194669 -0.00862203905736262 … -0.0038975113764370285 -0.008622039044467941])], 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.009532892475012794 -0.006919175286507595 … 0.0006640012195979587 -0.006919175287946836; -0.0069191752527551255 -0.004731678542027454 … 0.0014956970772529322 -0.004731678563905372; … ; 0.0006640011260718882 0.0014956969337929426 … 0.00356104640649859 0.0014956970026042968; -0.006919175322442843 -0.004731678635820403 … 0.001495697060654222 -0.004731678607875475;;; -0.006919175265228681 -0.004731678561241827 … 0.0014956971623162523 -0.004731678549687805; -0.004731678543619649 -0.0029911475557426216 … 0.0017333330922910586 -0.002991147550552825; … ; 0.0014956970385304668 0.0017333329334422564 … 0.002146715515036864 0.0017333329773431974; -0.00473167858573093 -0.0029911476119801553 … 0.0017333330505148285 -0.0029911475852207985;;; 0.0006640012993795203 0.0014956971776236873 … 0.003561046634722714 0.0014956971613890226; 0.0014956971375036845 0.00173333307806944 … 0.0021467156451544216 0.001733333072122832; … ; 0.0035610465342797 0.0021467155747351683 … 1.0229315853219068e-5 0.0021467155612816657; 0.0014956971493309034 0.001733333084490165 … 0.0021467156117474175 0.0017333330694787067;;; … ;;; -0.0009170700922242142 0.0007544783384472826 … 0.00414415825862401 0.0007544783448554422; 0.0007544783181966108 0.002952824485631264 … 0.00150774533546993 0.0029528245045034236; … ; 0.004144158238999421 0.0015077453106270419 … -0.0015401023617193383 0.001507745336233996; 0.0007544783624280731 0.0029528245447276416 … 0.0015077453290903995 0.002952824531035384;;; 0.0006640012675635902 0.0014956971210151881 … 0.0035610463297784368 0.0014956970675031947; 0.0014956971393459487 0.0017333330702268474 … 0.002146715416544975 0.0017333329853209435; … ; 0.003561046371690778 0.0021467154417022065 … 1.0229284568376296e-5 0.002146715484110277; 0.0014956970674373515 0.0017333329768302462 … 0.0021467154619291846 0.0017333329758508868;;; -0.0069191752890298 -0.004731678582089963 … 0.0014956969823612578 -0.004731678609405383; -0.004731678554314202 -0.0029911475549506995 … 0.0017333328982081899 -0.0029911476129549155; … ; 0.00149569695829368 0.001733332850888503 … 0.0021467154732363316 0.0017333329128643125; -0.004731678626532744 -0.0029911476558912367 … 0.0017333329266943209 -0.0029911476432071954]), 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.3471869794020585e-6 - 1.8865065307015032e-6im -0.006692607212294305 + 0.003945148294332172im … -0.004502602813795432 + 0.0026404927850408894im 0.00662980102825426 - 0.0039019423073063317im; -0.013644714376058187 - 0.01188039129848933im 0.0164729639966067 + 0.008595545845137362im … -0.004637504416483884 - 0.008513112760625465im 0.00866643875211517 + 0.013193768358705902im; … ; -0.016627545807367704 - 0.014505954851466728im 0.013269545997980345 + 0.013538866495953927im … -0.010952026688559328 - 0.007383570888397302im 0.0159491832238971 + 0.011960398140178557im; 0.013686634668000738 + 0.011844305875283049im -0.008647013428872558 - 0.013206642568995751im … 0.010814367036554142 + 0.0048971608607748865im -0.016453425861446276 - 0.008608523579487118im;;; 0.0019259145134999053 - 0.00199956264262585im 0.0027962656093892895 - 0.0011441678415030946im … 0.0034746428403351058 - 0.0022096449641122434im -0.005010534076094703 + 0.003453548671758636im; 0.01142349181415173 + 0.012081191594952505im -0.012616121787770053 - 0.008653625104430778im … 0.0039204513769290844 + 0.007323838118053837im -0.007589716544822172 - 0.011613742371238814im; … ; 0.015002090154615906 + 0.012341117059235332im -0.011864094622049636 - 0.011307445676131862im … 0.009678021263252266 + 0.006291245325341779im -0.01396760291845995 - 0.010068556102882476im; -0.013696318896289677 - 0.00972113941784105im 0.009019441541887913 + 0.01012386904350683im … -0.009130524624464006 - 0.004003004337338474im 0.014045721962613588 + 0.007163843612293067im;;; -0.0012946440936174688 + 0.0013427256455168252im -0.00046159648284827004 - 0.00012261497116774266im … -0.0019055280200709493 + 0.0013052341116975249im 0.002218253266190114 - 0.0017007459208715297im; -0.006840877884729551 - 0.007623644045371709im 0.006925571529108808 + 0.005720584051074432im … -0.002586283216208707 - 0.004718026931315701im 0.004821934205961777 + 0.006959360521686523im; … ; -0.009892701579044174 - 0.007809954131257158im 0.008008557866549524 + 0.007188961583465208im … -0.0067839024979511175 - 0.0042899953498943665im 0.009217613117858375 + 0.006476888352302467im; 0.008610873571083392 + 0.005786728764155931im -0.006125389563386364 - 0.005606282959879987im … 0.005773768586547504 + 0.002537679543630429im -0.008228983039093404 - 0.004367531200704728im;;; … ;;; -0.0006171436193243495 + 0.0006488479099820913im 0.0010514417211371942 - 0.0008620963673045305im … 0.00037434456718233265 - 5.589785036066392e-5im -1.5715956069403953e-6 - 0.00024165504696307705im; 0.0041828299174060096 + 0.002623467232205207im -0.004151401802798885 - 0.002063876190867362im … 0.0020039197383157356 + 0.002078065883956672im -0.0032375316486275836 - 0.0026022846695415254im; … ; 0.004803071619398364 + 0.004851507515306464im -0.004015981939568678 - 0.0045676399824462685im … 0.003780624426988762 + 0.0030329454177098992im -0.004659054292402919 - 0.004188708782158225im; -0.003133962430279394 - 0.0037266406631844056im 0.0023277609973702917 + 0.0035592914546256826im … -0.0026726314872550106 - 0.0019806033401399067im 0.0032415926466242712 + 0.00302081558468594im;;; 0.001365686354040948 - 0.0014349539422143906im -0.0022879756504130496 + 0.0017848759108985662im … -0.0010280975580858858 + 0.00038929369867854165im 0.0003915892767093238 + 0.00020619781820228498im; -0.00866716940028419 - 0.005716284382214478im 0.008288615950506877 + 0.004296868937199183im … -0.0033215429012122166 - 0.003950287340614103im 0.0061850713878171554 + 0.005535923192327532im; … ; -0.009015547533344423 - 0.00872603391926776im 0.0072734650951691935 + 0.00795682938218799im … -0.006304308621652584 - 0.004790870917589807im 0.00848252804776949 + 0.007244517728069626im; 0.006784521120612393 + 0.007693819028672373im -0.00476236873345387 - 0.007029875607173124im … 0.0050385352102714585 + 0.0033052515953689867im -0.006865862974719983 - 0.005790659208964542im;;; -0.0019760678374757482 + 0.0020735727668689677im 0.005082177948427774 - 0.0035312731728404743im … 0.00268998628293705 - 0.0013907051336137072im -0.002724036475772959 + 0.001067124502026375im; 0.013709154237201744 + 0.009694297950906629im -0.014087614672947248 - 0.007116686912154956im … 0.004536262681404142 + 0.006680910765039731im -0.009061405670166243 - 0.010077107690323809im; … ; 0.014217882884868409 + 0.013160153861299768im -0.011248511996183343 - 0.011950444009674537im … 0.009324314704889551 + 0.006660668301841114im -0.013352051212562938 - 0.010711325905490395im; -0.011410480588354428 - 0.01210771708183037im 0.007547949838837111 + 0.011660648830645135im … -0.008514788111560676 - 0.004645741801902477im 0.012574061919678326 + 0.00870004637145217im],)])]), [8.51987342155646e-7 4.6272648617277987e-7 … 2.417777760368781e-8 4.6272648631895314e-7; 4.627264821968486e-7 2.5014474465731265e-7 … 1.3582993517483898e-8 2.501447463333146e-7; … ; 2.4177779031590606e-8 1.3582994990025777e-8 … 1.5659213841907038e-9 1.3582994277952172e-8; 4.6272649038467164e-7 2.5014475185716274e-7 … 1.35829936813817e-8 2.501447496796628e-7;;; 4.627264835577263e-7 2.5014474606518554e-7 … 1.3582992629406904e-8 2.5014474515556575e-7; 2.5014474472189155e-7 1.389560693346136e-7 … 1.1279985656504618e-8 1.3895606905711215e-7; … ; 1.3582993901629646e-8 1.1279987092188116e-8 … 7.937752698039592e-9 1.1279986691681389e-8; 2.501447479246054e-7 1.3895607218978402e-7 … 1.1279986029351855e-8 1.3895607080329344e-7;;; 2.417777636268701e-8 1.3582992467966073e-8 … 1.5659208467743982e-9 1.3582992633316093e-8; 1.3582992884518032e-8 1.127998578252019e-8 … 7.937751772741361e-9 1.127998583462515e-8; … ; 1.5659210815292437e-9 7.937752269781491e-9 … 3.5718556695356175e-8 7.937752364244587e-9; 1.3582992755198815e-8 1.1279985717771973e-8 … 7.937752005905441e-9 1.12799858521817e-8;;; … ;;; 5.791843076738664e-8 2.281603050647156e-8 … 5.62036918049051e-10 2.2816030408784718e-8; 2.2816030808894075e-8 3.4892799900058307e-9 … 1.3458509201257682e-8 3.489279912485724e-9; … ; 5.620369404076324e-10 1.3458509448266668e-8 … 7.72831495027878e-8 1.3458509183493999e-8; 2.2816030146044193e-8 3.4892797470708307e-9 … 1.3458509259081917e-8 3.4892798017806846e-9;;; 2.417777689015884e-8 1.3582993078507807e-8 … 1.5659215689257249e-9 1.3582993628364238e-8; 1.358299289205342e-8 1.1279985877178447e-8 … 7.937753421865276e-9 1.1279986646149172e-8; … ; 1.5659214693960924e-9 7.937753237180616e-9 … 3.5718557376051267e-8 7.937752933155072e-9; 1.3582993626526426e-8 1.1279986718669475e-8 … 7.93775309298024e-9 1.1279986725619004e-8;;; 4.627264865355582e-7 2.501447477777285e-7 … 1.358299449934966e-8 2.501447498741957e-7; 2.501447456491443e-7 1.3895606936070496e-7 … 1.1279987430563616e-8 1.3895607233212963e-7; … ; 1.3582994742822475e-8 1.1279987853602295e-8 … 7.937753005553817e-9 1.1279987288954144e-8; 2.501447511802102e-7 1.3895607451407358e-7 … 1.1279987165422886e-8 1.389560738511061e-7;;;;], Matrix{ComplexF64}[[0.17077572027523508 - 0.196729823458917im 0.6718893769946265 - 0.6927801096253855im -1.7059750173342706e-7 + 1.3278243308386295e-7im -2.264125208325004e-6 + 3.618555329150347e-6im; -0.11586772385132839 + 0.1334770353495453im 0.029606109524587414 - 0.03052663365829313im -0.20905765791353864 - 0.40101196661743155im -0.3995778530344733 - 0.062283439899188414im; … ; 0.03681204165732976 - 0.04240665151112846im -0.010708134210833655 + 0.011041076726277462im 0.00016459231338910544 - 0.003369230015756768im -0.011589352735126936 - 0.001929018071398141im; -0.06600424280964806 + 0.0760354165579423im 0.02121365086904865 - 0.021873236548741345im -0.012138619088000242 - 0.0028191068188794556im 0.04292672009550956 + 0.007375447248499136im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944124718953, [[-0.5541104252599477, -0.010678399683842911, 0.14495656203278873, 0.1449567623115975]], 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.620693158885135e-7, -6.293615279978113e-7, -2.295120317017455e-7, 4.994026223045915e-7, 7.923312645962285e-7, 5.85107309147537e-7, 6.516971945112156e-7, 1.018930965954701e-6, 3.6385928534733087e-7, -2.201093250213291e-6  …  2.781018373727849e-6, 1.2016955777510656e-6, 2.4456377242459133e-8, 1.0078545026114585e-7, 5.535308004945035e-7, 4.5341727088266994e-7, 4.6318959789519464e-9, -1.043705480554327e-7, 1.1211096152149985e-7, 1.25508665260645e-7], resid_history = [0.24939209308744453, 0.0037664957815568842, 0.00028518990851962903, 4.681842689378919e-6, 1.9316490055703553e-8, 5.864063340046711e-11, 3.1037728899635125e-8, 8.714642054623401e-10], converged = true, n_iter = 8, residual_norm = 8.714642054623401e-10, maxiter = 100, tol = 1.0e-8, s = 0.935228794513672, residual = [6.829882032878631e-7, 4.489436898712773e-7, 9.017005025499718e-8, -1.2042712747103392e-7, -1.2562415337291328e-7, -6.456094804919713e-8, -5.153291086614667e-8, -5.708706942068914e-8, -1.3885230436729262e-8, 5.498797795512421e-8  …  8.47692023326427e-7, 3.744368140563278e-7, 7.962694651851577e-9, 3.458008252335335e-8, 2.0044219684512462e-7, 1.7669841328695825e-7, 2.0557798357636853e-9, -5.8575932795724846e-8, 9.680268059536153e-8, 3.0945491044991333e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [2.2753912017740848e-7, -3.066549578198418e-8, 8.714642054623401e-10, -0.0002043993745504295, 4.686806060479118e-6, -1.9519079626879922e-8, 5.864063340046711e-11, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], V = KrylovKit.OrthonormalBasis{Vector{Float64}}([[6.068306689518551e-7, 3.3856273284754585e-6, 1.1834473267887993e-6, -2.0374793608171078e-6, -2.549527005145084e-6, -1.6636506026005531e-6, -1.781666612682391e-6, -2.5875835187247004e-6, -7.876933224699688e-7, 4.075754730546686e-6  …  -6.070189201635461e-6, -2.822541761883617e-6, -6.596645091163732e-8, -2.9722935673746936e-7, -1.6431839750657744e-6, -1.334190679946154e-6, -1.4676841514986263e-8, 3.9755496707662554e-7, -5.546195596005674e-7, -1.016764804392137e-6], [7.52985765701341e-7, -1.2273862717107413e-6, -5.227547160682908e-7, 9.277736519011093e-7, 1.1661268317114373e-6, 7.700136109300128e-7, 8.339822862179181e-7, 1.204360743012807e-6, 3.5852817949781675e-7, -1.82356936560606e-6  …  4.504440426459086e-6, 2.108010051093568e-6, 4.905244332279415e-8, 2.2125144017721148e-7, 1.2452010975960882e-6, 1.0444916086695956e-6, 1.186745205109921e-8, -3.308823742903082e-7, 4.899497335497861e-7, 1.1240406714023116e-6]]), H = [1.0884538771478105 0.011856828191237321 … 0.0 0.0; 0.13769369908868165 1.0209007620110242 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.0971287059624433 0.13988988903953928 … 0.0 0.0; 0.0 1.0117394614002726 … 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.9257125142879803
Interacting polarizability:     1.7736548458357393

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

Manual solution of the Dyson equations

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

using KrylovKit

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

# Apply ``χ_0`` once to get non-interacting dipole
δρ_nointeract = apply_χ0(scfres, δVext).δρ

# Solve Dyson equation to get interacting dipole
δρ = linsolve(dielectric_operator, δρ_nointeract, verbosity=3)[1]

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

WARNING: using KrylovKit.basis in module Main conflicts with an existing identifier.
[ Info: GMRES linsolve starts with norm of residual = 4.19e+00
[ Info: GMRES linsolve in iteration 1; step 1: normres = 2.49e-01
[ Info: GMRES linsolve in iteration 1; step 2: normres = 3.77e-03
[ Info: GMRES linsolve in iteration 1; step 3: normres = 2.85e-04
[ Info: GMRES linsolve in iteration 1; step 4: normres = 4.69e-06
[ Info: GMRES linsolve in iteration 1; step 5: normres = 1.09e-08
[ Info: GMRES linsolve in iteration 1; step 6: normres = 6.27e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 5.00e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 6.95e-11
[ Info: GMRES linsolve in iteration 2; step 2: normres = 4.32e-12
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 2.73e-13
└ * number of operations = 12
Non-interacting polarizability: 1.9257125142878941
Interacting polarizability:     1.7736548428815262

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