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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457369536366235
Displaced dipole:  0.017612220814714412
Polarizability :   1.7746794510078072

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.0001496724660075052 - 0.003292493983731345im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.30894843935531424 - 0.016645026849794988im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.041014937987898833 - 0.002030198394749048im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.08086818419911301 + 0.003365809263059055im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206944016473276e-7 -2.5020155524827386e-7 … -4.880641647192035e-8 -2.5020146730892793e-7; -6.293616414315573e-7 -4.793036557390818e-7 … -1.207784344048352e-7 -4.793035914565236e-7; … ; 1.3626260545913106e-7 1.1356229380408188e-7 … 4.7119923778608454e-8 1.1356230853288762e-7; 1.3305472216922953e-7 1.3632634670582277e-7 … 5.200687174452448e-8 1.3632641055697867e-7;;; -2.5020148380674804e-7 -1.7404909128736362e-7 … -3.55649287230468e-8 -1.7404902715582457e-7; -4.793036035960695e-7 -3.6797319995374315e-7 … -1.1202894095561407e-7 -3.6797315252595344e-7; … ; 1.1356230551988708e-7 1.121109855670089e-7 … 1.0980686029961482e-7 1.1211099914278657e-7; 1.3632639854906768e-7 1.2550863627729439e-7 … 5.094910781821761e-8 1.2550868360042614e-7;;; -4.8806423899380695e-8 -3.556494616041935e-8 … -1.1412146546872104e-8 -3.5564931542479137e-8; -1.207784400522291e-7 -1.1202895664655195e-7 … -9.786224672891052e-8 -1.1202894357695461e-7; … ; 4.711992179772884e-8 1.0980684686596622e-7 … 2.493823467694333e-7 1.0980685837815589e-7; 5.200686619289786e-8 5.0949092004439126e-8 … 4.9106301682672527e-8 5.094910538526167e-8;;; … ;;; 6.475763107707572e-8 4.110335688799187e-8 … -6.681701867111032e-9 4.1103338198317104e-8; 1.689755672123565e-7 6.653284278304535e-8 … -1.32944520443899e-7 6.65328353490788e-8; … ; 3.267766380875319e-8 1.6017368206463468e-7 … 3.847036619214623e-7 1.6017369660802746e-7; -9.141988858591221e-8 -3.583032211033407e-8 … 7.037285835626493e-8 -3.5830329191735046e-8;;; -4.880644867960145e-8 -3.556496460247245e-8 … -1.1412152749892178e-8 -3.556495000520085e-8; -1.2077845837287062e-7 -1.1202897326554e-7 … -9.78622602248506e-8 -1.1202895999419272e-7; … ; 4.7119915649106656e-8 1.0980683310693304e-7 … 2.4938231700361796e-7 1.0980684424437498e-7; 5.200684811116946e-8 5.0949075641565183e-8 … 4.910628797212889e-8 5.0949088940821375e-8;;; -2.50201538676783e-7 -1.7404913124916276e-7 … -3.556493795073999e-8 -1.7404906721735035e-7; -4.79303643510176e-7 -3.679732295658574e-7 … -1.1202894915300941e-7 -3.6797318199486066e-7; … ; 1.1356229668749029e-7 1.1211097767980066e-7 … 1.0980685333502052e-7 1.1211099096690948e-7; 1.3632635869969678e-7 1.2550860693391371e-7 … 5.0949099682575357e-8 1.2550865412913497e-7;;;;], δHtotψ = Matrix{ComplexF64}[[-8.706253856480443e-5 + 0.001915202241803406im 0.00042656929282124513 - 0.09529607305683345im 0.19547851623364249 - 0.7316016053466236im 0.6760877326326311 + 0.699034212942552im; 0.1980071554188126 + 0.010805336451508859im -1.623640322686106 - 0.10230947235268854im 0.04301340621023282 - 0.16194557171457435im 0.15002104439639768 + 0.15430123143463306im; … ; 0.04895322989656348 + 0.002556254305343178im 0.020686308257177245 + 4.0639124490130615e-5im -0.009894169676067283 + 0.016969062671004634im -0.006408690910823912 - 0.0015196274907050369im; -0.07231168279110028 - 0.0030145076605742624im -0.025590481726628384 - 0.00017187900271725654im -0.0971798534895005 + 0.0400715974885468im 0.012070143922814117 + 0.11139405102923688im]], δVind = [0.0019654319793672246 0.0020810010577417106 … 0.0031536634491801732 0.002081000314742409; 0.00388913086447326 0.004779269812757575 … 0.010276668696026893 0.004779268981644526; … ; -0.006150846468141489 -0.008269600872676912 … -0.017302965573025665 -0.008269602285201814; 0.00023868248082289955 -0.0003986967148660145 … -0.003664695584213718 -0.00039869754499413396;;; 0.0020810004552188786 0.0022232285267306607 … 0.0034201174778602786 0.002223227694915675; 0.004779269139574708 0.0057349536989527525 … 0.010903594591668153 0.005734952772048645; … ; -0.008269602025576427 -0.009601005918705927 … -0.012623111577540987 -0.009601007369678956; -0.00039869739059920114 -0.0010730125678534768 … -0.004230681901403032 -0.001073013492214583;;; 0.003153663934825179 0.003420119173618557 … 0.004832500306062225 0.0034201177473385785; 0.010276669229047424 0.01090359635278685 … 0.012272137361300369 0.010903594872201801; … ; -0.017302964746475726 -0.012623109767632914 … -0.00962026639316645 -0.012623111294222791; -0.003664695053907667 -0.004230680166120153 … -0.005491110070151985 -0.004230681627580211;;; … ;;; -0.002336561544469723 -0.0028154652275764647 … 0.005693315307501369 -0.0028154639699547283; -0.012913095102514961 -0.017841186665956837 … 0.011498338843435107 -0.0178411848901277; … ; -0.025174449248614006 -0.012810882834919044 … -0.008395874664249106 -0.012810884214545577; 0.00759797794771703 0.010225412178092247 … -0.005449732867638555 0.010225413976016007;;; 0.0031536655363294556 0.003420120944406313 … 0.00483250292540612 0.0034201195238162765; 0.01027667100050369 0.010903598172857664 … 0.012272139312193495 0.010903596709463309; … ; -0.017302962130108637 -0.012623107813904291 … -0.009620264934559208 -0.012623109352884156; -0.0036646932531778113 -0.004230678313759964 … -0.005491108116942825 -0.004230679773524351;;; 0.002081000916701067 0.00222322904170736 … 0.003420118363522934 0.002223228211785619; 0.004779269653348232 0.005734954269133462 … 0.010903595516345618 0.005734953345921715; … ; -0.008269601128373836 -0.009601004992340473 … -0.012623110620787682 -0.009601006452205652; -0.00039869686938697607 -0.0010730119885267585 … -0.004230680979602124 -0.0010730129112984535;;;;], δρ0 = [-3.6021405581015113e-7 -2.488842983638021e-7 … -4.8528380385018494e-8 -2.488842977978008e-7; -5.97225084510115e-7 -4.606978705925631e-7 … -1.1905771095472612e-7 -4.606978700322507e-7; … ; 1.2735792015237132e-7 1.0836341712648912e-7 … 4.640720835512553e-8 1.083634172862554e-7; 1.0352008408137646e-7 1.1957962734743183e-7 … 5.0695074019797915e-8 1.1957962751693956e-7;;; -2.48884297716898e-7 -1.731067453418354e-7 … -3.5356988118451595e-8 -1.7310674494028263e-7; -4.6069787005256326e-7 -3.572171129560286e-7 … -1.1091442337499236e-7 -3.5721711269723964e-7; … ; 1.0836341704428603e-7 1.084649953294493e-7 … 1.086223988076953e-7 1.0846499567749624e-7; 1.1957962756634416e-7 1.1611148390180197e-7 … 5.020831886288916e-8 1.1611148408247984e-7;;; -4.852838036347224e-8 -3.535698820054518e-8 … -1.1341323553108181e-8 -3.5356988154912704e-8; -1.1905771104654104e-7 -1.1091442331803078e-7 … -9.77067180447766e-8 -1.1091442347388434e-7; … ; 4.640720823192467e-8 1.0862239845234437e-7 … 2.4845249266084274e-7 1.0862239897006096e-7; 5.0695074007212574e-8 5.020831878623738e-8 … 4.92686005204127e-8 5.020831894585875e-8;;; … ;;; 6.434570658204698e-8 4.0838622083541064e-8 … -6.637938485970738e-9 4.083862196243938e-8; 1.6975033718655048e-7 6.683082498893845e-8 … -1.334048532031141e-7 6.683082480638283e-8; … ; 3.271002155122388e-8 1.6023527901771565e-7 … 3.8390150030143725e-7 1.6023527933587846e-7; -9.271365843080212e-8 -3.633619545570837e-8 … 7.12631117449529e-8 -3.6336195278985066e-8;;; -4.8528380777711346e-8 -3.535698850222853e-8 … -1.1341323678926802e-8 -3.5356988467264706e-8; -1.1905771124271688e-7 -1.109144235138783e-7 … -9.770671790184807e-8 -1.1091442345571984e-7; … ; 4.6407208166377117e-8 1.0862239835423296e-7 … 2.4845249170330587e-7 1.0862239849577363e-7; 5.069507409026862e-8 5.0208318888689966e-8 … 4.92686004655412e-8 5.02083189689804e-8;;; -2.488842984393054e-7 -1.7310674586155765e-7 … -3.5356988316016207e-8 -1.7310674550495334e-7; -4.606978705488857e-7 -3.5721711341936943e-7 … -1.1091442338255308e-7 -3.5721711297651837e-7; … ; 1.083634172337455e-7 1.0846499563159787e-7 … 1.0862239866816506e-7 1.0846499570041348e-7; 1.1957962724283726e-7 1.161114837889687e-7 … 5.020831892895996e-8 1.1611148386817164e-7;;;;], δeigenvalues = [[0.0004950053137440722, 0.09679859346860589, 0.046574080937973494, 0.04299679981339548]], δ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.16039065806885477 0.16024715657380806 … 0.15981840036267847 0.16024715657380803; 0.16024715657380792 0.1601032154076877 … 0.15967315066597973 0.16010321540768768; … ; 0.15981840036267841 0.15967315066597976 … 0.15923916395256388 0.15967315066597976; 0.16024715657380792 0.1601032154076877 … 0.15967315066597973 0.16010321540768768;;; 0.16024715657380803 0.16010321540768782 … 0.15967315066597979 0.1601032154076878; 0.1601032154076877 0.1599588345026576 … 0.1595274471435323 0.15995883450265758; … ; 0.15967315066597973 0.15952744714353234 … 0.1590921094505264 0.15952744714353234; 0.1601032154076877 0.1599588345026576 … 0.1595274471435323 0.15995883450265758;;; 0.15981840036267841 0.15967315066597979 … 0.1592391639525639 0.15967315066597976; 0.15967315066597967 0.1595274471435323 … 0.15909210945052638 0.15952744714353229; … ; 0.15923916395256385 0.15909210945052643 … 0.15865272241895542 0.15909210945052643; 0.15967315066597967 0.1595274471435323 … 0.15909210945052638 0.15952744714353229;;; … ;;; 0.15910963392556976 0.15896221057980847 … 0.15852172729264627 0.15896221057980844; 0.15896221057980836 0.1588143192035732 … 0.1583724283891424 0.15881431920357317; … ; 0.15852172729264624 0.15837242838914245 … 0.15792633252253205 0.15837242838914245; 0.15896221057980836 0.1588143192035732 … 0.15837242838914237 0.15881431920357317;;; 0.15981840036267841 0.15967315066597979 … 0.1592391639525639 0.15967315066597976; 0.15967315066597967 0.1595274471435323 … 0.1590921094505264 0.15952744714353229; … ; 0.15923916395256385 0.1590921094505264 … 0.15865272241895542 0.1590921094505264; 0.15967315066597967 0.1595274471435323 … 0.15909210945052638 0.15952744714353229;;; 0.16024715657380803 0.16010321540768782 … 0.1596731506659798 0.1601032154076878; 0.1601032154076877 0.15995883450265758 … 0.15952744714353229 0.15995883450265755; … ; 0.1596731506659797 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; 0.1601032154076877 0.1599588345026576 … 0.1595274471435323 0.1599588345026576]), 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.1547596528221812 -0.1546162026976636 … -0.15418753099699867 -0.1546162026974089; -0.15461620269749005 -0.15447229889687267 … -0.15404229467399494 -0.15447229889661815; … ; -0.1541875309973329 -0.15404229467450556 … -0.15360833734507795 -0.15404229467424618; -0.1546162026975788 -0.15447229889696254 … -0.15404229467408287 -0.15447229889670627;;; -0.15461620269748022 -0.1544722988968623 … -0.15404229467398578 -0.15447229889660868; -0.15447229889669017 -0.15432794308808634 … -0.15389660282661755 -0.15432794308783296; … ; -0.15404229467431726 -0.153896602827124 … -0.15346128700871903 -0.15389660282686574; -0.15447229889677724 -0.15432794308817457 … -0.15389660282670375 -0.15432794308791942;;; -0.1541875309971395 -0.15404229467430897 … -0.15360833734488744 -0.15404229467405484; -0.1540422946741372 -0.1538966028269412 … -0.15346128700854114 -0.15389660282668727; … ; -0.15360833734521837 -0.15346128700904677 … -0.15302191313332886 -0.15346128700878783; -0.15404229467422306 -0.15389660282702833 … -0.15346128700862616 -0.15389660282677253;;; … ;;; -0.15347880585324186 -0.15333139196323353 … -0.15289092492584497 -0.15333139196296536; -0.15333139196304924 -0.153183505741076 … -0.15274162877752237 -0.153183505740808; … ; -0.15289092492620027 -0.15274162877806538 … -0.15229553958906342 -0.15274162877779202; -0.15333139196314596 -0.15318350574117398 … -0.15274162877761846 -0.1531835057409041;;; -0.15418753099736085 -0.15404229467453498 … -0.15360833734510368 -0.1540422946742729; -0.15404229467435535 -0.15389660282716405 … -0.15346128700875408 -0.15389660282690215; … ; -0.15360833734544987 -0.1534612870092831 … -0.1530219131335554 -0.15346128700901618; -0.15404229467444885 -0.15389660282725884 … -0.15346128700884698 -0.15389660282699508;;; -0.15461620269759146 -0.1544722988969759 … -0.15404229467409453 -0.15447229889671832; -0.15447229889679986 -0.15432794308819844 … -0.15389660282672468 -0.1543279430879411; … ; -0.15404229467443367 -0.15389660282724285 … -0.153461287008833 -0.15389660282698048; -0.15447229889689076 -0.15432794308829048 … -0.1538966028268148 -0.1543279430880313]), 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.015163898082546298 -0.012550129579739456 … -0.004966868888283622 -0.012550129575731046; -0.012550129585071672 -0.010362595578014784 … -0.004135159764042741 -0.010362595575287835; … ; -0.004966868895105948 -0.004135159758518388 … -0.0020697813842313734 -0.004135159761780271; -0.012550129576006066 -0.010362595569590821 … -0.0041351597511382135 -0.010362595566942348;;; -0.012550129573734387 -0.010362595567045902 … -0.004135159752837231 -0.01036259556506931; -0.010362595575103509 -0.00862203955723001 … -0.003897512197933259 -0.008622039557668645; … ; -0.004135159756013768 -0.003897512183918832 … -0.003484107822593036 -0.0038975121909222365; -0.010362595562771361 -0.008622039546054523 … -0.0038975121804759078 -0.008622039545509172;;; -0.004966868888180915 -0.004135159752178976 … -0.0020697813879092227 -0.0041351597553463985; -0.00413515976663213 -0.003897512193746488 … -0.0034841078338730654 -0.0038975121992346315; … ; -0.0020697813833329927 -0.0034841078143935033 … -0.005620580627502389 -0.0034841078239854096; -0.004135159747624179 -0.0038975121770159714 … -0.0034841078178545026 -0.003897512182138153;;; … ;;; -0.006547897437147471 -0.0048763394488651225 … -0.001486645517293664 -0.0048763394437427955; -0.004876339443337187 -0.0026779878183282955 … -0.004123055042677765 -0.0026779878144083907; … ; -0.0014866454850237824 -0.0041230550211619425 … -0.007170895787623102 -0.004123055026387042; -0.004876339448516639 -0.0026779878274176463 … -0.004123055036910993 -0.0026779878192071775;;; -0.00496686890092834 -0.004135159763888027 … -0.0020697813932429475 -0.004135159765007456; -0.004135159768050967 -0.0038975121958863468 … -0.003484107827443292 -0.0038975121955996286; … ; -0.0020697813830267446 -0.003484107816546827 … -0.005620580621713682 -0.0034841078193002905; -0.004135159764264463 -0.0038975121918208753 … -0.0034841078241408486 -0.0038975121945945784;;; -0.012550129581581852 -0.010362595574992917 … -0.004135159758924407 -0.01036259557176769; -0.01036259557844069 -0.008622039561552036 … -0.0038975121950547633 -0.00862203955846395; … ; -0.00413515976278369 -0.0038975121920502083 … -0.0034841078202908 -0.0038975121930971447; -0.010362595572954169 -0.008622039556369143 … -0.003897512188078392 -0.008622039554388013])], 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.009532892835872733 -0.006919175703595 … 0.0006640004773961804 -0.0069191756993319045; -0.006919175708753798 -0.004731679067199741 … 0.0014956962279420493 -0.004731679064218302; … ; 0.0006640004702395666 0.0014956962329558104 … 0.0035610452232545547 0.0014956962299533033; -0.006919175699776955 -0.004731679058865651 … 0.0014956962407586469 -0.004731679055960938;;; -0.0069191756974065766 -0.004731679056220368 … 0.0014956962391567737 -0.0047316790539902; -0.004731679064105973 -0.002991148142658745 … 0.0017333321189815093 -0.0029911481428440273; … ; 0.0014956962356486962 0.0017333321324895078 … 0.00214671461921434 0.001733332125744369; -0.004731679051860895 -0.0029911481315714923 … 0.0017333321363526516 -0.002991148130771013;;; 0.000664000477358 0.0014956962394918434 … 0.0035610452197672475 0.0014956962365785227; 0.0014956962252103568 0.0017333321228446223 … 0.0021467146081121683 0.001733332117610387; … ; 0.0035610452240124922 0.0021467146270861626 … 1.0228658124175963e-5 0.002146714617753188; 0.0014956962441324325 0.0017333321394880142 … 0.002146714624045716 0.0017333321346216;;; … ;;; -0.0009170693648195725 0.0007544791677098171 … 0.004144156849507637 0.0007544791731002908; 0.0007544791734219389 0.0029528256441688962 … 0.0015077445689422597 0.002952825648356781; … ; 0.004144156881422193 0.0015077445899151273 … -0.001540102854154471 0.001507744584963392; 0.0007544791681457583 0.0029528256349815683 … 0.0015077445746129138 0.0029528256434618767;;; 0.0006640004643892247 0.0014956962275567785 … 0.003561045214217279 0.0014956962266993896; 0.0014956962235733616 0.0017333321204819141 … 0.0021467146143290287 0.0017333321210305061; … ; 0.003561045224087231 0.0021467146246964726 … 1.0228663686341274e-5 0.0021467146222099343; 0.0014956962272663571 0.0017333321244526003 … 0.0021467146175385464 0.0017333321219426307;;; -0.006919175705365287 -0.004731679064280986 … 0.0014956962329608792 -0.0047316790607982145; -0.004731679067552844 -0.002991148147092903 … 0.0017333321217528406 -0.0029911481437474956; … ; 0.0014956962287623402 0.00173333212423931 … 0.0021467146214026117 0.001733332123454747; -0.004731679062157223 -0.0029911481420020203 … 0.0017333321286391173 -0.002991148139761708]), 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 = [-3.345989049939237e-8 + 2.53453787291149e-6im -0.007606803939953575 + 0.0019440969266971905im … -0.005190308900710656 + 0.0014081804129758833im 0.007608988735059226 - 0.0020248764813521193im; 0.02752585481697611 + 0.001278974636886397im -0.01769367727218372 - 0.0021580542639648703im … 0.015684885468476316 - 0.00042118426662799167im -0.02660978274016719 + 0.00016803414415376644im; … ; 0.028480095772834054 + 0.0013510510279069838im -0.025321836261697667 - 0.0007356291462162977im … 0.012958140558475502 + 0.0011607914329334653im -0.02226169479641158 - 0.0015340901540442607im; -0.027524829975026893 - 0.00122391896796762im 0.02660708859281256 - 0.00014337679491102553im … -0.008555763469133356 - 0.0015236704609449907im 0.01769171298089728 + 0.00218236250071477im;;; 0.004777087711410928 - 0.0025461341113880214im 0.001657183358599735 + 0.0003575578988695039im … 0.004525383088279227 - 0.0014858322896093218im -0.007258006984998893 + 0.0026830222690744842im; -0.024950669736132337 + 0.000472663512069795im 0.016807917853421193 + 0.0005859867444494942im … -0.01353917302044775 + 0.0006734325388648175im 0.022548012399411545 - 0.0009114079817356862im; … ; -0.02283119892317286 - 0.001639200217903764im 0.020467881154653506 + 0.001043277128377304im … -0.010857938711269381 - 0.0012259206874327883im 0.01806582793449528 + 0.0016700803480183392im; 0.019351031994633536 + 0.0024885152564426936im -0.01894063389397346 - 0.0010405461983935514im … 0.006708504083954634 + 0.0015946472634934346im -0.013201001869742986 - 0.002537959240676935im;;; -0.0032594521849103305 + 0.0017635423040411814im 0.0007090763778425054 - 0.0008268154180043352im … -0.0027886094192063437 + 0.0010574010934111923im 0.003769055390975851 - 0.001625080766651731im; 0.014359230175172802 - 0.0006647007107647132im -0.010594572577427542 + 0.00014651078637807315im … 0.008724417381347644 - 0.00060043738672129im -0.012996695699951148 + 0.0007730947505019775im; … ; 0.01359499193648261 + 0.001278143984689889im -0.012515398927070858 - 0.0009166827308731052im … 0.007245276989579701 + 0.0009929854311811833im -0.011134768883600755 - 0.0012769714500641578im; -0.009881515778753638 - 0.001768053242121173im 0.009653007165693244 + 0.0010679699154866247im … -0.004049382726216863 - 0.001181820228205678im 0.0072510281728986855 + 0.0016947395143475116im;;; … ;;; -0.0014358746615275894 + 0.0007085571274550926im 0.001813158019252606 - 0.00074934070229806im … 3.001319156886754e-5 + 0.00016458677249112658im 0.0006106723233302798 - 0.000435513384890377im; -0.004374220272149794 - 0.000847065367783083im 0.003368726371140312 + 0.0008658797832113434im … -0.00357930702227305 - 0.00030383470302822917im 0.004412260458899529 + 0.0005935517062836703im; … ; -0.008663557359427608 + 2.507631401242529e-5im 0.008172559822260199 - 0.00011141579039336961im … -0.0053411602066538045 - 0.00014616911669842005im 0.007438190755281315 + 8.019739825348224e-5im; 0.006798038368006517 - 0.00033899606137446805im -0.006498141937653169 + 0.00041715416894842663im … 0.003519032530148585 + 1.904794111110154e-5im -0.005454557778825078 + 0.00014481968586458466im;;; 0.003254480550103369 - 0.0016531028050272947im -0.003763616700305409 + 0.0015191460738413274im … -0.0006394495882781376 - 6.93718112857265e-5im -0.0007037662845643629 + 0.0007204936520007117im; 0.009885939831236076 + 0.0016804277573336218im -0.007255858422485224 - 0.0016039683424668553im … 0.0069230383194753674 + 0.0003437453472808146im -0.009657807074760237 - 0.000977086626870747im; … ; 0.01574394035042795 + 0.00015051417954236354im -0.014316475696094831 + 2.8300802634979893e-5im … 0.00842064883802121 + 0.0003764088879157513im -0.01293583713044745 - 0.00033196610368130047im; -0.01435406527199256 + 0.0005786817672566907im 0.012991318339908707 - 0.0006834332233653349im … -0.005850540624550475 - 0.00023722143698492814im 0.010589362142417687 - 5.656303627046334e-5im;;; -0.0047743570057886845 + 0.0024645707134015593im 0.007253607032305787 - 0.002578308495471632im … 0.0026042664579756384 - 0.00047862431695060994im -0.0016614466908805902 - 0.0002521643323261208im; -0.019353245811094204 - 0.0024617234027541277im 0.013204348740756513 + 0.002475329819919674im … -0.01203094677895623 - 0.00011707101625175332im 0.018944273623215546 + 0.0009775410952868266im; … ; -0.024752073355786015 - 0.0006308028316789471im 0.021975743651724482 + 0.0002519024526260282im … -0.011724699073384477 - 0.0007711696347312596im 0.01957364676836642 + 0.0008786661647413873im; 0.024947226455653435 - 0.00044873187281591694im -0.02254350666061158 + 0.000850177980080997im … 0.008216450654970254 + 0.0008037548022120007im -0.016803864963250682 - 0.0006474162890179615im],)])]), [8.519874101906337e-7 4.627265372544453e-7 … 2.417778920640775e-8 4.627265367786901e-7; 4.627265378872958e-7 2.5014478628026524e-7 … 1.3583002420422855e-8 2.5014478606934774e-7; … ; 2.4177789311087556e-8 1.3583002363431047e-8 … 1.5659241969563852e-9 1.3583002397085074e-8; 4.6272653681132685e-7 2.50144785628721e-7 … 1.3583002287292993e-8 2.5014478542387586e-7;;; 4.627265365417184e-7 2.501447854318847e-7 … 1.3583002304824153e-8 2.501447852790173e-7; 2.501447860551011e-7 1.389561002143158e-7 … 1.1279994624277798e-8 1.3895610023686623e-7; … ; 1.3583002337590975e-8 1.1279994497032639e-8 … 7.937759196798645e-9 1.1279994560624978e-8; 2.50144785101289e-7 1.3895609963974387e-7 … 1.1279994465775904e-8 1.3895609961169533e-7;;; 2.417778920483303e-8 1.3583002298033296e-8 … 1.5659242055935882e-9 1.3583002330709397e-8; 1.3583002447139523e-8 1.1279994586265752e-8 … 7.937759277253677e-9 1.1279994636098308e-8; … ; 1.5659241948466429e-9 7.937759138310417e-9 … 3.5718570226226214e-8 7.937759206729401e-9; 1.3583002251039256e-8 1.1279994434360195e-8 … 7.937759162997626e-9 1.1279994480866765e-8;;; … ;;; 5.7918410690558554e-8 2.2816018446515992e-8 … 5.620385777660674e-10 2.2816018370994928e-8; 2.281601836501197e-8 3.489275329974625e-9 … 1.3458517178471806e-8 3.489275314056605e-9; … ; 5.62038540077406e-10 1.3458516957902108e-8 … 7.728316679247412e-8 1.345851701146836e-8; 2.2816018441377877e-8 3.4892753668842266e-9 … 1.3458517119354857e-8 3.4892753335426655e-9;;; 2.4177789400427155e-8 1.3583002418826999e-8 … 1.5659242181184167e-9 1.3583002430379708e-8; 1.358300246177653e-8 1.1279994605692577e-8 … 7.93775923139262e-9 1.1279994603094004e-8; … ; 1.5659241941279067e-9 7.937759153673233e-9 … 3.571857010997485e-8 7.937759173309476e-9; 1.358300242271283e-8 1.1279994568784039e-8 … 7.937759207835069e-9 1.1279994593964904e-8;;; 4.6272653747310946e-7 2.501447860465482e-7 … 1.3583002367620659e-8 2.5014478579708853e-7; 2.50144786313211e-7 1.389561004365295e-7 … 1.1279994598144175e-8 1.3895610027775516e-7; … ; 1.3583002407436538e-8 1.1279994570865841e-8 … 7.93775918037387e-9 1.1279994580368283e-8; 2.501447858888557e-7 1.389561001700497e-7 … 1.127999453480135e-8 1.3895610006819392e-7;;;;], Matrix{ComplexF64}[[-0.011830350459565814 + 0.2602441413704729im 0.004319906826375916 - 0.9650704376433991im -1.7640553966198345e-7 + 2.4324793991994046e-6im 4.803274390733027e-7 - 2.182957185806007e-6im; 0.008026643229373302 - 0.17657016003870712im 0.00019035313397392414 - 0.0425247943949815im -0.23120101567614573 - 0.06190782645322442im 0.22128563544644175 - 0.21337722945614418im; … ; -0.0025501245238516985 + 0.05609765918457968im -6.884774602240572e-5 + 0.015380668660243266im -0.002115241027712815 + 0.003832974664605215im 0.009335783299455986 - 0.004120170082052694im; 0.004572390746279503 - 0.10058348715628561im 0.00013639310813450282 - 0.03047028304610094im -0.0006764503307284989 - 0.024654626362009446im -0.040023418113434316 + 0.01137623927707857im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944067155256, [[-0.5541104220653278, -0.010678391365723471, 0.14495659821467524, 0.1449566199260838]], 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.6206944016473276e-7, -6.293616414315573e-7, -2.295121136942102e-7, 4.994025393823447e-7, 7.923311477126468e-7, 5.85107166160578e-7, 6.516970399335604e-7, 1.0189307887749155e-6, 3.638591132312953e-7, -2.201093322643534e-6  …  2.7810181121663054e-6, 1.2016955098283806e-6, 2.4456429805641953e-8, 1.0078550147057368e-7, 5.535308262467282e-7, 4.5341731526535795e-7, 4.631965997408036e-9, -1.0437048556835826e-7, 1.1211099096690948e-7, 1.2550865412913497e-7], resid_history = [0.24939209630534534, 0.0037665490162731956, 0.00028528104996069863, 4.6803082971746775e-6, 6.352902536582029e-8, 4.281534148019413e-11, 4.6221166840841455e-8, 1.2128675432694028e-9], converged = true, n_iter = 8, residual_norm = 1.2128675432694028e-9, maxiter = 100, tol = 1.0e-8, s = 0.935658605989957, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.2135637876872867e-7, -4.5677977044401285e-8, 1.2128675432694028e-9, -0.00020446478374706693, 4.682645253803098e-6, -6.411189074179985e-8, 4.281534148019413e-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}}([[-8.841741415844185e-8, -1.8183512903753918e-6, -6.508984129969282e-7, 1.0475456216289511e-6, 1.2482101744565565e-6, 8.949138335764374e-7, 1.0996692483791363e-6, 1.6355903716531545e-6, 4.3506164097524674e-7, -1.8716633242383535e-6  …  2.8262302495452454e-6, 1.097000208413521e-6, 2.6938322321329227e-8, 1.3674470566958933e-7, 7.510544615137911e-7, 5.098708247950408e-7, 4.6051044229268235e-9, -1.308014214979867e-7, 2.224833375896107e-7, 4.6694088204157475e-7], [4.4642871309530793e-7, 1.563698553006022e-6, 5.315411968957501e-7, -8.645140491972833e-7, -1.0402034213817938e-6, -7.274024331489933e-7, -8.695421994381352e-7, -1.292199899779122e-6, -3.5302327038342945e-7, 1.562342617817628e-6  …  -1.6115622226051467e-6, -5.970589545617183e-7, -1.4979979379567083e-8, -7.858863366557008e-8, -4.2761244745329957e-7, -2.686577352129673e-7, -2.1210482928932924e-9, 5.761713968321495e-8, -9.72008970929311e-8, -1.3655751155960454e-7]]), H = [1.0981138844595142 0.011695967013092896 … 0.0 0.0; 0.1453034483999878 1.0215565959044652 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1076855128418406 0.14560017082480115 … 0.0 0.0; 0.0 1.0115432852242512 … 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.925712528804124
Interacting polarizability:     1.7736548588043435

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.28e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 1.47e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.61e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 1.35e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125288040082
Interacting polarizability:     1.7736548514856294

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