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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.000134575147629413
Displaced dipole:  0.017612224118337863
Polarizability :   1.7746799265967277

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.0026526492353688095 + 0.001956110611919117im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.18571109322841367 - 0.24746229334125985im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.02450518057298269 - 0.03295213124906242im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.04778689618890386 + 0.06532537353817289im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206936080073187e-7 -2.5020148521455667e-7 … -4.880642534626172e-8 -2.5020143390789837e-7; -6.293616355120162e-7 -4.793036435745483e-7 … -1.20778458088666e-7 -4.793035982576895e-7; … ; 1.362626453626837e-7 1.1356233574113315e-7 … 4.711993847838679e-8 1.1356232760902526e-7; 1.3305483220270281e-7 1.3632643617983086e-7 … 5.200688801587797e-8 1.3632647051664633e-7;;; -2.5020144746794643e-7 -1.7404905528320275e-7 … -3.5564941051435916e-8 -1.7404901891872305e-7; -4.793036143244381e-7 -3.6797320021831135e-7 … -1.1202896477710751e-7 -3.679731690264681e-7; … ; 1.1356234991543458e-7 1.1211103942044544e-7 … 1.098068935384654e-7 1.1211103056090211e-7; 1.3632646896810072e-7 1.2550869878864477e-7 … 5.094912640621202e-8 1.2550872290989106e-7;;; -4.8806430749961355e-8 -3.556494875673768e-8 … -1.1412159612802077e-8 -3.5564942698977235e-8; -1.2077845916747654e-7 -1.120289652768624e-7 … -9.786227227065497e-8 -1.1202896391011505e-7; … ; 4.711997163311872e-8 1.0980690723535186e-7 … 2.4938240403979267e-7 1.0980690960500457e-7; 5.2006892816613826e-8 5.094912018753948e-8 … 4.91063260894387e-8 5.09491285976058e-8;;; … ;;; 6.475758055068855e-8 4.110331406587604e-8 … -6.681717762917916e-9 4.1103302815513486e-8; 1.689754758907648e-7 6.653276840885413e-8 … -1.3294456371341592e-7 6.653276122173779e-8; … ; 3.267776494805288e-8 1.601737961436373e-7 … 3.847037658496372e-7 1.6017380149092463e-7; -9.141985205872903e-8 -3.5830279962103396e-8 … 7.03729048829368e-8 -3.5830285477537686e-8;;; -4.880644110434897e-8 -3.556495592611559e-8 … -1.141215926901821e-8 -3.5564948899284594e-8; -1.2077846048401832e-7 -1.1202896847252112e-7 … -9.786227283126609e-8 -1.1202896392289498e-7; … ; 4.711994295344706e-8 1.098068716966724e-7 … 2.493823676254928e-7 1.098068766534171e-7; 5.200687558440497e-8 5.0949101806222236e-8 … 4.9106308992026216e-8 5.094911124046173e-8;;; -2.502014695713384e-7 -1.7404907030591742e-7 … -3.5564942165824926e-8 -1.7404903378328142e-7; -4.793036189609506e-7 -3.679732019586368e-7 … -1.1202896295673528e-7 -3.679731710292754e-7; … ; 1.1356231863023022e-7 1.1211100356061307e-7 … 1.0980687205714964e-7 1.121110024604581e-7; 1.363264396787704e-7 1.255086721529596e-7 … 5.094911434520007e-8 1.255086995969166e-7;;;;], δHψ = Matrix{ComplexF64}[[-0.0015430126972570163 - 0.0011378449800318621im -0.0718213032642675 + 0.06263570479471743im 0.41786979901636 + 1.5720731229198022im -0.8250033238559891 - 0.2750248740494im; -0.11913769664666286 + 0.15852381126655246im 0.9958286886343172 + 1.2864681809252037im 0.09218476444231757 + 0.3477076066540677im -0.18358482162290032 - 0.06050493011065877im; … ; -0.029358684586634726 + 0.039255837284382575im -0.01363559911445534 - 0.015556229422865778im -0.004188798645160325 - 0.020129041240519122im 0.017945341600751966 - 0.004931596168256781im; 0.042734671814557715 - 0.05841073980762024im 0.01677691667339176 + 0.019324493825457457im 0.011989922265774432 + 0.047432076048862265im 0.08946249973116159 - 0.08017231986347925im]], δVind = [0.0019654317867564748 0.0020810007015092382 … 0.0031536635536445336 0.002081000280760611; 0.0038891314665029193 0.004779270316302038 … 0.010276670004366856 0.00477926987042108; … ; -0.006150849938906261 -0.008269604609381307 … -0.01730296900403344 -0.008269606044711168; 0.00023868139453795865 -0.0003986980925922367 … -0.0036646966993474635 -0.00039869857631230545;;; 0.0020810003578371084 0.0022232282457054775 … 0.0034201177006180556 0.002223227782342176; 0.004779269966297532 0.005734954419593584 … 0.010903596189845976 0.00573495392121607; … ; -0.008269605035289257 -0.009601008819015643 … -0.012623113779649944 -0.009601010344401282; -0.0003986984688478018 -0.0010730139336228544 … -0.004230682587213657 -0.0010730144712271835;;; 0.003153663862226204 0.0034201187641859404 … 0.004832499768515049 0.003420117892053816; 0.010276670572818683 0.010903597623974321 … 0.012272139131677476 0.010903596548685073; … ; -0.017302963197412966 -0.012623111572185437 … -0.00962026794789805 -0.01262311269028435; -0.003664696092293099 -0.004230681473674364 … -0.0054911099542692375 -0.004230682267536522;;; … ;;; -0.0023365614344967515 -0.0028154649916341096 … 0.005693309665856113 -0.0028154642394191075; -0.012913098092820667 -0.01784119372754753 … 0.011498338481776303 -0.017841192950626025; … ; -0.025174427607653295 -0.01281088428642222 … -0.008395876711159704 -0.012810885299377099; 0.007597981766331854 0.010225420522189249 … -0.005449733571608421 0.010225421536153596;;; 0.003153664724390733 0.0034201197743985554 … 0.004832501556092858 0.0034201188756595047; 0.010276671577525768 0.010903598607710235 … 0.012272140236795 0.010903597623552952; … ; -0.0173029655312882 -0.012623111843271122 … -0.009620267527183299 -0.012623112891185057; -0.0036646953598338285 -0.004230680794192003 … -0.0054911094373418165 -0.0042306816151693495;;; 0.002081000615593825 0.002223228554422875 … 0.0034201182570295892 0.002223228066288328; 0.004779270252725334 0.005734954766113007 … 0.010903596807151424 0.005734954245960775; … ; -0.00826960544046878 -0.00960100951983384 … -0.012623114159807698 -0.0096010107725226; -0.0003986982064510126 -0.001073013662487973 … -0.004230682411089262 -0.0010730142053188756;;;;], δρ0 = [-3.6021397983670894e-7 -2.488842485972921e-7 … -4.85283814595704e-8 -2.488842485542117e-7; -5.97225047831414e-7 -4.606978494356228e-7 … -1.1905772350661748e-7 -4.606978416590082e-7; … ; 1.2735794506629325e-7 1.0836344470220695e-7 … 4.640722145367663e-8 1.0836342867121045e-7; 1.0352015944575135e-7 1.1957967892533666e-7 … 5.0695090752283445e-8 1.195796758711762e-7;;; -2.488842527697492e-7 -1.731067160957187e-7 … -3.535699279706147e-8 -1.731067170622161e-7; -4.6069785103233826e-7 -3.572171013097519e-7 … -1.1091443567757726e-7 -3.572170977103749e-7; … ; 1.083634515752236e-7 1.0846503759056933e-7 … 1.0862243004385709e-7 1.0846502140585139e-7; 1.1957968094027054e-7 1.1611152445673719e-7 … 5.0208339645691304e-8 1.1611152088978999e-7;;; -4.852838265664472e-8 -3.5356990400310696e-8 … -1.1341333408701026e-8 -3.535699289238913e-8; -1.1905772162010387e-7 -1.1091442717899221e-7 … -9.770673302884241e-8 -1.109144334648286e-7; … ; 4.64072552014073e-8 1.086224512835074e-7 … 2.484525495562743e-7 1.0862244714651944e-7; 5.0695098616507124e-8 5.0208342684068194e-8 … 4.926862853355368e-8 5.0208343234683834e-8;;; … ;;; 6.434566714087435e-8 4.083859003474639e-8 … -6.637954173573408e-9 4.0838589858822386e-8; 1.6975024621990178e-7 6.683075213345401e-8 … -1.3340489187435105e-7 6.683074924462113e-8; … ; 3.271011981127181e-8 1.6023537580117955e-7 … 3.839015841129176e-7 1.6023537270007186e-7; -9.271360813524476e-8 -3.633614609291984e-8 … 7.126315486553071e-8 -3.633614729495034e-8;;; -4.852837822324172e-8 -3.535698653923727e-8 … -1.134132942284519e-8 -3.5356988147449933e-8; -1.1905771187494822e-7 -1.1091442027181233e-7 … -9.770672530045105e-8 -1.109144235069263e-7; … ; 4.640723055557599e-8 1.0862242423157704e-7 … 2.4845253052738955e-7 1.0862242266219944e-7; 5.069509246443818e-8 5.020833431471997e-8 … 4.926861957979453e-8 5.020833585512083e-8;;; -2.4888424222431505e-7 -1.7310670728500646e-7 … -3.5356988554108096e-8 -1.7310670831555274e-7; -4.6069783187609693e-7 -3.5721708536328733e-7 … -1.1091442895291588e-7 -3.572170822263454e-7; … ; 1.0836342682827283e-7 1.0846500758799862e-7 … 1.0862241278747026e-7 1.0846499891235085e-7; 1.1957967586576352e-7 1.1611151583697678e-7 … 5.020833252356362e-8 1.1611151534768416e-7;;;;], δeigenvalues = [[0.0004950066413438731, 0.09679867978374211, 0.021497254483929935, 0.0452818935344359]], δoccupation = [[0.0, 0.0, 0.0, 0.0]], δεF = 0.0, ε = DFTK.DielectricAdjoint{Array{Float64, 4}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}, Float64, Vector{Vector{Float64}}, StaticArraysCore.SVector{3, Float64}}(Hamiltonian(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), HamiltonianBlock[DFTK.DftHamiltonianBlock(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Any[DFTK.FourierMultiplication{Float64, Vector{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.0, 0.19739208802178715, 0.7895683520871486, 1.7765287921960844, 3.1582734083485944, 4.934802200544679, 7.106115168784338, 9.67221231306757, 12.633093633394378, 15.988759129764759  …  20.133992978222288, 16.38354330580833, 13.027877809437953, 10.066996489111146, 7.5008993448279115, 5.329586376588253, 3.5530575843921683, 2.1713129682396586, 1.184352528130723, 0.5921762640653614]), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.1603906580688546 0.1602471565738079 … 0.15981840036267833 0.16024715657380792; 0.16024715657380803 0.16010321540768782 … 0.15967315066597979 0.16010321540768782; … ; 0.15981840036267841 0.1596731506659797 … 0.15923916395256382 0.1596731506659797; 0.16024715657380809 0.16010321540768785 … 0.1596731506659798 0.16010321540768785;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.1599588345026577 … 0.15952744714353234 0.1599588345026577; … ; 0.15967315066597973 0.15952744714353229 … 0.15909210945052638 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772;;; 0.1598184003626783 0.15967315066597965 … 0.15923916395256382 0.15967315066597965; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256385 0.15909210945052635 … 0.15865272241895537 0.15909210945052635; 0.1596731506659798 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; … ;;; 0.15910963392556968 0.15896221057980836 … 0.1585217272926462 0.15896221057980836; 0.15896221057980842 0.15881431920357322 … 0.15837242838914242 0.15881431920357322; … ; 0.15852172729264621 0.15837242838914237 … 0.157926332522532 0.15837242838914237; 0.15896221057980844 0.15881431920357328 … 0.1583724283891425 0.15881431920357328;;; 0.1598184003626783 0.15967315066597962 … 0.15923916395256382 0.15967315066597962; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256382 0.15909210945052635 … 0.15865272241895534 0.15909210945052635; 0.15967315066597979 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.15995883450265772 … 0.15952744714353237 0.15995883450265772; … ; 0.1596731506659797 0.15952744714353229 … 0.15909210945052635 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772]), DFTK.NonlocalOperator{Float64, Matrix{ComplexF64}, Matrix{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Matrix{ComplexF64}(undef, 7809, 0), Matrix{Float64}(undef, 0, 0)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [-0.15475965282619192 -0.15461620270202853 … -0.15418753100020535 -0.15461620270102702; -0.1546162027004806 -0.1544722989002185 … -0.15404229467617217 -0.1544722988992146; … ; -0.15418753100344315 -0.15404229468098715 … -0.15360833735041635 -0.15404229467997263; -0.15461620270262993 -0.15447229890237457 … -0.1540422946783468 -0.1544722989013698;;; -0.15461620270154766 -0.15447229890128975 … -0.15404229467725372 -0.15447229890028505; -0.1544722988997383 -0.15432794309149456 … -0.153896602828854 -0.15432794309048714; … ; -0.15404229468049915 -0.1538966028336804 … -0.1534612870141286 -0.153896602832663; -0.1544722989018911 -0.15432794309365408 … -0.1538966028310331 -0.15432794309264627;;; -0.15418753100124216 -0.15404229467877958 … -0.1536083373481897 -0.15404229467776542; -0.154042294677214 -0.15389660283038506 … -0.1534612870108033 -0.1538966028293679; … ; -0.15360833735146381 -0.15346128701567252 … -0.15302191313879762 -0.1534612870146455; -0.15404229467938538 -0.15389660283256323 … -0.15346128701300116 -0.15389660283154577;;; … ;;; -0.15347880585691429 -0.15333139196727 … -0.15289092492870407 -0.15333139196624107; -0.15333139196567497 -0.15318350574406553 … -0.15274162877932296 -0.15318350574303397; … ; -0.15289092493203882 -0.15274162878428388 … -0.1522955395941201 -0.1527416287832407; -0.1533313919678913 -0.15318350574628936 … -0.1527416287815666 -0.15318350574525658;;; -0.15418753100115698 -0.15404229467869346 … -0.15360833734810111 -0.15404229467767921; -0.15404229467712385 -0.15389660283029377 … -0.15346128701070993 -0.15389660282927697; … ; -0.15360833735138338 -0.15346128701559225 … -0.15302191313871485 -0.15346128701456435; -0.15404229467930353 -0.15389660283248094 … -0.15346128701291634 -0.15389660283146292;;; -0.15461620270150722 -0.15447229890124844 … -0.15404229467720965 -0.15447229890024336; -0.15447229889969444 -0.1543279430914498 … -0.1538966028288073 -0.15432794309044243; … ; -0.15404229468045977 -0.15389660283364112 … -0.1534612870140872 -0.15389660283262285; -0.15447229890185227 -0.15432794309361486 … -0.1538966028309909 -0.15432794309260625]), 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.015163897242399138 -0.012550128941215296 … -0.0049668692271259245 -0.0125501289089789; -0.01255012888713324 -0.01036259507361962 … -0.004135159924350544 -0.010362594967977772; … ; -0.004966868670028893 -0.004135159815189204 … -0.0020697815215944856 -0.004135159428525858; -0.012550128785006864 -0.01036259499260944 … -0.004135159925762504 -0.010362594889026121;;; -0.012550129008192346 -0.010362595173618636 … -0.0041351602919268305 -0.010362595158264085; -0.01036259505210508 -0.008622039162173243 … -0.003897512347271458 -0.008622039102826157; … ; -0.004135159980116132 -0.0038975126271501752 … -0.003484108080597781 -0.003897512261034482; -0.010362595119231034 -0.008622039308789603 … -0.003897512618009398 -0.00862203920039679;;; -0.0049668692587544715 -0.004135160116514563 … -0.0020697826141986026 -0.004135160274404361; -0.004135159848564246 -0.0038975121453608576 … -0.003484108047203313 -0.003897512280737399; … ; -0.0020697824915633855 -0.0034841085152677027 … -0.005620580894745398 -0.00348410843254671; -0.0041351603329000675 -0.003897512755869281 … -0.0034841086500084642 -0.0038975127986950703;;; … ;;; -0.006547895260209623 -0.004876337307235974 … -0.0014866478959362976 -0.00487633729898515; -0.0048763376460476265 -0.002677985874930617 … -0.004123055715730382 -0.0026779857988987604; … ; -0.001486648220294076 -0.004123056203614081 … -0.007170896073380846 -0.0041230561436311495; -0.004876337166471078 -0.0026779851625507 … -0.004123056103602125 -0.0026779852114647996;;; -0.004966869124306767 -0.004135159941583723 … -0.002069782225538419 -0.004135160062093822; -0.0041351597052014105 -0.003897512052861042 … -0.0034841079259495212 -0.0038975121235282567; … ; -0.002069781783400368 -0.003484107961958493 … -0.005620580605749714 -0.003484107924826894; -0.004135159976643154 -0.0038975123225235137 … -0.0034841082071842317 -0.0038975123929460426;;; -0.012550128818813506 -0.01036259493731051 … -0.004135160042999094 -0.010362594968946233; -0.010362594882297339 -0.008622038992112885 … -0.0038975122274889167 -0.008622038947400695; … ; -0.004135159382210181 -0.003897511946321002 … -0.0034841077233684053 -0.003897511745243486; -0.010362594812140836 -0.008622038914628562 … -0.0038975122467766153 -0.008622038901003558])], 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.009532891999736448 -0.006919175069435931 … 0.0006640001353470544 -0.0069191750361980026; -0.0069191750138058086 -0.004731678566150289 … 0.0014956960654570707 -0.00473167845950455; … ; 0.0006640006892063705 0.0014956961698033462 … 0.0035610450805529906 0.0014956965574812135; -0.0069191749138287125 -0.004731678487296163 … 0.0014956960618705169 -0.0047316783827080645;;; -0.006919175135932144 -0.004731678667220726 … 0.0014956956967990941 -0.004731678650861451; -0.00473167854415555 -0.002991147751010113 … 0.0017333319674068771 -0.0029911476906556104; … ; 0.0014956960053644445 0.0017333316827017149 … 0.0021467143558000054 0.0017333320498348166; -0.004731678613434282 -0.0029911478997859675 … 0.0017333316944899298 -0.0029911477903853487;;; 0.0006640001026816701 0.0014956958706854985 … 0.003561043990175509 0.001495695713809862; 0.0014956961402014801 0.001733332167786452 … 0.002146714392519786 0.001733332033427069; … ; 0.0035610441095366508 0.0021467139195861244 … 1.0228385412347149e-5 0.002146714003334129; 0.0014956956536943677 0.0017333315550998544 … 0.002146713787516837 0.001733331513291529;;; … ;;; -0.0009170671915542314 0.0007544813053023888 … 0.004144154468005819 0.0007544813145821403; 0.0007544809680858217 0.002952827584577077 … 0.001507743894089083 0.0029528276616404975; … ; 0.004144154140313319 0.0015077434012444072 … -0.0015401031449689476 0.0015077434622705046; 0.0007544814454460598 0.0029528282947332177 … 0.001507743503973773 0.002952828246851903;;; 0.0006640002372145567 0.001495696045702437 … 0.0035610443789242884 0.001495695926206582; 0.0014956962836544657 0.0017333322603775556 … 0.0021467145138669473 0.0017333321907271387; … ; 0.003561044817780076 0.002146714472975603 … 1.0228674490769807e-5 0.002146714511135102; 0.001495696010033105 0.001733331988527917 … 0.0021467142304258907 0.001733331919123407;;; -0.006919174946512863 -0.0047316784308713 … 0.001495695945770907 -0.004731678461501911; -0.0047316783743039555 -0.0029911475809049574 … 0.0017333320872361586 -0.002991147535185406; … ; 0.0014956966033097523 0.0017333323635701624 … 0.002146714713070737 0.001733332565665947; -0.004731678306305255 -0.0029911475055857083 … 0.001733332065764901 -0.0029911474909520945]), 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.3598129937764474e-6 - 1.6783057301187582e-6im 0.008534855822536457 + 0.0019279215062177905im … 0.005732626680871501 + 0.0012424393045692961im -0.008436827017987614 - 0.0018449069892143885im; -0.011486808192503436 + 0.015516839130785727im 0.00530506832986178 - 0.015073287566117605im … -0.01015047371254714 + 0.007575687774973773im 0.015249644957531967 - 0.012862153322545063im; … ; -0.013615815860772566 + 0.018382866512791574im 0.01358597172449428 - 0.015646795379586666im … -0.0057409590295962164 + 0.01072123205978074im 0.010172322745573748 - 0.016406022116030906im; 0.011419223979299818 - 0.015545111237841636im -0.015271968027527394 + 0.012830504554969513im … 0.002279424288588841 - 0.009280091640553996im -0.005327314192806442 + 0.015041605725053246im;;; -0.0012700191601404073 + 0.0018177467934392027im -0.004258717419995159 - 0.002197696841800463im … -0.004187857236068072 - 0.00048326503771240305im 0.005685860679502422 + 1.3401624626523454e-5im; 0.011046995142152724 - 0.015013280063657378im -0.006035113697101563 + 0.013306818236839846im … 0.008374218176419247 - 0.00676235634006892im -0.012437710804785949 + 0.011883036472964723im; … ; 0.011616916642186968 - 0.015664377917743534im -0.011354640730418475 + 0.013196086490059218im … 0.00495458035939773 - 0.009165564844771205im -0.008674992426467831 + 0.01379206252975394im; -0.009529860093558322 + 0.012890493688985132im 0.011505764812676077 - 0.01047989768794601im … -0.0021072435647668235 + 0.007451255971060225im 0.0051031176073598776 - 0.011903660836829992im;;; 0.0008457486139847409 - 0.0012542858676570698im 0.0011392760690318057 + 0.0012527575678119096im … 0.002163722680285419 - 1.0039805712947417e-5im -0.0022743415485932217 + 0.0004935083473227086im; -0.006792258569709043 + 0.009294154419520412im 0.00432098959818638 - 0.008060513316225384im … -0.0051597226920626175 + 0.0045535920433002995im 0.007000597603605749 - 0.007464538370068565im; … ; -0.007356272784984583 + 0.009892649824807212im 0.007196727598508396 - 0.008463371636567902im … -0.003435645085448897 + 0.0062343081471368435im 0.005656506097212783 - 0.008805932067436722im; 0.005637559613486486 - 0.007561641180787509im -0.006160462125187935 + 0.0061530924631741195im … 0.0015395508905390264 - 0.00453112368799009im -0.0034808540500733337 + 0.006749110033619394im;;; … ;;; 0.00042345652839788627 - 0.000524169801763404im -0.001033024760895541 + 0.0002894476882670745im … -0.0006810668889896542 - 0.0004910565368266439im 0.0003085639655363751 + 0.0005878313168979001im; 0.002570062822632905 - 0.003535939743519392im -0.001697448441776149 + 0.003273464410863721im … 0.0024803927188792483 - 0.002101752513012018im -0.002861687386346646 + 0.003014512958458684im; … ; 0.004380376421394324 - 0.005887972725298491im -0.0043635557843199075 + 0.005218576559424717im … 0.0023772005604464266 - 0.004157624142647901im -0.0035442213229096827 + 0.005400676376425878im; -0.003292713991966947 + 0.00441426774125845im 0.003492167143638539 - 0.0037638959292913927im … -0.0012635282212573062 + 0.002975285193517128im 0.002327869383782804 - 0.0040228083984164605im;;; -0.0009215046878975858 + 0.0012199923327498967im 0.0023526386734288793 - 0.000446036792796765im … 0.001580848697677935 + 0.0008059865047900963im -0.001061036878001842 - 0.0012052534470403757im; -0.005579016684201482 + 0.007595408451883534im 0.0034154662034498974 - 0.006792636954693037im … -0.00467117597594946 + 0.0038696134213893362im 0.006095080205959 - 0.006196658878329666im; … ; -0.007939226368652103 + 0.010708682172906367im 0.007685283114780173 - 0.009147377927847075im … -0.003754414786287238 + 0.0066806078664788816im 0.006145065906138036 - 0.009489842809682254im; 0.006850985750356997 - 0.009260402887971106im -0.00706613744829472 + 0.007421011341531532im … 0.0020280959864664595 - 0.005215015391214458im -0.004386480203211227 + 0.008016964953684282im;;; 0.0013217784873133763 - 0.0018109234847926443im -0.005777108178071617 - 7.177198291673683e-5im … -0.003666827814390487 - 0.0012128372520098032im 0.004167529837996324 + 0.0021393230624910074im; 0.00952875539340968 - 0.012887385117484669im -0.0050577509295637 + 0.011938235751023115im … 0.007965244815165726 - 0.006189651528090546im -0.011460371715284766 + 0.010514458119232534im; … ; 0.012138080770542355 - 0.016393921595856023im -0.011763717633206574 + 0.013768795800566558im … 0.0051897068488593375 - 0.009494745598509186im -0.009084097040988528 + 0.01436474515705486im; -0.011048372631981038 + 0.015016410459622872im 0.012483307310416221 - 0.011848490259843084im … -0.0025162954448350204 + 0.008023952969590986im 0.006080629757985978 - 0.013272251732688728im],)])]), [8.519872575203825e-7 4.627264614697478e-7 … 2.4177794405552907e-8 4.627264576437038e-7; 4.627264550508849e-7 2.5014474726854846e-7 … 1.3583004074256083e-8 2.501447390978284e-7; … ; 2.417778585754355e-8 1.3583002948084073e-8 … 1.5659245195189177e-9 1.358299895903291e-8; 4.6272644292979953e-7 2.5014474100293103e-7 … 1.3583004088824331e-8 2.5014473299144264e-7;;; 4.627264694190749e-7 2.501447550028264e-7 … 1.3583007866397426e-8 2.501447538152544e-7; 2.5014474560454366e-7 1.3895607990281521e-7 … 1.1279995980210638e-8 1.389560768515441e-7; … ; 1.3583004649565779e-8 1.1279998521401413e-8 … 7.937761037063379e-9 1.1279995197213453e-8; 2.501447507962996e-7 1.389560874409742e-7 … 1.1279998438404535e-8 1.389560818680488e-7;;; 2.4177794890856096e-8 1.3583006056736706e-8 … 1.565927085210206e-9 1.3583007685621539e-8; 1.3583003292399788e-8 1.1279994146944348e-8 … 7.937760798869773e-9 1.1279995376109854e-8; … ; 1.5659267972350229e-9 7.937764137426403e-9 … 3.571857559314586e-8 7.93776354740447e-9; 1.3583008289102075e-8 1.1279999690120701e-8 … 7.937765098492583e-9 1.1280000078959785e-8;;; … ;;; 5.791834964078162e-8 2.2815986872554797e-8 … 5.620413558683402e-10 2.281598675091543e-8; 2.281599186763881e-8 3.4892674385969185e-9 … 1.345852407826381e-8 3.4892671298622237e-9; … ; 5.620417346989821e-10 1.3458529079799118e-8 … 7.728317657617618e-8 1.3458528464880515e-8; 2.2815984797266104e-8 3.4892645459047014e-9 … 1.3458528054525309e-8 3.489264744525487e-9;;; 2.4177792827910927e-8 1.3583004252047145e-8 … 1.5659261725440925e-9 1.358300549529913e-8; 1.3583001813383153e-8 1.1279993307084058e-8 … 7.937759934005646e-9 1.1279993948714192e-8; … ; 1.565925134299415e-9 7.937760190845163e-9 … 3.571856978937249e-8 7.937759925998686e-9; 1.3583004613740345e-8 1.1279995755510133e-8 … 7.93776193996498e-9 1.1279996394917529e-8;;; 4.627264469422035e-7 2.501447367259263e-7 … 1.3583005298304533e-8 2.50144739172743e-7; 2.501447324710068e-7 1.3895607115932016e-7 … 1.1279994892634179e-8 1.3895606886047962e-7; … ; 1.3582998481211735e-8 1.1279992339740208e-8 … 7.937758489059375e-9 1.1279990514038214e-8; 2.50144727044856e-7 1.3895606717552619e-7 … 1.1279995067756631e-8 1.389560664750175e-7;;;;], Matrix{ComplexF64}[[-0.20966975200961321 - 0.15461424809128965im -0.7273395035438232 + 0.6343160586743704im 5.215611218571905e-7 + 1.6122531184680274e-6im -2.270496430166571e-7 + 8.263144845147691e-7im; 0.14225650311209007 + 0.10490250525299848im -0.032049464386048616 + 0.027950496576054903im 0.4968615924718062 - 0.1327332324524069im -0.08727779861424632 + 0.25961615029075175im; … ; -0.04519595414740601 - 0.03332830982073979im 0.011591861930863662 - 0.010109313965719917im 0.01090470011438079 - 0.002796931212490305im -0.0053481053490199645 - 7.98975210214783e-5im; 0.08103665510832071 + 0.059757887662263184im -0.022964387937257563 + 0.020027347499178477im -0.03386081610785774 + 0.00844587925720101im 0.02502410212006623 + 0.014460704700362597im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944124183226, [[-0.5541104208624663, -0.010678403974178807, 0.14495661104247287, 0.14495827249111543]], 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.6206936080073187e-7, -6.293616355120162e-7, -2.2951215222556635e-7, 4.994023869209871e-7, 7.923307328904864e-7, 5.851065657394511e-7, 6.516966006032196e-7, 1.0189306847011097e-6, 3.6385909546541705e-7, -2.201093620405244e-6  …  2.7810188812433237e-6, 1.2016958552075658e-6, 2.445648949409408e-8, 1.007857108035993e-7, 5.535313616218756e-7, 4.534179065310575e-7, 4.632288696897457e-9, -1.0437042231472339e-7, 1.121110024604581e-7, 1.255086995969166e-7], resid_history = [0.24939209547385532, 0.003766548284935027, 0.0002852379782647066, 4.681648283525091e-6, 3.076425747029639e-8, 8.454870414863744e-11, 4.225020447012707e-8, 1.1686769767536589e-9], converged = true, n_iter = 8, residual_norm = 1.1686769767536589e-9, maxiter = 100, tol = 1.0e-8, s = 0.9354772332985191, residual = [9.316691593564759e-7, 6.147142996580742e-7, 1.2746256826952922e-7, -1.8403556488856371e-7, -2.1351571572619763e-7, -1.2201083292487228e-7, -1.090126416499476e-7, -1.4364761424453332e-7, -4.5993778939095215e-8, 2.563443347584118e-7  …  9.236602845970793e-7, 4.0934394409634677e-7, 8.897294535987042e-9, 3.9676288228197694e-8, 2.3363979995529634e-7, 2.076757574003618e-7, 2.456836754281465e-9, -7.230222634252675e-8, 1.242305652802141e-7, 4.100099165788321e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.3523792902611464e-7, -4.113177928674583e-8, 1.1686769767536589e-9, -0.0002044345494043764, 4.686294570539404e-6, -3.112003400211183e-8, 8.454870414863744e-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}}([[-9.816303754180385e-7, -2.4375955759694647e-6, -8.026051547649222e-7, 1.3371467400401934e-6, 1.6827451615705585e-6, 1.1918317120981137e-6, 1.4039501691751438e-6, 2.1037055390177912e-6, 6.122858404987716e-7, -2.993981358663981e-6  …  3.2426628486275657e-6, 1.2791749429970903e-6, 2.9657677446297494e-8, 1.4185073366306118e-7, 7.560392939110402e-7, 5.02899693159863e-7, 4.273120201919365e-9, -1.0664961999661448e-7, 1.5178994534074338e-7, 1.54859945072673e-7], [1.5896124863079968e-6, 3.0716549891845404e-6, 9.883605584261743e-7, -1.660303954213417e-6, -2.0985513168440015e-6, -1.4639221886847015e-6, -1.6969276718095186e-6, -2.5415539754968225e-6, -7.494098165139761e-7, 3.6922853652692726e-6  …  -3.132436881603706e-6, -1.2018795220857281e-6, -2.7833788485730605e-8, -1.3452457195400373e-7, -7.083561140452157e-7, -4.420173527659818e-7, -3.285597911570402e-9, 7.268937629138518e-8, -8.879772952869992e-8, 8.500490665999483e-8]]), H = [1.1377088405310873 0.01382823232717191 … 0.0 0.0; 0.12700150952699651 1.0343242308996667 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1447754317964405 0.12849091336745613 … 0.0 0.0; 0.0 1.0267982276637933 … 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.9257125333940879
Interacting polarizability:     1.7736548639098646

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

Manual solution of the Dyson equations

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

using KrylovKit

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

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

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

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

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

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