General

Ref: [95].

Perturbation theory describes the reaction of a system onto an external perturbation. At the time when the perturbation occurs, the system is in its ground state (or unperturbed state). The perturbation then changes slightly the potential energy surface and therefore also the state where the system's energy is minimum. As a consequence, the system tries to move towards that state of minimum energy. This movement or the new state often have properties which can be accessed experimentally. Example: An external electric field will slightly deform the electronic cloud, creating a dipole. That dipole can then be measured.

Assume that the magnitude of the perturbation is small compared to the strength of the forces acting in the unperturbed system. Then, the change in the minimum energy state will be small as well and perturbation theory can be applied to compute how the system reacts onto the perturbation. Generally, the Schrödinger equation is expanded in powers of the perturbation parameter (ex: the strength of the electric field), and the equations obtained for those powers are solved individually. At power zero, one refinds the equation of the unperturbed system:

$$(H^{(0)}- \varepsilon_k) \vert\varphi_{k}^{(0)}\rangle = 0.$$ (3)

For the part which is linear in the perturbation, the general format of the resulting equation is

$$(H^{(0)}- \varepsilon_k) \vert\varphi_{k}^{(1)}\rangle = - H^{(1)}\vert\varphi_{k}^{(0)}\rangle .$$ (4)

Grosso modo, this equation is solved during a linear response calculation through a wavefunction optimization process for $\vert\varphi_{k}^{(1)}\rangle$.

The presence of a first order perturbation correction for the wavefunctions, $\vert\varphi_{k}^{(\text{\sf tot})}\rangle = \vert\varphi_{k}^{(0)}\rangle + \vert\varphi_{k}^{(1)}\rangle$ implies that the total density of the perturbed system is no longer equal to the unperturbed one, $n^{(0)}$, but also contains a first order perturbation correction, $n^{(1)}$. That density is given by

$$n^{(1)}(\mathbf{r}) = \sum_k \langle \varphi^{(1)}_k \vert \mathbf{r}\rangle \; \langle \mathbf{r}\vert \varphi^{(0)}_k \rangle +$$ (5)

The Hamiltonian depends on the electronic density. Therefore, the first order density correction implies automatically an additional indirect perturbation hamiltonian coming from the expansion of the unperturbed Hamiltonian in the density. It has to be added to the explicit perturbation Hamiltonian determined by the type of the (external) perturbation. The contribution is

$$H^{(1)}_{\text{\sf indirect}}(\mathbf{r}) = \int d^3r^\prime\; \frac{\partial H^{(0)}(\mathbf{r})}{\partial n^{(0)}(\mathbf{r}^\prime)}\; n^{(1)}(\mathbf{r}^\prime)$$ (6)

The calculation of this indirect Hamiltonian represents almost 50% of the computational cost of the response calculation, especially in connection with xc-functionals. After several unsuccessful trials with analytic expressions for the derivative of the xc-potential with respect to the density, this is done numerically. That means that at each step, the xc-potential is calculated for the density $n^{(0)}+\epsilon n^{(1)}$ and for $n^{(0)}-\epsilon n^{(1)}$ (with an $\epsilon$ empirically set to 0.005), and the dervative needed in (6) is calculated as

$$\int d^3r^\prime\; n^{(1)}(\mathbf{r}^\prime)\; \frac{\partial v_{\text{xc}}}{\partial n^{(0)}(\mathbf{r}^\prime)} = \frac{v_{\text{xc}}[n^{(0)}+\epsilon n^{(1)}] - v_{\text{xc}}[n^{(0)}-\epsilon n^{(1)}]} {2\epsilon}.$$ (7)

In the case of the local density approximation, the derivative can be done analytically, in which case it only needs to be done once. This improves the performance of the optimization.