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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457334070088152
Displaced dipole:  0.017612221451565976
Polarizability :   1.774679479226686

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.003167111252571807 + 0.0009123169852574026im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.08314273693212132 + 0.2980159043707467im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.011207833753024844 + 0.03950609190964163im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.022701773000066656 - 0.07768926139984322im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206935819507673e-7 -2.5020144883852235e-7 … -4.8806420856905773e-8 -2.502014511677007e-7; -6.293615758346271e-7 -4.793035661848509e-7 … -1.207784270945501e-7 -4.793035720894514e-7; … ; 1.36262610073836e-7 1.1356229766911463e-7 … 4.7119913456000845e-8 1.1356230314281891e-7; 1.3305477883738748e-7 1.363264151558321e-7 … 5.200686257256468e-8 1.3632641822703525e-7;;; -2.502014232732569e-7 -1.7404901283912292e-7 … -3.5564930605308146e-8 -1.7404901420070456e-7; -4.793035511412489e-7 -3.6797312924020234e-7 … -1.1202893132798373e-7 -3.679731331343904e-7; … ; 1.1356230538929101e-7 1.1211098765682706e-7 … 1.0980684480695525e-7 1.1211099123562306e-7; 1.363264368425446e-7 1.255086844199831e-7 … 5.09490974887945e-8 1.2550868640545168e-7;;; -4.880640746618038e-8 -3.5564925124730955e-8 … -1.1412144132095132e-8 -3.556492531606963e-8; -1.2077841673088707e-7 -1.1202892506285658e-7 … -9.786223317437193e-8 -1.1202892608435106e-7; … ; 4.7119912257965935e-8 1.0980684425826015e-7 … 2.493823213542957e-7 1.0980684586513455e-7; 5.200686955502529e-8 5.094910003472818e-8 … 4.910628957808639e-8 5.0949100394280206e-8;;; … ;;; 6.475761445120175e-8 4.1103338834284416e-8 … -6.681699727642722e-9 4.110333789688493e-8; 1.689755661540984e-7 6.653284904123595e-8 … -1.3294450698802612e-7 6.653284720758384e-8; … ; 3.2677654264135206e-8 1.6017367951098876e-7 … 3.8470363120049237e-7 1.601736865072018e-7; -9.14199029574659e-8 -3.58303362398329e-8 … 7.037284625783117e-8 -3.583033346846695e-8;;; -4.880643144483182e-8 -3.556494234622333e-8 … -1.1412151040590606e-8 -3.556494402274389e-8; -1.2077843151928167e-7 -1.1202893730733569e-7 … -9.786224667410552e-8 -1.120289413991349e-7; … ; 4.711990423334037e-8 1.0980682569699669e-7 … 2.493822940607242e-7 1.0980683312400895e-7; 5.200685073171949e-8 5.0949081565619576e-8 … 4.910627788452013e-8 5.0949084991961656e-8;;; -2.502014768094365e-7 -1.7404905149980712e-7 … -3.556494018567913e-8 -1.7404905387753506e-7; -4.793035877670757e-7 -3.679731553129575e-7 … -1.1202893947395598e-7 -3.679731611895058e-7; … ; 1.1356229458089205e-7 1.1211097635422608e-7 … 1.0980683948841911e-7 1.1211098325969032e-7; 1.3632639625718991e-7 1.255086535802004e-7 … 5.0949090418060666e-8 1.2550865713050196e-7;;;;], δHψ = Matrix{ComplexF64}[[0.0018422686591328843 - 0.0005306832014437793im -0.08899191768757947 - 0.03408759440797666im 1.8638146139825933 - 0.006848741119974482im -0.12149994808984059 + 0.07041387074551109im; -0.053156553390340776 - 0.19104441745165318im -0.6696845591785694 + 1.4826320372374746im 0.4119263848568348 - 0.0021122926499424623im -0.026304785634817376 + 0.01340124541717945im; … ; -0.013250957480072365 - 0.04719497088623257im 0.007350947056952414 - 0.019336222044768703im -0.03127198777526625 - 0.00506234506826027im 0.0042679932842284345 - 0.004642440019983925im; 0.02029516985062343 + 0.06947066132131113im -0.009207399497038861 + 0.023877293040573736im -0.015274840677439542 - 0.06871406098652579im 0.039068647751702434 - 0.03168170326691602im]], δVind = [0.0019654316004924562 0.0020810002643447713 … 0.003153664049496301 0.002081000262563172; 0.0038891305494062316 0.0047792690413836925 … 0.010276669825176016 0.004779269053165421; … ; -0.006150847714095397 -0.008269603090043438 … -0.017302966279216805 -0.008269602998109398; 0.00023868195892140748 -0.0003986977013869539 … -0.0036646952222857496 -0.0003986977164078655;;; 0.0020810000357965983 0.00222322764237291 … 0.003420118194561934 0.0022232276445163066; 0.004779268794222479 0.005734952855481109 … 0.010903595870639715 0.005734952873725169; … ; -0.008269603496995416 -0.009601008236666958 … -0.012623111764682856 -0.009601008197011129; -0.00039869797067540024 -0.0010730136682563073 … -0.00423068159458216 -0.0010730136849704682;;; 0.0031536632742857405 0.0034201177574010463 … 0.004832501790567827 0.0034201177696905147; 0.010276668951161649 0.010903595401097768 … 0.01227213890487613 0.010903595432447495; … ; -0.01730296813226345 -0.012623112311289715 … -0.0096202663314962 -0.012623112338034152; -0.0036646961596222116 -0.004230682060970374 … -0.0054911098662789905 -0.004230682090312804;;; … ;;; -0.0023365608504439033 -0.0028154637163976598 … 0.005693317222314281 -0.002815463747941567; -0.01291309392962364 -0.01784118309082393 … 0.011498339919857556 -0.017841183358644947; … ; -0.025174453895199518 -0.01281088485750681 … -0.008395874401305032 -0.01281088471659498; 0.007597978558043644 0.010225413427329579 … -0.005449732468952928 0.010225413757762947;;; 0.003153664869637889 0.003420119556922329 … 0.004832504250802829 0.0034201195039944513; 0.010276670769640648 0.010903597322536405 … 0.012272140771464185 0.010903597231541266; … ; -0.017302965672816438 -0.012623110502626197 … -0.00962026476127287 -0.012623110235575004; -0.003664694357502666 -0.004230680250624433 … -0.0054911077827943045 -0.004230680169297083;;; 0.0020810004929449 0.002223228158718786 … 0.003420119044775534 0.0022232281517783137; 0.00477926930047787 0.0057349534280864665 … 0.010903596740218168 0.005734953429880688; … ; -0.008269602620208032 -0.009601007390338624 … -0.012623110653652531 -0.009601007216092219; -0.00039869744668695863 -0.001073013087409722 … -0.004230680613297094 -0.0010730130966404492;;;;], δρ0 = [-3.602139738196703e-7 -2.488842340294131e-7 … -4.8528365750317465e-8 -2.488842394079933e-7; -5.972249810274401e-7 -4.606977852614799e-7 … -1.1905768403655799e-7 -4.606977933365533e-7; … ; 1.2735790890095234e-7 1.0836340264229278e-7 … 4.640719917396917e-8 1.083634073960662e-7; 1.035201030151517e-7 1.1957963871907062e-7 … 5.0695073496326465e-8 1.1957963943289663e-7;;; -2.4888423863583974e-7 -1.7310669871708193e-7 … -3.5356976163586635e-8 -1.7310670236807345e-7; -4.606977921908019e-7 -3.5721704697983903e-7 … -1.1091439640052737e-7 -3.572170525521562e-7; … ; 1.083634051740312e-7 1.0846498068948824e-7 … 1.0862238564387588e-7 1.0846498370221963e-7; 1.195796382807163e-7 1.1611148962445772e-7 … 5.020831622253588e-8 1.161114898738952e-7;;; -4.852836460861767e-8 -3.5356974731638004e-8 … -1.1341316635881843e-8 -3.53569754542947e-8; -1.1905768280676213e-7 -1.1091439380826338e-7 … -9.770669055437926e-8 -1.1091439531024985e-7; … ; 4.640719496584248e-8 1.0862238208623373e-7 … 2.48452472293589e-7 1.0862238327895069e-7; 5.06950713838523e-8 5.020831511411991e-8 … 4.926859474243613e-8 5.020831499099865e-8;;; … ;;; 6.434568778746376e-8 4.08386116317251e-8 … -6.6379333630598095e-9 4.0838611360705484e-8; 1.6975032838308883e-7 6.683083183067177e-8 … -1.3340482046774388e-7 6.683083025135117e-8; … ; 3.271001014988464e-8 1.6023526069593386e-7 … 3.839014819278553e-7 1.6023526720263142e-7; -9.271366777569127e-8 -3.6336204240225764e-8 … 7.126310972906869e-8 -3.633620117595316e-8;;; -4.852836268697344e-8 -3.5356972642385474e-8 … -1.134131703353409e-8 -3.5356974815781576e-8; -1.1905767831415462e-7 -1.1091438855565752e-7 … -9.770668944861038e-8 -1.1091439308248668e-7; … ; 4.64071934689115e-8 1.0862237816803226e-7 … 2.484524758921459e-7 1.0862238514592863e-7; 5.069507190319251e-8 5.02083142072793e-8 … 4.92685976876536e-8 5.020831714182067e-8;;; -2.48884234834477e-7 -1.7310669551690635e-7 … -3.5356976085701486e-8 -1.7310670007923976e-7; -4.606977869815461e-7 -3.572170420794924e-7 … -1.1091439577905145e-7 -3.5721704951240777e-7; … ; 1.0836340406841556e-7 1.084649782267944e-7 … 1.0862238762655725e-7 1.0846498448023931e-7; 1.1957963964634414e-7 1.1611148990110516e-7 … 5.0208317914609216e-8 1.1611149169670618e-7;;;;], δeigenvalues = [[0.000495005246553015, 0.09679862290517006, 0.014274086459241303, 0.05219197103583797]], δ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.15475965283223664 -0.15461620270753784 … -0.1541875310074154 -0.15461620270764495; -0.15461620270785326 -0.15447229890705985 … -0.15404229468471722 -0.1544722989071598; … ; -0.15418753100673563 -0.15404229468371597 … -0.15360833735486176 -0.15404229468384067; -0.15461620270731194 -0.15447229890650926 … -0.1540422946841854 -0.15447229890662467;;; -0.15461620270752188 -0.15447229890672093 … -0.15404229468439593 -0.15447229890683434; -0.1544722989070329 -0.15432794309825137 … -0.15389660283732592 -0.1543279430983573; … ; -0.15404229468372194 -0.15389660283633416 … -0.1534612870185115 -0.1538966028364652; -0.15447229890650432 -0.1543279430977133 … -0.1538966028368073 -0.15432794309783512;;; -0.1541875310071692 -0.15404229468415323 … -0.1536083373552972 -0.15404229468427375; -0.15404229468446373 -0.1538966028370873 … -0.15346128701924378 -0.15389660283720014; … ; -0.1536083373546247 -0.15346128701825523 … -0.1530219131431335 -0.15346128701839376; -0.1540422946839447 -0.15389660283655873 … -0.15346128701873554 -0.15389660283668782;;; … ;;; -0.15347880586334486 -0.1533313919731565 … -0.15289092493631873 -0.15333139197325282; -0.15333139197349427 -0.15318350575134626 … -0.15274162878833517 -0.1531835057514348; … ; -0.15289092493559575 -0.15274162878726769 … -0.15229553959884337 -0.15274162878738112; -0.15333139197289794 -0.15318350575073983 … -0.15274162878774683 -0.1531835057508444;;; -0.15418753100744384 -0.15404229468443906 … -0.15360833735555163 -0.15404229468453784; -0.15404229468476754 -0.1538966028374024 … -0.1534612870195284 -0.1538966028374937; … ; -0.15360833735484727 -0.15346128701848913 … -0.1530219131433331 -0.15346128701860498; -0.1540422946841922 -0.1538966028368174 … -0.15346128701896147 -0.15389660283692438;;; -0.15461620270765997 -0.15447229890686473 … -0.1540422946845237 -0.154472298906967; -0.15447229890718572 -0.15432794309840978 … -0.15389660283746892 -0.15432794309850484; … ; -0.15404229468383385 -0.1538966028364517 … -0.1534612870186119 -0.15389660283657136; -0.15447229890662878 -0.15432794309784342 … -0.15389660283692083 -0.15432794309795392]), 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.015163897859021867 -0.012550129318221371 … -0.004966868668156137 -0.01255012937212336; -0.012550129347753548 -0.010362595296485283 … -0.0041351594660004745 -0.01036259534413283; … ; -0.004966868653717293 -0.004135159439572012 … -0.002069781136121717 -0.004135159548039826; -0.012550129358955376 -0.010362595307545842 … -0.00413515955041605 -0.010362595375626138;;; -0.012550129348355029 -0.010362595307625923 … -0.004135159460642779 -0.010362595337311138; -0.010362595334639006 -0.008622039273817602 … -0.0038975118573993175 -0.008622039304526305; … ; -0.004135159499252973 -0.0038975118795863146 … -0.0034841075747084407 -0.0038975119432673663; -0.010362595341531517 -0.008622039288103792 … -0.0038975119124265403 -0.008622039324848124;;; -0.004966868597064428 -0.004135159408628751 … -0.0020697808844313565 -0.004135159411486011; -0.004135159444407352 -0.0038975118251478105 … -0.0034841074286561264 -0.0038975118368786256; … ; -0.0020697810211400343 -0.003484107508979787 … -0.0056205803788800905 -0.0034841075300135085; -0.004135159451978214 -0.003897511852352783 … -0.00348410746596595 -0.0038975118560253567;;; … ;;; -0.006547897559879683 -0.004876339675750293 … -0.0014866450553359639 -0.004876339600014205; -0.004876339656501695 -0.0026779882094372515 … -0.004123054693357958 -0.0026779881543990357; … ; -0.0014866451612730305 -0.004123054744699518 … -0.007170895605088391 -0.004123054845445463; -0.004876339580259248 -0.002677988144468852 … -0.004123054819107781 -0.002677988000185623;;; -0.0049668685771842745 -0.004135159339061989 … -0.002069780982277237 -0.004135159458094277; -0.004135159382275713 -0.003897511741964957 … -0.0034841074323495327 -0.003897511821164161; … ; -0.0020697809742806854 -0.0034841074248971966 … -0.005620580421691724 -0.003484107562808429; -0.00413515944162228 -0.003897511778702122 … -0.0034841075693804457 -0.003897511926833844;;; -0.012550129338811677 -0.010362595278697439 … -0.0041351595056049635 -0.01036259535432161; -0.010362595312012234 -0.008622039239571277 … -0.0038975118618244453 -0.008622039300218631; … ; -0.004135159469117139 -0.003897511820893898 … -0.0034841076129313066 -0.0038975119579180866; -0.010362595333147777 -0.008622039253930583 … -0.0038975119768746777 -0.008622039351166138])], 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.009532892622403897 -0.0069191754519513215 … 0.0006640006871068034 -0.006919175505960391; -0.006919175481798774 -0.004731678795857314 … 0.0014956965152620868 -0.004731678843604809; … ; 0.0006640007022254927 0.001495696542691721 … 0.003561045461580343 0.0014956964340992005; -0.006919175492459229 -0.004731678806367257 … 0.001495696431378364 -0.0047316788745629616;;; -0.0069191754820690475 -0.004731678806659197 … 0.0014956965209409418 -0.004731678836457793; -0.004731678833984086 -0.002991147869411289 … 0.0017333324488070997 -0.0029911479002259073; … ; 0.0014956964830048197 0.0017333324276118094 … 0.002146714857306435 0.0017333323637997236; -0.004731678840347992 -0.002991147883159382 … 0.0017333323942985731 -0.0029911479200255327;;; 0.0006640007584446765 0.001495696573197665 … 0.003561045712835274 0.0014956965702198899; 0.0014956965371086452 0.0017333324812972492 … 0.002146715002626502 0.0017333324694536077; … ; 0.0035610455767991138 0.0021467149232913285 … 1.0228896941789217e-5 0.0021467149021190787; 0.0014956965300568963 0.0017333324546208537 … 0.002146714965824966 0.001733332450819189;;; … ;;; -0.00091706949765487 0.0007544789309015566 … 0.004144157300991494 0.0007544790065413323; 0.000754478949812452 0.002952825242789711 … 0.001507744907449299 0.0029528252977393866; … ; 0.004144157195777432 0.0015077448571751625 … -0.0015401026813997697 0.001507744756315781; 0.0007544790266512556 0.002952825308364598 … 0.0015077447822878943 0.002952825452543244;;; 0.000664000778050189 0.0014956966424785724 … 0.003561045614734958 0.001495696523347502; 0.001495696598936472 0.0017333325641650213 … 0.002146714998648462 0.0017333324848745014; … ; 0.003561045623435863 0.0021467150071400225 … 1.0228853930509783e-5 0.0021467148691129384; 0.001495696540165306 0.0017333325280128606 … 0.00214671486218454 0.0017333323797741411;;; -0.006919175472663779 -0.004731678777874514 … 0.001495696475850971 -0.004731678853600938; -0.0047316788115101364 -0.002991147835323338 … 0.001733332444239003 -0.002991147896065754; … ; 0.0014956965130287162 0.001733332486186681 … 0.0021467148189831496 0.0017333323490428382; -0.004731678832088708 -0.0029911478491162907 … 0.0017333323297369155 -0.0029911479464623403]), 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.715162548354328e-6 - 8.146685532828851e-7im 0.005211842376499838 - 0.00130318149535776im … 0.003642548889642904 - 0.0009089733763449742im -0.005303720725216851 + 0.0013224613261264106im; -0.0051990607834715066 - 0.0181624125061111im 0.0016645682064416944 + 0.017197201065480137im … -0.005381007988802941 - 0.009358928840629523im 0.007826456107790689 + 0.015658656774193418im; … ; -0.0062093390385489795 - 0.02174237343191929im 0.006480660968967204 + 0.01874336599334707im … -0.002299105448427878 - 0.012548929313949016im 0.004365336718217275 + 0.0192715215562985im; 0.0052620624078335365 + 0.01815294852165761im -0.007802033035899988 - 0.015666171692350822im … 0.0003713123078189677 + 0.010610520941739122im -0.0016401286050284599 - 0.017204721750391214im;;; 0.0043801311778836555 - 0.00525616933779526im -0.005617067321324217 + 0.003842098928508721im … -0.0016259972152825987 - 0.00043107301992031526im 0.0005448533395812799 + 0.002303544704352022im; 0.0021531204526842036 + 0.019513568654634113im -0.00042710515670011555 - 0.016528268825496947im … 0.003610392410346641 + 0.008785159058921452im -0.004394459111815099 - 0.015537682206742386im; … ; 0.006312087535259899 + 0.017948206577358326im -0.006192215684786396 - 0.015327760240418873im … 0.0024416846700090434 + 0.010455827841659139im -0.004531709098208487 - 0.015742356010227533im; -0.0073134122336644194 - 0.013349832234298153im 0.007669788980342013 + 0.011580437721330649im … -0.0012428457723660663 - 0.008063136145594438im 0.003702410226423176 + 0.012571032209072816im;;; -0.002961954808392183 + 0.0035775072660681165im 0.0030767441252359047 - 0.0027166602933101154im … 0.00023408793390858855 + 0.0008772623600408637im 0.0009614076028147371 - 0.0021885005718443673im; -0.0008468217408018231 - 0.01243968483162715im 9.409109302140558e-5 + 0.010460495413513331im … -0.0017527276849165852 - 0.006122109423100522im 0.0017545868681274925 + 0.010045902474058196im; … ; -0.004476993305615572 - 0.01106866719341296im 0.004343109404954384 + 0.009573982399969656im … -0.0019932178603529705 - 0.006937525124809697im 0.003388646823230976 + 0.009812286505895558im; 0.0049054682154860115 + 0.0075297751802993656im -0.004759148253013962 - 0.006387797556519628im … 0.0013514905647336504 + 0.004654362424287242im -0.003098635805213527 - 0.0068023912749536255im;;; … ;;; -0.0013720488367085922 + 0.0015975838344232415im 0.0007482824743657319 - 0.0012465564190256846im … -0.001359937690836156 + 0.0010434442262161076im 0.0015796009278792288 - 0.0014541186582917717im; 0.0025215279510881986 + 0.003365253930592109im -0.001848018681206349 - 0.003100953147754361im … 0.002086601923262329 + 0.0020796703209225763im -0.0025694680828321076 - 0.002920828598881372im; … ; 0.0011322088244097974 + 0.007547241875876506im -0.0012808733176611301 - 0.006740715219568973im … 0.0003966103866752564 + 0.0052469161518546425im -0.0007731829099501874 - 0.006867470598154983im; -0.00019514457511755112 - 0.006065743881498713im 0.0005570124167942903 + 0.005248177851199454im … 0.00035123658486876247 - 0.003944702527616522im -0.0001644368474167253 + 0.005428309487085703im;;; 0.0030614805430227597 - 0.003604538666151148im -0.0010594632935037366 + 0.002215136181699706im … 0.002221495294314396 - 0.0014924303735768685im -0.0031747618039681767 + 0.002743277879826644im; -0.004983045464871413 - 0.0075078873224417204im 0.0031814749815973794 + 0.0067792932687945244im … -0.003418403654912287 - 0.0041360434835601145im 0.004841941035585639 + 0.0063647148575628545im; … ; -0.002489593878235056 - 0.013438349031925999im 0.0026774425659955537 + 0.011560042356126489im … -0.0009062710579142301 - 0.008233565362437328im 0.0017229991046490423 + 0.011798340840944183im; 0.0007693170177608922 + 0.012461530084919886im -0.0016718144819021742 - 0.010068970845035316im … -0.0003141694157695452 + 0.006640422287064053im -1.133884206342666e-5 - 0.010483552816547824im;;; -0.004452129109937385 + 0.005275023549732615im -0.00044150631960200717 - 0.002329001284805816im … -0.003402696357501043 + 0.0016873943684161583im 0.00572036426333002 - 0.0038675303121333037im; 0.007328706757820484 + 0.013342469806142749im -0.0037594511631271242 - 0.012554990024329828im … 0.005005078547217977 + 0.0071222060766761176im -0.007726758254392289 - 0.011564423800855633im; … ; 0.004535406726718545 + 0.02006664704383866im -0.00479753715308765 - 0.016990694242417084im … 0.0016400389387383448 + 0.011411681456700646im -0.0031370572285060513 - 0.01740528272787495im; -0.0021379381479255397 - 0.0195208561910358im 0.004337519868275938 + 0.015553668360625859im … 0.00015181576418432213 - 0.009726072621421084im 0.00037019195103632876 + 0.016544241915675134im],)])]), [8.519873695721162e-7 4.6272650621557724e-7 … 2.4177785828805752e-8 4.6272651261304043e-7; 4.627265097206593e-7 2.5014476450577424e-7 … 1.3582999345642097e-8 2.501447681910097e-7; … ; 2.4177785607264172e-8 1.3582999072988776e-8 … 1.5659236143368096e-9 1.3583000192009354e-8; 4.6272651105017374e-7 2.5014476536122833e-7 … 1.3583000216526114e-8 2.501447706268041e-7;;; 4.6272650979205104e-7 2.5014476536742957e-7 … 1.3582999290369693e-8 2.5014476766339326e-7; 2.501447674567168e-7 1.389560856429189e-7 … 1.1279991532366889e-8 1.389560872217757e-7; … ; 1.3582999688694016e-8 1.127999173381636e-8 … 7.937757428712554e-9 1.1279992312016792e-8; 2.5014476798980667e-7 1.3895608637742908e-7 … 1.1279992031994878e-8 1.3895608826659906e-7;;; 2.417778473798572e-8 1.3582998753759937e-8 … 1.5659230233107748e-9 1.3582998783240444e-8; 1.3582999122873979e-8 1.1279991239534498e-8 … 7.937756386971943e-9 1.127999134604712e-8; … ; 1.565923344332878e-9 7.937756959892892e-9 … 3.571856523325892e-8 7.93775710992196e-9; 1.3582999200981132e-8 1.1279991486547341e-8 … 7.937756653088228e-9 1.1279991519891018e-8;;; … ;;; 5.79184141324532e-8 2.28160217914727e-8 … 5.620380382297653e-10 2.281602067489915e-8; 2.2816021507695532e-8 3.4892769181167328e-9 … 1.3458513597422145e-8 3.48927669462957e-9; … ; 5.620381619581988e-10 1.3458514123753461e-8 … 7.728316054289915e-8 1.3458515156549507e-8; 2.281602038365158e-8 3.489276654306092e-9 … 1.3458514886545461e-8 3.4892760684263417e-9;;; 2.4177784432949334e-8 1.3582998036067794e-8 … 1.5659232530757434e-9 1.3582999264078229e-8; 1.3582998481888294e-8 1.127999048426897e-8 … 7.937756413311993e-9 1.1279991203368142e-8; … ; 1.565923234296727e-9 7.93775636015945e-9 … 3.5718566093021765e-8 7.937757343833558e-9; 1.358299909414446e-8 1.1279990817828375e-8 … 7.937757390711055e-9 1.1279992162804502e-8;;; 4.627265086593735e-7 2.5014476313000014e-7 … 1.3582999754226995e-8 2.5014476897903533e-7; 2.501447657066786e-7 1.3895608388216954e-7 … 1.1279991572544118e-8 1.3895608700030304e-7; … ; 1.3582999377797984e-8 1.1279991200910388e-8 … 7.937757701347689e-9 1.1279992445037703e-8; 2.501447673413739e-7 1.3895608462043896e-7 … 1.1279992617156533e-8 1.3895608961972252e-7;;;;], Matrix{ComplexF64}[[0.25033371715220154 - 0.07211102642860902im -0.9012275585325155 - 0.34520791014977303im 5.634425651394212e-7 - 4.6744730740186286e-7im -8.627601978300489e-7 - 2.9843592975951836e-6im; -0.16984614655024263 + 0.048925810376722174im -0.03971163235658334 - 0.015211215142257387im -0.0019236819598030042 - 0.5896453691300418im 0.024874785263649935 + 0.03876595229996352im; … ; 0.05396138989290808 - 0.015544095503460522im 0.014363178753564862 + 0.005501700461053128im -0.0034881164626761646 - 0.009447384568215336im -0.001168192830027646 - 0.0013352134150369259im; -0.0967531415843154 + 0.02787066970382591im -0.02845456931741789 - 0.010899296953642007im 0.019250424931755286 + 0.020820057984214286im 0.007587266484331693 + 0.00945272657527118im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944105667999, [[-0.554110423735072, -0.01067839739852787, 0.14495658413841336, 0.14495754664422988]], 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.6206935819507673e-7, -6.293615758346271e-7, -2.2951208148773962e-7, 4.994025479553151e-7, 7.923311479847341e-7, 5.851071678620778e-7, 6.516970638724782e-7, 1.0189308665680189e-6, 3.638592433445862e-7, -2.201093195264199e-6  …  2.781018166205132e-6, 1.2016954817120018e-6, 2.4456377733713487e-8, 1.0078549131129596e-7, 5.535308574450898e-7, 4.5341734218407676e-7, 4.6319673101851055e-9, -1.0437049701333233e-7, 1.1211098325969032e-7, 1.2550865713050196e-7], resid_history = [0.24939209401325296, 0.0037665477895093076, 0.0002852324846299492, 4.6817332683989165e-6, 2.8544714066296475e-8, 9.274637266863058e-11, 4.2415821540481e-8, 1.1629336716400006e-9], converged = true, n_iter = 8, residual_norm = 1.1629336716400006e-9, maxiter = 100, tol = 1.0e-8, s = 0.9354406777336488, residual = [9.348064546109092e-7, 6.135109896923486e-7, 1.2652270314574828e-7, -1.8186597146946774e-7, -2.1011128783302067e-7, -1.1935992143874029e-7, -1.0580731308322246e-7, -1.3851856138530742e-7, -4.4263732813127316e-8, 2.46487555412934e-7  …  9.333733958295629e-7, 4.130235911135222e-7, 8.967903022130056e-9, 3.995573317837445e-8, 2.3494088734493948e-7, 2.0836514543202642e-7, 2.4587250350998396e-9, -7.220055590611446e-8, 1.2378343305220105e-7, 4.07299244468478e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.371392924161941e-7, -4.126805340286087e-8, 1.1629336716400006e-9, -0.0002044305738735834, 4.686493289380406e-6, -2.8907820851188503e-8, 9.274637266863058e-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.0417345522592815e-6, -2.4698560034960018e-6, -8.094946007088119e-7, 1.350618131407852e-6, 1.7029325017384699e-6, 1.2043005407657634e-6, 1.4149897268949532e-6, 2.1206770629098967e-6, 6.198359454290168e-7, -3.0476893042987784e-6  …  3.3024114275189576e-6, 1.3019291064390632e-6, 3.008616940294739e-8, 1.4372356305105904e-7, 7.685169632007942e-7, 5.152289326531379e-7, 4.4119476284030205e-9, -1.0973497764006995e-7, 1.5532777043850508e-7, 1.6402781169512018e-7], [1.697783894719068e-6, 3.188746908540312e-6, 1.0222149370831558e-6, -1.7194103478543475e-6, -2.177054237538432e-6, -1.5170557353033527e-6, -1.7547855396363213e-6, -2.628779650879854e-6, -7.780408619516801e-7, 3.852333323536785e-6  …  -3.310550213609898e-6, -1.2681147355811215e-6, -2.9137965036773278e-8, -1.402217725587966e-7, -7.403964617491734e-7, -4.6661869186157375e-7, -3.5166615701804555e-9, 7.799954667992295e-8, -9.630458894354173e-8, 6.730522528684276e-8]]), H = [1.1412540236769346 0.014113511617765582 … 0.0 0.0; 0.1267755767760091 1.0349249594951158 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1482738320739017 0.12828848526122866 … 0.0 0.0; 0.0 1.0274261252752503 … 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.9257125192668594
Interacting polarizability:     1.7736548502832914

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.08e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.02e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 4.10e-12
└ * number of operations = 11
Non-interacting polarizability: 1.925712519266769
Interacting polarizability:     1.7736548463630086

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