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

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.017612220545455907
polarizability = (με - μref) / ε

println("Reference dipole:  $μref")
println("Displaced dipole:  $με")
println("Polarizability :   $polarizability")
Reference dipole:  -0.00013457288162335064
Displaced dipole:  0.017612220545455907
Polarizability :   1.7746793427079257

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.003018027583490829 - 0.0013245479308976087im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; -0.12671425107063783 + 0.282258201192861im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; … ; -0.016654641572072866 + 0.03753624576726076im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.03224316393285483 - 0.07423860314006282im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]], δρ = [-3.620693914661668e-7 -2.5020146477788914e-7 … -4.8806423877730155e-8 -2.5020148147514036e-7; -6.293615359232363e-7 -4.793035310189686e-7 … -1.2077840448512214e-7 -4.7930354210812e-7; … ; 1.3626256018115486e-7 1.1356226012277712e-7 … 4.711986856888496e-8 1.1356225176477916e-7; 1.3305468976077914e-7 1.3632635716646362e-7 … 5.2006833266628025e-8 1.3632634273862993e-7;;; -2.5020141798290207e-7 -1.740490026791217e-7 … -3.556492691017343e-8 -1.7404901440106454e-7; -4.793034959978433e-7 -3.679730831133496e-7 … -1.1202890375774285e-7 -3.679730900014889e-7; … ; 1.1356226761352829e-7 1.1211095484666173e-7 … 1.0980679627210997e-7 1.1211094776014105e-7; 1.363263910467783e-7 1.2550865478030473e-7 … 5.094907489301714e-8 1.2550864408729738e-7;;; -4.880639521715101e-8 -3.5564913839172213e-8 … -1.1412137630910528e-8 -3.556491603738274e-8; -1.2077838097455783e-7 -1.1202889276012608e-7 … -9.78622024165076e-8 -1.120288920453475e-7; … ; 4.7119884752868616e-8 1.0980681320061355e-7 … 2.493822676728719e-7 1.0980680908243083e-7; 5.2006857318075074e-8 5.0949089235842944e-8 … 4.910627294118881e-8 5.09490860479703e-8;;; … ;;; 6.47576611411934e-8 4.1103368452898543e-8 … -6.681696635396171e-9 4.1103372244446086e-8; 1.6897561658675537e-7 6.653288465683312e-8 … -1.3294449496346966e-7 6.653288725386871e-8; … ; 3.267761204622742e-8 1.6017361796104965e-7 … 3.847035118701604e-7 1.6017361548195393e-7; -9.141988762753495e-8 -3.583033941764204e-8 … 7.037280120966565e-8 -3.583033816361983e-8;;; -4.880644255630164e-8 -3.5564949053070255e-8 … -1.1412149699356774e-8 -3.5564951501017085e-8; -1.2077841277019955e-7 -1.1202892140918576e-7 … -9.786223130887977e-8 -1.1202892298520128e-7; … ; 4.7119865825729196e-8 1.0980678045268047e-7 … 2.493822053309665e-7 1.0980677568389021e-7; 5.200681710855587e-8 5.094905327854586e-8 … 4.910624349923659e-8 5.094904967475311e-8;;; -2.502015267550469e-7 -1.740490819214958e-7 … -3.556494534013614e-8 -1.740490940201323e-7; -4.793035739638129e-7 -3.6797314051360837e-7 … -1.1202892208522856e-7 -3.6797314875080314e-7; … ; 1.1356224538886904e-7 1.1211093506099133e-7 … 1.0980677982894178e-7 1.1211092755424099e-7; 1.363263092201741e-7 1.2550859416238408e-7 … 5.094905735737615e-8 1.2550858315548387e-7;;;;], δHψ = Matrix{ComplexF64}[[0.0017555483330531224 + 0.000770472675267789im 0.005650995191559602 - 0.09512932703670933im 0.947438381851101 - 1.0281767061018257im -0.0947911322142993 - 0.13476560935964155im; 0.08134024338262161 - 0.18085174428079448im -1.6155882986371501 - 0.19117948738891674im 0.20874958628049706 - 0.22754796551731798im -0.020711796768309584 - 0.02984416883424892im; … ; 0.020002286105170318 - 0.044753342931952385im 0.020652977955976405 + 0.0011748056208364453im -0.01870488463030095 + 0.020993718489030234im 0.005230628153949815 + 0.009804633203923344im; -0.028836083540459613 + 0.06638182727163522im -0.025542545789253343 - 0.001574742100423585im -0.067641146185728 + 0.043048003925208626im 0.014293639979893418 + 0.0975956195371247im]], δVind = [0.001965431840084975 0.002081000463228322 … 0.0031536649263605872 0.0020810006144254926; 0.003889130421976107 0.004779268862654625 … 0.010276670583828827 0.004779269027732959; … ; -0.0061508461588164285 -0.008269601684856957 … -0.017302968210714954 -0.008269601568938875; 0.00023868267358570072 -0.000398697015771846 … -0.0036646940093250126 -0.0003986968486570379;;; 0.002081000069674367 0.002223227609620393 … 0.0034201187515563243 0.0022232277758657965; 0.004779268426276814 0.005734952414797517 … 0.010903596314575846 0.005734952603684982; … ; -0.008269602443198274 -0.009601007399411238 … -0.012623111267017421 -0.009601007267443017; -0.00039869745541574824 -0.0010730132218190413 … -0.004230680845782261 -0.0010730130415742686;;; 0.003153663030589515 0.0034201174073180385 … 0.004832502510190884 0.003420117694898789; 0.010276668588167172 0.010903594918030302 … 0.012272139170601446 0.010903595314178319; … ; -0.01730296995859482 -0.012623112393037877 … -0.009620265637681555 -0.012623112172414616; -0.0036646960318308364 -0.0042306821484654755 … -0.005491109701985602 -0.004230681884278659;;; … ;;; -0.002336562066608367 -0.0028154650146734817 … 0.005693323262190713 -0.002815465249570582; -0.012913093984560552 -0.017841181643067017 … 0.011498342315788222 -0.01784118186030215; … ; -0.025174459058365475 -0.012810882931608876 … -0.008395872234151232 -0.012810882625405471; 0.007597975796544561 0.010225409004848916 … -0.0054497302544481195 0.01022540862875346;;; 0.003153666303169269 0.003420121029800067 … 0.004832507842416596 0.003420121322161829; 0.010276672307327525 0.01090359875908359 … 0.012272142995612577 0.010903599072140065; … ; -0.01730296568943181 -0.012623108738716735 … -0.009620262761783556 -0.012623108537124495; -0.003664692553210889 -0.004230678553652355 … -0.005491105867368094 -0.004230678292505614;;; 0.002081001003361789 0.0022232286526785517 … 0.003420120539446012 0.0022232288216079274; 0.004779269468503995 0.005734953577223548 … 0.010903598081947033 0.005734953760197581; … ; -0.008269600781076292 -0.009601005673304108 … -0.012623109426751964 -0.009601005551124203; -0.000398696414778477 -0.0010730120670950937 … -0.004230679022860937 -0.0010730118840853958;;;;], δρ0 = [-3.60214007179124e-7 -2.488842600277741e-7 … -4.852836440255759e-8 -2.488842599339568e-7; -5.972249785398561e-7 -4.6069778362770493e-7 … -1.1905766374807451e-7 -4.606977825030619e-7; … ; 1.2735787580841967e-7 1.0836337544284322e-7 … 4.640715944584714e-8 1.0836337012494087e-7; 1.0352005095833281e-7 1.195795993704947e-7 … 5.0695053066802476e-8 1.1957959731377725e-7;;; -2.488842602422664e-7 -1.7310671540916288e-7 … -3.535697362223488e-8 -1.7310671488921438e-7; -4.606977829394807e-7 -3.5721703945526266e-7 … -1.1091437472291086e-7 -3.5721703731247455e-7; … ; 1.0836337511487355e-7 1.0846495286051596e-7 … 1.086223443259564e-7 1.0846494845235983e-7; 1.1957959900133727e-7 1.1611145873754825e-7 … 5.0208297429503e-8 1.161114571674643e-7;;; -4.8528364427203535e-8 -3.535697408161302e-8 … -1.1341312028927177e-8 -3.53569734863178e-8; -1.1905766164509212e-7 -1.1091437598452161e-7 … -9.770666757833346e-8 -1.1091437275284868e-7; … ; 4.640716834198351e-8 1.0862235103639704e-7 … 2.48452425367198e-7 1.0862234904906301e-7; 5.069505570363807e-8 5.0208299495588795e-8 … 4.9268577480581404e-8 5.020829884985156e-8;;; … ;;; 6.434570741769006e-8 4.083862644156599e-8 … -6.637927028419204e-9 4.083862664779896e-8; 1.6975036819994534e-7 6.683086415680411e-8 … -1.3340479652847296e-7 6.683086535603076e-8; … ; 3.270997245016314e-8 1.6023521963216232e-7 … 3.8390142245221587e-7 1.6023521988798355e-7; -9.271367562126017e-8 -3.633621569434487e-8 … 7.126308479592329e-8 -3.6336215839917355e-8;;; -4.8528361871261764e-8 -3.535697214118057e-8 … -1.1341311653614493e-8 -3.535697178586948e-8; -1.1905765633845142e-7 -1.1091437098263643e-7 … -9.770666852474308e-8 -1.1091437002035774e-7; … ; 4.6407161956106784e-8 1.0862234627015346e-7 … 2.484524219435689e-7 1.086223436449167e-7; 5.0695052739857136e-8 5.0208297253550525e-8 … 4.926857596785752e-8 5.020829620721665e-8;;; -2.4888425805410553e-7 -1.7310671377173186e-7 … -3.5356973570791606e-8 -1.7310671362457682e-7; -4.606977800323583e-7 -3.572170369982567e-7 … -1.1091437625341745e-7 -3.572170361955513e-7; … ; 1.083633715418461e-7 1.084649499406117e-7 … 1.0862234185651352e-7 1.084649451313759e-7; 1.1957959817505192e-7 1.1611145800831e-7 … 5.020829667340888e-8 1.161114561535842e-7;;;;], δeigenvalues = [[0.0004950045450071426, 0.09679858596811622, 0.031506491579619746, 0.052500089328695286]], δoccupation = [[0.0, 0.0, 0.0, 0.0]], δεF = 0.0, ε = DFTK.DielectricAdjoint{Array{Float64, 4}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}, Float64, Vector{Vector{Float64}}, StaticArraysCore.SVector{3, Float64}}(Hamiltonian(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), HamiltonianBlock[DFTK.DftHamiltonianBlock(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Any[DFTK.FourierMultiplication{Float64, Vector{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.0, 0.19739208802178715, 0.7895683520871486, 1.7765287921960844, 3.1582734083485944, 4.934802200544679, 7.106115168784338, 9.67221231306757, 12.633093633394378, 15.988759129764759  …  20.133992978222288, 16.38354330580833, 13.027877809437953, 10.066996489111146, 7.5008993448279115, 5.329586376588253, 3.5530575843921683, 2.1713129682396586, 1.184352528130723, 0.5921762640653614]), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [0.1603906580688546 0.1602471565738079 … 0.15981840036267833 0.16024715657380792; 0.16024715657380803 0.16010321540768782 … 0.15967315066597979 0.16010321540768782; … ; 0.15981840036267841 0.1596731506659797 … 0.15923916395256382 0.1596731506659797; 0.16024715657380809 0.16010321540768785 … 0.1596731506659798 0.16010321540768785;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.1599588345026577 … 0.15952744714353234 0.1599588345026577; … ; 0.15967315066597973 0.15952744714353229 … 0.15909210945052638 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772;;; 0.1598184003626783 0.15967315066597965 … 0.15923916395256382 0.15967315066597965; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256385 0.15909210945052635 … 0.15865272241895537 0.15909210945052635; 0.1596731506659798 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; … ;;; 0.15910963392556968 0.15896221057980836 … 0.1585217272926462 0.15896221057980836; 0.15896221057980842 0.15881431920357322 … 0.15837242838914242 0.15881431920357322; … ; 0.15852172729264621 0.15837242838914237 … 0.157926332522532 0.15837242838914237; 0.15896221057980844 0.15881431920357328 … 0.1583724283891425 0.15881431920357328;;; 0.1598184003626783 0.15967315066597962 … 0.15923916395256382 0.15967315066597962; 0.15967315066597973 0.15952744714353237 … 0.1590921094505264 0.15952744714353237; … ; 0.15923916395256382 0.15909210945052635 … 0.15865272241895534 0.15909210945052635; 0.15967315066597979 0.15952744714353237 … 0.15909210945052646 0.15952744714353237;;; 0.16024715657380786 0.16010321540768765 … 0.15967315066597965 0.16010321540768768; 0.16010321540768782 0.15995883450265772 … 0.15952744714353237 0.15995883450265772; … ; 0.1596731506659797 0.15952744714353229 … 0.15909210945052635 0.15952744714353229; 0.16010321540768785 0.15995883450265772 … 0.15952744714353242 0.15995883450265772]), DFTK.NonlocalOperator{Float64, Matrix{ComplexF64}, Matrix{Float64}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), Matrix{ComplexF64}(undef, 7809, 0), Matrix{Float64}(undef, 0, 0)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.NoopOperator{Float64}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809)), DFTK.RealSpaceMultiplication{Float64, Array{Float64, 3}}(PlaneWaveBasis(model = Model(lda_x+lda_c_pw, spin_polarization = :none), Ecut = 30.0 Ha, kgrid = MonkhorstPack([1, 1, 1])), KPoint([     0,      0,      0], spin = 1, num. G vectors =  7809), [-0.15475965283607007 -0.15461620271151652 … -0.15418753101093005 -0.15461620271132298; -0.15461620271200818 -0.15447229891135802 … -0.15404229468856584 -0.154472298911164; … ; -0.1541875310099598 -0.15404229468709532 … -0.15360833735774423 -0.15404229468689856; -0.15461620271083346 -0.15447229891018008 … -0.15404229468737657 -0.1544722989099857;;; -0.154616202711674 -0.15447229891102168 … -0.15404229468823033 -0.1544722989108303; -0.15447229891151715 -0.1543279431028819 … -0.153896602841505 -0.15432794310269013; … ; -0.15404229468725258 -0.153896602840023 … -0.15346128702169998 -0.1538966028398285; -0.15447229891033698 -0.15432794310169862 … -0.1538966028403101 -0.15432794310150627;;; -0.15418753101162594 -0.15404229468876351 … -0.1536083373594366 -0.15404229468857292; -0.1540422946892655 -0.1538966028420402 … -0.15346128702374098 -0.15389660284184936; … ; -0.15360833735844567 -0.15346128702223863 … -0.15302191314660915 -0.153461287022045; -0.1540422946880731 -0.15389660284084447 … -0.1534612870225331 -0.15389660284065285;;; … ;;; -0.1534788058661709 -0.15333139197612985 … -0.15289092493882134 -0.15333139197592188; -0.15333139197662776 -0.15318350575462397 … -0.15274162879115574 -0.15318350575441572; … ; -0.15289092493783776 -0.15274162878966513 … -0.1522955396007325 -0.15274162878945333; -0.15333139197542606 -0.15318350575341894 … -0.15274162878993805 -0.15318350575320971;;; -0.15418753101060878 -0.15404229468774905 … -0.15360833735839619 -0.154042294687547; -0.1540422946882418 -0.15389660284101922 … -0.1534612870226937 -0.15389660284081677; … ; -0.15360833735742302 -0.15346128702121897 … -0.1530219131455629 -0.15346128702101344; -0.1540422946870563 -0.15389660283983045 … -0.153461287021493 -0.1538966028396274;;; -0.1546162027111629 -0.15447229891051173 … -0.15404229468770717 -0.1544722989103146; -0.15447229891100245 -0.15432794310236853 … -0.15389660284097828 -0.15432794310217093; … ; -0.1540422946867384 -0.15389660283951026 … -0.15346128702117381 -0.15389660283930984; -0.15447229890982603 -0.15432794310118886 … -0.15389660283978712 -0.1543279431009908]), 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.015163897692835581 -0.012550129163876513 … -0.004966868185714476 -0.012550129137203832; -0.012550129138391016 -0.010362595096620398 … -0.004135158985506569 -0.01036259508067117; … ; -0.004966868216839828 -0.004135159051038787 … -0.0020697801114536634 -0.004135158930208081; -0.012550129160826096 -0.010362595130139506 … -0.004135158930146163 -0.010362595077619309;;; -0.012550129158078958 -0.010362595116191423 … -0.00413515895705667 -0.01036259508679132; -0.010362595083205334 -0.008622039030410974 … -0.003897511349181779 -0.008622039000769745; … ; -0.004135159045697608 -0.0038975114732484224 … -0.0034841069103046314 -0.003897511378229447; -0.010362595129886618 -0.008622039076461262 … -0.0038975113503124017 -0.008622039031444243;;; -0.004966868202993382 -0.004135158998646883 … -0.0020697801849372706 -0.004135158957530548; -0.004135158938566746 -0.0038975113759897066 … -0.003484106889005676 -0.003897511306483103; … ; -0.002069780362067157 -0.0034841070419369836 … -0.005620579967930453 -0.003484107004312104; -0.00413515904892745 -0.0038975114403257856 … -0.003484106954265361 -0.0038975114121695843;;; … ;;; -0.006547898128928834 -0.004876340241869103 … -0.0014866439829344456 -0.0048763402594246864; -0.004876340313776559 -0.002677988991585561 … -0.00412305416997549 -0.0026779890239097244; … ; -0.0014866441172452596 -0.004123054246222489 … -0.00717089522134362 -0.004123054254023492; -0.004876340205914222 -0.0026779888857407083 … -0.004123054210854703 -0.0026779888953839476;;; -0.004966868073178311 -0.004135158886258834 … -0.0020697801475577193 -0.004135158855933645; -0.004135158835459459 -0.0038975112765041115 … -0.0034841069164787243 -0.003897511256136369; … ; -0.002069780177173867 -0.003484106942508771 … -0.005620579919031622 -0.0034841068932716934; -0.0041351588871996895 -0.00389751132037243 … -0.0034841068835737804 -0.0038975112764765792;;; -0.012550129121377762 -0.01036259508499209 … -0.004135158958175901 -0.010362595060115986; -0.010362595060071943 -0.008622039007211571 … -0.0038975113864812188 -0.008622038995070528; … ; -0.0041351589594788895 -0.0038975114063967025 … -0.003484106860939839 -0.00389751130316804; -0.010362595075009785 -0.008622039032137834 … -0.003897511314248385 -0.008622038981783766])], 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.009532892460051045 -0.006919175301585142 … 0.0006640011660338039 -0.006919175274718893; -0.006919175276591168 -0.0047316786002906015 … 0.001495696991907377 -0.004731678584147362; … ; 0.0006640011358787865 0.0014956969278455932 … 0.0035610464833659245 0.001495697048873059; -0.006919175297851465 -0.004731678632631735 … 0.001495697048457081 -0.004731678579917166;;; -0.0069191752959451 -0.004731678619525451 … 0.0014956970206926462 -0.004731678589933945; -0.00473167858703466 -0.002991147630635179 … 0.0017333329528455644 -0.002991147600802187; … ; 0.0014956969330295368 0.0017333328302608746 … 0.002146715518521767 0.0017333329254743332; -0.004731678632535749 -0.002991147675502164 … 0.0017333329529099304 -0.0029911476302927984;;; 0.0006640011480589826 0.0014956969785692486 … 0.003561046408189949 0.0014956970198761813; 0.0014956970381474815 0.001733332925502454 … 0.002146715537779744 0.0017333329951999049; … ; 0.00356104623205102 0.002146715386350735 … 1.0229304415762183e-5 0.0021467154241692375; 0.001495696928979268 0.001733332862362113 … 0.0021467154737280095 0.0017333328907099387;;; … ;;; -0.0009170700695300656 0.000754478361809403 … 0.004144158370890403 0.0007544783444617918; 0.0007544782894040947 0.00295282445736369 … 0.0015077454280111894 0.002952824425247777; … ; 0.0041441582375631916 0.0015077453532547488 … -0.0015401022995441263 0.0015077453456655487; 0.000754478398468165 0.0029528245644136345 … 0.001507745388349753 0.0029528245549796167;;; 0.0006640012788912122 0.0014956970919717365 … 0.003561046446609918 0.0014956971224989856; 0.001495697142278477 0.001733333026009038 … 0.002146715511353969 0.0017333330465792297; … ; 0.0035610464179669362 0.002146715486798604 … 1.0229354360812015e-5 0.0021467155362412116; 0.0014956970917237985 0.0017333329833294905 … 0.002146715545459675 0.0017333330274284012;;; -0.006919175258732785 -0.004731678587816165 … 0.0014956970200965796 -0.004731678562742913; -0.00473167856338657 -0.0029911476069223817 … 0.00173333291607287 -0.0029911475945837464; … ; 0.0014956970197623999 0.0017333328976253232 … 0.0021467155684126945 0.0017333330010544088; -0.004731678577147964 -0.0029911476306689767 … 0.0017333329894969178 -0.0029911475801168447]), 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.4547116913129674e-6 + 1.11343544141163e-6im 0.0017639168249389083 - 0.0008235088952157476im … 0.0012785032216856129 - 0.0005250393208436853im -0.0018388409803593107 + 0.0007854512361497888im; 0.01028296686018038 - 0.02335810894113493im -0.00947797530649062 + 0.019691195150494im … 0.0038295398636852282 - 0.010377044596057157im -0.00736690664315591 + 0.018748445963096133im; … ; 0.01095263398995263 - 0.024824222385082213im -0.008874916906650872 + 0.02076061300180889im … 0.006165875429504142 - 0.013193935995311244im -0.009599866167869509 + 0.02108426213147122im; -0.010233077327933026 + 0.023389192479630977im 0.007388871790675231 - 0.018740315379133655im … -0.005594047388107118 + 0.011095246645484821im 0.009499769401770143 - 0.019683122620049814im;;; 0.005781255546585961 - 0.007971819148486898im -0.004417385074012895 + 0.005154313954526323im … 0.0002789449072852562 - 0.001245044288273663im -0.0023061173591984964 + 0.004212138916347542im; -0.011831824335251151 + 0.02387987653028381im 0.009730463126653169 - 0.018999204666593696im … -0.004192570249292148 + 0.009979191127860191im 0.008371132298200176 - 0.018392766491233576im; … ; -0.008067164320230833 + 0.019323045792547785im 0.006611422613665799 - 0.01630392198886442im … -0.004818536566750759 + 0.010729114183121499im 0.007180648673772083 - 0.01655787061828714im; 0.0050348139634926615 - 0.014552135804313727im -0.004035692818540974 + 0.012367839644078279im … 0.003695176728132999 - 0.00799372687681613im -0.005395038389241358 + 0.012974523370732087im;;; -0.003918770609264249 + 0.00544667226305504im 0.003040058640533109 - 0.003907396906516419im … -0.0008613390735113166 + 0.0016275502827007551im 0.002314945484118964 - 0.0035842151626050657im; 0.007398236927269843 - 0.014477676649605371im -0.006222300873780808 + 0.01191092270508208im … 0.003212873325984704 - 0.0070516128773537605im -0.005653104519194851 + 0.011657261048509338im; … ; 0.004451176152796514 - 0.011223795346509835im -0.0037456642805664677 + 0.009726733693744233im … 0.0029706953886124146 - 0.006915524136495651im -0.0040730255247504715 + 0.009872624768675384im; -0.0020253410495209593 + 0.006994858097098243im 0.0016653366167541151 - 0.006062010421291359im … -0.001777690106375219 + 0.004319922890603775im 0.0022346274253629184 - 0.006315750683863362im;;; … ;;; -0.0017911745188747543 + 0.002377041533304009im 0.0013748000486319464 - 0.0019409831446132353im … -0.0012151027233304728 + 0.0014148375652866117im 0.0016599586068440889 - 0.0020682517323840637im; -0.0006659092038099105 + 0.0029685708239578735im 0.000742171043737319 - 0.0027259855802603083im … -0.0003002534991783254 + 0.0019022667255128864im 0.0004946341537237981 - 0.0026155094931431385im; … ; -0.004089354479046211 + 0.008327428238200502im 0.0035945419091626078 - 0.007458395041851847im … -0.002914698043544204 + 0.0057186077602146175im 0.0037687547508015867 - 0.007536097526732975im; 0.0036996086327181262 - 0.006978423297611663im -0.003115236608079255 + 0.006063232180312104im … 0.0024765209133199577 - 0.004424772886079214im -0.00336273351767498 + 0.006173666293179756im;;; 0.004013781199755766 - 0.005401423814010802im -0.0024072992932664032 + 0.0035409954075219623im … 0.0017559529771547298 - 0.0019516096176414727im -0.0031321401498232613 + 0.003864858725508363im; 0.0019510764546635726 - 0.0070287087608554835im -0.0021564714840399462 + 0.006351382343020748im … 0.0010192280300191235 - 0.004051791486506811im -0.0015871226996045833 + 0.006097050934961534im; … ; 0.007068643655247799 - 0.014802337964521423im -0.0059393931579494845 + 0.012725844777194188im … 0.004402196654276604 - 0.008872870154165283im -0.0062667458563027325 + 0.012871935334508684im; -0.007472708078829453 + 0.014443227279416125im 0.00573145714637272 - 0.011621180318474714im … -0.003971319252295527 + 0.007319508586693399im 0.0063006144961628395 - 0.011875399874384861im;;; -0.005851003196517622 + 0.007936329406805531im 0.0023994435829859926 - 0.00416617055359781im … -0.0020609860045501544 + 0.0019546687966577347im 0.0045101186659985155 - 0.0051093535020192375im; -0.005015389626313508 + 0.014558507417833716im 0.0053415750230204105 - 0.012998567439207477im … -0.002355649923166596 + 0.007467361749660438im 0.003982143592956282 - 0.012391108750589932im; … ; -0.010407352444937654 + 0.022522029476959936im 0.008448404376078479 - 0.01881496533589451im … -0.005874347044671381 + 0.012172679182085934im 0.00901760579104513 - 0.019069092763586516im; 0.011851561901597843 - 0.02387254927595226im -0.008424882164271608 + 0.01836804760064438im … 0.0055320197332019766 - 0.010505279741041672im -0.009783960755957825 + 0.018975309647638005im],)])]), [8.519873393729828e-7 4.6272648789680347e-7 … 2.4177778426301967e-8 4.6272648473109194e-7; 4.6272648487199466e-7 2.50144749047512e-7 … 1.3582994388583538e-8 2.5014474781394853e-7; … ; 2.417777890388238e-8 1.3582995064649269e-8 … 1.565921208178181e-9 1.3582993818087362e-8; 4.6272648753475766e-7 2.501447516400033e-7 … 1.3582993817448407e-8 2.5014474757789374e-7;;; 4.6272648720869217e-7 2.50144750561202e-7 … 1.3582994095074262e-8 2.501447482872928e-7; 2.501447480099425e-7 1.3895607312837326e-7 … 1.1279986917954648e-8 1.3895607160439603e-7; … ; 1.358299500954607e-8 1.1279988044428169e-8 … 7.937752689734035e-9 1.127998718169819e-8; 2.501447516204434e-7 1.3895607549601312e-7 … 1.1279986928220525e-8 1.389560731815041e-7;;; 2.417777869142452e-8 1.3582994524146636e-8 … 1.5659213807342083e-9 1.3582994099964016e-8; 1.3582993904320187e-8 1.1279987161362266e-8 … 7.93775253781514e-9 1.1279986530268762e-8; … ; 1.5659217966779213e-9 7.93775362862608e-9 … 3.5718556980342554e-8 7.937753360260934e-9; 1.3582995042870352e-8 1.1279987745507077e-8 … 7.937753003291519e-9 1.1279987489860291e-8;;; … ;;; 5.791843009079921e-8 2.2816030137747522e-8 … 5.620367857373003e-10 2.2816030396570073e-8; 2.2816031197878825e-8 3.489280094118428e-9 … 1.3458508231979833e-8 3.489280225376382e-9; … ; 5.620369426034033e-10 1.3458509013624483e-8 … 7.728314740434435e-8 1.345850909359623e-8; 2.2816029607669162e-8 3.4892796643249063e-9 … 1.3458508651052231e-8 3.489279703483361e-9;;; 2.4177776699565283e-8 1.3582993364679847e-8 … 1.5659212929603715e-9 1.358299305182979e-8; 1.3582992840605399e-8 1.1279986258070054e-8 … 7.937752733771675e-9 1.127998607314225e-8; … ; 1.5659213625057418e-9 7.937752919436111e-9 … 3.571855599832791e-8 7.937752568243574e-9; 1.3582993374390129e-8 1.12799866563757e-8 … 7.937752499073012e-9 1.1279986257821478e-8;;; 4.6272648285273217e-7 2.5014474814813813e-7 … 1.358299410662334e-8 2.5014474622412386e-7; 2.5014474622072837e-7 1.3895607193560113e-7 … 1.127998725662014e-8 1.389560713113781e-7; … ; 1.3582994120062071e-8 1.1279987437443265e-8 … 7.93775233763142e-9 1.1279986500167936e-8; 2.5014474737607395e-7 1.3895607321716095e-7 … 1.1279986600773791e-8 1.3895607062826054e-7;;;;], Matrix{ComplexF64}[[0.23854987474125183 + 0.10469444726363472im 0.05722805213239462 - 0.9633818369011435im -6.408112863656749e-7 - 6.961850681997812e-7im 8.844600815646789e-7 + 4.269362328223694e-7im; -0.16185105802952476 - 0.07103297391622033im 0.002521697692097065 - 0.042450408715733585im -0.3251993912990726 - 0.30023703281678216im -0.04255408747736807 + 0.030313597003734816im; … ; 0.05142129054284411 + 0.0225677066584641im -0.0009120632903527198 + 0.015353753248063566im -0.0035341392147232747 - 0.0017008923947156062im 0.004151669736961439 - 6.532309160020469e-5im; -0.09219872607288077 - 0.0404640525886742im 0.001806868216751336 - 0.03041697446411511im 0.0021790474667771624 - 0.006646994179880246im -0.025308615881804823 + 0.0019032147457154285im]], [[2.0, 0.0, 0.0, 0.0]], -0.28239441081885897, [[-0.5541104250034026, -0.010678396634315357, 0.1449566048291054, 0.1449566478186669]], 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.620693914661668e-7, -6.293615359232363e-7, -2.295120363607812e-7, 4.994025627956864e-7, 7.923311605795349e-7, 5.851072014664433e-7, 6.516970821347131e-7, 1.0189308279882245e-6, 3.638592037990827e-7, -2.2010931109495195e-6  …  2.781017817131195e-6, 1.201695316410186e-6, 2.4456329841319945e-8, 1.007854069640385e-7, 5.535306972608179e-7, 4.5341719465691167e-7, 4.63189448314093e-9, -1.0437053109484774e-7, 1.1211092755424099e-7, 1.2550858315548387e-7], resid_history = [0.24939209480034968, 0.0037665486281898575, 0.00028527459980011026, 4.680615527282139e-6, 6.998988103738359e-8, 4.397347887014917e-11, 4.6069331589019985e-8, 1.2460091934708189e-9], converged = true, n_iter = 8, residual_norm = 1.2460091934708189e-9, maxiter = 100, tol = 1.0e-8, s = 0.9356383556466076, residual = [7.361027840123217e-7, 5.241753664162266e-7, 1.167888821798706e-7, -1.7824177369234426e-7, -2.175222161035164e-7, -1.329192292825965e-7, -1.2899275016830414e-7, -1.811250598693351e-7, -5.912523150634532e-8, 3.323953277472511e-7  …  5.432774966364475e-7, 2.442560828548005e-7, 5.362266650399169e-9, 2.4193236510694444e-8, 1.4551819281859016e-7, 1.333221940362406e-7, 1.6284380096962246e-9, -4.937102067535835e-8, 8.757983075159966e-8, 3.0144608793030463e-7], restart_history = [6], stage = :finalize, krylovdim = 20, y = [3.3728764083494196e-7, -4.554251121592558e-8, 1.2460091934708189e-9, -0.00020446038061474315, 4.6836354521620965e-6, -7.065662195590398e-8, 4.397347887014917e-11, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], V = KrylovKit.OrthonormalBasis{Vector{Float64}}([[-2.39326082603623e-7, -1.848441604006086e-6, -6.506217501486605e-7, 1.058340783022244e-6, 1.2807765419445868e-6, 9.158334361787903e-7, 1.1136259233281684e-6, 1.6615035199386447e-6, 4.5386890521754974e-7, -2.033110521248343e-6  …  3.054060835123104e-6, 1.2022639931871255e-6, 2.8785119654949834e-8, 1.426145257814723e-7, 7.823742203486296e-7, 5.421031807264948e-7, 4.986294344383606e-9, -1.3834141376376794e-7, 2.276977888687668e-7, 4.760368905576313e-7], [6.057802096926078e-7, 1.9118203331976073e-6, 6.463323925443234e-7, -1.0615766399601442e-6, -1.295494952470666e-6, -9.080596919232752e-7, -1.0797136373995616e-6, -1.6101144126445375e-6, -4.495249751326008e-7, 2.0568382586798795e-6  …  -2.337720521060516e-6, -8.878353583014889e-7, -2.118651321268581e-8, -1.062027479431166e-7, -5.760702699703716e-7, -3.766242920603828e-7, -3.14315869691133e-9, 8.383503428682981e-8, -1.3725717664470867e-7, -2.2662068707619364e-7]]), H = [1.1166064302066039 0.009953614676965525 … 0.0 0.0; 0.13778800419979734 1.0197229635733132 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], R = [1.1250757548183583 0.13476413609503723 … 0.0 0.0; 0.0 1.0111976096714306 … 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.9257125155267059
Interacting polarizability:     1.7736548468557054

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

Manual solution of the Dyson equations

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

using KrylovKit

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

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

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

println("Non-interacting polarizability: $(dipole(basis, δρ_nointeract))")
println("Interacting polarizability:     $(dipole(basis, δρ))")
WARNING: using KrylovKit.basis in module Main conflicts with an existing identifier.
[ Info: GMRES linsolve starts with norm of residual = 4.19e+00
[ Info: GMRES linsolve in iteration 1; step 1: normres = 2.49e-01
[ Info: GMRES linsolve in iteration 1; step 2: normres = 3.77e-03
[ Info: GMRES linsolve in iteration 1; step 3: normres = 2.85e-04
[ Info: GMRES linsolve in iteration 1; step 4: normres = 4.69e-06
[ Info: GMRES linsolve in iteration 1; step 5: normres = 1.09e-08
[ Info: GMRES linsolve in iteration 1; step 6: normres = 6.28e-11
[ Info: GMRES linsolve in iteration 1; step 7: normres = 1.45e-09
[ Info: GMRES linsolve in iteration 2; step 1: normres = 7.42e-11
┌ Info: GMRES linsolve converged at iteration 2, step 2:
* norm of residual = 2.10e-12
* number of operations = 11
Non-interacting polarizability: 1.9257125155265888
Interacting polarizability:     1.7736548398177938

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