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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457365822168132
Displaced dipole:  0.01761222167359714
Polarizability :   1.774679533181882

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.0014417020541519184 + 0.0029638510702778356im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.2770803053770981 + 0.13766880482401492im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.03685535491278676 + 0.018111588811510725im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.07291920476780081 - 0.035125225373962024im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206946466389894e-7 -2.502015495062854e-7 … -4.8806433145687373e-8 -2.5020150951636163e-7; -6.293617201899408e-7 -4.79303699129044e-7 … -1.2077846194660482e-7 -4.793036673479019e-7; … ; 1.3626263469499964e-7 1.1356232358998504e-7 … 4.711993741440094e-8 1.1356232846783813e-7; 1.3305476587264856e-7 1.363263988165075e-7 … 5.2006874586545226e-8 1.3632642507588833e-7;;; -2.502014687941553e-7 -1.7404906268146645e-7 … -3.5564935517887826e-8 -1.740490339769943e-7; -4.793036381669765e-7 -3.679732148437475e-7 … -1.1202895969022217e-7 -3.679731921468519e-7; … ; 1.135623353084028e-7 1.1211101574672019e-7 … 1.0980688162755369e-7 1.1211102497432037e-7; 1.3632645375440264e-7 1.2550869192810544e-7 … 5.0949115468137536e-8 1.2550871231930768e-7;;; -4.880641472375375e-8 -3.556493530033505e-8 … -1.1412147021875857e-8 -3.5564928855766874e-8; -1.20778451238881e-7 -1.1202896214351606e-7 … -9.786225511184987e-8 -1.1202895607291024e-7; … ; 4.7119939792511876e-8 1.0980687360396023e-7 … 2.4938238154129014e-7 1.0980688523954776e-7; 5.200688684688516e-8 5.094911418761449e-8 … 4.9106310810673935e-8 5.094912078138653e-8;;; … ;;; 6.475766601969223e-8 4.110337087419082e-8 … -6.681710856134469e-9 4.110336464201868e-8; 1.6897558849661238e-7 6.65328357896669e-8 … -1.32944567572812e-7 6.653284009284887e-8; … ; 3.267767675663589e-8 1.601737023871901e-7 … 3.847036617058033e-7 1.601737057295102e-7; -9.141987105676797e-8 -3.583031102941208e-8 … 7.037284936835255e-8 -3.583031716409826e-8;;; -4.8806471528906505e-8 -3.556498017084166e-8 … -1.1412163623977932e-8 -3.55649711603082e-8; -1.207784972694775e-7 -1.1202901031879663e-7 … -9.786229561124532e-8 -1.1202899698797895e-7; … ; 4.711993022270418e-8 1.0980685698716675e-7 … 2.4938232675920653e-7 1.0980685633330171e-7; 5.2006850325572865e-8 5.094908474985267e-8 … 4.910628544670289e-8 5.0949086759293333e-8;;; -2.502015911657956e-7 -1.7404915372581242e-7 … -3.5564957123949684e-8 -1.7404912328692714e-7; -4.793037311265183e-7 -3.6797328714844966e-7 … -1.1202898189687782e-7 -3.6797326056234364e-7; … ; 1.1356231692864333e-7 1.1211100433828262e-7 … 1.0980686588581776e-7 1.1211100508004006e-7; 1.363263695804909e-7 1.2550863150881047e-7 … 5.094909837603628e-8 1.255086491960004e-7;;;;], δHtotψ = Matrix{ComplexF64}[[0.000838619874955755 - 0.001724034775611319im -0.015280385930828862 + 0.09406397326841119im 1.1419430006353972 - 1.477738668956793im 0.11967559500691208 - 0.17865273852724425im; -0.17752821259112767 - 0.08835905618065189im 1.5878288459299215 + 0.35422409611770567im 0.25227509483325505 - 0.32697987136310747im 0.0273241104054483 - 0.03949673201083111im; … ; -0.04393584169843898 - 0.021739249143916062im -0.020426976449156177 - 0.0032656588644510093im -0.01932885295722618 + 0.02332014131926558im 0.005250841956419154 + 0.012488524661252172im; 0.06520185653158028 + 0.031413128272275236im 0.025250671207166054 + 0.004159965009289708im -0.02653677176246382 + 0.0018899047516692385im 0.043915060894825046 + 0.13078433833769235im]], δVind = [0.0019654320359982222 0.002081000897975547 … 0.0031536642100935615 0.002081000584723232; 0.0038891312759765425 0.004779269990777073 … 0.01027666990068616 0.004779269647750099; … ; -0.006150847421260135 -0.0082696021479234 … -0.017302964221928736 -0.008269602807746463; 0.00023868216873182174 -0.00039869727732231804 … -0.0036646951449654072 -0.00039869763910052873;;; 0.002081000233312073 0.0022232280301415176 … 0.0034201177627262514 0.0022232276732723913; 0.004779269218914012 0.005734953476291238 … 0.010903595227946119 0.00573495308403043; … ; -0.008269603374833031 -0.009601007710653524 … -0.012623111754962876 -0.009601008200922087; -0.0003986980058139915 -0.0010730134936634804 … -0.0042306819409996255 -0.0010730138938079981;;; 0.00315366299066147 0.003420117700578997 … 0.004832499871755047 0.003420117085719033; 0.0102766683218853 0.010903594998793526 … 0.012272137428292397 0.010903594366752575; … ; -0.01730296627651172 -0.01262311213805441 … -0.009620267147870723 -0.012623112456624395; -0.0036646964406413613 -0.00423068201221677 … -0.005491110666580564 -0.004230682613512416;;; … ;;; -0.002336562893805392 -0.002815466539820155 … 0.005693318655559607 -0.0028154659459599684; -0.012913097632727764 -0.017841190912242507 … 0.011498341498017034 -0.017841189456849017; … ; -0.025174444640780367 -0.012810882308269576 … -0.008395873608726524 -0.012810883008513647; 0.0075979771825799266 0.010225412161204548 … -0.0054497307299647794 0.01022541265982526;;; 0.0031536665128180275 0.003420121519753454 … 0.004832505299468562 0.003420121005857212; 0.010276672332083608 0.010903598874769709 … 0.012272141531173385 0.010903598548393758; … ; -0.017302960414795655 -0.012623107587538616 … -0.009620264100335979 -0.012623108542049564; -0.0036646925639915984 -0.004230677862681215 … -0.005491106303211289 -0.0042306786452850865;;; 0.002081001253467061 0.0022232291644167564 … 0.0034201197135882307 0.0022232288228859697; 0.004779270410535266 0.005734954778363076 … 0.010903597248778398 0.005734954423249297; … ; -0.008269601560324169 -0.009601005622620834 … -0.012623109856335211 -0.009601006455639025; -0.0003986969049347936 -0.0010730122532101175 … -0.004230679951123509 -0.0010730126677196044;;;;], δρ0 = [-3.60214077206049e-7 -2.4888431656478524e-7 … -4.8528384752620365e-8 -2.4888431149109907e-7; -5.972251186829224e-7 -4.606979017000101e-7 … -1.1905772317704206e-7 -4.606978955301798e-7; … ; 1.273579315079702e-7 1.0836342960110463e-7 … 4.640722085953487e-8 1.0836342866435517e-7; 1.0352008789177038e-7 1.19579632388056e-7 … 5.0695079838499647e-8 1.1957963323710152e-7;;; -2.4888431115478083e-7 -1.7310675671593554e-7 … -3.535699098065992e-8 -1.7310675352343692e-7; -4.606978958135185e-7 -3.572171363571786e-7 … -1.1091443320491518e-7 -3.572171326103614e-7; … ; 1.0836342868430503e-7 1.0846500556576653e-7 … 1.086224140859867e-7 1.0846500934233595e-7; 1.195796323883088e-7 1.1611148888720248e-7 … 5.02083246690304e-8 1.1611149034815306e-7;;; -4.852838454528426e-8 -3.535699169838105e-8 … -1.1341324489955684e-8 -3.535699111611201e-8; -1.1905772602505539e-7 -1.109144364754402e-7 … -9.770672342355154e-8 -1.109144357351909e-7; … ; 4.640721869186477e-8 1.0862240546426951e-7 … 2.4845250751029174e-7 1.0862241258182867e-7; 5.0695078533427174e-8 5.02083225722667e-8 … 4.92686045193596e-8 5.0208323825435555e-8;;; … ;;; 6.434570878463043e-8 4.083862112634817e-8 … -6.637942523639627e-9 4.083862239411733e-8; 1.6975033098504098e-7 6.683080943174024e-8 … -1.3340487140786474e-7 6.683081670425423e-8; … ; 3.2710035504763524e-8 1.6023530011651188e-7 … 3.8390153309067466e-7 1.6023529769080142e-7; -9.271365442344131e-8 -3.633618745909816e-8 … 7.126312206746848e-8 -3.6336190698304705e-8;;; -4.852838911731943e-8 -3.53569975668155e-8 … -1.1341327961294044e-8 -3.535699443678523e-8; -1.1905773301629824e-7 -1.109144491508691e-7 … -9.770673442649131e-8 -1.1091444121880515e-7; … ; 4.640722218266781e-8 1.08622418176336e-7 … 2.4845251491825015e-7 1.0862241312764145e-7; 5.0695080910595144e-8 5.020832843234343e-8 … 4.9268608620152456e-8 5.020832515771898e-8;;; -2.4888431781113715e-7 -1.7310676315018347e-7 … -3.5356993113689574e-8 -1.7310675823379476e-7; -4.606979041226378e-7 -3.5721714578099867e-7 … -1.1091443768150873e-7 -3.5721713825073564e-7; … ; 1.0836342974674209e-7 1.0846501165986576e-7 … 1.0862241310268657e-7 1.0846500722031268e-7; 1.1957963270043386e-7 1.1611149102253658e-7 … 5.0208325247322976e-8 1.1611148995014051e-7;;;;], δeigenvalues = [[0.0004950052280080835, 0.09679857512815478, 0.01337792141174204, 0.052123284280652]], δ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.1547596528204824 -0.15461620269660387 … -0.15418753099400023 -0.15461620269506252; -0.15461620269583698 -0.15447229889586023 … -0.15404229467103683 -0.15447229889431482; … ; -0.15418753099556012 -0.15404229467337818 … -0.15360833734198 -0.15404229467181418; -0.15461620269584125 -0.15447229889586542 … -0.15404229467103792 -0.15447229889431766;;; -0.1546162026956405 -0.1544722988956664 … -0.15404229467083022 -0.15447229889411407; -0.15447229889489358 -0.1543279430869351 … -0.15389660282350068 -0.15432794308537925; … ; -0.15404229467240155 -0.1538966028258593 … -0.1534612870054624 -0.15389660282428358; -0.15447229889489855 -0.15432794308694148 … -0.15389660282350168 -0.15432794308538234;;; -0.1541875309951783 -0.15404229467299937 … -0.15360833734158508 -0.15404229467142697; -0.15404229467221617 -0.15389660282567347 … -0.15346128700527564 -0.1538966028240979; … ; -0.15360833734317736 -0.15346128700766618 … -0.15302191312992294 -0.1534612870060691; -0.15404229467222197 -0.15389660282568116 … -0.1534612870052769 -0.15389660282410153;;; … ;;; -0.15347880585206758 -0.1533313919627034 … -0.15289092492335724 -0.15333139196113588; -0.15333139196192364 -0.1531835057405971 … -0.1527416287750791 -0.15318350573902487; … ; -0.15289092492494227 -0.1527416287774585 … -0.1522955395864735 -0.15274162877587022; -0.15333139196192735 -0.153183505740601 … -0.1527416287750815 -0.15318350573902811;;; -0.15418753099600094 -0.1540422946738134 … -0.1536083373424424 -0.15404229467226407; -0.15404229467304284 -0.15389660282649226 … -0.15346128700613573 -0.1538966028249384; … ; -0.15360833734400914 -0.15346128700848757 … -0.15302191313079305 -0.1534612870069171; -0.15404229467304636 -0.15389660282649606 … -0.15346128700613743 -0.1538966028249411;;; -0.15461620269605428 -0.15447229889607567 … -0.1540422946712619 -0.15447229889453537; -0.15447229889530953 -0.15432794308734699 … -0.15389660282393372 -0.15432794308580233; … ; -0.1540422946728202 -0.1538966028262727 … -0.15346128700590062 -0.15389660282471054; -0.1544722988953133 -0.15432794308735132 … -0.153896602823935 -0.154327943085805]), 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.01516389835838712 -0.012550129880773113 … -0.004966869127477155 -0.012550129817161818; -0.012550129853059155 -0.010362595865705704 … -0.0041351600044780635 -0.010362595814506619; … ; -0.004966869138454531 -0.004135160037933102 … -0.002069781732168968 -0.004135160005446075; -0.012550129829063514 -0.010362595853334855 … -0.004135160019251783 -0.01036259579525912;;; -0.012550129837509446 -0.010362595833816582 … -0.004135159964476365 -0.010362595803774056; -0.010362595852743762 -0.008622039836399507 … -0.0038975124078381145 -0.008622039809448687; … ; -0.004135160022268412 -0.0038975124073793886 … -0.0034841081176927103 -0.00389751248442142; -0.010362595807838763 -0.008622039789536157 … -0.0038975124362299906 -0.008622039784262719;;; -0.0049668691400836535 -0.004135159987072918 … -0.0020697815174015484 -0.004135159981472081; -0.0041351600927360226 -0.0038975124809594863 … -0.0034841079655474347 -0.003897512475556978; … ; -0.0020697816815244885 -0.003484107958933384 … -0.0056205808278941595 -0.00348410809720775; -0.004135159979042714 -0.0038975123680062597 … -0.003484107976568893 -0.0038975124083207306;;; … ;;; -0.0065478972997986915 -0.00487633915501105 … -0.0014866460554919239 -0.004876339293183909; -0.004876339230208703 -0.0026779873610448457 … -0.00412305532930835 -0.0026779875537952125; … ; -0.0014866458814668916 -0.004123055327157621 … -0.007170896089246638 -0.004123055281471987; -0.004876339296671969 -0.0026779874419722806 … -0.0041230553853456664 -0.0026779875859363183;;; -0.004966869271684463 -0.004135160283259626 … -0.0020697818481256814 -0.004135160083415682; -0.004135160167084422 -0.003897512690909148 … -0.003484108152045893 -0.003897512523560605; … ; -0.0020697817683256073 -0.0034841081924740632 … -0.005620580893552987 -0.00348410808232026; -0.0041351601056665846 -0.003897512659267607 … -0.003484108165954971 -0.0038975124785430885;;; -0.012550129870509177 -0.010362595925681832 … -0.00413516004500482 -0.010362595817513788; -0.010362595861606427 -0.008622039893985519 … -0.00389751247032293 -0.008622039805627184; … ; -0.004135160027705914 -0.0038975125288706943 … -0.0034841080850410494 -0.0038975124189095647; -0.010362595836776338 -0.00862203989113484 … -0.0038975124703139437 -0.00862203977908474])], 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.009532893110014914 -0.006919176003569087 … 0.0006640002412009434 -0.006919175938416414; -0.006919175975088104 -0.00473167935387811 … 0.001495695990464888 -0.004731679301133622; … ; 0.0006640002286637595 0.0014956959546684172 … 0.0035610448784148432 0.0014956959887194431; -0.006919175951096682 -0.004731679341512423 … 0.0014956959756901138 -0.004731679281888926;;; -0.006919175959342079 -0.004731679321795338 … 0.00149569603067306 -0.004731679290200443; -0.004731679339949524 -0.0029911484206769125 … 0.0017333319121935493 -0.0029911483921702536; … ; 0.0014956959713097705 0.0017333319102936064 … 0.002146714327371255 0.0017333318348272866; -0.0047316792950494665 -0.0029911483738199185 … 0.0017333318838007573 -0.002991148366987339;;; 0.0006640002274163597 0.001495696005907353 … 0.003561045093577197 0.0014956960130805995; 0.001495695901027539 0.0017333318368994155 … 0.0021467144797033258 0.0017333318438774968; … ; 0.0035610449278620027 0.002146714483926789 … 1.0228461138264006e-5 0.002146714347249506; 0.0014956960147151296 0.0017333319498449538 … 0.002146714468680674 0.0017333319111101082;;; … ;;; -0.0009170692262965938 0.0007544794620939098 … 0.004144156313797027 0.0007544793254885751; 0.0007544793876760776 0.0029528261019312685 … 0.0015077442847549697 0.002952825910753144; … ; 0.004144156486237049 0.0015077442845262412 … -0.0015401031531881342 0.0015077443318001604; 0.0007544793212091207 0.0029528260210000034 … 0.0015077442287153497 0.0029528258786088463;;; 0.000664000094992903 0.0014956957089066018 … 0.0035610447619957497 0.0014956959102998624; 0.0014956958258524672 0.0017333316261309642 … 0.002146714292344778 0.0017333317950333754; … ; 0.003561044840229077 0.002146714249564711 … 1.0228394609299377e-5 0.00214671436128898; 0.0014956958872668355 0.0017333316577687026 … 0.002146714278434062 0.001733331840048172;;; -0.006919175992755591 -0.0047316794140698434 … 0.001495695949712923 -0.004731679304361476; -0.004731679349228135 -0.00299114847867479 … 0.001733331849275719 -0.002991148388771802; … ; 0.0014956959654535756 0.0017333317883888814 … 0.0021467143595846833 0.001733331899912178; -0.004731679324401793 -0.00299114847582844 … 0.001733331849283484 -0.002991148362232021]), 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.1682595838539101e-6 - 2.526311752018857e-6im 0.0013685742976595435 + 0.0024716992448842526im … 0.000969976122351745 + 0.0016075132207258034im -0.0014049067066797642 - 0.0023918014613803973im; -0.01682001483456825 - 0.008219670164626224im 0.01447618495679229 + 0.006000506967890838im … -0.009980078833516636 - 0.005629997697788154im 0.016101412033439667 + 0.008850448961197875im; … ; -0.020198327585306312 - 0.00988066806938396im 0.017955009057322086 + 0.009142007540216511im … -0.011087111079059193 - 0.0050479259311638595im 0.01739707400200877 + 0.008163637673684913im; 0.016845205389874954 + 0.008156160387631124im -0.016090833119359965 - 0.008868296900444383im … 0.008643481061109319 + 0.0034291251872022852im -0.014465610611274424 - 0.00601836253308639im;;; -0.005619839057354987 - 0.0059204623592462105im 0.0024931830280830376 + 0.0020178034382191293im … -0.0017995288978833097 - 0.002289784227150528im 0.0041184105227400035 + 0.004867759569200585im; 0.018571326328897825 + 0.010918013314773253im -0.014695835223713557 - 0.007672788858186952im … 0.009247258144056032 + 0.005573868445686823im -0.0157422551385822 - 0.009507750851135605im; … ; 0.016549146867744183 + 0.007453967663072929im -0.014560104749887599 - 0.0069021579153933504im … 0.009169332911163793 + 0.003883957749195996im -0.014122132108922627 - 0.00613415396004521im; -0.011995694612045575 - 0.0039507962453641954im 0.011488322583576401 + 0.005064055917715759im … -0.006469277287556107 - 0.0020711767705039407im 0.010441907678821756 + 0.0032290999691609048im;;; 0.0038444586942728213 + 0.004014007423534098im -0.002367555068231028 - 0.0022541005648875276im … 0.001593263040295256 + 0.001846278576399653im -0.0029254927704432518 - 0.0032324812189299893im; -0.011953969855779532 - 0.007282364033768733im 0.009617046602761754 + 0.005482087372526805im … -0.006365384104454255 - 0.003989093379102089im 0.01005502045537415 + 0.006250098351773951im; … ; -0.010138641434942416 - 0.004254062244393795im 0.009008748987255207 + 0.003980263058962931im … -0.006039581550663291 - 0.002368222229098267im 0.008756992434214496 + 0.003538816045229185im; 0.006669589440872849 + 0.0017767559186941046im -0.006099485345854379 - 0.0021629531813896637im … 0.0036904751092789964 + 0.0009024077117587544im -0.005661510998659719 - 0.0013949461098989755im;;; … ;;; 0.0016673407574613315 + 0.001834918490558246im -0.0015148293042943845 - 0.001746458155983902im … 0.0007760627949297107 + 0.0007749016408400773im -0.0012955483492569592 - 0.001361962430672197im; 0.0029947798965970617 + 0.0005869030110194473im -0.0025929362296558643 - 0.00038636035567591155im … 0.002120906920201341 + 0.0006739371406977149im -0.0027832279195090666 - 0.0007200411901118837im; … ; 0.0071318448389279055 + 0.004058584598231411im -0.006493593153035302 - 0.0037730331478724993im … 0.004756896550384611 + 0.0025858503775493106im -0.006359684936000393 - 0.003538210973295182im; -0.005805705155237238 - 0.0036939484483329847im 0.005198456535908195 + 0.003403636242663277im … -0.003501413883507094 - 0.002060947668415538im 0.005008159925599595 + 0.00306994935758257im;;; -0.003796447913199273 - 0.004114201156631688im 0.00287932973509144 + 0.0033273760270400897im … -0.0009278574916724907 - 0.0008356101784212906im 0.0023213933865592748 + 0.002349017845235296im; -0.006707081519815004 - 0.0017008788402737849im 0.0057006544476843665 + 0.001315956542133939im … -0.004252400335672395 - 0.0017413754867269776im 0.006138623294469654 + 0.002083953203032069im; … ; -0.012659757562742222 - 0.006935940842141191im 0.011121721096954262 + 0.006227976699927961im … -0.007418446932633473 - 0.00383500746111805im 0.010869977931139117 + 0.005786527585034793im; 0.011916475818883478 + 0.0073582434019136365im -0.010015877481488317 - 0.006329087241269389im … 0.005803454108118635 + 0.0031501260260698863im -0.009577905451229747 - 0.005561089545710342im;;; 0.005584155635232919 + 0.005998105302228787im -0.00407215685480088 - 0.004966210651628954im … 0.0004543043846945166 + 0.0001077742221873362im -0.002446936144860872 - 0.002116292722787363im; 0.012005981337392089 + 0.003933982186665625im -0.01046872206569298 - 0.00317612483453258im … 0.007478043430932765 + 0.0036918313304693936im -0.011515135634600406 - 0.005011063712584279im; … ; 0.018802978628871446 + 0.009851511603593256im -0.016329318068682312 - 0.008784180873912327im … 0.01018627584223375 + 0.004965743420765661im -0.015891347570585673 - 0.00801618135510817im; -0.018561038609690553 - 0.01093482664576867im 0.015715437792488338 + 0.009560719285683552im … -0.008238490480163492 - 0.003953208233831466im 0.014669024284395125 + 0.00772578318753737im],)])]), [8.519874603160076e-7 4.627265729833096e-7 … 2.4177792876559738e-8 4.6272656543345317e-7; 4.627265696940199e-7 2.501448085312968e-7 … 1.3583004900899008e-8 2.5014480457138857e-7; … ; 2.4177793044990804e-8 1.3583005246042033e-8 … 1.565925013996935e-9 1.358300491088764e-8; 4.6272656684604134e-7 2.5014480757449195e-7 … 1.3583005053315324e-8 2.5014480308270816e-7;;; 4.6272656784846184e-7 2.501448060648798e-7 … 1.3583004488219125e-8 2.5014480374128734e-7; 2.501448075287815e-7 1.389561145675701e-7 … 1.1279996530134304e-8 1.389561131819204e-7; … ; 1.3583005084438855e-8 1.1279996525969689e-8 … 7.937761301647545e-9 1.1279997225478508e-8; 2.5014480405567013e-7 1.3895611215813036e-7 … 1.127999678792096e-8 1.3895611188700988e-7;;; 2.4177793069989366e-8 1.3583004721338982e-8 … 1.5659245096713234e-9 1.3583004663555648e-8; 1.3583005811423015e-8 1.1279997194048417e-8 … 7.93776021644537e-9 1.1279997144992855e-8; … ; 1.565924895072189e-9 7.937760169268732e-9 … 3.571857425060436e-8 7.937761155536913e-9; 1.3583004638494059e-8 1.127999616847335e-8 … 7.937760295055956e-9 1.12799965345154e-8;;; … ;;; 5.791840683876558e-8 2.281601411423034e-8 … 5.620392063460848e-10 2.281601615131283e-8; 2.281601522286529e-8 3.4892734731251593e-9 … 1.345852011685971e-8 3.489274255810386e-9; … ; 5.620390030965982e-10 1.345852009481452e-8 … 7.728317711938626e-8 1.3458519626465016e-8; 2.2816016202731614e-8 3.489273801741216e-9 … 1.3458520691326746e-8 3.489274386320073e-9;;; 2.4177795089250363e-8 1.35830077769795e-8 … 1.5659252862898008e-9 1.3583005715269941e-8; 1.3583006578445622e-8 1.1279999100308034e-8 … 7.937761546681762e-9 1.1279997580845217e-8; … ; 1.5659250989004037e-9 7.937761835042522e-9 … 3.571857556919823e-8 7.937761049348101e-9; 1.3583005944820254e-8 1.127999881301321e-8 … 7.937761645889338e-9 1.1279997172107043e-8;;; 4.627265717651046e-7 2.5014481317006625e-7 … 1.3583005318996626e-8 2.501448048039625e-7; 2.5014480821424657e-7 1.3895611752830466e-7 … 1.12799970974728e-8 1.3895611298544155e-7; … ; 1.3583005140535035e-8 1.1279997629058703e-8 … 7.937761068757304e-9 1.1279996630658946e-8; 2.501448062938054e-7 1.38956117381745e-7 … 1.127999709739004e-8 1.389561116207883e-7;;;;], Matrix{ComplexF64}[[0.11395450117843973 - 0.23426767190632833im -0.15474559855874706 + 0.9525929930527698im -9.120399871466109e-8 - 6.785033880498543e-7im -2.6470159199581057e-6 + 2.319092074742018e-6im; -0.07731572508229423 + 0.15894567337931048im -0.006818701436552305 + 0.04197498524020271im -0.46729354848912485 - 0.36104119093894776im -0.0564223783647533 - 0.03678496360920907im; … ; 0.02456378356813884 - 0.050498228047451216im 0.0024662352374016736 - 0.015181811507408993im -0.00796675250122037 - 0.004867150522644461im 0.005593537055783114 - 0.002652325049068665im; -0.044043032186177306 + 0.09054366871580155im -0.004885800371948571 + 0.030076331373802608im 0.019145691759324908 + 0.007629715555357706im -0.03412065995441429 + 0.012682551842466104im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944063911946, [[-0.5541104223294051, -0.01067839045298409, 0.14495656498248474, 0.14495678842904713]], 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.6206946466389894e-7, -6.293617201899408e-7, -2.2951215151941526e-7, 4.99402589199648e-7, 7.923312430822712e-7, 5.851072492716465e-7, 6.516971206189764e-7, 1.0189308997406666e-6, 3.638591716132273e-7, -2.201093488836659e-6  …  2.7810182699022314e-6, 1.201695554333045e-6, 2.4456416146130962e-8, 1.0078551446353204e-7, 5.535308754933765e-7, 4.5341735129397143e-7, 4.6319719630417455e-9, -1.0437048274483227e-7, 1.1211100508004006e-7, 1.255086491960004e-7], resid_history = [0.24939209480867475, 0.0037665477416317958, 0.0002852302424692876, 4.681812439277999e-6, 5.786539001864306e-8, 4.931908817973217e-11, 4.4564725665270646e-8, 1.329739787947426e-9], converged = true, n_iter = 8, residual_norm = 1.329739787947426e-9, maxiter = 100, tol = 1.0e-8, s = 0.9355286696768914, residual = [8.174742419820129e-7, 5.35922843196094e-7, 1.105495944322714e-7, -1.591541912711447e-7, -1.841410868004234e-7, -1.0463909284574704e-7, -9.27392352940967e-8, -1.215373607254159e-7, -3.8929431031636005e-8, 2.170693161681071e-7  …  8.133625871779562e-7, 3.5988429960992015e-7, 7.81859991615216e-9, 3.485026087465538e-8, 2.048759022869878e-7, 1.8160638507222826e-7, 2.142593713469613e-9, -6.293328922990845e-8, 1.0790671238153603e-7, 3.548710509427922e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.741051473432191e-7, -4.3660847431272255e-8, 1.329739787947426e-9, -0.00020442898092665857, 4.686830405871563e-6, -5.848072696499532e-8, 4.931908817973217e-11, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], V = KrylovKit.OrthonormalBasis{Vector{Float64}}([[-9.558782541594615e-7, -2.2187822835549015e-6, -7.256028942733618e-7, 1.2117169920247319e-6, 1.5296658828122104e-6, 1.0812684734088287e-6, 1.2691372180786967e-6, 1.9030422097708721e-6, 5.577734614053608e-7, -2.7519512148114448e-6  …  2.982426136586195e-6, 1.1773738688517835e-6, 2.713883030239599e-8, 1.2930245508531666e-7, 6.913875996072704e-7, 4.6477277669116444e-7, 3.9901107532589e-9, -9.896937845214147e-8, 1.3930683401159797e-7, 1.4654888296827827e-7], [1.8065475657176837e-6, 3.3801831355147654e-6, 1.0825958336998009e-6, -1.8226445205419703e-6, -2.312369911368237e-6, -1.6132288097621165e-6, -1.8661262453960576e-6, -2.7979209831691236e-6, -8.307945046634683e-7, 4.133912091192013e-6  …  -3.6375034869651187e-6, -1.4010334879502613e-6, -3.203120170633038e-8, -1.531178580482709e-7, -8.087804296484737e-7, -5.152752808816524e-7, -3.961482798201076e-9, 8.86627189886437e-8, -1.1010744219147133e-7, 5.079530013836444e-8]]), H = [1.1585084585459395 0.009266418634471479 … 0.0 0.0; 0.11987026507445601 1.0261966777245952 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1646934055671085 0.11483339006346763 … 0.0 0.0; 0.0 1.0202477784851378 … 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.9257125299024138
Interacting polarizability:     1.7736548603371818

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

Manual solution of the Dyson equations

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

using KrylovKit

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

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

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

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

[ Info: GMRES linsolve starts with norm of residual = 4.19e+00
[ Info: GMRES linsolve in iteration 1; step 1: normres = 2.49e-01
[ Info: GMRES linsolve in iteration 1; step 2: normres = 3.77e-03
[ Info: GMRES linsolve in iteration 1; step 3: normres = 2.85e-04
[ Info: GMRES linsolve in iteration 1; step 4: normres = 4.69e-06
[ Info: GMRES linsolve in iteration 1; step 5: normres = 1.09e-08
[ Info: GMRES linsolve in iteration 1; step 6: normres = 6.27e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 9.11e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.12e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 4.18e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125299023214
Interacting polarizability:     1.7736548566282875

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