We compute the polarizability of a Helium atom. The polarizability is defined as the change in dipole moment μ=∫rρ(r)dr
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;
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)
res = self_consistent_field(basis; tol)
μref = dipole(basis, res.ρ)
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -2.770417161655 -0.52 9.0 200ms 2 -2.771692104137 -2.89 -1.32 1.0 96.4ms 3 -2.771714268741 -4.65 -2.43 1.0 508ms 4 -2.771714637491 -6.43 -3.13 1.0 110ms 5 -2.771714715159 -7.11 -4.49 2.0 144ms 6 -2.771714715221 -10.21 -4.65 1.0 130ms 7 -2.771714715246 -10.59 -5.55 1.0 103ms 8 -2.771714715250 -11.45 -6.30 2.0 115ms 9 -2.771714715250 -13.24 -7.24 1.0 101ms 10 -2.771714715250 + -15.05 -8.04 2.0 117ms
-0.0001345736247947825
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_ε.ρ)
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -2.770499848770 -0.52 9.0 223ms 2 -2.771779455241 -2.89 -1.32 1.0 150ms 3 -2.771801607532 -4.65 -2.43 1.0 93.7ms 4 -2.771802003727 -6.40 -3.14 1.0 97.3ms 5 -2.771802074301 -7.15 -4.48 2.0 517ms 6 -2.771802074425 -9.91 -4.64 1.0 96.0ms 7 -2.771802074474 -10.30 -5.55 1.0 98.8ms 8 -2.771802074476 -11.94 -5.79 1.0 119ms 9 -2.771802074476 -12.30 -6.48 2.0 115ms 10 -2.771802074476 + -13.64 -6.75 1.0 103ms 11 -2.771802074476 -13.62 -7.06 1.0 101ms 12 -2.771802074476 + -13.89 -7.60 1.0 370ms 13 -2.771802074476 + -Inf -7.83 1.0 126ms 14 -2.771802074476 -14.12 -8.64 1.0 117ms
0.017612221070465728
polarizability = (με - μref) / ε
println("Reference dipole: $μref")
println("Displaced dipole: $με")
println("Polarizability : $polarizability")
Reference dipole: -0.0001345736247947825 Displaced dipole: 0.017612221070465728 Polarizability : 1.7746794695260508
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 (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 δVext an external perturbing potential (like in this case the uniform electric field). Then: δρ=χ0δV=χ0(δVext+Kδρ),
using KrylovKit
# Apply $(1- χ_0 K)$
function dielectric_operator(δρ)
δV = apply_kernel(basis, δρ; res.ρ)
χ0δV = apply_χ0(res, δV)
δρ - χ0δV
end
# `δ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)
# Apply $χ_0$ once to get non-interacting dipole
δρ_nointeract = apply_χ0(res, δ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 ##281 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 = 2.19e-08 [ Info: GMRES linsolve in iteration 2; step 1: normres = 1.76e-09 [ Info: GMRES linsolve in iteration 2; step 2: normres = 3.21e-11 [ Info: GMRES linsolve in iteration 2; step 3: normres = 4.92e-12 ┌ Info: GMRES linsolve converged at iteration 2, step 4: │ * norm of residual = 6.38e-14 └ * number of operations = 13 Non-interacting polarizability: 1.9257125373920252 Interacting polarizability: 1.7736548592813037
As expected, the interacting polarizability matches the finite difference result. The non-interacting polarizability is higher.