Sternheimer¶
Preparing the ground state¶
In this tutorial, we are going to calculate the dielectric constant of diamond-structure silicon, using the linear-response Sternheimer equation \(^4\) (aka density-functional perturbation theory \(^2\)) and the quantum theory of polarization \(^3\). Like in the tutorial on dynamic polarizabilities of a molecule, we will use an electric field as a perturbation, but unlike in the case of a finite system, we will need to use a different form for the electric field due to the peculiar properties of electric fields in periodic systems. An electric field \(\vec{\mathcal{E}}\) gives rise to a term in the Hamiltonian \(\vec{\mathcal{E}} \cdot \vec{r}\), but such a perturbation is not periodic. If we use such a perturbation, we lose the ability to use Bloch’s Theorem ( \(\psi_{\vec{k}} \left( \vec{r} \right) = e^{i \vec{k} \cdot \vec{r}} u_{\vec{k}} \left( \vec{r} \right)\) ) and all the other convenient related formalism for periodic systems. Instead, we can express the Hamiltonian term as \(\vec{\mathcal{E}} \cdot i \vec{\nabla}_{\vec{k}}\), which is now no longer a local potential but instead an operator that applies to the periodic part of the wavefunction \(u_{\vec{k}} \left( \vec{r} \right)\). The general theory is laid out in Refs. \(^4\) and \(^5\).
While some codes such as Quantum ESPRESSO use this formalism to calculate dielectric constants, or Raman intensities (related to \(\frac{\partial \epsilon}{\partial R}\) for an atomic position \(R\)), they evaluate the derivative with respect to \(k\) by finite differences, i.e. \(\frac{\partial}{\partial k_i} u_{\vec{k}} \approx \left( u_{\vec{k} + \Delta \vec{k}_i} - u_{\vec{k}} \right) / \Delta k_i \). In Octopus, instead we use linear response to evaluate that derivative, via the same Sternheimer equation we would use for the electric field perturbation. The perturbation in this case, to the effective Bloch Hamiltonian \(H_{\vec{k}} = \frac{\hbar^2}{2m} (i \nabla + \vec{k})^2 + V(\vec{r}, \vec{r'})\) is \(\frac{\partial H_{\vec k}}{\partial k_i} = -i \frac{\partial}{\partial k_i} + k_i + \left[ V_{\rm nl}, r_i \right] \). The main contribution is from the kinetic energy term, but the potential can contribute as well insofar as it is nonlocal, in which case it is effectively a function of \(k\) too. The Hartree potential is local. The exchange-correlation potential is local in the Kohn-Sham framework, i.e. LDA, GGA, and meta-GGA functionals, but it is nonlocal when hybrid functionals are used. Moreover, when pseudopotentials are used, as we do in Octopus, these are almost always non-local (angular-momentum-dependent) and therefore there is a non-zero commutator to include. This perturbation is referred to as \(k \cdot p\) in Octopus, as it is the same one used in the classic \(k \cdot p\) perturbation theory for semiconductor effective masses. Indeed, from our \(k \cdot p\) run, we will be able to obtain group velocities and effective masses – albeit with some complication in both cases whenever there is a degeneracy, in which case a degenerate perturbation theory needs to be used to consider properly the way that bands split when the \(k\)-point changes from the symmetric point. One simplification of using this perturbation compared to other ones is that since this is a non-physical perturbation, i.e. it does not constitute a change of the physical system, but only of an internal parameter we use to describe states, it should not be solved self-consistently, and so the calculation is faster.
Once we have obtained \(\frac{u_{n \vec{k}}}{\partial k_i}\) for each band \(n\), k-point \(\vec{k}\), and Cartesian direction \(i\), we can now calculate the electric field response as we would for a finite system, and obtain a polarizability \(\alpha_{ij}\). We generally prefer to think in periodic systems of the intensive quantity \(\epsilon_{ij} = 1 + \frac{4 \pi}{V} \alpha_{ij}\). In the simplest case, we can find the static dielectric constant, i.e. the response at frequency \(\omega = 0\). Note that this includes only the electronic contribution, not any ionic contributions, and thus is referred to as \(\epsilon_\infty\) since the frequency is infinite with respect to phonons (but is very small as far as the electrons are concerned). We can also extend this approach straightforwardly to the dielectric function at nonzero frequencies \(\omega\), in which case we are indeed using TDDFT. To prevent divergences at resonant frequencies, and allow us to obtain the imaginary part \(\epsilon_2\) (related to optical absorption) as well as the real part \(\epsilon_1\), we introduce a small imaginary part of the frequency \(i \eta\), which will broaden out these resonances into a Lorentzian lineshape.
Now that we have laid out the theory, let’s do some calculations. We begin with a ground-state SCF run as usual, with the input file below. We ask for somewhat higher precision than usual, which helps with Sternheimer calculations’ convergence later. We use FilterPotentials as this is sometimes helpful for convergence in Sternheimer.
mkdir -p 7-Sternheimer
%%writefile 7-Sternheimer/inp
stdout = 'stdout_gs.txt'
stderr = 'stderr_gs.txt'
FromScratch = yes
CalculationMode = gs
PeriodicDimensions = 3
a = 10.2
BoxShape = parallelepiped
%LatticeParameters
a | a | a
90 | 90 | 90
%
Spacing = a/14
%ReducedCoordinates
"Si" | 0.0 | 0.0 | 0.0
"Si" | 1/2 | 1/2 | 0.0
"Si" | 1/2 | 0.0 | 1/2
"Si" | 0.0 | 1/2 | 1/2
"Si" | 1/4 | 1/4 | 1/4
"Si" | 1/4 + 1/2 | 1/4 + 1/2 | 1/4
"Si" | 1/4 + 1/2 | 1/4 | 1/4 + 1/2
"Si" | 1/4 | 1/4 + 1/2 | 1/4 + 1/2
%
%KPointsGrid
4 | 4 | 4
1/2 | 1/2 | 1/2
%
KPointsUseSymmetries = yes
ExperimentalFeatures = yes
SymmetrizeDensity = yes
Eigensolver = chebyshev_filter
ExtraStates = 4
ConvRelDens = 1e-7
%Output
dos
%
FilterPotentials = filter_none
MixField = density
Writing 7-Sternheimer/inp
!cd 7-Sternheimer && mpirun -n 4 octopus
Performing a k dot calculation¶
Next, we perform the kdotp run, based on the ground-state results. We need the Calculationmode = kdotp. The ExtraStates is not needed but we will be able to see the velocity and effective mass of the lowest conduction band in this way.
%%writefile 7-Sternheimer/inp
stdout = 'stdout_kdot.txt'
stderr = 'stderr_kdot.txt'
FromScratch = yes
CalculationMode = kdotp
ExperimentalFeatures = yes
PeriodicDimensions = 3
a = 10.2
BoxShape = parallelepiped
%LatticeParameters
a | a | a
%
Spacing = a/14
%ReducedCoordinates
"Si" | 0.0 | 0.0 | 0.0
"Si" | 1/2 | 1/2 | 0.0
"Si" | 1/2 | 0.0 | 1/2
"Si" | 0.0 | 1/2 | 1/2
"Si" | 1/4 | 1/4 | 1/4
"Si" | 1/4 + 1/2 | 1/4 + 1/2 | 1/4
"Si" | 1/4 + 1/2 | 1/4 | 1/4 + 1/2
"Si" | 1/4 | 1/4 + 1/2 | 1/4 + 1/2
%
%KPointsGrid
4 | 4 | 4
1/2 | 1/2 | 1/2
%
KPointsUseSymmetries = yes
SymmetrizeDensity = yes
Eigensolver = chebyshev_filter
ExtraStates = 4
FilterPotentials = filter_none
MixField = density
Overwriting 7-Sternheimer/inp
!cd 7-Sternheimer && mpirun -n 4 octopus
Look at the file kdotp/velocity to find the band velocities for each band and k-point. This is not the typical approach to finding velocities. How else might you obtain them? They play a key role in Boltzmann transport calculations of electric conductivity for example, which you can do in the BoltzTraP code.
!head -n 47 "7-Sternheimer/kdotp/velocity" && echo "..."
# Band velocities
# spin = 1, k-point = 1
# state energy vg(x) vg(y) vg(z)
# [H] [bH(2pi/h)] [bH(2pi/h)] [bH(2pi/h)]
1 -0.28036 -0.06523 -0.06523 -0.06523
2 -0.16639 0.11101 0.11101 0.11101
3 -0.16639 -0.02579 0.19774 0.16109
4 -0.16639 0.24782 0.02428 0.06094
5 -0.09649 -0.14116 -0.14116 -0.14116
6 -0.09649 -0.18955 -0.24124 0.00733
7 -0.09649 -0.09276 -0.04107 -0.28964
8 0.02227 0.12240 0.12240 0.12240
9 0.02227 0.07759 0.06080 0.22881
10 0.02227 0.16721 0.18400 0.01599
11 0.05121 -0.05264 -0.05026 -0.05457
12 0.05121 -0.05234 -0.05472 -0.05040
13 0.05121 -0.05249 -0.05249 -0.05249
14 0.09276 0.36403 0.36403 0.36403
15 0.13943 0.15782 -0.03556 0.10297
16 0.13943 -0.00766 0.18571 0.04718
17 0.17298 -0.07293 -0.07293 -0.07293
18 0.17298 -0.07556 -0.07359 -0.06963
19 0.17298 -0.07029 -0.07226 -0.07622
20 0.23956 0.00324 0.00324 0.00324
# spin = 1, k-point = 2
# state energy vg(x) vg(y) vg(z)
# [H] [bH(2pi/h)] [bH(2pi/h)] [bH(2pi/h)]
1 -0.26041 -0.19283 -0.06232 -0.06232
2 -0.22149 0.31010 -0.05416 -0.05416
3 -0.15957 -0.04015 0.17306 0.17306
4 -0.15957 -0.04015 0.17306 0.17306
5 -0.10013 0.02859 -0.11912 -0.11912
6 -0.10013 0.02859 -0.11912 -0.11912
7 -0.02789 -0.42680 0.15952 0.15952
8 -0.00140 0.06466 -0.02245 -0.02245
9 -0.00140 0.06466 -0.02245 -0.02245
10 0.03390 0.38577 0.32121 0.32121
11 0.03955 0.10174 -0.20850 -0.20850
12 0.03955 0.10174 -0.20850 -0.20850
13 0.05106 -0.28279 0.09502 0.09502
14 0.06878 -0.16812 -0.04396 -0.04396
15 0.10165 0.25030 -0.02688 -0.02688
16 0.10452 0.37171 -0.04626 -0.04626
17 0.19233 -0.17917 -0.13254 -0.13254
18 0.22500 0.22647 -0.16822 -0.16822
19 0.22593 -0.51925 0.12296 0.12296
20 0.22593 -0.51925 0.12296 0.12296
...
Then, look at the effective masses in the kdotp/kpoint_1_1 etc.
!head -n 25 "7-Sternheimer/kdotp/kpoint_1_1" && echo "..."
# spin index = 1
# k-point index = 1
# k-point coordinates = -0.07699982 -0.07699982 -0.07699982
# Inverse effective-mass tensors
State #1, Energy = -0.28035994 H
0.884358 -0.006741 -0.006741
-0.006741 0.884358 -0.006741
-0.006741 -0.006741 0.884358
Isotropic average 0.884358
State #2, Energy = -0.16638868 H
0.309220 0.114836 0.114836
0.114836 0.309220 0.114836
0.114836 0.114836 0.309220
Isotropic average 0.309220
State #3, Energy = -0.16638869 H
0.399703 0.057993 0.012963
0.057993 0.300196 0.277229
0.012963 0.277229 0.316513
Isotropic average 0.338804
State #4, Energy = -0.16638869 H
...
These are \(3 \times 3\) tensors defined from the band energies \(\epsilon_{n \vec{k}}\). The inverse effective mass is \(m^{-1}\_{ij} \left( n, \vec{k} \right) = \frac{1}{\hbar^2} \frac{\partial^2 \epsilon\_{n \vec{k}}}{\partial k\_i \partial k\_j}\). The effective mass tensor itself is obtained by matrix inversion of \(m^{-1}\). Note that the complication of degenerate perturbation theory is not dealt with currently in the code (as pointed out in a warning), so focus on either the non-degenerate states, or look at invariants of a degenerate subspace, in which you trace over the degenerate states. The end of the main output file classifies the degenerate subspaces for your convenience.
Calculating the linear reponse¶
Finally, we do a em_resp run.
%%writefile 7-Sternheimer/inp
stdout = 'stdout_em_resp.txt'
stderr = 'stderr_em_resp.txt'
FromScratch = yes
CalculationMode = em_resp
ExperimentalFeatures = yes
PeriodicDimensions = 3
a = 10.2
BoxShape = parallelepiped
%LatticeParameters
a | a | a
%
Spacing = a/14
%ReducedCoordinates
"Si" | 0.0 | 0.0 | 0.0
"Si" | 1/2 | 1/2 | 0.0
"Si" | 1/2 | 0.0 | 1/2
"Si" | 0.0 | 1/2 | 1/2
"Si" | 1/4 | 1/4 | 1/4
"Si" | 1/4 + 1/2 | 1/4 + 1/2 | 1/4
"Si" | 1/4 + 1/2 | 1/4 | 1/4 + 1/2
"Si" | 1/4 | 1/4 + 1/2 | 1/4 + 1/2
%
%KPointsGrid
4 | 4 | 4
1/2 | 1/2 | 1/2
%
KPointsUseSymmetries = yes
SymmetrizeDensity = yes
Eigensolver = chebyshev_filter
ExtraStates = 4
%EMFreqs
2 | 0.0 | 0.1
%
EMEta = 0.01
FilterPotentials = filter_none
MixField = density
RestartFixedOccupations = no
Overwriting 7-Sternheimer/inp
!cd 7-Sternheimer && mpirun -n 4 octopus
Look in em_resp/freq_0.0000/epsilon to find the static dielectric constant.
!cat "7-Sternheimer/em_resp/freq_0.0000/epsilon"
# Real part of dielectric constant
14.298823 -0.000000 0.000000
-0.000000 14.298823 -0.000000
0.000000 0.000000 14.298823
Isotropic average 14.298823
# Imaginary part of dielectric constant
0.000000 0.000000 0.000000
0.000000 0.000000 0.000000
0.000000 0.000000 0.000000
Isotropic average 0.000000
You should obtain a diagonal isotropic tensor (\(\epsilon\_{xx} = \epsilon\_{yy} = \epsilon\_{zz}\)) as required by the \(T\_d\) point group of silicon. The value obtained is somewhat overestimated compared to experiment. Sometimes you will hear people attribute this discrepancy to the “band-gap” problem of Kohn-Sham DFT, thinking of a simple perturbation theory expression for the dielectric constant in which a too-small energy denominator will give a too-large result. This reasoning is a misunderstanding however, since the eigenvalues were not used in the calculation of the dielectric constant, and moreover the dielectric constant is a ground-state property which should be calculable according to the Hohenberg-Kohn Theorems regardless of whether electronic excitation energies are computed correctly. We should look instead for the source of any discrepancies in other deficiencies of LDA.
We can also look at the variation of \(\epsilon\) with frequency by checking em_resp/freq_0.1000/epsilon, which is the result for frequency \(\hbar \omega\) = 0.1 Ha = 2.72 eV.
!cat "7-Sternheimer/em_resp/freq_0.1000/epsilon"
# Real part of dielectric constant
29.065849 0.000000 0.000000
0.000000 29.065849 0.000000
0.000000 0.000000 29.065849
Isotropic average 29.065849
# Imaginary part of dielectric constant
15.274348 0.000000 0.000000
0.000000 15.274348 -0.000000
-0.000000 -0.000000 15.274348
Isotropic average 15.274348
Is symmetry still preserved? How big are the real and imaginary parts?
You will see in the iterations that some lines list a negative number for the state index, i.e. “-2”: this means the negative frequency for state 2, as we need both \(\omega\) and \(-\omega\) for the dielectric constant calculations. As an exercise, you can calculate values from 0 to 0.2 Ha, with 10 points in between (using the EMFreqs block), and make plots of the real and imaginary parts. Compare to results from time-propagation if you have done that tutorial – the results are identical in principle, but can differ due to the use of \(i \eta\) here, and the fact that that time-propagation is not linear response and depends somewhat on the strength of the perturbation used.
Two final notes: We can study partially periodic systems (e.g. monolayer hexagonal boron nitride) through this same approach, where we use the finite scheme for finite directions, and this \(k \cdot p\) scheme for periodic directions. We can also extend this scheme to the calculation of magnetic susceptibilities, and even to nonlinear properties like magneto-optical susceptibilities \(^6\).
References¶
X. Andrade, S. Botti, M. A. L. Marques, and A. Rubio, “Time-dependent density functional theory scheme for efficient calculations of dynamic (hyper)polarizabilities,” J. Chem. Phys. 126, 184106 (2007)
S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Gianozzi, “Phonons and related crystal properties from density-functional perturbation theory,” Rev. Mod. Phys. 73, 515 (2001).
R. Resta, “Macroscopic polarization in crystalline dielectrics: the geometric phase approach,” Rev. Mod. Phys. 66, 899 (1994)
D. A. Strubbe, L. Lehtovaara, A. Rubio, M. A. L Marques, and S. G. Louie, “Response Functions in TDDFT: Concepts and Implementation,” in Fundamentals of Time-Dependent Density Functional Theory, Lecture Notes in Physics, Springer (2012).
X. Andrade, D. A. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodríguez, A. Varas, I. Theophilou, N. Helbig, M. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and Á. Rubio, “Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems,” Phys. Chem. Chem. Phys. 17, 31371-31396 (2015).
I. V. Lebedeva, D. A. Strubbe, I. V. Tokatly, and A. Rubio, “Orbital magneto-optical response of periodic insulators from first principles,” npj Comput. Mater. 5, 32 (2019).
Tutorial Validation Checks¶
from postopus import Run
import numpy as np
velocity = np.loadtxt("7-Sternheimer/kdotp/velocity")
np.testing.assert_allclose(velocity[0], [1. , -0.28036, -0.06523, -0.06523, -0.06523], rtol = 0.01)
np.testing.assert_allclose(velocity[20], [ 1. , -0.26041, -0.19283, -0.06232, -0.06232], rtol = 0.01)