/
index.html
80 lines (70 loc) · 14.4 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Polarizability by linear response · DFTK.jl</title><link rel="canonical" href="https://juliamolsim.github.io/DFTK.jl/stable/examples/polarizability/"/><link href="https://fonts.googleapis.com/css?family=Lato|Roboto+Mono" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.0/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.0/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.0/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.11.1/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script><link href="../../assets/favicon.ico" rel="icon" type="image/x-icon"/></head><body><div id="documenter"><nav class="docs-sidebar"><a class="docs-logo" href="../../"><img src="../../assets/logo.png" alt="DFTK.jl logo"/></a><div class="docs-package-name"><span class="docs-autofit">DFTK.jl</span></div><form class="docs-search" action="../../search/"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="../../">Home</a></li><li><span class="tocitem">Getting started</span><ul><li><a class="tocitem" href="../../guide/installation/">Installation</a></li><li><a class="tocitem" href="../../guide/tutorial/">Tutorial</a></li><li><a class="tocitem" href="../../guide/parallelization/">Timings and parallelization</a></li><li><a class="tocitem" href="../../guide/density_functional_theory/">Density-functional theory</a></li></ul></li><li><span class="tocitem">Examples</span><ul><li><a class="tocitem" href="../metallic_systems/">Temperature and metallic systems</a></li><li><a class="tocitem" href="../pymatgen/">Creating supercells with pymatgen</a></li><li><a class="tocitem" href="../ase/">Creating slabs with ASE</a></li><li><a class="tocitem" href="../collinear_magnetism/">Collinear spin and magnetic systems</a></li><li><a class="tocitem" href="../geometry_optimization/">Geometry optimization</a></li><li><a class="tocitem" href="../scf_callbacks/">Monitoring self-consistent field calculations</a></li><li><a class="tocitem" href="../scf_checkpoints/">Saving SCF results on disk and SCF checkpoints</a></li><li class="is-active"><a class="tocitem" href>Polarizability by linear response</a><ul class="internal"><li><a class="tocitem" href="#Polarizability-by-finite-differences"><span>Polarizability by finite differences</span></a></li><li><a class="tocitem" href="#Polarizability-by-linear-response-2"><span>Polarizability by linear response</span></a></li></ul></li><li><a class="tocitem" href="../gross_pitaevskii/">Gross-Pitaevskii equation in one dimension</a></li><li><a class="tocitem" href="../gross_pitaevskii_2D/">Gross-Pitaevskii equation with magnetism</a></li><li><a class="tocitem" href="../cohen_bergstresser/">Cohen-Bergstresser model</a></li><li><a class="tocitem" href="../arbitrary_floattype/">Arbitrary floating-point types</a></li><li><a class="tocitem" href="../custom_solvers/">Custom solvers</a></li><li><a class="tocitem" href="../custom_potential/">Custom potential</a></li></ul></li><li><span class="tocitem">Advanced topics</span><ul><li><a class="tocitem" href="../../advanced/conventions/">Notation and conventions</a></li><li><a class="tocitem" href="../../advanced/data_structures/">Data structures</a></li><li><a class="tocitem" href="../../advanced/useful_formulas/">Useful formulas</a></li><li><a class="tocitem" href="../../advanced/symmetries/">Crystal symmetries</a></li></ul></li><li><a class="tocitem" href="../../api/">API reference</a></li><li><a class="tocitem" href="../../publications/">Publications</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">Examples</a></li><li class="is-active"><a href>Polarizability by linear response</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>Polarizability by linear response</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/polarizability.jl" title="Edit on GitHub"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitHub</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Polarizability-by-linear-response"><a class="docs-heading-anchor" href="#Polarizability-by-linear-response">Polarizability by linear response</a><a id="Polarizability-by-linear-response-1"></a><a class="docs-heading-anchor-permalink" href="#Polarizability-by-linear-response" title="Permalink"></a></h1><p><a href="https://mybinder.org/v2/gh/JuliaMolSim/DFTK.jl/gh-pages?filepath=v0.2.6/examples/polarizability.ipynb"><img src="https://mybinder.org/badge_logo.svg" alt/></a> <a href="https://nbviewer.jupyter.org/github/JuliaMolSim/DFTK.jl/blob/gh-pages/v0.2.6/examples/polarizability.ipynb"><img src="https://img.shields.io/badge/show-nbviewer-579ACA.svg" alt/></a></p><p>We compute the polarizability of a Helium atom. The polarizability is defined as the change in dipole moment</p><p class="math-container">\[\mu = \int r ρ(r) dr\]</p><p>with respect to a small uniform electric field <span>$E = -x$</span>.</p><p>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).</p><p>As in other tests, this is not fully converged, convergence parameters were simply selected for fast execution on CI,</p><pre><code class="language-julia">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;</code></pre><h2 id="Polarizability-by-finite-differences"><a class="docs-heading-anchor" href="#Polarizability-by-finite-differences">Polarizability by finite differences</a><a id="Polarizability-by-finite-differences-1"></a><a class="docs-heading-anchor-permalink" href="#Polarizability-by-finite-differences" title="Permalink"></a></h2><p>We first compute the polarizability by finite differences. First compute the dipole moment at rest:</p><pre><code class="language-julia">model = model_LDA(lattice, atoms; symmetries=false)
basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
res = self_consistent_field(basis, tol=tol)
μref = dipole(res.ρ)</code></pre><pre class="documenter-example-output">-0.00013528987861585455</pre><p>Then in a small uniform field:</p><pre><code class="language-julia">ε = .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_ε.ρ)</code></pre><pre class="documenter-example-output">0.01761502802925802</pre><pre><code class="language-julia">polarizability = (με - μref) / ε
println("Reference dipole: $μref")
println("Displaced dipole: $με")
println("Polarizability : $polarizability")</code></pre><pre class="documenter-example-output">Reference dipole: -0.00013528987861585455
Displaced dipole: 0.01761502802925802
Polarizability : 1.7750317907873874</pre><p>The result on more converged grids is very close to published results. For example <a href="https://doi.org/10.1039/C8CP03569E">DOI 10.1039/C8CP03569E</a> quotes <strong>1.65</strong> with LSDA and <strong>1.38</strong> with CCSD(T).</p><h2 id="Polarizability-by-linear-response-2"><a class="docs-heading-anchor" href="#Polarizability-by-linear-response-2">Polarizability by linear response</a><a class="docs-heading-anchor-permalink" href="#Polarizability-by-linear-response-2" title="Permalink"></a></h2><p>Now we use linear response to compute this analytically; we refer to standard textbooks for the formalism. In the following, <span>$\chi_0$</span> is the independent-particle polarizability, and <span>$K$</span> the Hartree-exchange-correlation kernel. We denote with <span>$dV_{\rm ext}$</span> an external perturbing potential (like in this case the uniform electric field). Then:</p><p class="math-container">\[d\rho = \chi_0 dV = \chi_0 (dV_{\rm ext} + K d\rho),\]</p><p>which implies</p><p class="math-container">\[d\rho = (1-\chi_0 K)^-1 \chi_0 dV_{\rm ext}.\]</p><p>From this we identify the polarizability operator to be <span>$\chi = (1-\chi_0 K)^-1 \chi_0$</span>. Numerically, we apply <span>$\chi$</span> to <span>$dV = -x$</span> by solving a linear equation (the Dyson equation) iteratively.</p><pre><code class="language-julia">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ρ))")</code></pre><pre class="documenter-example-output">WARNING: using KrylovKit.basis in module ex-polarizability conflicts with an existing identifier.
[ Info: GMRES linsolve in iter 1; step 1: normres = 2.491181012722e-01
[ Info: GMRES linsolve in iter 1; step 2: normres = 3.729196710082e-03
[ Info: GMRES linsolve in iter 1; step 3: normres = 1.931096354930e-04
[ Info: GMRES linsolve in iter 1; step 4: normres = 5.336595719728e-06
[ Info: GMRES linsolve in iter 1; step 5: normres = 4.522772099275e-07
[ Info: GMRES linsolve in iter 1; step 6: normres = 1.856475240484e-09
[ Info: GMRES linsolve in iter 1; step 7: normres = 1.136033591874e-11
[ Info: GMRES linsolve in iter 1; step 8: normres = 2.302135562488e-13
[ Info: GMRES linsolve in iter 1; finised at step 8: normres = 2.302135562488e-13
[ Info: GMRES linsolve in iter 2; step 1: normres = 1.218980046467e-06
[ Info: GMRES linsolve in iter 2; step 2: normres = 6.920576273222e-08
[ Info: GMRES linsolve in iter 2; step 3: normres = 2.817025246420e-09
[ Info: GMRES linsolve in iter 2; step 4: normres = 2.180354369470e-11
[ Info: GMRES linsolve in iter 2; step 5: normres = 6.788532734724e-13
[ Info: GMRES linsolve in iter 2; finised at step 5: normres = 6.788532734724e-13
[ Info: GMRES linsolve in iter 3; step 1: normres = 2.313511677849e-12
[ Info: GMRES linsolve in iter 3; finised at step 1: normres = 2.313511677849e-12
┌ Info: GMRES linsolve converged at iteration 3, step 1:
│ * norm of residual = 2.3146282361139995e-12
└ * number of operations = 16
Non-interacting polarizability: 1.9262816188532756
Interacting polarizability: 1.7740221748801412</pre><p>As expected, the interacting polarizability matches the finite difference result. The non-interacting polarizability is higher.</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../scf_checkpoints/">« Saving SCF results on disk and SCF checkpoints</a><a class="docs-footer-nextpage" href="../gross_pitaevskii/">Gross-Pitaevskii equation in one dimension »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> on <span class="colophon-date" title="Monday 1 February 2021 15:21">Monday 1 February 2021</span>. Using Julia version 1.5.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>