Sternheimer linear response

The Sternheimer approach to perturbation theory allows efficient calculations of linear and non-linear response properties. ¹

The basis of this method, just as in standard perturbation theory, is to calculate the variation of the wave-functions $\psi^{1}$ under a given perturbing potential. The advantage of the method is that the variations are obtained by solving the linear equation

$$ (H^0-\epsilon^0 + \omega )|\psi^{1}>=-P_{\rm c} H^{1}|\psi^{0}>\ , $$

that only depends on the occupied states instead of requiring an (infinite) sum over unoccupied states. In the case of (time-dependent) density functional theory the variation of the Hamiltonian includes a term that depends on the variation of the density, so this equation must be solved self-consistently.

To run a Sternheimer calculation with Octopus, the only previous calculation you need is a ground-state calculation. For this tutorial we will use a water molecule, with this basic input file for the ground state:

CalculationMode = gs

%Coordinates
 'O'  |  0.000000  | -0.553586  |  0.000000
 'H'  |  1.429937  |  0.553586  |  0.000000
 'H'  | -1.429937  |  0.553586  |  0.000000
%

Radius = 10
Spacing = 0.435
ConvRelDens = 1e-6

 CalculationMode = gs

 %Coordinates
  'O'  |  0.000000  | -0.553586  |  0.000000
  'H'  |  1.429937  |  0.553586  |  0.000000
  'H'  | -1.429937  |  0.553586  |  0.000000
 %

 Radius = 10
 Spacing = 0.435
 ConvRelDens = 1e-6

We use a tighter setting on SCF convergence (ConvRelDens) which will help the ability of the Sternheimer calculation to converge numerically, and we increase a bit the size of the box as response calculations tend to require more space around the molecule than ground-state calculations to be converged. ²

After the ground-state calculation is finished, we change the run mode to em_resp, to run a calculation of the electric-dipole response:

 CalculationMode = em_resp

Next, to specify the frequency of the response we use the EMFreqs block; in this case we will use three values 0.00, 0.15 and 0.30 {{< Hartree > }}:

 %EMFreqs
 3 | 0.0 | 0.3
 %

and we will also specify a small imaginary part to the frequency of 0.1 eV , which avoids divergence on resonance:

 EMEta = 0.1*eV

and finally we add a specification of the linear solver, which will greatly speed things up compared to the default:

 LinearSolver = qmr_dotp
 ExperimentalFeatures = yes

In the run, you will see calculations for each frequency for the ‘‘x’’, ‘‘y’’, and ‘‘z’’ directions, showing SCF iterations, each having linear-solver iterations for the individual states’ $\psi^{1}$, labelled by the k-point/spin (ik) and state (ist). The norm of $\psi^{1}$, the number of linear-solver iterations (iter), and the residual $\left|(H^0-\epsilon^0 + \omega )|\psi^{1}>+P_{\rm c} H^{1}|\psi^{0}>\right|$ are shown for each. First we see the static response:

****************** Linear-Response Polarizabilities ******************
Wavefunctions type: Complex
Calculating response for   3 frequencies.
**********************************************************************

Info: EM Resp. restart information will be written to 'restart/em_resp'.
Info: EM Resp. restart information will be read from 'restart/em_resp'.
Info: Finished reading information from 'restart/em_resp'.
Info: Calculating response for the x-direction and frequency 0.0000.
--------------------------------------------
LR SCF Iteration:   1
Frequency:             0.000000 Eta :             0.003675
   ik  ist                norm   iters            residual
    1    1            0.161255      17        0.343555E-03
    1    2            1.700485      17        0.777896E-02
    1    3            1.530852      17        0.852600E-02
    1    4            1.301154      17        0.229271E-02

             Info: Writing states. 2024/10/21 at 12:15:59


        Info: Finished writing states. 2024/10/21 at 12:15:59

SCF Residual:     0.125284E+01 (abs),     0.156605E+00 (rel)

Later will come the dynamical response. The negative state indices listed indicate response for $-\omega$. For each frequency, the code will try to use a saved response density from the closest previously calculated frequency.

Info: Calculating response for the x-direction and frequency 0.1500.
--------------------------------------------
LR SCF Iteration:   1
Frequency:             0.150000 Eta :             0.003675
   ik  ist                norm   iters            residual
    1    1            0.151339       7        0.111851E-02
    1   -1            0.197593      17        0.828131E-04
    1    2            1.006174       7        0.534947E-02
    1   -2            1.840154      17        0.234366E-02
    1    3            0.808163       7        0.706163E-02
    1   -3            1.593366      17        0.938379E-02
    1    4            0.792556       7        0.516886E-02
    1   -4            1.612798      17        0.822362E-03

             Info: Writing states. 2024/10/21 at 12:16:02


        Info: Finished writing states. 2024/10/21 at 12:16:02


             Info: Writing states. 2024/10/21 at 12:16:02


        Info: Finished writing states. 2024/10/21 at 12:16:02

SCF Residual:     0.763103E-01 (abs),     0.953879E-02 (rel)
--------------------------------------------

At the end, you will have a directory called em_resp containing a subdirectory for each frequency calculated, each in turn containing eta (listing $\eta$ = 0.1 eV ), alpha (containing the real part of the polarizability tensor), and cross_section (containing the cross-section for absorption, based on the imaginary part of the polarizability).

For example, em_resp/freq_0.0000/alpha says

# Polarizability tensor [b^3]
           10.713755           -0.000000           -0.000000
           -0.000000           10.734819           -0.000000
            0.000000           -0.000000           10.874435
Isotropic average           10.774336

Exercise: compare results for polarizability or cross-section to a calculation from time-propagation or the Casida approach.