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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457351817738536
Displaced dipole:  0.017612223482183775
Polarizability :   1.774679700036116

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.002043649858895685 + 0.0025858099568802275im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.24111911907802272 + 0.19387565047307198im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.032114962205700065 + 0.025592500068971017im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.06369223720486188 - 0.04994287628082848im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206944963558417e-7 -2.5020149344415525e-7 … -4.8806445823699076e-8 -2.502015415204659e-7; -6.293616529919869e-7 -4.793036129247023e-7 … -1.2077845164548465e-7 -4.793036483060312e-7; … ; 1.362626005217756e-7 1.1356229985766626e-7 … 4.711991331759172e-8 1.1356229384785804e-7; 1.3305472086555614e-7 1.3632639438209699e-7 … 5.200684629202653e-8 1.3632636007442427e-7;;; -2.5020154456040154e-7 -1.740490853959549e-7 … -3.556495940295816e-8 -1.74049120652331e-7; -4.793036490450718e-7 -3.679731956468649e-7 … -1.1202896551006439e-7 -3.679732222613391e-7; … ; 1.1356229134527004e-7 1.1211098350492528e-7 … 1.0980683571394391e-7 1.121109776408999e-7; 1.363263562124869e-7 1.255086410644279e-7 … 5.094907632071343e-8 1.2550861547389677e-7;;; -4.880644656853952e-8 -3.556495155806907e-8 … -1.1412158412823873e-8 -3.556495944439894e-8; -1.207784503577636e-7 -1.1202895694374842e-7 … -9.786227374488792e-8 -1.120289642681737e-7; … ; 4.711990908733336e-8 1.0980683832279848e-7 … 2.4938228432666233e-7 1.0980683267296906e-7; 5.200684335794968e-8 5.094908210209789e-8 … 4.910626799887311e-8 5.094907465953819e-8;;; … ;;; 6.47575817621355e-8 4.110331212276221e-8 … -6.6816985836548745e-9 4.1103323087479235e-8; 1.6897555069493426e-7 6.653283533301556e-8 … -1.3294451103338255e-7 6.653284188403755e-8; … ; 3.26776608351967e-8 1.6017371802269696e-7 … 3.8470367736248134e-7 1.6017370730379124e-7; -9.14199319415699e-8 -3.5830346875209394e-8 … 7.037285942470678e-8 -3.583034430846868e-8;;; -4.88064239281344e-8 -3.556493515058816e-8 … -1.1412153043286533e-8 -3.55649425820173e-8; -1.2077843805741306e-7 -1.120289467046966e-7 … -9.786226137371104e-8 -1.1202895181734854e-7; … ; 4.711991883816427e-8 1.0980685557968333e-7 … 2.4938231479052087e-7 1.0980684916781526e-7; 5.2006864851617076e-8 5.0949101695264844e-8 … 4.9106283354305554e-8 5.094909393033878e-8;;; -2.502014905704263e-7 -1.7404904645895584e-7 … -3.55649505879716e-8 -1.7404908114647082e-7; -4.793036122766712e-7 -3.6797316956984876e-7 … -1.1202895757879666e-7 -3.679731945956554e-7; … ; 1.1356230211060495e-7 1.1211099294665953e-7 … 1.0980684310049888e-7 1.12110987224103e-7; 1.3632639816737992e-7 1.255086722269713e-7 … 5.0949085435699685e-8 1.2550864675627496e-7;;;;], δHtotψ = Matrix{ComplexF64}[[0.001188765148431528 - 0.001504132905429854im 0.024786962296416887 - 0.09201699104956489im 1.0295457657563143 + 0.18434583233461865im 0.6526222716632946 - 0.3312792602018435im; -0.15445465404018316 - 0.12436779829363585im -1.5434625904641206 - 0.5141967674938639im 0.22797165755904272 + 0.04032840543674128im 0.14307957206355246 - 0.07312954040724592im; … ; -0.03825298266181418 - 0.03065391441781911im 0.019987778632834213 + 0.00533037818726905im -0.01789702396478541 + 0.007296388384451789im -0.00752175733093603 - 0.002209226485308992im; 0.056950269764068244 + 0.0446624396520621im -0.024695276397768612 - 0.006711622355493257im -0.058816403808501054 + 0.1011408610605425im 0.05123615789398882 - 0.07168984838752386im]], δVind = [0.0019654320678083835 0.0020810005795624583 … 0.0031536656644984007 0.002081000987608346; 0.003889131048196411 0.004779269374438765 … 0.010276671550214314 0.004779269843657327; … ; -0.006150846699920111 -0.008269602290167509 … -0.017302962965764473 -0.00826960146091466; 0.00023868250578509115 -0.000398697336155623 … -0.003664693469674991 -0.0003986968860921823;;; 0.0020810010202288766 0.0022232285113982543 … 0.0034201208447667995 0.002223228971783093; 0.004779269873028705 0.005734953797319418 … 0.010903598447093196 0.0057349543229240005; … ; -0.008269601438844901 -0.009601006437190378 … -0.012623108461300086 -0.009601005594039331; -0.0003986968449207068 -0.0010730126851323783 … -0.004230678800607363 -0.0010730121817634425;;; 0.003153665786689519 0.0034201201180268505 … 0.004832506821898248 0.0034201209190871545; 0.010276671703283237 0.010903597702767227 … 0.012272142447957532 0.010903598547874867; … ; -0.017302963289182696 -0.012623109339246014 … -0.009620263247328082 -0.012623108484733084; -0.0036646933721410447 -0.004230679555532398 … -0.005491105841308446 -0.0042306787616963595;;; … ;;; -0.0023365595685055133 -0.0028154618963680098 … 0.005693315524798805 -0.0028154625743024055; -0.01291309262849748 -0.0178411813648253 … 0.011498339067926721 -0.01784118211734156; … ; -0.025174456825099507 -0.012810886878942433 … -0.008395875175259221 -0.012810886212547517; 0.007597980187489696 0.010225416527076186 … -0.005449733493410404 0.010225415406206988;;; 0.0031536640989391854 0.003420118242200249 … 0.004832504001312894 0.003420119043074551; 0.010276669701290216 0.010903595616639668 … 0.012272140399311838 0.010903596523425575; … ; -0.01730296544782974 -0.01262311119837492 … -0.00962026474092498 -0.01262311037071294; -0.0036646951082386593 -0.004230681336057684 … -0.005491107761220431 -0.004230680555759589;;; 0.0020810005455138988 0.0022232279828665255 … 0.003420119922606208 0.0022232284387232668; 0.004779269340515063 0.005734953202367067 … 0.010903597501940506 0.005734953727548452; … ; -0.008269602311234458 -0.009601007348139477 … -0.012623109451829564 -0.009601006501806956; -0.00039869737411096524 -0.0010730132703099416 … -0.004230679716031576 -0.0010730127694230315;;;;], δρ0 = [-3.6021403735454094e-7 -2.488842834822279e-7 … -4.852837474549302e-8 -2.4888428413558246e-7; -5.972250725712862e-7 -4.6069785995551797e-7 … -1.1905770179663526e-7 -4.606978606748283e-7; … ; 1.2735791911749103e-7 1.0836341486637764e-7 … 4.6407206412844023e-8 1.0836341666589922e-7; 1.0352010034217185e-7 1.1957963895507997e-7 … 5.069507426809743e-8 1.1957963908389915e-7;;; -2.488842832157922e-7 -1.7310673378593953e-7 … -3.535698361469045e-8 -1.7310673436237817e-7; -4.6069785879269314e-7 -3.572171032073065e-7 … -1.1091441572646058e-7 -3.5721710410835197e-7; … ; 1.0836341489956278e-7 1.084649939030712e-7 … 1.0862240016364019e-7 1.0846499519887839e-7; 1.1957963828595967e-7 1.161114916392605e-7 … 5.0208318952204936e-8 1.1611149164984779e-7;;; -4.8528373795902425e-8 -3.535698309314686e-8 … -1.1341320559305475e-8 -3.535698303708741e-8; -1.1905769935553972e-7 -1.109144138672739e-7 … -9.770671143156333e-8 -1.1091441397675023e-7; … ; 4.6407202678164976e-8 1.0862239729092599e-7 … 2.4845249954756365e-7 1.0862239764947436e-7; 5.069507269662625e-8 5.020831814019061e-8 … 4.9268599964815384e-8 5.020831790486882e-8;;; … ;;; 6.434570601737412e-8 4.083862213172683e-8 … -6.637935475749492e-9 4.083862299026115e-8; 1.6975035076916578e-7 6.683083490272734e-8 … -1.3340484467484688e-7 6.68308375124894e-8; … ; 3.271001369101233e-8 1.6023528013713078e-7 … 3.83901511321887e-7 1.602352771772076e-7; -9.271367098981626e-8 -3.63362026464593e-8 … 7.126311122439589e-8 -3.6336204046129964e-8;;; -4.85283760908956e-8 -3.5356985217348564e-8 … -1.1341321387377194e-8 -3.535698473673852e-8; -1.1905770546374491e-7 -1.1091442031162122e-7 … -9.770671294965531e-8 -1.1091441824610831e-7; … ; 4.640720618536375e-8 1.0862240063662589e-7 … 2.4845250065509803e-7 1.0862240019959966e-7; 5.069507554666904e-8 5.020832089301527e-8 … 4.926860139104892e-8 5.020832030461217e-8;;; -2.4888428481514323e-7 -1.731067352295022e-7 … -3.535698402219259e-8 -1.7310673536905758e-7; -4.6069786216243706e-7 -3.572171067552813e-7 … -1.1091441612908503e-7 -3.572171061724402e-7; … ; 1.0836341640452984e-7 1.0846499499523944e-7 … 1.0862240059356007e-7 1.0846499638008856e-7; 1.1957963960353342e-7 1.1611149286071307e-7 … 5.0208319678356925e-8 1.1611149286575472e-7;;;;], δeigenvalues = [[0.0004950046928189986, 0.0967985718836501, 0.04101130545771005, 0.04645452058398241]], δ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.15475965279917514 -0.15461620267470647 … -0.15418753097387886 -0.1546162026743482; -0.15461620267464868 -0.15447229887408295 … -0.15404229465103814 -0.15447229887372066; … ; -0.15418753097400273 -0.15404229465122093 … -0.15360833732163895 -0.15404229465086436; -0.15461620267440884 -0.1544722988738396 … -0.1540422946508016 -0.15447229887348318;;; -0.15461620267401033 -0.15447229887343925 … -0.15404229465039843 -0.15447229887308367; -0.15447229887338204 -0.15432794306482794 … -0.15389660280319015 -0.1543279430644682; … ; -0.15404229465052133 -0.15389660280337125 … -0.15346128698481024 -0.15389660280301762; -0.15447229887314376 -0.15432794306458594 … -0.1538966028029553 -0.1543279430642324;;; -0.15418753097321117 -0.15404229465042615 … -0.15360833732083876 -0.15404229465007108; -0.1540422946503692 -0.15389660280322146 … -0.15346128698465034 -0.15389660280286208; … ; -0.1536083373209616 -0.15346128698483105 … -0.15302191310895497 -0.15346128698447825; -0.15404229465013114 -0.15389660280297948 … -0.15346128698441586 -0.15389660280262663;;; … ;;; -0.15347880583160525 -0.15333139194165496 … -0.15289092490410894 -0.15333139194127648; -0.15333139194159312 -0.15318350571968045 … -0.15274162875596436 -0.15318350571929779; … ; -0.15289092490424058 -0.15274162875615901 … -0.15229553956700534 -0.1527416287557818; -0.15333139194134124 -0.15318350571942474 … -0.1527416287557153 -0.1531835057190482;;; -0.1541875309752763 -0.154042294652505 … -0.1536083373229151 -0.15404229465213534; -0.15404229465244487 -0.1538966028053108 … -0.15346128698673756 -0.15389660280493708; … ; -0.15360833732304344 -0.15346128698692724 … -0.15302191311104746 -0.15346128698655884; -0.15404229465219832 -0.15389660280506062 … -0.15346128698649397 -0.15389660280469278;;; -0.15461620267504866 -0.1544722988744847 … -0.15404229465144287 -0.15447229887412178; -0.1544722988744259 -0.15432794306587874 … -0.15389660280424008 -0.15432794306551184; … ; -0.15404229465156855 -0.15389660280442571 … -0.1534612869858629 -0.15389660280406423; -0.15447229887418334 -0.15432794306563266 … -0.1538966028040006 -0.15432794306527156]), 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.01516389800508557 -0.012550129495291891 … -0.004966868684268351 -0.012550129491753246; -0.012550129492594432 -0.010362595488520383 … -0.004135159571184296 -0.010362595480908055; … ; -0.004966868741668568 -0.004135159607110429 … -0.0020697812919488277 -0.004135159649716507; -0.012550129506194808 -0.010362595489743536 … -0.004135159604536237 -0.010362595502408496;;; -0.012550129473816006 -0.010362595467901612 … -0.004135159546074701 -0.010362595461923628; -0.010362595466207923 -0.008622039457051119 … -0.003897512030144569 -0.0086220394536058; … ; -0.004135159611314215 -0.003897512068587917 … -0.0034841077824499652 -0.0038975120972179404; -0.010362595478586312 -0.008622039463407772 … -0.003897512050342273 -0.008622039469698327;;; -0.004966868631682749 -0.004135159526578763 … -0.0020697810985331323 -0.004135159507165407; -0.004135159525814774 -0.003897512000759624 … -0.003484107681298683 -0.0038975119953723724; … ; -0.0020697811884812243 -0.003484107728664173 … -0.005620580624065506 -0.0034841077350505944; -0.004135159533364977 -0.0038975120127380056 … -0.0034841076926308734 -0.003897512003853404;;; … ;;; -0.006547897709253449 -0.004876339682866692 … -0.0014866451445176845 -0.004876339734483083; -0.00487633971015832 -0.002677988084757916 … -0.004123054874090842 -0.002677988150585752; … ; -0.0014866452282905466 -0.004123054951464722 … -0.0071708957578688086 -0.004123054904908341; -0.004876339710407078 -0.0026779880791345675 … -0.004123054884194647 -0.002677988148927687;;; -0.004966868781162469 -0.0041351596827179414 … -0.00206978120498848 -0.004135159638924663; -0.004135159657553795 -0.0038975121405734284 … -0.0034841077179350757 -0.003897512090460095; … ; -0.002069781285763359 -0.00348410778816268 … -0.005620580629881568 -0.0034841077804672726; -0.004135159666947181 -0.003897512135755075 … -0.0034841077485399875 -0.003897512111067293;;; -0.01255012952073327 -0.010362595515334014 … -0.004135159582471784 -0.010362595503522618; -0.010362595507520911 -0.008622039503679515 … -0.0038975120403297154 -0.008622039483648127; … ; -0.004135159641479149 -0.0038975120878078163 … -0.003484107787445006 -0.0038975121191833577; -0.010362595519440819 -0.008622039500380636 … -0.003897512079572495 -0.008622039506419784])], 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.009532892735406102 -0.006919175596190465 … 0.0006640007045311241 -0.006919175592293523; -0.00691917559343508 -0.004731678954915514 … 0.0014956964437573524 -0.004731678946940893; … ; 0.0006640006470071177 0.0014956964076483382 … 0.0035610453389760457 0.0014956963653988357; -0.006919175606795565 -0.004731678955895279 … 0.0014956964106419722 -0.004731678968203829;;; -0.0069191755740184695 -0.004731678933653203 … 0.001495696469506519 -0.004731678927319615; -0.004731678931902144 -0.0029911480192213753 … 0.001733332310197626 -0.0029911480154163166; … ; 0.0014956964041441872 0.001733332271573117 … 0.0021467146832661746 0.0017333322432967276; -0.0047316789440422224 -0.0029911480253359997 … 0.001733332290234845 -0.002991148031273004;;; 0.000664000757784387 0.00149569648897473 … 0.0035610455331919343 0.0014956965087431625; 0.0014956964897957564 0.0017333323395512905 … 0.0021467147845773845 0.001733332345297921; … ; 0.003561045443121024 0.0021467147370311224 … 1.0228685934894968e-5 0.002146714730997502; 0.0014956964824836957 0.0017333323278148818 … 0.0021467147734797287 0.0017333323370523402;;; … ;;; -0.000917069615289025 0.0007544789552867085 … 0.004144157244019564 0.0007544789040487927; 0.0007544789280569762 0.0029528253991348554 … 0.001507744759087215 0.0029528253336896858; … ; 0.004144157160115085 0.0015077446815186298 … -0.0015401028023421552 0.001507744728452237; 0.0007544789280601273 0.0029528254050139716 … 0.001507744749232572 0.002952825335597384;;; 0.0006640006062395409 0.0014956963307566863 … 0.0035610454246602476 0.001495696374919613; 0.0014956963559810627 0.0017333321976481294 … 0.0021467147458537725 0.0017333322481351918; … ; 0.003561045343757027 0.0021467146754364313 … 1.0228678026312325e-5 0.0021467146835002383; 0.0014956963468342847 0.0017333322027166718 … 0.002146714715492499 0.0017333322277722982;;; -0.006919175621974069 -0.004731678982131046 … 0.0014956964320649944 -0.004731678969956719; -0.004731678974258991 -0.0029911480669005414 … 0.0017333322989625693 -0.0029911480465022525; … ; 0.0014956963729320077 0.0017333322512987559 … 0.002146714677218448 0.001733332220284703; -0.004731678985936315 -0.002991148063355582 … 0.0017333322599593203 -0.0029911480690336303]), 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.6319720535619662e-6 - 1.8783517041011398e-6im -0.002922137232303737 - 0.0025825469464472233im … -0.0019426327381058938 - 0.0018259909590054694im 0.002872094364352342 + 0.002646679851763197im; -0.02120771495613228 - 0.016799829874263844im 0.01887067888397588 + 0.015092921240999063im … -0.008107332992042263 - 0.006202457214476045im 0.015473981131729219 + 0.012025878607368569im; … ; -0.02207547204717162 - 0.017513511873005705im 0.01789236912524908 + 0.014137368384163039im … -0.012059776474184127 - 0.009660910436977909im 0.019058806449862332 + 0.015191856556664543im; 0.021253753979027786 + 0.01675891993668819im -0.01546638942331127 - 0.012044595701881704im … 0.010771745288458541 + 0.00871840544276808im -0.01886256505932171 - 0.015111892539413543im;;; -0.0011721974927561703 + 0.008845881416013154im 0.0024013918894902637 - 0.003674447244696512im … 0.0010955231481700946 + 0.0030413490118176003im -0.000994925551930093 - 0.006742127208070357im; 0.017843226854149612 + 0.008391802835545727im -0.015432412018809724 - 0.008614893872095193im … 0.007082149264826468 + 0.0038921351503442953im -0.013244973579443101 - 0.006638443628157942im; … ; 0.01824550806331231 + 0.016440430603275977im -0.014921879354275311 - 0.013337681875545213im … 0.010338598612732642 + 0.009170130940411234im -0.015837647582250643 - 0.014165494192354322im; -0.016493209894782077 - 0.018746380714357757im 0.012348298368209724 + 0.01334153564057012im … -0.008694860620708449 - 0.008577023000487367im 0.014535501987490025 + 0.015317944458192354im;;; 0.0008159004770920708 - 0.00605412064230176im -0.001151823470037425 + 0.003642538813115619im … -0.0003667824928619774 - 0.002610131098730453im 1.4769367202263552e-5 + 0.004697288446347321im; -0.010001222665357655 - 0.0032981572232309996im 0.008773151937572301 + 0.003531810108000059im … -0.004675330109884572 - 0.0017448072504120304im 0.00785738244436889 + 0.0027038100948441183im; … ; -0.011149498086543617 - 0.011057248712071938im 0.00941628429299895 + 0.009301177889387035im … -0.007031528835881614 - 0.006844684604519023im 0.009942648747536809 + 0.009777127774092963im; 0.008877775526469771 + 0.011622955551888298im -0.007003831626897723 - 0.008937568210557193im … 0.005284604423929545 + 0.006126965697870407im -0.007919527735002325 - 0.009765434728701113im;;; … ;;; 0.00031611441565179184 - 0.002625378536528394im -3.461256578821199e-5 + 0.002419506494974381im … 0.0004866665732790581 - 0.0011735450097731387im -0.0004928409061150935 + 0.0020053798205529im; 0.004091591871882453 + 0.005655229729903692im -0.003821535060909318 - 0.00512135119475058im … 0.0023350163810057873 + 0.003230988380671055im -0.0034239762288761095 - 0.004761974430119211im; … ; 0.006298221165688797 + 0.003383064435173234im -0.005536618163048431 - 0.002941133377308934im … 0.004505987783438031 + 0.0025401858539748083im -0.005816715214433332 - 0.003194274605242903im; -0.004620036200028494 - 0.0012296956428326509im 0.0038728545967524454 + 0.0009597344253661508im … -0.0032014340664078492 - 0.0011445643450579412im 0.004270675999265408 + 0.0013192700117471423im;;; -0.000740982903107916 + 0.0059776253283399024im -8.226632017711074e-5 - 0.0046216572113695325im … -0.0008810324085000299 + 0.0013624301694548377im 0.0010841125286124966 - 0.003567659766483498im; -0.008932096672016342 - 0.011562571957983957im 0.00797450978607391 + 0.009701116961446446im … -0.004244312740864806 - 0.005074034038199678im 0.007059062810188381 + 0.008873643066680245im; … ; -0.011664006253993227 - 0.0070845181544455455im 0.009847412948076017 + 0.005971633152774641im … -0.007312909590647456 - 0.004672027636411386im 0.010373829441002414 + 0.006447602810955503im; 0.009947294958039303 + 0.0033579357432417375im -0.007802367751620334 - 0.0027676803460362585im … 0.0057156332267416215 + 0.0027975878344627223im -0.008717908487568439 - 0.0035952731330323837im;;; 0.0011094686453887448 - 0.00879107554302191im 0.001064245259147795 + 0.006664684365796092im … 0.0015548865663940993 - 0.0005108703543950328im -0.0023322653830115878 + 0.003598222928110517im; 0.016506697959556536 + 0.018731342933275177im -0.014571213020232858 - 0.015273569299106126im … 0.006721461846153587 + 0.006680229507876501im -0.012383973126753715 - 0.013297831604973206im; … ; 0.018705096079448487 + 0.012888225450000875im -0.015282808331352607 - 0.010549524199938314im … 0.010545935529064812 + 0.007567514647526816im -0.016198304765886806 - 0.011377234850137266im; -0.017830677137777796 - 0.00840574486342403im 0.013209618443340326 + 0.006682120754263633im … -0.009055379413734545 - 0.005788747946959619im 0.015396569546369426 + 0.008657948385280015im],)])]), [8.519873961145749e-7 4.627265272315892e-7 … 2.4177786076028676e-8 4.6272652681160094e-7; 4.627265269114269e-7 2.50144779358446e-7 … 1.3583000430784655e-8 2.50144778769688e-7; … ; 2.4177786956771255e-8 1.358300080142017e-8 … 1.565923980255984e-9 1.3583001240967334e-8; 4.627265285256188e-7 2.5014477945304887e-7 … 1.3583000774858054e-8 2.5014478043260134e-7;;; 4.6272652468267224e-7 2.501447777637172e-7 … 1.3583000171737179e-8 2.501447773013626e-7; 2.50144777632723e-7 1.3895609506371103e-7 … 1.1279993100826357e-8 1.389560948865703e-7; … ; 1.3583000844784462e-8 1.1279993449874379e-8 … 7.9377589104707e-9 1.127999370982276e-8; 2.501447785901107e-7 1.3895609539053086e-7 … 1.127999328421121e-8 1.389560957139479e-7;;; 2.41777852691668e-8 1.3582999970603027e-8 … 1.5659235260703667e-9 1.3582999770324538e-8; 1.3582999962723058e-8 1.127999283402404e-8 … 7.93775818898798e-9 1.1279992785107978e-8; … ; 1.5659237372892996e-9 7.937758526831531e-9 … 3.571857015720468e-8 7.937758572384564e-9; 1.3583000040613661e-8 1.12799929427792e-8 … 7.937758269816033e-9 1.1279992862110602e-8;;; … ;;; 5.791841832146883e-8 2.281602189639107e-8 … 5.620381423884736e-10 2.281602265737376e-8; 2.2816022298749554e-8 3.4892764118412045e-9 … 1.3458515450203035e-8 3.4892766791426834e-9; … ; 5.620382402287424e-10 1.3458516243400312e-8 … 7.728316577375167e-8 1.3458515766129913e-8; 2.2816022302415476e-8 3.4892763890070996e-9 … 1.3458515553782773e-8 3.489276672410769e-9;;; 2.417778756275758e-8 1.3583001581429974e-8 … 1.565923776053699e-9 1.358300112963158e-8; 1.3583001321821372e-8 1.127999410347347e-8 … 7.937758450306829e-9 1.1279993648467489e-8; … ; 1.5659239657306343e-9 7.93775895121467e-9 … 3.57185702740032e-8 7.937758896328432e-9; 1.358300141873151e-8 1.1279994059726036e-8 … 7.937758668601174e-9 1.1279993835569395e-8;;; 4.627265302511583e-7 2.501447814323142e-7 … 1.3583000547227785e-8 2.501447805187746e-7; 2.501447808280178e-7 1.3895609746106462e-7 … 1.1279993193302e-8 1.3895609643116467e-7; … ; 1.3583001155985033e-8 1.1279993624383101e-8 … 7.937758946097697e-9 1.1279993909262238e-8; 2.501447817499458e-7 1.389560972914608e-7 … 1.1279993549610915e-8 1.389560976019532e-7;;;;], Matrix{ComplexF64}[[0.16153343770061024 - 0.20438669059244002im 0.2510194517770099 - 0.9318631056173453im -1.0151283593777154e-6 + 9.96481315099383e-6im 1.3292530471076212e-6 - 8.984424900501638e-6im; -0.10959702961500994 + 0.13867205771034063im 0.0110609114062081 - 0.041061567252333814im 0.058605108294519946 - 0.3252509726852907im -0.10524035300648718 - 0.20777517981488147im; … ; 0.034819795143106735 - 0.0440571487898115im -0.004000584924074114 + 0.014851428590783176im 0.00592344174462786 - 0.0025114306796592597im -0.0053198115821744395 - 0.0062389256452762346im; -0.06243213127283946 + 0.0789947696625772im 0.007925467584059494 - 0.02942182692930724im -0.029845229821150544 - 0.003496534064700936im 0.023809252612611756 + 0.023368550153006994im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944074271533, [[-0.5541104183003411, -0.010678396553965495, 0.1449566538563966, 0.1449574295134834]], 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.6206944963558417e-7, -6.293616529919869e-7, -2.2951210243856136e-7, 4.994026037702256e-7, 7.92331278104313e-7, 5.85107335965742e-7, 6.516971974327513e-7, 1.0189308925544904e-6, 3.6385911101467483e-7, -2.201093498514542e-6  …  2.7810183706958605e-6, 1.2016956713479164e-6, 2.4456484079956677e-8, 1.0078547851849846e-7, 5.535307523785669e-7, 4.534172265764083e-7, 4.631906147723457e-9, -1.0437050596864155e-7, 1.12110987224103e-7, 1.2550864675627496e-7], resid_history = [0.24939209826314732, 0.003766548946086814, 0.000285279619812531, 4.68035105917777e-6, 4.333160693476078e-8, 5.962101384796747e-11, 3.836930258469233e-8, 9.56273696252421e-10], converged = true, n_iter = 8, residual_norm = 9.56273696252421e-10, maxiter = 100, tol = 1.0e-8, s = 0.935655221451883, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.0602649753393907e-7, -3.7603229961143326e-8, 9.56273696252421e-10, -0.00020446384472413035, 4.6826790536163395e-6, -4.375554238411569e-8, 5.962101384796747e-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.7158478448791546e-7, -2.024556475762929e-6, -7.174036448395366e-7, 1.1543716818781536e-6, 1.3772254640428783e-6, 9.85192465683872e-7, 1.2072723773347473e-6, 1.7973470042168798e-6, 4.812454260339672e-7, -2.0913062273015624e-6  …  3.130141460927302e-6, 1.2183651373598098e-6, 2.977028428411418e-8, 1.504123269825033e-7, 8.257681536375138e-7, 5.622358598869874e-7, 5.081663678826598e-9, -1.428319849532449e-7, 2.3990970392820593e-7, 4.93896300405046e-7], [1.3153785128876743e-6, 2.112007210844192e-6, 6.457079568391798e-7, -1.045891192035443e-6, -1.277276435540771e-6, -8.73258143142321e-7, -1.0056839546259035e-6, -1.489857743274655e-6, -4.2790052947125676e-7, 2.04579362036329e-6  …  -1.1477360019636294e-6, -3.9037061300556153e-7, -9.994921227044838e-9, -5.375099204336607e-8, -2.742997442372636e-7, -1.317472283508459e-7, -4.5070374119307214e-10, 3.958223957456159e-9, 4.1678938578471595e-9, 2.175160678898427e-7]]), H = [1.0685813713009964 0.014701106340132726 … 0.0 0.0; 0.12646674942274808 1.028595375993464 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.0760390261515025 0.1354899204419901 … 0.0 0.0; 0.0 1.0200555703280212 … 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.925712549441308
Interacting polarizability:     1.7736548768108873

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.07e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.20e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 1.40e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125494411924
Interacting polarizability:     1.7736548702991637

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