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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457110688041708
Displaced dipole:  0.017612220597271775
Polarizability :   1.7746791704152192

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.0001300622756650903 + 0.003293323908661281im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.309042086291714 + 0.01480478858982288im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.04102630145178337 + 0.0017858979906723188im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.0808867945134239 - 0.002884139330803384im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620691009248387e-7 -2.502013222481772e-7 … -4.88063242951367e-8 -2.5020118415171907e-7; -6.293614327936376e-7 -4.793035550831579e-7 … -1.207783509925278e-7 -4.793034152689834e-7; … ; 1.3626266062955704e-7 1.1356234927017055e-7 … 4.7119881464455594e-8 1.1356233707607028e-7; 1.330550263468599e-7 1.3632657693335648e-7 … 5.200691155717496e-8 1.3632665226475522e-7;;; -2.502012452828285e-7 -1.7404893021476063e-7 … -3.556486619376001e-8 -1.7404883075367323e-7; -4.793034819925475e-7 -3.67973152843101e-7 … -1.1202888088111831e-7 -3.6797304689538494e-7; … ; 1.1356237848082694e-7 1.1211105899785133e-7 … 1.0980691502821735e-7 1.1211107290706948e-7; 1.3632664281384627e-7 1.255088225368567e-7 … 5.094916094665801e-8 1.255088824230281e-7;;; -4.880636102345668e-8 -3.556490994811674e-8 … -1.1412121217284202e-8 -3.556487846812743e-8; -1.2077837286464612e-7 -1.1202893755481506e-7 … -9.786216973095213e-8 -1.1202886653535568e-7; … ; 4.7119983396564725e-8 1.0980693703945832e-7 … 2.493825620694407e-7 1.098069786484083e-7; 5.2006939835043005e-8 5.094916257088871e-8 … 4.910636584558315e-8 5.09491798970489e-8;;; … ;;; 6.475757123762854e-8 4.110334515933093e-8 … -6.6816215574639485e-9 4.110332492045757e-8; 1.6897590448881023e-7 6.653315068664922e-8 … -1.3294414231080083e-7 6.65332041082383e-8; … ; 3.26774181427356e-8 1.6017355497622626e-7 … 3.847037098894135e-7 1.6017354435503696e-7; -9.142015045759385e-8 -3.583048997793581e-8 … 7.037283210047709e-8 -3.583051635546038e-8;;; -4.8806322284256664e-8 -3.5564880301981456e-8 … -1.1412096492881685e-8 -3.5564844929034066e-8; -1.2077822899571156e-7 -1.1202882025845397e-7 … -9.786205412748124e-8 -1.1202873262976628e-7; … ; 4.7119835350974466e-8 1.0980682149193036e-7 … 2.4938236706318795e-7 1.0980680477772042e-7; 5.2006871296964595e-8 5.094910491757461e-8 … 4.910629931040298e-8 5.094911107933759e-8;;; -2.502012315106286e-7 -1.740489267676595e-7 … -3.556485383585957e-8 -1.7404881582003662e-7; -4.793034059147636e-7 -3.679731073475699e-7 … -1.1202882318129858e-7 -3.679729730084767e-7; … ; 1.1356231129800179e-7 1.1211101777107356e-7 … 1.0980683281624093e-7 1.1211099518819681e-7; 1.3632659551256604e-7 1.2550878846489934e-7 … 5.094913047834061e-8 1.2550883990676604e-7;;;;], δHtotψ = Matrix{ComplexF64}[[7.565561068790795e-5 - 0.0019156844595492612im 0.025125808681902943 + 0.09192504635405856im 0.6296935285572467 - 0.03107982361071445im -1.400054424910947 + 0.21364689205707588im; -0.19806799405030848 - 0.009625913368240925im 1.5916741325971944 - 0.33652449708187687im 0.1415180574011218 - 0.004168504932632261im -0.3082507681929349 + 0.04915516097433612im; … ; -0.048967585839442045 - 0.0022646678502760585im -0.019940713885241993 + 0.005504271919374539im -0.005962554671201801 - 0.005771007229611854im 0.016587167803069174 - 0.011278766631552826im; 0.07232835312623873 + 0.002583801850986982im 0.024700556127235106 - 0.006691921412758699im 0.06129357291358831 - 0.05509638963642876im -0.02239966887050532 - 0.10085809522552716im]], δVind = [0.00196542947198519 0.002080998303990544 … 0.0031536579163285333 0.0020809972823198036; 0.0038891271267861233 0.004779265542469399 … 0.010276661271505227 0.004779264469538862; … ; -0.006150849647709454 -0.008269603855072732 … -0.017302980336438383 -0.008269606505345948; 0.0002386803719073207 -0.000398698958176791 … -0.003664701062163653 -0.00039870004060813436;;; 0.0020809976952467432 0.00222322548784354 … 0.003420111495214895 0.002223224373123978; 0.004779265125868168 0.005734949063741287 … 0.010903587274581766 0.005734948019279725; … ; -0.008269604921749057 -0.009601008488224755 … -0.012623115436232558 -0.009601010339088912; -0.00039869975518912605 -0.0010730150358832004 … -0.004230686982807206 -0.0010730162591804128;;; 0.003153658954845148 0.003420113588922626 … 0.004832492249222735 0.0034201119988852622; 0.01027666361787658 0.010903589204851329 … 0.012272131738711464 0.010903589360775337; … ; -0.017302964258778204 -0.012623111526426003 … -0.009620267959561937 -0.012623112125079775; -0.0036646989913733984 -0.004230683749998675 … -0.00549111476758536 -0.00423068578885161;;; … ;;; -0.0023365538153419774 -0.0028154562175054863 … 0.005693320063955787 -0.002815453815414405; -0.012913073691856288 -0.017841142500943083 … 0.01149833866529316 -0.017841132841068957; … ; -0.025174533964421018 -0.012810892373207594 … -0.008395878281601887 -0.012810895001128274; 0.007597978236457642 0.010225406256908478 … -0.005449741876196515 0.01022540722105768;;; 0.0031536600429990028 0.003420115119581586 … 0.004832496113759609 0.0034201132418067983; 0.01027666736624342 0.010903593156487753 … 0.012272135738821358 0.010903593011026672; … ; -0.017302981945472605 -0.012623113849631174 … -0.009620268713227607 -0.012623117151241928; -0.003664699446413653 -0.004230684051612966 … -0.005491116011773451 -0.004230686485438451;;; 0.0020809977847192616 0.002223225644735936 … 0.0034201119228713796 0.00222322445401617; 0.004779265052029054 0.005734949112908379 … 0.010903588595533016 0.0057349479432035715; … ; -0.008269605472003 -0.009601008373247594 … -0.012623117737937891 -0.009601011552035927; -0.0003986993688654312 -0.001073014611935742 … -0.004230687164261296 -0.0010730158727558963;;;;], δρ0 = [-3.602137173660373e-7 -2.4888407002573475e-7 … -4.8528286310340896e-8 -2.488840108338602e-7; -5.97224859940837e-7 -4.606977616350681e-7 … -1.1905762515533071e-7 -4.6069768095406344e-7; … ; 1.273579719079709e-7 1.0836346876141815e-7 … 4.640716665824974e-8 1.0836344477046291e-7; 1.0352037301733598e-7 1.1957984433123272e-7 … 5.069511453218586e-8 1.1957986299739775e-7;;; -2.4888405148455967e-7 -1.7310658159916922e-7 … -3.535692287535654e-8 -1.7310653988188728e-7; -4.6069773227689784e-7 -3.5721705623521367e-7 … -1.1091435987298175e-7 -3.5721699407608004e-7; … ; 1.0836348764476317e-7 1.084650661176294e-7 … 1.0862245377420846e-7 1.0846506814961448e-7; 1.1957986729868251e-7 1.1611166533631036e-7 … 5.0208373447288055e-8 1.1611168271618451e-7;;; -4.852831380389102e-8 -3.535694988652972e-8 … -1.1341296920008362e-8 -3.535693173508486e-8; -1.1905764017376795e-7 -1.1091440114432316e-7 … -9.770663737191484e-8 -1.109143426653816e-7; … ; 4.640726931545327e-8 1.0862248791920829e-7 … 2.48452710524482e-7 1.086225192468339e-7; 5.069514868849957e-8 5.0208389839941445e-8 … 4.926866700457753e-8 5.020839513766277e-8;;; … ;;; 6.434565435191809e-8 4.083861590247206e-8 … -6.63785916444786e-9 4.0838612487698e-8; 1.6975068052131522e-7 6.683113636685423e-8 … -1.3340447609938256e-7 6.683119670990057e-8; … ; 3.270977438466902e-8 1.6023514505039018e-7 … 3.839015372540291e-7 1.6023512162886004e-7; -9.271391636040785e-8 -3.633636238021544e-8 … 7.12630807108567e-8 -3.6336382429951765e-8;;; -4.8528258321380194e-8 -3.5356907795937634e-8 … -1.1341268113524985e-8 -3.535688573747425e-8; -1.1905748548370545e-7 -1.1091427343614242e-7 … -9.770651247078718e-8 -1.1091419848253245e-7; … ; 4.6407127583150285e-8 1.0862238285177296e-7 … 2.4845253628355875e-7 1.086223565521428e-7; 5.069509406406028e-8 5.020834420335209e-8 … 4.9268609834804565e-8 5.0208338575058805e-8;;; -2.4888399994679807e-7 -1.7310655065440058e-7 … -3.535690428150309e-8 -1.7310649740240242e-7; -4.606976303750973e-7 -3.5721699123673963e-7 … -1.1091429698613059e-7 -3.5721690118129877e-7; … ; 1.083634291557871e-7 1.0846503152583948e-7 … 1.0862237713589864e-7 1.0846499826310043e-7; 1.1957984895650838e-7 1.1611165217842968e-7 … 5.0208349068925484e-8 1.1611166183171457e-7;;;;], δeigenvalues = [[0.0004950021205721699, 0.09679852747209985, 0.04583339399446872, 0.03408697520749402]], δ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.15475965285824703 -0.15461620273258836 … -0.15418753103541474 -0.15461620273465718; -0.1546162027357511 -0.15447229893399028 … -0.15404229471461817 -0.15447229893606193; … ; -0.1541875310291283 -0.15404229470511482 … -0.1536083373791564 -0.1540422947072054; -0.15461620273149523 -0.15447229892971978 … -0.1540422947103123 -0.15447229893179376;;; -0.15461620273349297 -0.15447229893173547 … -0.15404229471232042 -0.15447229893378941; -0.15447229893489148 -0.1543279431251478 … -0.15389660286715726 -0.15432794312720483; … ; -0.1540422947060397 -0.15389660285766568 … -0.15346128704270653 -0.1538966028597431; -0.15447229893063116 -0.15432794312087295 … -0.153896602862848 -0.15432794312293246;;; -0.15418753103314278 -0.15404229470916336 … -0.1536083373831862 -0.15404229471121386; -0.15404229471233044 -0.15389660286398885 … -0.15346128704905446 -0.15389660286604281; … ; -0.1536083373768759 -0.15346128703951958 … -0.1530219131672378 -0.15346128704159573; -0.1540422947080423 -0.15389660285968504 … -0.15346128704471595 -0.1538966028617424;;; … ;;; -0.15347880588972856 -0.15333139199851278 … -0.1528909249647689 -0.15333139200068016; -0.15333139200179696 -0.1531835057786173 … -0.152741628818729 -0.15318350578079007; … ; -0.15289092495825352 -0.15274162880886621 … -0.15229553962348968 -0.15274162881106126; -0.153331391997402 -0.15318350577420506 … -0.15274162881427827 -0.15318350577638057;;; -0.15418753103371458 -0.15404229470969763 … -0.15360833738382812 -0.1540422947118219; -0.1540422947129273 -0.15389660286454782 … -0.1534612870497231 -0.15389660286667683; … ; -0.15360833737742088 -0.1534612870400289 … -0.1530219131678564 -0.15346128704217804; -0.1540422947085975 -0.15389660286020346 … -0.15346128704534293 -0.15389660286233442;;; -0.15461620273378887 -0.15447229893200917 … -0.1540422947126454 -0.1544722989341016; -0.15447229893519834 -0.15432794312543277 … -0.1538966028674953 -0.15432794312752884; … ; -0.1540422947063157 -0.15389660285792264 … -0.1534612870430175 -0.15389660286003726; -0.15447229893091644 -0.1543279431211371 … -0.15389660286316478 -0.15432794312323495]), 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.015163900551401144 -0.012550131920650475 … -0.004966868614911828 -0.012550131539080519; -0.012550131791118344 -0.010362597889206541 … -0.004135159674646105 -0.01036259733659619; … ; -0.004966868988482503 -0.004135160165686667 … -0.002069779785716749 -0.0041351594996493065; -0.012550131476982883 -0.010362597313839777 … -0.004135159108160475 -0.010362597033941452;;; -0.012550131810131806 -0.01036259773335901 … -0.004135159761743356 -0.010362597375321301; -0.010362597491351287 -0.008622041641724126 … -0.0038975123173394417 -0.00862204108775414; … ; -0.004135160473040334 -0.0038975129840104327 … -0.003484108218610476 -0.0038975129377305045; -0.010362597464082308 -0.00862204127304365 … -0.003897512404613672 -0.008622041183504047;;; -0.004966869589023598 -0.0041351609554781225 … -0.002069780858258685 -0.004135160142757069; -0.004135159637586658 -0.0038975129492544724 … -0.003484107236187389 -0.003897511780171601; … ; -0.002069782681460044 -0.0034841089550633575 … -0.0056205825160234625 -0.003484109506702279; -0.004135160935463388 -0.0038975135028717594 … -0.0034841088125883626 -0.003897513461592234;;; … ;;; -0.0065479045551425415 -0.0048763459929381785 … -0.0014866364426558535 -0.004876346823774738; -0.004876348100355551 -0.002677997004063403 … -0.004123049655779689 -0.0026779987524263023; … ; -0.001486638027245987 -0.00412305233461336 … -0.007170895128988122 -0.004123051887280039; -0.004876345803992769 -0.002677994562093325 … -0.004123051798814892 -0.002677995348817551;;; -0.004966866994693167 -0.004135158511527783 … -0.002069777965969636 -0.00413515764910193; -0.004135156744662625 -0.0038975104780838788 … -0.0034841048486152354 -0.003897509085166373; … ; -0.0020697787512065913 -0.003484107003840805 … -0.005620580343715435 -0.003484106409620014; -0.004135158238143536 -0.003897511358896011 … -0.0034841063218317304 -0.003897510794222273;;; -0.012550131269777345 -0.010362597320764881 … -0.004135158860057234 -0.010362596812982908; -0.010362596809936285 -0.008622041134323907 … -0.003897511193019936 -0.008622040348010142; … ; -0.004135159244381525 -0.0038975123357980097 … -0.0034841067620695085 -0.0038975115066825217; -0.010362596864750152 -0.008622040889077014 … -0.0038975112644773872 -0.008622040477511125])], 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.009532895340793563 -0.0069191780794309425 … 0.0006640007123517597 -0.006919177699929776; -0.00691917795306142 -0.004731681415509002 … 0.0014956962767155137 -0.004731680864970299; … ; 0.0006640003450676122 0.0014956957951782106 … 0.0035610467876906614 0.0014956964591249936; -0.00691917763467003 -0.004731680835871709 … 0.0014956968475070328 -0.004731680558047364;;; -0.006919177969816911 -0.00473168125740683 … 0.001495696191915869 -0.004731680901423034; -0.004731681018554944 -0.002991150264214242 … 0.0017333319590356437 -0.0029911497123012772; … ; 0.0014956954868997003 0.0017333313018561714 … 0.002146714189209365 0.001733331346058678; -0.004731680987025623 -0.002991149891258881 … 0.001733331876070744 -0.0029911498037787972;;; 0.0006639997405119228 0.0014956950013381605 … 0.003561045711118949 0.0014956958120087155; 0.0014956963160626304 0.0017333313302890507 … 0.0021467151652845577 0.001733332497317954; … ; 0.0035610438942279034 0.0021467134559434065 … 1.0226735694115335e-5 0.00214671290222834; 0.0014956950224741371 0.001733330780975571 … 0.0021467135932221413 0.0017333308201977421;;; … ;;; -0.0009170765193014293 0.0007544725883573983 … 0.004144165885221446 0.0007544717553534617; 0.0007544704776559028 0.002952816420892513 … 0.0015077499146337412 0.002952814670356852; … ; 0.004144164307146711 0.0015077472456627924 … -0.001540102229945814 0.0015077476908010633; 0.0007544727784136694 0.002952818867274895 … 0.0015077477760493396 0.002952818078375159;;; 0.0006640023342705612 0.0014956974447542055 … 0.003561048602766067 0.0014956983050557855; 0.0014956992083898077 0.0017333338009006749 … 0.0021467175521880794 0.0017333351916891614; … ; 0.0035610478239363474 0.002146715406656644 … 1.0228907383499214e-5 0.002146715998728293; 0.0014956977192387387 0.0017333329244329006 … 0.002146716083351803 0.001733333486975677;;; -0.0069191774297583525 -0.004731680845086399 … 0.0014956970932770014 -0.0047316803393968355; -0.004731680337446807 -0.002991149757098962 … 0.0017333330830171419 -0.002991148972881271; … ; 0.0014956967152824803 0.0017333319498116333 … 0.0021467156454393313 0.0017333327768125073; -0.004731680387978738 -0.0029911495075563947 … 0.0017333330158902545 -0.0029911490980883554]), 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.847876873938608e-7 - 2.7121431032667815e-6im 0.006369491908318879 + 0.0048796877688371525im … 0.004412617479315005 + 0.003273071909419798im -0.006465630125832545 - 0.004837667242841034im; -0.02366480004947328 - 0.0009615608033160235im 0.016153953830467076 - 0.0020641368734447513im … -0.01437201132585457 - 0.0026888587708678597im 0.023674017387306463 + 0.0036294078489466463im; … ; -0.025852459119394974 - 0.001076287545167242im 0.023275218866264068 + 0.001895903384275705im … -0.012355025641086011 + 0.0005045404723187438im 0.020694499667707204 - 5.8260029136644596e-5im; 0.023667916226796155 + 0.0009077488542302589im -0.02362917651674399 - 0.003635738158465921im … 0.008320675962417326 - 0.001792295428786336im -0.016109467658466298 + 0.002057732125187207im;;; -0.0033391477529355045 + 0.0052260967618284445im -0.0017746100018438328 - 0.005925481458389446im … -0.00370310432732232 - 0.0011904560014445328im 0.005745449557788041 - 0.00023198268016641873im; 0.021854202540908228 - 0.0022621332773045105im -0.015279104192153146 + 0.0031626986594338576im … 0.01230513707363167 + 0.0012184180939499556im -0.02012055133411151 - 0.0005027066312899735im; … ; 0.021311352565900197 + 0.0019605789165049684im -0.019154054537700256 - 0.0023692635038687323im … 0.010480198892590185 + 6.924500625272976e-5im -0.017128121093765127 - 0.0008351786148320428im; -0.017929286695242058 - 0.003833619466820045im 0.01756288996569865 + 0.004460268097075356im … -0.00672773967188218 + 0.00046664984258801116im 0.012721829932849128 + 0.0007947156605783282im;;; 0.0022780971719274304 - 0.0035905504372874316im -0.000282855456698254 + 0.003453573630814847im … 0.002207900673168407 - 0.00011494659387926921im -0.0028639643823150165 + 0.0014993954508303819im; -0.012916014443022695 + 0.0020200580663568457im 0.00968015117299376 - 0.0022437100911275373im … -0.007888714225657262 - 0.00013739400540682553im 0.011706499768089842 - 0.0007097363749890046im; … ; -0.013047628752481263 - 0.00175508492133887im 0.01194727097475772 + 0.0019341424110354208im … -0.0070896156484387355 - 0.00040549667147921206im 0.01078272743604005 + 0.0010522488026018764im; 0.009776544633997079 + 0.0029164193121069834im -0.009352097050334403 - 0.0029955898488518732im … 0.0041869817167241525 + 0.00033942805541233836im -0.007326238629316179 - 0.0014614808718574441im;;; … ;;; 0.0010244294450435115 - 0.0015077277128711327im -0.0013701413085439192 + 0.0008839419496079814im … -0.00013305046565502036 - 0.0013646508638842651im -0.00035605837346371885 + 0.0016519500770971374im; 0.004418211066854896 + 0.0014769578729244958im -0.003487300898502288 - 0.0009384461210774814im … 0.0034852979994991033 + 0.001348920691941161im -0.004367407253070835 - 0.0016049974648176062im; … ; 0.008009363190934284 - 0.0005378646672745035im -0.007529256398909586 + 0.00023940181331418795im … 0.004990627102789121 - 0.0006986841335461244im -0.006909975953722149 + 0.0007084758202396284im; -0.006145823720499658 + 0.0010595941528857714im 0.005851796095566392 - 0.0005694783569081173im … -0.003241794899284785 + 0.0010831463719052528im 0.004971643955744927 - 0.0012359779825958268im;;; -0.0023343636547310436 + 0.0034568173842970187im 0.0028841324456860057 - 0.0013849931791696918im … 0.0006858770310318561 + 0.002209512784585361im 0.0003034471960948372 - 0.0033391840361443642im; -0.009748607679262335 - 0.0028187952106937895im 0.007316259196194197 + 0.0013676847240871657im … -0.006612898563267673 - 0.0020856543490056544im 0.009342061102410582 + 0.002901839131452166im; … ; -0.014568887112255865 + 0.0005696559108819429im 0.013222290544363 - 1.4195309739571342e-5im … -0.007921901116347534 + 0.0008657811269753332im 0.012057897241308301 - 0.0008960782098554361im; 0.012943608939457778 - 0.001922270215484505im -0.011716210550525685 + 0.0006158332869384395im … 0.005462532876363097 - 0.0016086091495725128im -0.00969022589150103 + 0.0021498476522729877im;;; 0.0034244108471312503 - 0.0051091561167554275im -0.005737549609458813 + 0.00012961146528512333im … -0.0023425249083611367 - 0.003268501276508736im 0.0017821684168030489 + 0.005823060624070182im; 0.017890909169192292 + 0.0037930149153665464im -0.012727581353199124 - 0.0007353116242777736im … 0.011236814416613623 + 0.0028497393292793343im -0.017568728909845975 - 0.004400833075801608im; … ; 0.022671513376202206 - 0.0001174381588835679im -0.020221593997558907 - 0.0007379646572105252im … 0.01109401767417795 - 0.0008683038426445628im -0.018195818509861447 + 0.000796091604168374im; -0.021892279576605528 + 0.002221508166749978im 0.020114391203812667 + 0.0005621119528528463im … -0.007795636701705957 + 0.0020977876727621293im 0.015273326394642333 - 0.0031033334505290748im],)])]), [8.519878588272769e-7 4.627268150908774e-7 … 2.417778501183767e-8 4.627267698033545e-7; 4.627267997170436e-7 2.501449650360879e-7 … 1.3583001498155801e-8 2.501449222952302e-7; … ; 2.4177790743843035e-8 1.358300656402328e-8 … 1.5659204432721904e-9 1.358299969278236e-8; 4.6272676243315547e-7 2.5014492053516743e-7 … 1.358299565395126e-8 2.501448988868159e-7;;; 4.6272680197371335e-7 2.501449529822843e-7 … 1.3583002396704131e-8 2.5014492529035555e-7; 2.5014493426454066e-7 1.389562073867799e-7 … 1.1279995708440676e-8 1.3895617890488532e-7; … ; 1.3583009734871157e-8 1.1280001761552598e-8 … 7.937762021463414e-9 1.1280001341346718e-8; 2.5014493215546023e-7 1.3895618843138677e-7 … 1.1279996500855306e-8 1.3895618382778163e-7;;; 2.4177799958456775e-8 1.3583014711989443e-8 … 1.5659229618513142e-9 1.3583006327469545e-8; 1.3583001115827545e-8 1.1280001445980334e-8 … 7.937755014152835e-9 1.1279990831169536e-8; … ; 1.5659272431567302e-9 7.937767274349867e-9 … 3.5718608152544014e-8 7.937771209023076e-9; 1.3583014505504533e-8 1.1280006472609457e-8 … 7.937766258123192e-9 1.1280006097807787e-8;;; … ;;; 5.791861030705827e-8 2.2816114925729207e-8 … 5.620279792536664e-10 2.281612717476225e-8; 2.281614599542388e-8 3.489312629815841e-9 … 1.3458461954839364e-8 3.4893197292962264e-9; … ; 5.620298299299906e-10 1.345848941681029e-8 … 7.728314424230052e-8 1.3458484830984829e-8; 2.2816112140101284e-8 3.4893027138593313e-9 … 1.3458483924084837e-8 3.4893059084601243e-9;;; 2.4177760151470746e-8 1.3582989498733262e-8 … 1.5659161700892197e-9 1.3582980601441629e-8; 1.358297127071193e-8 1.1279979008735218e-8 … 7.937737984372661e-9 1.1279966361617168e-8; … ; 1.5659180140058394e-9 7.937753356895988e-9 … 3.571856452706106e-8 7.937749118511135e-9; 1.3582986678340243e-8 1.1279987006157961e-8 … 7.937748492344588e-9 1.1279981879147673e-8;;; 4.627267378404908e-7 2.5014492107077447e-7 … 1.3582993094372777e-8 2.501448817970984e-7; 2.501448815614596e-7 1.3895618129923407e-7 … 1.127998550006833e-8 1.389561408715852e-7; … ; 1.3582997059291207e-8 1.1279995876036521e-8 … 7.937751632421025e-9 1.1279988348000314e-8; 2.5014488580095433e-7 1.3895616869006688e-7 … 1.127998614887279e-8 1.3895614752976383e-7;;;;], Matrix{ComplexF64}[[0.010280257132370955 - 0.2603099809817432im 0.25445113424507304 + 0.9309319154438631im -4.581457300194096e-7 + 1.1539040283817674e-6im 1.5664981640021083e-6 + 5.181284523352124e-6im; -0.006974936980080534 + 0.1766148305562558im 0.01121209797209585 + 0.04102051618413605im -0.012842532311674877 - 0.19612001486297445im 0.06584738472998518 + 0.44448165811076973im; … ; 0.002215989834733825 - 0.056111851700982214im -0.004055278200933457 - 0.014836590095198682im -0.0027987586649585875 - 0.006579563484499017im -0.0038705434153125406 + 0.008546481203129846im; -0.003973284863574507 + 0.10060893451607693im 0.008033816311598971 + 0.029392422271242724im 0.015086472741873976 + 0.025732514811891868im 0.025141227803455634 - 0.023723368265199142im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239440774457514, [[-0.5541104258857869, -0.01067838960336339, 0.14495843385116483, 0.14495918577680109]], 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.620691009248387e-7, -6.293614327936376e-7, -2.29511975071825e-7, 4.994030278404241e-7, 7.923324532964158e-7, 5.851090654475132e-7, 6.516986606471215e-7, 1.0189316672990082e-6, 3.638596019893565e-7, -2.2010926690492566e-6  …  2.781018128136559e-6, 1.2016955964073707e-6, 2.4456482936506525e-8, 1.0078525725889711e-7, 5.535301208076914e-7, 4.5341637191377127e-7, 4.63122929853081e-9, -1.0437080643218767e-7, 1.1211099518819681e-7, 1.2550883990676604e-7], resid_history = [0.24939209326459938, 0.0037665477113969434, 0.0002852656250282847, 4.680875232429533e-6, 5.2820389435238034e-8, 6.125133998173548e-11, 4.3424120281138895e-8, 1.1507357448314428e-9], converged = true, n_iter = 8, residual_norm = 1.1507357448314428e-9, maxiter = 100, tol = 1.0e-8, s = 0.9356154493506436, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.262487398668544e-7, -4.262032186022687e-8, 1.1507357448314428e-9, -0.00020445403807277536, 4.6843132060968455e-6, -5.343419726241531e-8, 6.125133998173548e-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}}([[-4.3007286251137565e-7, -2.0736671764719033e-6, -7.159138418705164e-7, 1.1713305369789087e-6, 1.432454691859507e-6, 1.021954860978648e-6, 1.2332552222180863e-6, 1.8457601107652827e-6, 5.148941621538037e-7, -2.3752999437454126e-6  …  3.084655520712571e-6, 1.2082034979063185e-6, 2.886697966694855e-8, 1.428351035567598e-7, 7.783833711770448e-7, 5.290732627169653e-7, 4.6981616467727254e-9, -1.262734378473572e-7, 2.002659727411946e-7, 3.527020272759043e-7], [8.025590036402861e-7, 2.1153332963368857e-6, 7.041801599227136e-7, -1.1633929925063884e-6, -1.4331038702084643e-6, -1.0024419091356947e-6, -1.1836869344309945e-6, -1.7700092045371e-6, -5.032581974354943e-7, 2.359012213405403e-6  …  -2.267251263367979e-6, -8.588866792523855e-7, -2.070484194928812e-8, -1.0455991613207412e-7, -5.647254677554579e-7, -3.6047203417740403e-7, -2.845063959249351e-9, 7.121337846385389e-8, -1.0686048114537226e-7, -8.832151493564542e-8]]), H = [1.112810736945858 0.013126650568420738 … 0.0 0.0; 0.1342722330549962 1.0271126567076048 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1208821386887913 0.13607156581714566 … 0.0 0.0; 0.0 1.0185017030779882 … 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.9257125000894753
Interacting polarizability:     1.7736548327113497

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.09e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.10e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
│ * norm of residual = 2.77e-12
└ * number of operations = 11
Non-interacting polarizability: 1.9257125000893618
Interacting polarizability:     1.7736548260219025

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