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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457365182548757
Displaced dipole:  0.01761222136715441
Polarizability :   1.7746795018979897

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.0013432920150998988 - 0.0030097321478734377im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.2835823733712559 + 0.12372240472565905im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.037566882339560424 + 0.016585421327890263im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.0737839293216131 - 0.03327046265750991im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206945992167514e-7 -2.502015242841601e-7 … -4.8806438746072914e-8 -2.502015261612165e-7; -6.293617172387298e-7 -4.793036722973558e-7 … -1.20778453903354e-7 -4.7930367527009e-7; … ; 1.3626262490396954e-7 1.1356231137893779e-7 … 4.7119918788512005e-8 1.1356231414928808e-7; 1.3305477132684266e-7 1.3632641090458082e-7 … 5.200686169427702e-8 1.3632641136246134e-7;;; -2.502015468183431e-7 -1.7404910379300736e-7 … -3.556495292393309e-8 -1.7404910537205156e-7; -4.793036901789528e-7 -3.679732348270739e-7 … -1.1202896482849184e-7 -3.679732373956012e-7; … ; 1.1356230838136365e-7 1.1211099191191796e-7 … 1.0980684740619285e-7 1.1211099402461453e-7; 1.3632639516949524e-7 1.2550865479068416e-7 … 5.094909006715297e-8 1.2550865500971087e-7;;; -4.8806447683886647e-8 -3.556495556434352e-8 … -1.1412156063623115e-8 -3.556495627148116e-8; -1.207784606548701e-7 -1.1202896588644016e-7 … -9.786226999719196e-8 -1.1202896792492928e-7; … ; 4.7119912695137275e-8 1.0980684071609404e-7 … 2.493823111190104e-7 1.098068435621518e-7; 5.200685323362737e-8 5.0949085607839736e-8 … 4.9106278636633006e-8 5.0949086280079586e-8;;; … ;;; 6.475758950006499e-8 4.110332249751959e-8 … -6.681699393677697e-9 4.1103321784383334e-8; 1.6897555717574396e-7 6.653284045394773e-8 … -1.3294451521006488e-7 6.653283845209282e-8; … ; 3.267766510544614e-8 1.6017371660876286e-7 … 3.847037114751066e-7 1.6017372306190105e-7; -9.141992911595918e-8 -3.5830344633607625e-8 … 7.037286966111584e-8 -3.58303422210287e-8;;; -4.8806428841110405e-8 -3.556494142950822e-8 … -1.1412151965062323e-8 -3.556494260492935e-8; -1.207784466615189e-7 -1.1202895305697118e-7 … -9.786226011132327e-8 -1.1202895580077346e-7; … ; 4.7119920163145456e-8 1.0980685227901358e-7 … 2.4938233841960184e-7 1.09806859015424e-7; 5.200686819547855e-8 5.0949098456830855e-8 … 4.9106291206519206e-8 5.094910081453319e-8;;; -2.502015042390406e-7 -1.740490727992395e-7 … -3.556494613747939e-8 -1.7404907439522513e-7; -4.793036586072519e-7 -3.6797321109208596e-7 … -1.1202895818955006e-7 -3.679732141127413e-7; … ; 1.1356231687509527e-7 1.1211099840211504e-7 … 1.098068557171799e-7 1.1211100315703381e-7; 1.363264266967605e-7 1.2550867759226726e-7 … 5.0949097598458056e-8 1.255086789720135e-7;;;;], δHtotψ = Matrix{ComplexF64}[[0.0007813759335248623 + 0.001750723029789715im -0.06565028578280961 - 0.06907653723964724im 0.18150054406779023 - 1.9105640233770935im 0.10211926350120865 - 0.58036345736914im; 0.1818115219416958 - 0.07917170947268846im -1.2426975888518152 + 1.049941887222828im 0.039908306353043044 - 0.42213044826034263im 0.022407967994391598 - 0.128269896157212im; … ; 0.04489753431235969 - 0.01967649786213845im 0.014958809399169635 - 0.01428852982392358im -0.007596757681120127 + 0.026020751372040293im 0.001632736235256345 + 0.01606009122380993im; -0.06597916458596788 + 0.029745866394809365im -0.018589223792885398 + 0.017588072177326868im -0.04474212869477096 - 0.04977999524656663im 0.010538070433858883 + 0.08986970005593456im]], δVind = [0.0019654321204911245 0.0020810008414760494 … 0.003153665059087006 0.0020810008494415795; 0.003889131624004337 0.004779270223505149 … 0.010276671677136329 0.004779270236208289; … ; -0.0061508482408141356 -0.008269603579872546 … -0.01730296653039434 -0.008269603487425784; 0.00023868205354685753 -0.00039869758451732554 … -0.003664694851067146 -0.0003986975788218062;;; 0.002081001027664919 0.0022232287464950496 … 0.003420120078691899 0.002223228759953414; 0.004779270427719571 0.0057349546479498106 … 0.010903598491266035 0.005734954666538541; … ; -0.008269603205420448 -0.009601007878083472 … -0.012623111182792031 -0.009601007800240395; -0.00039869737681847285 -0.001073012995533217 … -0.004230680327118428 -0.0010730129869507441;;; 0.0031536656881303313 0.003420120415042535 … 0.004832505528460608 0.003420120430662942; 0.010276672343754147 0.010903598855339428 … 0.012272142376872852 0.010903598838087114; … ; -0.017302966084678652 -0.012623110990596358 … -0.009620265236303119 -0.012623110861922268; -0.0036646941950820005 -0.004230680024259988 … -0.005491107524549023 -0.004230679987196038;;; … ;;; -0.002336560009249448 -0.0028154627424314444 … 0.005693314635886586 -0.002815462768777948; -0.012913093900772474 -0.01784118355102107 … 0.011498339288968889 -0.01784118376758667; … ; -0.025174457888908368 -0.012810887580253122 … -0.008395876319721487 -0.012810887356514531; 0.007597980434740037 0.010225416282347586 … -0.005449734305816252 0.010225416504910557;;; 0.0031536644397257975 0.0034201190539621296 … 0.004832503376738135 0.00342011902034988; 0.010276670981836222 0.010903597480256303 … 0.012272140832167681 0.010903597393589325; … ; -0.017302967858336907 -0.012623112520869928 … -0.009620266313661651 -0.012623112185293637; -0.0036646955742098528 -0.004230681466049006 … -0.005491108961445142 -0.004230681374709451;;; 0.0020810006682639125 0.0022232283506439412 … 0.003420119367089527 0.0022232283523158673; 0.004779270033322158 0.005734954219240667 … 0.01090359778230182 0.005734954220182985; … ; -0.008269603871909621 -0.00960100861372976 … -0.012623111809029435 -0.009601008429061374; -0.0003986977834699514 -0.0010730134474601144 … -0.004230681009048244 -0.0010730134380491447;;;;], δρ0 = [-3.602140499929466e-7 -2.488842925267328e-7 … -4.852837889543878e-8 -2.4888429441288775e-7; -5.97225065014871e-7 -4.606978519032999e-7 … -1.1905770175903855e-7 -4.6069785492679586e-7; … ; 1.2735791128564377e-7 1.0836340716963843e-7 … 4.6407202065071643e-8 1.0836340987676092e-7; 1.0352007508151378e-7 1.1957961953542333e-7 … 5.069507064786281e-8 1.1957962004814767e-7;;; -2.4888429412431463e-7 -1.7310674133137036e-7 … -3.5356986606721235e-8 -1.731067427425829e-7; -4.6069785452649206e-7 -3.572170975351168e-7 … -1.1091441425924387e-7 -3.572171000132808e-7; … ; 1.0836340764373949e-7 1.0846498456336554e-7 … 1.0862238930107033e-7 1.0846498666701457e-7; 1.195796190643063e-7 1.1611147649641079e-7 … 5.020831474493602e-8 1.1611147685800875e-7;;; -4.852837838973657e-8 -3.535698575579334e-8 … -1.1341321758139981e-8 -3.535698641253758e-8; -1.190577014928407e-7 -1.1091441217750037e-7 … -9.77067084338314e-8 -1.1091441417131531e-7; … ; 4.640719835819132e-8 1.0862238525270766e-7 … 2.4845247853071614e-7 1.0862238810348001e-7; 5.0695069227014005e-8 5.020831346252808e-8 … 4.9268595418591284e-8 5.0208314172100064e-8;;; … ;;; 6.434570755452527e-8 4.0838623568356714e-8 … -6.6379372792493665e-9 4.0838622930004845e-8; 1.697503413037899e-7 6.683083192095439e-8 … -1.3340484518800983e-7 6.683082993075683e-8; … ; 3.2710014663562915e-8 1.6023526670891823e-7 … 3.8390149047418243e-7 1.6023527321022764e-7; -9.271365986072388e-8 -3.6336199242155724e-8 … 7.126310801654966e-8 -3.633619679002058e-8;;; -4.852837830711667e-8 -3.535698554681433e-8 … -1.1341322363762558e-8 -3.535698677334338e-8; -1.1905770151497815e-7 -1.1091441201448026e-7 … -9.770670913930686e-8 -1.109144147883053e-7; … ; 4.640720111173347e-8 1.0862238631209298e-7 … 2.484524836677423e-7 1.0862239300820619e-7; 5.069507023339197e-8 5.0208313701224615e-8 … 4.926859747108766e-8 5.0208316014040124e-8;;; -2.4888429323423224e-7 -1.7310674056167967e-7 … -3.535698682601622e-8 -1.731067423359695e-7; -4.606978534755119e-7 -3.5721709627332273e-7 … -1.1091441400102002e-7 -3.5721709942242627e-7; … ; 1.0836340915469452e-7 1.0846498477815928e-7 … 1.0862239234896612e-7 1.0846498942438078e-7; 1.1957962028457145e-7 1.1611147698908084e-7 … 5.020831594422291e-8 1.1611147827205062e-7;;;;], δeigenvalues = [[0.0004950050373552736, 0.09679862707521988, 0.010545849030840932, 0.04819124015917561]], δ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.15475965281399257 -0.15461620268929516 … -0.1541875309891694 -0.15461620268940018; -0.1546162026893418 -0.1544722988885458 … -0.15404229466620556 -0.15447229888864947; … ; -0.1541875309890598 -0.1540422946660484 … -0.1536083373371759 -0.15404229466615874; -0.15461620268934795 -0.15447229888855024 … -0.1540422946662158 -0.15447229888865743;;; -0.15461620268959977 -0.15447229888880426 … -0.15404229466646752 -0.15447229888890854; -0.15447229888885022 -0.15432794308007 … -0.15389660281913917 -0.15432794308017292; … ; -0.15404229466635924 -0.15389660281898404 … -0.15346128700113432 -0.15389660281909348; -0.154472298888857 -0.15432794308007508 … -0.15389660281914985 -0.15432794308018152;;; -0.15418753098955887 -0.15404229466655212 … -0.1536083373376747 -0.15404229466665623; -0.15404229466659763 -0.1538966028192265 … -0.15346128700136866 -0.15389660281932918; … ; -0.15360833733756726 -0.15346128700121484 … -0.15302191312605587 -0.1534612870013239; -0.1540422946666051 -0.15389660281923226 … -0.15346128700137993 -0.15389660281933848;;; … ;;; -0.1534788058441136 -0.15333139195391415 … -0.15289092491708814 -0.15333139195402548; -0.15333139195396367 -0.15318350573180028 … -0.15274162876880715 -0.15318350573191014; … ; -0.15289092491697148 -0.1527416287686399 … -0.15229553958021533 -0.15274162876875658; -0.15333139195396986 -0.1531835057318046 … -0.1527416287688174 -0.15318350573191822;;; -0.15418753098854296 -0.15404229466553088 … -0.153608337336651 -0.15404229466563948; -0.15404229466557923 -0.15389660281820272 … -0.15346128700034237 -0.15389660281830994; … ; -0.15360833733653723 -0.15346128700017916 … -0.15302191312501787 -0.153461287000293; -0.15404229466558525 -0.15389660281820697 … -0.1534612870003523 -0.1538966028183178;;; -0.15461620268908868 -0.15447229888829056 … -0.15404229466595246 -0.15447229888839692; -0.15447229888833794 -0.1543279430795551 … -0.15389660281862283 -0.1543279430796602; … ; -0.15404229466584096 -0.153896602818463 … -0.15346128700061204 -0.15389660281857478; -0.15447229888834385 -0.1543279430795593 … -0.1538966028186328 -0.15432794307966788]), 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.01516389797555193 -0.01255012945842314 … -0.004966868777475462 -0.012550129478969111; -0.012550129452251691 -0.010362595425268442 … -0.004135159585338706 -0.010362595444497209; … ; -0.004966868761875513 -0.004135159582159268 … -0.002069781233416358 -0.004135159644959849; -0.012550129478788169 -0.010362595452787855 … -0.004135159655844733 -0.010362595483705342;;; -0.012550129469913145 -0.010362595448414991 … -0.004135159614496233 -0.01036259546257879; -0.010362595446004923 -0.0086220394071127 … -0.003897512011838582 -0.008622039424126047; … ; -0.004135159596629266 -0.003897511997580895 … -0.0034841076719392414 -0.0038975120430575637; -0.010362595462235856 -0.008622039427879325 … -0.00389751205403153 -0.008622039448530665;;; -0.004966868745112358 -0.004135159570310172 … -0.0020697811543622497 -0.004135159599644326; -0.004135159583592221 -0.003897511975693534 … -0.0034841076268981217 -0.003897512011098698; … ; -0.0020697811287179327 -0.003484107593305543 … -0.005620580495139003 -0.0034841076488909403; -0.0041351595929303864 -0.003897511993101137 … -0.003484107653224158 -0.003897512028924464;;; … ;;; -0.006547897511241156 -0.004876339576905196 … -0.0014866453384181427 -0.004876339516801416; -0.004876339547488885 -0.0026779880140040836 … -0.004123054900895361 -0.0026779879596435163; … ; -0.001486645294554737 -0.004123054852442956 … -0.007170895724148087 -0.00412305495670786; -0.004876339543888884 -0.0026779880325768744 … -0.004123054942305688 -0.0026779879163594776;;; -0.004966868757717146 -0.004135159563719204 … -0.0020697812354554725 -0.004135159645599936; -0.004135159584037149 -0.0038975119697258077 … -0.0034841076453160673 -0.003897512026681945; … ; -0.002069781197335965 -0.0034841076031852554 … -0.005620580554299915 -0.0034841077380067245; -0.004135159621044321 -0.003897511985450482 … -0.0034841077305501635 -0.0038975120985645477;;; -0.012550129467859241 -0.010362595435381688 … -0.004135159643557364 -0.010362595474350807; -0.010362595434586865 -0.008622039388416557 … -0.003897512008371058 -0.008622039419819618; … ; -0.004135159624163085 -0.003897511994749357 … -0.003484107728862443 -0.0038975120984710856; -0.010362595468633756 -0.008622039416735328 … -0.0038975121020418608 -0.008622039474742])], 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.009532892720689887 -0.006919175573910405 … 0.000664000596033467 -0.0069191755945613755; -0.006919175567785447 -0.004731678906126422 … 0.00149569641443552 -0.004731678925458855; … ; 0.0006640006117431055 0.0014956964177720265 … 0.0035610453819715547 0.0014956963548611174; -0.006919175594328031 -0.004731678933650247 … 0.0014956963439192788 -0.004731678964674926;;; -0.006919175585705055 -0.004731678929531597 … 0.0014956963850158971 -0.004731678943799646; -0.004731678927167325 -0.002991147984525013 … 0.001733332312554593 -0.002991148001641278; … ; 0.0014956964029912215 0.0017333323269673505 … 0.0021467147774528175 0.0017333322813812414; -0.004731678943405003 -0.002991148005296691 … 0.0017333322703510425 -0.0029911480260544736;;; 0.0006640006280070766 0.0014956964291173592 … 0.0035610454605268676 0.001495696399679094; 0.0014956964157898742 0.001733332348612322 … 0.0021467148222596246 0.00173333231310449; … ; 0.003561045486278654 0.0021467148560059693 … 1.0228797760494096e-5 0.0021467148003115203; 0.0014956964064443256 0.0017333323311989734 … 0.002146714795922375 0.0017333322952694259;;; … ;;; -0.0009170694297850873 0.0007544790489890114 … 0.004144157037139905 0.0007544791089814634; 0.0007544790783558615 0.002952825457768855 … 0.001507744719439914 0.0029528255120195657; … ; 0.004144157081119995 0.0015077447680595189 … -0.0015401027818314223 0.0015077446636779304; 0.0007544790819497011 0.0029528254391918175 … 0.001507744678019429 0.002952825555295583;;; 0.0006640006164181983 0.0014956964367295374 … 0.0035610453804573537 0.0014956963547401985; 0.0014956964163633539 0.0017333323556038404 … 0.0021467148048679692 0.0017333322985404833; … ; 0.0035610454186906313 0.0021467148471619284 … 1.0228739637557077e-5 0.0021467147122266337; 0.001495696379350214 0.0017333323398749196 … 0.002146714719623992 0.0017333322266500258;;; -0.00691917558314006 -0.004731678915984593 … 0.0014956963564698262 -0.0047316789550600445; -0.004731678915236983 -0.00299114796531395 … 0.0017333323165384818 -0.002991147996822093; … ; 0.0014956963759756546 0.0017333323303199161 … 0.002146714721051865 0.0017333322264864157; -0.0047316789492897575 -0.002991147993636912 … 0.0017333322228577702 -0.002991148051752164]), DFTK.NonlocalOperator{Float64, Matrix{ComplexF64}, Matrix{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Matrix{ComplexF64}(undef, 7809, 0), Matrix{Float64}(undef, 0, 0)), nothing, @NamedTuple{ψ_reals::Array{ComplexF64, 3}}[(ψ_reals = [1.91434524965993e-6 + 3.504914918775587e-6im 9.509086019608679e-5 + 0.00027165827618728515im … 9.367290530305409e-5 + 0.0002645507627690392im -0.00012323675973045475 - 0.00035133569579273575im; 0.015099071790614028 - 0.006701166916232437im -0.014400885614008738 + 0.006227180800300447im … 0.008984544492562542 - 0.004225122796054904im -0.014276891598105863 + 0.0065925476584534235im; … ; 0.019190342260505124 - 0.0085095193914752im -0.01703167363417511 + 0.00762008796610269im … 0.01118738480403887 - 0.0048599126836937475im -0.01707556663042074 + 0.0074926504405564524im; -0.015078991724645493 + 0.006767927675283541im 0.014283141092746046 - 0.006577232085425589im … -0.009118116462112497 + 0.0038527134879732415im 0.014411716291881805 - 0.006209252024723135im;;; 0.000603519131303308 + 0.006756904863596504im -0.0004053832815939182 - 0.004117646727152369im … 5.450584759609781e-5 + 0.00117290014857356im -0.0002809605763952817 - 0.00375102738825315im; -0.01470074681540355 + 0.002425722950225447im 0.012796063426148098 - 0.0029376912961204013im … -0.0077049492354895 + 0.002495752313537568im 0.012718226633420565 - 0.003173944738512586im; … ; -0.016930004974957627 + 0.00892533873722412im 0.014687188327140315 - 0.0076773622340524095im … -0.009717095269480695 + 0.0048611098328173175im 0.01472142969034672 - 0.00757726544282354im; 0.01398632675980166 - 0.010377648621169532im -0.012274526680987947 + 0.008253042829714567im … 0.007619832118420282 - 0.0043427650522926825im -0.012355430933841983 + 0.008015257908380672im;;; -0.00039796829111153555 - 0.004584964292777355im 0.00028729597819171314 + 0.003199369094671698im … -9.56061180542694e-5 - 0.0014164341798486673im 0.0002464681634798419 + 0.003072383603532329im; 0.009320865864155509 - 0.0008940156929962067im -0.007873913344992732 + 0.00103330835636489im … 0.0050006132370672834 - 0.0010141860249473756im -0.007842018223657281 + 0.0011326309081461552im; … ; 0.0110250237097648 - 0.00646228601830471im -0.009627323393818483 + 0.005615377210883847im … 0.006732262582575259 - 0.0037852083106139174im -0.00964685514119985 + 0.005557604336328336im; -0.008780533806561703 + 0.007184950203542392im 0.007449187105001765 - 0.005806431151800842im … -0.004845587632712676 + 0.0033796268521964806im 0.007482075729148275 - 0.005706030654467679im;;; … ;;; -0.0002082774840181169 - 0.0020930181476815648im 0.00016980167506820402 + 0.0017471403342445722im … -0.00013453689722947365 - 0.0012215041173107358im 0.00018652101822162626 + 0.0017974842913756584im; -0.004131080722171527 + 0.003694979650464384im 0.003669501597474799 - 0.00321210632027655im … -0.0026359786693836776 + 0.0023239334936747793im 0.003655316840240865 - 0.0032554365671528624im; … ; -0.006201675518004896 + 0.001526333547921522im 0.005594571982532088 - 0.001403223810342994im … -0.0042589033621546375 + 0.0010404987734039376im 0.00560495959227617 - 0.0013723985637165074im; 0.0044870664712467236 - 0.00015131829608545574im -0.00396574077905236 + 0.00019590502979949343im … 0.0028791457733077013 - 0.00013723017967008026im -0.003980426440751533 + 0.00015206194412089235im;;; 0.0004501223078787406 + 0.004694204934384397im -0.00029342113943107756 - 0.0031729717811763824im … 0.0001841211132306913 + 0.0016465309000449487im -0.00033608129182909 - 0.003300978902285429im; 0.008739392844385623 - 0.007267237625303789im -0.007440876043853543 + 0.0057906604908564445im … 0.004766647483498972 - 0.0035815187120750826im -0.007407917073035659 + 0.0058904307130260795im; … ; 0.011303664510227712 - 0.0033990266257257463im -0.009860673731675457 + 0.003047810921849931im … 0.0068846558976564986 - 0.0021097670494057987im -0.009880137481832935 + 0.0029900562129192im; -0.009361825088267807 + 0.0008117991384568447im 0.00788215771475476 - 0.0010490652173448774im … -0.005079610176284988 + 0.0008125637128505749im 0.007915866942090528 - 0.0009483560651142094im;;; -0.0006436054651454982 - 0.0068441799982263424im 0.0003233432470574338 + 0.00385320913410145im … -0.00019622488693495403 - 0.0015651953474728757im 0.0004502502529549467 + 0.004221207020564392im; -0.013970889345094464 + 0.010397841766268057im 0.012327611985597876 - 0.008071300269682145im … -0.007508506447990377 + 0.004645196335561822im 0.01224815792057834 - 0.008307956757843848im; … ; -0.01717912205631109 + 0.0061871946801730845im 0.014882536261485979 - 0.005527613517523412im … -0.009829545069228455 + 0.0036254949461631703im 0.014916779298226171 - 0.005427677646468227im; 0.01471551412752405 - 0.002405779039075958im -0.012742578078212458 + 0.0031196010370532386im … 0.007816006927584674 - 0.002193295743867103im -0.01282481607004681 + 0.0028812537583123706im],)])]), [8.51987390747774e-7 4.6272652285573097e-7 … 2.4177787506186894e-8 4.6272652529427953e-7; 4.627265221232728e-7 2.5014477446631556e-7 … 1.3583000576809113e-8 2.501447759535398e-7; … ; 2.4177787266826596e-8 1.3583000544008151e-8 … 1.5659238428090814e-9 1.3583001191894818e-8; 4.6272652527280703e-7 2.501447765947678e-7 … 1.3583001304192003e-8 2.5014477898604313e-7;;; 4.6272652421944565e-7 2.501447762565522e-7 … 1.3583000877613669e-8 2.5014477735204e-7; 2.501447760701516e-7 1.389560924961642e-7 … 1.1279992934616442e-8 1.3895609337088747e-7; … ; 1.3583000693289582e-8 1.1279992805162575e-8 … 7.937758122230635e-9 1.127999321807213e-8; 2.5014477732551254e-7 1.3895609356386352e-7 … 1.127999331770853e-8 1.3895609462563668e-7;;; 2.417778700961401e-8 1.3583000421763382e-8 … 1.5659236571709828e-9 1.3583000724393828e-8; 1.3583000558789425e-8 1.1279992606431192e-8 … 7.937757800967757e-9 1.1279992927898417e-8; … ; 1.5659235969515317e-9 7.937757561361247e-9 … 3.571856756802794e-8 7.937757957834315e-9; 1.3583000655126115e-8 1.1279992764486197e-8 … 7.937757988740696e-9 1.1279993089748074e-8;;; … ;;; 5.791841276843879e-8 2.281602033420745e-8 … 5.620383688508929e-10 2.281601944809923e-8; 2.281601990051948e-8 3.489276124536905e-9 … 1.3458515724993602e-8 3.4892759038019695e-9; … ; 5.620383176219526e-10 1.3458515228283068e-8 … 7.728316461923232e-8 1.3458516297152556e-8; 2.2816019847448106e-8 3.4892761999561824e-9 … 1.3458516149510159e-8 3.489275728040106e-9;;; 2.417778720302071e-8 1.3583000353767282e-8 … 1.5659238475958236e-9 1.3583001198495696e-8; 1.3583000563377406e-8 1.1279992552248764e-8 … 7.937757932334553e-9 1.1279993069385047e-8; … ; 1.5659237580828558e-9 7.937757631831655e-9 … 3.5718568756135484e-8 7.937758593468233e-9; 1.3583000945167193e-8 1.1279992695021387e-8 … 7.93775854028549e-9 1.1279993722052965e-8;;; 4.627265239756778e-7 2.5014477524851947e-7 … 1.358300117742396e-8 2.50144778262525e-7; 2.5014477518704526e-7 1.3895609153492145e-7 … 1.1279992903131661e-8 1.3895609314948094e-7; … ; 1.358300097734248e-8 1.1279992779448549e-8 … 7.93775852824657e-9 1.1279993721202989e-8; 2.501447778203522e-7 1.389560929909101e-7 … 1.12799937536245e-8 1.389560959732653e-7;;;;], Matrix{ComplexF64}[[0.10617599081511188 + 0.23789415600189662im -0.6648451602975075 - 0.6995430842668947im -1.9666711036315254e-7 - 3.0734229748465398e-6im 1.8024890745727587e-7 - 8.210334907962834e-7im; -0.07203816992744243 - 0.1614061663331617im -0.029295707431919095 - 0.030824635595717808im -0.6045201343324907 - 0.05717257993309149im -0.18363522375845223 - 0.03238600501341309im; … ; 0.022887064826794985 + 0.051279945033161835im 0.010595869315826643 + 0.01114886209120458im -0.012925864365207741 + 0.0011595297700987552im 0.0013368744055087952 - 0.0009428500698914718im; -0.04103666402912093 - 0.09194529272093527im -0.02099123893929264 - 0.022086761413020427im 0.03928401199179141 - 0.009631697149915035im -0.01725901142259408 + 0.0035736710398487654im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239440999501625, [[-0.5541104208731916, -0.010678399116840877, 0.14495668357373737, 0.1449567771304171]], 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.6206945992167514e-7, -6.293617172387298e-7, -2.295121316948226e-7, 4.994026279787937e-7, 7.923312810999209e-7, 5.85107287977264e-7, 6.516971956225319e-7, 1.0189310161087874e-6, 3.638592286229113e-7, -2.2010936348137644e-6  …  2.7810185894970577e-6, 1.2016957030498723e-6, 2.4456418054721007e-8, 1.0078549150506978e-7, 5.535308836549616e-7, 4.5341736028854977e-7, 4.631956708119154e-9, -1.0437050548660454e-7, 1.1211100315703381e-7, 1.255086789720135e-7], resid_history = [0.2493920944987058, 0.0037664991303411146, 0.0002851522401517166, 4.682282659013154e-6, 1.791669609923396e-8, 4.8608943972239946e-11, 3.1006843891938975e-8, 9.898228694437727e-10], converged = true, n_iter = 8, residual_norm = 9.898228694437727e-10, maxiter = 100, tol = 1.0e-8, s = 0.93496943637638, restart_history = [6], stage = :finalize, krylovdim = 20, y = [2.4965234710465765e-7, -3.0583447873640333e-8, 9.898228694437727e-10, -0.00020437208781900752, 4.687720133922055e-6, -1.8099320856996816e-8, 4.8608943972239946e-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}}([[2.169988971540852e-7, 1.593267439087628e-6, 5.147922752957199e-7, -7.521519071315905e-7, -7.222796409894316e-7, -3.323660215147143e-7, -2.2131702231104031e-7, -7.616693783483093e-8, 8.43281240560826e-8, -1.0011101151775546e-6  …  4.707123047520352e-7, 7.1796924647737195e-9, -6.289315254041142e-9, -4.372578119863343e-8, -2.5260981337497026e-7, -2.107803773116925e-7, -3.0304423650033367e-9, 1.1657074544824217e-7, -2.0935633026108003e-7, -4.6339060577515134e-7], [1.0574080790163052e-6, -2.790747716616036e-7, -1.5762701982214622e-7, 1.7177314124097978e-7, 1.2820689267008422e-8, -1.0772334594479368e-7, -2.2201552920805232e-7, -5.478247915267034e-7, -2.744009582476184e-7, 1.980171244746004e-6  …  -1.2084703830697115e-7, 1.446320213700593e-7, 8.982913353331202e-9, 5.467104053064588e-8, 3.337148936739743e-7, 3.0879137221857124e-7, 4.424175970550275e-9, -1.6150396346203396e-7, 2.990608145724521e-7, 8.599213633409852e-7]]), H = [1.1047856408016759 0.007786289277877716 … 0.0 0.0; 0.12510020477598174 1.0201649821191656 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1118459305841624 0.12252149781773115 … 0.0 0.0; 0.0 1.0133272439583796 … 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.9257125358461429
Interacting polarizability:     1.773654865112382

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 = 9.37e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 8.58e-10
[ Info: GMRES linsolve in iteration 2; step 2: normres = 3.18e-11
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 3.95e-12
└ * number of operations = 12
Non-interacting polarizability: 1.9257125358460743
Interacting polarizability:     1.773654864070938

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