Approach for Higher Order Response Properties

Question

There are three (at least that I'm aware of) commonly used approaches to obtain linear response properties (e.g electric polarizability, optical rotation, NMR shielding tensors)

• Sum over states: The properties can formally be written as a sum of matrix elements of the perturbations $A$ and $B$ over all excited states. In practice, compute enough excited states to converge the property. Tends to converge slowly with the number of states.
• Response functions: The properties can also be written in terms of response relations, which leads us instead to compute the perturbed density with respect to $A$ (or $B$) and contract it with $B$ (or $A$) to compute the property.
• Derivatives: These properties are also derivatives of the energy with respect to these perturbations. One can derive analytic formulas or compute numerical derivatives. Analytical formulas are complex and don't (directly) apply to frequency dependent properties. Numerical derivatives require repeated calculations and tuning the step size.

For optical rotation, and I believe most other linear properties, the second approach above has won out as the best way to do the computation in general. But I'm curious if this holds for computing nonlinear properties (e.g. $n^{\text{th}}$ hyperpolarizability, the Kerr Effect) as well or whether the cost/benefit analysis of these methods changes.

Answer 1

Options 2 and 3 appear to be the same: the response is almost always the response of the energy, since the wave function is determined by the energy principle.

To clarify: in many methods (e.g. Hartree-Fock or CC) one computes the derivative of the energy functional with respect to the property (e.g. polarizability or NMR shielding constants); this turns out to lead to (generalized) response densities that you need to solve from the Schrödinger equation, and at the end you get your property by contracting the densities.

There may also be several ways in which to choose your perturbation, e.g. for NMR shielding constants your variables are the external magnetic field and the nuclear shieldings. I forget the details, but the idea is this: you can first perturb with the shieldings, and then contract with the magnetic field response, but this obviously becomes horribly slow for many nuclei. Instead, you can first perturb with the magnetic field, solving the response equations for the 3 components of the field, after which you can get the shieldings just by contracting the fixed response to the external field with the relevant matrices for the individual nuclei; this approach scales to large numbers of atoms.