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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457361827735504
Displaced dipole:  0.01761222066323871
Polarizability :   1.7746794281516063

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.0008344683711474262 - 0.0031885078848987725im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.29996298304891134 + 0.07581822951616257im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.03976876263892862 + 0.010236814094783908im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.07822200763002254 - 0.020792051207303226im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620694128621567e-7 -2.5020144239902345e-7 … -4.880644308519737e-8 -2.502015347237721e-7; -6.293616532988726e-7 -4.793035850396903e-7 … -1.2077843309461443e-7 -4.793036542780553e-7; … ; 1.3626259380538864e-7 1.1356228551129811e-7 … 4.711987195096428e-8 1.1356227228959363e-7; 1.3305477156798873e-7 1.3632644031996244e-7 … 5.2006830903406433e-8 1.3632637335053588e-7;;; -2.5020149109891366e-7 -1.7404902821137423e-7 … -3.556495091602286e-8 -1.7404909538960845e-7; -4.79303620735684e-7 -3.679731520577084e-7 … -1.1202894035322196e-7 -3.6797320308830753e-7; … ; 1.1356227775392454e-7 1.1211096573928462e-7 … 1.0980679794742975e-7 1.1211095423762528e-7; 1.3632640337349177e-7 1.2550868116048037e-7 … 5.094906291954732e-8 1.255086315242019e-7;;; -4.880642217383082e-8 -3.556492738136056e-8 … -1.1412148467020925e-8 -3.5564942677539784e-8; -1.2077841363489512e-7 -1.1202891579909329e-7 … -9.786224170934236e-8 -1.120289301609373e-7; … ; 4.711987718718228e-8 1.0980681424876284e-7 … 2.4938225554814176e-7 1.0980680288895281e-7; 5.2006844027838294e-8 5.0949082891133593e-8 … 4.9106254091677815e-8 5.094906845900774e-8;;; … ;;; 6.475760883663698e-8 4.110332783813597e-8 … -6.681694663721811e-9 4.110334746719124e-8; 1.6897560005852163e-7 6.6532876051512e-8 … -1.329445046617303e-7 6.65328791791559e-8; … ; 3.2677617129039286e-8 1.6017367475421804e-7 … 3.847036072194342e-7 1.6017365720495197e-7; -9.14199420765444e-8 -3.5830365218720865e-8 … 7.037282597221586e-8 -3.583035748485416e-8;;; -4.880642068279276e-8 -3.556492572696998e-8 … -1.1412148596697747e-8 -3.5564942096890765e-8; -1.2077841347092074e-7 -1.1202891323047923e-7 … -9.786224294202761e-8 -1.120289316727456e-7; … ; 4.711988090185771e-8 1.0980681965307195e-7 … 2.493822570142497e-7 1.0980680600132478e-7; 5.2006849158352107e-8 5.094908747404095e-8 … 4.9106256182812215e-8 5.094907276551542e-8;;; -2.50201484997855e-7 -1.7404902317933252e-7 … -3.5564950254586266e-8 -1.7404909147054354e-7; -4.793036164520203e-7 -3.6797314767597174e-7 … -1.1202893981900683e-7 -3.6797320074635805e-7; … ; 1.1356228194915098e-7 1.1211097060915574e-7 … 1.0980679712701475e-7 1.1211095679460067e-7; 1.3632641070886028e-7 1.2550868723655046e-7 … 5.0949064194866615e-8 1.2550863688575066e-7;;;;], δHtotψ = Matrix{ComplexF64}[[0.00048539971356621503 + 0.001854714687591015im -0.044117363465032464 - 0.08447000818689679im -0.7581554004795341 - 1.4498478213097141im -0.9415601123749514 + 0.7284834921132761im; 0.19228911088234216 - 0.04846118802194806im -1.483563188508676 + 0.667619216390576im -0.1682146218088955 - 0.3210161163642553im -0.2067825765866299 + 0.16058685515906107im; … ; 0.047505378753311775 - 0.012090994436158688im 0.018312008779003765 - 0.009622693428859408im 0.006879165938335011 + 0.02791576153518592im 0.009647425859586198 - 0.016070867805853605im; -0.06994692609981641 + 0.018587475374648978im -0.022710035893215307 + 0.011796420988186332im -0.07249550906990325 + 0.04942209199703852im -0.032084411926274314 - 0.06974385558614696im]], δVind = [0.0019654319470004612 0.002081000271608863 … 0.0031536662958991705 0.0020810010527418316; 0.003889131368683844 0.004779269541737564 … 0.010276673610232493 0.004779270422212057; … ; -0.006150848825766344 -0.00826960532919077 … -0.017302970534888253 -0.008269603795369989; 0.0002386819736535839 -0.00039869811965826795 … -0.0036646940523034025 -0.0003986972575115502;;; 0.0020810006958273475 0.002223227960304936 … 0.0034201212392988185 0.0022232288343991304; 0.004779270021074752 0.005734953771053686 … 0.010903600305426949 0.005734954750846453; … ; -0.008269604514758977 -0.009601010206973458 … -0.012623111781637077 -0.009601008597238192; -0.000398697655833486 -0.0010730138115941784 … -0.00423067988370088 -0.0010730128507997648;;; 0.0031536650685674594 0.003420119056756331 … 0.004832507716094614 0.003420120569865227; 0.010276672362687183 0.010903598150431102 … 0.012272144280485654 0.010903599708408889; … ; -0.017302972451905353 -0.012623114158231374 … -0.009620265119366922 -0.012623112520098134; -0.003664695471882479 -0.0042306821553082475 … -0.0054911075840294084 -0.0042306806296713365;;; … ;;; -0.002336560194007866 -0.0028154621946642863 … 0.005693320820964392 -0.0028154635939395315; -0.012913093073777906 -0.017841179354313317 … 0.011498342083524673 -0.01784118177493239; … ; -0.025174470019341408 -0.01281088886130536 … -0.00839587513322708 -0.012810887473031331; 0.00759797920119544 0.010225414453569594 … -0.005449732876474202 0.010225412512498543;;; 0.00315366481473464 0.003420118770811658 … 0.00483250726303242 0.003420120286538157; 0.010276672080319493 0.010903597930341685 … 0.012272143915872582 0.010903599363119333; … ; -0.017302972443579395 -0.01262311431196008 … -0.00962026539855183 -0.012623112783374574; -0.003664695646932366 -0.00423068235336088 … -0.005491107831518884 -0.004230680809207923;;; 0.0020810006280667244 0.0022232278796215565 … 0.0034201211130928237 0.0022232287616793475; 0.004779269954722676 0.0057349536969961154 … 0.010903600179097367 0.005734954675135696; … ; -0.008269604566302552 -0.009601010217811394 … -0.012623112019966805 -0.00960100869452057; -0.00039869772888083193 -0.0010730138937386157 … -0.004230679998561512 -0.0010730129260190012;;;;], δρ0 = [-3.6021402869511814e-7 -2.488842729874399e-7 … -4.852836744475806e-8 -2.488842775867312e-7; -5.972250196677871e-7 -4.6069781031098766e-7 … -1.1905766944752907e-7 -4.6069781542827696e-7; … ; 1.273578776470401e-7 1.0836337214300637e-7 … 4.6407159756541705e-8 1.0836337361883635e-7; 1.0352005728766926e-7 1.1957960405813408e-7 … 5.069505202373407e-8 1.1957960136315656e-7;;; -2.488842731510065e-7 -1.7310672290118804e-7 … -3.535697553161963e-8 -1.7310672588062378e-7; -4.6069781058765445e-7 -3.572170573625418e-7 … -1.1091438172046172e-7 -3.5721706081100297e-7; … ; 1.0836337256742252e-7 1.0846494746490128e-7 … 1.0862234704045265e-7 1.0846494935896047e-7; 1.1957960280458156e-7 1.1611146113967079e-7 … 5.020829608744432e-8 1.1611145924305903e-7;;; -4.852836419895823e-8 -3.535697336962606e-8 … -1.13413120471175e-8 -3.53569738781128e-8; -1.1905766319484267e-7 -1.1091437657312963e-7 … -9.770667526563206e-8 -1.1091437751505192e-7; … ; 4.640716054470566e-8 1.086223471279561e-7 … 2.484524332116455e-7 1.0862234702865635e-7; 5.0695051970192224e-8 5.020829661265072e-8 … 4.926857567333981e-8 5.020829566455859e-8;;; … ;;; 6.434571641833495e-8 4.083863188522474e-8 … -6.637928639503536e-9 4.083863255309159e-8; 1.6975037763736764e-7 6.683086880668972e-8 … -1.3340481512327457e-7 6.683086451165832e-8; … ; 3.2709967593439964e-8 1.602352222734345e-7 … 3.839014322035797e-7 1.6023521924115116e-7; -9.271367942504482e-8 -3.6336218859976205e-8 … 7.126308338704205e-8 -3.633621854464091e-8;;; -4.852836698202595e-8 -3.535697490798299e-8 … -1.1341313251661341e-8 -3.535697649679937e-8; -1.1905766615470795e-7 -1.1091437688957257e-7 … -9.770667893534349e-8 -1.1091438187993725e-7; … ; 4.640716309959557e-8 1.0862235004293466e-7 … 2.484524294855982e-7 1.0862234760484242e-7; 5.069505376858035e-8 5.020829822721957e-8 … 4.926857533335821e-8 5.0208296970107685e-8;;; -2.4888427636475477e-7 -1.7310672467418952e-7 … -3.535697645109398e-8 -1.7310672881864528e-7; -4.6069781291661596e-7 -3.572170579540555e-7 … -1.1091438262467136e-7 -3.5721706342785006e-7; … ; 1.083633749685205e-7 1.0846495074656013e-7 … 1.086223450139152e-7 1.0846495035057116e-7; 1.1957960291492097e-7 1.1611146194018703e-7 … 5.0208295913124916e-8 1.1611145926853986e-7;;;;], δeigenvalues = [[0.0004950051582409561, 0.09679860386579647, 0.023012472097985917, 0.037025908405228375]], δoccupation = [[0.0, 0.0, 0.0, 0.0]], δεF = 0.0, ε_adj = DFTK.DielectricAdjoint{Array{Float64, 4}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}, Float64, Vector{Vector{Float64}}, StaticArraysCore.SVector{3, Float64}}(Hamiltonian(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), HamiltonianBlock[DFTK.DftHamiltonianBlock(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Any[DFTK.FourierMultiplication{Float64, Vector{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.0, 0.19739208802178715, 0.7895683520871486, 1.7765287921960844, 3.1582734083485944, 4.934802200544679, 7.106115168784338, 9.67221231306757, 12.633093633394378, 15.988759129764759  …  20.133992978222288, 16.38354330580833, 13.027877809437953, 10.066996489111146, 7.5008993448279115, 5.329586376588253, 3.5530575843921683, 2.1713129682396586, 1.184352528130723, 0.5921762640653614]), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.1603906580688546 0.1602471565738079 … 0.15981840036267833 0.16024715657380792; 0.16024715657380803 0.16010321540768782 … 0.15967315066597979 0.16010321540768782; … ; 0.15981840036267841 0.1596731506659797 … 0.15923916395256382 0.1596731506659797; 0.16024715657380809 0.16010321540768785 … 0.1596731506659798 0.16010321540768785;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.1599588345026577 … 0.15952744714353234 0.1599588345026577; … ; 0.15967315066597973 0.15952744714353229 … 0.15909210945052638 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772;;; 0.1598184003626783 0.15967315066597965 … 0.15923916395256382 0.15967315066597965; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256385 0.15909210945052635 … 0.15865272241895537 0.15909210945052635; 0.1596731506659798 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; … ;;; 0.15910963392556968 0.15896221057980836 … 0.1585217272926462 0.15896221057980836; 0.15896221057980842 0.15881431920357322 … 0.15837242838914242 0.15881431920357322; … ; 0.15852172729264621 0.15837242838914237 … 0.157926332522532 0.15837242838914237; 0.15896221057980844 0.15881431920357328 … 0.1583724283891425 0.15881431920357328;;; 0.1598184003626783 0.15967315066597962 … 0.15923916395256382 0.15967315066597962; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256382 0.15909210945052635 … 0.15865272241895534 0.15909210945052635; 0.15967315066597979 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.15995883450265772 … 0.15952744714353237 0.15995883450265772; … ; 0.1596731506659797 0.15952744714353229 … 0.15909210945052635 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772]), DFTK.NonlocalOperator{Float64, Matrix{ComplexF64}, Matrix{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Matrix{ComplexF64}(undef, 7809, 0), Matrix{Float64}(undef, 0, 0)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [-0.15475965282527468 -0.1546162027006785 … -0.1541875310002283 -0.15461620270057314; -0.1546162027006379 -0.15447229889994435 … -0.15404229467728142 -0.15447229889983705; … ; -0.15418753100028507 -0.15404229467737904 … -0.15360833734818213 -0.15404229467727537; -0.15461620270060578 -0.15447229889991035 … -0.15404229467725136 -0.1544722988998063;;; -0.15461620270062262 -0.15447229889992717 … -0.1540422946772681 -0.15447229889982297; -0.15447229889988692 -0.15432794309120784 … -0.1538966028299562 -0.15432794309110187; … ; -0.1540422946773241 -0.15389660283005277 … -0.15346128701187972 -0.15389660282995019; -0.15447229889985525 -0.15432794309117442 … -0.15389660282992645 -0.1543279430910714;;; -0.15418753100032387 -0.1540422946774178 … -0.15360833734822116 -0.15404229467731426; -0.15404229467737798 -0.15389660283010811 … -0.1534612870119317 -0.15389660283000295; … ; -0.1536083373482766 -0.15346128701202735 … -0.15302191313654562 -0.15346128701192546; -0.1540422946773461 -0.15389660283007478 … -0.15346128701190198 -0.1538966028299724;;; … ;;; -0.15347880585615248 -0.15333139196606746 … -0.152890924928915 -0.15333139196595566; -0.15333139196602483 -0.15318350574397616 … -0.15274162878065578 -0.15318350574386233; … ; -0.15289092492897519 -0.1527416287807593 … -0.15229553959200254 -0.15274162878064965; -0.15333139196599008 -0.1531835057439396 … -0.15274162878062414 -0.15318350574382925;;; -0.15418753100033306 -0.15404229467743014 … -0.15360833734822568 -0.15404229467732086; -0.15404229467738842 -0.15389660283012183 … -0.15346128701193706 -0.15389660283001058; … ; -0.1536083373482843 -0.1534612870120379 … -0.15302191313654911 -0.1534612870119308; -0.1540422946773544 -0.15389660283008605 … -0.153461287011906 -0.15389660282997825;;; -0.15461620270062681 -0.15447229889993322 … -0.1540422946772704 -0.15447229889982608; -0.15447229889989217 -0.15432794309121475 … -0.15389660282995907 -0.1543279430911057; … ; -0.15404229467732788 -0.153896602830058 … -0.15346128701188158 -0.1538966028299528; -0.15447229889985908 -0.15432794309117984 … -0.1538966028299285 -0.15432794309107414]), 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.015163897687932588 -0.012550129139234832 … -0.004966868092273858 -0.012550129132105698; -0.01255012915958327 -0.010362595094462222 … -0.004135158940458497 -0.010362595105935074; … ; -0.004966868092438635 -0.004135158843898883 … -0.0020697800680422895 -0.004135158887657753; -0.012550129119179312 -0.010362595060491312 … -0.004135158827661517 -0.010362595049115143;;; -0.012550129103134609 -0.010362595051644014 … -0.004135158861893536 -0.010362595041415935; -0.010362595081313045 -0.008622039010878484 … -0.0038975113519852497 -0.008622039018126325; … ; -0.004135158854713663 -0.0038975112364698647 … -0.0034841068798180047 -0.0038975112865900457; -0.010362595023103726 -0.008622038955087558 … -0.00389751123554975 -0.008622038946749074;;; -0.004966867994462365 -0.004135158806592128 … -0.0020697800643012334 -0.004135158790026914; -0.004135158835920263 -0.0038975112713771198 … -0.0034841069339109944 -0.003897511274950395; … ; -0.002069780083301647 -0.0034841068703191495 … -0.005620579947789233 -0.0034841068750413517; -0.004135158771898132 -0.0038975112126485734 … -0.0034841068303656303 -0.003897511191870458;;; … ;;; -0.006547898161351713 -0.0048763403085934805 … -0.001486644073828944 -0.004876340253424819; -0.004876340277968705 -0.0026779890067000107 … -0.004123054339351043 -0.0026779888801952436; … ; -0.0014866439501642964 -0.004123054198302542 … -0.007170895187420022 -0.00412305415867932; -0.004876340307728178 -0.0026779890134670248 … -0.004123054155482461 -0.002677989018747458;;; -0.004966868107025449 -0.004135158876305238 … -0.0020697801664493025 -0.004135158916079621; -0.004135158879484935 -0.003897511262164789 … -0.0034841069974100095 -0.0038975113499216014; … ; -0.002069780160715882 -0.0034841069339304077 … -0.005620579905011129 -0.0034841068932732985; -0.004135158894303356 -0.0038975113114657297 … -0.0034841068411389587 -0.0038975112868485156;;; -0.01255012914645224 -0.01036259508387561 … -0.004135158896354732 -0.01036259508691446; -0.010362595089515154 -0.008622039004782221 … -0.0038975113618665374 -0.008622039034517907; … ; -0.004135158916478373 -0.0038975113134099666 … -0.0034841068426058296 -0.003897511313580316; -0.010362595080671057 -0.008622039012623114 … -0.003897511243806566 -0.00862203899316659])], 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.00953289244435266 -0.006919175266105436 … 0.000664001270176167 -0.006919175258870914; -0.006919175286413129 -0.0047316785867187496 … 0.001495697048239866 -0.004731678598084299; … ; 0.0006640012699547131 0.0014956971447017808 … 0.0035610465363393995 0.001495697101046578; -0.006919175245977002 -0.004731678552713812 … 0.001495697161066933 -0.004731678541233587;;; -0.006919175229949368 -0.004731678543883528 … 0.0014956971268180116 -0.004731678533551227; -0.004731678573512147 -0.0029911475994286382 … 0.0017333329615908835 -0.0029911476065705087; … ; 0.0014956971339419572 0.0017333330770096514 … 0.0021467155588286496 0.001733333026992055; -0.004731678515271131 -0.002991147543604267 … 0.0017333330780562204 -0.002991147535162754;;; 0.0006640013678920701 0.0014956971819697056 … 0.0035610465400414313 0.0014956971986384762; 0.0014956971526814836 0.001733333042047135 … 0.0021467155046837015 0.0017333330385790255; … ; 0.0035610465209855897 0.002146715568179845 … 1.0229334620515812e-5 0.0021467155635595337; 0.0014956972167355885 0.0017333331008090157 … 0.0021467156082588474 0.0017333331216895214;;; … ;;; -0.0009170700919345146 0.0007544783051474206 … 0.004144158289902231 0.0007544783604278814; 0.000754478335814884 0.0029528244528970564 … 0.0015077452691355988 0.002952824579515649; … ; 0.004144158413506732 0.001507745410080516 … -0.0015401022568905676 0.0015077454498134005; 0.000754478306090189 0.002952824446166652 … 0.0015077454530359052 0.002952824440996575;;; 0.0006640012553197983 0.001495697112244244 … 0.003561046437888838 0.0014956970725791352; 0.0014956971091063748 0.0017333330512457543 … 0.0021467154411793296 0.0017333329636001864; … ; 0.003561046443563639 0.00214671550455804 … 1.022937739509483e-5 0.002146715545322258; 0.0014956970943220379 0.0017333330019805907 … 0.0021467155974814944 0.0017333330267056074;;; -0.00691917527327119 -0.004731678576121175 … 0.0014956970923545119 -0.004731678579052861; -0.004731678581719502 -0.0029911475933392592 … 0.0017333329517067648 -0.0029911476229658937; … ; 0.0014956970721734441 0.0017333330000643314 … 0.0021467155960389374 0.0017333329999991757; -0.004731678572842292 -0.0029911476011452356 … 0.0017333330697973507 -0.0029911475815830174]), 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.4069020887073947e-6 + 2.678435464353936e-6im 0.005669910646374065 + 0.0031410289510296245im … 0.003854490121926139 + 0.0022357325006100355im -0.005647902716322147 - 0.003240217487638827im; 0.015741280757933337 - 0.004088188691120504im -0.018375896101985263 + 0.002084885414214953im … 0.006892052270435325 - 0.004035143232048272im -0.0117437752485086 + 0.0058225975318755225im; … ; 0.020147822919092673 - 0.005217557813477504im -0.016799339430051245 + 0.005282139576427013im … 0.013028507267078477 - 0.002310998062361622im -0.01907441161540767 + 0.0039992955246200244im; -0.015728786725065218 + 0.004147798071203848im 0.011728126968162984 - 0.0057923449337308405im … -0.012188292003121416 + 0.0009571834736201565im 0.018358340716204975 - 0.0020538757011991787im;;; -0.000441623781374039 + 0.0003279596749061963im -0.0030634909535273607 - 0.0020298635633170118im … -0.002734696033521933 - 0.0014742622725179796im 0.0035671801992688923 + 0.0017089043783437257im; -0.014799951471160467 + 0.0037328813773305192im 0.015164506588735745 - 0.002139959743737249im … -0.006034006876993016 + 0.003216192139016648im 0.010895262083714682 - 0.0045464020587000775im; … ; -0.018022425616305255 + 0.004715682650101448im 0.014753877130086863 - 0.004588225384961526im … -0.011375582617295179 + 0.002073857886461485im 0.016539873361000504 - 0.0035812832027394im; 0.015302729735383836 - 0.004143353082366989im -0.011214452010560571 + 0.0047619460360525286im … 0.010124375436186748 - 0.0010116739826146128im -0.015482829669767842 + 0.002355085787137292im;;; 0.00030489921360785674 - 0.00019875019504532481im 0.0009300254613599885 + 0.0007628586444671936im … 0.0013725353550962682 + 0.0006811499042600871im -0.001345546958724087 - 0.0005202289587547123im; 0.009333548335515909 - 0.002364556866110736im -0.00881809324458847 + 0.0015250522351308421im … 0.004042854960104291 - 0.0019346152768969435im -0.007032012233261888 + 0.002532136215975929im; … ; 0.011853054843014181 - 0.0031130212915403206im -0.009838709463678777 + 0.003001526604969435im … 0.007902642796925444 - 0.0014833657521834817im -0.010865341593563019 + 0.0024227031842288323im; -0.00974702759713628 + 0.002627727519619161im 0.007340185796100949 - 0.002702982084166741im … -0.0063444784367880785 + 0.0007833511586551631im 0.009126597299431709 - 0.0016956481433689002im;;; … ;;; 0.00011349097448437793 - 0.00015307171932767876im -0.0005426065145698948 - 0.00011658842480217582im … -0.00055753896167863 - 0.00045104293292578716im 0.00035149094119610205 + 0.0003873672475132806im; -0.004641595343264521 + 0.001324251409393436im 0.00449790745211682 - 0.0009567011696296923im … -0.0024093540493149296 + 0.0011687117538626726im 0.0037220500715934605 - 0.0013940272643725622im; … ; -0.00635160093907231 + 0.0015774321317922257im 0.0054632519047066634 - 0.001575915456977225im … -0.004770995318020335 + 0.0008391991443605865im 0.0060091592358174385 - 0.0012681469666761536im; 0.004450798782439444 - 0.0010550646651796497im -0.0035585061522156814 + 0.00115151604477im … 0.0034100157246777433 - 0.00035400489842953036im -0.004334383757591953 + 0.0007140899023474444im;;; -0.000259503664807665 + 0.00029595455418994206im 0.0013171590948868207 + 0.0004210561697780972im … 0.001187281248774151 + 0.0008466491466696612im -0.0009589769319355615 - 0.0008615398163906174im; 0.009720756922984476 - 0.0027062868720623572im -0.009107048538906619 + 0.0017809751751104266im … 0.004198121326206606 - 0.0020733922459626984im -0.007320403231465183 + 0.002787853769966396im; … ; 0.011667748764255517 - 0.0029475806297168854im -0.009683541530868114 + 0.0028628581325940853im … 0.007801413394486946 - 0.0013928632405048249im -0.010710030872357897 + 0.002284131662518937im; -0.009360422366140457 + 0.0022863030431819545im 0.0070519304333134265 - 0.0024470591891357493im … -0.006189231799692124 + 0.0006446169211604337im 0.008838313936893513 - 0.0014401921908472758im;;; 0.00038771936213090644 - 0.00039357046359088377im -0.003548923013138216 - 0.0016035948096710729im … -0.0025693633175348403 - 0.001621967266095854im 0.0030830645002000246 + 0.002134299490051117im; -0.01528510910915113 + 0.0041591873102286075im 0.015476635259558707 - 0.002416001237645741im … -0.006163828698587209 + 0.0033322382580533204im 0.011206669946715927 - 0.00482196752717514im; … ; -0.017856960005679833 + 0.00456805948303482im 0.01462395950911692 - 0.004472413373927479im … -0.011301036094615548 + 0.0020071462766102386im 0.016409824645761593 - 0.0034653558102608747im; 0.014818451068864404 - 0.0037179082400260596im -0.010902845876355734 + 0.004486369200535762im … 0.00999462942434319 - 0.0008957549575500077im -0.015171777751455355 + 0.0020799923257551417im],)])]), [8.519873384820376e-7 4.6272648497213973e-7 … 2.4177776992562027e-8 4.627264841259989e-7; 4.6272648738724174e-7 2.5014474888059483e-7 … 1.3582993923838872e-8 2.5014474976793946e-7; … ; 2.4177776995091317e-8 1.3582992927669876e-8 … 1.5659211062383528e-9 1.358299337911397e-8; 4.627264825918074e-7 2.5014474625316835e-7 … 1.3582992760157849e-8 2.50144745373283e-7;;; 4.6272648068750446e-7 2.5014474556888045e-7 … 1.3582993113313766e-8 2.501447447778108e-7; 2.5014474786358704e-7 1.3895607212413055e-7 … 1.1279986943411378e-8 1.3895607249677565e-7; … ; 1.3582993039244801e-8 1.127998589457562e-8 … 7.937752472285001e-9 1.12799863496474e-8; 2.5014474336148003e-7 1.3895606925568735e-7 … 1.1279985886221212e-8 1.3895606882697195e-7;;; 2.417777549175683e-8 1.3582992542792787e-8 … 1.5659210974526541e-9 1.35829923718944e-8; 1.3582992845360422e-8 1.12799862115215e-8 … 7.937752858112024e-9 1.1279986243965559e-8; … ; 1.5659211420703355e-9 7.937752404531557e-9 … 3.5718556575856674e-8 7.937752438211071e-9; 1.3582992184869901e-8 1.1279985678290767e-8 … 7.937752119555041e-9 1.1279985489631628e-8;;; … ;;; 5.791843100006389e-8 2.281603112146438e-8 … 5.620368918961917e-10 2.281603030811356e-8; 2.2816030669964866e-8 3.4892801554946926e-9 … 1.3458509968328938e-8 3.4892796418063083e-9; … ; 5.620367474637727e-10 1.3458508522373108e-8 … 7.728314624288465e-8 1.345850811617799e-8; 2.2816031108709883e-8 3.489280182970717e-9 … 1.3458508083405005e-8 3.489280204413635e-9;;; 2.417777721890612e-8 1.3582993261993959e-8 … 1.5659213373201082e-9 1.3582993672328738e-8; 1.3582993294799682e-8 1.1279986127877133e-8 … 7.937753311029858e-9 1.1279986924671679e-8; … ; 1.5659213238573775e-9 7.937752858249005e-9 … 3.571855571676457e-8 7.937752568256162e-9; 1.3582993447672724e-8 1.127998657551072e-8 … 7.937752196399026e-9 1.1279986351994831e-8;;; 4.6272648582876676e-7 2.5014474806179085e-7 … 1.3582993468839105e-8 2.5014474829681586e-7; 2.501447484979642e-7 1.3895607181069956e-7 … 1.1279987033127326e-8 1.3895607333953108e-7; … ; 1.3582993676442352e-8 1.12799865931595e-8 … 7.937752206861283e-9 1.127998659470761e-8; 2.5014474781393376e-7 1.3895607221382956e-7 … 1.1279985961188993e-8 1.3895607121348818e-7;;;;], Matrix{ComplexF64}[[0.06595774992873271 + 0.25202489068606115im -0.44677981300601605 - 0.8554340488025821im 4.747448387216451e-6 - 2.094590733773433e-6im 9.102935462107696e-6 + 1.2579695019788496e-5im; -0.04475094193813239 - 0.1709935717816952im -0.01968688198875269 - 0.037693800945132805im -0.4578183108194919 + 0.23933017017279043im 0.23133833325702954 + 0.2991877146460759im; … ; 0.014217708645738911 + 0.054325935457546715im 0.007120486141792157 + 0.01363335193272926im -0.005517282085845693 + 0.007805389502985617im 0.00045280456393618216 + 0.006600295598872903im; -0.025492448962491608 - 0.09740677437952747im -0.014106234465517061 - 0.027008725194495643im 0.0059801012249373646 - 0.030451574968832266im 0.009926357563758146 - 0.020646539334027547im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944095310946, [[-0.5541104220550804, -0.010678397007108771, 0.14495669188132615, 0.14495724828586318]], 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.620694128621567e-7, -6.293616532988726e-7, -2.2951208412740465e-7, 4.994026735549289e-7, 7.923313409224658e-7, 5.851073326181161e-7, 6.516971826399777e-7, 1.018930927607833e-6, 3.6385903520349244e-7, -2.2010939802897998e-6  …  2.78101884204095e-6, 1.201695891837675e-6, 2.4456535829463852e-8, 1.0078553675162869e-7, 5.5353085664096e-7, 4.5341729037545303e-7, 4.6318894678299865e-9, -1.0437055487142915e-7, 1.1211095679460067e-7, 1.2550863688575066e-7], resid_history = [0.24939209356054567, 0.0037664998521441256, 0.0002851433508296768, 4.682374740254381e-6, 3.4770430404878206e-8, 5.342101655083252e-11, 3.926916157648659e-8, 1.1267074367739535e-9], converged = true, n_iter = 8, residual_norm = 1.1267074367739535e-9, maxiter = 100, tol = 1.0e-8, s = 0.9353684800229439, residual = [5.582210699199959e-7, 3.697932688954942e-7, 7.3251223708572e-8, -9.293459923715711e-8, -8.82342981047293e-8, -3.9474392567417206e-8, -2.537839856022625e-8, -1.746594926419902e-8, 2.8671044661845296e-10, -3.1587878798566694e-8  …  7.680716640855302e-7, 3.3587197986394955e-7, 7.041273467732312e-9, 3.0198502811960137e-8, 1.7333293344706696e-7, 1.5098658214386036e-7, 1.7250889569901261e-9, -4.805646571282736e-8, 7.776017613831698e-8, 2.44767384999368e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.478233291149013e-7, -3.895907644099683e-8, 1.1267074367739535e-9, -0.00020436567490577157, 4.687936576887906e-6, -3.51516064192728e-8, 5.342101655083252e-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}}([[1.329815808155449e-7, 8.975171310052382e-7, 2.7874478305851564e-7, -3.543314882677803e-7, -2.2688996306855992e-7, -2.5427463489419616e-8, 6.075134824228102e-8, 2.443812069130355e-7, 1.5730847987448593e-7, -1.3173668867259752e-6  …  1.3000047339092936e-6, 4.2957346855859746e-7, 4.900397331743088e-9, 9.055599912839974e-9, 2.377927069828585e-8, -9.128226790404141e-9, -8.496971797368109e-10, 4.860831765500522e-8, -9.479360107655859e-8, -1.984893243590437e-7], [1.392748622445856e-7, -1.1085976554190104e-6, -3.5561702340615855e-7, 4.18206930345036e-7, 1.8837331213450005e-7, -6.015072609846688e-8, -1.9525156655106753e-7, -5.296555550217216e-7, -2.923537215483265e-7, 2.286363138803886e-6  …  -1.74972725406501e-6, -5.268560194821888e-7, -4.144780447056799e-9, 1.9578062423581344e-9, 5.761021442594791e-8, 1.015873018164382e-7, 2.420504477134441e-9, -1.0819769362732161e-7, 1.9908001089934233e-7, 4.610777499465032e-7]]), H = [1.1755128126945384 0.00015628139143408152 … 0.0 0.0; 0.11342490812784702 1.0118220143423482 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.180972303905911 0.09733465333297822 … 0.0 0.0; 0.0 1.0075442781889778 … 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.9257125271817925
Interacting polarizability:     1.7736548575393751

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

Manual solution of the Dyson equations

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

using KrylovKit

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

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

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

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

[ Info: GMRES linsolve starts with norm of residual = 4.19e+00
[ Info: GMRES linsolve in iteration 1; step 1: normres = 2.49e-01
[ Info: GMRES linsolve in iteration 1; step 2: normres = 3.77e-03
[ Info: GMRES linsolve in iteration 1; step 3: normres = 2.85e-04
[ Info: GMRES linsolve in iteration 1; step 4: normres = 4.69e-06
[ Info: GMRES linsolve in iteration 1; step 5: normres = 1.09e-08
[ Info: GMRES linsolve in iteration 1; step 6: normres = 6.27e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 1.43e-08
[ Info: GMRES linsolve in iteration 2; step 1: normres = 1.11e-09
[ Info: GMRES linsolve in iteration 2; step 2: normres = 3.25e-11
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 4.00e-12
└ * number of operations = 12
Non-interacting polarizability: 1.92571252718173
Interacting polarizability:     1.7736548565753363

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