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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.0001345734947229307
Displaced dipole:  0.01761222236507526
Polarizability :   1.774679585979819

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.0032601876935385516 - 0.000483835589758676im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.04798741276059172 + 0.3056524241345873im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.006192037700326695 + 0.04059563422008327im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.01157491449660487 - 0.0801062622724581im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620694560114834e-7 -2.5020149240236825e-7 … -4.880645220376899e-8 -2.5020155330874163e-7; -6.293616678640967e-7 -4.793036209884051e-7 … -1.207784602709649e-7 -4.793036636838991e-7; … ; 1.362626087673592e-7 1.1356230795651641e-7 … 4.711991439569727e-8 1.1356229815300408e-7; 1.3305472626090228e-7 1.3632640381752839e-7 … 5.20068470461721e-8 1.3632635904136105e-7;;; -2.502015060431426e-7 -1.7404905282260126e-7 … -3.5564956679175806e-8 -1.7404909760396915e-7; -4.793036278145252e-7 -3.6797317675862585e-7 … -1.120289639127228e-7 -3.679732091880301e-7; … ; 1.1356230474711308e-7 1.121109954028269e-7 … 1.098068407897366e-7 1.1211098738197972e-7; 1.3632639169954166e-7 1.2550867091512623e-7 … 5.094908313034436e-8 1.2550863818317627e-7;;; -4.880643023140157e-8 -3.55649378453091e-8 … -1.1412155444149483e-8 -3.5564948489290684e-8; -1.2077844224122078e-7 -1.1202894606944305e-7 … -9.786226482146372e-8 -1.1202895646638309e-7; … ; 4.711991722249784e-8 1.0980685177395132e-7 … 2.493823017661941e-7 1.0980684406388078e-7; 5.2006860653001686e-8 5.0949097935354194e-8 … 4.910627810556613e-8 5.0949088942410555e-8;;; … ;;; 6.475762421867717e-8 4.1103336417597575e-8 … -6.6817062483181314e-9 4.110334984728831e-8; 1.689755686063917e-7 6.653283600879737e-8 … -1.3294454432647222e-7 6.653284303855075e-8; … ; 3.267766829666675e-8 1.60173710864764e-7 … 3.8470363372230533e-7 1.6017369690268222e-7; -9.141989798933609e-8 -3.583033098685499e-8 … 7.037284205963186e-8 -3.583032676932583e-8;;; -4.880644618013588e-8 -3.556495057055079e-8 … -1.1412160361303384e-8 -3.5564960448945026e-8; -1.2077845861066985e-7 -1.1202896354219815e-7 … -9.786227803192841e-8 -1.1202897071493586e-7; … ; 4.711991851724542e-8 1.0980685154707877e-7 … 2.493822902146976e-7 1.0980684014270168e-7; 5.200685308013871e-8 5.094909209160091e-8 … 4.910627285128055e-8 5.0949081914167466e-8;;; -2.502015392046025e-7 -1.740490775153437e-7 … -3.5564962750219426e-8 -1.7404912184260242e-7; -4.793036558980868e-7 -3.6797319934723336e-7 … -1.1202897171164066e-7 -3.679732302006474e-7; … ; 1.1356230141544315e-7 1.121109939974327e-7 … 1.0980683616391116e-7 1.1211098307216906e-7; 1.363263711739822e-7 1.2550865654557703e-7 … 5.0949079134938904e-8 1.2550862289195576e-7;;;;], δHψ = Matrix{ComplexF64}[[0.0018964098863253024 + 0.0002814409583850352im -0.09529503361773352 - 0.0006205013561020329im 0.31373045134017913 + 1.3147765380722911im 0.22111405318020225 - 0.5644642514156277im; 0.03089217166851318 - 0.19588073642600587im -0.1056134095519578 + 1.6234287738209006im 0.06946866825999032 + 0.2905356644359117im 0.04840103332762556 - 0.12505268289990454im; … ; 0.007522911872117552 - 0.04843922945739803im 8.273267639969379e-5 - 0.020686212744996626im 0.00015554783840901238 - 0.027734318277013184im 0.0015472012028530232 + 0.012168230414788945im; -0.010354998776564229 + 0.07162988702877354im -0.000223966367704451 + 0.025590053176405694im 0.08146470942135142 - 0.06862989070354909im 0.039476227925813574 + 0.05558617491185445im]], δVind = [0.001965432061070104 0.0020810005266787275 … 0.003153665810701963 0.002081001034756785; 0.003889131085541421 0.004779269373645941 … 0.010276671653028532 0.0047792699304992445; … ; -0.006150846785833428 -0.008269602480127182 … -0.017302962707465987 -0.008269601425114374; 0.000238682447871684 -0.0003986974613267814 … -0.0036646932144373595 -0.0003986968720888662;;; 0.00208100064761508 0.0022232280460995862 … 0.0034201202993643146 0.0022232286134224776; 0.004779269497933485 0.0057349533404501985 … 0.01090359790095588 0.005734953955231184; … ; -0.008269602190246094 -0.009601007389051123 … -0.012623109102187486 -0.009601006268118495; -0.00039869730784356525 -0.0010730132693383843 … -0.0042306792686548815 -0.0010730126107935589;;; 0.0031536644818859635 0.0034201185968143638 … 0.004832504913841424 0.003420119565897487; 0.010276670213315584 0.010903596195269876 … 0.012272141170267174 0.010903597132442945; … ; -0.017302965210344246 -0.012623111151979846 … -0.009620264351086863 -0.01262311002282398; -0.003664694761182717 -0.004230681136407304 … -0.005491107136003556 -0.004230680077875585;;; … ;;; -0.0023365612907192816 -0.0028154638610394245 … 0.005693318933233568 -0.0028154647350704808; -0.012913095074384998 -0.01784118523548361 … 0.011498341197699252 -0.017841186281821096; … ; -0.025174450353634522 -0.01281088486572172 … -0.008395873642328956 -0.012810883979794371; 0.007597978564991949 0.01022541469607441 … -0.005449731139755414 0.010225413322404883;;; 0.0031536653740749346 0.0034201195521111256 … 0.004832506240474633 0.003420120554785284; 0.010276671154816432 0.010903597055922911 … 0.012272142160513628 0.010903598131182492; … ; -0.017302963217169502 -0.0126231098868765 … -0.009620263603155432 -0.01262310893525825; -0.0036646937421559757 -0.004230680050512416 … -0.005491105997726857 -0.004230679032253152;;; 0.002081000912530815 0.0022232283379791584 … 0.0034201207967628194 0.002223228911662699; 0.004779269806724293 0.005734953668844531 … 0.01090359837858846 0.0057349543002981045; … ; -0.008269601717348755 -0.00960100684021439 … -0.012623108689813391 -0.009601005830444603; -0.0003986970284451112 -0.0010730129553121687 … -0.004230678758044027 -0.001073012298655486;;;;], δρ0 = [-3.602140733396248e-7 -2.488843112185686e-7 … -4.852838279003601e-8 -2.4888431042751013e-7; -5.972251070807275e-7 -4.606978875139971e-7 … -1.1905771133941031e-7 -4.6069788527188105e-7; … ; 1.2735792113748433e-7 1.0836341718036093e-7 … 4.640720682815988e-8 1.0836341765078925e-7; 1.0352008196081539e-7 1.1957962578632091e-7 … 5.0695073428439117e-8 1.1957962599367903e-7;;; -2.4888430907594413e-7 -1.7310675331468147e-7 … -3.535698925007273e-8 -1.7310675305561766e-7; -4.606978829020448e-7 -3.572171217210354e-7 … -1.1091442144041534e-7 -3.57217120821937e-7; … ; 1.0836341660945503e-7 1.084649942225171e-7 … 1.0862239842393086e-7 1.0846499552655686e-7; 1.195796252024441e-7 1.1611148174270317e-7 … 5.020831787429935e-8 1.1611148233047679e-7;;; -4.8528381529944286e-8 -3.535698850632313e-8 … -1.13413227980377e-8 -3.535698881055491e-8; -1.1905770874343257e-7 -1.1091441999705029e-7 … -9.770671397701895e-8 -1.1091442099863478e-7; … ; 4.640720451401974e-8 1.0862239580334809e-7 … 2.484524912713599e-7 1.0862239590528548e-7; 5.069507163485594e-8 5.020831628772185e-8 … 4.926859813545702e-8 5.0208316673948564e-8;;; … ;;; 6.434571235965564e-8 4.083862543357568e-8 … -6.637939142414499e-9 4.0838625656993286e-8; 1.6975034565240977e-7 6.683082718604759e-8 … -1.334048582842473e-7 6.683082907295794e-8; … ; 3.271002452920189e-8 1.6023528826815757e-7 … 3.83901512634124e-7 1.6023528440315322e-7; -9.271365905172915e-8 -3.633619468230296e-8 … 7.126311306615549e-8 -3.6336195637131276e-8;;; -4.8528383959930345e-8 -3.535699116170571e-8 … -1.1341324334909445e-8 -3.535699070246657e-8; -1.1905771498500105e-7 -1.1091442827578899e-7 … -9.77067196138768e-8 -1.1091442609550408e-7; … ; 4.6407209123800864e-8 1.0862240309655948e-7 … 2.4845249564162923e-7 1.0862239952256369e-7; 5.069507408050951e-8 5.020831952612598e-8 … 4.926860045443597e-8 5.0208318743319345e-8;;; -2.4888431227147926e-7 -1.7310675617747363e-7 … -3.535699031152448e-8 -1.7310675543630757e-7; -4.6069788904243715e-7 -3.572171280369565e-7 … -1.1091442467214424e-7 -3.5721712560539765e-7; … ; 1.0836341825019004e-7 1.0846499729521679e-7 … 1.0862239757849786e-7 1.0846499576793867e-7; 1.1957962648119053e-7 1.1611148346881601e-7 … 5.020831841157968e-8 1.1611148318854115e-7;;;;], δeigenvalues = [[0.0004950054073907791, 0.0967986350167122, 0.031936916568209725, 0.04909189304957477]], δ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.15475965281456855 -0.15461620268951634 … -0.15418753099043525 -0.15461620269032628; -0.15461620269004994 -0.15447229888889774 … -0.15404229466760713 -0.1544722988897093; … ; -0.15418753098937626 -0.15404229466600505 … -0.1536083373381833 -0.154042294666827; -0.15461620268979334 -0.15447229888863942 … -0.1540422946673498 -0.1544722988894529;;; -0.1546162026903181 -0.15447229888916636 … -0.154042294667879 -0.15447229888997835; -0.15447229888970151 -0.15432794308056375 … -0.1538966028206874 -0.15432794308137732; … ; -0.15404229466681668 -0.15389660281908035 … -0.15346128700228698 -0.15389660281990447; -0.15447229888944372 -0.15432794308030426 … -0.15389660282042902 -0.15432794308111986;;; -0.1541875309904267 -0.1540422946670599 … -0.15360833733924317 -0.15404229466787908; -0.15404229466759997 -0.15389660281986772 … -0.15346128700307576 -0.1538966028206883; … ; -0.15360833733817097 -0.15346128700145364 … -0.15302191312736352 -0.15346128700228523; -0.15404229466733946 -0.15389660281960538 … -0.15346128700281467 -0.15389660282042822;;; … ;;; -0.1534788058443277 -0.15333139195376252 … -0.152890924918001 -0.15333139195459572; -0.15333139195431122 -0.1531835057317809 … -0.15274162876985806 -0.15318350573261574; … ; -0.15289092491691264 -0.15274162876821068 … -0.1522955395808587 -0.1527416287690566; -0.15333139195404785 -0.15318350573151576 … -0.15274162876959368 -0.15318350573235265;;; -0.15418753098886703 -0.1540422946654949 … -0.15360833733766663 -0.1540422946663156; -0.15404229466603545 -0.1538966028182977 … -0.15346128700149417 -0.15389660281911985; … ; -0.15360833733659446 -0.15346128699987166 … -0.15302191312576927 -0.15346128700070463; -0.15404229466577585 -0.1538966028180364 … -0.15346128700123385 -0.1538966028188607;;; -0.1546162026895331 -0.1544722988883788 … -0.154042294667086 -0.15447229888919167; -0.1544722988889142 -0.1543279430797738 … -0.1538966028198919 -0.15432794308058823; … ; -0.15404229466602368 -0.15389660281828457 … -0.15346128700148523 -0.15389660281910955; -0.15447229888865693 -0.15432794307951486 … -0.15389660281963394 -0.15432794308033124]), 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.015163898086781414 -0.012550129565757923 … -0.004966868865329068 -0.0125501295834188; -0.012550129576290585 -0.010362595560912586 … -0.004135159689693553 -0.0103625955584894; … ; -0.004966868866976378 -0.004135159716731425 … -0.0020697813286110313 -0.00413515973467914; -0.012550129576167758 -0.010362595551166417 … -0.004135159739407669 -0.010362595574279458;;; -0.012550129550970734 -0.010362595519887724 … -0.004135159691693056 -0.010362595547262765; -0.010362595534238957 -0.008622039500352579 … -0.003897512082927491 -0.008622039509918645; … ; -0.0041351597129842405 -0.003897512126480954 … -0.0034841077938407967 -0.0038975121610290405; -0.010362595541732631 -0.00862203950234922 … -0.003897512150311783 -0.008622039534824013;;; -0.004966868789548189 -0.0041351596218669 … -0.0020697812237959494 -0.004135159657549155; -0.004135159638978164 -0.0038975120445094934 … -0.0034841076831167877 -0.003897512075210454; … ; -0.002069781265488538 -0.0034841077404388398 … -0.005620580583362458 -0.0034841077453487985; -0.0041351596537649036 -0.003897512074779034 … -0.0034841077284884533 -0.003897512096199692;;; … ;;; -0.006547897457205231 -0.004876339438447022 … -0.0014866455326756262 -0.004876339466839021; -0.004876339469796313 -0.002677987801262429 … -0.004123055038015199 -0.0026779878536109685; … ; -0.0014866455540631244 -0.004123055115924706 … -0.007170895841960516 -0.004123055057561409; -0.004876339455852266 -0.0026779877797620524 … -0.00412305507395133 -0.0026779878276621933;;; -0.0049668688964105276 -0.004135159777152184 … -0.002069781378324723 -0.004135159752969389; -0.004135159744699751 -0.003897512200547423 … -0.0034841077851195572 -0.003897512158990993; … ; -0.0020697813892706623 -0.003484107870139032 … -0.005620580618828328 -0.0034841078018014505; -0.004135159761386072 -0.0038975122167437004 … -0.0034841078226870085 -0.00389751218316168;;; -0.012550129591314657 -0.010362595577969941 … -0.004135159742682374 -0.010362595580069847; -0.01036259557547994 -0.008622039561090918 … -0.003897512138891098 -0.0086220395449006; … ; -0.004135159738908513 -0.0038975121840390146 … -0.003484107768870857 -0.003897512155480201; -0.01036259557658952 -0.00862203955617269 … -0.0038975121727797034 -0.008622039557025445])], 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.009532892832495355 -0.0069191756814663725 … 0.0006640005069140122 -0.0069191756999371565; -0.006919175692532496 -0.004731679042122506 … 0.0014956963086790991 -0.0047316790405108645; … ; 0.0006640005063257718 0.0014956962832432235 … 0.003561045285769493 0.0014956962644735554; -0.006919175692153013 -0.004731679032117988 … 0.0014956962592223332 -0.004731679056044517;;; -0.006919175667480959 -0.004731679001366429 … 0.0014956963064075915 -0.004731679029553432; -0.004731679016252651 -0.0029911480782586364 … 0.0017333322399174502 -0.0029911480886382746; … ; 0.0014956962861788076 0.0017333321979709796 … 0.002146714654398601 0.0017333321625987745; -0.004731679023488504 -0.0029911480799957625 … 0.001733332172791618 -0.002991148113286154;;; 0.0006640005827034121 0.0014956963770528424 … 0.0035610453895247004 0.0014956963405514095; 0.001495696359401594 0.0017333322791551532 … 0.002146714764333852 0.001733332247633627; … ; 0.0035610453489043374 0.0021467147086338635 … 1.0228708229390211e-5 0.00214671470289232; 0.0014956963448754515 0.001733332249147958 … 0.0021467147192233385 0.0017333322269044585;;; … ;;; -0.000917069375963241 0.0007544791875988138 … 0.004144156841969568 0.0007544791583736204; 0.0007544791557008781 0.002952825670529883 … 0.0015077445812691666 0.0029528256173465112; … ; 0.004144156821670449 0.0015077445050069808 … -0.0015401029002872256 0.001507744562524371; 0.0007544791699083259 0.002952825692295464 … 0.001507744545597491 0.002952825643558437;;; 0.0006640004774007427 0.001495696223332529 … 0.0035610452365724713 0.001495696246694619; 0.0014956962552445333 0.0017333321246872455 … 0.0021467146639126784 0.0017333321654215278; … ; 0.003561045226698702 0.0021467145805156556 … 1.022867435774514e-5 0.0021467146480202647; 0.001495696238817866 0.0017333321087522591 … 0.002146714626605602 0.0017333321415099944;;; -0.0069191757070398985 -0.004731679058661082 … 0.0014956962562112782 -0.004731679061573837; -0.004731679056706319 -0.002991148138206997 … 0.0017333321847493735 -0.0029911481228311106; … ; 0.001495696261047512 0.001733332141208699 … 0.0021467146801702605 0.0017333321689425338; -0.0047316790575586045 -0.002991148133029836 … 0.001733332151118784 -0.0029911481346989663]), 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.4704006982476143e-6 + 8.384347134425424e-7im 0.008891708909667582 - 0.0040556144562367774im … 0.00614331853688317 - 0.0027594077292231204im -0.008969123992741482 + 0.004050243657331141im; 0.003640437547362466 - 0.024317250093041178im -0.008222683220813282 + 0.022588382269478402im … -0.002529869288518604 - 0.009516808073626376im 0.0022423814606623567 + 0.01783877682697934im; … ; 0.003941292400832486 - 0.026200561337552133im -0.0015470013064881924 + 0.021368402687791277im … 0.0041129583316256105 - 0.014775547069141812im -0.005139517469091412 + 0.022998927007717385im; -0.0035792170089255634 + 0.024331273015009536im -0.002222667445674261 - 0.017842346184558933im … -0.005917104458840073 + 0.013305847895097568im 0.008242887423408302 - 0.022592057118779117im;;; 0.0004856927379664961 - 0.0025443703365198144im -0.005489466831004798 + 0.003867533389234344im … -0.004127606705951728 + 0.0013823983029828445im 0.004975812207604471 - 0.0008819795022966339im; -0.0032974572594897844 + 0.021702383363143438im 0.0060002058648073045 - 0.019081015955293476im … 0.00166310449050219 + 0.00854947132373719im -0.0007374238641286604 - 0.016023061826098566im; … ; -0.0032376808109144737 + 0.021672969779035902im 0.001345742228909383 - 0.017756699385340524im … -0.0034687368626499142 + 0.01257951758593798im 0.004165843797393235 - 0.019036622972385422im; 0.002702891090308306 - 0.018728376038940613im 0.0010753951793023037 + 0.01410423470291506im … 0.0045011062583906305 - 0.01057301441501212im -0.005662574804932362 + 0.017162204032366884im;;; -0.0003036433202829175 + 0.0017386726336724111im 0.001989946870026749 - 0.002010213702129865im … 0.0019375087333373617 - 0.0003357519514286426im -0.0016023612471982034 - 0.00037979084736263196im; 0.001897095227317647 - 0.012604918454040792im -0.002968379641109203 + 0.011093292579045896im … -0.0007296127290668658 - 0.005777001421612439im -0.00014846640583900868 + 0.009813409886988946im; … ; 0.0019855903990929945 - 0.0132921483890875im -0.0009229015297442892 + 0.011198816796449068im … 0.0023101216418524653 - 0.00850613859529405im -0.002543913551984054 + 0.011934509452237218im; -0.001490035725229114 + 0.010217698456458467im -0.00013051878409667402 - 0.008029165507478021im … -0.002527676568326672 + 0.006338816635314241im 0.002689501542350912 - 0.009309064634954777im;;; … ;;; -0.00020459988277800892 + 0.0007574666949689566im -0.000527730325236769 - 0.000318624629321744im … -0.001108276113173471 + 0.0008677012038829557im 0.0008841877433728106 - 0.0009594728471978554im; -0.0006088387214375742 + 0.004660331841556745im 0.0011499059561497188 - 0.004412124132866043im … 0.0004955091380758332 + 0.002597268625700401im -7.53803789265489e-5 - 0.0038560096911191832im; … ; -0.001235869722019051 + 0.007884058400050414im 0.0006856866318310493 - 0.006912046648315118im … -0.0015080508687094446 + 0.00567715479027993im 0.0015479753640398482 - 0.007303462286596903im; 0.0009640726368252014 - 0.005938507614748685im -0.00024026119946241815 + 0.004955304831791699im … 0.0014965638353545913 - 0.004148137522745118im -0.001465576629264362 + 0.005511510746141519im;;; 0.00041839296113100406 - 0.0017179575894408486im 0.0014949798326686597 + 0.0003640404457004107im … 0.0021748291146135 - 0.0014770116274792913im -0.002097345895281563 + 0.001994517976456516im; 0.0014022092563051383 - 0.010230727899854877im -0.0025993547836584275 + 0.009320850726450041im … -0.000928417590683632 - 0.004820452482522714im 0.0002205833693326525 + 0.008040925690352918im; … ; 0.002222886001185799 - 0.01443339852065732im -0.0011217084330299102 + 0.012155334446257273im … 0.0024398126632215614 - 0.009130397130072722im -0.0027427698733467434 + 0.012891063979795843im; -0.0019850053718307396 + 0.012591970632903121im 0.00023856516505910337 - 0.009801647711401934im … -0.002726532173865947 + 0.007295341842729284im 0.003058620218909 - 0.01108158709180447im;;; -0.0005750323237881607 + 0.0025225558208675546im -0.004869204472701868 + 0.0008973282672694037im … -0.004339742117507208 + 0.0024027248559874963im 0.005596003363307541 - 0.003852314680156407im; -0.0026774524991011036 + 0.018732349979217792im 0.005601633243232021 - 0.01716818967278697im … 0.00182956608028826 + 0.007748531141773324im -0.0011360150029701504 - 0.01411017895578074im; … ; -0.0034497820653563295 + 0.022693244919075303im 0.0015122242107875347 - 0.018557619196421012im … -0.003564393409277766 + 0.013039941082943918im 0.00433235843160523 - 0.019837565502752265im; 0.00332305727427596 - 0.021698603281979152im 0.0006767114692474885 + 0.016017150318461338im … 0.004667622406277629 - 0.011373957776307472im -0.006061276413535384 + 0.019075167183832543im],)])]), [8.519874109602262e-7 4.6272653559500543e-7 … 2.4177788854197226e-8 4.6272653769113156e-7; 4.6272653684510536e-7 2.5014478495751636e-7 … 1.3583001653395607e-8 2.5014478477009925e-7; … ; 2.417778887947748e-8 1.3583001932335574e-8 … 1.5659240663470964e-9 1.3583002117491099e-8; 4.627265368305343e-7 2.501447842037076e-7 … 1.3583002166276997e-8 2.5014478599136237e-7;;; 4.6272653383995e-7 2.501447817845095e-7 … 1.3583001674021272e-8 2.5014478390179106e-7; 2.501447828944894e-7 1.389560972900138e-7 … 1.1279993580071466e-8 1.389560977818486e-7; … ; 1.3583001893677808e-8 1.1279993975520361e-8 … 7.937758991714879e-9 1.1279994289205692e-8; 2.5014478347406986e-7 1.3895609739266912e-7 … 1.1279994191897116e-8 1.3895609906233086e-7;;; 2.4177787691429214e-8 1.3583000953653754e-8 … 1.565923820216187e-9 1.3583001321774496e-8; 1.3583001130182984e-8 1.1279993231254656e-8 … 7.937758201955204e-9 1.1279993510006965e-8; … ; 1.5659239181213303e-9 7.93775861081484e-9 … 3.571856933978437e-8 7.937758645838586e-9; 1.3583001282732986e-8 1.1279993506087249e-8 … 7.937758525577727e-9 1.1279993700578617e-8;;; … ;;; 5.791841125306325e-8 2.2816018292919868e-8 … 5.620385957305773e-10 2.281601871150298e-8; 2.2816018755102817e-8 3.4892752606783197e-9 … 1.3458517130674676e-8 3.4892754732447393e-9; … ; 5.620386207102303e-10 1.3458517929360121e-8 … 7.728316865286693e-8 1.3458517331053235e-8; 2.2816018549525975e-8 3.4892751733724534e-9 … 1.3458517499074577e-8 3.4892753678750433e-9;;; 2.4177789331104695e-8 1.3583002555668169e-8 … 1.565924183086615e-9 1.358300230618781e-8; 1.3583002220871339e-8 1.127999464801395e-8 … 7.937758929511209e-9 1.1279994270701252e-8; … ; 1.5659242087904415e-9 7.93775953592542e-9 … 3.571857005202501e-8 7.93775904849922e-9; 1.358300239301673e-8 1.1279994795070095e-8 … 7.937759197467941e-9 1.1279994490160535e-8;;; 4.6272653862827176e-7 2.5014478627679145e-7 … 1.3583002200061123e-8 2.5014478643920414e-7; 2.5014478608420935e-7 1.3895610041282714e-7 … 1.1279994088200833e-8 1.3895609958041027e-7; … ; 1.358300216112417e-8 1.1279994498127035e-8 … 7.93775881361539e-9 1.1279994238820624e-8; 2.5014478617003004e-7 1.3895610015996017e-7 … 1.1279994395896257e-8 1.3895610020379516e-7;;;;], Matrix{ComplexF64}[[0.25769057894241004 + 0.038243118094542526im -0.9650596494429409 - 0.006283830417132192im 1.4153179017837785e-6 - 1.2448737508262171e-6im 1.181712565019374e-5 + 7.3060084739670335e-6im; -0.17483762142685813 - 0.025947148813431183im -0.04252434378755687 - 0.0002768886848837247im 0.4158959819866187 - 0.09941978162079151im -0.17830785398072507 - 0.0703851505123255im; … ; 0.055547218738583135 + 0.008243603064630045im 0.015380492445980424 + 0.00010014736130319048im 0.00383427282794567 - 0.005855934131293783im -0.0002327097800676821 - 0.0030616347175833206im; -0.0995965436762178 - 0.014780836754315841im -0.030469947676683186 - 0.00019840014226432im 0.0010341685268693607 + 0.027202156617930493im -0.008248696729166801 + 0.01327725416381594im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239440942919525, [[-0.5541104206301244, -0.01067839822826614, 0.1449565724267698, 0.14495681751543554]], 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.620694560114834e-7, -6.293616678640967e-7, -2.2951211937386145e-7, 4.994025682673312e-7, 7.923311878395275e-7, 5.851071883408568e-7, 6.516970495654654e-7, 1.018930801509401e-6, 3.63859089404586e-7, -2.2010934655835484e-6  …  2.781018248807814e-6, 1.2016955983383577e-6, 2.4456485006131473e-8, 1.0078555030273109e-7, 5.535308672523145e-7, 4.534173333421849e-7, 4.631965990323311e-9, -1.0437048800205498e-7, 1.1211098307216906e-7, 1.2550862289195576e-7], resid_history = [0.2493920967672644, 0.003766548805124366, 0.00028527170407797336, 4.680718305681144e-6, 4.279271223305402e-8, 6.231962654551307e-11, 4.2425196509141486e-8, 1.0966616866608587e-9], converged = true, n_iter = 8, residual_norm = 1.0966616866608587e-9, maxiter = 100, tol = 1.0e-8, s = 0.9356264600267401, residual = [9.008134500303847e-7, 6.34042941776252e-7, 1.394374702663698e-7, -2.1001180256876008e-7, -2.530992193104225e-7, -1.5278963032161168e-7, -1.4647503991861862e-7, -2.0333912842264862e-7, -6.577701109648806e-8, 3.6724210971731533e-7  …  7.258324063458872e-7, 3.251948787129497e-7, 7.10222309929461e-9, 3.186069455219409e-8, 1.9054729661141563e-7, 1.734759235634811e-7, 2.1027553642178534e-9, -6.31958113856203e-8, 1.110849380945202e-7, 3.788116488246548e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.116737528781672e-7, -4.1511227459321406e-8, 1.0966616866608587e-9, -0.00020445835354098208, 4.683805472453688e-6, -4.321987529528894e-8, 6.231962654551307e-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}}([[-3.5080434132295774e-7, -2.1163678583619976e-6, -7.394696717261253e-7, 1.2091833931285502e-6, 1.474003468372022e-6, 1.0526007321846102e-6, 1.2730137462491734e-6, 1.8997294044022262e-6, 5.240269383390109e-7, -2.3818197440568057e-6  …  3.3935624746519325e-6, 1.3358562390832629e-6, 3.1860783622680055e-8, 1.5720113548345988e-7, 8.593747447968109e-7, 5.922309287808352e-7, 5.391440265578436e-9, -1.4796954385184193e-7, 2.4035720666749515e-7, 4.79375910670607e-7], [6.762771199981194e-7, 1.8499718848024097e-6, 6.192828893811782e-7, -1.0225766595925744e-6, -1.2551680304079674e-6, -8.760464710340522e-7, -1.034001576583425e-6, -1.5415330586563523e-6, -4.342071075198545e-7, 2.0078148746582763e-6  …  -2.0668924576529194e-6, -7.782582839010102e-7, -1.856535612677035e-8, -9.329485649303176e-8, -5.035742711072756e-7, -3.222099378442676e-7, -2.575415442404478e-9, 6.677605127677403e-8, -1.06991720053955e-7, -1.4132911199485922e-7]]), H = [1.0964555824073257 0.01454657569308105 … 0.0 0.0; 0.13742697018598388 1.031148783163715 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1050343959925804 0.1426718731801371 … 0.0 0.0; 0.0 1.0216758878793406 … 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.9257125359484086
Interacting polarizability:     1.7736548654078528

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 = 9.79e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.04e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 2.46e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125359482932
Interacting polarizability:     1.7736548584994447

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