Optical spectra from sternheimer¶
We have just seen how to calculate optical spectra in the time domain with a finite perturbation, and in a frequency-domain, linear-response matrix formulation with the Casida equation. Now we will try a third approach, which is in the frequency domain and linear response but rather than using a pseudo-eigenvalue equation as in Casida, uses a self-consistent linear equation, the Sternheimer equation. This approach is also known as density-functional perturbation theory. It has superior scaling, is more efficient for dense spectra, and is more applicable to nonlinear response. One disadvantage is that one needs to proceed one frequency point at a time, rather than getting the whole spectrum at once. We will find we can obtain equivalent results with this approach for the optical spectra as for time propagation and Casida, by calculating the polarizability and taking the imaginary part.
import os
import math
!mkdir -p 4_Optical_spectra_from_Sternheimer
Ground state¶
Before doing linear response, we need to obtain the ground state of the system, for which we can use the same input file as for Optical spectra from Casida, but we will use a tighter numerical tolerance, which helps the Sternheimer equation to be solved more rapidly. Unlike for Casida, no unoccupied states are required. If they are present, they won’t be used anyway.
%%writefile 4_Optical_spectra_from_Sternheimer/inp
stdout = 'stdout_gs.txt'
stderr = 'stderr_gs.txt'
CalculationMode = gs
UnitsOutput = eV_angstrom
Eigensolver = chebyshev_filter
ExtraStates = 4
Radius = 6.5*angstrom
Spacing = 0.24*angstrom
CH = 1.097*angstrom
%Coordinates
"C" | 0 | 0 | 0
"H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3)
"H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3)
"H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3)
"H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3)
%
ConvRelDens = 1e-5
Writing 4_Optical_spectra_from_Sternheimer/inp
!cd 4_Optical_spectra_from_Sternheimer && octopus
Linear response¶
Change the CalculationMode to em_resp and add
the lines starting with %EMFreqs to the input file.
%%writefile 4_Optical_spectra_from_Sternheimer/inp
stdout = 'stdout_em_resp.txt'
stderr = 'stderr_em_resp.txt'
CalculationMode = em_resp
UnitsOutput = eV_angstrom
Radius = 6.5*angstrom
Spacing = 0.24*angstrom
CH = 1.097*angstrom
%Coordinates
"C" | 0 | 0 | 0
"H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3)
"H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3)
"H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3)
"H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3)
%
ConvRelDens = 1e-7
%EMFreqs
5 | 0*eV | 8*eV
9 | 10*eV | 12*eV
%
EMEta = 0.1*eV
Preconditioner = no
LinearSolver = qmr_dotp
ExperimentalFeatures = yes
Overwriting 4_Optical_spectra_from_Sternheimer/inp
The frequencies of interest must be specified, and we choose them based on the what we have seen from the Casida spectrum. The block above specifies 5 frequencies spanning the range 0 to 8 eV (below the resonances) and 9 frequencies spanning the range 10 to 12 eV (where there are peaks). If we didn’t know where to look, then looking at a coarse frequency grid and then sampling more points in the region that seems to have a peak (including looking for signs of resonances in the real part of the polarizability) would be a reasonable approach. We must add a small imaginary part (ηη) to the frequency in order to be able to obtain the imaginary part of the response, and to avoid divergence at resonances. The resonances are broadened into Lorentzians with this width. The larger the ηη, the easier the SCF convergence is, but the lower the resolution of the spectrum.
To help in the numerical solution, we turn off the preconditioner (which sometimes causes trouble here), and use a linear solver that is experimental but will give convergence much faster than the default one.
!cd 4_Optical_spectra_from_Sternheimer && mpirun -n 8 octopus
Spectrum¶
With Postopus we can exctract the needed info out of the files in the em_resp directory:
from postopus import Run
import pandas as pd
import math
run = Run("4_Optical_spectra_from_Sternheimer")
# As cross_section and alpha use both 'freq' as index name we can (horizontally) concat them via `.join`
df = run.em_resp.cross_section().join(run.em_resp.alpha())
c = 137.035999679
df['Re alpha'] = df['Isotropic average']
df['Im alpha'] = df['(1/3)*Tr[sigma]'] * c / (4 * math.pi * df['Energy'])
# Set 'Im alpha' per hand for energy = 0
df.loc[0, 'Im alpha'] = 0
Our result is
df[["Energy", "Re alpha", "Im alpha", "(1/3)*Tr[sigma]"]]
| Energy | Re alpha | Im alpha | (1/3)*Tr[sigma] | ||
|---|---|---|---|---|---|
| freq | |||||
| 0.00 | 0 | 0.00 | 2.639931 | 0.000000 | 0.000000 |
| 2.00 | 0 | 2.00 | 2.698269 | 0.000418 | 0.000077 |
| 4.00 | 0 | 4.00 | 2.897126 | 0.001018 | 0.000373 |
| 6.00 | 0 | 6.00 | 3.348183 | 0.002343 | 0.001289 |
| 8.00 | 0 | 8.00 | 4.691392 | 0.009967 | 0.007312 |
| 10.00 | 0 | 10.00 | 3.549525 | 0.043533 | 0.039921 |
| 10.25 | 0 | 10.25 | 4.202297 | 0.119447 | 0.112273 |
| 10.50 | 0 | 10.50 | 4.828190 | 0.050175 | 0.048312 |
| 10.75 | 0 | 10.75 | 6.725430 | 0.071974 | 0.070951 |
| 11.00 | 0 | 11.00 | 10.372729 | 0.250986 | 0.253173 |
| 11.25 | 0 | 11.25 | 4.029440 | 1.057495 | 1.090952 |
| 11.50 | 0 | 11.50 | -0.608458 | 0.212971 | 0.224592 |
| 11.75 | 0 | 11.75 | 4.216867 | 0.269344 | 0.290215 |
| 12.00 | 0 | 12.00 | -1.133445 | 0.221212 | 0.243425 |
See also¶
Tutorial Validation Checks¶
import numpy as np
np.testing.assert_allclose(df["Im alpha"][8.00], 0.009967, rtol=0.01)
np.testing.assert_allclose(df["Im alpha"][11.25], 1.05750, rtol=0.01)