We compute the polarizability of a Helium atom. The polarizability
is defined as the change in dipole moment
$$
\mu = \int 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
a = 10.
lattice = a * I(3) # cube of ``a`` bohrs
He = ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))
atoms = [He => [[1/2; 1/2; 1/2]]] # Helium at the center of the box
kgrid = [1, 1, 1] # no kpoint sampling for an isolated system
Ecut = 30
tol = 1e-8
# dipole moment of a given density (assuming the current geometry)
function dipole(ρ)
basis = ρ.basis
rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
d = sum(rr .* ρ.real) * basis.model.unit_cell_volume / prod(basis.fft_size)
end;
We first compute the polarizability by finite differences. First compute the dipole moment at rest:
model = model_LDA(lattice, atoms; symmetries=false)
basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
res = self_consistent_field(basis, tol=tol)
μref = dipole(res.ρ)
n Energy Eₙ-Eₙ₋₁ ρout-ρin Diag --- --------------- --------- -------- ---- 1 -2.769852765115 NaN 2.97e-01 8.0 2 -2.771174027634 -1.32e-03 4.98e-02 1.0 3 -2.771211881160 -3.79e-05 2.03e-03 1.0 4 -2.771212976493 -1.10e-06 2.33e-04 2.0 5 -2.771212981528 -5.03e-09 6.89e-05 1.0
-0.00013465564185778244
Then in a small uniform field:
ε = .01
model_ε = model_LDA(lattice, atoms; extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
symmetries=false)
basis_ε = PlaneWaveBasis(model_ε, Ecut; kgrid=kgrid)
res_ε = self_consistent_field(basis_ε, tol=tol)
με = dipole(res_ε.ρ)
n Energy Eₙ-Eₙ₋₁ ρout-ρin Diag --- --------------- --------- -------- ---- 1 -2.769998039447 NaN 3.00e-01 9.0 2 -2.771276504507 -1.28e-03 4.72e-02 1.0 3 -2.771299805056 -2.33e-05 3.86e-03 1.0 4 -2.771300276868 -4.72e-07 7.94e-04 1.0 5 -2.771300362241 -8.54e-08 2.22e-05 2.0 6 -2.771300362351 -1.10e-10 1.52e-05 2.0
0.01761899891902222
polarizability = (με - μref) / ε
println("Reference dipole: $μref")
println("Displaced dipole: $με")
println("Polarizability : $polarizability")
Reference dipole: -0.00013465564185778244 Displaced dipole: 0.01761899891902222 Polarizability : 1.7753654560880003
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).
Now we use linear response to compute this analytically; we refer to standard
textbooks for the formalism. In the following, \chi_0
is the
independent-particle polarizability, and K
the
Hartree-exchange-correlation kernel. We denote with dV_{\rm ext}
an external
perturbing potential (like in this case the uniform electric field). Then:
$$
d\rho = \chi_0 dV = \chi_0 (dV_{\rm ext} + K d\rho),
$$
which implies
$$
d\rho = (1-\chi_0 K)^-1 \chi_0 dV_{\rm ext}.
$$
From this we identify the polarizability operator to be \chi = (1-\chi_0 K)^-1 \chi_0
.
Numerically, we apply \chi
to dV = -x
by solving a linear equation
(the Dyson equation) iteratively.
using KrylovKit
# KrylovKit cannot deal with the density as a 3D array, so we need `vec` to vectorize
# and `devec` to "unvectorize"
devec(arr) = from_real(basis, reshape(arr, size(res.ρ.real)))
# Apply (1- χ0 K) to a vectorized dρ
function dielectric_operator(dρ)
dρ = devec(dρ)
dv = apply_kernel(basis, dρ; ρ=res.ρ)
χ0dv = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dv...)[1]
vec((dρ - χ0dv).real)
end
# dVext is the potential from a uniform field interacting with the dielectric dipole
# of the density.
dVext = from_real(basis, [-a * (r[1] - 1/2) for r in r_vectors(basis)])
# Apply χ0 once to get non-interacting dipole
dρ_nointeract = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dVext)[1]
# Solve Dyson equation to get interacting dipole
dρ = devec(linsolve(dielectric_operator, vec(dρ_nointeract.real), verbosity=3)[1])
println("Non-interacting polarizability: $(dipole(dρ_nointeract))")
println("Interacting polarizability: $(dipole(dρ))")
WARNING: using KrylovKit.basis in module ##276 conflicts with an existing identifier. ┌ Info: GMRES linsolve in iter 1; step 1: normres = 2.500933252371e-01 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:55 ┌ Info: GMRES linsolve in iter 1; step 2: normres = 3.838695181580e-03 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; step 3: normres = 4.151617778933e-04 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; step 4: normres = 7.148429096267e-06 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; step 5: normres = 6.602563423436e-07 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; step 6: normres = 1.926005643034e-09 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; step 7: normres = 1.720440212648e-11 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; step 8: normres = 6.114992993067e-13 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 1; finised at step 8: normres = 6.114992993067e-13 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:96 ┌ Info: GMRES linsolve in iter 2; step 1: normres = 6.907935889496e-07 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:55 ┌ Info: GMRES linsolve in iter 2; step 2: normres = 6.509721767133e-08 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 2; step 3: normres = 1.069002068840e-09 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 2; step 4: normres = 2.456703619498e-11 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 2; step 5: normres = 7.378397421897e-13 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 2; finised at step 5: normres = 7.378397421897e-13 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:96 ┌ Info: GMRES linsolve in iter 3; step 1: normres = 3.835875029306e-11 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:55 ┌ Info: GMRES linsolve in iter 3; step 2: normres = 6.087075747925e-12 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 3; step 3: normres = 8.337624120713e-14 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:89 ┌ Info: GMRES linsolve in iter 3; finised at step 3: normres = 8.337624120713e-14 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:96 ┌ Info: GMRES linsolve converged at iteration 3, step 3: │ * norm of residual = 9.834653476749866e-14 │ * number of operations = 18 └ @ KrylovKit /home/runner/.julia/packages/KrylovKit/OLgKs/src/linsolve/gmres.jl:127 Non-interacting polarizability: 1.9326533851427157 Interacting polarizability: 1.7795297437405155
As expected, the interacting polarizability matches the finite difference result. The non-interacting polarizability is higher.