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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.0001345730118606202
Displaced dipole:  0.017612221122203797
Polarizability :   1.7746794134064414

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.0025656222343223858 - 0.002068937494571041im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.19623348934645776 + 0.2392041172225694im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.02590655522786026 + 0.03186215976653495im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.05056576737577731 - 0.06319885238446904im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620695240717607e-7 -2.502015753920234e-7 … -4.880645870141571e-8 -2.502015758747951e-7; -6.293617719244565e-7 -4.793037305057789e-7 … -1.207785034449411e-7 -4.79303730442531e-7; … ; 1.3626264025736363e-7 1.1356233749497088e-7 … 4.7119967607444894e-8 1.1356234095725818e-7; 1.330547260604224e-7 1.3632638610652903e-7 … 5.2006873906701336e-8 1.3632638743324875e-7;;; -2.5020161375434576e-7 -1.7404915735367293e-7 … -3.5564974756989314e-8 -1.7404915723845506e-7; -4.793037613128154e-7 -3.6797330577597425e-7 … -1.1202902162499563e-7 -3.6797330468542917e-7; … ; 1.1356233412112562e-7 1.1211102399149343e-7 … 1.0980690252199178e-7 1.1211102997449602e-7; 1.363263611080853e-7 1.2550863704273612e-7 … 5.0949105795253604e-8 1.255086392753082e-7;;; -4.880647700607202e-8 -3.556498178988368e-8 … -1.1412174376055802e-8 -3.556498178817332e-8; -1.2077852354078193e-7 -1.1202903293313639e-7 … -9.78623354264248e-8 -1.1202903084771214e-7; … ; 4.7119964416349896e-8 1.0980689081567377e-7 … 2.4938237302966643e-7 1.0980689880830323e-7; 5.200686336964346e-8 5.094909918619428e-8 … 4.910630227888875e-8 5.094910137445983e-8;;; … ;;; 6.475757502465928e-8 4.110330954876388e-8 … -6.681712405502822e-9 4.1103308250334155e-8; 1.6897551395449424e-7 6.653279754108549e-8 … -1.329445741961441e-7 6.653278917503242e-8; … ; 3.2677721669656814e-8 1.6017377966983567e-7 … 3.847037877278427e-7 1.6017378449943093e-7; -9.141991288728745e-8 -3.583032387596101e-8 … 7.037290269361153e-8 -3.583032031903658e-8;;; -4.880644033081824e-8 -3.556495326821887e-8 … -1.141216476601388e-8 -3.5564954979871695e-8; -1.2077848708225777e-7 -1.120289967571679e-7 … -9.786231431860196e-8 -1.120290021425832e-7; … ; 4.7119966171183516e-8 1.0980690416234258e-7 … 2.4938239688203874e-7 1.09806906988414e-7; 5.200688382784571e-8 5.0949117517948024e-8 … 4.9106315676379844e-8 5.094911907132477e-8;;; -2.5020153717475626e-7 -1.7404910060308427e-7 … -3.556496130005863e-8 -1.7404910203155936e-7; -4.793036998419908e-7 -3.679732578381184e-7 … -1.1202900703084606e-7 -3.679732599834123e-7; … ; 1.1356234407313618e-7 1.1211103409539741e-7 … 1.0980690314663184e-7 1.121110361332953e-7; 1.3632641276895463e-7 1.255086754129652e-7 … 5.094911375198497e-8 1.2550867609972797e-7;;;;], δHtotψ = Matrix{ComplexF64}[[0.0014923899293459554 + 0.0012034746908181945im 0.0020376358028757475 + 0.09527524416683675im 0.6140592301577062 + 0.8772285059700216im 0.6601957285234983 + 0.4039578673075058im; 0.12587816185816894 - 0.15322622679553366im 1.6257428847822983 + 0.06029298703949802im 0.1358116589052283 + 0.1933664367055941im 0.14526424709921137 + 0.08998012705432291im; … ; 0.031027985976558207 - 0.03795019431345833im -0.02068046574183312 + 0.0004942590949884933im -0.005337310275370396 - 0.012191129413750503im -0.0045055468402591305 - 0.015606519632836314im; -0.045219396526381614 + 0.05650905127710346im 0.02558635662289357 - 0.0004898638123510527im 0.04073402046544503 + 0.024567745726972116im 0.0921683969123859 - 0.08538684819684071im]], δVind = [0.0019654323211393612 0.0020810010628937894 … 0.0031536651887929135 0.002081001058023448; 0.0038891315300771277 0.004779270093836865 … 0.010276670254518096 0.004779270089583015; … ; -0.006150846358890565 -0.008269601088330011 … -0.01730295720194659 -0.008269600996092552; 0.00023868251796291848 -0.0003986970652156321 … -0.003664693539319784 -0.00039869706787939757;;; 0.0020810013688979004 0.0022232291257175364 … 0.0034201203008264584 0.002223229111514088; 0.004779270429345063 0.005734954593335434 … 0.010903596959986309 0.005734954584861928; … ; -0.00826960044385633 -0.009601004777719212 … -0.01262310680563131 -0.009601004592807454; -0.00039869672337615155 -0.0010730122458791995 … -0.0042306785404898324 -0.001073012253311089;;; 0.0031536662187207315 0.003420120895757842 … 0.004832504638390685 0.0034201208581631413; 0.010276671180013271 0.010903597369404877 … 0.012272140321544016 0.010903597421015007; … ; -0.017302955409195457 -0.012623106462515498 … -0.009620262865808478 -0.012623106132833318; -0.0036646923115438287 -0.004230677929386547 … -0.00549110490673323 -0.004230677896366822;;; … ;;; -0.0023365598574564277 -0.002815462656529535 … 0.0056933100360609615 -0.0028154627503065763; -0.012913094390644798 -0.01784118661670158 … 0.01149833673717663 -0.017841187489563565; … ; -0.02517444036773945 -0.012810885350577027 … -0.008395875591217614 -0.012810885162485956; 0.007597981146075368 0.010225418642692509 … -0.00544973351354307 0.010225418920706162;;; 0.003153664121320735 0.003420118582515869 … 0.004832501414855772 0.003420118562910322; 0.010276669135925369 0.010903595354119328 … 0.012272137857892936 0.010903595185760447; … ; -0.017302959729469887 -0.012623109257453273 … -0.009620264951990514 -0.012623109106883046; -0.003664694805816735 -0.0042306805398938515 … -0.005491107633707176 -0.0042306804526047025;;; 0.002081000746961772 0.002223228424483108 … 0.003420119178158734 0.0022232284273266912; 0.004779269744456735 0.005734953836652361 … 0.01090359586656355 0.0057349538331814475; … ; -0.00826960166129485 -0.009601006008400895 … -0.012623108435680368 -0.009601005937091856; -0.0003986974142175671 -0.0010730130243850613 … -0.004230679824882257 -0.0010730130164478822;;;;], δρ0 = [-3.602141404759742e-7 -2.4888436359920893e-7 … -4.852840235148489e-8 -2.4888436181108183e-7; -5.972252129068397e-7 -4.60697976333733e-7 … -1.1905776430434598e-7 -4.6069797454418944e-7; … ; 1.2735795301695086e-7 1.0836345183304915e-7 … 4.640725648180965e-8 1.0836345553590897e-7; 1.0352008560414102e-7 1.195796324875563e-7 … 5.0695090748929896e-8 1.195796353780187e-7;;; -2.488843645516079e-7 -1.7310679789020764e-7 … -3.535700844639113e-8 -1.7310679613972633e-7; -4.6069797956706274e-7 -3.572172071802939e-7 … -1.1091448026233643e-7 -3.5721720485553127e-7; … ; 1.0836345458358018e-7 1.0846503529727015e-7 … 1.086224596365719e-7 1.0846504145453384e-7; 1.1957963455738356e-7 1.1611149347885874e-7 … 5.020833972392998e-8 1.1611149683808637e-7;;; -4.852840486278319e-8 -3.535700994470156e-8 … -1.1341339132899801e-8 -3.535700957782509e-8; -1.1905777248921081e-7 -1.1091448649567287e-7 … -9.770677926511359e-8 -1.1091448410041059e-7; … ; 4.6407257218559626e-8 1.0862245214208126e-7 … 2.484525760582459e-7 1.0862246034726935e-7; 5.06950919030032e-8 5.0208338121572585e-8 … 4.926862846945001e-8 5.0208340625004275e-8;;; … ;;; 6.43457156471766e-8 4.083862494558328e-8 … -6.637952281351597e-9 4.0838623205112524e-8; 1.6975032430184692e-7 6.683079909391769e-8 … -1.3340492340695082e-7 6.683079053129287e-8; … ; 3.271007310686911e-8 1.6023533831105177e-7 … 3.8390158688373163e-7 1.602353435341406e-7; -9.271364171926995e-8 -3.633617736564847e-8 … 7.126314011414286e-8 -3.633617395694926e-8;;; -4.852840098219803e-8 -3.535700584793276e-8 … -1.1341337702542178e-8 -3.535700719221227e-8; -1.1905776064813391e-7 -1.1091447261151147e-7 … -9.770677656460723e-8 -1.1091447763857408e-7; … ; 4.640725090395599e-8 1.0862244722411038e-7 … 2.4845256134347705e-7 1.0862245030663364e-7; 5.069508809351522e-8 5.020833443066776e-8 … 4.926862355836607e-8 5.020833630455894e-8;;; -2.4888436074987274e-7 -1.731067941928126e-7 … -3.5357007124754456e-8 -1.7310679399582962e-7; -4.60697971555184e-7 -3.5721719880635437e-7 … -1.1091447676222393e-7 -3.572171996948527e-7; … ; 1.0836345252457748e-7 1.0846503445675543e-7 … 1.0862245121791528e-7 1.0846503677498748e-7; 1.1957963347361485e-7 1.1611149281107412e-7 … 5.02083367555067e-8 1.1611149466086583e-7;;;;], δeigenvalues = [[0.0004950047443470925, 0.09679859376823714, 0.040617789059406766, 0.044788695354550115]], δ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.15475965285951745 -0.15461620273497145 … -0.15418753103441035 -0.154616202734781; -0.15461620273537008 -0.15447229893472475 … -0.15404229471195346 -0.15447229893453468; … ; -0.15418753103362884 -0.15404229471076558 … -0.1536083373814264 -0.15404229471057188; -0.15461620273438634 -0.15447229893373862 … -0.1540422947109559 -0.15447229893354694;;; -0.1546162027345099 -0.15447229893386272 … -0.15404229471108152 -0.1544722989336711; -0.15447229893426218 -0.1543279431256297 … -0.15389660286425924 -0.15432794312543863; … ; -0.15404229471029887 -0.15389660286306955 … -0.15346128704474468 -0.15389660286287418; -0.15447229893327574 -0.15432794312464104 … -0.15389660286325868 -0.154327943124448;;; -0.15418753103387178 -0.15404229471100928 … -0.15360833738167182 -0.15404229471081543; -0.15404229471141223 -0.15389660286418524 … -0.1534612870458734 -0.1538966028639919; … ; -0.15360833738088323 -0.1534612870446745 … -0.1530219131690276 -0.15346128704447615; -0.1540422947104167 -0.1538966028631877 … -0.15346128704486345 -0.15389660286299214;;; … ;;; -0.15347880589151075 -0.15333139200146606 … -0.15289092496419207 -0.15333139200127233; -0.153331392001878 -0.15318350577987022 … -0.1527416288164353 -0.1531835057796761; … ; -0.1528909249633858 -0.1527416288152097 … -0.15229553962631506 -0.15274162881501294; -0.1533313920008646 -0.15318350577885403 … -0.1527416288154089 -0.15318350577865927;;; -0.15418753103532598 -0.15404229471246653 … -0.1536083373831458 -0.15404229471227512; -0.15404229471287176 -0.15389660286564819 … -0.15346128704735235 -0.15389660286545664; … ; -0.15360833738235244 -0.15346128704614642 … -0.15302191317051933 -0.15346128704595213; -0.1540422947118741 -0.15389660286464807 … -0.1534612870463418 -0.15389660286445564;;; -0.15461620273524232 -0.15447229893459644 … -0.15404229471182324 -0.1544722989344061; -0.15447229893499712 -0.1543279431263662 … -0.1538966028650034 -0.15432794312617598; … ; -0.1540422947110381 -0.15389660286381013 … -0.15346128704549516 -0.1538966028636168; -0.15447229893400946 -0.15432794312537607 … -0.1538966028640025 -0.15432794312518464]), 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.015163898501150713 -0.012550130022991247 … -0.004966869621910768 -0.012550130046951654; -0.012550130070060651 -0.010362596110682144 … -0.004135160608459811 -0.010362596113330135; … ; -0.004966869517126129 -0.004135160464678091 … -0.0020697827020692665 -0.0041351605390730455; -0.012550130011466893 -0.01036259604009183 … -0.004135160612633599 -0.010362596078075055;;; -0.012550130072138769 -0.010362596104478568 … -0.0041351606689440205 -0.010362596130949676; -0.010362596158787258 -0.008622040205137746 … -0.003897513157400548 -0.008622040198255613; … ; -0.004135160519650765 -0.003897512976178148 … -0.003484108952113004 -0.0038975130997078774; -0.010362596085722873 -0.008622040110623922 … -0.003897513195157544 -0.008622040167899837;;; -0.004966869740846884 -0.0041351607177476084 … -0.0020697829017452514 -0.00413516073858733; -0.004135160776051334 -0.003897513280047875 … -0.0034841089577210206 -0.0038975132388489823; … ; -0.0020697827230636485 -0.003484108815277477 … -0.005620581637331288 -0.0034841089680045126; -0.004135160690083908 -0.0038975131571697317 … -0.0034841090557383024 -0.003897513245416921;;; … ;;; -0.006547897036994424 -0.004876338936777043 … -0.0014866472841467788 -0.00487633879020478; -0.004876338871031613 -0.002677986967007895 … -0.004123056128143732 -0.002677986742251628; … ; -0.0014866469626001325 -0.004123055883100536 … -0.007170896597643737 -0.0041230559636351635; -0.004876338878135916 -0.002677986965697801 … -0.004123056065982123 -0.0026779868015212012;;; -0.004966869526773483 -0.0041351604642784374 … -0.002069782726073584 -0.004135160567357273; -0.004135160531600001 -0.0038975129908377024 … -0.0034841088985540693 -0.003897513102418129; … ; -0.0020697825493848883 -0.003484108718542724 … -0.005620581451990349 -0.0034841087761197415; -0.00413516049785916 -0.003897512968259129 … -0.0034841088571009853 -0.0038975130514668404;;; -0.012550130003703305 -0.010362596032077264 … -0.004135160566596436 -0.01036259606455915; -0.010362596074455945 -0.008622040108150807 … -0.0038975130817330717 -0.008622040132799972; … ; -0.004135160476144558 -0.003897512959308679 … -0.0034841087899448843 -0.003897513002594926; -0.010362596032050947 -0.008622040063329712 … -0.0038975130645415067 -0.008622040097163618])], 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.00953289329181356 -0.0069191761841548086 … 0.0006639997063572108 -0.006919176207924729; -0.006919176231622699 -0.004731679637719073 … 0.0014956953455665094 -0.004731679640176994; … ; 0.0006639998119234471 0.0014956954905360293 … 0.003561043869068154 0.001495695416334781; -0.006919176172045145 -0.004731679566142604 … 0.0014956953423903123 -0.004731679603934148;;; -0.00691917623284081 -0.004731679630653631 … 0.0014956952859541024 -0.004731679656933087; -0.004731679685361612 -0.002991148828109755 … 0.0017333311218725527 -0.002991148821036552; … ; 0.0014956954360300931 0.0017333313042845844 … 0.0021467134536686896 0.0017333311809502264; -0.004731679611310767 -0.00299114873260725 … 0.0017333310851162008 -0.002991148789690125;;; 0.000663999587959636 0.0014956952372227623 … 0.0035610436691467543 0.0014956952165768858; 0.0014956951785161645 0.0017333309992992496 … 0.0021467134469319925 0.0017333310406914878; … ; 0.0035610438486169765 0.0021467135905743825 … 1.0227612596471619e-5 0.0021467134380456884; 0.0014956952654792107 0.0017333311231749284 … 0.0021467133499247083 0.0017333310351233047;;; … ;;; -0.0009170690029354969 0.0007544796415652573 … 0.0041441550443073365 0.0007544797883312542; 0.0007544797068987944 0.0029528264566951066 … 0.001507743444563392 0.002952826681645496; … ; 0.004144155366660282 0.0015077436908321352 … -0.0015401037014268074 0.0015077436104942669; 0.0007544797008079314 0.0029528264590214432 … 0.0015077435077514847 0.002952826623392804;;; 0.0006639998005788395 0.0014956954892346545 … 0.003561043843344434 0.0014956953863472216; 0.001495695421507971 0.0017333312870464815 … 0.002146713504619988 0.001733331175657596; … ; 0.0035610440208264952 0.002146713685837202 … 1.0227796445659229e-5 0.0021467136284544735; 0.0014956954562465137 0.0017333313106251717 … 0.00214671354708368 0.0017333312276098896;;; -0.00691917616513776 -0.004731679558986046 … 0.0014956953875599746 -0.0047316795912775555; -0.004731679601765239 -0.0029911487318592823 … 0.001733331196795902 -0.002991148756318238; … ; 0.0014956954787970298 0.001733331320413479 … 0.0021467136150862985 0.0017333312773205496; -0.004731679558372559 -0.002991148686048063 … 0.0017333312149884162 -0.002991148719690539]), 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.9275421529004636e-6 + 1.8998927907995878e-6im 0.005571081251263506 + 0.004249226073827045im … 0.003857892373196457 + 0.0029457328841398633im -0.005626197498794316 - 0.004295855055158282im; 0.01669983316427083 - 0.020648241328289154im -0.01682883006902971 + 0.014305935349428726im … 0.004837693070317264 - 0.011314219600428586im -0.010268313249698052 + 0.019312834835212296im; … ; 0.017514495858538294 - 0.021643412609183482im -0.013574434424674816 + 0.019012586711740647im … 0.010374932336724672 - 0.010271484716525378im -0.01582592396820303 + 0.01729386990871655im; -0.016647559696455434 + 0.02070499584487141im 0.010280128073177386 - 0.019310727342970603im … -0.010144328171924741 + 0.007264342782323098im 0.016840335705044247 - 0.01430376225764925im;;; -0.001923915870243538 + 0.006570615114404246im -0.00213369914572184 - 0.0063379935295723965im … -0.0030403312785604575 - 0.0007041188058128514im 0.004426590357232303 - 0.001331253591376376im; -0.012438964309452051 + 0.012944545948879553im 0.012374151069526903 - 0.00955123648167574im … -0.00397506302929234 + 0.0082438104360796im 0.008150674483290422 - 0.012774730621710294im; … ; -0.015079174228398424 + 0.01948081845003074im 0.011816785702275276 - 0.0170230061888349im … -0.009141987720223246 + 0.009540156852298275im 0.013584170461214531 - 0.01567384067871125im; 0.014668622433053004 - 0.02067228865715851im -0.009619635903339148 + 0.017723175306354505im … 0.008551421350044944 - 0.007289966294588013im -0.013842958282409957 + 0.01449964437948276im;;; 0.001333377650586355 - 0.004463755084663746im 0.0001997048967151051 + 0.003916481743399727im … 0.0017204594543631115 - 0.0004953974108148137im -0.002051837544351367 + 0.0021979139716294437im; 0.0065723392261553305 - 0.00622538703896245im -0.006452037332058434 + 0.004812530349319783im … 0.0025111512544810166 - 0.00449418809812228im -0.004684604490728873 + 0.006161626053492406im; … ; 0.009540078988175674 - 0.012717467860897353im -0.007763353895613268 + 0.011378854313215468im … 0.006387823867943586 - 0.007002088486349757im -0.008779237324040333 + 0.010603356302504807im; -0.008409071200111223 + 0.012353302550516976im 0.006074268514319541 - 0.010721292922043291im … -0.005402460316921766 + 0.005319306806176943im 0.007841665951550125 - 0.009372195336847137im;;; … ;;; 0.0005292253458978259 - 0.0020190062974425275im -0.0008848698705846048 + 0.0013699919043565407im … -0.0003284522869615868 - 0.001608563252095983im 4.26112129312846e-9 + 0.002045472486716323im; -0.003904067473079204 + 0.00599908561921319im 0.0038187641595643063 - 0.004970675879343939im … -0.0019020550154317016 + 0.004167319821701633im 0.0030508458995844303 - 0.0055568709794348425im; … ; -0.0045851453730328525 + 0.004900175511424408im 0.0038706781254013137 - 0.0046351100816443315im … -0.003557119417720264 + 0.003051054067295263im 0.00441103813870896 - 0.004222555986829998im; 0.003018228102324779 - 0.0025846252364107635im -0.0022961551217194863 + 0.0026119257820094556im … 0.002498384084990499 - 0.0012892973092673743im -0.003064058258612568 + 0.0020256917833083154im;;; -0.0012345970345237888 + 0.0045479268791739375im 0.001962613109793105 - 0.002270702898100215im … 0.0008737806172589133 + 0.0024768724700570564im -0.00028905172079115524 - 0.003989246350788077im; 0.008335879980131609 - 0.012412438726342644im -0.007768109061838934 + 0.009430271832885466im … 0.003220892120536731 - 0.006985205240221455im -0.006000588902506003 + 0.010779356003654728im; … ; 0.008693836045145821 - 0.00974514563075628im -0.007054047666949127 + 0.008887801088867376im … 0.005924933361811694 - 0.005376568428905771im -0.008069969282010994 + 0.008112260220420862im; -0.006646926235514617 + 0.006166027401140267im 0.004759044481241216 - 0.006103427528432603im … -0.004692994937532553 + 0.002828239263599581im 0.006526527505934985 - 0.004754329501128472im;;; 0.0018427046554668207 - 0.006645614886049029im -0.0043400413998372096 + 0.0014048771280842044im … -0.002283495097896578 - 0.003361219956617971im 0.002220454129588472 + 0.006411554681854336im; -0.014646421916298332 + 0.020687170809263924im 0.01379494727587403 - 0.014535634365067749im … -0.0045693492143475765 + 0.010329537457750133im 0.009571341188972293 - 0.01775909377605392im; … ; -0.014322973029848921 + 0.016823610581150582im 0.011223032686243651 - 0.014937181981568029im … -0.008800648632164274 + 0.00834131117090891im 0.01299046002409662 - 0.013588016605375754im; 0.01246361445415772 - 0.012929265749179561im -0.008200148905586938 + 0.01273856001514486im … 0.007957474663508476 - 0.005204194060274156im -0.012423600218334418 + 0.009515069650527225im],)])]), [8.519874862588016e-7 4.6272658986279377e-7 … 2.4177800463071518e-8 4.627265927065865e-7; 4.627265954493353e-7 2.50144827478646e-7 … 1.3583011131939169e-8 2.501448276834548e-7; … ; 2.417779885526974e-8 1.3583009648601248e-8 … 1.5659272915511383e-9 1.3583010416105425e-8; 4.627265884950007e-7 2.5014482201894337e-7 … 1.358301117500175e-8 2.5014482495669557e-7;;; 4.627265956959821e-7 2.501448269988503e-7 … 1.3583011755929725e-8 2.5014482904621834e-7; 2.5014483119927603e-7 1.3895613352592807e-7 … 1.1280003335865162e-8 1.3895613317209052e-7; … ; 1.3583010215729878e-8 1.1280001690437678e-8 … 7.937767253308257e-9 1.128000281203773e-8; 2.5014482554821457e-7 1.3895612866657444e-7 … 1.1280003678685761e-8 1.3895613161137286e-7;;; 2.4177802288006945e-8 1.3583012259418674e-8 … 1.5659277604387957e-9 1.3583012474415223e-8; 1.3583012860916104e-8 1.128000444945267e-8 … 7.937767293309417e-9 1.1280004075382943e-8; … ; 1.5659273408529466e-9 7.937766277302116e-9 … 3.571859050615851e-8 7.937767366654798e-9; 1.3583011974024557e-8 1.1280003333768284e-8 … 7.937767992435792e-9 1.128000413501697e-8;;; … ;;; 5.791839946871687e-8 2.281601089681719e-8 … 5.620406413358369e-10 2.281600873590214e-8; 2.281600992752956e-8 3.4892718730949084e-9 … 1.3458528306111434e-8 3.4892709604462458e-9; … ; 5.6204026579103e-10 1.3458525794056998e-8 … 7.728319452575639e-8 1.3458526619655235e-8; 2.2816010032270573e-8 3.4892718677734264e-9 … 1.3458527668866821e-8 3.489271201118605e-9;;; 2.4177799003299347e-8 1.3583009644480668e-8 … 1.5659273479193988e-9 1.358301070790413e-8; 1.3583010339009048e-8 1.1280001823538984e-8 … 7.93776687129045e-9 1.1280002836648475e-8; … ; 1.565926933013036e-9 7.937765587321743e-9 … 3.5718586784039224e-8 7.937765998003973e-9; 1.358300999091643e-8 1.1280001618536924e-8 … 7.937766575615819e-9 1.128000237402634e-8;;; 4.6272658757356113e-7 2.501448213990724e-7 … 1.3583010700053256e-8 2.501448239113377e-7; 2.50144824676786e-7 1.389561285394308e-7 … 1.1280002648832488e-8 1.3895612980674058e-7; … ; 1.35830097668944e-8 1.1280001537268172e-8 … 7.937766096615031e-9 1.1280001930291048e-8; 2.5014482139702754e-7 1.3895612623498536e-7 … 1.1280002492742564e-8 1.3895612797453502e-7;;;;], Matrix{ComplexF64}[[0.20279098385777336 + 0.163532218505866im 0.02063525386925563 + 0.9648594707271072im -3.880282342330608e-6 + 1.4783068089943856e-7im 1.2753259893319336e-6 + 1.4612391823782973e-6im; -0.13758940496896693 - 0.11095316070143463im 0.0009092683991275696 + 0.042515497284834815im 0.27795725698330614 - 0.19441481849046074im 0.126732754735988 - 0.20961762921012883im; … ; 0.043713182219409584 + 0.035250648407479726im -0.00032887128932697417 - 0.015377306806543983im 0.006021155388029182 - 0.005323693328752199im -0.00191539874186398 - 0.008280504784264613im; -0.07837803501259233 - 0.06320465390018702im 0.0006515191608355281 + 0.030463622403391684im -0.018588117227432307 + 0.019153993782387593im 0.017549704885263227 + 0.03479005645164065im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239440856413445, [[-0.5541104268349304, -0.010678390293338511, 0.14495675724731785, 0.14495685555286572]], 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.620695240717607e-7, -6.293617719244565e-7, -2.2951220053798912e-7, 4.99402515358235e-7, 7.923311447798266e-7, 5.851071589310742e-7, 6.516970298698118e-7, 1.0189307697392273e-6, 3.638590751047232e-7, -2.2010933524082127e-6  …  2.78101813785013e-6, 1.2016955582628994e-6, 2.445647907422291e-8, 1.0078555342890346e-7, 5.535308896044322e-7, 4.534173776210058e-7, 4.632013822303556e-9, -1.0437044339613065e-7, 1.121110361332953e-7, 1.2550867609972797e-7], resid_history = [0.2493920935251157, 0.003766548708641359, 0.00028527788931128234, 4.680425504998838e-6, 4.5587246009379846e-8, 5.737547508899833e-11, 4.33293637997915e-8, 1.0617159958132852e-9], converged = true, n_iter = 8, residual_norm = 1.0617159958132852e-9, maxiter = 100, tol = 1.0e-8, s = 0.9356498748619111, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.068755583981175e-7, -4.2480190601058876e-8, 1.0617159958132852e-9, -0.00020446260764942544, 4.6829275755848155e-6, -4.608177256798144e-8, 5.737547508899833e-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.6396864750610062e-7, -1.9978082283790358e-6, -7.113073964254031e-7, 1.1530947991674596e-6, 1.3870140341221156e-6, 9.931340748571954e-7, 1.2123392117402048e-6, 1.8033159651722776e-6, 4.854816506040023e-7, -2.1306713204117308e-6  …  3.0332883301073225e-6, 1.17920766745241e-6, 2.8833813807335063e-8, 1.4564400314271125e-7, 7.966393978248516e-7, 5.376121422562417e-7, 4.805876277378232e-9, -1.349625354660228e-7, 2.2685433625146863e-7, 4.548747481580518e-7], [5.095550067713171e-7, 1.5416399737183805e-6, 5.198225047986137e-7, -8.514583526346963e-7, -1.0327431633284178e-6, -7.196016866362847e-7, -8.533439360437146e-7, -1.2679001924304982e-6, -3.5033738328041013e-7, 1.5753979560496196e-6  …  -1.4560826438044295e-6, -5.365051222476521e-7, -1.352708015303e-8, -7.128577278208728e-8, -3.8608188822024725e-7, -2.3730235064124366e-7, -1.787904632432747e-9, 4.697187825449598e-8, -7.695926964703844e-8, -7.391394881747844e-8]]), H = [1.0848521211758975 0.015207844729790474 … 0.0 0.0; 0.1426867057117202 1.0301570775719122 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.094195421671425 0.14941387903839146 … 0.0 0.0; 0.0 1.0196836082172813 … 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.9257124985341523
Interacting polarizability:     1.7736548313333713

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

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