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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457354746192448
Displaced dipole:  0.017612221178745682
Polarizability :   1.7746794726207606

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.003294819394920525 + 8.417581150154338e-5im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.010498668095058412 - 0.3092183246796922im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.0012142643496665633 - 0.041047197241035276im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.001757177253114789 + 0.08091912124575983im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206943699603454e-7 -2.502014762372775e-7 … -4.8806445887262506e-8 -2.502015397818735e-7; -6.293616273253034e-7 -4.793035829244113e-7 … -1.2077844127508434e-7 -4.793036312762556e-7; … ; 1.3626258578465176e-7 1.1356228412368777e-7 … 4.711989220624146e-8 1.135622750394276e-7; 1.3305471313376658e-7 1.3632639047960015e-7 … 5.200683657813012e-8 1.3632634557029233e-7;;; -2.5020152264831935e-7 -1.7404906325454686e-7 … -3.556495636990095e-8 -1.7404910956871557e-7; -4.793036170616378e-7 -3.6797316224554894e-7 … -1.1202895076847021e-7 -3.679731981987992e-7; … ; 1.1356227686970853e-7 1.1211096639188945e-7 … 1.098068109340484e-7 1.1211095865895555e-7; 1.3632635682637976e-7 1.255086418455815e-7 … 5.0949067849363836e-8 1.2550860890531146e-7;;; -4.880643766285464e-8 -3.556494233698811e-8 … -1.1412154694020107e-8 -3.5564953160715895e-8; -1.2077843229680405e-7 -1.1202893532224215e-7 … -9.786225542187304e-8 -1.1202894632702673e-7; … ; 4.7119892041373026e-8 1.0980681771695376e-7 … 2.4938225476729787e-7 1.0980681189660575e-7; 5.200684100511533e-8 5.094907846551649e-8 … 4.910625988223814e-8 5.094906978072527e-8;;; … ;;; 6.475759436827287e-8 4.1103320078452804e-8 … -6.681698369603759e-9 4.110333298565856e-8; 1.6897556142684223e-7 6.653284733317013e-8 … -1.329445062493903e-7 6.653285086667975e-8; … ; 3.267764330265347e-8 1.6017368861042858e-7 … 3.8470362437271135e-7 1.6017367920855025e-7; -9.141992451007457e-8 -3.58303490371646e-8 … 7.03728425206195e-8 -3.58303430908439e-8;;; -4.880642452229104e-8 -3.556493254488209e-8 … -1.1412151855861249e-8 -3.5564943632215995e-8; -1.2077842430308414e-7 -1.12028927899195e-7 … -9.786225022689861e-8 -1.120289396684717e-7; … ; 4.71198965200909e-8 1.0980682767284506e-7 … 2.493822724106979e-7 1.0980682065004667e-7; 5.200685195152133e-8 5.094908863670205e-8 … 4.910626837300205e-8 5.0949079813586975e-8;;; -2.5020149307092415e-7 -1.740490415679018e-7 … -3.5564951707328285e-8 -1.7404908819080148e-7; -4.793035968114709e-7 -3.6797314716376674e-7 … -1.1202894754912071e-7 -3.679731837256571e-7; … ; 1.1356228241859947e-7 1.1211097204853449e-7 … 1.0980681528710844e-7 1.1211096362917724e-7; 1.3632637949554566e-7 1.255086588708608e-7 … 5.094907280249868e-8 1.2550862569726943e-7;;;;], δHψ = Matrix{ComplexF64}[[-0.0019165546320700246 - 4.896395673447647e-5im 0.09229405768099608 + 0.023734636317210966im -0.11722713994804024 + 1.0694284998319465im -1.2191646133450282 - 0.753603661733303im; -0.006866066162432228 + 0.198182859603568im 0.49653891300763187 - 1.5492334277005184im -0.02605559460452746 + 0.23628133219602826im -0.2693217195508793 - 0.1651929461837481im; … ; -0.001582372892788701 + 0.048994379780028435im -0.005101792941008429 + 0.020047348523826402im -0.008390309135806935 - 0.02159064635193532im 0.013325989757253194 + 0.013310885936930418im; 0.0015760808441059749 - 0.0723573263890528im 0.006429207456644777 - 0.02477030171825194im -0.09160132782825747 - 0.09290576144697553im -0.06318722454621571 + 0.0037746504562639997im]], δVind = [0.00196543202264238 0.0020810004819443394 … 0.003153665883218339 0.0020810010105154746; 0.0038891309812523816 0.004779269255638565 … 0.010276671911465833 0.004779269847617985; … ; -0.006150846816910301 -0.008269602678829839 … -0.01730296506088918 -0.008269601636859449; 0.00023868250375307164 -0.0003986973936925412 … -0.0036646932768469934 -0.000398696804939906;;; 0.002081000871091273 0.0022232282960440593 … 0.0034201209184413434 0.0022232288860100724; 0.004779269689572324 0.005734953557773422 … 0.010903598718995489 0.005734954214332322; … ; -0.008269601904869937 -0.009601007172781128 … -0.012623109123948314 -0.009601006085501207; -0.00039869695503111954 -0.001073012874610854 … -0.0042306788522621305 -0.0010730122198166477;;; 0.003153665417238921 0.003420119667240372 … 0.004832506960805647 0.003420120663181697; 0.010276671466721902 0.010903597505159157 … 0.012272142684417179 0.010903598503087116; … ; -0.017302966129926194 -0.012623110627798912 … -0.0096202635045942 -0.012623109450568413; -0.0036646938198071025 -0.004230680194096114 … -0.00549110619383528 -0.004230679138156599;;; … ;;; -0.0023365599961535717 -0.0028154622640882533 … 0.005693317662151627 -0.002815463179384841; -0.012913092773704365 -0.017841180671824518 … 0.011498339959853865 -0.017841182110581527; … ; -0.025174460038961413 -0.012810886725142788 … -0.008395874486709066 -0.012810885703140734; 0.007597979347101109 0.010225415074486014 … -0.0054497326446017146 0.010225413865274844;;; 0.0031536645191404594 0.003420118665327365 … 0.004832505448534677 0.0034201196643173777; 0.01027667044536604 0.010903596456400141 … 0.01227214152904857 0.010903597426032988; … ; -0.017302967470730814 -0.012623111670225605 … -0.009620264319587447 -0.012623110531304474; -0.003664694817862242 -0.004230681221310562 … -0.005491107241525485 -0.004230680151988445;;; 0.002081000618550169 0.0022232280117758384 … 0.0034201204284942006 0.00222322860333356; 0.004779269411142311 0.005734953246509008 … 0.010903598187153554 0.005734953900096239; … ; -0.008269602415349011 -0.009601007682850032 … -0.012623109661572077 -0.009601006608980243; -0.00039869724648649066 -0.0010730132008256986 … -0.004230679350822466 -0.0010730125417667147;;;;], δρ0 = [-3.6021405504704793e-7 -2.4888429543471653e-7 … -4.85283762508075e-8 -2.4888429668434837e-7; -5.972250446376914e-7 -4.606978341420489e-7 … -1.1905768893471389e-7 -4.606978368926576e-7; … ; 1.27357889374516e-7 1.0836338707860495e-7 … 4.640718283941661e-8 1.0836338838441204e-7; 1.0352004621231411e-7 1.1957959671327364e-7 … 5.069506002286607e-8 1.195795972591304e-7;;; -2.488842953869884e-7 -1.7310674159142742e-7 … -3.535698359396123e-8 -1.7310674231984595e-7; -4.606978343206985e-7 -3.572170796588769e-7 … -1.109144005405588e-7 -3.5721708175289404e-7; … ; 1.0836338761882973e-7 1.084649642752216e-7 … 1.0862236809028593e-7 1.0846496599858558e-7; 1.195795970419076e-7 1.1611145810258343e-7 … 5.0208304669979675e-8 1.161114589182698e-7;;; -4.85283751054717e-8 -3.5356982648778895e-8 … -1.1341318490229018e-8 -3.5356983031555155e-8; -1.1905768529989093e-7 -1.1091439703620315e-7 … -9.770669425571401e-8 -1.1091439853264973e-7; … ; 4.640718093702079e-8 1.086223650092901e-7 … 2.484524516733218e-7 1.0862236707180705e-7; 5.069505921166827e-8 5.020830342485114e-8 … 4.926858567587365e-8 5.020830421648152e-8;;; … ;;; 6.43457139978957e-8 4.083862858340948e-8 … -6.637934280466052e-9 4.083862834495895e-8; 1.6975035434813635e-7 6.683084473571667e-8 … -1.334048320436968e-7 6.683084312915395e-8; … ; 3.2709995265060506e-8 1.6023524540376228e-7 … 3.8390145490419827e-7 1.6023524608762418e-7; -9.271366227544136e-8 -3.633620384819202e-8 … 7.12630961010578e-8 -3.633620303723262e-8;;; -4.8528376017003654e-8 -3.535698327378871e-8 … -1.1341319165765158e-8 -3.535698400793438e-8; -1.1905768774146744e-7 -1.1091439904120376e-7 … -9.770669697762254e-8 -1.1091440138307816e-7; … ; 4.6407181869498506e-8 1.086223670594617e-7 … 2.484524525017619e-7 1.0862236784768614e-7; 5.0695059634034794e-8 5.020830410881509e-8 … 4.926858624350397e-8 5.020830467507977e-8;;; -2.488842966941918e-7 -1.7310674230386177e-7 … -3.535698412268757e-8 -1.7310674359986745e-7; -4.6069783660410324e-7 -3.5721708118377356e-7 … -1.1091440204815149e-7 -3.5721708406484003e-7; … ; 1.0836338793242769e-7 1.0846496519726437e-7 … 1.086223684629062e-7 1.0846496618881955e-7; 1.1957959703893475e-7 1.1611145845206005e-7 … 5.020830488055792e-8 1.161114588234175e-7;;;;], δeigenvalues = [[0.0004950053047365892, 0.09679861763749684, 0.039294942276273066, 0.03128372403892443]], δ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.15475965283263715 -0.15461620270795334 … -0.15418753100775606 -0.15461620270802062; -0.15461620270806575 -0.15447229890728215 … -0.1540422946848768 -0.15447229890735137; … ; -0.15418753100755428 -0.15404229468456104 … -0.1536083373556055 -0.1540422946846254; -0.1546162027079156 -0.15447229890713338 … -0.15404229468472158 -0.1544722989071991;;; -0.15461620270805107 -0.15447229890727035 … -0.1540422946848562 -0.15447229890733377; -0.15447229890738076 -0.1543279430986141 … -0.15389660283761208 -0.15432794309867942; … ; -0.15404229468465877 -0.15389660283730294 … -0.15346128701936354 -0.15389660283736298; -0.1544722989072308 -0.15432794309846565 … -0.153896602837457 -0.15432794309852735;;; -0.15418753100782018 -0.1540422946848298 … -0.15360833735586996 -0.1540422946848899; -0.15404229468493905 -0.15389660283758297 … -0.15346128701964862 -0.153896602837645; … ; -0.15360833735567542 -0.15346128701934367 … -0.15302191314408856 -0.15346128701939993; -0.15404229468478825 -0.15389660283743364 … -0.1534612870194924 -0.15389660283749187;;; … ;;; -0.15347880586333515 -0.1533313919731483 … -0.15289092493627185 -0.15333139197323029; -0.15333139197327175 -0.15318350575111953 … -0.15274162878808178 -0.15318350575120343; … ; -0.15289092493604825 -0.1527416287877328 … -0.1522955395992443 -0.15274162878781208; -0.1533313919731141 -0.15318350575096318 … -0.15274162878791894 -0.15318350575104364;;; -0.15418753100757376 -0.15404229468457442 … -0.15360833735563736 -0.15404229468465108; -0.15404229468469338 -0.15389660283732837 … -0.1534612870194167 -0.1538966028374069; … ; -0.15360833735542243 -0.15346128701908113 … -0.1530219131438502 -0.15346128701915507; -0.15404229468453942 -0.15389660283717563 … -0.15346128701925765 -0.15389660283725082;;; -0.154616202707927 -0.1544722989071418 … -0.1540422946847391 -0.1544722989072135; -0.15447229890725706 -0.1543279430984859 … -0.15389660283749532 -0.15432794309855952; … ; -0.1540422946845315 -0.15389660283717077 … -0.1534612870192436 -0.15389660283723974; -0.1544722989071055 -0.15432794309833572 … -0.1538966028373389 -0.15432794309840595]), 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.015163897874290287 -0.012550129332159096 … -0.00496686851914464 -0.012550129352886033; -0.012550129340839909 -0.010362595290917842 … -0.004135159348546722 -0.010362595316586105; … ; -0.004966868478150281 -0.004135159283590737 … -0.002069780754996231 -0.004135159311499834; -0.012550129342232891 -0.010362595294566786 … -0.004135159333587489 -0.010362595313407;;; -0.012550129339056748 -0.010362595290337805 … -0.004135159332989152 -0.010362595310679078; -0.01036259529704125 -0.008622039236564665 … -0.0038975117530562677 -0.008622039260382327; … ; -0.004135159295783551 -0.0038975116929669874 … -0.0034841073526990307 -0.003897511727318003; -0.010362595302055671 -0.008622039241128244 … -0.0038975117448183452 -0.008622039262085372;;; -0.0049668684684335135 -0.004135159271741604 … -0.002069780748910149 -0.00413515930186208; -0.004135159278864836 -0.0038975116814590526 … -0.0034841073502076994 -0.0038975117129821666; … ; -0.002069780699558538 -0.0034841072945943044 … -0.00562058027718331 -0.0034841073325335777; -0.004135159290531307 -0.00389751168960508 … -0.003484107350415818 -0.003897511721312917;;; … ;;; -0.0065478977528394995 -0.004876339830619176 … -0.0014866448561614127 -0.004876339797701459; -0.004876339822684903 -0.002677988363620894 … -0.004123054675292554 -0.002677988319972083; … ; -0.001486644769031863 -0.004123054626038313 … -0.007170895524077771 -0.00412305463873191; -0.004876339820581906 -0.0026779883570348523 … -0.004123054670898395 -0.0026779883208717034;;; -0.004966868497829093 -0.004135159300594682 … -0.002069780809774848 -0.00413515934320853; -0.004135159319053273 -0.0038975117152206355 … -0.0034841074005796465 -0.0038975117632256005; … ; -0.0020697807271775463 -0.003484107331906102 … -0.005620580290565615 -0.0034841073480577837; -0.0041351593147739545 -0.003897511719838474 … -0.003484107384231166 -0.0038975117509540257;;; -0.012550129346849122 -0.010362595297515844 … -0.004135159353861018 -0.010362595322810196; -0.010362595309217418 -0.008622039246109851 … -0.0038975117789705666 -0.008622039277987886; … ; -0.004135159300624447 -0.003897511709567313 … -0.0034841073604029976 -0.0038975117308883266; -0.010362595307501232 -0.008622039249660915 … -0.003897511760089807 -0.00862203927019286])], 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.00953289263807283 -0.006919175466304547 … 0.0006640008357776291 -0.006919175487098736; -0.006919175475097631 -0.0047316787905121675 … 0.0014956966325562725 -0.0047316788162496525; … ; 0.0006640008769738543 0.0014956966978279217 … 0.0035610458419620904 0.0014956966698544596; -0.006919175476340401 -0.004731678794012314 … 0.001495696647670743 -0.004731678812918253;;; -0.006919175473299954 -0.0047316787899205 … 0.0014956966481342983 -0.004731678810325168; -0.00473167879673419 -0.00299114783252109 … 0.0017333325528639896 -0.00299114785640406; … ; 0.0014956966855374084 0.001733332613262356 … 0.0021467150784638045 0.0017333325788513052; -0.00473167880159862 -0.002991147836936176 … 0.0017333325612570658 -0.002991147857955005;;; 0.0006640008864246122 0.0014956967094082417 … 0.0035610458477837175 0.0014956966792276753; 0.0014956967021758473 0.001733332624490348 … 0.002146715080670086 0.0017333325929052004; … ; 0.003561045897329895 0.002146715136588376 … 1.0228997683500198e-5 0.0021467150985928422; 0.001495696690660256 0.0017333326164936458 … 0.0021467150806182313 0.0017333325847275773;;; … ;;; -0.000917069690604972 0.0007544787760408889 … 0.004144157500212924 0.0007544788088766161; 0.0007544787838517607 0.002952825088832804 … 0.0015077449257680838 0.0029528251323977096; … ; 0.004144157587566101 0.0015077449753712396 … -0.001540102600790079 0.0015077449625983726; 0.0007544787861124376 0.0029528250955752483 … 0.0015077449303251677 0.0029528251316579338;;; 0.0006640008572754465 0.001495696680810515 … 0.00356104578715161 0.001495696638120006; 0.0014956966622330752 0.001733332590983367 … 0.0021467150305300646 0.0017333325428998815; … ; 0.003561045869963851 0.0021467150995391183 … 1.0228984539532512e-5 0.002146715083313496; 0.001495696666666409 0.0017333325865182673 … 0.00214671504703764 0.001733332555327526;;; -0.006919175480968261 -0.004731678796969976 … 0.0014956966273795332 -0.00473167882233602; -0.004731678808786652 -0.002991147841938045 … 0.0017333325270664861 -0.0029911478738896876; … ; 0.001495696680823744 0.0017333325967942024 … 0.0021467150708797417 0.0017333325754042162; -0.004731678806918892 -0.0029911478453389234 … 0.0017333325461037319 -0.0029911478659410892]), 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 = [-2.1954580547599716e-6 - 7.942571276270282e-7im -0.002714339752766786 + 0.004457574393585723im … -0.001914701286957682 + 0.0030189793493351274im 0.0027721718843983292 - 0.004426454573666907im; -0.0005423441505205929 + 0.020101037187925948im 0.00205462121984363 - 0.020414886586626664im … 0.001048505691283423 + 0.008555273072738864im -0.0011572928745988659 - 0.015205574334170332im; … ; -0.0006493946389686831 + 0.023443585732506773im 1.5293604226280864e-5 - 0.019437421645252168im … -0.0009900465011914291 + 0.014008596398464024im 0.0011174102867018603 - 0.02122826318255849im; 0.00048447196618244697 - 0.020128781250417567im 0.0011468936396956402 + 0.01520416505359299im … 0.00159058089311804 - 0.012708987636832573im -0.0020664171467647545 + 0.02041497510540219im;;; 0.00755420486603717 - 0.002074902708672208im -0.0028402675838087605 - 0.0013977410537169082im … 0.00284057296752711 - 0.002493108147174081im -0.00605603165812853 + 0.0038142566186923693im; -0.003960664175937427 - 0.01658960039837047im 0.0014417987258587193 + 0.015978108298652705im … -0.001954956988156236 - 0.007077462145515598im 0.003510635495805317 + 0.012621249785443746im; … ; 0.002076914929889475 - 0.02075116772306814im -0.0012426917468235108 + 0.01706416056561656im … 0.0015276339610996485 - 0.012360397123141622im -0.0021081934817497204 + 0.018470054675513173im; -0.004865527935479189 + 0.019032818117800228im 0.002208608559044376 - 0.014184454193193564im … -0.0024111715977362547 + 0.010834715608772048im 0.0042781014651749796 - 0.017542175688016467im;;; -0.0051771171637634665 + 0.001414076639374732im 0.003015465476669737 - 7.360952802175724e-5im … -0.0023497405283509345 + 0.0014653505144955636im 0.0041191249138719885 - 0.0018650890071885702im; 0.003288678240472961 + 0.009661361509709566im -0.002006258752933348 - 0.008935364564042102im … 0.0018281553788691888 + 0.00448779000102459im -0.00287226318691519 - 0.007529510441227027im; … ; -0.0020638239881847396 + 0.01347583846586752im 0.0014971786203355644 - 0.011321022774438916im … -0.001499068048680865 + 0.008663602899885223im 0.0019947056654460588 - 0.012129101191651032im; 0.00382817900189556 - 0.011604267379957948im -0.0024612467620160113 + 0.008977122800932331im … 0.002119020545993161 - 0.006944379393675547im -0.003327369894493319 + 0.010383272287205278im;;; … ;;; -0.0022152083247688643 + 0.0006251504280534762im 0.002081676778847785 - 0.0008795570771338517im … -0.0009283150089663 - 0.00013640162635466147im 0.0016487360655069493 - 0.00017661672984092336im; 0.001992924648124963 - 0.005529598747734448im -0.0019019660079658957 + 0.005186042074059319im … 0.0009049061341025804 - 0.0030687943703066956im -0.0015262156652095624 + 0.004575861156373753im; … ; -0.0010492404760928187 - 0.0066992250855681355im 0.001063584556322118 + 0.005834917579416472im … -0.00048155072892188723 - 0.004919143316362276im 0.0007991357391427978 + 0.0062644070521350095im; 0.0017384187799845568 + 0.004474405263405716im -0.0016767151432278476 - 0.0036669951644348227im … 0.0007425127947723732 + 0.003316089094302748im -0.0013010000371339195 - 0.004277117411471311im;;; 0.005057541688128858 - 0.0014145638855940135im -0.004011424846969682 + 0.0018665601051631967im … 0.00102808722085422 + 0.0005316732090948853im -0.0029102517951089084 + 7.810696483493997e-5im; -0.003740209109705174 + 0.01160256419594452im 0.0032393181373232835 - 0.010383241782739153im … -0.0010026580593725637 + 0.005270078546071712im 0.002374788464695778 - 0.008979292980907386im; … ; 0.0013130648311725788 + 0.012543477863698586im -0.0013327394692336401 - 0.010539939014433243im … 0.0003479811214488289 + 0.008153341282684134im -0.0008356123551648976 - 0.01134730021157538im; -0.003199934011370252 - 0.009662882611914461im 0.0027838456389769756 + 0.007529267091659346im … -0.0007117158898287382 - 0.006162091552992575im 0.0019192191164123064 + 0.008933452380618904im;;; -0.0074537236315736005 + 0.0020743856163057245im 0.005952413691417558 - 0.0038253957483151612im … -0.00017988930420501168 - 0.0016579781059296735im 0.0027410180310193807 + 0.0013813259417352407im; 0.004835737364829652 - 0.019018962402628684im -0.004220398314899631 + 0.01754046400105115im … 0.00041572131151014323 - 0.007732738125118096im -0.002153890075058569 + 0.014186553752155034im; … ; -0.0009424497847904972 - 0.01991792776259014im 0.001126988981563872 + 0.01641037217670818im … 0.00016532175213139998 - 0.01198398803017534im 0.0002619838578396212 + 0.017815438867029198im; 0.003929485521252549 + 0.01660321745357892im -0.0034529031834982737 - 0.0126220330777489im … -4.057089864665938e-5 + 0.010179443003113988im -0.001385690319506886 - 0.01597673275300286im],)])]), [8.519873723466656e-7 4.6272650786980516e-7 … 2.4177783542400463e-8 4.6272651032983056e-7; 4.6272650890010877e-7 2.5014476407516196e-7 … 1.3582998133917778e-8 2.5014476606043587e-7; … ; 2.4177782913388704e-8 1.3582997463797483e-8 … 1.565922719365119e-9 1.3582997751724493e-8; 4.627265090654222e-7 2.5014476435738825e-7 … 1.3582997979589149e-8 2.5014476581456423e-7;;; 4.627265086884672e-7 2.5014476403030643e-7 … 1.358299797341969e-8 2.5014476560357325e-7; 2.501447645487714e-7 1.3895608372758553e-7 … 1.1279990584971698e-8 1.3895608495214704e-7; … ; 1.3582997589585805e-8 1.1279990039389483e-8 … 7.9377558451931e-9 1.1279990351282167e-8; 2.50144764936611e-7 1.3895608396222636e-7 … 1.1279990510175345e-8 1.3895608503970904e-7;;; 2.417778276429958e-8 1.3582997341552057e-8 … 1.5659227050753852e-9 1.3582997652290524e-8; 1.3582997415040742e-8 1.1279989934899009e-8 … 7.937755827422976e-9 1.1279990221118723e-8; … ; 1.5659225891852023e-9 7.937755430750918e-9 … 3.571856319092454e-8 7.937755701356566e-9; 1.358299753539733e-8 1.1279990008861797e-8 … 7.93775582890834e-9 1.1279990296754722e-8;;; … ;;; 5.79184195437979e-8 2.2816024074704526e-8 … 5.620378056081732e-10 2.2816023589399305e-8; 2.2816023957727267e-8 3.4892775441968547e-9 … 1.3458513412228814e-8 3.489277366954755e-9; … ; 5.620377038451093e-10 1.3458512907301173e-8 … 7.728315776927803e-8 1.345851303742527e-8; 2.2816023926721497e-8 3.4892775174535082e-9 … 1.3458513367181787e-8 3.489277370609534e-9;;; 2.41777832153391e-8 1.3582997639217982e-8 … 1.5659228479985999e-9 1.3582998078850323e-8; 1.3582997829645727e-8 1.1279990241443017e-8 … 7.937756186707611e-9 1.127999067730784e-8; … ; 1.5659226540421936e-9 7.937755696884506e-9 … 3.5718563459680175e-8 7.937755812087909e-9; 1.3582997785497824e-8 1.1279990283369045e-8 … 7.937756070100044e-9 1.127999056588502e-8;;; 4.6272650961331994e-7 2.501447645854788e-7 … 1.3582998188747671e-8 2.5014476654184017e-7; 2.5014476549051203e-7 1.3895608421834363e-7 … 1.1279990820265589e-8 1.3895608585732483e-7; … ; 1.358299763952509e-8 1.1279990190110248e-8 … 7.93775590014109e-9 1.1279990383696324e-8; 2.50144765357785e-7 1.3895608440091904e-7 … 1.1279990648836807e-8 1.3895608545655648e-7;;;;], Matrix{ComplexF64}[[-0.2604279224516406 - 0.006653369978697369im 0.9346686695287589 + 0.24036241407309183im 3.2787696157889994e-7 - 2.2103757918591557e-6im 3.76254870550395e-6 - 4.525002574719907e-6im; 0.17669485168403712 + 0.004514171145108177im 0.04118517793140929 + 0.010591311637893739im 0.3383603536061764 + 0.036644696437715306im -0.2396182100318619 + 0.3861335781193139im; … ; -0.056137274749016475 - 0.001434185927088187im -0.01489614472782681 - 0.003830740646223038im 0.0010508024486663623 + 0.0050539228030417604im -0.003989986387603345 + 0.009755602243098752im; 0.10065451806841835 + 0.0025715051808347135im 0.029510404783953195 + 0.007588991418135962im 0.01230816740412005 - 0.026077472090355394im 0.009370310899597294 - 0.03351820824149862im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944080518443, [[-0.5541104231603817, -0.010678392943306858, 0.14495664680089407, 0.1449571327954985]], 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.6206943699603454e-7, -6.293616273253034e-7, -2.2951208303254293e-7, 4.994025945138364e-7, 7.923312220149517e-7, 5.85107244774948e-7, 6.516971234871358e-7, 1.0189308872592872e-6, 3.6385918909347967e-7, -2.2010933530305837e-6  …  2.7810181295515867e-6, 1.2016954763006736e-6, 2.4456372249359826e-8, 1.0078544626664469e-7, 5.535307741858382e-7, 4.5341726226036947e-7, 4.631921742511642e-9, -1.0437051647012074e-7, 1.1211096362917724e-7, 1.2550862569726943e-7], resid_history = [0.24939209463765716, 0.0037665480943702154, 0.00028524525683628975, 4.681486810589749e-6, 3.6891566933631397e-8, 9.195082612607024e-11, 4.220213170991931e-8, 1.164941594961988e-9], converged = true, n_iter = 8, residual_norm = 1.164941594961988e-9, maxiter = 100, tol = 1.0e-8, s = 0.9355297314016175, residual = [9.153265668247358e-7, 6.109618226781955e-7, 1.280977938052816e-7, -1.8650335328035676e-7, -2.1806083563650274e-7, -1.2606059527711034e-7, -1.1443888386542069e-7, -1.5278564255321924e-7, -4.910885432848885e-8, 2.7412384947536453e-7  …  8.779486581073341e-7, 3.8959748293109823e-7, 8.473227110409207e-9, 3.78214123781313e-8, 2.2326400622662523e-7, 1.9921965259291744e-7, 2.366465750467434e-9, -6.991865510438209e-8, 1.2065387830539037e-7, 4.0059527218153036e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.3342546149299873e-7, -4.115407468240017e-8, 1.164941594961988e-9, -0.00020443974185030956, 4.6859087492185165e-6, -3.733640199981251e-8, 9.195082612607024e-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}}([[-7.977255513861701e-7, -2.3176910903751055e-6, -7.785194317328451e-7, 1.3025514603855441e-6, 1.6421041641137817e-6, 1.1679755728406838e-6, 1.3796793898680597e-6, 2.060334970240509e-6, 5.92954353295809e-7, -2.8612254840395958e-6  …  3.2904065973139405e-6, 1.2962130059743854e-6, 3.0104928514822444e-8, 1.445393905427006e-7, 7.74965844583113e-7, 5.223812115317826e-7, 4.563001697503056e-9, -1.1785984761788033e-7, 1.7571846852665585e-7, 2.452763907869858e-7], [1.3332626299669617e-6, 2.8168628371455756e-6, 9.21029690439344e-7, -1.5527325898741525e-6, -1.9651801735880046e-6, -1.375386702129762e-6, -1.5978358551305993e-6, -2.3857742709127357e-6, -6.965863951223261e-7, 3.392705160758014e-6  …  -3.05191033637864e-6, -1.1673858108116994e-6, -2.7022536733669473e-8, -1.309746465377561e-7, -6.937639244122917e-7, -4.397104261462351e-7, -3.4042502153946528e-9, 8.0247826298681e-8, -1.0923434008855273e-7, -1.492144247958016e-8]]), H = [1.1320961450853155 0.013608632227362078 … 0.0 0.0; 0.12756039479403575 1.0327028163609602 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.139259994925233 0.1291524846828782 … 0.0 0.0; 0.0 1.025075901501192 … 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.925712521526251
Interacting polarizability:     1.7736548523318685

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

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