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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457465073954994
Displaced dipole:  0.017612222381127928
Polarizability :   1.7746797031867478

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.0007127342072702089 + 0.0032179060604667167im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.30150317189688436 + 0.06944086712835165im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.040057354245743096 + 0.009041859406477136im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.0790895194887945 - 0.01719999574108593im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620692457147452e-7 -2.502013383645491e-7 … -4.880640398397794e-8 -2.5020139103565894e-7; -6.293614724885022e-7 -4.79303468483966e-7 … -1.2077838445717075e-7 -4.793035011600566e-7; … ; 1.3626259536824603e-7 1.1356228511267841e-7 … 4.711986484617805e-8 1.135622712830793e-7; 1.3305484154092766e-7 1.3632648048035677e-7 … 5.2006845966664556e-8 1.3632643637322542e-7;;; -2.5020132846279223e-7 -1.740489227659538e-7 … -3.5564914621223936e-8 -1.740489608227779e-7; -4.793034593643464e-7 -3.6797304394550773e-7 … -1.1202888957292676e-7 -3.679730669313745e-7; … ; 1.1356228719693379e-7 1.1211097424068844e-7 … 1.0980679978504924e-7 1.1211096040749215e-7; 1.3632648619467527e-7 1.255087346850273e-7 … 5.094908266076652e-8 1.2550870128035138e-7;;; -4.8806376605669176e-8 -3.556489641774088e-8 … -1.141213254175177e-8 -3.55649043059613e-8; -1.207783682539512e-7 -1.1202887770941325e-7 … -9.786219296746327e-8 -1.1202888176494863e-7; … ; 4.7119880198475075e-8 1.0980682487090732e-7 … 2.4938228341867615e-7 1.0980681227584752e-7; 5.200687084818837e-8 5.0949104936225836e-8 … 4.910627491166115e-8 5.094909390588003e-8;;; … ;;; 6.475760588548557e-8 4.11033293866055e-8 … -6.681690006390529e-9 4.110334267208938e-8; 1.6897558806921943e-7 6.653287156490303e-8 … -1.329444704150151e-7 6.653288097744065e-8; … ; 3.267761009311453e-8 1.601736543014593e-7 … 3.8470357589678755e-7 1.6017363853120763e-7; -9.141993113979529e-8 -3.583035810142674e-8 … 7.037282122193081e-8 -3.5830356374590135e-8;;; -4.880640911178721e-8 -3.55649204113034e-8 … -1.1412139526172226e-8 -3.5564927943488905e-8; -1.207783904223335e-7 -1.1202889708971269e-7 … -9.786220571501056e-8 -1.1202889989720436e-7; … ; 4.7119865264429665e-8 1.0980679921732012e-7 … 2.4938223883586084e-7 1.0980678713231665e-7; 5.20068445595874e-8 5.0949080884429625e-8 … 4.9106254124008924e-8 5.0949069495468303e-8;;; -2.5020140074122773e-7 -1.7404897511230986e-7 … -3.5564925852274354e-8 -1.7404901324486756e-7; -4.793035099642744e-7 -3.679730802392024e-7 … -1.1202889692526717e-7 -3.679731032329048e-7; … ; 1.135622697097756e-7 1.1211095745236272e-7 … 1.0980678751950577e-7 1.1211094475771304e-7; 1.363264307162556e-7 1.255086931276646e-7 … 5.0949070004456516e-8 1.2550865974637848e-7;;;;], δHtotψ = Matrix{ComplexF64}[[0.00041458856058274106 - 0.0018718160380930494im 0.07690917942234693 + 0.056271611999682385im -0.5056034840494343 + 1.1317457749527968im 1.0824045393977377 - 0.9833869269668637im; -0.19321183769144842 - 0.04464050866665708im 1.0357012928315545 - 1.2545907779859584im -0.11210331519324214 + 0.25022829516457057im 0.23929404007266813 - 0.21639563396715458im; … ; -0.047787454208052846 - 0.01092302008755368im -0.012173070173234381 + 0.016725525957374553im 0.011844469643221512 - 0.019885184385621136im -0.011059120898456419 + 0.02041354790681395im; 0.0707203859215055 + 0.01538485471469039im 0.015157384882057324 - 0.020619264713648514im 0.055720214987927466 - 0.014957699732668604im 0.0373950733471086 + 0.0761733063447436im]], δVind = [0.001965431082119373 0.0020809994512471885 … 0.003153664051738269 0.0020809999204517374; 0.0038891300046648876 0.0047792682183312426 … 0.01027667074791956 0.004779268731389801; … ; -0.006150849385310499 -0.008269605753453901 … -0.017302973892658297 -0.008269604876506579; 0.0002386813870797744 -0.00039869864991833384 … -0.003664696011251066 -0.00039869810450187965;;; 0.0020809993731896724 0.0022232266435325077 … 0.0034201180904721647 0.0022232271714784685; 0.004779268122220703 0.005734951878895964 … 0.010903596705954774 0.005734952461016542; … ; -0.008269605824721448 -0.009601011280970961 … -0.012623114377607634 -0.009601010432166582; -0.0003986987257046397 -0.0010730148553308927 … -0.004230682625357031 -0.0010730142541845206;;; 0.0031536621216637453 0.003420116079222781 … 0.004832502255293254 0.0034201170184639937; 0.01027666852793751 0.010903594549583775 … 0.012272139848555458 0.010903595575683826; … ; -0.017302975933528537 -0.012623116223090473 … -0.009620267456011874 -0.01262311538562038; -0.0036646980755979967 -0.0042306845830892776 … -0.005491111343507101 -0.004230683694742655;;; … ;;; -0.0023365601559070083 -0.0028154624502761357 … 0.0056933205420412825 -0.002815463214374956; -0.012913092446630406 -0.01784117907826762 … 0.011498341369404275 -0.01784117969070215; … ; -0.0251744692837814 -0.012810887500531387 … -0.008395874712488639 -0.01281088670673496; 0.007597978739507818 0.01022541359201534 … -0.005449733027414053 0.010225412085990123;;; 0.0031536642937590657 0.0034201184931730755 … 0.004832506006368729 0.003420119449344792; 0.010276670993157638 0.010903597122160389 … 0.012272142743407625 0.010903598197117626; … ; -0.017302973461279555 -0.01262311392619906 … -0.009620265534537148 -0.012623113038359337; -0.0036646957030961195 -0.004230682156938105 … -0.005491108811495252 -0.004230681269321244;;; 0.0020809999961069266 0.002223227342305162 … 0.0034201193193709382 0.0022232278706468623; 0.004779268823044374 0.005734952669232945 … 0.010903598079814763 0.005734953249503538; … ; -0.008269604773797462 -0.009601010243715688 … -0.012623113187439822 -0.009601009336882554; -0.0003986980271687957 -0.0010730140838400737 … -0.004230681435817127 -0.001073013478313243;;;;], δρ0 = [-3.6021385946222735e-7 -2.4888415261141847e-7 … -4.8528335418718674e-8 -2.4888414825501955e-7; -5.972248655863765e-7 -4.6069770138249886e-7 … -1.1905763019863056e-7 -4.606976926415739e-7; … ; 1.2735789114871873e-7 1.0836338283178987e-7 … 4.6407153597643186e-8 1.0836337872679337e-7; 1.035201559063892e-7 1.1957967516776343e-7 … 5.0695065777730104e-8 1.1957967292483328e-7;;; -2.488841498346476e-7 -1.731066352107077e-7 … -3.53569513410847e-8 -1.7310663166997837e-7; -4.606976976175333e-7 -3.572169763199648e-7 … -1.1091434549576716e-7 -3.57216968597007e-7; … ; 1.0836338368181699e-7 1.084649595376558e-7 … 1.0862234548502361e-7 1.084649545018694e-7; 1.1957967574533243e-7 1.1611151499105328e-7 … 5.0208308306212457e-8 1.1611151253081055e-7;;; -4.852833687398599e-8 -3.53569534232317e-8 … -1.1341302516519589e-8 -3.535695178045271e-8; -1.190576354557469e-7 -1.1091435201201293e-7 … -9.77066435108207e-8 -1.109143474241833e-7; … ; 4.640716161418666e-8 1.0862235519358215e-7 … 2.484524422348936e-7 1.0862234983498356e-7; 5.069506909350403e-8 5.020831213012756e-8 … 4.926858477596564e-8 5.020830977795683e-8;;; … ;;; 6.434567228056484e-8 4.0838604271016224e-8 … -6.637921030740612e-9 4.083860533215142e-8; 1.6975034435173496e-7 6.683085462722084e-8 … -1.3340476936881363e-7 6.683085929905087e-8; … ; 3.2709966037771306e-8 1.602352310646875e-7 … 3.839014517377262e-7 1.6023522469859063e-7; -9.27136987894434e-8 -3.63362248136181e-8 … 7.126309625484648e-8 -3.6336227920624663e-8;;; -4.852833723301782e-8 -3.535695345150498e-8 … -1.1341301442462866e-8 -3.5356951419665066e-8; -1.1905763364742509e-7 -1.1091434961693473e-7 … -9.770663813995872e-8 -1.109143438034831e-7; … ; 4.6407154843472745e-8 1.0862234767851486e-7 … 2.4845243595688367e-7 1.0862234291462476e-7; 5.069506681889274e-8 5.0208309824449545e-8 … 4.926858210367698e-8 5.020830719441316e-8;;; -2.48884150758517e-7 -1.7310663546725746e-7 … -3.53569505842685e-8 -1.731066319030248e-7; -4.606976960694748e-7 -3.5721697394898257e-7 … -1.1091434197876108e-7 -3.572169662549465e-7; … ; 1.0836337840387641e-7 1.0846495377155683e-7 … 1.0862234232897692e-7 1.0846494990430504e-7; 1.1957967250955214e-7 1.1611151207250826e-7 … 5.020830656012459e-8 1.1611150979546706e-7;;;;], δeigenvalues = [[0.0004950055064465038, 0.09679853196627522, 0.03457865064865978, 0.02575503504466709]], δ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.15475965280045265 -0.15461620267745244 … -0.15418753097220486 -0.1546162026741451; -0.15461620267505652 -0.15447229887596287 … -0.1540422946484779 -0.15447229887264158; … ; -0.15418753097695603 -0.1540422946556699 … -0.1536083373216377 -0.15404229465232797; -0.15461620267653708 -0.15447229887744457 … -0.1540422946499839 -0.15447229887413022;;; -0.15461620267565598 -0.1544722988765626 … -0.15404229464908817 -0.15447229887324448; -0.1544722988741597 -0.15432794306708641 … -0.1538966028009944 -0.15432794306375439; … ; -0.15404229465385325 -0.15389660280820744 … -0.15346128698518047 -0.15389660280485462; -0.15447229887564354 -0.15432794306857153 … -0.1538966028025038 -0.15432794306524641;;; -0.1541875309752239 -0.15404229465393607 … -0.1536083373198763 -0.15404229465058653; -0.15404229465151106 -0.1538966028058633 … -0.15346128698279796 -0.15389660280249964; … ; -0.1536083373246846 -0.15346128699007708 … -0.1530219131096923 -0.1534612869866924; -0.1540422946530075 -0.15389660280736103 … -0.1534612869843205 -0.15389660280400438;;; … ;;; -0.15347880583180262 -0.1533313919433587 … -0.15289092490128836 -0.1533313919399595; -0.15333139194089437 -0.15318350572049086 … -0.1527416287522263 -0.15318350571707695; … ; -0.15289092490617426 -0.15274162875962408 … -0.1522955395659082 -0.15274162875618905; -0.15333139194241965 -0.15318350572201725 … -0.15274162875377914 -0.15318350571861103;;; -0.15418753097581897 -0.15404229465453242 … -0.1536083373204801 -0.1540422946511841; -0.15404229465210562 -0.1538966028064592 … -0.15346128698340117 -0.1538966028030966; … ; -0.15360833732529203 -0.1534612869906858 … -0.1530219131103093 -0.15346128698730263; -0.15404229465360691 -0.15389660280796166 … -0.15346128698492886 -0.1538966028046064;;; -0.15461620267595597 -0.15447229887686317 … -0.15404229464939218 -0.15447229887354566; -0.15447229887445924 -0.15432794306738665 … -0.1538966028012981 -0.1543279430640551; … ; -0.1540422946541591 -0.15389660280851394 … -0.153461286985491 -0.15389660280516185; -0.15447229887594563 -0.15432794306887418 … -0.1538966028028101 -0.1543279430655497]), 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.015163897554553812 -0.012550129031123902 … -0.004966867884390715 -0.012550128964479611; -0.012550128994675409 -0.01036259495906387 … -0.004135158637963965 -0.01036259489123395; … ; -0.0049668679643341315 -0.004135158756435617 … -0.002069779755349279 -0.004135158675712673; -0.012550128996084672 -0.010362594958202801 … -0.004135158627741086 -0.010362594891487501;;; -0.012550129016777796 -0.010362594986834088 … -0.004135158664993191 -0.010362594913288708; -0.010362594942643834 -0.008622038897830465 … -0.0038975110602573718 -0.00862203882443778; … ; -0.004135158779170937 -0.0038975112104864683 … -0.0034841066293482244 -0.00389751111161924; -0.010362594953169189 -0.008622038906829256 … -0.0038975110439426522 -0.008622038828634182;;; -0.004966867988663603 -0.004135158801669738 … -0.002069779801685584 -0.004135158704869956; -0.004135158744069953 -0.00389751118468738 … -0.0034841066503056603 -0.0038975110998243486; … ; -0.0020697799743110267 -0.0034841068055072905 … -0.005620579742620944 -0.0034841067076730073; -0.0041351587600181455 -0.0038975112038861226 … -0.0034841066260165253 -0.003897511098557551;;; … ;;; -0.006547898054144287 -0.0048763401744908695 … -0.0014866436173651353 -0.0048763402846209725; -0.004876340240894248 -0.002677988942940639 … -0.0041230539530618645 -0.0026779890736115513; … ; -0.00148664381547438 -0.004123054115760949 … -0.00717089503203369 -0.004123054024715012; -0.004876340216568444 -0.002677988890731685 … -0.004123053976738384 -0.0026779890445703877;;; -0.004966867950694195 -0.004135158753100006 … -0.0020697796457211583 -0.004135158627730152; -0.004135158692443075 -0.0038975111206085897 … -0.0034841065194985707 -0.003897511008448328; … ; -0.0020697797850250804 -0.0034841066685249176 … -0.005620579680482571 -0.0034841065835273894; -0.004135158692985443 -0.00389751113182956 … -0.0034841065357567036 -0.0038975110135427904;;; -0.012550128980807625 -0.01036259494222554 … -0.004135158587321438 -0.010362594869882651; -0.010362594907939598 -0.00862203885207757 … -0.003897510977666579 -0.00862203877977381; … ; -0.004135158667062228 -0.003897511091573859 … -0.0034841065735365536 -0.0038975110195613462; -0.010362594902578435 -0.008622038847446983 … -0.003897510977958892 -0.008622038776970434])], 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.00953289228615185 -0.006919175134768446 … 0.0006640015060827547 -0.006919175064816801; -0.006919175095923897 -0.004731678427338914 … 0.001495697379537913 -0.004731678356187707; … ; 0.0006640014213882498 0.001495697253874181 … 0.0035610468755768432 0.0014956973379390628; -0.00691917509881367 -0.0047316784279595216 … 0.00149569738825483 -0.004731678357929873;;; -0.006919175118625916 -0.004731678455709041 … 0.0014956973518982872 -0.00473167837884551; -0.00473167840911571 -0.0029911474622591923 … 0.0017333332822805666 -0.002991147385534479; … ; 0.001495697232955547 0.0017333331248383789 … 0.002146715835997684 0.0017333332270584251; -0.004731678421124878 -0.0029911474727430734 … 0.0017333332970859657 -0.0029911473912228818;;; 0.0006640013987908094 0.0014956972103738397 … 0.003561046831001935 0.0014956973105231652; 0.0014956972703987128 0.001733333152981703 … 0.0021467158174227866 0.0017333332412083772; … ; 0.0035610466535682273 0.002146715654941979 … 1.0229566642130407e-5 0.0021467157561609438; 0.0014956972529541618 0.0017333331322852141 … 0.002146715840189445 0.0017333332409704338;;; … ;;; -0.000917069960377233 0.0007544784619587831 … 0.00414415877399269 0.0007544783552278778; 0.0007544783980197945 0.0029528245401417253 … 0.0015077456838542584 0.002952824412884721; … ; 0.004144158570997576 0.0015077455137573418 … -0.0015401020754098865 0.0015077456082383082; 0.0007544784208203467 0.0029528245908243445 … 0.0015077456586249805 0.0029528244403918616;;; 0.0006640014361651379 0.001495697258347188 … 0.003561046986362566 0.001495697387065364; 0.0014956973214310387 0.0017333332164645812 … 0.002146715947626664 0.0017333333319874307; … ; 0.003561046842246715 0.0021467157913156168 … 1.0229628163468778e-5 0.0021467158796963275; 0.001495697319387427 0.001733333203741146 … 0.002146715929840892 0.0017333333253831762;;; -0.006919175082955727 -0.004731678411401059 … 0.0014956974292660331 -0.004731678335740629; -0.004731678374711012 -0.002991147416806502 … 0.0017333333645676856 -0.0029911473411711846; … ; 0.0014956973447583614 0.0017333332434444834 … 0.0021467158914987975 0.0017333333188090925; -0.0047316783708362155 -0.002991147413663447 … 0.001733333362763443 -0.002991147339862419]), 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 = [8.212796226349338e-7 - 2.2244433728735773e-6im -0.0047949016615373265 + 0.005063854914659877im … -0.003234001082633893 + 0.003354887268613884im 0.004746894252256578 - 0.004955960255190342im; -0.017650251513856862 - 0.003961870188668333im 0.019039946193136814 + 0.0006508603231633221im … -0.007763626662045071 - 0.00451381639445706im 0.01344248431056469 + 0.006523038624127486im; … ; -0.021536503280155116 - 0.004807042035852245im 0.017961982919208155 + 0.005240431895649267im … -0.013329138790146274 - 0.0016292288706140036im 0.019886384534403598 + 0.0032257812244439786im; 0.017682134087742395 + 0.0038635471061370173im -0.013434238764782372 - 0.00653722197707043im … 0.012212510630336876 - 9.155065931891812e-5im -0.01903125851180896 - 0.0006672661413657608im;;; -0.004175496290915751 - 0.0015379444358222499im 0.005259736953576711 - 0.0020670150754855377im … 0.0013808003220546483 - 0.002610455600374261im -0.0003361401403957129 + 0.0038048158706924192im; 0.018676629633090518 + 0.004518348693022379im -0.0173683585466164 - 0.0017748761678451385im … 0.007436646808090988 + 0.003825596864159024im -0.013762783868606797 - 0.005555677745516735im; … ; 0.01808474670100823 + 0.003913553457080991im -0.01491975254659591 - 0.004200522108127385im … 0.011212303893313668 + 0.0013202024763668643im -0.016430016769149174 - 0.002618812807267672im; -0.013795521774697249 - 0.002669858620015765im 0.010600939841067428 + 0.004415552800962709im … -0.009435778261945684 + 9.055836743219777e-5im 0.01420637964802849 + 0.0006358323590074466im;;; 0.002853411922818387 + 0.001024961387776171im -0.0029261696000975105 + 0.0003177378212158242im … -0.0001256140062758396 + 0.0014392670556674942im -0.0010029626412354042 - 0.0016982726851369504im; -0.011947782199394087 - 0.002912562213538406im 0.010657907141285401 + 0.0015882032079863394im … -0.005299332544787646 - 0.0023987421264790114im 0.009148162758399896 + 0.003170692999929459im; … ; -0.011314847940276429 - 0.0023984943599640735im 0.009469086291432127 + 0.0025639787348713584im … -0.007502735528867157 - 0.0008716070111523457im 0.01033707083420891 + 0.00165455182298296im; 0.008027262096845142 + 0.0015083451477100125im -0.006214057059042124 - 0.0021461042155017764im … 0.0055087692302681025 - 6.188099329234469e-6im -0.00772389085394536 - 0.0005638660468243429im;;; … ;;; 0.0012507574480228672 + 0.0005192043873377922im -0.0006769124500654408 - 0.0008399352671463588im … 0.001222448570805775 - 0.00024485879539014245im -0.0014328925280147598 - 4.7396990013893956e-5im; 0.0036744523497838303 + 0.0006045818718617934im -0.0036089535292077066 - 0.00019765571535257763im … 0.0019311586928927949 + 0.0008866928943243737im -0.002952933967434764 - 0.0008854333231054854im; … ; 0.007368357230011615 + 0.0017799705563966953im -0.006412858608887526 - 0.0018487029147651365im … 0.0053805977869227795 + 0.00085277325089116im -0.0068745975733764016 - 0.001364774188744023im; -0.005785282632512409 - 0.0014898350136125778im 0.00476778133520421 + 0.0016569166368294239im … -0.004115235311663977 - 0.0004519326378529261im 0.005423888228155888 + 0.0009691920909182286im;;; -0.002840405541659186 - 0.0011340187013420971im 0.000985138626558922 + 0.0018008750304456835im … -0.00200655855141373 + 0.0007263400522711763im 0.002908924883422826 - 0.0002154462040541492im; -0.008037380455898042 - 0.0014303105156183795im 0.0077382817001564 + 0.00048141382830355756im … -0.0037231951698999044 - 0.001801470963745037im 0.006228363210637071 + 0.002064098376564632im; … ; -0.01319600589791739 - 0.0031123397020255293im 0.011045720123322817 + 0.0031622633878569647im … -0.008531858833033137 - 0.0012619786197628083im 0.011913764863335839 + 0.002252789887944609im; 0.011939563388385355 + 0.002992713612554341im -0.009134789765833514 - 0.003254275992427117im … 0.007085421279241373 + 0.0005916464569708315im -0.010644849578859602 - 0.0016718191950958455im;;; 0.0041708733035173105 + 0.0016258028548459373im 0.0003662419559011745 - 0.00392173188485791im … 0.0030625208025882407 - 0.001973417611695526im -0.005230178407473484 + 0.0019505630134744314im; 0.013785294656383822 + 0.0026650833272225363im -0.014217641088154001 - 0.0005810021334396006im … 0.00611687601665589 + 0.00332567972022166im -0.010611830988813367 - 0.004361995083128369im; … ; 0.01976670426770611 + 0.004551721849251716im -0.016239970136111132 - 0.004701518286782176im … 0.0119713728711602 + 0.001608180469658738im -0.017750334722813323 - 0.0031198049295768264im; -0.018690045942507998 - 0.004526360274286354im 0.013753350797358914 + 0.005611283866954491im … -0.01075605217377722 - 0.0004099859368880596im 0.017359048348811017 + 0.001831347076053398im],)])]), [8.519873142446504e-7 4.6272647214074165e-7 … 2.417777380283628e-8 4.6272646423091214e-7; 4.6272646781478133e-7 2.5014473840840273e-7 … 1.358299080312497e-8 2.501447331621936e-7; … ; 2.4177775029475603e-8 1.3582992025351121e-8 … 1.5659203719618175e-9 1.3582991192564156e-8; 4.6272646798204085e-7 2.5014473834181133e-7 … 1.3582990697658859e-8 2.5014473318181303e-7;;; 4.627264704380604e-7 2.501447405562501e-7 … 1.3582991081972956e-8 2.501447348679952e-7; 2.5014473713841267e-7 1.3895606631187236e-7 … 1.1279984294636909e-8 1.389560625384559e-7; … ; 1.3582992259898271e-8 1.127998565865606e-8 … 7.937750685762435e-9 1.1279984760984006e-8; 2.5014473795248317e-7 1.3895606677453362e-7 … 1.1279984146505756e-8 1.3895606275420322e-7;;; 2.4177775402787835e-8 1.3582992492008313e-8 … 1.5659204807714555e-9 1.3582991493364012e-8; 1.3582991897777583e-8 1.1279985424412042e-8 … 7.937750835244925e-9 1.1279984653891129e-8; … ; 1.5659208861351505e-9 7.937751942247075e-9 … 3.571855245555526e-8 7.937751244429174e-9; 1.3582992062307001e-8 1.1279985598730499e-8 … 7.937750661997673e-9 1.1279984642388484e-8;;; … ;;; 5.791842799354914e-8 2.2816029144395212e-8 … 5.620363587766437e-10 2.2816030768042046e-8; 2.2816030123379776e-8 3.489279896591041e-9 … 1.3458506008293471e-8 3.4892804271947847e-9; … ; 5.620365901553871e-10 1.3458507676199028e-8 … 7.728314092280541e-8 1.3458506742841116e-8; 2.2816029764739832e-8 3.489279684590265e-9 … 1.3458506251013745e-8 3.4892803092707358e-9;;; 2.4177774820185167e-8 1.3582991990937067e-8 … 1.5659201145319964e-9 1.358299069754425e-8; 1.3582991365164487e-8 1.1279984842604585e-8 … 7.937749902242109e-9 1.12799838242316e-8; … ; 1.5659204416495228e-9 7.937750965199239e-9 … 3.5718551207664735e-8 7.937750358939113e-9; 1.3582991370756348e-8 1.1279984944483365e-8 … 7.937750018202777e-9 1.1279983870486514e-8;;; 4.627264661688552e-7 2.5014473710605934e-7 … 1.3582990280666452e-8 2.501447315108065e-7; 2.501447344542772e-7 1.3895606395952939e-7 … 1.1279983544745098e-8 1.3895606024209737e-7; … ; 1.3582991103319994e-8 1.1279984578980971e-8 … 7.937750287675493e-9 1.1279983925133386e-8; 2.5014473403962755e-7 1.3895606372144928e-7 … 1.1279983547398445e-8 1.3895606009796065e-7;;;;], Matrix{ComplexF64}[[0.05633569820119774 - 0.25434869724796977im 0.778865303651872 + 0.5698670480755981im 2.865954991979768e-8 - 1.925833116284221e-6im 2.4770117691397414e-6 - 1.1456723108640971e-7im; -0.03822258281152441 + 0.17257022545800954im 0.03431987819532846 + 0.025110603093877453im 0.35796462386690947 + 0.15925728963326746im -0.31241318766035586 - 0.3419224040871587im; … ; 0.012143600107410562 - 0.054826849903903324im -0.012413047643244251 - 0.009082169906620384im 0.005474208774649323 - 5.965468885519997e-5im -0.001476975884091325 - 0.007614855801547005im; -0.021773558227437765 + 0.09830491751419553im 0.0245912132055402 + 0.017992483036387532im -0.011200978052833493 + 0.00888042239841821im -0.008365313187383731 + 0.023854832855714606im]], [[2.0, 0.0, 0.0, 0.0]], -0.2823944087959322, [[-0.5541104196781454, -0.010678397913718978, 0.14495665121455126, 0.1449571543504763]], 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.620692457147452e-7, -6.293614724885022e-7, -2.2951201694369196e-7, 4.994026040359059e-7, 7.923312341219848e-7, 5.851072830318005e-7, 6.516971822157196e-7, 1.018930964211686e-6, 3.6385927624140693e-7, -2.2010932896884465e-6  …  2.781018273666848e-6, 1.2016954812407125e-6, 2.4456291491075683e-8, 1.0078538775476933e-7, 5.535307672364371e-7, 4.534172567004623e-7, 4.631882759242121e-9, -1.0437056657621727e-7, 1.1211094475771304e-7, 1.2550865974637848e-7], resid_history = [0.24939209612074564, 0.0037664958040930004, 0.0002851899306825758, 4.681853336031085e-6, 4.159726024216522e-8, 8.816776557979112e-11, 4.006773727999861e-8, 1.0153928567942804e-9], converged = true, n_iter = 8, residual_norm = 1.0153928567942804e-9, maxiter = 100, tol = 1.0e-8, s = 0.9354784264705803, restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.26113046468934e-7, -3.989131060644207e-8, 1.0153928567942804e-9, -0.00020439947444484167, 4.686880828708941e-6, -4.213469377261484e-8, 8.816776557979112e-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}}([[5.185571693150068e-7, 2.2355512485244417e-6, 7.70948588539463e-7, -1.338156906081932e-6, -1.696688117001198e-6, -1.1093102306348478e-6, -1.17593515344875e-6, -1.6887026791815817e-6, -5.096225445087793e-7, 2.609023206560252e-6  …  -3.827636630226312e-6, -1.78768123258602e-6, -4.1682806084139516e-8, -1.8614321264056285e-7, -1.0199888574303905e-6, -8.291756482138824e-7, -9.238653828487978e-9, 2.518941169027369e-7, -3.441092734726008e-7, -5.74436991395967e-7], [-4.283309022828966e-7, -2.8943937752740894e-6, -1.0206908659254195e-6, 1.7855217550276289e-6, 2.2827051809463644e-6, 1.5042420637813256e-6, 1.600495011365326e-6, 2.3023468422340524e-6, 6.985404708136821e-7, -3.615831611855128e-6  …  5.792608294902931e-6, 2.7149981753590464e-6, 6.32410235692311e-8, 2.8224002120239973e-7, 1.5532283746110053e-6, 1.2736057800449825e-6, 1.4278448026583806e-8, -3.895517263703689e-7, 5.346690288445391e-7, 9.420361210957087e-7]]), H = [1.172838815200667 0.004443401494593711 … 0.0 0.0; 0.1236013757651565 1.0098046494237365 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.179333789277804 0.11025261791814642 … 0.0 0.0; 0.0 1.0041001062449362 … 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.9257125454269264
Interacting polarizability:     1.7736548739025253

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 = 6.95e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.03e-11
[ Info: GMRES linsolve in iteration 2; step 2: normres = 4.28e-12
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 2.79e-13
└ * number of operations = 12
Non-interacting polarizability: 1.92571254542684
Interacting polarizability:     1.7736548709006719

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