First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance

First-order nonadiabatic coupling matrix elements (NACMEs) are key for phenomena such as nonradiative transitions and excited-state decay, yet a consistent and practical first principles treatment has been elusive for molecules with more than a few heavy atoms. Here we present theory, implementation using Gaussian basis sets, and benchmarks of first-order NACMEs between ground and excited states in the framework of time-dependent hybrid density functional theory (TDDFT). A time-dependent response approach to NACMEs which avoids explicit computation of excited-state wave functions is outlined. In contrast to previous approaches, the present treatment produces exact analytical derivative couplings between time-dependent Kohn-Sham (TDKS) determinants in a finite atom-centered basis set. As in analytical gradient theory, derivative molecular orbital coefficients can be eliminated, making the computational cost independent of the number of nuclear degrees of freedom. Our expression reduces to the exact Chernyak-Mukamel formula for first-order NACMEs in the complete basis-set limit, but greatly improves basis-set convergence in finite atom-centered basis sets due to additional Pulay type terms. The Chernyak-Mukamel formula is shown to be equivalent to the Hellmann-Feynman contribution in analytical gradient theory. Our formalism may be implemented in TDDFT analytical excited-state gradient codes with minor modifications. Tests for systems with up to 147 atoms show that evaluation of first-order NACMEs causes total computation times to increase by an insignificant 10% on average. The resolution-of-the-identity approximation for the Coulomb energy (RI-J) reduces the computational cost by an order of magnitude for nonhybrid functionals, while errors are insignificant with standard auxiliary basis sets. We compare the computed NACMEs to full configuration interaction (FCI) in benchmark results for diatomic molecules; hybrid TDDFT and FCI are found to be in agreement for regions of the potential energy curve where the Kohn-Sham ground-state reference is stable and the character of the excitation is properly captured by the present functionals. With these developments, nonadiabatic molecular dynamics simulations of molecular systems in the 100 atoms regime are within reach.