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

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

polarizability = (με - μref) / ε

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

Reference dipole:  -0.00013457275193736801
Displaced dipole:  0.01761222157846419
Polarizability :   1.774679433040156

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.0028563942650291505 + 0.0016443638030818706im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.15660795016698123 - 0.2668335533473186im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; 0.02063123121103346 - 0.035506325444359714im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.04011211554539665 + 0.07029943181407411im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.6206944835273916e-7 -2.502015178089527e-7 … -4.880642270061191e-8 -2.50201509849736e-7; -6.293616501653438e-7 -4.793036099705986e-7 … -1.207783956936092e-7 -4.793036063767416e-7; … ; 1.3626256726898847e-7 1.1356223333058855e-7 … 4.7119845123356586e-8 1.1356225555574146e-7; 1.3305471337377948e-7 1.363263495027412e-7 … 5.200682583374215e-8 1.3632636193726367e-7;;; -2.5020152567221374e-7 -1.7404908869023667e-7 … -3.5564931388564535e-8 -1.7404907957140212e-7; -4.793036188822208e-7 -3.679731703144727e-7 … -1.1202889082118578e-7 -3.679731582116844e-7; … ; 1.135622454823148e-7 1.1211090371593626e-7 … 1.0980676151495e-7 1.1211092635302942e-7; 1.3632634538224837e-7 1.2550859866398806e-7 … 5.0949049831918744e-8 1.2550860964789683e-7;;; -4.880642860023271e-8 -3.556493866131312e-8 … -1.1412128845416636e-8 -3.556493105671611e-8; -1.2077840523404252e-7 -1.1202891402799512e-7 … -9.786217397922557e-8 -1.1202889108888014e-7; … ; 4.711984486273318e-8 1.0980674984811817e-7 … 2.4938218700323907e-7 1.098067536161508e-7; 5.200681628693593e-8 5.094904295322132e-8 … 4.91062225099018e-8 5.094904128460284e-8;;; … ;;; 6.475761218724441e-8 4.11033390980644e-8 … -6.681687544341521e-9 4.1103340831190754e-8; 1.6897559465857775e-7 6.653287782297146e-8 … -1.3294446459466066e-7 6.653288064221034e-8; … ; 3.26775944469049e-8 1.601736359069805e-7 … 3.8470357238617813e-7 1.6017362222271214e-7; -9.141993098083017e-8 -3.5830354509771e-8 … 7.03728189275165e-8 -3.583036445383233e-8;;; -4.880641962984084e-8 -3.5564932450054875e-8 … -1.1412134532636058e-8 -3.556492947750182e-8; -1.2077838899826062e-7 -1.1202889500526232e-7 … -9.786219334734835e-8 -1.1202889417274958e-7; … ; 4.711984930979262e-8 1.0980676868878467e-7 … 2.493822089394341e-7 1.0980676773567187e-7; 5.200683069889287e-8 5.094906021121758e-8 … 4.9106239533798565e-8 5.094905812140229e-8;;; -2.502015009361876e-7 -1.7404907044621123e-7 … -3.556493015660691e-8 -1.740490649569175e-7; -4.793035966554729e-7 -3.6797315045925393e-7 … -1.1202889732930468e-7 -3.679731501957504e-7; … ; 1.1356225448306612e-7 1.1211092032648347e-7 … 1.0980675930587745e-7 1.1211093275684252e-7; 1.3632636870016454e-7 1.2550861988719365e-7 … 5.094905426761161e-8 1.2550862677626904e-7;;;;], δHtotψ = Matrix{ComplexF64}[[-0.001661528337849327 - 0.0009565055612739707im 0.08916863064043115 + 0.03362259215045882im 1.774706296860182 - 0.4523589775294354im 0.7794062408538293 - 0.6348815880935819im; -0.10049319106943794 + 0.1709523543624295im 0.6619355554443459 - 1.4861078297150658im 0.3916556357179464 - 0.10067231966318037im 0.17336672353454136 - 0.13977339950819234im; … ; -0.024742595165231338 + 0.042317338531383523im -0.00724990486121883 + 0.019374335296735073im -0.025656631896508478 + 0.009941470424409272im -0.01793868889322853 + 0.001432696991154198im; 0.03587222166659097 - 0.06285897284001413im 0.009082622073403067 - 0.023925040908638864im 0.002450023895697015 + 0.035842519175763025im -0.03914328831454175 - 0.11112797466187972im]], δVind = [0.0019654323107566973 0.002081001145389207 … 0.003153665694978549 0.0020810010443293395; 0.0038891318171413564 0.004779270589994403 … 0.010276673549128004 0.004779270480965198; … ; -0.0061508483715722335 -0.008269604786327699 … -0.017302975812361867 -0.008269604207530853; 0.0002386823672659293 -0.0003986972085233793 … -0.003664695193063266 -0.0003986972896345901;;; 0.002081001201814322 0.002223229057845647 … 0.0034201208961170686 0.0022232289379720627; 0.004779270635404839 0.0057349550674706835 … 0.010903600997039622 0.005734954994211987; … ; -0.008269604104582587 -0.009601009871877914 … -0.012623114144323547 -0.009601009191128369; -0.0003986971184974219 -0.0010730127103708133 … -0.004230681155764242 -0.0010730128089792583;;; 0.003153666251793423 0.003420121318062836 … 0.0048325095560995705 0.003420121311562108; 0.010276673932530705 0.010903600711646666 … 0.012272146557370857 0.010903601383754563; … ; -0.017302974527795733 -0.012623114007994603 … -0.009620266447695447 -0.012623114043008966; -0.0036646946950073325 -0.004230680723338377 … -0.005491109671161525 -0.004230681060463461;;; … ;;; -0.002336560576669855 -0.0028154634255335865 … 0.005693320134353359 -0.002815463259159096; -0.012913093554751643 -0.017841181110095347 … 0.011498341655143528 -0.01784118069308052; … ; -0.025174475864885193 -0.012810888317976376 … -0.008395876436301023 -0.012810888975897302; 0.007597978807788234 0.010225413062008266 … -0.0054497345320220015 0.010225412334046368;;; 0.0031536654206559185 0.0034201203335094458 … 0.0048325071694286835 0.003420120268786234; 0.010276673435515337 0.010903600351190796 … 0.012272144549100123 0.010903600203443645; … ; -0.017302975924629573 -0.01262311464535541 … -0.0096202672172203 -0.012623114878361784; -0.0036646954648854885 -0.004230681397960147 … -0.005491110092094317 -0.00423068165252879;;; 0.002081000975038589 0.0022232287850325098 … 0.003420120458912357 0.002223228694139178; 0.004779270429477005 0.005734954853486725 … 0.010903600291426476 0.005734954729557276; … ; -0.00826960451900002 -0.009601009984996302 … -0.01262311504042522 -0.009601009674451896; -0.00039869738610623265 -0.0010730130056154413 … -0.004230681573953782 -0.0010730130836023866;;;;], δρ0 = [-3.602140311500052e-7 -2.4888427504292746e-7 … -4.8528364114215074e-8 -2.488842780478926e-7; -5.972249656264886e-7 -4.606977629209352e-7 … -1.1905764085589991e-7 -4.6069776719825414e-7; … ; 1.2735784478475158e-7 1.0836332474523977e-7 … 4.640712626735626e-8 1.0836334416160513e-7; 1.0351999578255325e-7 1.195795481785842e-7 … 5.069503027345849e-8 1.1957955198573838e-7;;; -2.488842786908837e-7 -1.7310672606537437e-7 … -3.535696793244721e-8 -1.7310672509656242e-7; -4.60697768613364e-7 -3.5721701967377723e-7 … -1.1091433960729436e-7 -3.5721701372182537e-7; … ; 1.0836333695753406e-7 1.0846489082538628e-7 … 1.0862229587839199e-7 1.0846491106581255e-7; 1.1957954687646981e-7 1.1611140796993357e-7 … 5.0208268725278805e-8 1.1611141254648533e-7;;; -4.8528363200251615e-8 -3.535697074144881e-8 … -1.1341296255799234e-8 -3.535696508180108e-8; -1.1905764526035728e-7 -1.1091435861710139e-7 … -9.770661310118579e-8 -1.1091433760399327e-7; … ; 4.640712770844444e-8 1.0862228768886127e-7 … 2.484523349604893e-7 1.0862228997845634e-7; 5.0695025936349186e-8 5.0208265948396385e-8 … 4.92685327402969e-8 5.020826258381429e-8;;; … ;;; 6.434571743843957e-8 4.083863087277368e-8 … -6.6379246253809795e-9 4.083863508781042e-8; 1.6975036591305612e-7 6.683086391337644e-8 … -1.3340478652818128e-7 6.683086772980085e-8; … ; 3.270994361864006e-8 1.602351855509229e-7 … 3.8390134434161436e-7 1.6023516993093988e-7; -9.271366507406184e-8 -3.6336211109104084e-8 … 7.12630576532083e-8 -3.6336220207382166e-8;;; -4.8528363103785904e-8 -3.5356971150767223e-8 … -1.1341304114499152e-8 -3.535697009402216e-8; -1.1905763577573542e-7 -1.1091434571621141e-7 … -9.77066374439527e-8 -1.1091434662768787e-7; … ; 4.640712991169665e-8 1.0862230140520539e-7 … 2.484523462855609e-7 1.0862229901254709e-7; 5.069503359668511e-8 5.02082770611123e-8 … 4.926854480467535e-8 5.020827328216883e-8;;; -2.488842737098151e-7 -1.7310672227601512e-7 … -3.535696996689515e-8 -1.7310672485588473e-7; -4.6069776098774996e-7 -3.5721701065855787e-7 … -1.1091434901742841e-7 -3.5721701629599825e-7; … ; 1.0836334249630778e-7 1.0846490406394135e-7 … 1.0862229124209246e-7 1.0846491442001112e-7; 1.1957955554057216e-7 1.16111418143862e-7 … 5.0208270174667176e-8 1.1611141886099457e-7;;;;], δeigenvalues = [[0.0004950042707630642, 0.09679859211892729, 0.01604350814921992, 0.03955344254635487]], δ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.15475965280801438 -0.15461620268277657 … -0.1541875309841963 -0.15461620268393644; -0.15461620268399698 -0.1544722988826657 … -0.1540422946618817 -0.15447229888382372; … ; -0.1541875309817478 -0.15404229465818708 … -0.15360833733089324 -0.15404229465937463; -0.15461620268270837 -0.1544722988813668 … -0.15404229466058644 -0.1544722988825362;;; -0.1546162026831436 -0.1544722988818038 … -0.1540422946610266 -0.1544722988829727; -0.15447229888302652 -0.15432794307370648 … -0.1538966028143453 -0.15432794307487355; … ; -0.1540422946585721 -0.15389660281064246 … -0.15346128699438247 -0.15389660281183834; -0.1544722988817416 -0.15432794307241143 … -0.15389660281305328 -0.15432794307358946;;; -0.15418753098263205 -0.15404229465907587 … -0.15360833733178775 -0.154042294660261; -0.15404229466030916 -0.15389660281239229 … -0.15346128699613146 -0.15389660281357523; … ; -0.15360833732931167 -0.15346128699239658 … -0.15302191311885285 -0.15346128699360953; -0.15404229465901914 -0.15389660281109174 … -0.1534612869948355 -0.15389660281228684;;; … ;;; -0.15347880583954795 -0.15333139194881346 … -0.15289092491358752 -0.1533313919499952; -0.15333139195007436 -0.15318350572737996 … -0.1527416287659924 -0.15318350572855835; … ; -0.15289092491106576 -0.15274162876218328 … -0.15229553957538405 -0.15274162876339561; -0.15333139194872905 -0.1531835057260224 … -0.1527416287646427 -0.15318350572721579;;; -0.15418753098350388 -0.15404229465995886 … -0.15360833733265217 -0.15404229466112643; -0.15404229466119942 -0.1538966028132942 … -0.15346128699701345 -0.1538966028144586; … ; -0.153608337330169 -0.1534612869932644 … -0.15302191311970345 -0.15346128699446077; -0.15404229465988 -0.15389660281196316 … -0.1534612869956893 -0.15389660281314155;;; -0.1546162026835818 -0.1544722988822477 … -0.15404229466146166 -0.15447229888340763; -0.15447229888347463 -0.1543279430741606 … -0.1538966028147892 -0.15432794307531794; … ; -0.15404229465900326 -0.15389660281107875 … -0.15346128699481099 -0.15389660281226672; -0.15447229888217412 -0.1543279430728491 … -0.1538966028134833 -0.1543279430740193]), 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.015163897140633736 -0.01255012852423038 … -0.004966867572787583 -0.012550128612661603; -0.012550128659220157 -0.010362594550526735 … -0.004135158357479136 -0.010362594596950169; … ; -0.004966867759839658 -0.004135158055610828 … -0.002069779254219755 -0.004135158531072377; -0.012550128580077105 -0.010362594334631568 … -0.004135158203860031 -0.010362594546507845;;; -0.01255012854898762 -0.010362594412831837 … -0.00413515800393173 -0.0103625944155316; -0.010362594602535816 -0.008622038499086056 … -0.003897510451270147 -0.008622038410458095; … ; -0.0041351583791280865 -0.00389751028291731 … -0.0034841060913954943 -0.0038975107476841215; -0.010362594432877544 -0.00862203813290725 … -0.0038975103580856712 -0.00862203831748072;;; -0.004966867449127986 -0.004135158161534622 … -0.0020697780158851135 -0.0041351577554054825; -0.004135158446064222 -0.0038975108676772216 … -0.0034841055425784482 -0.0038975103989508076; … ; -0.0020697793113400587 -0.0034841059334129376 … -0.005620578960616341 -0.0034841059854291137; -0.004135158041476662 -0.003897510253095298 … -0.00348410529369065 -0.003897510085492599;;; … ;;; -0.006547897812521703 -0.004876339926465684 … -0.0014866434059704636 -0.004876340155779107; -0.004876340111757948 -0.0026779889054695826 … -0.00412305382699679 -0.002677988990679738; … ; -0.0014866433112189596 -0.004123053799468103 … -0.007170894509268705 -0.004123053556056119; -0.004876340181965358 -0.0026779888380340545 … -0.00412305350092018 -0.002677989289948622;;; -0.004966867626688624 -0.004135158374240435 … -0.0020697789824914332 -0.004135158249782578; -0.004135158270065262 -0.0038975106136083953 … -0.00348410610084671 -0.0038975106217137264; … ; -0.002069779335507503 -0.0034841061444292037 … -0.005620579027494084 -0.0034841061058811097; -0.004135158279441742 -0.0038975106479328648 … -0.0034841057565326962 -0.003897510491984141;;; -0.012550128598391067 -0.010362594489697972 … -0.004135158180550383 -0.010362594529463174; -0.010362594541602169 -0.008622038405507983 … -0.0038975106597012314 -0.008622038469778841; … ; -0.004135158458747522 -0.0038975105341561536 … -0.003484105964373853 -0.003897510780144567; -0.010362594521587071 -0.00862203832006621 … -0.003897510380668729 -0.008622038410464047])], 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.009532891879793504 -0.006919174633199054 … 0.0006640018056944513 -0.0069191747227901265; -0.006919174769409104 -0.004731678025504612 … 0.001495697646618937 -0.0047316780730860646; … ; 0.0006640016210909452 0.001495697952181789 … 0.0035610473674508266 0.0014956974755326903; -0.006919174688977389 -0.004731677808310512 … 0.0014956978015333387 -0.004731678021356187;;; -0.006919174658323361 -0.004731677886947987 … 0.0014956980010213256 -0.00473167789081662; -0.0047316780778745154 -0.002991147070134849 … 0.0017333338779168877 -0.0029911469826739544; … ; 0.0014956976282795333 0.0017333340499725127 … 0.0021467163647484137 0.0017333335840098245; -0.004731677906931299 -0.0029911467026609677 … 0.0017333339723934688 -0.0029911468884124687;;; 0.0006640019309182687 0.0014956978453691504 … 0.0035610486048909614 0.001495698250313155; 0.0014956975596063429 0.001733333463462862 … 0.0021467169118164976 0.0017333339310063333; … ; 0.003561047311912119 0.0021467165247168263 … 1.0230339486171774e-5 0.002146716471487704; 0.0014956979654840094 0.0017333340793453285 … 0.002146717162000315 0.001733334245752928;;; … ;;; -0.0009170697264999812 0.0007544787045292208 … 0.0041441589730882005 0.0007544784740340484; 0.0007544785179761041 0.002952824570723681 … 0.0015077457961532332 0.002952824484335135; … ; 0.004144159070361492 0.001507745827490985 … -0.001540101562120766 0.0015077460696906325; 0.0007544784491140344 0.0029528246395168177 … 0.0015077461235796247 0.0029528241864088716;;; 0.0006640017524858002 0.0014956976317803215 … 0.003561047637420222 0.0014956977550706126; 0.001495697734715043 0.0017333337166297708 … 0.002146716352666247 0.0017333337073600378; … ; 0.0035610472868873324 0.0021467163128327544 … 1.0230271757803722e-5 0.002146716350184472; 0.0014956977266580348 0.0017333336836363477 … 0.0021467166983044515 0.001733333838406681;;; -0.006919174708165013 -0.004731677964258017 … 0.001495697823967604 -0.004731678005183124; -0.004731678017388982 -0.0029911469770108577 … 0.001733333669041936 -0.0029911470424390674; … ; 0.0014956975482289152 0.0017333337982973793 … 0.002146716491341509 0.0017333335511209996; -0.004731677996073342 -0.0029911468902576056 … 0.0017333339493803938 -0.002991146981825645]), 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.440399432011966e-6 - 1.0527229499590173e-6im 0.0008058753634624159 - 0.00011047464424904298im … 0.00048224730889024633 - 0.00010006740221627286im -0.0007420719810677877 + 0.00012348151955380833im; -0.007548720528511245 + 0.013040222070590586im 0.006955291458536355 - 0.012803130925674061im … -0.005126923518028306 + 0.0084108359521579im 0.007859866464213251 - 0.012942455216106392im; … ; -0.010019675835671264 + 0.017305886032972614im 0.0091623690012577 - 0.015573709826993716im … -0.005791037427409589 + 0.010241823991377692im 0.008853502436683602 - 0.015525054471027677im; 0.007487780543715317 - 0.013079281046440593im -0.00787024413841263 + 0.01294567419606834im … 0.004473735146277639 - 0.008262858990102196im -0.006966307837486409 + 0.01280721402509136im;;; -0.0034954556854979234 - 0.0030311004432965103im 0.0015745123905849259 + 0.001841018329958227im … -0.0010283338182752897 - 0.0005389192405832398im 0.002478076836178343 + 0.0017023798303982536im; 0.009475591233532588 - 0.011087321726911131im -0.007596153007948777 + 0.01023386272754942im … 0.004863067569823226 - 0.0066510355221691im -0.008176924684354062 + 0.010324458136802436im; … ; 0.008297299023404522 - 0.01616624332262066im -0.0074262078675886205 + 0.01406684488047917im … 0.004753068504658969 - 0.00923349981793723im -0.007183513275591491 + 0.014028512213723652im; -0.005350431251632868 + 0.014662511492288137im 0.00555779642245074 - 0.012611175486789403im … -0.0032817140188374167 + 0.00749510878375009im 0.004976856781620086 - 0.012521127880745778im;;; 0.0023604951460976823 + 0.002052596079510728im -0.001455053654745071 - 0.0014296792905891197im … 0.0009285569069565308 + 0.0006390092832232046im -0.0017636474613518124 - 0.0013809182971346408im; -0.00641833717528869 + 0.006927580858579766im 0.005148896664501105 - 0.006006652328297245im … -0.003398650962266207 + 0.004025149846118606im 0.005391206453207554 - 0.006045161228673214im; … ; -0.005169159291879701 + 0.010958254812874355im 0.0046223401499835975 - 0.009566564304227823im … -0.0031185005517163165 + 0.0066181891742811555im 0.004482895945612467 - 0.009544383382045928im; 0.003182135082084794 - 0.009746220313842506im -0.002996089019771678 + 0.008139430137433535im … 0.0018498869884580474 - 0.005090934792125753im -0.002753519907646361 + 0.008101041047058243im;;; … ;;; 0.0011049735402472159 + 0.0009559623590999851im -0.0009981394323082288 - 0.0008001473738802039im … 0.000547212546742522 + 0.0005565530666121477im -0.0008768347034903214 - 0.0008194964146553284im; 0.0013605645904556064 - 0.004800088953238581im -0.0011733101123324826 + 0.004228312801929758im … 0.0009809176250802423 - 0.003054880444264765im -0.0012785181908047142 + 0.0042450517476140355im; … ; 0.00387315908451502 - 0.005088085846130583im -0.00353069681238648 + 0.0046052225683011015im … 0.0025925065823031754 - 0.003486921760984308im -0.0034566412554093216 + 0.004593344903876163im; -0.0032343563858531106 + 0.0031807050115322663im 0.0028998259228736937 - 0.0028461772855664412im … -0.0019626574247689705 + 0.0020577738066392274im 0.0027945387586532605 - 0.0028293192378664837im;;; -0.0024652355018066406 - 0.0021416484339200677im 0.0018585352550804049 + 0.0014520988110678398im … -0.0006631569835589242 - 0.0007469662408303725im 0.0015495259786568723 + 0.0015012344161905965im; -0.00310486654343199 + 0.009809753818380965im 0.0026760282217927894 - 0.00815892548876882im … -0.0020647491104290308 + 0.005186923917247599im 0.0029183896752938802 - 0.008197317294687163im; … ; -0.006760956251044638 + 0.009572798367709988im 0.005956415865649947 - 0.00840538574075064im … -0.0039891707584111535 + 0.0058602831518499935im 0.00581717940075604 - 0.008383144167068087im; 0.006495308829479254 - 0.006864789549446493im -0.005469102925265895 + 0.005988022455594591im … 0.003184068036749297 - 0.003929498731232236im -0.005226630213656902 + 0.00594942740802868im;;; 0.0035822592114612174 + 0.003115407942931677im -0.00257292989992441 - 0.0017635035287055535im … 0.00039461086168533495 + 0.0007004970273268565im -0.0016682183369112849 - 0.0019029758468474241im; 0.005328308408337333 - 0.01469235518304452im -0.004926608968105298 + 0.012556723739023726im … 0.003746252845035762 - 0.007623975411822459im -0.005507603886362039 + 0.01264719905466123im; … ; 0.009720234350994218 - 0.01492761113730533im -0.00854318527995842 + 0.013094527575142064im … 0.005395294097180952 - 0.008674472410609532im -0.008300836249951773 + 0.013056169169464536im; -0.009496565770130406 + 0.011058524364511518im 0.008227390386106828 - 0.010289549847165253im … -0.004398992026439685 + 0.006522484282840939im 0.007646536856797231 - 0.010199102306959359im],)])]), [8.519872390277141e-7 4.6272641197894527e-7 … 2.4177769021654052e-8 4.627264224746026e-7; 4.627264280005028e-7 2.5014470681069435e-7 … 1.3582987909477004e-8 2.5014471040125483e-7; … ; 2.4177771891744408e-8 1.3582984795225674e-8 … 1.5659191951954648e-9 1.3582989700368004e-8; 4.6272641860723357e-7 2.5014469011259716e-7 … 1.3582986324650893e-8 2.5014470649985893e-7;;; 4.6272641491731157e-7 2.5014469616087896e-7 … 1.3582984262075362e-8 2.5014469636969153e-7; 2.5014471083325365e-7 1.3895604581078393e-7 … 1.127997876527865e-8 1.389560412540472e-7; … ; 1.358298813281994e-8 1.1279977236701608e-8 … 7.937746848719362e-9 1.1279981456599074e-8; 2.501446977112876e-7 1.3895602698402317e-7 … 1.1279977919201132e-8 1.3895603647370923e-7;;; 2.4177767124242262e-8 1.3582985887997636e-8 … 1.5659162873031971e-9 1.3582981698131375e-8; 1.3582988823372962e-8 1.127998254608662e-8 … 7.93774293418714e-9 1.1279978290241496e-8; … ; 1.56591932932747e-9 7.937745721881405e-9 … 3.5718536750926076e-8 7.937746092894689e-9; 1.358298464941029e-8 1.1279976965933372e-8 … 7.937741158949245e-9 1.1279975444168911e-8;;; … ;;; 5.7918421217516156e-8 2.2816025487765324e-8 … 5.6203611188313e-10 2.2816028868526024e-8; 2.2816028219524946e-8 3.489279744435589e-9 … 1.3458504715938294e-8 3.4892800904419928e-9; … ; 5.620360012196619e-10 1.345850443372757e-8 … 7.728312302452429e-8 1.3458501938394675e-8; 2.2816029254587012e-8 3.4892794706048958e-9 … 1.3458501373169718e-8 3.4892813056571107e-9;;; 2.417776984869775e-8 1.3582988082395634e-8 … 1.5659185571146235e-9 1.3582986798413348e-8; 1.3582987007663187e-8 1.1279980239246393e-8 … 7.937746916133736e-9 1.1279980312837009e-8; … ; 1.5659193860772836e-9 7.937747226991693e-9 … 3.571853809399687e-8 7.937746952041498e-9; 1.3582987104396873e-8 1.1279980550897241e-8 … 7.937744460251819e-9 1.1279979134947492e-8;;; 4.6272642078086676e-7 2.501447021059815e-7 … 1.3582986084173987e-8 2.501447051815719e-7; 2.501447061204388e-7 1.3895604099954907e-7 … 1.1279980657748845e-8 1.3895604430397838e-7; … ; 1.3582988954221861e-8 1.1279979517852598e-8 … 7.937745942717574e-9 1.127998175132622e-8; 2.501447045723945e-7 1.3895603660663135e-7 … 1.1279978124248585e-8 1.3895604125435344e-7;;;;], Matrix{ComplexF64}[[-0.2257740797936502 - 0.1299732099394576im 0.9030173873807836 + 0.3404984711083705im -4.829051456401186e-7 - 2.5505835810026777e-6im -3.1001089205362336e-6 + 1.0433555830047976e-6im; 0.15318295057596387 + 0.08818408125169264im 0.03979049176267401 + 0.015003700972793462im -0.14278651516490737 - 0.5620534578322365im -0.20148783779732385 - 0.24504410223239964im; … ; -0.04866736783809272 - 0.02801674145401735im -0.01439170755215324 - 0.005426644432830482im -0.0006795827009406535 - 0.009693914756885965im -0.009036230621207666 - 0.002460613677934267im; 0.08726092404241509 + 0.05023420937071913im 0.028511075096699742 + 0.010750597929534882im -0.003918098423357815 + 0.023689700578904804im 0.039388999633826915 + 0.000325052232637029im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239440634084206, [[-0.5541104220105867, -0.010678390671097465, 0.14495664238927924, 0.14495740236351545]], 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.6206944835273916e-7, -6.293616501653438e-7, -2.295120638745721e-7, 4.994026567536676e-7, 7.923312203248644e-7, 5.851071287857125e-7, 6.516970439143957e-7, 1.0189309702321188e-6, 3.6385921474589005e-7, -2.201093797749147e-6  …  2.781018412706599e-6, 1.201695520218961e-6, 2.4456283568633687e-8, 1.0078541361642917e-7, 5.535308561698935e-7, 4.5341736537957554e-7, 4.631957379143081e-9, -1.0437054342704381e-7, 1.1211093275684252e-7, 1.2550862677626904e-7], resid_history = [0.24939209329334172, 0.0037665038925986344, 0.0002850905963752577, 4.682848657253062e-6, 1.1501519588959642e-8, 5.887277005460825e-11, 3.1923296852185617e-8, 1.0854262029532574e-9], converged = true, n_iter = 8, residual_norm = 1.0854262029532574e-9, maxiter = 100, tol = 1.0e-8, s = 0.933314850479688, restart_history = [6], stage = :finalize, krylovdim = 20, y = [2.812742988903544e-7, -3.1303930715391465e-8, 1.0854262029532574e-9, -0.00020432737645477688, 4.6887272768325e-6, -1.1621343899525634e-8, 5.887277005460825e-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}}([[-3.031057683281672e-8, -2.95755257285452e-7, -2.3392361423250184e-7, 7.635761350976216e-7, 1.5124523967018935e-6, 1.2861774088499085e-6, 1.595205678321612e-6, 2.7085475736226544e-6, 1.0233419754081228e-6, -6.442730250003669e-6  …  7.601757151168454e-6, 3.0853919869654207e-6, 5.891071222974204e-8, 2.3175318842924286e-7, 1.2178281338908371e-6, 9.292979766896342e-7, 8.49842021275598e-9, -1.6721310208407652e-7, 1.624943180098513e-7, 2.182407785839372e-7], [1.4454514354445848e-6, 1.5498949104990334e-6, 5.857192399237746e-7, -1.4161172880145786e-6, -2.4419759882844877e-6, -1.937019567117589e-6, -2.313879225878581e-6, -3.83292673744047e-6, -1.4181251627686784e-6, 8.775408411780198e-6  …  -8.541428557146118e-6, -3.42552585195418e-6, -6.442667613939244e-8, -2.493233602294274e-7, -1.2805604171792503e-6, -9.364677193581594e-7, -7.839422347995916e-9, 1.247161483546042e-7, -5.231129133162107e-8, 2.7142673972508447e-7]]), H = [1.127850165081253 0.007228192656270403 … 0.0 0.0; 0.1141552016236386 1.0245424621322063 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1336125462130098 0.11036325421347828 … 0.0 0.0; 0.0 1.0191959272026467 … 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.9257125357500082
Interacting polarizability:     1.773654864824322

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 = 3.28e-10
[ Info: GMRES linsolve in iteration 2; step 1: normres = 4.87e-11
[ Info: GMRES linsolve in iteration 2; step 2: normres = 5.07e-12
┌ Info: GMRES linsolve converged at iteration 2, step 3:
│ * norm of residual = 1.85e-13
└ * number of operations = 12
Non-interacting polarizability: 1.9257125357499647
Interacting polarizability:     1.7736548666171674

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