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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457355029187126
Displaced dipole:  0.0176122211136519
Polarizability :   1.774679466394377

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.0007087076610060837 + 0.0032187968589272318im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.30270707982020495 - 0.06398920385512469im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.04013982828848882 - 0.008668392453412829im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.07897764692720566 + 0.017706583962228073im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206934551233504e-7 -2.5020144143339927e-7 … -4.88064076806022e-8 -2.502014374884117e-7; -6.293614574675233e-7 -4.793034678635184e-7 … -1.207783657328929e-7 -4.79303470106466e-7; … ; 1.362625352937342e-7 1.135622227529126e-7 … 4.711984326058594e-8 1.1356223019418198e-7; 1.3305467766748252e-7 1.3632633046753252e-7 … 5.200682184877714e-8 1.363263362648332e-7;;; -2.502014775556741e-7 -1.7404905179836998e-7 … -3.556492830381981e-8 -1.7404904950420742e-7; -4.793034978244236e-7 -3.679730745071002e-7 … -1.1202887877253437e-7 -3.679730792649148e-7; … ; 1.1356222616793677e-7 1.1211090201400023e-7 … 1.0980675628820132e-7 1.1211090995870957e-7; 1.36326309071432e-7 1.255085797406517e-7 … 5.094904819726284e-8 1.2550858443300744e-7;;; -4.880642840473476e-8 -3.556493683559404e-8 … -1.1412139945253544e-8 -3.556493773369821e-8; -1.2077838079088618e-7 -1.1202888120807433e-7 … -9.786219172703815e-8 -1.1202888925729672e-7; … ; 4.71198486549902e-8 1.0980675036169091e-7 … 2.493821851254604e-7 1.0980675895609356e-7; 5.200681469107704e-8 5.0949044003030256e-8 … 4.9106236839863044e-8 5.094904717409302e-8;;; … ;;; 6.475755994244902e-8 4.1103307833300786e-8 … -6.681689173248004e-9 4.110330463232156e-8; 1.6897552351614207e-7 6.653284482474921e-8 … -1.329444466446963e-7 6.653283919278941e-8; … ; 3.2677620365921245e-8 1.6017364393079015e-7 … 3.84703597642227e-7 1.601736520597896e-7; -9.141992904153999e-8 -3.583035497471537e-8 … 7.037284645278755e-8 -3.583035077609235e-8;;; -4.880639407393567e-8 -3.556491092538843e-8 … -1.1412132764596148e-8 -3.5564912524138574e-8; -1.2077835800314826e-7 -1.1202886291842704e-7 … -9.786217849728836e-8 -1.120288681826194e-7; … ; 4.711984642492442e-8 1.0980675880263174e-7 … 2.4938221830212413e-7 1.0980676835746113e-7; 5.20068344795596e-8 5.0949061052946353e-8 … 4.910625483989305e-8 5.094906607167655e-8;;; -2.502014027620189e-7 -1.7404899694661928e-7 … -3.556491641245054e-8 -1.740489950539808e-7; -4.793034426143412e-7 -3.679730349448537e-7 … -1.1202886941680992e-7 -3.679730378620693e-7; … ; 1.1356222986997796e-7 1.1211090485488117e-7 … 1.0980675977213485e-7 1.1211091320698905e-7; 1.3632635898630356e-7 1.2550861550387462e-7 … 5.094905812063943e-8 1.2550862103452879e-7;;;;], δHtotψ = Matrix{ComplexF64}[[-0.0004122462051604816 - 0.0018723334299527736im -0.07682274423352094 + 0.056389635916744736im -0.8555089678829582 - 0.48733085854526303im 0.29314456974080705 - 0.811565361402355im; -0.19404264060533033 + 0.04087838539543679im 0.8843939622961978 + 1.3654751567500856im -0.1890845826463585 - 0.10693413005317165im 0.06481062427011285 - 0.17980179191798842im; … ; -0.047943244330159486 + 0.010217557757919826im -0.012282473706433724 - 0.016645278033678882im 0.001388485000650599 + 0.0109917334441005im -0.002396068521840158 + 0.015791120394479635im; 0.0706224307938937 - 0.015828422596831707im 0.015096617755695029 + 0.02066385274608719im -0.13878105806113414 + 0.0012706573115913505im 0.009394154811344852 + 0.044467672348596im]], δVind = [0.0019654318168764373 0.0020810005657302746 … 0.003153664827314846 0.0020810005180819614; 0.0038891304853428147 0.0047792691324768845 … 0.010276671420719969 0.00477926907615338; … ; -0.006150846823126199 -0.008269602782240041 … -0.017302973073770572 -0.008269602693618882; 0.00023868260649067934 -0.0003986969998518964 … -0.003664694914613332 -0.0003986970571512138;;; 0.002081000842202971 0.002223228615474248 … 0.00342012021945501 0.002223228563809603; 0.0047792694403354635 0.005734953791651438 … 0.010903598890182626 0.005734953713583632; … ; -0.008269601956339289 -0.009601007194868758 … -0.012623111765864533 -0.009601007061258304; -0.00039869667960246245 -0.0010730122713709222 … -0.004230680332916991 -0.001073012332111579;;; 0.003153665891866854 0.0034201208870755936 … 0.004832507513956689 0.0034201207582355143; 0.01027667259966492 0.010903599744478304 … 0.012272143707787647 0.010903599358981325; … ; -0.01730296960041118 -0.012623110917730202 … -0.009620264319933796 -0.012623110677596014; -0.0036646934280900726 -0.004230679456196238 … -0.0054911076156857 -0.004230679472207044;;; … ;;; -0.0023365590499518024 -0.002815461542631361 … 0.005693314433895699 -0.002815461524218554; -0.012913091199092447 -0.017841178204274735 … 0.011498337828917135 -0.017841178617783976; … ; -0.025174467376002777 -0.012810887706630341 … -0.00839587598668004 -0.012810887567867392; 0.007597980046333556 0.010225414753698949 … -0.005449734997274402 0.010225415315543534;;; 0.0031536637984946863 0.0034201185860882046 … 0.004832503867677121 0.0034201184130374074; 0.01027667015985917 0.01090359714514668 … 0.012272140824775816 0.010903596846016776; … ; -0.017302974727603716 -0.012623114003456677 … -0.009620266314218816 -0.01262311367701288; -0.0036646959469882794 -0.004230682093600774 … -0.005491110199366373 -0.004230682031015422;;; 0.002081000237563378 0.0022232279424849416 … 0.0034201190323724408 0.002223227883263387; 0.004779268758885506 0.005734953028165974 … 0.010903597589443755 0.005734952953159771; … ; -0.008269603420583664 -0.009601008737849987 … -0.01262311331608912 -0.00960100856953694; -0.0003986973747450589 -0.0010730130561303886 … -0.004230681578535173 -0.0010730131087382882;;;;], δρ0 = [-3.602139624223367e-7 -2.4888422610162603e-7 … -4.852835354293527e-8 -2.4888422824827724e-7; -5.972249057826274e-7 -4.606977163231385e-7 … -1.1905762969582989e-7 -4.606977227451427e-7; … ; 1.2735785332052234e-7 1.0836334212814984e-7 … 4.6407133427548495e-8 1.0836334828742695e-7; 1.0352004309181913e-7 1.1957958406852548e-7 … 5.069503840196445e-8 1.1957958531164758e-7;;; -2.488842280718266e-7 -1.7310668932270873e-7 … -3.5356963555290886e-8 -1.731066914125053e-7; -4.6069772116890476e-7 -3.5721697770454624e-7 … -1.1091433988480975e-7 -3.5721698548830077e-7; … ; 1.0836335083119315e-7 1.0846491742282528e-7 … 1.0862231381231172e-7 1.0846492423786068e-7; 1.1957958725592804e-7 1.1611144180562549e-7 … 5.020828168523834e-8 1.1611144315777601e-7;;; -4.8528356145722765e-8 -3.5356964317984154e-8 … -1.1341305306496435e-8 -3.53569662045512e-8; -1.1905763124662742e-7 -1.109143353272573e-7 … -9.770663677221105e-8 -1.10914344201327e-7; … ; 4.640714317603692e-8 1.086223136955629e-7 … 2.4845238677346057e-7 1.0862232145741936e-7; 5.069504456131437e-8 5.020828447375945e-8 … 4.92685619113413e-8 5.0208286710852615e-8;;; … ;;; 6.434570113896078e-8 4.083862270987045e-8 … -6.63792956535274e-9 4.0838620897570166e-8; 1.6975033457794372e-7 6.683084652604283e-8 … -1.3340479736239736e-7 6.683084139904928e-8; … ; 3.270997163207316e-8 1.6023520279367487e-7 … 3.8390139273909634e-7 1.6023520996771044e-7; -9.271365781604922e-8 -3.633620911987599e-8 … 7.126308135563971e-8 -3.633620433781663e-8;;; -4.8528355262500884e-8 -3.535696330931561e-8 … -1.1341306487200992e-8 -3.535696594386328e-8; -1.1905763334708927e-7 -1.1091433957101537e-7 … -9.770664229660095e-8 -1.1091434571804413e-7; … ; 4.6407132760187733e-8 1.0862230349900088e-7 … 2.484523804789051e-7 1.0862231222972607e-7; 5.069503955826893e-8 5.020827901856562e-8 … 4.926856114515253e-8 5.020828305664257e-8;;; -2.488842275991188e-7 -1.7310668864589467e-7 … -3.535696410241429e-8 -1.731066912486531e-7; -4.606977202257194e-7 -3.572169782365773e-7 … -1.1091434178757022e-7 -3.5721698429045907e-7; … ; 1.0836334252260272e-7 1.0846490926659021e-7 … 1.0862230800471371e-7 1.0846491649230817e-7; 1.1957958331646624e-7 1.1611143769489873e-7 … 5.020828034935399e-8 1.161114397982376e-7;;;;], δeigenvalues = [[0.0004950052655174405, 0.09679860002714546, 0.042775343889577884, 0.04455212620892875]], δ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.1547596528350744 -0.15461620271043983 … -0.15418753101004426 -0.15461620271038742; -0.15461620271048698 -0.1544722989097573 … -0.15404229468715333 -0.1544722989097032; … ; -0.15418753100996052 -0.1540422946870212 … -0.15360833735789092 -0.1540422946869722; -0.15461620271034354 -0.1544722989096119 … -0.1540422946870118 -0.1544722989095612;;; -0.1546162027104848 -0.15447229890975273 … -0.1540422946871576 -0.1544722989097042; -0.15447229890980366 -0.15432794310108958 … -0.1538966028399049 -0.15432794310103926; … ; -0.1540422946870661 -0.1538966028397613 … -0.1534612870216613 -0.15389660283971593; -0.15447229890965639 -0.15432794310094022 … -0.1538966028397595 -0.1543279431008934;;; -0.1541875310102404 -0.15404229468730019 … -0.15360833735817903 -0.15404229468725517; -0.15404229468735492 -0.15389660284005177 … -0.15346128702195258 -0.15389660284000506; … ; -0.15360833735807958 -0.15346128702179723 … -0.1530219131463961 -0.15346128702175513; -0.15404229468720357 -0.15389660283989837 … -0.1534612870218027 -0.1538966028398549;;; … ;;; -0.1534788058656963 -0.15333139197558132 … -0.15289092493848314 -0.15333139197551646; -0.15333139197561965 -0.15318350575354206 … -0.15274162879027814 -0.15318350575347514; … ; -0.15289092493841935 -0.15274162879017547 … -0.15229553960147713 -0.15274162879011324; -0.15333139197548223 -0.15318350575340245 … -0.15274162879014266 -0.1531835057533391;;; -0.15418753100995874 -0.15404229468702502 … -0.1536083373578767 -0.1540422946869644; -0.15404229468706584 -0.15389660283976916 … -0.15346128702164274 -0.1538966028397069; … ; -0.1536083373578068 -0.15346128702153106 … -0.1530219131461034 -0.15346128702147357; -0.15404229468692746 -0.1538966028396287 … -0.1534612870215062 -0.15389660283956974;;; -0.1546162027103439 -0.15447229890961486 … -0.15404229468700562 -0.15447229890955835; -0.15447229890965858 -0.15432794310094755 … -0.153896602839749 -0.15432794310088946; … ; -0.15404229468692915 -0.15389660283962756 … -0.1534612870215142 -0.1538966028395745; -0.15447229890951825 -0.15432794310080516 … -0.15389660283961049 -0.1543279431007504]), 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.01516389687206731 -0.012550128333847894 … -0.004966867487447216 -0.012550128371274904; -0.01255012835600298 -0.010362594268545371 … -0.004135158231384717 -0.010362594332345793; … ; -0.004966867630031686 -0.004135158231332596 … -0.002069779350447105 -0.0041351583722479615; -0.012550128370434907 -0.010362594288203162 … -0.004135158170340533 -0.010362594335577822;;; -0.012550128411423313 -0.010362594328328301 … -0.004135158210420445 -0.010362594373171451; -0.010362594329243227 -0.00862203818240214 … -0.0038975105664565083 -0.008622038275088457; … ; -0.0041351584323064405 -0.0038975106517830944 … -0.0034841062678511526 -0.00389751079949212; -0.010362594399650104 -0.008622038275195563 … -0.0038975105709447054 -0.008622038324384247;;; -0.004966867695480716 -0.004135158300010092 … -0.002069779375590659 -0.004135158419737325; -0.0041351582778241785 -0.0038975104808799395 … -0.0034841061743194094 -0.003897510667198231; … ; -0.0020697796311959886 -0.0034841062694600682 … -0.005620579388640203 -0.003484106421684703; -0.0041351584701358555 -0.003897510703941701 … -0.003484106254764219 -0.003897510811394971;;; … ;;; -0.006547897393911114 -0.004876339745291039 … -0.0014866442813321569 -0.00487633958546182; -0.004876339627816806 -0.0026779885054932866 … -0.0041230541324993 -0.002677988367558333; … ; -0.0014866440874947862 -0.004123053974335116 … -0.007170894864565618 -0.004123054088215409; -0.004876339676713747 -0.0026779886100602546 … -0.004123054125142302 -0.0026779883829591965;;; -0.004966867587601996 -0.004135158179739608 … -0.002069779492623859 -0.004135158361736859; -0.0041351583146084885 -0.003897510565697415 … -0.003484106295102313 -0.003897510694421596; … ; -0.0020697793356719042 -0.0034841060697611117 … -0.005620579318907369 -0.003484106245257013; -0.004135158244272819 -0.003897510461199908 … -0.003484106242064248 -0.003897510659530085;;; -0.012550128335647576 -0.010362594247217185 … -0.00413515822783535 -0.01036259431840832; -0.010362594302667987 -0.008622038177481374 … -0.0038975106056106297 -0.00862203824901626; … ; -0.004135158240588582 -0.0038975104721473723 … -0.0034841061563571805 -0.003897510631506059; -0.010362594276347649 -0.008622038137792203 … -0.0038975105227543544 -0.008622038225975437])], 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.009532891638287098 -0.006919174470479829 … 0.0006640018651868553 -0.006919174507854408; -0.006919174492681932 -0.004731677770614855 … 0.0014956977474417377 -0.004731677834361181; … ; 0.0006640017226862069 0.0014956977476259193 … 0.003561047244225795 0.0014956976067595426; -0.006919174506970363 -0.004731677790127207 … 0.0014956978086274747 -0.004731677837451185;;; -0.006919174548100239 -0.004731677830393372 … 0.0014956977684015965 -0.004731677875187978; -0.004731677831359063 -0.0029911467808340295 … 0.0017333337371709442 -0.0029911468734700256; … ; 0.0014956975466071942 0.0017333336519878822 … 0.0021467161610139096 0.001733333504324237; -0.004731677901618641 -0.002991146873478071 … 0.001733333732828214 -0.002991146922619932;;; 0.0006640016569571788 0.0014956976786693688 … 0.0035610472187941377 0.0014956975589871559; 0.0014956977008006318 0.0017333338226006583 … 0.0021467162542544134 0.0017333336363290794; … ; 0.003561046963288284 0.002146716159269047 … 1.0229883919064299e-5 0.0021467160070865175; 0.0014956975086403892 0.001733333599692302 … 0.0021467161739595394 0.001733333492282497;;; … ;;; -0.0009170693340377258 0.0007544788589359998 … 0.0041441580728308934 0.0007544790188300839; 0.0007544789763719584 0.002952824944537877 … 0.0015077454663649835 0.0029528250825397494; … ; 0.004144158266732074 0.00150774562463178 … -0.0015401019435107507 0.0015077455108137147; 0.0007544789276124628 0.002952824840110575 … 0.0015077454738575396 0.0029528250672749713;;; 0.0006640017651175622 0.0014956977992149944 … 0.0035610471020632514 0.001495697617278361; 0.001495697664305396 0.00173333373806579 … 0.0021467161337813453 0.0017333336094038647; … ; 0.0035610472590851226 0.0021467163592341796 … 1.0229953944580934e-5 0.00214671618379576; 0.0014956977347795105 0.001733333842703768 … 0.0021467161869560235 0.0017333336444325439;;; -0.006919174472183615 -0.004731677749144394 … 0.0014956977511386806 -0.004731677820278991; -0.004731677804638745 -0.0029911467757712095 … 0.001733333698172726 -0.0029911468472480027; … ; 0.001495697738461971 0.0017333338317573584 … 0.0021467162726549584 0.0017333336724517402; -0.004731677778178046 -0.002991146735939653 … 0.0017333337811675847 -0.0029911468240681246]), 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 = [-5.71030574239101e-7 - 2.740572974743045e-6im -0.0009384827574427641 - 0.006714919758282853im … -0.000649778839034912 - 0.004651977644675404im 0.0009450531447669113 + 0.006783735731471921im; -0.02607971031646504 + 0.0057119299858305675im 0.02171649229341023 - 0.0007131916554475707im … -0.011233226466663701 + 0.005771400719744923im 0.02061377599624029 - 0.008623294242622773im; … ; -0.02729120032683095 + 0.005954586329459861im 0.02269462399069434 - 0.006345251283592694im … -0.014377296046161368 + 0.00155443087126803im 0.023072664495140532 - 0.003629270814855424im; 0.026071364757277265 - 0.00577061498139104im -0.0206157761099268 + 0.008607731979809085im … 0.012128968777854077 + 0.0006277756364902687im -0.02171863108596045 + 0.0006971221263789417im;;; 0.006929215515984259 + 0.001416173976296772im -0.0035122295977856217 + 0.00310210293431359im … 0.00184142065802334 + 0.003485365980861806im -0.004614866367728607 - 0.004809324230992834im; 0.017108308700604483 - 0.005476389610858608im -0.014898158764267373 + 0.0017767294594479217im … 0.008335711263159802 - 0.004767347117080024im -0.014188564419858 + 0.006870482213954595im; … ; 0.024275274520553357 - 0.00471121361275806im -0.020238859538501405 + 0.005048925206322313im … 0.013100526231531476 - 0.0012655190642219144im -0.02053563537890791 + 0.002916878310541546im; -0.02522396125707703 + 0.003844902531575576im 0.019420664480752955 - 0.005779940982169178im … -0.011194934800747867 - 0.0004671853658920497im 0.020130331457661934 - 0.0006859963341112204im;;; -0.004731523573839024 - 0.0009912999108730447im 0.003062614215293507 - 0.0006646803324760022im … -0.0017778738729372374 - 0.0018721072530125749im 0.0034405051219784533 + 0.0020518843782130316im; -0.008431233161774696 + 0.003256427416463604im 0.007485537488148557 - 0.0016062604717889074im … -0.004676897685315281 + 0.0029315744610048998im 0.0071888102475011935 - 0.003738565679455907im; … ; -0.01572297022162788 + 0.002764275942529412im 0.013515107962458788 - 0.003016784092743633im … -0.009447961715121545 + 0.0008271755381411439im 0.013685636431325155 - 0.0017913016054539874im; 0.014930836211958001 - 0.0018878023650064858im -0.012045001828881954 + 0.0026941485858897525im … 0.007659228256223877 + 0.00021543610234834706im -0.012341722564576855 + 0.0005618852720309061im;;; … ;;; -0.0020725146400467557 - 0.00036838591530550754im 0.0018228167579880843 + 0.000839415781565319im … -0.0010559065830523483 + 0.0005474707512303793im 0.0016744998801048769 - 0.0002279161231955741im; 0.007143852259196896 - 0.0008815222121434456im -0.006334390894906975 + 0.00032718653178611214im … 0.004370317309329363 - 0.001249706614284059im -0.0062056730911215045 + 0.0012534363457417843im; … ; 0.006433621233787961 - 0.001864619594568749im -0.0057692562110205365 + 0.0020027227443255467im … 0.004489552195581828 - 0.0007644995239425411im -0.005859862117660163 + 0.0013508612881469781im; -0.003646545064390129 + 0.0014941655279621549im 0.0031971288561893525 - 0.0017713832022352802im … -0.0024885258195134433 + 0.0002603635741708617im 0.003325849556656327 - 0.0008450745168374916im;;; 0.004702360698851745 + 0.0008775930796825339im -0.0034150789194046646 - 0.001947136119908902im … 0.0013349655826653997 - 0.0012555789342110692im -0.0030377142729240644 + 0.0007685240476371887im; -0.01490918942004002 + 0.001973027000499261im 0.012320603311490572 - 0.0006487903262460989im … -0.007285816072668955 + 0.0024148845734465975im 0.012024341221112303 - 0.002780591371900519im; … ; -0.012610754495088744 + 0.0033805759780036596im 0.01090670035218786 - 0.0035332753781839015im … -0.007745597254302361 + 0.0011641605717747257im 0.01107711104829909 - 0.002307812067356914im; 0.008453331963956347 - 0.003170622538941005im -0.0072102080140879245 + 0.0036513615893555883im … 0.005050392569701081 - 0.00030115946568520216im -0.0075065553485516474 + 0.0015194791086530122im;;; -0.006904208309723039 - 0.0013248302556654523im 0.004593271058883691 + 0.004706725229022159im … -0.0009414705916592984 + 0.0029341902589413733im 0.003491361636894753 - 0.0032031925135194506im; 0.02521390265719036 - 0.0038702249443406302im -0.020116662683596048 + 0.0007433973810121648im … 0.0105202303689935 - 0.004334665716501833im -0.01940765784140929 + 0.005836414814601102im; … ; 0.02149332282077429 - 0.005261938869522437im -0.01805504156362652 + 0.005481257475602252im … 0.011844989044578126 - 0.0015140292760170725im -0.018351657550856486 + 0.003349249385382269im; -0.01711887876556742 + 0.005450087635514852im 0.014202515532981451 - 0.006812851473355221im … -0.0090104917139806 - 3.4644377933023475e-5im 0.014911641038568513 - 0.0017195610967864782im],)])]), [8.519871902242509e-7 4.627263893829586e-7 … 2.417776771220606e-8 4.627263938250776e-7; 4.627263920125027e-7 2.501446850012664e-7 … 1.358298660861245e-8 2.5014468993581656e-7; … ; 2.4177769899995938e-8 1.358298660807682e-8 … 1.5659194211598137e-9 1.358298806184157e-8; 4.627263937253858e-7 2.501446865216554e-7 … 1.3582985978841817e-8 2.5014469018578917e-7;;; 4.6272639859018315e-7 2.501446896250832e-7 … 1.3582986392334908e-8 2.5014469309340907e-7; 2.5014468969584387e-7 1.3895602952876197e-7 … 1.1279979811125999e-8 1.389560342941429e-7; … ; 1.3582988681442446e-8 1.1279980585853499e-8 … 7.937748107319565e-9 1.1279981926996356e-8; 2.501446951413599e-7 1.389560342996538e-7 … 1.1279979851876198e-8 1.389560368286463e-7;;; 2.4177770904233836e-8 1.3582987316593641e-8 … 1.5659194802029346e-9 1.3582988551770825e-8; 1.3582987087710185e-8 1.1279979034125747e-8 … 7.937747440190853e-9 1.127998072581927e-8; … ; 1.5659200804233515e-9 7.937748118795706e-9 … 3.5718545346729795e-8 7.937749204566223e-9; 1.3582989071712266e-8 1.1279981059435663e-8 … 7.937748013977083e-9 1.1279982035066315e-8;;; … ;;; 5.791840947804276e-8 2.2816022816714737e-8 … 5.620371342462896e-10 2.2816020460355532e-8; 2.2816021084790672e-8 3.489278120285528e-9 … 1.3458507847791068e-8 3.4892775601839878e-9; … ; 5.620369078562448e-10 1.3458506226374323e-8 … 7.728313518907605e-8 1.3458507393814362e-8; 2.281602180568141e-8 3.4892785448925933e-9 … 1.3458507772370532e-8 3.4892776227213385e-9;;; 2.417776924896124e-8 1.3582986075808405e-8 … 1.5659197550234743e-9 1.3582987953400718e-8; 1.358298746719829e-8 1.1279979804231366e-8 … 7.937748301694053e-9 1.1279980972995173e-8; … ; 1.5659193864646024e-9 7.937746694410347e-9 … 3.571854394631829e-8 7.937747946164575e-9; 1.3582986741574662e-8 1.1279978855438817e-8 … 7.937747923389801e-9 1.1279980656196834e-8;;; 4.6272638959656006e-7 2.501446833516668e-7 … 1.3582986571997336e-8 2.501446888578391e-7; 2.501446876404193e-7 1.3895602927575855e-7 … 1.1279980166626455e-8 1.3895603295367013e-7; … ; 1.3582986703565225e-8 1.1279978954834227e-8 … 7.93774731206926e-9 1.127998040174908e-8; 2.5014468560471073e-7 1.3895602723517958e-7 … 1.1279979414327768e-8 1.3895603176904712e-7;;;;], Matrix{ComplexF64}[[-0.05601742296753455 - 0.2544189821402997im -0.7779895163836476 + 0.5710621014196525im -6.688085598231486e-7 + 4.0806740568522814e-7im 1.9140288301763984e-6 - 1.2260760875744884e-6im; 0.03800663979126473 + 0.17261791232464604im -0.03428127768311829 + 0.0251632399581247im -0.15484314948270478 + 0.27083964568599234im -0.25628885067766194 - 0.09257063151962361im; … ; -0.012074993471931967 - 0.05484200069996459im 0.012399093932972055 - 0.00910121879533637im -0.0027645094049240622 + 0.011544750247276858im -0.002168836389936929 - 0.0019797809973419265im; 0.02165054595692233 + 0.09833208280696994im -0.024563554906577904 + 0.018030211419908078im 0.0071039546456298985 - 0.04957282947891544im -0.0016799480999072317 + 0.006053972017137212im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944072563396, [[-0.5541104233253091, -0.010678391187370071, 0.1449566437702413, 0.14495671494987442]], 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.6206934551233504e-7, -6.293614574675233e-7, -2.2951199367130265e-7, 4.994024966012261e-7, 7.923309330941329e-7, 5.851068867451924e-7, 6.516968763333041e-7, 1.0189308328453568e-6, 3.638593212946138e-7, -2.2010930597420304e-6  …  2.781018034816846e-6, 1.201695383210051e-6, 2.445626820066791e-8, 1.0078538820872157e-7, 5.53530830383239e-7, 4.534174030570173e-7, 4.6320248452955776e-9, -1.0437051582762379e-7, 1.1211091320698905e-7, 1.2550862103452879e-7], resid_history = [0.2493920954036487, 0.003766548907361272, 0.00028527869884337144, 4.680406918931959e-6, 4.669272969168979e-8, 5.559147333039899e-11, 4.351459883818028e-8, 1.0705926738908872e-9], converged = true, n_iter = 8, residual_norm = 1.0705926738908872e-9, maxiter = 100, tol = 1.0e-8, s = 0.9356513533233597, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.0730515288685187e-7, -4.268557240079202e-8, 1.0705926738908872e-9, -0.00020446320540815543, 4.682871590761997e-6, -4.7194231087725525e-8, 5.559147333039899e-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.8129463556290142e-7, -2.0014298711124666e-6, -7.091650777333924e-7, 1.144690159716854e-6, 1.3714333120170686e-6, 9.81950662439638e-7, 1.2016815497264806e-6, 1.7899248019693495e-6, 4.813482930051916e-7, -2.105211177465382e-6  …  3.004228545996525e-6, 1.165799041065157e-6, 2.853141193755598e-8, 1.4437565582817398e-7, 7.906431770403964e-7, 5.33690511937378e-7, 4.76573701560616e-9, -1.334446015625424e-7, 2.2350672437437035e-7, 4.4280887272999334e-7], [5.176512336993764e-7, 1.540636689221586e-6, 5.175203540018141e-7, -8.446528192831736e-7, -1.0211839575490376e-6, -7.114220613164028e-7, -8.452807161043015e-7, -1.2573018565928618e-6, -3.470797393811539e-7, 1.5564477564336568e-6  …  -1.4380701616857972e-6, -5.2889454112344e-7, -1.3368716942725271e-8, -7.065625967134764e-8, -3.8335480300138295e-7, -2.3587535258016708e-7, -1.7767631848270283e-9, 4.644122535298808e-8, -7.546014583120775e-8, -6.783267995349647e-8]]), H = [1.0849693767695294 0.015050680853800518 … 0.0 0.0; 0.14310121142662643 1.0296131333704979 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.094365800470496 0.14955549089858503 … 0.0 0.0; 0.0 1.0191131210208924 … 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.925712519602586
Interacting polarizability:     1.7736548505637184

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

Manual solution of the Dyson equations

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

using KrylovKit

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

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

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

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

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

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