Numerical gradients

The form of vibronic coupling is essentially the derivative of the wave function. Each component of the vibronic coupling vector can be calculated with numerical differentiation methods using wave functions at displaced geometries. This is the procedure used in MOLPRO.^[8]

First order accuracy can be achieved with forward difference formula:

${\displaystyle (\mathbf {f} _{k'k})_{l}\approx {\frac {1}{d}}\left[\gamma ^{k'k}(\mathbf {R} |\mathbf {R} +d\mathbf {e} _{l})-\gamma ^{k'k}(\mathbf {R} |\mathbf {R} )\right]}$

Second order accuracy can be achieved with central difference formula:

${\displaystyle (\mathbf {f} _{k'k})_{l}\approx {\frac {1}{2d}}\left[\gamma ^{k'k}(\mathbf {R} |\mathbf {R} +d\mathbf {e} _{l})-\gamma ^{k'k}(\mathbf {R} |\mathbf {R} -d\mathbf {e} _{l})\right]}$

Here, ${\displaystyle \mathbf {e} _{l}}$ is a unit vector along direction ${\displaystyle l}$. ${\displaystyle \gamma ^{k'k}}$ is the transition density between the two electronic states.

${\displaystyle \gamma ^{k'k}(\mathbf {R} _{1}|\mathbf {R} _{2})=\langle \chi _{k'}(\mathbf {r}$ ;\mathbf {R} _{1})\,|\,\chi _{k}(\mathbf {r} ;\mathbf {R} _{2})\rangle _{(\mathbf {r} )}}

Evaluation of electronic wave functions for both electronic states are required at N displacement geometries for first order accuracy and 2*N displacements to achieve second order accuracy, where N is the number of nuclear degrees of freedom. This can be extremely computationally demanding for large molecules.

As with other numerical differentiation methods, the evaluation of nonadiabatic coupling vector with this method is numerically unstable, limiting the accuracy of the result. Moreover, the calculation of the two transition densities in the numerator are not straightforward. The wave functions of both electronic states are expanded with Slater determinants or configuration state functions (CSF). The contribution from the change of CSF basis is too demanding to evaluate using numerical method, and is usually ignored by employing an approximate diabatic CSF basis. This will also cause further inaccuracy of the calculated coupling vector, although this error is usually tolerable.

Analytic gradient methods

Evaluating derivative couplings with analytic gradient methods has the advantage of high accuracy and very low cost, usually much cheaper than one single point calculation. This means an acceleration factor of 2N. However, the process involves intense mathematical treatment and programming. As a result, few programs have currently implemented analytic evaluation of vibronic couplings at wave function theory levels. Details about this method can be found in ref.^[9] For the implementation for SA-MCSCF and MRCI in COLUMBUS, please see ref.^[10]

TDDFT-based methods

The computational cost of evaluating the vibronic coupling using (multireference) wave function theory has led to the idea of evaluating them at the TDDFT level, which indirectly describes the excited states of a system without describing its excited state wave functions. However, the derivation of the TDDFT vibronic coupling theory is not trivial, since there are no electronic wave functions in TDDFT that are available for plugging into the defining equation of the vibronic coupling.^[5]

In 2000, Chernyak and Mukamel^[11] showed that in the complete basis set (CBS) limit, knowledge of the reduced transition density matrix between a pair of states (both at the unperturbed geometry) suffices to determine the vibronic couplings between them. The vibronic couplings between two electronic states are given by contracting their reduced transition density matrix with the geometric derivatives of the nuclear attraction operator, followed by dividing by the energy difference of the two electronic states:

${\displaystyle (\mathbf {f} _{k'k})_{l}={\frac {1}{E_{k}-E_{k'}}}\sum _{pq}\langle \psi _{p}|{\frac {\partial }{\partial \mathbf {e} _{l}}}{\hat {V}}_{\rm {ne}}|\psi _{q}\rangle (\gamma ^{k'k}(\mathbf {R} |\mathbf {R} ))_{pq}}$

This enables one to calculate the vibronic couplings at the TDDFT level, since although TDDFT does not give excited state wave functions, it does give reduced transition density matrices, not only between the ground state and an excited state, but also between two excited states. The proof of the Chernyak-Mukamel formula is straightforward and involves the Hellmann-Feynman theorem. While the formula provides useful accuracy for a plane-wave basis (see e.g. ref.^[12]), it converges extremely slowly with respect to the basis set if an atomic orbital basis set is used, due to the neglect of the Pulay force. Therefore, modern implementations in molecular codes typically use expressions that include the Pulay force contributions, derived from the Lagrangian formalism.^[5][6][7] They are more expensive than the Chernyak-Mukamel formula, but still much cheaper than the vibronic couplings at wave function theory levels (more specifically, they are roughly as expensive as the SCF gradient for ground state-excited state vibronic couplings, and as expensive as the TDDFT gradient for excited state-excited state vibronic couplings). Moreover, they are much more accurate than the Chernyak-Mukamel formula for realistically sized atomic orbital basis sets.^[5]