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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457313408023864
Displaced dipole:  0.01761222205238181
Polarizability :   1.774679518646205

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.002467216834847583 + 0.0021853505471924367im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.20319412521628807 + 0.23332025557209896im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.027104210099627317 + 0.030849775332797973im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.05389796009883935 - 0.06038213087817989im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206926214952306e-7 -2.502013314926999e-7 … -4.880640798818143e-8 -2.5020141986923576e-7; -6.293614328137214e-7 -4.793034087659015e-7 … -1.207783494382052e-7 -4.793034750193999e-7; … ; 1.3626254591047698e-7 1.135622329314865e-7 … 4.711981476529985e-8 1.1356222057436178e-7; 1.3305477870114172e-7 1.3632643446350523e-7 … 5.200681326395224e-8 1.3632637216023979e-7;;; -2.5020134686608893e-7 -1.7404892013936826e-7 … -3.556491548368024e-8 -1.7404898544057982e-7; -4.79303419466157e-7 -3.67972985163544e-7 … -1.1202884979122167e-7 -3.679730361106563e-7; … ; 1.1356223339891925e-7 1.1211091564851661e-7 … 1.0980673604728312e-7 1.12110903097412e-7; 1.3632642515317448e-7 1.2550868771817885e-7 … 5.094904769302627e-8 1.2550864118944107e-7;;; -4.880637557017864e-8 -3.5564888172130375e-8 … -1.1412126618780038e-8 -3.5564904255881656e-8; -1.207783231403595e-7 -1.1202882197256031e-7 … -9.786215303000602e-8 -1.1202884076880992e-7; … ; 4.711982461717503e-8 1.0980675700125704e-7 … 2.493821919383357e-7 1.0980674585279018e-7; 5.200683904510058e-8 5.094907270508683e-8 … 4.9106238688371586e-8 5.0949060056318646e-8;;; … ;;; 6.475761760290571e-8 4.11033358918e-8 … -6.6816872371526686e-9 4.11033548269975e-8; 1.689756057035282e-7 6.653289213936434e-8 … -1.3294445342591565e-7 6.653290066050484e-8; … ; 3.267757715487977e-8 1.601736078220928e-7 … 3.847034980943069e-7 1.6017359260529122e-7; -9.141992861047504e-8 -3.5830367522586226e-8 … 7.037279405786434e-8 -3.5830360264824506e-8;;; -4.8806402560186255e-8 -3.5564908604458856e-8 … -1.1412134467189935e-8 -3.55649246589851e-8; -1.2077834909628334e-7 -1.1202884831100895e-7 … -9.786217548076247e-8 -1.1202886456400613e-7; … ; 4.7119816326769835e-8 1.098067386990375e-7 … 2.4938216020789553e-7 1.0980673083976088e-7; 5.2006818303025506e-8 5.094905300739029e-8 … 4.910622402757257e-8 5.094904174943804e-8;;; -2.502014051905758e-7 -1.7404896299720267e-7 … -3.556492612477921e-8 -1.7404902827157268e-7; -4.793034646301048e-7 -3.6797302026205556e-7 … -1.1202886381098793e-7 -3.6797307065710884e-7; … ; 1.1356221952953347e-7 1.1211090087947815e-7 … 1.0980672614997049e-7 1.12110891356485e-7; 1.3632638161636208e-7 1.255086540836764e-7 … 5.0949038374384136e-8 1.2550860890906588e-7;;;;], δHψ = Matrix{ComplexF64}[[0.0014351488916828667 - 0.0012711908643124744im 0.07065025959088048 - 0.06395361033933736im 0.7589630691256612 + 0.1255285833274998im 1.1212531208515237 - 1.469019398172775im; -0.13012986762614698 - 0.14963223660712174im -1.0194574174190727 - 1.2678256202214158im 0.16740264018964737 + 0.027627012191184273im 0.24699526869954502 - 0.3247827874249935im; … ; -0.03225459361784249 - 0.03691333574635146im 0.01392107657080608 + 0.015301338629778121im -0.013639811148891411 - 0.005846711588109218im -0.01875080104853562 + 0.024947160980834396im; 0.048191659507867156 + 0.05399658000470088im -0.01713152713532486 - 0.019010818599760303im -0.009892569877193893 - 0.0514087390840469im -0.030113907621550397 + 0.022618632018685934im]], δVind = [0.0019654314051930277 0.0020809997167082553 … 0.0031536654716153526 0.002081000439962274; 0.003889130472868457 0.004779268681742362 … 0.01027667357636769 0.004779269485883886; … ; -0.00615084945229806 -0.00826960660968187 … -0.017302978718445664 -0.00826960509648909; 0.0002386817419311462 -0.000398698402553905 … -0.0036646952562609165 -0.00039869757407480955;;; 0.002080999827142217 0.00222322704219365 … 0.003420119933681728 0.0022232278622604617; 0.004779268802789069 0.00573495260211164 … 0.010903600121748227 0.005734953503184803; … ; -0.008269606228871813 -0.009601012358057045 … -0.01262311521854714 -0.00960101084398861; -0.0003986982548833246 -0.0010730145234610633 … -0.004230681836531332 -0.001073013593303677;;; 0.003153663323460456 0.0034201173211777887 … 0.004832506311262136 0.0034201187209681153; 0.010276671248320018 0.010903597531294475 … 0.012272144069843534 0.010903598816813097; … ; -0.017302981876416247 -0.012623118022362058 … -0.00962026709006051 -0.012623116450187228; -0.00366469755772473 -0.004230684506822743 … -0.00549111057790132 -0.004230682989210164;;; … ;;; -0.002336560543983182 -0.002815462606801059 … 0.005693323971634141 -0.0028154638816460563; -0.012913092502824035 -0.017841177423168667 … 0.011498343362726622 -0.017841179127702438; … ; -0.025174478325854895 -0.01281088865058785 … -0.008395874532633924 -0.012810887250034285; 0.007597977979873798 0.010225411994376156 … -0.005449732471047464 0.010225410158956378;;; 0.003153665060202277 0.0034201192604327254 … 0.004832508978757043 0.00342012062402132; 0.010276673018055736 0.010903599274645724 … 0.012272145821050226 0.010903600606292307; … ; -0.017302979537459576 -0.012623116278553797 … -0.009620265628351502 -0.012623114524977306; -0.0036646957414570146 -0.004230682656113082 … -0.005491108527942951 -0.004230681094715478;;; 0.0020810003314506225 0.002223227616011382 … 0.003420120873392324 0.0022232284214970693; 0.004779269370512913 0.005734953237187474 … 0.010903600929177444 0.00573495411802131; … ; -0.008269605496865846 -0.009601011684693746 … -0.012623114364851801 -0.00960101004793753; -0.00039869771973443887 -0.0010730139299334092 … -0.0042306808855285686 -0.0010730129970424609;;;;], δρ0 = [-3.6021387770137076e-7 -2.4888416131748054e-7 … -4.852833246809323e-8 -2.4888416304600755e-7; -5.972248099667835e-7 -4.6069764152377787e-7 … -1.1905758764522787e-7 -4.6069764396903546e-7; … ; 1.2735783533948231e-7 1.0836332387292345e-7 … 4.640710401180662e-8 1.0836332603589599e-7; 1.035200747274222e-7 1.1957960601601333e-7 … 5.069503581948337e-8 1.1957960687016943e-7;;; -2.4888416161595754e-7 -1.7310663849076e-7 … -3.5356945740375805e-8 -1.7310664015332652e-7; -4.6069764091315346e-7 -3.572169137811541e-7 … -1.1091429768876904e-7 -3.572169171635765e-7; … ; 1.0836332673731784e-7 1.0846489580194975e-7 … 1.0862228290555585e-7 1.0846489655322729e-7; 1.1957960742248483e-7 1.1611145510657438e-7 … 5.0208276939542435e-8 1.1611145564839629e-7;;; -4.8528332487615185e-8 -3.535694528152402e-8 … -1.1341294034596587e-8 -3.535694665387956e-8; -1.1905758558712935e-7 -1.1091429433618031e-7 … -9.770659577108661e-8 -1.1091429965764504e-7; … ; 4.6407105657320275e-8 1.086222834980046e-7 … 2.484523545955027e-7 1.086222834984111e-7; 5.069503729561579e-8 5.020827754375422e-8 … 4.9268552328993366e-8 5.020827822727528e-8;;; … ;;; 6.434569075058248e-8 4.083861838687541e-8 … -6.63791793919787e-9 4.0838618454221104e-8; 1.6975036307880674e-7 6.683087706437452e-8 … -1.3340474828402763e-7 6.683087820902361e-8; … ; 3.270993127210181e-8 1.6023517329329225e-7 … 3.839013655999643e-7 1.6023517254046318e-7; -9.271368866430966e-8 -3.633623011185007e-8 … 7.126306817687732e-8 -3.633623022595984e-8;;; -4.852833396060331e-8 -3.535694664251219e-8 … -1.1341295486552297e-8 -3.535694805135577e-8; -1.190575924725591e-7 -1.109143033112387e-7 … -9.770660390152925e-8 -1.1091430620549871e-7; … ; 4.6407103810993015e-8 1.0862227965301045e-7 … 2.4845235325726883e-7 1.0862228289091563e-7; 5.069503568220954e-8 5.0208275217613855e-8 … 4.926855205953221e-8 5.0208277240062095e-8;;; -2.488841635976602e-7 -1.731066400039976e-7 … -3.535694694011968e-8 -1.7310664194960942e-7; -4.6069764493743865e-7 -3.572169181838355e-7 … -1.1091430311638075e-7 -3.572169213191934e-7; … ; 1.0836332256520462e-7 1.0846488987458913e-7 … 1.0862228022019496e-7 1.0846489352825352e-7; 1.1957960537593131e-7 1.1611145234303793e-7 … 5.0208276235616466e-8 1.1611145400788803e-7;;;;], δeigenvalues = [[0.0004950046351643645, 0.09679858894808525, 0.04609019477915319, 0.014835723926235429]], δ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.15475965283644935 -0.15461620271125306 … -0.15418753101253363 -0.15461620271232335; -0.1546162027121936 -0.154472298910901 … -0.1540422946899837 -0.1544722989119743; … ; -0.15418753101066424 -0.15404229468715627 … -0.15360833735971768 -0.15404229468824224; -0.1546162027113828 -0.15447229891008737 … -0.1540422946891654 -0.15447229891116201;;; -0.15461620271180357 -0.15447229891051006 … -0.15404229468958955 -0.1544722989115833; -0.1544722989114531 -0.15432794310217657 … -0.15389660284267684 -0.15432794310325254; … ; -0.15404229468771546 -0.15389660283984227 … -0.15346128702343226 -0.15389660284093112; -0.1544722989106403 -0.1543279431013607 … -0.15389660284185655 -0.1543279431024383;;; -0.15418753101149862 -0.1540422946879947 … -0.15360833736055898 -0.15404229468907757; -0.15404229468894556 -0.1538966028410779 … -0.153461287024674 -0.1538966028421633; … ; -0.153608337358669 -0.15346128702181525 … -0.15302191314811353 -0.15346128702291378; -0.15404229468812672 -0.1538966028402556 … -0.15346128702384695 -0.15389660284134285;;; … ;;; -0.153478805867214 -0.15333139197651818 … -0.15289092494117187 -0.1533313919776211; -0.1533313919774861 -0.15318350575482603 … -0.15274162879332437 -0.153183505755932; … ; -0.15289092493924777 -0.15274162879041375 … -0.1522955396034788 -0.15274162879153233; -0.15333139197665321 -0.15318350575398992 … -0.15274162879248238 -0.15318350575509718;;; -0.15418753101143337 -0.15404229468792777 … -0.15360833736049528 -0.15404229468901331; -0.15404229468888084 -0.15389660284101125 … -0.15346128702461112 -0.15389660284209974; … ; -0.15360833735860085 -0.15346128702174575 … -0.15302191314804717 -0.15346128702284662; -0.15404229468806038 -0.1538966028401877 … -0.15346128702378214 -0.15389660284127746;;; -0.15461620271177104 -0.1544722989104767 … -0.15404229468955752 -0.1544722989115511; -0.1544722989114205 -0.15432794310214304 … -0.1538966028426452 -0.1543279431032205; … ; -0.1540422946876814 -0.15389660283980758 … -0.15346128702339906 -0.15389660284089757; -0.15447229891060743 -0.15432794310132705 … -0.15389660284182405 -0.1543279431024057]), 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.015163896764803859 -0.012550128167764654 … -0.0049668670504268075 -0.012550128226448888; -0.012550128166154039 -0.010362594044423738 … -0.004135157623243734 -0.010362594084379038; … ; -0.004966867127397037 -0.004135157734725494 … -0.0020697784841564237 -0.004135157789445986; -0.012550128228563377 -0.010362594096358764 … -0.004135157718413482 -0.010362594154883775;;; -0.012550128218538485 -0.010362594088263628 … -0.0041351576869918015 -0.010362594139918888; -0.010362594057451179 -0.00862203789104819 … -0.0038975099204877194 -0.008622037940367006; … ; -0.004135157808835711 -0.0038975101138846724 … -0.0034841056038599633 -0.003897510133055318; -0.010362594168307294 -0.008622038012464611 … -0.0038975100640603684 -0.00862203805142151;;; -0.004966867074845692 -0.004135157666831908 … -0.0020697784791744367 -0.004135157755693555; -0.004135157586868401 -0.003897509850913097 … -0.003484105491902523 -0.0038975099659881688; … ; -0.002069778534734037 -0.003484105618125984 … -0.005620578881644318 -0.003484105619006244; -0.0041351577939833795 -0.0038975100917481163 … -0.0034841056152598027 -0.003897510131608481;;; … ;;; -0.006547898094058736 -0.0048763404049400275 … -0.001486642900389162 -0.004876340411927338; -0.004876340460708797 -0.0026779894566212265 … -0.004123053418484977 -0.002677989486742156; … ; -0.0014866428752781163 -0.0041230534261586585 … -0.0071708944363810334 -0.004123053415816005; -0.004876340422131535 -0.002677989445201732 … -0.004123053400216202 -0.002677989452159098;;; -0.004966867091837561 -0.004135157690514368 … -0.00206977861508636 -0.004135157797916562; -0.004135157709113758 -0.003897510020692047 … -0.003484105657080583 -0.0038975100879726247; … ; -0.0020697784766563476 -0.00348410553137694 … -0.00562057885677268 -0.00348410559992857; -0.00413515771965261 -0.003897509988603786 … -0.003484105613620333 -0.0038975100919131866;;; -0.012550128186139368 -0.010362594048977748 … -0.004135157738339333 -0.010362594130067723; -0.010362594070119118 -0.008622037917310203 … -0.0038975100285193985 -0.008622037971805638; … ; -0.004135157702196374 -0.0038975099742122373 … -0.0034841055475466056 -0.003897510060170789; -0.010362594093181706 -0.00862203791451709 … -0.0038975100395534554 -0.008622037999852901])], 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.009532891532398603 -0.006919174305209826 … 0.0006640022997178941 -0.006919174364964316; -0.0069191743045395986 -0.0047316775476369186 … 0.0014956983527523457 -0.004731677588665526; … ; 0.0006640022246171412 0.0014956982440979347 … 0.003561048108689715 0.001495698188291478; -0.006919174366138085 -0.00473167759875829 … 0.0014956982584009432 -0.004731677658357941;;; -0.0069191743565341895 -0.004731677591086038 … 0.0014956982893982963 -0.004731677643814495; -0.004731677561216445 -0.0029911464905670697 … 0.0017333343803677837 -0.0029911465409618587; … ; 0.0014956981694285547 0.0017333341898053412 … 0.0021467168232341544 0.0017333341695458444; -0.004731677671259739 -0.002991146611167589 … 0.0017333342376155063 -0.0029911466512020975;;; 0.0006640022763339873 0.0014956983111530254 … 0.003561048112830403 0.0014956982212085222; 0.0014956983901657653 0.0017333344515413772 … 0.0021467169339498935 0.0017333343353808958; … ; 0.0035610480591608184 0.002146716810585118 … 1.0230389197517925e-5 0.00214671680860632; 0.0014956981838697148 0.0017333342115286444 … 0.002146716811419702 0.0017333341705810384;;; … ;;; -0.0009170700357030506 0.0007544781983501499 … 0.00414415945108515 0.0007544781902599166; 0.0007544781416135155 0.002952823992125964 … 0.001507746177333076 0.002952823960899058; … ; 0.004144159478120323 0.0015077461725699556 … -0.0015401015173278432 0.001507746181794032; 0.0007544781810236946 0.0029528240043816233 … 0.0015077461964439279 0.002952823996317004;;; 0.0006640022594073718 0.0014956982875374836 … 0.0035610479769821787 0.0014956981790497413; 0.0014956982679851342 0.0017333342818290683 … 0.0021467167688347 0.0017333342134600001; … ; 0.00356104811730662 0.0021467168974036617 … 1.0230414135491486e-5 0.0021467168277511627; 0.0014956982582667927 0.0017333343147408926 … 0.0021467168131239812 0.001733334210341725;;; -0.006919174324102544 -0.0047316775517667956 … 0.001495698238082795 -0.004731677633931133; -0.0047316775738517995 -0.0029911465167955266 … 0.0017333342723677737 -0.002991146572368434; … ; 0.0014956982761019193 0.0017333343295124707 … 0.00214671687958068 0.00173333424246393; -0.004731677596101288 -0.0029911465131864275 … 0.001733334262154921 -0.002991146599600876]), 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.4910459250445243e-6 - 1.794024580103606e-6im -0.0020546787816517822 - 0.0046129326767872575im … -0.0013544828641458924 - 0.0031877791554452994im 0.002011817454710546 + 0.00465009197805859im; -0.011007548509196212 - 0.012474498767394775im 0.011655551934265169 + 0.01452209491994011im … -0.005662603272718448 - 0.00524721093612796im 0.00927225963872892 + 0.009094781499432502im; … ; -0.013944406760209768 - 0.015818255245796324im 0.0119676601698963 + 0.013112149332979373im … -0.008540007487415315 - 0.010228843971576019im 0.012786096318680306 + 0.01497498912221483im; 0.011071513287532083 + 0.012451676514095004im -0.009275454832753622 - 0.00910846304162913im … 0.0075178387714213805 + 0.00963324990500961im -0.011658893067896942 - 0.014535851886101459im;;; 0.0005662601688911331 - 0.004702050512800067im 0.0008730613467397136 + 0.005454112288798267im … 0.0010419517759643955 + 0.0012480424868150117im -0.001510139403125465 + 2.6703816573597048e-5im; 0.010115869106831469 + 0.014575199526008716im -0.009691094963708522 - 0.013826141608280485im … 0.004798079202107716 + 0.005398473916151739im -0.008156460215402422 - 0.010332000494115269im; … ; 0.012498763883677136 + 0.013093298841576697im -0.01050757087086342 - 0.010706726304253025im … 0.007495022607935614 + 0.008470752292401186im -0.011149991456805727 - 0.012168992091877872im; -0.010815282047733205 - 0.009055344161883496im 0.00858096403307267 + 0.006802543158515217im … -0.006357463700828566 - 0.00719592405204023im 0.010115647737270014 + 0.01029671597391859im;;; -0.00036804748327143983 + 0.0031894220732951034im -0.00016355271763746903 - 0.003112577432417294im … -0.0005632381350265075 - 6.9491048110069845e-6im 0.0006547419449147248 - 0.0012497559282545873im; -0.006337141668055483 - 0.009626682727087839im 0.0057171956191170645 + 0.008652711140005938im … -0.003113697456690317 - 0.004035995356232698im 0.005074852102580922 + 0.0071904639315509im; … ; -0.008217952879146978 - 0.00812420815227145im 0.006985394334729273 + 0.006715465035556093im … -0.005230338087694841 - 0.005581940447694157im 0.00735464372320841 + 0.007555886724151383im; 0.006843294359183549 + 0.005253757451169386im -0.0054383485774983364 - 0.0039416603225203755im … 0.004052057779719726 + 0.0040540493681223075im -0.006080691997220998 - 0.005403909896025489im;;; … ;;; -0.0002191364764583167 + 0.0014524916987430758im 0.0003492372303290969 - 0.0008637239979353963im … 9.103881470658696e-5 + 0.001330869063236704im 2.7665398892212542e-5 - 0.00159573906197469im; 0.0033241070337279104 + 0.002311035436560609im -0.0030829388528189495 - 0.0023928819646513973im … 0.0019268614513788863 + 0.0010737983105321021im -0.002803854592573358 - 0.0017576463976281516im; … ; 0.004292136742819556 + 0.005817951712495471im -0.0037764869150250582 - 0.00502279776827043im … 0.0030874268174014177 + 0.004314328478365307im -0.0039729469413338735 - 0.005469738537009555im; -0.0029481251539885027 - 0.004770187463494183im 0.002473885805170683 + 0.0038807009657749448im … -0.002086352466327686 - 0.0034571162653964857im 0.0027529714323608545 + 0.00451592009493524im;;; 0.000469414228172406 - 0.003259979304533156im -0.0007390270391475918 + 0.0013156225142305442im … -0.00028648087042568106 - 0.0021344271175564424im 7.927605713675941e-5 + 0.0031785093417071334im; -0.0069126149543601075 - 0.005198427677746647im 0.006146912474346616 + 0.005347582018691082im … -0.0033456635170879246 - 0.00225293749120157im 0.005504563386910105 + 0.00388531234288462im; … ; -0.007941202173365487 - 0.010251665132600393im 0.00675344764001348 + 0.008498492339339709im … -0.005079036183078743 - 0.006745429790445832im 0.007122673648663376 + 0.009338870924848485im; 0.006267838971432782 + 0.009681979907465034im -0.005008639174693818 - 0.007246777298353047im … 0.0038201274614238128 + 0.005837074386009192im -0.005650990520274377 - 0.00870904762368536im;;; -0.0006610529751135909 + 0.004755702452773691im 0.0015929091335310395 - 8.719418385136293e-5im … 0.000794489079737239 + 0.003149940646092333im -0.0007903483409558291 - 0.005514684995476559im; 0.01083567022913854 + 0.00903386280769888im -0.01015484278755811 - 0.010258708119088493im … 0.004992353774428699 + 0.0039055193502837366im -0.008620203319582322 - 0.006764526214262304im; … ; 0.012251289277317574 + 0.014995151643001275im -0.010313311595579555 - 0.012199641912574876im … 0.00738332998293332 + 0.009328838525567106im -0.010955720847647366 - 0.013661899679637685im; -0.010095464894659683 - 0.014596685682476905im 0.00811718975167612 + 0.010370001666530187im … -0.006163189770880883 - 0.008688849082790214im 0.00965190919993469 + 0.013864208256195469im],)])]), [8.519871707324948e-7 4.6272636967099523e-7 … 2.4177761006633075e-8 4.627263766360816e-7; 4.627263694798359e-7 2.501446676668972e-7 … 1.3582980334672564e-8 2.501446707571937e-7; … ; 2.417776218765559e-8 1.3582981484785027e-8 … 1.565917386909722e-9 1.3582982049312315e-8; 4.6272637688703836e-7 2.501446716837471e-7 … 1.358298131649618e-8 2.501446762102682e-7;;; 4.627263756972045e-7 2.5014467105763804e-7 … 1.3582980992332885e-8 2.501446750528369e-7; 2.501446686744854e-7 1.3895601454905982e-7 … 1.1279973945992647e-8 1.389560170847411e-7; … ; 1.3582982249348221e-8 1.1279975701958351e-8 … 7.93774337128655e-9 1.1279975876018919e-8; 2.501446772484867e-7 1.3895602079157678e-7 … 1.1279975249571099e-8 1.3895602279450837e-7;;; 2.4177761381312054e-8 1.3582980784351412e-8 … 1.5659173752125398e-9 1.3582981701103364e-8; 1.3582979959399805e-8 1.1279973314282567e-8 … 7.937742572732193e-9 1.1279974359116273e-8; … ; 1.5659175056788207e-9 7.937743473040078e-9 … 3.5718535164964555e-8 7.937743479322411e-9; 1.3582982096120745e-8 1.1279975500967286e-8 … 7.937743452598136e-9 1.1279975862883024e-8;;; … ;;; 5.791842911290629e-8 2.2816032541900258e-8 … 5.620355214002189e-10 2.281603264491464e-8; 2.281603336409534e-8 3.489281982450257e-9 … 1.3458500528087926e-8 3.4892821047617355e-9; … ; 5.620354920723096e-10 1.3458500606757681e-8 … 7.72831205290164e-8 1.3458500500728673e-8; 2.2816032795350796e-8 3.489281936079489e-9 … 1.3458500340806368e-8 3.489281964333312e-9;;; 2.4177761642035327e-8 1.358298102867622e-8 … 1.565917694364477e-9 1.3582982136699778e-8; 1.3582981220559044e-8 1.1279974855802831e-8 … 7.937743750892056e-9 1.1279975466683254e-8; … ; 1.565917369297651e-9 7.937742854291384e-9 … 3.571853466547925e-8 7.93774334324606e-9; 1.3582981329279463e-8 1.1279974564457766e-8 … 7.937743440906117e-9 1.1279975502465744e-8;;; 4.627263718518504e-7 2.5014466801913197e-7 … 1.3582981522066835e-8 2.5014467429091167e-7; 2.501446696542664e-7 1.389560158992989e-7 … 1.127997492687358e-8 1.3895601870112507e-7; … ; 1.3582981149193442e-8 1.127997443378682e-8 … 7.937742969621252e-9 1.1279975214255584e-8; 2.501446714380085e-7 1.3895601575569566e-7 … 1.1279975027061554e-8 1.3895602014315903e-7;;;;], Matrix{ComplexF64}[[0.1950128963475652 - 0.17273372668977574im 0.7154804470010004 - 0.6476629850099624im -2.7092986439601672e-6 - 4.547363384418189e-7im 4.230805747267199e-6 - 1.929277327005261e-6im; -0.1323121363665184 + 0.11719618982489477im 0.031526885404762624 - 0.028538577284238666im 0.03961347336513628 - 0.24051887491579235im -0.464948434751489 - 0.355282529086912im; … ; 0.04203655445915677 - 0.03723410525795472im -0.011402863936730685 + 0.010322033109662648im -0.0018963914732424243 - 0.0036927563140412845im -0.006881788862272542 - 0.004537954143719923im; -0.07537183004952629 + 0.06676100576004217im 0.022589949834022804 - 0.020448740763796564im 0.0126659240950111 + 0.007636630097880789im 0.013285578239071536 + 0.00618786798655781im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944082676792, [[-0.5541104233724997, -0.010678393162858714, 0.1449565955408893, 0.14495679606198464]], 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.6206926214952306e-7, -6.293614328137214e-7, -2.2951195984500963e-7, 4.994026231299731e-7, 7.923311649999066e-7, 5.851071649400082e-7, 6.516971465195703e-7, 1.018931060280713e-6, 3.6385935769288854e-7, -2.201093444426784e-6  …  2.781018579136238e-6, 1.201695646754714e-6, 2.4456301976076343e-8, 1.0078535591772482e-7, 5.535307815328355e-7, 4.534173060317613e-7, 4.631908793337421e-9, -1.0437058895987745e-7, 1.12110891356485e-7, 1.2550860890906588e-7], resid_history = [0.24939209303134624, 0.003766498634357822, 0.00028515516768058577, 4.682253246613109e-6, 4.2116942507398584e-8, 6.759355310585287e-11, 4.085756751510658e-8, 1.0805512268590958e-9], converged = true, n_iter = 8, residual_norm = 1.0805512268590958e-9, maxiter = 100, tol = 1.0e-8, s = 0.935430828276514, residual = [4.892175756242443e-7, 3.2255555587600885e-7, 6.372694717360468e-8, -8.102385547881673e-8, -7.71151524043138e-8, -3.3992192317372796e-8, -2.0476618651151164e-8, -1.1120134579373909e-8, 1.8593681596542696e-9, -3.653950304762912e-8  …  6.838875885631663e-7, 2.982868101380664e-7, 6.229710397963413e-9, 2.6682588277110216e-8, 1.53526609339652e-7, 1.3410139198748997e-7, 1.5297497138245971e-9, -4.233158291926219e-8, 6.793760295633816e-8, 2.1239389406396338e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.782532536239574e-7, -4.060927076472382e-8, 1.0805512268590958e-9, -0.00020437430744746764, 4.687720003748128e-6, -4.260319684116641e-8, 6.759355310585287e-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}}([[-2.746565972010076e-8, 9.393828019451606e-7, 3.4205305380195235e-7, -5.71433993327967e-7, -6.387178657073302e-7, -3.253827346853226e-7, -2.2993842765228564e-7, -1.5099797868294386e-7, 1.7070962659716176e-8, -3.88472958840584e-7  …  1.6150402442758265e-7, -8.901474937348896e-8, -6.655288490112794e-9, -3.8349101946933843e-8, -1.9377733333555326e-7, -1.4207053828969154e-7, -2.037719580072281e-9, 8.250675558442523e-8, -1.508312863770193e-7, -3.352278927376836e-7], [3.584655518961875e-7, -1.3347101286977557e-6, -5.10565960469374e-7, 8.521819174725657e-7, 9.40826842052847e-7, 4.666612449341984e-7, 3.045549750799161e-7, 1.3067890375877888e-7, -6.803106295836787e-8, 8.468928883705709e-7  …  -2.3047466777179772e-8, 2.871748103737639e-7, 1.4766379078914756e-8, 8.14566572980949e-8, 4.2655235663850586e-7, 3.3148624871995784e-7, 4.5808643330236155e-9, -1.7224787003145703e-7, 3.0383071622435604e-7, 6.970777226103159e-7]]), H = [1.183305802697928 -0.00033148458338621905 … 0.0 0.0; 0.10846537316242187 1.0095951424612053 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1882665357039435 0.09182608695982246 … 0.0 0.0; 0.0 1.0057623713664414 … 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.9257125186089437
Interacting polarizability:     1.7736548496648366

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 = 1.65e-08
[ Info: GMRES linsolve in iteration 2; step 1: normres = 1.20e-09
[ Info: GMRES linsolve in iteration 2; step 2: normres = 3.47e-11
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 4.16e-12
└ * number of operations = 12
Non-interacting polarizability: 1.925712518608875
Interacting polarizability:     1.77365484813328

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