# On the Mass of Atoms in Molecules: Beyond the Born-Oppenheimer Approximation

Arne Scherrer, Federica Agostini, Daniel Sebastiani, E. k. u. Gross, Rodolphe Vuilleumier

► To cite this version:

Arne Scherrer, Federica Agostini, Daniel Sebastiani, E. k. u. Gross, Rodolphe Vuilleumier. On the Mass of Atoms in Molecules: Beyond the Born-Oppenheimer Approximation. Physical Review X, 2017, 7 (3), pp.031035. 10.1103/PhysRevX.7.031035. [hal-01585233](https://hal.sorbonne-universite.fr/hal-01585233)

HAL Id: hal-01585233

https://hal.sorbonne-universite.fr/hal-01585233

Submitted on 11 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Mass of Atoms in Molecules: Beyond the Born-Oppenheimer Approximation

Arne Scherrer,<sup>1,2,3</sup> Federica Agostini,<sup>4,5,*</sup> Daniel Sebastiani,<sup>1</sup> E. K. U. Gross,<sup>4,6</sup> and Rodolphe Vuilleumier<sup>2,3,†</sup>

<sup>1</sup>Martin-Luther-Universität Halle-Wittenberg, von-Danckelmann-Platz 4, D-06120 Halle, Germany

<sup>2</sup>UMR 8640 ENS-CNRS-UPMC, Département de Chimie,
24 rue Lhomond, École Normale Supérieure, 75005 Paris, France

<sup>3</sup>UPMC Université Paris 06, 4, Place Jussieu, 75005 Paris, France

<sup>4</sup>Max-Planck-Institut für Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany

<sup>5</sup>Laboratoire de Chimie Physique, UMR 8000 CNRS/University Paris-Sud,
University Paris-Saclay, 91405 Orsay, France

<sup>6</sup>Fritz Haber Center for Molecular Dynamics, Institute of Chemistry,
The Hebrew University of Jerusalem, Jerusalem 91904, Israel

(Received 23 May 2017; published 25 August 2017)

Describing the dynamics of nuclei in molecules requires a potential energy surface, which is traditionally provided by the Born-Oppenheimer or adiabatic approximation. However, we also need to assign masses to the nuclei. There, the Born-Oppenheimer picture does not account for the inertia of the electrons, and only bare nuclear masses are considered. Nowadays, experimental accuracy challenges the theoretical predictions of rotational and vibrational spectra and requires the participation of electrons in the internal motion of the molecule. More than 80 years after the original work of Born and Oppenheimer, this issue has still not been solved, in general. Here, we present a theoretical and numerical framework to address this problem in a general and rigorous way. Starting from the exact factorization of the electron-nuclear wave function, we include electronic effects beyond the Born-Oppenheimer regime in a perturbative way via position-dependent corrections to the bare nuclear masses. This maintains an adiabaticlike point of view: The nuclear degrees of freedom feel the presence of the electrons via a single potential energy surface, whereas the inertia of electrons is accounted for and the total mass of the system is recovered. This constitutes a general framework for describing the mass acquired by slow degrees of freedom due to the inertia of light, bounded particles; thus, it is applicable not only in electron-nuclear systems but in light-heavy nuclei or ions as well. We illustrate this idea with a model of proton transfer, where the light particle is the proton and the heavy particles are the oxygen atoms to which the proton is bounded. Inclusion of the light-particle inertia allows us to gain orders of magnitude in accuracy. The electron-nuclear perspective is adopted, instead, to calculate position-dependent mass corrections using density functional theory for a few polyatomic molecules at their equilibrium geometry. These data can serve as input for the computation of high-precision molecular spectra.

DOI: 10.1103/PhysRevX.7.031035

Subject Areas: Atomic and Molecular Physics,
Chemical Physics,
Computational Physics

## I. INTRODUCTION

The Born-Oppenheimer (BO) [1], or adiabatic, treatment of the coupled motion of electrons and nuclei in molecular systems is among the most fundamental approximations in condensed matter physics and chemical physics. Based on the hypothesis that part of the system, usually electrons or protons, evolves on a much shorter time scale than the rest, i.e., (heavy) nuclei, the BO approximation allows one to visualize molecules as a set of nuclei moving on a single

potential energy surface that represents the effect of the electrons in a given eigenstate. Yet, it is an approximation, yielding the correct dynamics only in the limit of infinite nuclear masses. For instance, when compared to highly accurate molecular spectroscopy measurements, theoretical predictions might deviate from experimentally observed behavior.

Often in the literature, in relation to calculations of rotational and vibrational spectra of light molecules, for instance, hydrogen-based [2–9], the question of which masses [10–13] are to be considered is addressed, aiming to provide accurate numerical predictions. However, beyond the simple interest in spectroscopy applications, this question carries a very fundamental significance, which can be summarized as follows: How is the inertia of the electrons accounted for in the nuclear motion, when the BO approximation is employed? In fact, within the BO approximation, the electrons appear only implicitly in the dynamics, as a

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*federica.agostini@u-psud.fr</footnote>

<footnote>†rodolphe.vuilleumier@ens.fr</footnote>
potential energy contribution to the Hamiltonian driving the motion of the nuclei. The kinetic energy arising from the molecular motion then involves only the bare nuclear masses. Indeed, ideas to improve the agreement between numerical and experimental spectroscopy results have been proposed, such as the use of atomic masses [14] or fractional masses [15]. Such propositions do not unambiguously answer the above question and, therefore, might still be of little interest to a broader community.

In addition, performing a full nonadiabatic treatment of the coupled electron-nuclear problem, with a numerical cost that is much larger than a BO calculation, does not help from a fundamental point of view. The following question is still unanswered: What is the mechanism by which the inertia of the electrons affects the mass of the heavy degrees of freedom beyond the dependence of the potential energy surface on the electronic kinetic energy? An alternative approach, pioneered by Bunker and Moss [2,3,16], is to treat nonadiabatic effects perturbatively, but applications are still limited to di- and triatomic molecules. In connection with the perturbation idea of Bunker and Moss, accurate numerical calculations have been performed on small molecules, like H₂, D₂, HD, H₃⁺ [17–21]. However, despite the effort to push forth the applications, it seems that the basic formalism still represents a major obstacle for the treatment of molecular systems comprising more than three atoms. The main reason for this is found in the use of internal coordinates, obtained after separation of the rotational and translational degrees of freedom of the center of mass of the molecule, as a starting point for the application of the perturbation approach. Several complications arise when using internal coordinates: (i) Internal coordinates are not uniquely defined; thus, calculations of nuclear displacements, which are crucial to determine nonadiabatic corrections, have to be reformulated for each choice of coordinates, either analytically or numerically. (ii) As the size of the molecular system of interest increases, derivative operators in internal coordinates become more and more involved, and their implementation becomes less and less desirable [22]. (iii) Internal coordinates are defined in the nuclear center of mass [4], which depends on the mass of the nuclei; thus, electronic structure calculations to determine the nonadiabatic corrections have to be repeated when considering different isotopes. (iv) Derivatives with respect to nuclear displacements in internal coordinates cannot take advantage of analytic methods, as is the case of derivatives in Cartesian coordinates [23].

In the present paper, we examine this problem in the framework of the exact factorization of the electron-nuclear wave function [24]. This (nonadiabatic) reformulation of the quantum-mechanical problem is used as a starting point to develop a procedure that settles the issue described above in a rigorous way. The key point in the exact factorization is that the electronic effect on the nuclear system is taken into account by time-dependent vector and scalar potentials. These concepts are the generalization of similar, but static,

quantities appearing also within the BO approximation. We show that nonadiabatic effects can be accounted for, by formulating a theory that treats these effects as a perturbation to the BO problem. Such a framework has been discussed in previous work [23] to derive the nuclear velocity perturbation theory [25] for vibrational circular dichroism [26]. As we will show below, here we propose a new perspective on the nuclear velocity perturbation theory, which will allow us to access a broader class of both static, e.g., energetics, and dynamical, e.g., vibrational spectra, problems in quantum mechanics. Within nuclear velocity perturbation theory, nonadiabatic effects can be included by taking into account corrections to the BO approximation up to within linear order in the classical nuclear velocity. We show here that this is equivalent to a perturbation approach where the small parameter is the electron-nuclear mass ratio [27].

The major achievement of such a formulation is presented in this paper: Electronic nonadiabatic effects appear as a position-dependent mass correction to the bare nuclear mass, up to within linear order in the perturbation. From a fundamental perspective, we prove that it is possible to recover an adiabaticlike structure of the Hamiltonian governing the dynamics of the heavy degrees of freedom, with a kinetic energy contribution and a separate potential energy term. Since the mass correction can be fully identified with the electronic mass, totally missing in the BO approximation, we propose a theory able to restore a fundamental property, often overlooked, of the dynamical problem: the translational invariance of an isolated system with its physical mass, i.e., nuclear *and* electronic. If in the BO approximation the nuclear masses are made position dependent in the way proposed in this paper, the center of mass can be separated from rotations and internal vibrations, and it evolves as a free particle with mass equal to the total mass of the system (expected from the Galilean invariance of the problem [13]). This property enables us to apply the perturbation approach *before* moving to the molecular center-of-mass reference frame, with the formal advantage of a very simple and intuitive theory. From an algorithmic perspective, the corrections to the mass involve only ground-state properties and can be calculated as a response to the nuclear motion, within standard perturbation theory [28–30]. Therefore, we are able to perform numerical studies of molecular systems, easily pushing the applications beyond di- and triatomic molecules. The experimental implications are clear: The approach proposed here has the potential to predict and to describe rovibrational spectroscopic data for a large class of molecular systems when high accuracy is required. The theory developed here will make the nuclear velocity perturbation scheme of Refs. [23,31] consistent. In previous work, we have focused only on the effect of nuclear motion on electronic properties, neglecting the “backreaction” of electrons. The advantage of taking this effect into account as done in the present work is twofold: (i) From a fundamental perspective, even a weak coupling to electronic excited states has an observable effect on nuclear properties, and (ii) from a
practical perspective, the translational motion of the center of mass of the *total* system can now be correctly separated from the internal rotations and vibrations.

The paper is organized as follows. First, we show how, starting from the exact factorization, nonadiabatic effects are included by constructing a perturbative scheme based on the BO approach. Then, we prove that the vector potential of the theory can be expressed as a position-dependent correction to the bare nuclear mass. In the nuclear Hamiltonian, nonadiabatic effects are taken into account in an adiabaticlike picture, if the nuclear masses are corrected for the electronic contribution. We prove that (i) the position-dependent corrections sum up to the total electronic mass of the complete system, and (ii) the Hamiltonian with position-dependent dressed masses is appropriate to compute rotational and vibrational spectra, as it is possible to exactly separate the center-of-mass motion. Results are first presented for a model of a hydrogen bond where the light particle is the shared proton. Then, we present the computation of the position-dependent mass corrections using density functional theory (DFT) for some molecular systems, up to six atoms, at their equilibrium geometry. To illustrate the effect on vibrational spectra, the order of magnitude of this nonadiabatic correction is evaluated in the harmonic approximation.

## II. BEYOND THE BORN-OPPENHEIMER APPROXIMATION

### A. Exact factorization of the electron-nuclear wave function

The exact factorization of the electron-nuclear wave function has been presented [24] and discussed [32,33] in previous work. Therefore, we only introduce here the basic formalism, and we refer to the above references for a detailed presentation.

A system of interacting particles, which will be taken as electrons of mass $m_e$ and nuclei of masses $M_ν$, is described by the Hamiltonian $\hat{H} = \hat{T}_n + \hat{H}_{BO}$, with $\hat{T}_n$ the nuclear kinetic energy and $\hat{H}_{BO}$ the standard BO Hamiltonian. The evolution of the electron-nuclear wave function $\Psi(\mathbf{r}, \mathbf{R}, t)$, in the absence of an external time-dependent field, is described by the time-dependent Schrödinger equation $\hat{H}\Psi = i\hbar\partial_t\Psi$. The symbols $\mathbf{r}, \mathbf{R}$ collectively indicate the Cartesian coordinates of $N_{el}$ electrons and $N_n$ nuclei, respectively, in a fixed laboratory frame. When the exact factorization is employed, the solution of the time-dependent Schrödinger equation is written as the product $\Psi(\mathbf{r}, \mathbf{R}, t) = \Phi_R(\mathbf{r}, t)\chi(\mathbf{R}, t)$, where $\chi(\mathbf{R}, t)$ is the nuclear wave function and $\Phi_R(\mathbf{r}, t)$ is an electronic factor that parametrically depends on the nuclear configuration $\mathbf{R}$. Here, $\Phi_R(\mathbf{r}, t)$ satisfies the partial normalization condition $\int d\mathbf{r}|\Phi_R(\mathbf{r}, t)|^2 = 1 \forall \mathbf{R}, t$, which makes the factorization unique up to a gauge transformation. Starting from the time-dependent Schrödinger equation for $\Psi(\mathbf{r}, \mathbf{R}, t)$,

Frenkel's action principle [34–36] and the partial normalization condition yield the evolution equations for $\Phi_R(\mathbf{r}, t)$ and $\chi(\mathbf{R}, t)$,

$$[\hat{H}_{el} - \epsilon(\mathbf{R}, t)]\Phi_R = i\hbar\partial_t\Phi_R \quad \text{and} \quad \hat{H}_n\chi = i\hbar\partial_t\chi. \quad (1)$$

Here, the electronic and nuclear Hamiltonians are $\hat{H}_{el} = \hat{H}_{BO} + \hat{U}_{en}[\Phi_R, \chi]$ and $\hat{H}_n = \sum_ν[-i\hbar\nabla_ν + \mathbf{A}_ν(\mathbf{R}, t)]^2/(2M_ν) + \epsilon(\mathbf{R}, t)$, respectively. The index $ν$ is used to label the nuclei. The electron-nuclear coupling operator (ENCO),

$$\hat{U}_{en}[\Phi_R, \chi] = \sum_{\nu=1}^{N_n} \frac{1}{M_\nu} \left[ \frac{[-i\hbar\nabla_\nu - \mathbf{A}_\nu(\mathbf{R}, t)]^2}{2} + \left( \frac{-i\hbar\nabla_\nu\chi(\mathbf{R}, t)}{\chi(\mathbf{R}, t)} + \mathbf{A}_\nu(\mathbf{R}, t) \right) \times (-i\hbar\nabla_\nu - \mathbf{A}_\nu(\mathbf{R}, t)) \right], \quad (2)$$

the time-dependent vector potential (TDVP),

$$\mathbf{A}_\nu(\mathbf{R}, t) = \langle \Phi_R(t) | -i\hbar\nabla_\nu \Phi_R(t) \rangle_\mathbf{r}, \quad (3)$$

and the time-dependent potential energy surface (TDPES),

$$\epsilon(\mathbf{R}, t) = \langle \Phi_R(t) | \hat{H}_{BO} + \hat{U}_{en} - i\hbar\partial_t | \Phi_R(t) \rangle_\mathbf{r}, \quad (4)$$

mediate the *exact* coupling between the two subsystems; thus, they include all effects beyond BO. The symbol $\langle \dots \rangle_\mathbf{r}$ indicates integration over the electronic coordinates. The TDVP and TDPES transform [24] as standard gauge potentials when the electronic and nuclear wave functions transform with a phase $\theta(\mathbf{R}, t)$. The gauge, the only freedom in the definition of the electronic and nuclear wave functions, will be fixed below.

### B. Large nuclear mass limit

Starting from the exact factorization described above, we now consider the limit of large nuclear masses. The ENCO is inversely proportional to the nuclear masses $M_ν$; then, the BO limit [37] corresponds to the solution of Eq. (1), setting the ENCO to zero [23,27]. Formally, however, approaching this limit of large but finite nuclear masses depends on the physical situation considered [38]. In the time-dependent case, keeping the kinetic energy fixed, it has been shown [38] that the BO limit is recovered asymptotically in terms of a small expansion parameter $\mu^4$ used to scale the nuclear mass, $M \rightarrow M^{(\mu)} \equiv M/\mu^4$. Making $\mu$ approach zero corresponds to the ratio of the nuclear mass over the electron mass $M^{(\mu)}/m_e$ going to infinity. This scaling factor will be used only to perturbatively estimate the order of the terms in the electronic equation and will be set equal to unity to recover the values of the physical masses. Since the nuclear mass is made larger, the nuclear dynamics is slower such that the time variable must then be scaled as well, by a
factor $\mu^2$, i.e., $t \to t/\mu^2$ [38], increasing the separation of time scales between the light and heavy particles. Similarly, following a simple scaling argument, the nuclear momentum behaves as $\mu^{-2}$ in the semiclassical limit (see Appendix A). Then, the ENCO from Eq. (2) scales with $\mu^4$ as

$$ \hat{U}_{\text{en},\mu}[\Phi_{\mathbf{R}},\chi] = \sum_{\nu=1}^{N_n} \left[ \frac{\mu^4 [-i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R},t)]^2}{2} + \frac{\mu^2}{M_{\nu}} (\lambda_{\nu}(\mathbf{R},t) + \mu^2 \mathbf{A}_{\nu}(\mathbf{R},t)) \times (-i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R},t)) \right], \quad (5) $$

where $\lambda_{\nu}(\mathbf{R},t) = \mu^2\{[-i\hbar\nabla_{\nu}\chi(\mathbf{R},t)]/\chi(\mathbf{R},t)\}$. Note that $\lambda_{\nu}(\mathbf{R},t)$ tends towards a quantity independent of $\mu$ in the limit of small $\mu$ since $-i\hbar\nabla_{\nu}\chi/\chi$ is related to the nuclear momentum [23,33] and thus scales as $\mu^{-2}$.

Using the definition in Eq. (4), the scaled TDPES is

$$ \epsilon_{\mu}(\mathbf{R},t) = \langle \Phi_{\mathbf{R}}(t) | \hat{H}_{\text{BO}} | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}} + \mu^2 \langle \Phi_{\mathbf{R}}(t) | -i\hbar \partial_t | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}} + \mu^4 \sum_{\nu=1}^{N_n} \frac{1}{2M_{\nu}} \langle \Phi_{\mathbf{R}}(t) | [-i\hbar \nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R},t)]^2 | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}}, \quad (6) $$

noting that the second term in Eq. (5) does not contribute to the TDPES. This statement can be easily verified in Eq. (2), by computing the expectation value of the second term on the right-hand side with the electronic wave function $\Phi_{\mathbf{R}}(\mathbf{r},t)$ and by recalling the definition of the TDVP given in Eq. (3): The contribution arising from $-i\hbar\nabla_{\nu}$ exactly cancels the contribution from $\mathbf{A}_{\nu}(\mathbf{R},t)$.

### C. Perturbative expansion

The scaled ENCO, from Eq. (5), and the scaled TDPES, from Eq. (6), are used to rewrite the electronic equation (1) such that its dependence on $\mu^4$ is explicit. To this end, the time derivative of the right-hand side of Eq. (1) also has to be scaled with $\mu^2$, as indicated above. We thus obtain

$$ [\hat{H}_{\text{BO}} + \hat{U}_{\text{en},\mu}[\Phi_{\mathbf{R}},\chi] - \epsilon_{\mu}(\mathbf{R},t)]\Phi_{\mathbf{R}} = i\hbar\mu^2\partial_t\Phi_{\mathbf{R}}, \quad (7) $$

which can be solved perturbatively in powers of $\mu^2$, with its solution of the form $\Phi_{\mathbf{R}}(\mathbf{r},t) = \Phi_{\mathbf{R}}^{(0)}(\mathbf{r},t) + \mu^2\Phi_{\mathbf{R}}^{(1)}(\mathbf{r},t) + \cdots [1,39]$.

The time dependence appears only at order $\mu^2$, as it is clear from Eqs. (6) and (7). Therefore, the time dependence of $\Phi_{\mathbf{R}}^{(0)}(\mathbf{r},t) = \phi_{\mathbf{R}}^{(0)}(\mathbf{r})$ can be dropped, and $\phi_{\mathbf{R}}^{(0)}(\mathbf{r})$ satisfies the zeroth-order equation

$$ [\hat{H}_{\text{BO}} - \epsilon^{(0)}(\mathbf{R})]\phi_{\mathbf{R}}^{(0)} = 0, \quad (8) $$

with $\epsilon^{(0)}(\mathbf{R})$ the first term on the right-hand side of Eq. (6). Here, $\phi_{\mathbf{R}}^{(0)}(\mathbf{r})$ is an eigenstate of the BO Hamiltonian with eigenvalue $\epsilon^{(0)}(\mathbf{R}) = \epsilon_{\text{BO}}^{(0)}(\mathbf{R})$, chosen to be the ground state.

At the zeroth order, (i) the TDVP identically vanishes, $\mathbf{A}_{\nu}^{(0)}(\mathbf{R},t) = 0$, as in the absence of a magnetic field $\phi_{\mathbf{R}}^{(0)}(\mathbf{r})$ can be real; (ii) the evolution of the nuclear wave function is determined by the usual BO equation; (iii) in order to uniquely determine the nuclear, $\chi(\mathbf{R},t)$, and electronic, $\Phi_{\mathbf{R}}(\mathbf{r},t)$, wave functions via Eq. (1), one has to make a choice of gauge, and for the derivations presented here, we impose $\langle\phi_{\mathbf{R}}^{(0)}|\Phi_{\mathbf{R}}(t)\rangle \in \mathbb{R}$.

The electronic equation at the next order yields

$$ [\hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R})]\Phi_{\mathbf{R}}^{(1)} = i \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R},t) \cdot (\hbar\nabla_{\nu}\phi_{\mathbf{R}}^{(0)}), \quad (9) $$

where $\lambda'_{\nu}(\mathbf{R},t) = [\lambda_{\nu}(\mathbf{R},t) + \mu^2\mathbf{A}_{\nu}(\mathbf{R},t)]/M_{\nu}$ from Eq. (2). We neglected the TDVP from the term in parentheses since $\mathbf{A}_{\nu}(\mathbf{R},t)$ is $O(\mu^2)$. However, in $\lambda'_{\nu}$, we include a term $O(\mu^2)$, which will be analyzed below along with the TDVP. Appendix B presents the connection between Eq. (9) and nuclear velocity perturbation theory, thus providing a numerical scheme [23] to compute $\Phi_{\mathbf{R}}^{(1)}(\mathbf{r},t)$ within perturbation theory [25]. To this end, we write $N_n$ Sternheimer equations

$$ [\hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R})]\phi_{\mathbf{R},\nu}^{(1)}(\mathbf{r}) = \hbar\nabla_{\nu}\phi_{\mathbf{R}}^{(0)}(\mathbf{r}), \quad (10) $$

which provide an expression for $\Phi_{\mathbf{R}}^{(1)}(\mathbf{r})$ decomposed in $N_n$ contributions, labeled by the nuclear index $\nu$, i.e., $\phi_{\mathbf{R},\nu}^{(1)}(\mathbf{r})$, independent of $\lambda'_{\nu}(\mathbf{R},t)$. Therefore, the electronic wave function up to within $O(\mu^2)$ is

$$ \Phi_{\mathbf{R}}(\mathbf{r},t) = \phi_{\mathbf{R}}^{(0)}(\mathbf{r}) + \mu^2 i \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R},t) \cdot \phi_{\mathbf{R},\nu}^{(1)}(\mathbf{r}), \quad (11) $$

where $\phi_{\mathbf{R},\nu}^{(1)}(\mathbf{r})$ is now given in Eq. (10). Equation (11) is also valid as an initial condition; i.e., the correction is included if at the initial time the nuclear velocity [the classical limit of $\lambda'_{\nu}(\mathbf{R},t)$] is nonzero [40].

Here, $\Phi_{\mathbf{R}}(\mathbf{r},t)$ is complex and can thus sustain an electronic current density [31,41] induced by the nuclear motion. The crucial point is that this current influences the nuclear motion through the TDVP.

It is worth noting here that Eq. (11) is a perturbative solution of the electronic equation (1). For small values of $\mu^2$, or equivalently of $\lambda'_{\nu}$, the effect of the ENCO is small; therefore, nonadiabatic effects can be treated as a perturbation to the ground-state dynamics. The theory is applicable in regions of weak nonadiabatic coupling, far away
from degeneracies between adiabatic potential energy surfaces.

### D. Expression of the time-dependent vector potential

The TDVP becomes nonzero when inserting Eq. (11) into Eq. (3). As described in more detail in Appendix C, up to within the linear order in $\mu^2$, the expression of the TDVP is

$$A_{\nu}(\mathbf{R}, t) = -2\hbar\mu^2 \int d\mathbf{r} \sum_{\nu'=1}^{N_n} [\lambda_{\nu'}(\mathbf{R}, t) \cdot \varphi_{\mathbf{R},\nu'}^{(1)}(\mathbf{r})] \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}). \quad (12)$$

Notice here that $\varphi_{\mathbf{R}}^{(0)}(\mathbf{r})$ and $\varphi_{\mathbf{R},\nu'}^{(1)}(\mathbf{r})$, as defined in Eq. (B7), are both real quantities, thus yielding a real TDVP. Using Eq. (10) and replacing the integral sign with $\langle \dots \rangle_{\mathbf{r}}$, we can define the quantity

$$\underline{\underline{A}}(\mathbf{R}) = 2\langle\varphi_{\mathbf{R}}^{(1)}|\hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R})|\varphi_{\mathbf{R}}^{(1)}\rangle_{\mathbf{r}}, \quad (13)$$

such that $\underline{\underline{A}}(\mathbf{R}, t) = -\mu^2 \underline{\underline{A}}(\mathbf{R}) \lambda'(\mathbf{R}, t)$. Here, the singly underlined symbols $\underline{\underline{A}}(\mathbf{R}, t)$, $\lambda'(\mathbf{R}, t)$, and $\varphi_{\mathbf{R}}^{(1)}$ indicate $(3N_n)$-dimensional vectors, whereas $\underline{\underline{A}}(\mathbf{R})$ is a $(3N_n \times 3N_n)$-dimensional matrix.

Since $\lambda'(\mathbf{R}, t)$ depends on $\underline{\underline{A}}(\mathbf{R}, t)$, we find $\underline{\underline{A}}(\mathbf{R}, t)$ self-consistently, which amounts to including an infinite number of terms of order $\mu^{2n}$. Recalling that $\underline{\underline{\lambda}} = \mu^2[-i\hbar\nabla\chi/\chi]$, the TDVP becomes

$$\underline{\underline{A}} = -\mu^2 \underline{\underline{A}} \underline{\underline{M}}^{-1} \underline{\underline{\lambda}}, \quad \text{with} \quad \underline{\underline{M}}(\mathbf{R}) = \underline{\underline{M}} + \mu^4 \underline{\underline{A}}(\mathbf{R}). \quad (14)$$

Here, $\underline{\underline{M}} = M_{\nu} \delta_{\nu,\nu'}$ is the $(3N_n \times 3N_n)$ diagonal mass matrix. If $\mu^4 = 1$, expressions where the physical masses appear are recovered. From Eq. (13), it is evident that $\underline{\underline{A}}(\mathbf{R})$ is a purely electronic quantity, which affects the nuclear momentum through the TDVP. Such correction, however, also appears in the nuclear evolution equation (1).

### E. Nuclear time-dependent Schrödinger equation

Having derived the electronic wave function up to order $\underline{\underline{O}}(\mu^2)$ and then the TDVP up to order $\underline{\underline{O}}(\mu^4)$, we can now determine an effective form for the nuclear Hamiltonian. On one side, the kinetic energy of the nuclei is reexpressed in a canonical form using the fact that the TDVP is shown to be proportional to $\underline{\underline{\lambda}}$, that is, to the nuclear momentum itself. On the other side, this involves determining the TDPES from the expansion of the electronic wave function, leading mainly to a term that describes the kinetic energy of the electrons induced by the nuclear motion. The detailed derivation is given in Appendix D. We then arrive at the nuclear time-dependent Schrödinger equation, written in matrix form,

$$\left[ \frac{1}{2} (-i\hbar\nabla)^T \underline{\underline{M}}^{-1}(\mathbf{R}) (-i\hbar\nabla) + E(\mathbf{R}) \right] \underline{\underline{\chi}} = i\hbar \partial_t \underline{\underline{\chi}}, \quad (15)$$

where the superscript T indicates the transpose vector and

$$E(\mathbf{R}) = \epsilon_{\text{BO}}^{(0)}(\mathbf{R}) + \sum_{\nu=1}^{N_n} \frac{\hbar^2}{2M_{\nu}} \langle \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} | \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}}. \quad (16)$$

The second term is the diagonal BO correction (DBOC). The kinetic energy term in Eq. (15) now involves dressed nuclear masses. It is important to notice that such a canonical form of the nuclear time-dependent Schrödinger equation arises from the self-consistent solution for $\underline{\underline{A}}$.

To better understand the meaning of this effective Hamiltonian, it is instructive to take the classical limit. The corresponding classical Hamiltonian [13] is simply $H_n = \underline{\underline{P}}^T \underline{\underline{M}}^{-1}(\mathbf{R}) \underline{\underline{P}}/2 + E(\mathbf{R})$, with nuclear velocity $\underline{\underline{R}} = \underline{\underline{M}}^{-1}(\mathbf{R}) \underline{\underline{P}}$. This Hamiltonian contains both the nuclear and electronic contributions to the kinetic energy, in the forms $\underline{\underline{R}}^T \underline{\underline{M}} \underline{\underline{R}}/2$ and $\underline{\underline{R}}^T \underline{\underline{A}}(\mathbf{R}) \underline{\underline{R}}/2$, respectively.

The key result of the paper is encoded in Eq. (15), where $\underline{\underline{M}}(\mathbf{R}) = \underline{\underline{M}} + \underline{\underline{A}}(\mathbf{R})$ since we have taken $\mu^4 = 1$. Even in the presence of (weak) nonadiabatic effects, the dynamical problem can be expressed in terms of nuclei moving on a single, static, potential energy surface—the electronic ground state (plus DBOC)—with masses that are corrected by the presence of the electrons. We have shown how, in a very simple and intuitive way, the electrons are carried along by the nuclei: $\underline{\underline{A}}(\mathbf{R})$, the $\underline{\underline{A}}$ matrix, is a position-dependent mass correction that dresses the bare nuclear masses $\underline{\underline{M}}$. The $\underline{\underline{A}}$ matrix is a purely electronic quantity, obtained by considering the lowest-order corrections $\underline{\underline{O}}(\mu^2)$ to the BO electronic wave function, and it appears both in the definition of the TDVP and in the nuclear Hamiltonian. The nuclear equation (15), together with the electronic equations (8) and (9), represents the limit of the exact factorization equations (1) up to $\underline{\underline{O}}(\mu^4)$ (see Appendix D). The $\underline{\underline{A}}$ matrix is the new and fundamental quantity introduced in this study, for which we are able to provide a rigorous derivation, in the context of the exact factorization; an intuitive interpretation, in terms of electronic mass carried along by the motion of the nuclei; and an efficient computation scheme, based on perturbation theory [23].

In the pioneering work of Bunker and Moss [2] on $H_2$ and $D_2$, the nuclear kinetic energy associated with the one-dimensional vibrational motion appears in the form $(2M_{\text{red}})^{-1}P(1+\beta(R))P$, with $M_{\text{red}}$ the reduced nuclear mass. Starting from the expression (15), it is straightforward algebra to rewrite $\underline{\underline{M}}^{-1}(\mathbf{R}) = (M + \underline{\underline{A}}(\mathbf{R}))^{-1}$ in the form $\underline{\underline{M}}^{-1}(1+\beta(R))$ by assuming that $\underline{\underline{A}}(\mathbf{R}) \ll \underline{\underline{M}}$. Further discussions on the connection between the
approach proposed in this work and previous analytical studies [12,13] are provided in Appendix C.

As we will now see, applying the correction to the mass matrix, through a resummation of an infinite number of terms $\mu^{2n}$, implies that the nuclear Hamiltonian has Galilean invariance.

### III. PROPERTIES OF THE DRESSED POSITION-DEPENDENT MASS CORRECTION

When Cartesian coordinates are employed as done here, the A matrix has the property of yielding the total electronic mass of the system when summed over all nuclei,

$$ \sum_{\nu, \nu'=1}^{N_n} \mathcal{A}_{\nu\nu'}^{ij}(\mathbf{R}) = m_e N_{el} \delta_{ij} \quad \forall \mathbf{R}, \quad (17) $$

supporting its interpretation as a correction term to the nuclear mass (indices $\nu$ and $\nu'$ run over the nuclei, and $i$ and $j$ over the three spatial dimensions). Here, $m_e$ is the electronic mass. It should also be noticed that the A matrix is positive definite in a ground-state dynamics [13]. The proof of Eq. (17) uses the property of the BO electronic wave function of being invariant under a translation of the reference system [13,42], and Eq. (9) (see Appendix E). This leads to

$$ \sum_{\nu, \nu'=1}^{N_n} [\underline{\underline{\mathcal{A}(\mathbf{R})}}]_{\nu\nu'} = \frac{m_e}{e} \sum_{\nu=1}^{N_n} [\underline{\underline{\mathcal{P}(\mathbf{R})}}]_{\nu} = m_e N_{el} \underline{\underline{\mathcal{I}}}^{(3)}, \quad (18) $$

where $[\underline{\underline{\mathcal{A}(\mathbf{R})}}]_{\nu\nu'}$ and $[\underline{\underline{\mathcal{P}(\mathbf{R})}}]_{\nu}$ are $(3 \times 3)$ matrices (in Cartesian components) and $\underline{\underline{\mathcal{I}}}^{(3)}$ is the identity matrix. Here, $[\underline{\underline{\mathcal{P}(\mathbf{R})}}]_{\nu} = \nabla_{\nu} \langle \hat{\mu}^{(\mathrm{el})}(\mathbf{R}) \rangle_{\mathrm{BO}}$ is the electronic contribution to the atomic polar tensor, defined as the variation with respect to nuclear positions of the electronic dipole moment (here averaged over the BO state) [43]. The second equality in Eq. (18) is obtained using the known property of the atomic polar tensor of yielding the total electronic charge of the system when summed over all nuclei [42,44].

It is common [2] to separate the c.m. motion before introducing the BO approximation. Calculations are then performed in internal coordinates. In this case, thus working within the molecular frame, the procedure presented here can be straightforwardly applied, by choosing coordinates in which the kinetic energy operator is the sum of two separated terms, i.e., nuclear and electronic. However, using the above sum rule, Eq. (18), it is possible to perform the calculations in Cartesian coordinates and to separate the c.m. motion *a posteriori*, still recovering the full mass of the system. In Cartesian coordinates, nuclear derivatives can be computed analytically via perturbation theory.

Starting from Cartesian coordinates, we make the following change of coordinates:

$$ \begin{gathered} \mathbf{R}'_1 = M_{\mathrm{tot}}^{-1} \left( \sum_{\nu=1}^{N_n} M_{\nu} \mathbf{R}_{\nu} + m_e \sum_{k=1}^{N_{\mathrm{el}}} \langle \hat{\mathbf{r}}_k \rangle_{\mathrm{BO}} \right), \\ \mathbf{R}'_{\nu} = \mathbf{R}_{\nu} - \mathbf{R}_1 \quad \text{with } \nu \ge 2, \end{gathered} \quad (19) $$

with $M_{\mathrm{tot}} = \sum_{\nu} M_{\nu} + m_e N_{\mathrm{el}}$. From the sum rule (18), the nuclear Hamiltonian of Eq. (15) becomes

$$ \hat{H}_n = \frac{\hat{P}_{\mathrm{c.m.}}^2}{2M_{\mathrm{tot}}} + \frac{1}{2} (-i\hbar\underline{\underline{\nabla}}')^T \underline{\underline{\underline{\mathcal{M}}}}^{-1}(\mathbf{R}') (-i\hbar\underline{\underline{\nabla}}') \\ + E(\mathbf{R}'). \quad (20) $$

Note that $\hat{P}_{\mathrm{c.m.}}$ is the momentum (operator) associated with the c.m. coordinate in Eq. (19); thus, the first term accounts for the motion of the c.m. as a free particle. The mass associated with the c.m. is, correctly, the total mass of the system, i.e., nuclei and electrons, rather than the nuclear mass only, as in the BO approximation. The following terms in Eq. (20) are the kinetic and potential energies corresponding to the internal, rotational and vibrational, degrees of freedom. Notice that, due to the fact that the Hamiltonian contains the gradient operator, the Jacobian of the transformation from Cartesian to internal coordinates has to be calculated. Therefore, the kinetic energy term in Eq. (15), $(-i\hbar\underline{\underline{\nabla}}')^T \underline{\underline{\underline{\mathcal{M}}}}^{-1}(-i\hbar\underline{\underline{\nabla}}')$, becomes $(-i\hbar\underline{\underline{\nabla}}')^T (\underline{\underline{\underline{J}}} \underline{\underline{\underline{\mathcal{M}}}}^{-1} \underline{\underline{\underline{J}}}) (-i\hbar\underline{\underline{\nabla}}')$. Here, $\underline{\underline{\underline{J}}} \underline{\underline{\underline{\mathcal{M}}}}^{-1} \underline{\underline{\underline{J}}}^T$ can be separated into a contribution from the c.m. [the first term on the right-hand side of Eq. (20)] and a contribution depending on the position-dependent mass matrix in internal coordinates [the second term of Eq. (20)]. Appendix F proves that off-diagonal contributions identically vanish with this choice of internal coordinates.

### IV. APPLICATIONS

The formalism introduced above is employed to construct a numerical procedure that (i) is fundamentally adiabatic, namely, only a single (static) potential energy surface is explicitly involved, but (ii) able to account for effects beyond BO via the position-dependent corrections to the bare nuclear masses.

The key quantity in the examples reported below is the nuclear Hamiltonian of Eq. (15). In the first numerical study of Sec. IVA, quantum mechanically, this Hamiltonian will be used to compute the spectrum of a model of a proton involved in a one-dimensional hydrogen bond [45]. In the same model system, interpreted classically, the Hamiltonian will be employed as the generator of the classical evolution of the oxygen atoms in the presence of a quantum proton. Furthermore, we have quantum mechanically computed the first four vibrational states of the whole system and compared them with the first four vibrational states of the reduced Hamiltonian. This model
system allows us to study the mass dependence of the error made by each approximation.

In Sec. IV B, the second numerical study presents results for some molecular systems, i.e., H₂, H₂O, NH₃, H₃O⁺, CH₄, and CH₃OH. The **A** matrix is computed for such molecules, proving that our approach can easily handle systems comprising more than three atoms. Furthermore, transforming the nuclear Hamiltonian of Eq. (15) to internal coordinates and within the harmonic approximation, position-dependent corrections are included in the calculation of the vibrational frequencies. Numerical details are given in Appendix G.

### A. Proton transfer

As a first application, we consider a model of a proton involved in a one-dimensional hydrogen bond O-H-O [45], in which nonadiabatic (vibrational) effects are known to be important [46]. The light particle is the proton, assumed to be in its vibrational ground state. In this case, everything that has been discussed on the electronic nature of the **A** matrix has to be adapted, as the **A** matrix now describes a mass correction to the heavy nuclei, the oxygens, because of the coupling with the proton. Hence, the **A** matrix is a purely protonic quantity. We use an asymmetric potential mimicking a strong hydrogen bond (as shown in Fig. 1): The proton is bonded to the oxygen atom O⁻ at large distances (we denote O⁻ the oxygen atom that is located on the left and O⁺ the one on the right), whereas at short distances, it is shared by the two oxygen atoms and is localized around the center of the O-O bond. The proton density corresponding to the ground state is shown in Fig. 2. At large distances, we expect the effective mass of O⁻ to be close to 17 a.m.u. as it carries the proton along. This is clear in Fig. 3, where it is shown that the element $A_{O^-O^-}(R)$ of the **A** matrix tends to a constant (equal to 1 a.m.u., the mass of the proton) at $R > 3$ Å, whereas all other components are zero, as

<figure><img src="image_2.png" /><figcaption>FIG. 1. Potential of the hydrogen bond model as a function of the O-O distance (R) and of the proton position (r) measured from the center of the O-O bond. The energy is measured in kcal/mol.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Proton density corresponding to the BO ground state, showing that in its vibrational ground state and for large O-O distances, the proton tends to be localized around the oxygen atom O⁻ because of the asymmetry of the potential.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Elements of the **A** matrix as functions of R. As expected, for a large O-O distance, the proton mass associated with the oxygen O⁻ tends to unity, while all other contributions tend to zero since the proton localizes around O⁻.</figcaption></figure>

expected from the sum rule of Eq. (18) (adapted to the present case—thus, with $m_e N_{el}$ replaced by $m_p$, the total proton mass of the system). We show this schematically in Fig. 4, where we plot the proton density along the O-O bond. We also report an estimate of the amount of proton mass associated with each oxygen, as the sum over the columns of the **A** matrix, e.g., $M_{O^-}(R) = M_{O^-} + [A_{O^-O^-}(R) + A_{O^+O^-}(R)]$. At short distances, the proton is instead shared by the oxygens: The elements of the **A** matrix are nonzero, but the O⁻ diagonal contribution remains dominant. Notice that it is not surprising that the off-diagonal elements of the **A** matrix are negative, as only two conditions are physically relevant: The diagonal elements must be non-negative, in a ground-state dynamics, and the sum of the elements must yield the proton mass, in
<figure><img src="image_5.png" /><figcaption>FIG. 4. Proton density at two values of the O—O distance (2.0 Å black and 2.8 Å red), where the masses of the oxygens (sum of columns of the matrix *M*, see text) *M*<sub>O+</sub> and *M*<sub>O-</sub> indicate the *A* matrix effect.</figcaption></figure>

a translationally invariant system. As seen in Fig. 4, the two oxygens then have similar masses at very short distances.

Figure 5 shows the classical trajectories of the two oxygen atoms starting from a compressed O—O distance and zero velocity. Calculations have been performed both in the standard adiabatic approximation (BO) and with position-dependent corrections to the oxygen masses

<figure><img src="image_6.png" /><figcaption>FIG. 5. Distance of the oxygens and position of the proton during the final steps of the dynamics: BO approximation with $M_{O^+} = M_{O^-} = 16$ a.m.u. (blue lines), BO approximation adding the proton mass $M_H = 1$ a.m.u. to the mass of $M_{O^-}$ (red lines), the BO approximation corrected by the position-dependent dressed mass (orange lines), and Ehrenfest dynamics (black lines). When position-dependent mass corrections are included in the BO approximation, the dynamics approaches reference results (Ehrenfest) that explicitly take into account nonadiabatic effects. As expected, the mass corrections yield a slower evolution than BO dynamics with bare oxygen masses.</figcaption></figure>

(BO + M). The two sets of calculations are compared with Ehrenfest dynamics, where nonadiabatic effects are included explicitly. Ehrenfest-type simulations, being explicitly nonadiabatic, require calculations of excited-state quantities, limiting not only the size of the accessible systems but also the time scales. BO + M calculations, based on the perturbation to the ground state, are instead easily affordable. The distance of the oxygens is plotted along with the mean position of the proton at the final steps of the dynamics. The masses are $M_{O^+} = M_{O^-} = 16$ a.m.u. and $M_{H^+} = 1$ a.m.u. In Fig. 6, it is shown that the c.m. of the system is perfectly fixed when position-dependent masses are employed, in contrast to the BO approximation. BO dynamics is faster than the Ehrenfest dynamics because the heavy atoms only have the bare nuclear mass. We have tested an *ad hoc* correction to the mass of the oxygen O<sup>-</sup>, i.e., $M_{O^-} = 17$ a.m.u. This improves the conservation of the c.m. but does not fix it completely. Changing $M_{O^-}$ to 17 a.m.u. improves the result, but only including the position-dependent dressed mass leads to a systematic convergence to the Ehrenfest results. We have further compared the error with respect to Ehrenfest dynamics, of BO and BO + M dynamics, as a function of the inverse mass ratio $\mu^{-4} = M_O/M_{H^+}$. This is shown in Fig. 7 as the root-mean-square deviation (RMSD) computed on the O—O distance along the trajectory with respect to the reference Ehrenfest trajectory. The position-dependent dressed mass greatly improves the precision of the dynamics even at small values of $\mu^{-4}$ (= 4 is the smallest value used) and leads to an error 4 orders of magnitude smaller than BO at large mass ratios. While at the physical value of the mass ratio, i.e., $\mu^{-4} = 16$, the population of the vibrational ground state along the Ehrenfest dynamics never decreases below 0.95, it reaches values as small as 0.8 for

<figure><img src="image_7.png" /><figcaption>FIG. 6. Position of the c.m. as a function of time. The color code is the same as that used in Fig. 5. Here, we show how the position of the center of mass is fixed when the position-dependent mass corrections are included in the BO approximation.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

FIG. 7. Root-mean-squared deviations of the O—O distance between Ehrenfest results and the BO approximation (black) or the BO approximation corrected by the position-dependent dressed mass (red). The results are shown as functions of the inverse mass ratio $\mu^{-4}$. The comparison shows that orders of magnitude in accuracy are gained with respect to the BO approximation, when the position-dependent mass correction is taken into account.

$\mu^{-4} = 4$. In this latter situation, nonadiabatic effects become substantial, but they are correctly captured by the BO + M.

Compared to Ehrenfest dynamics, the BO + M dynamics is much less computationally expensive, having a cost similar to the BO dynamics itself. Also, when interested in calculating electronic current densities induced by nuclear motion in molecules, an approach that takes into account excited-state effects is required to actually observe (in the simulation) a nonzero current [23,31]. The Ehrenfest approach is indeed suitable for such calculations, provided that all electronic excited states are accessible. The BO + M would indeed be preferred, as it yields (i) the correct electronic current density (employing the NVPT idea) and (ii) a correct description of the electronic effect on nuclear dynamics (via the *A* matrix), without the need for explicitly including all excited states, which might not be easy to compute anyway.

The proposed Hamiltonian formulation allows for a full quantum treatment of the nuclear dynamics: It maintains the simplicity of the BO approximation, making calculations of large systems feasible; at the same time, it gains accuracy.

To illustrate this, we have computed the four lowest eigenstates of the full quantum Hamiltonian at different values of $\mu^{-4}$. The diagonalization of the full Hamiltonian is compared to three approximations: BO, BO + DBOC, and BO + DBOC + M (where we also include the position-dependent correction). Figure 8 shows the error on the eigenvalues (the exact lowest eigenvalue is -4127.08527 cm⁻¹ at M<0xE2><0x82><0x95> = M<0xE2><0x82><0x95> = 16 a.m.u.). At small

<figure><img src="image_9.png" /></figure>

FIG. 8. Error between the four lowest eigenvalues of the full Hamiltonian and the BO approximation (solid lines with circles); the BO approximation corrected by adding the diagonal BO correction, BO + DBOC (dashed lines with crosses); and the BO approximation corrected by adding the diagonal BO correction and the position-dependent mass, BO + DBOC + M (dotted lines with squares). The results are shown as functions of the inverse mass ratio $\mu^{-4}$. As in Fig. 7, high precision is achieved when position-dependent mass corrections are taken into account.

$\mu^{-4}$, the BO approximation is expected to fail: The mass corrections allow us to gain 1 order of magnitude in the eigenvalues, even if compared to the case where the DBOC is included. Overall, in the static situation, the mass correction also leads to highly accurate results. At a mass ratio $\mu^{-4} = 1600$, an accuracy on the eigenvalues of about 10⁻⁵ cm⁻¹ is reached, whereas it is only 0.5 cm⁻¹ using the BO approximation. However, such an important gain of accuracy keeps the same favorable scaling and reduction of complexity as BO with respect to a full nonadiabatic treatment.

## B. Molecules: H₂, H₂O, NH₃, H₃O⁺, CH₄, CH₃OH

The major numerical result of this work is the possibility of computing from first principles the *A* matrix, a position-dependent mass correction to add to the bare nuclear masses in order to account for nonadiabatic effects. These position-dependent mass corrections were computed at the DFT level using analytical perturbation theory. The accuracy of DFT for this type of calculation can be asserted from the rotational *g* factors whose computation requires many similar terms. It has been shown that at equilibrium geometries, DFT gives very accurate results for the rotational *g* factors [47]. For the smaller molecules, H₂ and
$H_2O$, the results presented here are an improvement on the only available data in the literature [4] obtained at the self-consistent field (SCF) level with single excitations for excited states and numerical derivatives.

Once the *A* matrix is accessible, properties such as, but not limited to, corrections to harmonic vibrational frequencies can be determined. Therefore, we show in Tables I and II the *A* matrix for the molecules $H_2$ and $H_2O$, verifying the sum rule (17).

For larger molecules, the optimized geometries are listed in Appendix G 2, whereas the *A* matrix will be reported in Appendix H for readability reasons. A few comments can be made about these results. First, it appears that these matrices have significant off-diagonal contributions. The electrons give rise to inertial coupling between the nuclear motions that should not be neglected if we are aiming at high-accuracy molecular spectra. These off-diagonal terms are also rather long range in $CH_3OH$, as inertial coupling is found between the H(C) and O or even H(O) (see Table VI). Nevertheless, if we inspect the diagonal terms for the hydrogen atoms, we observe from $H_3O^+$ and $H_2O$, to $NH_3$, and then $CH_4$ or $CH_3OH$ a clear correlation with the electronegativity of the atom it is bonded to. While hydrogens bonded to O have an additional mass of about one half electron, the signature of the ionicity of the OH bond, it is nearly one full electron mass in the case of CH bonds.

<table><caption>TABLE I. Elements of the symmetric *A* matrix for the $H_2$ molecule at its equilibrium geometry, with the H–H bond length set to 0.743 Å and oriented along the *z* axis. The sum rule (17) yields (1.998, 1.998, 1.998) for *i*, *j* = *xx*, *yy*, *zz*, as the system contains two electrons.</caption><thead><tr><th></th><th colspan="2">Hydrogen 1</th><th colspan="2">Hydrogen 2</th></tr></thead><tbody><tr><td>Hydrogen 1</td><td>0.553</td><td></td><td>0.446</td><td></td></tr><tr><td></td><td></td><td>0.553</td><td></td><td>0.446</td></tr><tr><td></td><td></td><td></td><td>0.868</td><td>0.131</td></tr><tr><td>Hydrogen 2</td><td></td><td></td><td>0.553</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>0.553</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.868</td></tr></tbody></table>

In order to illustrate the effect of these position-dependent mass corrections on the vibrational frequencies, in Table III, we further report the values of the nonadiabatic corrections derived from the *A* matrix on vibrational frequencies for all molecules listed above. The corrected frequencies $\nu + \Delta\nu$ have been computed by diagonalizing the matrix $[M^{-1}(R_0)K(R_0)]$ at the equilibrium geometry $R_0$, where $K$ is the Hessian computed from the ground-state adiabatic potential energy surface. Negative frequency shifts will be computed: Nonadiabatic effects perturbing the ground-state dynamics tend to induce excitations of the light particles, and the energy necessary for the transition is “removed” from the heavy particles. Comparison with other theoretical predictions from the literature, when available, shows that the theory is capable of predicting accurate nonadiabatic corrections, even if working within the harmonic approximation (concerning the vibrational analysis) and with the generalized gradient approximation to DFT. Reference results [4] are determined at the SCF level, with excited states considered within a single-excitation CI approach and with an accurate quantum-mechanical calculation of anharmonic frequencies with the nuclear Hamiltonian.

What the approach developed in this study contributes to the field is the possibility of easily extending the numerical applications beyond di- and triatomic molecules [17–21,48]. To prove this, we provide the first predictions to the nonadiabatic corrections of vibrational frequencies of $NH_3$, $H_3O^+$, $CH_4$, and $CH_3OH$. It can be seen that the shifts of the N–H stretch frequencies of $NH_3$ are larger than those of the O–H stretch frequencies of $H_3O^+$, due to the fact that the N–H bonds are less ionic; as a result, the mass carried by the protons is larger in $NH_3$ than in $H_3O^+$ as observed above. The inertial coupling through the electronic motion also gives rise to different corrections for each vibrational mode. In particular, we observe that the stretch is more affected than bending. From these harmonic frequency shifts, it can be seen that correctly taking into account the electronic mass can lead to shifts of about 1 cm$^{-1}$. The mass matrices reported here should now allow for very accurate computations of molecular spectra.

<table><caption>TABLE II. The *A* matrix for the molecule $H_2O$ at the equilibrium geometry [$r_{OH} = 97.2$ pm and $\angle(\text{HOH}) = 104.5^\circ$]. The sum rule (17) yields (10.0004, 9.9997, 9.9999) for *i*, *j* = *xx*, *yy*, *zz*, as the system contains 10 electrons.</caption><thead><tr><th rowspan="2">A</th><th colspan="2">Oxygen</th><th colspan="2">Hydrogen 1</th><th colspan="2">Hydrogen 2</th></tr><tr><th>Oxygen</th><th></th><th>Oxygen</th><th></th><th>Oxygen</th><th></th></tr></thead><tbody><tr><td>Oxygen</td><td>8.0054</td><td></td><td>8.3000</td><td></td><td>0.3122</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>7.8025</td><td></td><td></td></tr><tr><td>Hydrogen 1</td><td></td><td></td><td></td><td></td><td>0.3365</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>0.3365</td></tr><tr><td>Hydrogen 2</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>TABLE III. Harmonic frequencies $\nu$ (in cm$^{-1}$) and their nonadiabatic corrections $\Delta\nu$. Benchmark values are taken from Ref. [4] when indicated. They have been calculated using a SCF approach for the ground state and including the nonadiabatic corrections in a single-excitation CI framework. Differences are observed in the predicted frequencies ($\nu$), which might be ascribed to the use of a different electronic structure approach and to the fact that in Ref. [4] the diagonal BO correction term has been included in the calculation of the ground-state potential for the vibrational analysis.</caption><thead><tr><th>Molecule</th><th>H<sub>2</sub></th><th>H<sub>2</sub> [4]</th><th>H<sub>2</sub>O</th><th>H<sub>2</sub>O [4]</th><th>NH<sub>3</sub></th><th>H<sub>3</sub>O<sup>+</sup></th><th>CH<sub>4</sub></th><th>CH<sub>3</sub>OH</th></tr></thead><tbody><tr><td>ν, Δν</td><td>4343.28, -0.89 -0.74</td><td>1594.93, -0.09 1597.60, -0.07</td><td>3656.19, -0.96 3661.00, -0.69</td><td>3757.77, -0.81 3758.63, -0.77</td><td>1016.73, -0.08</td><td>837.27, -0.05 1306.81, -0.12</td><td>297.67, -0.03</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>1628.30, -0.12</td><td>1639.25, -0.07 1522.48, -0.18</td><td>973.53, -0.14</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>3358.91, -0.97</td><td>3438.80, -0.47 2964.31, -0.68</td><td>1043.18, -0.16</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>3471.93, -0.89</td><td>3522.10, -0.37 3057.48, -0.88</td><td>1131.23, -0.19</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1331.39, -0.11</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1433.73, -0.15</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1458.72, -0.20</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1469.75, -0.20</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>2918.76, -0.87</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>2959.41, -1.12</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>3034.70, -0.91</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>3677.79, -1.06</td><td></td></tr></tbody></table>

### C. Methods

In the static calculations of the O-H-O model, the eigenvalues of the full Hamiltonian are determined using a Gaussian quadrature method with 20 points for R (the O-O distance) and 34 for r (the position of the proton from the center of the O-O bond). In the dynamics, we use the three coordinates, i.e., $R_O^+$, $R_O^-$, and $r = r_H$, in order to test the conservation of the position of the c.m.. The velocity-Verlet algorithm is used to integrate the classical nuclear equations, with a time step of 1 fs; the Crank-Nicolson [49] algorithm for the proton (quantum) equation in Ehrenfest dynamics, with a time step of $10^{-4}$ fs; the Euler algorithm with a time step of 0.0625 fs for BO + M calculations, where the force depends on the velocity. For the vibrational spectra, the A matrix has been computed using density functional perturbation theory [23,25,50], and it has been checked that the sum rule of Eq. (18) is satisfied. The numerical scheme has been implemented in the electronic structure package CPMD [51]. Calculations have been performed using Troullier-Martins [52] pseudopotentials in the Becke-Lee-Yang-Parr [53,54] (BLYP) approximation of the exchange-correlation kernel. The equilibrium molecular geometry is determined at the BLYP level, employing the aug-cc-pVTZ basis set [55] in the Gaussian electronic structure program [56].

### V. CONCLUSIONS

This work provides a rigorous theory to include the effect of light-particle motion on heavy-particle dynamics in molecules, within the adiabatic framework. The theory derived in this work has been developed by referring to coupled electron-nuclear systems, and we will adopt such an electron-nuclear perspective in the following discussion as well. However, we stress once again that the theory can

be straightforwardly applied to a more general situation, as proven by the numerical study of Sec. IVA.

Nuclear masses are dressed by position-dependent corrections that are purely electronic quantities and a consequence of the fact that electrons do not rigidly follow the motion of the nuclei. Various applications are discussed, yielding, in all cases, striking agreement with the benchmarks, either exact or highly accurate quantum-mechanical calculations. The idea of perturbatively including non-adiabatic electronic effects on the nuclear motion has been introduced previously [2,5,17,23,26,31], mainly as a tool to resolve some of the issues encountered in the context of theoretical vibrational spectroscopy when working in the BO approximation. Similarly, the idea of accounting for corrections to the nuclear masses has been proposed [2,5,12,13,57] to cure some fundamental inconsistencies of the BO treatment. The novelty of the present study is thus to be found in the overall picture that our work conveys: The theory is developed based on a rigorous starting point, the exact factorization of the molecular wave function; the perturbation treatment is justified in terms of the electron-nuclear mass ratio, as in the seminal paper of Born and Oppenheimer; the algebraic procedure is very simple, easily allowing for applications not restricted to di- and triatomic molecules; the proposed numerical scheme requires standard electronic structure calculations to determine the mass corrections, as the expressions of such corrections are explicitly given in terms of electronic properties. We expect that the theory will be able to provide solid information to predict and interpret highly accurate spectroscopy experiments in a large class of molecular systems.

Conceptually, we have resolved a well-known [58] fundamental inconsistency of the BO approximation. In a translationally invariant problem, the center of mass moves as a free particle with mass that equals the total mass of the
systems, i.e., nuclei and electrons, not only the nuclear mass. This feature is naturally built in the theory and corrects for a deficiency of the BO approximation, providing exactly the missing mass of the electrons. From a more practical point of view, our approach is very general and can be applied whenever a “factorization” of the underlying physical problem is possible, e.g., in the case of proton and oxygen atoms or in the case of electrons and nuclei.

Further applications are indeed envisaged since the perturbative incorporation of nonadiabatic effects greatly reduces the complexity of the fully coupled problem. For instance, the approximations can be applied to nuclear wave-packet methods for the calculation of highly accurate vibrational spectra beyond the BO approximation. The position-dependent mass is also shown to be related to the ionicity of the bonds and may serve as a proxy to access electronic properties.

## ACKNOWLEDGMENTS

The authors are grateful to Ari P. Seitsonen for providing them with local pseudopotentials for C, O, and N.

## APPENDIX A: ADIABATIC LIMIT OF THE EXACT FACTORIZATION

In this section, we argue that (i) the correct scaling of the time variable is $\mu^2$, when the parameter $\mu^4$ is used to scale the nuclear mass $M_\nu$, and (ii) the term in the electron-nuclear coupling operator of Eq. (2) containing the nuclear wave function scales with $\mu^2$ as well.

Statement (i) is obtained by taking the large mass limit, or small $\mu^4$ limit, as in Ref. [38]. In this situation, the dynamics of the heavy nuclei becomes semiclassical, and our scaling argument will make the nuclear kinetic energy tend towards a constant. In the classical limit, it is easy to see that at different values of $\mu^4$ the trajectories of the nuclei can be superimposed if the physical time $s$ is rescaled to a common time $t = \mu^2 s$ [27]. At each configuration $R(s)$ along the dynamics, the scaling of the time variable has the effect of yielding a kinetic energy that is a constant of $\mu^4$. In other words, the velocities $V(s)$ scale as $\mu^{-2}$. Notice that this is possible because we do not scale the positions with $\mu^4$, and therefore, the potential energy is not affected by the scaling. Using the common rescaled time $t$ to describe the nuclear trajectory, it then becomes possible to make a convergence statement about the nuclear dynamics.

Following Ref. [38], the nuclear wave packet can be considered to be a Gaussian wave packet localized at the position $\underline{R}(t)$, with momentum $\underline{P}(t)$:

$$ \chi(\underline{R}, t) = \pi^{-3N_n/4} \mu^{-3N_n/2} (\det \underline{\sigma}(t))^{1/4} \times e^{-\{\frac{-(\underline{R}-\underline{R}(t))^T \underline{\sigma}(t)(\underline{R}-\underline{R}(t))}{2\mu^2} + (i/\hbar)\{\underline{P}(t)\}/\mu^2\} \cdot (\underline{R}-\underline{R}(t))}, \quad (A1) $$

with $\underline{\sigma}(t)$ a $(3N_n \times 3N_n)$ symmetric matrix yielding the spatial extension of the wave packet. Notice here that both the amplitude and phase of the nuclear wave function are assumed to depend on $\mu^4$, and this dependence is inserted as $\mu^{-2}$ in the real and in the complex exponents. From this expression, we see that statement (ii) holds: $-i\hbar\nabla_\nu\chi/\chi$ scales as $\mu^{-2}$; thus, $\lambda_\nu(\underline{R}, t) = \mu^2\{[-i\hbar\nabla_\nu\chi(\underline{R}, t)]/\chi(\underline{R}, t)\}$ tends towards a quantity independent of $\mu$.

## APPENDIX B: NUCLEAR VELOCITY PERTURBATION THEORY

In this section, we show the relation between the $\mu^4$ expansion proposed in the paper and the nuclear velocity perturbation theory (NVPT) of Ref. [23]. We recall here the definition of $\lambda'_\nu(\underline{R}, t)$,

$$ \lambda'_{\nu}(\underline{R}, t) = \frac{1}{M_{\nu}} \left( \mu^2 \frac{-i\hbar\nabla_{\nu}\chi(\underline{R}, t)}{\chi(\underline{R}, t)} + \mu^2 \underline{A}_{\nu}(\underline{R}, t) \right) \quad (B1) $$

$$ = \frac{1}{M_\nu} (\lambda_\nu(\underline{R}, t) + \mu^2 \underline{A}_\nu(\underline{R}, t)). \quad (B2) $$

In the framework of NVPT, we have used $\lambda_\nu(\underline{R}, t)/M_\nu$ as the perturbation parameter that controls the degree of nonadiabaticity of the problem. The electronic equation (7) can be written using $\lambda'_\nu(\underline{R}, t)$ as

$$ \begin{aligned} & [\hat{H}_{\text{BO}} - e_{\text{BO}}^{(0)}(\mathbf{R})](\varphi_{\mathbf{R}}^{(0)}(\mathbf{r}) + \mu^2 \Phi_{\mathbf{R}}^{(1)}(\mathbf{r}, t)) \\ & = \mu^2 \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R}, t) \cdot [i\hbar\nabla_{\nu} + \underline{A}_{\nu}(\mathbf{R}, t)]\varphi_{\mathbf{R}}^{(0)}(\mathbf{r}). \end{aligned} \quad (B3) $$

Also, as will be proved in Appendix C, the TDVP is itself $O(\mu^2)$; thus, it will be neglected from the term in square brackets on the right-hand side. If we solve this equation order by order, Eqs. (8) and (9) are easily obtained. In particular, we recall here Eq. (9) whose solution yields $\Phi_{\mathbf{R}}^{(1)}(\mathbf{r}, t)$,

$$ \begin{aligned} & [\hat{H}_{\text{BO}} - e_{\text{BO}}^{(0)}(\mathbf{R})]\Phi_{\mathbf{R}}^{(1)}(\mathbf{r}, t) \\ & = i \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R}, t) \cdot (\hbar \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r})). \end{aligned} \quad (B4) $$

We have assumed that $\lambda'_\nu(\underline{R}, t)$ is small because of its dependence on $\mu^2$, the small parameter adopted in the perturbation approach. The correction $\Phi_{\mathbf{R}}^{(1)}(\mathbf{r}, t)$ to the electronic wave function can therefore be computed directly by inverting Eq. (B4) if $\lambda'_\nu(\underline{R}, t)$ itself is known. Moreover, the “resolution” of the perturbation in terms of $N_n$ (nuclear) contributions is somehow missing, which would make the numerical computation straightforward. To this end, we use below the same derivation proposed in Ref. [23], but we replace the small parameter used there,
i.e., $\lambda_{\nu}(\mathbf{R}, t)/M_{\nu}$, with the parameter introduced in this work, i.e., $\lambda'_{\nu}(\mathbf{R}, t)$: We assume that $\lambda'_{\nu}(\mathbf{R}, t)$ is itself a small quantity. The formal derivation is exactly the same as in Ref. [23], but the final result reported in Eq. (B9) will allow us to (i) decompose $\Phi_{\mathbf{R}}^{(1)}(\mathbf{r}, t)$ as a sum of $N_n$ independent contributions and (ii) compute each contribution without explicit knowledge of $\lambda'_{\nu}(\mathbf{R}, t)$.

We start from the electronic Hamiltonian of the form used in Ref. [23], replacing $\lambda_{\nu}(\mathbf{R}, t)/M_{\nu}$ with $\lambda'_{\nu}(\mathbf{R}, t)$,

$$ \hat{H}_{\text{el}} = \hat{H}_{\text{BO}} + \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R}, t) \cdot (-i\hbar\nabla_{\nu}), \quad (\text{B5}) $$

and we have solved it perturbatively, using $\hat{H}_{\text{BO}}$ as the unperturbed Hamiltonian. It is clear, as stated above, that $\lambda'_{\nu}(\mathbf{R}, t)$ is the small parameter that controls the strength of the perturbation and that $-i\hbar\nabla_{\nu}$ is the (nonadiabatic) perturbation. We have looked for the eigenstates of $\hat{H}_{\text{el}}$ in the form

$$ \Phi_{\mathbf{R}}(\mathbf{r}, t) = \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}) + \sum_{e \neq 0} \frac{\langle \varphi_{\mathbf{R}}^{(e)} | -i\hbar \sum_{\nu} \lambda'_{\nu}(\mathbf{R}, t) \cdot \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}}}{\epsilon_{\text{BO}}^{(0)}(\mathbf{R}) - \epsilon_{\text{BO}}^{(e)}(\mathbf{R})} \times \varphi_{\mathbf{R}}^{(e)}(\mathbf{r}), \quad (\text{B6}) $$

as straightforwardly follows from the application of standard time-independent perturbation theory. The first-order perturbation to the BO ground state can be written as

$$ i\varphi_{\mathbf{R},\nu}^{(1)}(\mathbf{r}) = i\sum_{e \neq 0} \frac{\mathbf{d}_{\nu,e0}(\mathbf{R})}{\omega_{e0}(\mathbf{R})} \varphi_{\mathbf{R}}^{(e)}(\mathbf{r}) \quad (\text{B7}) $$

with $\omega_{e0}(\mathbf{R}) = [\epsilon_{\text{BO}}^{(e)}(\mathbf{R}) - \epsilon_{\text{BO}}^{(0)}(\mathbf{R})]/\hbar$ and $\mathbf{d}_{\nu,e0}(\mathbf{R}) = \langle \varphi_{\mathbf{R}}^{(e)} | \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}}$, the nonadiabatic coupling vectors. This leads to a new expression of $\Phi_{\mathbf{R}}(\mathbf{r}, t)$,

$$ \Phi_{\mathbf{R}}(\mathbf{r}, t) = \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}) + i \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R}, t) \cdot \varphi_{\mathbf{R}}^{(1)}(\mathbf{r}), \quad (\text{B8}) $$

which is exactly Eq. (11) when setting $\mu^2 = 1$, to obtain the physical nuclear mass.

In the framework of NVPT, the perturbation parameter has been interpreted classically as the nuclear velocity [23,33,59]. It is worth mentioning here that, when performing a numerical simulation, such dependence on the nuclear velocity shall also be correctly accounted for in the preparation of the initial electronic state. When using NVPT to perform the calculations, the electronic evolution is not explicit, in the sense that at each time the electronic wave function is simply reconstructed using ground-state properties that are then inserted in Eq. (B8). However, when NVPT results are (or can be) compared with quantum-mechanical fully nonadiabatic results, the initial

electronic state cannot be simply prepared in the ground state, unless the initial nuclear velocity is zero. If this is not the case, then the first-order contribution in Eq. (B8), proportional to the finite value of the initial nuclear velocity, has to be included in the initial condition. Then, NVPT and nonadiabatic results can be directly compared, as the same initial conditions are used in both.

Equating the first-order corrections to the BO eigenstate, from the $\mu^4$ and the $\lambda'_{\nu}(\mathbf{R}, t)$ expansion, yields

$$ \Phi_{\mathbf{R}}^{(1)}(\mathbf{r}, t) = i \sum_{\nu=1}^{N_n} \lambda'_{\nu}(\mathbf{R}, t) \cdot \varphi_{\mathbf{R}}^{(1)}(\mathbf{r}). \quad (\text{B9}) $$

The comparison between the $\mu^4$ expansion and NVPT allows us, first of all, to derive an explicit expression of $\varphi_{\mathbf{R}}^{(1)}(\mathbf{r})$ as given in Eq. (B7) and, second, to decompose the perturbed state as a sum of independent (linear) responses to the nonadiabatic perturbations, thus leading to

$$ [\hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R})] \varphi_{\mathbf{R},\nu\alpha}^{(1)}(\mathbf{r}) = \hbar \partial_{\nu\alpha} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}). \quad (\text{B10}) $$

As above, the index $\nu$ is used to label the nuclei, and $\alpha$ labels the Cartesian components of the gradient. This equation can now be easily solved by employing density functional perturbation theory as described in Ref. [23].

## APPENDIX C: ANALYSIS OF THE PERTURBATION PARAMETER

The TDVP, defined in Eq. (3), is written using Eq. (11) as

$$ \mathbf{A}_{\nu}(\mathbf{R}, t) = \left\langle \varphi_{\mathbf{R}}^{(0)} + i\mu^2 \sum_{\nu'=1}^{N_n} \lambda'_{\nu'}(\mathbf{R}, t) \cdot \varphi_{\mathbf{R},\nu'}^{(1)} \right| \\ -i\hbar\nabla_{\nu}\varphi_{\mathbf{R}}^{(0)} + \mu^2\hbar\nabla_{\nu}\sum_{\nu'=1}^{N_n} \lambda'_{\nu'}(\mathbf{R}, t) \cdot \varphi_{\mathbf{R},\nu'}^{(1)} \rangle_{\mathbf{r}}. \quad (\text{C1}) $$

Up to within the linear order in $\mu^2$ [or, more precisely, $\mu^2\lambda'_{\nu}(\mathbf{R}, t)$], this expression is

$$ \mathbf{A}_{\nu}(\mathbf{R}, t) = -2\hbar\mu^2 \int d\mathbf{r} \sum_{\nu'=1}^{N_n} [\lambda'_{\nu'}(\mathbf{R}, t) \cdot \varphi_{\mathbf{R},\nu'}^{(1)}(\mathbf{r})] \nabla_{\nu'} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}), \quad (\text{C2}) $$

where we can use Eq. (B10) to identify the $\mathcal{A}$ matrix,

$$ \underline{\underline{\mathcal{A}}}(\mathbf{R}) = 2 \langle \varphi_{\mathbf{R}}^{(1)} | \hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R}) | \varphi_{\mathbf{R}}^{(1)} \rangle_{\mathbf{r}}. \quad (\text{C3}) $$

We derive the following expression of the TDVP, namely,

$$ \underline{\underline{\mathcal{A}}}(\mathbf{R}, t) = -\mu^2 \underline{\underline{\mathcal{A}}}(\mathbf{R}) \underline{\underline{\lambda}}'(\mathbf{R}, t). \quad (\text{C4}) $$
Once again, we keep the term $O(\mu^2)$ in $\lambda'$, but we show below how it will be included in the definition of the small parameter $\lambda$. Here, $\underline{\underline{A}}(\mathbf{R})$ is a matrix, thus the double- underlined notation, with $(3N_n \times 3N_n)$ elements, whereas $\varphi_{\mathbf{R}}^{(1)}(\mathbf{r})$ is a vector with $(3N_n)$ components (and each component is a function of $3N_e$ electronic coordinates, $\mathbf{r}$, with parametric dependence on $3N_n$ nuclear coordinates, $\mathbf{R}$). We have also written the TDVP and the parameter in matrix notation, with $\underline{\underline{A}}(\mathbf{R}, t)$ and $\underline{\underline{\lambda}}'(\mathbf{R}, t)$ $(3N_n)$- dimensional vectors. The elements of the $\mathcal{A}$ matrix are

$$ \underline{\underline{\mathcal{A}}}_{\nu\nu}^{ij}(\mathbf{R}) = \langle \varphi_{\mathbf{R},\nu i}^{(1)} | \hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R}) | \varphi_{\mathbf{R},\nu j}^{(1)} \rangle_{\mathbf{r}}, \quad (\text{C5}) $$

with i, j labeling the Cartesian components and $\nu'$, $\nu$ the nuclei. When using Eq. (B7), the elements of the $\mathcal{A}$ matrix can be written in terms of the nonadiabatic coupling vectors and of the BO eigenvalues as

$$ \underline{\underline{\mathcal{A}}}_{\nu\nu}^{ij}(\mathbf{R}) = 2\hbar \sum_{e \neq 0} \frac{d_{\nu',e0}(\mathbf{R}) d_{\nu,j,e0}(\mathbf{R})}{\omega_{e0}(\mathbf{R})} \quad (\text{C6}) $$

from which it follows that the $\mathcal{A}$ matrix is symmetric. The $\mathcal{A}$ matrix is also positive definite (i.e., for all nonzero real vectors $\mathbf{v}$, the relation $\mathbf{v}^T \underline{\underline{\mathcal{A}}} \mathbf{v} \ge 0$ holds) with non-negative diagonal elements, i.e.,

$$ \underline{\underline{\mathcal{A}}}_{\nu\nu}^{ii}(\mathbf{R}) = 2\hbar \sum_{e \neq 0} \frac{|\underline{\underline{d}}_{\nu,i,e0}(\mathbf{R})|^2}{\omega_{e0}(\mathbf{R})} \ge 0. \quad (\text{C7}) $$

This property is essential for the interpretation of the $\mathcal{A}$ matrix as a position-dependent mass correction. The components of the TDVP can be expressed in terms of the components of the $\mathcal{A}$ matrix,

$$ A_{\nu i}(\mathbf{R}, t) = -\mu^2 \sum_{\nu'=1}^{N_n} \sum_{j=x,y,z} \underline{\underline{\mathcal{A}}}_{\nu\nu'}^{ij}(\mathbf{R}) \lambda'_{\nu'j}(\mathbf{R}, t). \quad (\text{C8}) $$

This expression is used in the definition of the parameter $\lambda'_{\nu i}(\mathbf{R}, t)$, given in Eq. (B1),

$$ \lambda'_{\nu i}(\mathbf{R}, t) = M_{\nu}^{-1} \lambda_{\nu i}(\mathbf{R}, t) - \mu^4 M_{\nu}^{-1} \sum_{\nu', j} \underline{\underline{\mathcal{A}}}_{\nu\nu'}^{ij}(\mathbf{R}) \lambda'_{\nu'j}(\mathbf{R}, t), \quad (\text{C9}) $$

where

$$ \lambda_{\nu i}(\mathbf{R}, t) = \mu^2 \frac{-i\hbar\partial_{\nu i}\chi(\mathbf{R}, t)}{\chi(\mathbf{R}, t)}, \quad (\text{C10}) $$

which, we recall, tends towards a quantity independent of $\mu$ if $\mu \to 0$.

Writing Eq. (C9) in matrix form and solving for $\underline{\underline{\lambda}}(\mathbf{R}, t)$, we obtain

$$ \underline{\underline{\lambda}}(\mathbf{R}, t) = [\underline{\underline{M}} + \mu^4 \underline{\underline{\underline{A}}}(\mathbf{R})] \underline{\underline{\lambda}}'(\mathbf{R}, t) = \underline{\underline{\underline{M}}} (\mathbf{R}) \underline{\underline{\lambda}}'(\mathbf{R}, t), \quad (\text{C11}) $$

where $\underline{\underline{M}}$ is a diagonal $(3N_n \times 3N_n)$ matrix containing the masses of the nuclei and we have defined a position- dependent mass matrix $\underline{\underline{M}}(\mathbf{R})$. This equation can be inverted to obtain

$$ \underline{\underline{\lambda}}'(\mathbf{R}, t) = \underline{\underline{\underline{M}}}^{-1}(\mathbf{R}) \underline{\underline{\lambda}}(\mathbf{R}, t), \quad (\text{C12}) $$

yielding the TDVP in the form given in Eq. (14),

$$ \underline{\underline{\underline{A}}}(\mathbf{R}, t) = -\mu^2 \underline{\underline{\underline{A}}}(\mathbf{R}) \underline{\underline{\underline{M}}}^{-1}(\mathbf{R}) \underline{\underline{\lambda}}(\mathbf{R}, t), \quad (\text{C13}) $$

with $\mu^4 = 1$, where only $\underline{\underline{\lambda}}$ appears.

Equation (C8) shows that the TDVP is at least first order in the perturbation parameter, and this is the reason why it is not considered in the definition of the perturbed electronic Hamiltonian in Eq. (B5). Because of the explicit dependence of $\underline{\underline{A}}_{\nu}(\mathbf{R}, t)$ on $\underline{\underline{\lambda}}'_{\nu}(\mathbf{R}, t)$, which is known via the $\mathcal{A}$ matrix, we have been able to isolate the “actual” small parameter, i.e., $\underline{\underline{\lambda}}(\mathbf{R}, t)$. In all expressions, however, we find $\underline{\underline{\lambda}}'(\mathbf{R}, t)$, the matrix product of $\underline{\underline{\underline{M}}}^{-1}(\mathbf{R})$ and $\underline{\underline{\lambda}}(\mathbf{R}, t)$, which is a gauge-invariant quantity.

As anticipated in Sec. II E, we now draw connections between the position-dependent mass corrections derived in this work and Refs. [12,13]. Starting from a perturbed Hamiltonian of the type introduced in Ref. [23], namely, where the perturbation is driven by the nuclear velocity, Goldhaber [13] identifies a position-dependent inertia tensor, also discussed there as a metric tensor, exactly as in Eq. (C6). The tensor diverges as the third power of the energy difference between the electronic ground state and the excited states, preventing an adiabaticlike dynamics to drive nuclear motion too close to a region of degeneracy. Somehow different are the numerical results proposed by Kutzelnigg [12]. While the idea of assigning position- dependent mass corrections to the vibrational and rotational nuclear degrees of freedom of a diatomic molecule, i.e., $H_2^+$ and $H_2$, is similar to the one presented in this work (by considering the variation of the electronic ground state with respect to the nuclear position), the numerical results are not in accordance with our findings. The reason for this is probably because of the level of electronic structure theory considered in Ref. [12], i.e., LCAO with noninteracting electrons.

## APPENDIX D: NUCLEAR HAMILTONIAN

In this section, we describe the procedure leading to the appearance of the position-dependent mass $\underline{\underline{M}}(\mathbf{R})$ in the nuclear evolution equation (1) of the exact factorization. Taking into account the $\mu^4$ scaling, the action of the kinetic energy operator $\hat{T}_n = \sum_{\nu} \mu^4 [-i\hbar\nabla_{\nu} + \underline{\underline{A}}_{\nu}]^2 / (2M_{\nu})$ on the nuclear wave function $\chi(\mathbf{R}, t)$ can be written in matrix form as
$$ \hat{\tilde{T}}_n \chi = \frac{1}{2} [-i\hbar\mu^2\underline{\nabla} + \mu^2\underline{A}]^T \underline{M}^{-1} [-i\hbar\mu^2\underline{\nabla} + \mu^2\underline{A}] \chi. \quad (D1) $$

As discussed previously, the terms $-i\hbar\mu^2\underline{\nabla}$ are $O(1)$. From Eq. (C13) for the TDVP, the term $\mu^2\underline{A} = -\mu^4\underline{A}\underline{M}^{-1}\underline{\lambda}$ is of order $O(\mu^4)$. Using expression (C13), we then identify the following terms:

$$ \hat{\tilde{T}}_n \chi = \frac{1}{2} [(-i\hbar\mu^2\underline{\nabla})^T \underline{M}^{-1} (\underline{\underline{T}} - \mu^4\underline{A}\underline{M}^{-1})(-i\hbar\mu^2\underline{\nabla}) \\ - (\mu^4\underline{A}\underline{M}^{-1}\underline{\lambda})^T \underline{M}^{-1} \underline{M}(\underline{M}^{-1}\underline{\lambda}) \\ + (\mu^4\underline{A}\underline{M}^{-1}\underline{\lambda})^T \underline{M}^{-1}(\mu^4\underline{A}\underline{M}^{-1}\underline{\lambda})] \chi. \quad (D2) $$

In the second term on the right-hand side, we have again used the definition of $\underline{\lambda}$ to write $-i\hbar\mu^2\underline{\nabla}\chi = \underline{\lambda}\chi$, and we have inserted the definition of the identity matrix in the form $\underline{T} = \underline{M}^{-1}\underline{M}$. We recall the expression of the position-dependent mass matrix, $\underline{M} = \underline{M} + \mu^4\underline{A}$, leading to the kinetic energy operator in the nuclear Hamiltonian (1),

$$ \hat{\tilde{T}}_n \chi = \frac{1}{2} (-i\hbar\mu^2\underline{\nabla})^T \underline{M}^{-1} (-i\hbar\mu^2\underline{\nabla}) \chi \\ - \frac{1}{2} \mu^4 (\underline{M}^{-1}\underline{\lambda})^T \underline{A} (\underline{M}^{-1}\underline{\lambda}) \chi, \quad (D3) $$

where only the position-dependent mass (and mass correction) appears. In the second term on the right-hand side, we have used the property of the $\mathcal{A}$ matrix of being symmetric; thus, $\underline{A}^T = \underline{A}$. It is interesting to note that no approximation has been invoked in deriving Eq. (D3) from Eq. (D1). The second term in Eq. (D3) is $O(\mu^4)$, and there is no additional contribution of the same order in $\hat{\tilde{T}}_n$.

We now turn to the TDPES, $\epsilon(\mathbf{R}, t)$, appearing in the nuclear Hamiltonian from the exact factorization, in Eq. (1). We shall see that in its expansion up to order $\mu^4$, a kineticlike contribution, which can be identified as the induced electronic kinetic energy due to the nuclear motion, exactly balances the second term in Eq. (D3).

We first write the expression of $\langle \Phi_R(t) | \hat{H}_{BO} | \Phi_R(t) \rangle_r$ up to terms in $\mu^4$, when the electronic wave function is expanded as

$$ \Phi_R(\mathbf{r}, t) = \varphi_R^{(0)}(\mathbf{r}) + \mu^2 \lambda'(t) \varphi_R^{(1)}(\mathbf{r}) \\ + \mu^4 \lambda'^2(t) \varphi_R^{(2)}(\mathbf{r}). \quad (D4) $$

Here, we use a simplified notation, also using the property that the only time dependence in the electronic wave function appears via $\lambda'_\nu(\mathbf{R}, t)$. Using this form of the electronic wave function, we write

$$ \begin{aligned} \langle \Phi_R(t) | \hat{H}_{BO} | \Phi_R(t) \rangle_r \\ = \epsilon_{BO}^{(0)}(\mathbf{R}) + \mu^4 \lambda'^2(t) \langle \varphi_R^{(1)} | \hat{H}_{BO} | \varphi_R^{(1)} \rangle_r \\ + \mu^4 \lambda'^2(t) \epsilon_{BO}^{(0)}(\mathbf{R}) [\langle \varphi_R^{(2)} | \varphi_R^{(0)} \rangle_r + \langle \varphi_R^{(0)} | \varphi_R^{(2)} \rangle_r] \\ + O(\mu^6), \end{aligned} \quad (D5) $$

and, by using the partial normalization condition up to within second order,

$$ \langle \varphi_R^{(0)} | \varphi_R^{(0)} \rangle_r + \mu^4 \lambda'^2(t) \langle \varphi_R^{(1)} | \varphi_R^{(1)} \rangle_r + \mu^4 \lambda'^2(t) \langle \varphi_R^{(2)} | \varphi_R^{(0)} \rangle_r \\ + \mu^4 \lambda'^2(t) \langle \varphi_R^{(0)} | \varphi_R^{(2)} \rangle_r = 1, \quad (D6) $$

we find

$$ \langle \varphi_R^{(2)} | \varphi_R^{(0)} \rangle_r + \langle \varphi_R^{(0)} | \varphi_R^{(2)} \rangle_r = - \langle \varphi_R^{(1)} | \varphi_R^{(1)} \rangle_r \quad (D7) $$

since the normalization condition is already satisfied at zeroth order. We insert this result in Eq. (D5) to obtain

$$ \begin{aligned} \langle \Phi_R(t) | \hat{H}_{BO} | \Phi_R(t) \rangle_r \\ = \epsilon_{BO}^{(0)}(\mathbf{R}) + \mu^4 \lambda'^2(t) \langle \varphi_R^{(1)} | \hat{H}_{BO} - \epsilon_{BO}^{(0)}(\mathbf{R}) | \varphi_R^{(1)} \rangle_r \\ + O(\mu^6). \end{aligned} \quad (D8) $$

In the second term on the right-hand side, we identify the $\mathcal{A}$ matrix, and we thus write

$$ \begin{aligned} \langle \Phi_R(t) | \hat{H}_{BO} | \Phi_R(t) \rangle_r \\ = \epsilon_{BO}^{(0)}(\mathbf{R}) + \mu^4 \sum_{\nu, \nu'} \sum_{i, j} \frac{1}{2} \lambda'_{\nu i}(\mathbf{R}, t) A_{\nu \nu'}^{ij}(\mathbf{R}) \lambda'_{\nu' j}(\mathbf{R}, t) \end{aligned} \quad (D9) $$

$$ = \epsilon_{BO}^{(0)}(\mathbf{R}) + \mu^4 \frac{1}{2} \lambda'^T(\mathbf{R}, t) \underline{A}(\mathbf{R}) \underline{\lambda'}(\mathbf{R}, t), \quad (D10) $$

where Eq. (D10) is Eq. (D9) in matrix form. Inserting the expression of $\underline{\lambda'}(\mathbf{R}, t)$ in terms of $\underline{\lambda}(\mathbf{R}, t)$ given in Eq. (C12), we can express the second term of Eq. (D10) as

$$ \begin{aligned} \underline{\lambda'}^T(\mathbf{R}, t) \underline{A}(\mathbf{R}) \underline{\lambda'}(\mathbf{R}, t) \\ = [\underline{M}^{-1}(\mathbf{R}) \underline{\lambda}(\mathbf{R}, t)]^T \underline{A}(\mathbf{R}) [\underline{M}^{-1}(\mathbf{R}) \underline{\lambda}(\mathbf{R}, t)], \end{aligned} $$

which exactly cancels the second term on the right-hand side of Eq. (D3).

The TDPES for the nuclei contains a further term of order $\mu^4$ [see Eq. (6)] referred to as the diagonal BO correction (DBOC). Altogether, the nuclear Hamiltonian of Eq. (15) is thus derived (setting $\mu^4 = 1$),

$$ \hat{H}_n = \frac{1}{2} (-i\hbar\underline{\nabla})^T \underline{M}^{-1}(\mathbf{R}) (-i\hbar\underline{\nabla}) + E(\mathbf{R}). \quad (D11) $$
The potential energy is time independent and contains the BO energy, from the first term in Eq. (D10), and the DBOC, according to

$$E(\mathbf{R}) = \epsilon_{\text{BO}}^{(0)}(\mathbf{R}) + \sum_{\nu=1}^{N_n} \frac{\hbar^2}{2M_{\nu}} \langle \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} | \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}}. \quad (\text{D12})$$

It is worth noting that the first-order contribution to the time-dependent potential $\epsilon(\mathbf{R}, t)$ is zero; thus, only $\epsilon^{(0)}(\mathbf{R})$, the zeroth-order term, appears as potential energy in the nuclear Hamiltonian of Eq. (15). This statement has already been proven in Ref. [23] using the definition in Eq. (4) and the expression of the electronic wave function up to within first-order terms in the perturbation. In between the DBOC and the position-dependent correction to the mass, all terms of order $\mu^4$ are included in the nuclear Hamiltonian given in Eq. (D11).

The correspondence principle of quantum mechanics enables us to determine the classical nuclear Hamiltonian as

$$H_n = \frac{1}{2} P^T \underline{M}^{-1}(\mathbf{R}) P + E(\mathbf{R}), \quad (\text{D13})$$

where $P = \underline{M}(\mathbf{R})\dot{\mathbf{R}}$ is the nuclear momentum.

## APPENDIX E: ELECTRONIC MASS AND THE A MATRIX

The derivation of Eq. (17) uses the property of the BO electronic wave function of being invariant under a translation of the coordinate reference system, namely, $\varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}') = \varphi_{\mathbf{R}}^{(0)}(\mathbf{r})$, with $\mathbf{R}' = \mathbf{R}'_1, \dots, \mathbf{R}'_{N_n} = \mathbf{R}_1 + \eta\Delta, \dots, \mathbf{R}_{N_n} + \eta\Delta$, and analogously for $\mathbf{r}'$. Notice that $\Delta$ is a three-dimensional vector and that all positions, electronic and nuclear, are shifted of the same amount of $\eta\Delta$. Translational invariance [42] means

$$\begin{aligned} 0 &= \frac{\partial \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}')}{\partial \eta} \\ &= \sum_{i=x,y,z} \left[ \sum_{\nu=1}^{N_n} \frac{\partial \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}')}{\partial R'_{\nu i}} \frac{\partial R'_{\nu i}}{\partial \eta} + \sum_{k=1}^{N_{\text{el}}} \frac{\partial \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}')}{\partial r'_{ki}} \frac{\partial r'_{ki}}{\partial \eta} \right] \\ &= \sum_{i=x,y,z} \Delta_i \left[ \sum_{\nu=1}^{N_n} \frac{\partial \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}')}{\partial R'_{\nu i}} + \sum_{k=1}^{N_{\text{el}}} \frac{\partial \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}')}{\partial r'_{ki}} \right] \\ &= \mathbf{\Delta} \cdot \left[ \sum_{\nu=1}^{N_n} \nabla_{\nu} \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}') + \sum_{k=1}^{N_{\text{el}}} \nabla_k \varphi_{\mathbf{R}'}^{(0)}(\mathbf{r}') \right], \end{aligned} \quad (\text{E1})$$

which is valid for all values of $\mathbf{\Delta}$. Identifying $\nabla_k$ as the position representation of the momentum operator $\hat{\mathbf{p}}_k$ corresponding to the $k$th electron (divided by $-i\hbar$), which can also be written as

$$\hat{\mathbf{p}}_k = \frac{im}{\hbar} [\hat{H}, \hat{\mathbf{r}}_k] = \frac{im}{\hbar} [\hat{H}_{\text{BO}}, \hat{\mathbf{r}}_k], \quad (\text{E2})$$

and projecting the two terms in square brackets in Eq. (E1) onto $\varphi_{\mathbf{R},\nu i}^{(1)}(\mathbf{r})$, from Eq. (B7),

$$\sum_{\nu=1}^{N_n} \langle \varphi_{\mathbf{R},\nu i}^{(1)} | \hbar \nabla_{\nu} \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}} = \frac{m}{\hbar} \sum_{k=1}^{N_{\text{el}}} \langle \varphi_{\mathbf{R},\nu i}^{(1)} | [\hat{H}_{\text{BO}}, \hat{\mathbf{r}}_k] | \varphi_{\mathbf{R}}^{(0)} \rangle, \quad (\text{E3})$$

we identify the A matrix on the left-hand side and, for each Cartesian component $j$, we write

$$\sum_{\nu=1}^{N_n} \mathcal{A}_{\nu\nu'}^{ij}(\mathbf{R}) = -\frac{m}{\hbar} \sum_{k=1}^{N_{\text{el}}} \langle \varphi_{\mathbf{R},\nu i}^{(1)} | [\hat{H}_{\text{BO}}, \hat{\mathbf{r}}_{kj}] | \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}}. \quad (\text{E4})$$

From the term on the right-hand side, we derive the expression of the atomic polar tensor (APT). First, we write the commutator explicitly, and we use Eq. (B10) to obtain

$$\begin{aligned} & \langle \varphi_{\mathbf{R},\nu i}^{(1)} | [\hat{H}_{\text{BO}}, \hat{\mathbf{r}}_{kj}] | \varphi_{\mathbf{R}}^{(0)} \rangle_{\mathbf{r}} \\ &= \int d\mathbf{r} \varphi_{\mathbf{R},\nu i}^{(1)}(\mathbf{r}) [\hat{H}_{\text{BO}} - \epsilon_{\text{BO}}^{(0)}(\mathbf{R})] r_{kj} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}) \end{aligned} \quad (\text{E5})$$

$$= -\hbar \int d\mathbf{r} (\partial_{\nu i} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r})) r_{kj} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}). \quad (\text{E6})$$

Then, we identify the expectation value of the electronic dipole moment operator over the BO wave function in the following expression:

$$\partial_{\nu i} \sum_{k=1}^{N_{\text{el}}} \int d\mathbf{r} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}) r_{kj} \varphi_{\mathbf{R}}^{(0)}(\mathbf{r}) = \frac{1}{e} \partial_{\nu i} \langle \hat{\mu}_j^{(\text{el})}(\mathbf{R}) \rangle_{\text{BO}}. \quad (\text{E7})$$

The derivative with respect to the $i$th Cartesian component, relative to the $\nu$th nucleus, of the $j$th Cartesian component of the electronic dipole moment is the definition of the electronic contribution to the APT [43] $\mathcal{P}_{ij}^{\nu}(\mathbf{R})$. This leads to the relation [42,44]

$$\sum_{\nu, \nu'=1}^{N_n} \mathcal{A}_{\nu\nu'}^{ij}(\mathbf{R}) = \sum_{\nu=1}^{N_n} \frac{m}{e} \mathcal{P}_{ij}^{\nu}(\mathbf{R}) = m N_{\text{el}} \delta_{ij}, \quad (\text{E8})$$

when we further sum over the index $\nu$. This result states that when the A matrix is summed up over all nuclei, it yields the total electronic mass of the complete system. In Eqs. (15) and (D11), this means that the mass effect of the electrons is completely taken into account by the position-dependent mass corrections to the nuclear masses within the order of the perturbation considered here.
APPENDIX F: SEPARATION OF THE CENTER OF MASS

We introduce the coordinate transformation

$$ \mathbf{R}'_1 = \mathcal{R}_{\text{c.m.}} = \frac{1}{M_{\text{tot}}} \left[ \sum_{\nu=1}^{N_n} M_{\nu} \mathbf{R}_{\nu} + m \sum_{k=1}^{N_{\text{el}}} \langle \hat{\mathbf{r}}_k \rangle_{\text{BO}} \right] \quad (\text{F1}) $$

$$ \mathbf{R}'_{\nu} = \mathbf{R}_{\nu} - \mathbf{R}_{1} \quad \text{with} \quad \nu \ge 2, \qquad (\text{F2}) $$

with the position of the c.m. defined in Eq. (19) and $M_{\text{tot}} = \sum_{\nu} M_{\nu} + mN_{\text{el}}$ the total mass of the system. Such coordinate transformation is applied to the kinetic and potential energy terms in the nuclear Hamiltonian (D11). Since we need to evaluate the gradient of $\chi$, we have to compute the Jacobian matrix of the transformation from Cartesian to internal coordinates. The Jacobian is a $(3N_n \times 3N_n)$ matrix, whose elements are

$$ J_{\nu\nu'}^{ij} = \frac{\partial R'_{\nu i}}{\partial R_{\nu j}} = \begin{cases} \frac{1}{M_{\text{tot}}} \left( M_{\nu'} \delta_{ij} + \frac{m}{e} P_{ji}^{\nu'} \right) & \text{if } \nu = 1 \\ -\delta_{1\nu'} \delta_{ij} + \delta_{\nu\nu'} \delta_{ij} & \text{if } \nu \ge 2, \end{cases} \quad (\text{F3}) $$

with $P_{ji}^{\nu'}$ the electronic APT of Eq. (E7). It can be proved with some simple, but tedious, algebra that the determinant of the Jacobian is unity. In Eq. (D11), we replace $\nabla$ with $\nabla'$ according to

$$ \begin{aligned} (-i\hbar\nabla)^T \underline{\underline{\underline{\underline{\underline{M}}}}}^{-1} (-i\hbar\nabla) &= [\underline{\underline{\underline{\underline{\underline{J}}}}}^T (-i\hbar\nabla')]^T \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1} [\underline{\underline{\underline{\underline{\underline{J}}}}}^T (-i\hbar\nabla')] \\ &= (-i\hbar\nabla')^T (\underline{\underline{\underline{\underline{\underline{J}}}}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1} \underline{\underline{\underline{\underline{\underline{\underline{J}}}}}}^T) (-i\hbar\nabla'), \end{aligned} \quad (\text{F4}) $$

where the position-dependent mass in the last term on the right-hand side depends on $\mathbf{R}'$, namely,

$$ \underline{\underline{\underline{\underline{\underline{M}}}}}^{-1}(\mathbf{R}') = \underline{\underline{\underline{\underline{\underline{J}}}}} \underline{\underline{\underline{\underline{\underline{M}}}}}^{-1}(\mathbf{R}) \underline{\underline{\underline{\underline{\underline{J}}}}}^T. \quad (\text{F5}) $$

We rewrite the Jacobian matrix as the sum of two terms, $\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}$ and $\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}}$: The first three rows of $\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}$ are the same as $\underline{\underline{\underline{\underline{\underline{J}}}}}$, and thus given by Eq. (F3) for $\nu = 1$, i.e., $(\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}})_{\nu\nu'}^{ij} = \delta_{\nu 1} \underline{\underline{\underline{\underline{\underline{J}}}}}^{ij}$, with each row composed of $3N_n$ entries; all other elements of $\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}$ are zeros. The first three rows of $\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}}$ are zero, and the remaining $3(N_n - 1)$ rows are the same as $\underline{\underline{\underline{\underline{\underline{J}}}}}$, and thus given by the second expression in Eq. (F3). We now introduce the operator $\underline{\underline{\underline{\underline{\underline{T}}}}}$, defined as $\underline{\underline{\underline{\underline{T}}}}_{\nu\nu'}^{ij} = \delta_{ij}\delta_{\nu 1}$, and we notice that the product of the position-dependent mass matrix and $\underline{\underline{\underline{\underline{T}}}}$ yields

$$ \underline{\underline{\underline{\underline{\underline{M}}}}}(\mathbf{R})\underline{\underline{\underline{\underline{T}}}} = M_{\text{tot}}[\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}]^T, \quad (\text{F6}) $$

as we will now prove. First, we recall the expression of the position-dependent mass matrix,

$$ \underline{\underline{\underline{\underline{\underline{M}}}}}^{ij}_{\nu\nu'}(\mathbf{R}) = M_{\nu}\delta_{\nu\nu'}\delta_{ij} + \underline{\underline{\underline{\underline{\underline{A}}}}}^{ij}_{\nu\nu'}(\mathbf{R}), \quad (\text{F7}) $$

and then we write the matrix product with $\underline{\underline{\underline{T}}}$ as the sum of their components, namely,

$$ \begin{aligned} \sum_{j=x,y,z} \sum_{\nu=1}^{N_n} \underline{\underline{\underline{\underline{\underline{M}}}}}^{ij}_{\nu\nu'}(\mathbf{R}) \underline{\underline{\underline{\underline{T}}}}^{jk}_{\nu'\nu''} &= (M_{\nu}\delta_{ik} + \frac{m}{e}\underline{\underline{\underline{\underline{\underline{P}}}}}^{\nu}_{ik}(\mathbf{R}))\delta_{\nu''1} \\ &= M_{\text{tot}}[\delta_{\nu''1}\underline{\underline{\underline{\underline{\underline{J}}}}}^{ik}_{\nu\nu''}]^T, \end{aligned} \quad (\text{F8}) $$

where we used the sum rule of Eq. (18) in the first equality and Eq. (F3) in the second. We identify the term in square brackets in the last equality as $\underline{\underline{\underline{\underline{J}}}^{\text{c.m.}}}$. Further relations that will be used below are

$$ \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}\underline{\underline{\underline{\underline{T}}}} = \begin{pmatrix} \underline{\underline{\underline{\underline{I}}}}^{(3)} & \underline{\underline{\underline{\underline{\underline{0}}}}} \\ \underline{\underline{\underline{\underline{\underline{0}}}}} & \underline{\underline{\underline{\underline{\underline{0}}}}} \end{pmatrix}, \quad (\text{F9}) $$

$$ \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}}\underline{\underline{\underline{\underline{T}}}} = \underline{\underline{\underline{\underline{0}}}}, \quad (\text{F10}) $$

Equation (F5) is written by introducing the two components, c.m. and int, of the Jacobian as

$$ \begin{aligned} \underline{\underline{\underline{\underline{\underline{M}}}}}^{-1}(\mathbf{R}') = \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}]^T &+ \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}}]^T \\ &+ \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}]^T + \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{int}}]^T. \end{aligned} \quad (\text{F11}) $$

Using Eq. (F6), the first term on the right-hand side can be rewritten as

$$ \begin{aligned} \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}]^T &= \frac{1}{M_{\text{tot}}} \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}(\mathbf{R})\underline{\underline{\underline{\underline{\underline{T}}}}} \\ &= \frac{1}{M_{\text{tot}}} \underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}\underline{\underline{\underline{\underline{T}}}}, \end{aligned} \quad (\text{F12}) $$

and from Eq. (F9), we obtain

$$ \frac{1}{2} (-i\hbar\nabla')^T [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}} \underline{\underline{\underline{\underline{\underline{\underline{M}}}}}}^{-1}(\mathbf{R}) [\underline{\underline{\underline{\underline{\underline{J}}}}}^{\text{c.m.}}]^T] (-i\hbar\nabla') = \frac{\hat{P}_{\text{c.m.}}^2}{2M_{\text{tot}}}. \quad (\text{F13}) $$

A similar procedure, which uses Eq. (F10), is employed to show that the cross terms (second and third terms on the right-hand side) in Eq. (F11) do not contribute to the kinetic energy. Therefore, the final result reads

$$ \hat{H}_n = \frac{\hat{P}_{\text{c.m.}}^2}{2M_{\text{tot}}} + \frac{1}{2} (-i\hbar\nabla')^T \underline{\underline{\underline{\underline{\underline{M}}}}}^{-1}(\mathbf{R}') (-i\hbar\nabla') + E(\mathbf{R}'). \quad (\text{F14}) $$
# APPENDIX G: NUMERICAL DETAILS

## 1. Proton-transfer model

A model of a proton involved in a one-dimensional hydrogen bond like O-H-O is considered [45], with potential

$$V(r, R) = D[e^{-2a[(R/2)+r-d]} - 2e^{-a[(R/2)+r-d]} + 1] \\ + Dc^2[e^{-(2a/c)[(R/2)-r-d]} - 2e^{-(a/c)[(R/2)-r-d]} \\ + Ae^{-BR} - \frac{C}{R^6}. \quad (G1)$$

Here, *r* indicates the position of the proton measured from the center of the O-O bond, and *R* stands for the O-O distance. The chosen parameters of the Morse potential are *D*=60kcal/mol, *d*=0.95Å, *a*= 2.52 Å⁻¹; *c*= 0.707 makes the potential for the proton asymmetric, mimicking a strong O-H-O bond. The other parameters are *A*= 2.32 × 10⁵ kcal/mol, *B*= 3.15 Å⁻¹, and *C*= 2.31 × 10⁴ kcal/mol/Å⁶. The full Hamiltonian of the system involves *V(r, R)* and the kinetic energies of the oxygen atoms and of the proton, namely,

$$\hat{H}(r, R_{O^{-}}, R_{O^{+}}) = \sum_{\nu=+,-} \frac{-\hbar^2 \nabla_{O^{\nu}}^2}{2M_{O^{\nu}}} + \frac{-\hbar^2 \nabla_r^2}{2M_H} \\ + \hat{V}(r, R_{O^{-}}, R_{O^{+}}) \quad (G2) \\ = \sum_{\nu=+,-} \frac{-\hbar^2 \nabla_{O^{\nu}}^2}{2M_O} + \hat{H}_{BO}(r, R_{O^{-}}, R_{O^{+}}), \quad (G3)$$

where $\hat{V}$, according to Eq. (G1), depends only on the distance between the oxygen atoms, $R = |R_{O^{-}} - R_{O^{+}}|$.

In the static calculations, the adiabatic states have been computed by diagonalizing the BO Hamiltonian in Eq. (G3) on a 400 × 400 spatial grid. The eigenvalues of the full Hamiltonian in Eq. (G2) are determined using a Gaussian quadrature method with 20 points for *R*, the distance between the two heavy atoms, and 34 for *r*, the displacement of the proton from the c.m. of the heavy atoms. When the Hamiltonian with position-dependent dressed masses is used for computing the eigenvalues, *R* is again the distance between the two heavy atoms. In this case, as described in the text, the BO approximation has been introduced before separating the c.m. motion, and the eigenvalues of the Hamiltonian in internal coordinates [indicated by the prime symbols in Eq. (20)] have been computed.

In the dynamics, we use the three coordinates, i.e., $R_{O^{+}}$, $R_{O^{-}}$, and $r = r_{H}$, in order to test the conservation of the position of the c.m. The results in the paper are shown for

the same number of periods in all cases, using the following: the velocity-Verlet algorithm to integrate the classical equations, with a time step of 1 fs; the Crank-Nicolson [49] algorithm for the proton (quantum) equation in Ehrenfest, with a time step of 10⁻⁴ fs; and the Euler algorithm if the force depends on the velocity [see Eq. (G7)], with a time step of 0.0625 fs, where the stability of the integration has been tested based on the energy conservation. The position of the proton is estimated as the expectation value of the position operator on the proton wave function at the instantaneous O-O geometry.

## 2. Molecular geometries

The position-dependent mass corrections, *A* matrix, and the vibrational spectra of the molecular systems of Sec. IV B have been performed on geometries that have been optimized with the Gaussian electronic structure package [56] at the BLYP/aug-cc-pVTZ level. Bond lengths and angles for the molecules H₂O, H₃O⁺, NH₃, and CH₄ are listed in Table IV, whereas Cartesian coordinates of CH₃OH are given in Table V.

## 3. Calculation of the *A* matrix

We have computed the *A* matrix using density functional perturbation theory [25,28–30,50,60] as described in Ref. [23] and checked that the sum rule of Eq. (18) is satisfied. The numerical scheme has been implemented in the electronic structure package CPMD [51]. Calculations

<table><caption>TABLE IV. Optimized geometries of H₂O, H₃O⁺, NH₃, and CH₄ given in internal coordinates, i.e., bonds lengths and angles.</caption><thead><tr><th>Molecule</th><th>Coordinate</th><th>Value</th></tr></thead><tbody><tr><td rowspan="2">H<sub>2</sub>O</td><td>r<sub>OH</sub></td><td>97.2 pm</td></tr><tr><td>∠(HOH)</td><td>104.5°</td></tr><tr><td rowspan="2">H<sub>3</sub>O<sup>+</sup></td><td>r<sub>OH</sub></td><td>99.0 pm</td></tr><tr><td>∠(HOH)</td><td>112.0°</td></tr><tr><td rowspan="2">NH<sub>3</sub></td><td>r<sub>NH</sub></td><td>102.2 pm</td></tr><tr><td>∠(HNH)</td><td>106.7°</td></tr><tr><td rowspan="2">CH<sub>4</sub></td><td>r<sub>CH</sub></td><td>109.4 pm</td></tr><tr><td>∠(HCH)</td><td>109.5°</td></tr></tbody></table>

<table><caption>TABLE V. Optimized geometry of CH₃OH given in Cartesian coordinates.</caption><thead><tr><th>Molecule</th><th>Atom</th><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td rowspan="6">CH<sub>3</sub>OH</td><td>O</td><td>-0.6937</td><td>-0.0071</td><td>0.1741</td></tr><tr><td>C</td><td>0.7177</td><td>0.0086</td><td>-0.1244</td></tr><tr><td>H</td><td>1.2358</td><td>0.0367</td><td>0.8388</td></tr><tr><td>H</td><td>1.0055</td><td>0.8983</td><td>-0.7048</td></tr><tr><td>H</td><td>1.0335</td><td>-0.8950</td><td>-0.6675</td></tr><tr><td>H</td><td>-1.1832</td><td>-0.0322</td><td>-0.6644</td></tr></tbody></table>
have been performed using Troullier-Martins [52] pseudopotentials in the Becke-Lee-Yang-Parr [53,54] (BLYP) approximation of the exchange-correlation kernel. As stated above, the molecular geometry is the equilibrium geometry at the BLYP level, employing the aug-cc-pVTZ basis set [55] in the Gaussian electronic structure program [56].

Remember that the A matrix is a (3N_n × 3N_n) matrix. For the H_2 molecule shown in Table I, the blocks of the matrix are

$$ \left( \frac{(A_{H_1 H_1})^{ij}}{(A_{H_2 H_1})^{ij}} \middle| \frac{(A_{H_1 H_2})^{ij}}{(A_{H_2 H_2})^{ij}} \right) \qquad (G4) $$

and indices i, j running over the Cartesian components x, y, z. Each block is then a (3 × 3) matrix. The sum rule in Eq. (18) reads, in this case,

$$ \begin{aligned} \sum_{\nu, \nu'=1}^{N_n} A_{\nu\nu'}^{xx}(\mathbf{R}) &= [(A_{H_1 H_1})^{xx} + (A_{H_1 H_2})^{xx} \\ &\quad + (A_{H_2 H_1})^{xx} + (A_{H_2 H_2})^{xx}] \\ &= 1.998 \approx 2 \end{aligned} \qquad (G5) $$

and similarly for the other Cartesian components. This result is obtained by summing the entries of the matrix in Table I, and we indeed find the total electronic mass ($m_e = 1, N_{el} = 2$) of the system as expected from Eq. (18).

### 4. Normal mode analysis

It is easy to prove that, given a Lagrangian of the form

$$ \mathcal{L}(\dot{\mathbf{R}}, \mathbf{R}) = \frac{1}{2} \dot{\mathbf{R}}^T \underline{\underline{\mathcal{M}(\mathbf{R})\dot{\mathbf{R}}}} - E(\mathbf{R}), \qquad (G6) $$

<table><caption>TABLE VI. A matrix for H_3O⁺. The sum rule yields (9.9997, 9.9997, 9.9997) for i, j = xx, yy, zz, as 10 electrons are present in the system.</caption><thead><tr><th>A</th><th>Oxygen</th><th>Hydrogen 1</th><th>Hydrogen 2</th><th>Hydrogen 3</th></tr></thead><tbody><tr><td rowspan="3">Oxygen</td><td rowspan="2">7.7280</td><td>0.4013</td><td>0.1174</td><td>0.1174</td></tr><tr><td>7.7280</td><td>0.0227 0.0140</td><td>-0.1639 0.3067</td><td>-0.0121 0.1639</td></tr><tr><td>7.4688</td><td>0.1232 0.2546</td><td>-0.1067 -0.0616</td><td>0.2546 0.1067</td><td>-0.0616 0.2546</td></tr><tr><td rowspan="3">Hydrogen 1</td><td rowspan="3">0.2427</td><td></td><td>0.0073</td><td>0.0073</td></tr><tr><td>0.3900</td><td>-0.0692 -0.0353</td><td>0.0096 0.0297</td><td>0.0353 0.0096</td></tr><tr><td>0.2740</td><td>0.0162 0.0314</td><td>0.0302 -0.0162</td><td>0.0314 0.0302</td></tr><tr><td rowspan="3">Hydrogen 2</td><td rowspan="3"></td><td></td><td>0.3532</td><td>0.0638 0.0599</td><td>0.0107 0.0373</td></tr><tr><td>0.2795</td><td>0.0346 -0.0373</td><td>0.0062 -0.0017</td><td>0.0353 -0.0017</td></tr><tr><td>0.2740</td><td>0.0353 -0.0353</td><td>0.0302 0.0302</td><td>0.0353 -0.0638</td></tr><tr><td rowspan="3">Hydrogen 3</td><td rowspan="3"></td><td></td><td></td><td>0.3532</td><td>-0.0638 0.2795</td><td>-0.0599 0.0346</td></tr><tr><td>0.2740</td><td>0.0353 -0.0353</td><td>0.0302 0.0302</td><td>0.0353 -0.0638</td></tr><tr><td>0.2740</td><td>0.0353 -0.0353</td><td>0.0302 0.0302</td><td>0.0353 -0.0638</td></tr></tbody></table>

<table><caption>TABLE VII. A matrix for NH_3. The sum rule yields (9.9965, 9.9965, 9.9991) for i, j = xx, yy, zz, as 10 electrons are present in the system.</caption><thead><tr><th>A</th><th>Nitrogen</th><th>Hydrogen 1</th><th>Hydrogen 2</th><th>Hydrogen 3</th></tr></thead><tbody><tr><td rowspan="3">Nitrogen</td><td rowspan="2">6.8077</td><td>0.4125</td><td>0.0517 0.2084</td><td>0.0517 -0.2083</td></tr><tr><td>6.8077</td><td>-0.0685 -0.0125 0.2083 0.2922</td><td>-0.2083 0.2922</td><td>0.0062 0.0062</td></tr><tr><td>6.6408</td><td>0.2059 0.2826</td><td>0.1783 -0.1029 0.2826</td><td>-0.1783 -0.1029</td><td>0.2826</td></tr><tr><td rowspan="3">Hydrogen 1</td><td rowspan="3">0.4186</td><td></td><td>0.0121 -0.0663</td><td>0.0121 0.0663</td></tr><tr><td>0.9147 -0.2215</td><td>0.0175 0.0403</td><td>-0.0492 -0.0175</td><td>0.0403 -0.0492</td></tr><tr><td>0.4932</td><td>-0.0259 0.0535</td><td>0.0305 0.0259</td><td>0.0535 0.0305</td></tr><tr><td rowspan="3">Hydrogen 2</td><td rowspan="3"></td><td></td><td>0.7907 -0.2148</td><td>0.0543 -0.0419</td><td>-0.0593 -0.0419</td></tr><tr><td>0.5426</td><td>0.1107 0.0419</td><td>-0.0020 -0.0043</td><td>0.0419 -0.0020</td></tr><tr><td>0.4932</td><td>0.0593 -0.0043</td><td>0.0305 0.0305</td><td>0.0593 -0.0043</td></tr><tr><td rowspan="3">Hydrogen 3</td><td rowspan="3"></td><td></td><td></td><td>0.7907 0.2148</td><td>0.1919 0.5426</td></tr><tr><td>0.5426</td><td>0.1107 0.0419</td><td>-0.0020 -0.0043</td><td>0.1107 0.0419</td></tr><tr><td>0.4932</td><td>0.0593 -0.0043</td><td>0.0305 0.0305</td><td>0.4932 0.0305</td></tr></tbody></table>
<table><caption>TABLE VIII. A matrix for CH₄. The sum rule yields (9.9977, 9.9977, 9.9977) for i, j = xx, yy, zz, as 10 electrons are present in the system.</caption><thead><tr><th>A</th><th>Carbon</th><th>Hydrogen 1</th><th>Hydrogen 2</th><th>Hydrogen 3</th><th>Hydrogen 4</th></tr></thead><tbody><tr><td>Carbon</td><td>5.0406</td><td>0.2806<br>-0.0341</td><td>-0.1981<br>-0.0603</td><td>0.3757<br>-0.0677</td><td>0.0557<br>0.1207</td><td>0.3477<br>0.2238</td><td>0.1247<br>0.0453</td><td>0.0602<br>-0.1676</td><td>-0.0643<br>-0.0230</td><td>0.0230<br>0.3914</td><td>0.0822<br>0.0041</td></tr><tr><td></td><td>5.0407</td><td>-0.0341<br>-0.1981</td><td>0.3821<br>-0.0603</td><td>0.0420<br>0.0557</td><td>0.2238<br>0.2082</td><td>0.0602<br>-0.1675</td><td>0.3116<br>0.0821</td><td>0.0821<br>0.0041</td><td>0.0041<br>-0.0202</td><td>0.3777<br>-0.0081</td><td>0.0051<br>-0.0202</td></tr><tr><td>Hydrogen 1</td><td></td><td>0.6688</td><td>0.0352<br>0.5638</td><td>0.2046<br>0.0623</td><td>-0.0191<br>-0.0364</td><td>-0.0165<br>-0.0170</td><td>0.0019<br>-0.0695</td><td>-0.0269<br>0.0384</td><td>0.0114<br>0.0216</td><td>0.0112<br>0.0877</td><td>0.0183<br>-0.0081</td><td>0.0884<br>0.0051</td></tr><tr><td></td><td></td><td></td><td></td><td>0.0290<br>0.9154</td><td>-0.0260<br>0.0895</td><td>0.0134<br>0.0482</td><td>0.0134<br>0.0482</td><td>0.0587<br>0.0341</td><td>0.0011<br>0.0011</td><td>0.0011<br>-0.0202</td><td>0.0011<br>-0.0202</td><td></td></tr><tr><td>Hydrogen 2</td><td></td><td></td><td></td><td></td><td>0.5705<br>0.8340</td><td>0.0699<br>-0.2312</td><td>-0.0575<br>0.7436</td><td>-0.0163<br>0.0131</td><td>-0.0188<br>0.0061</td><td>0.0198<br>0.1093</td><td>0.0344<br>0.0385</td><td>0.0814<br>0.0053</td><td>0.0704<br>-0.0171</td></tr><tr><td>Hydrogen 3</td><td></td><td></td><td></td><td></td><td></td><td>0.5994<br>0.9120</td><td>-0.1288<br>0.6367</td><td>-0.0622<br>0.1732</td><td>0.0530<br>-0.0408</td><td>-0.0908<br>-0.0042</td><td>0.0030<br>0.0104</td><td>0.0439<br>0.0204</td><td>0.0004<br>0.0047</td></tr><tr><td>Hydrogen 4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1.0253<br>0.0237</td><td>0.0849<br>0.5543</td><td>-0.0293<br>-0.0043</td><td>0.0204<br>-0.0043</td><td>0.0047<br>0.5684</td><td></td></tr></tbody></table>

<table><caption>TABLE IX. A matrix for CH₃OH. The sum rule yields (18.0001, 18.0005, 18.0001) for i, j = xx, yy, zz, as 18 electrons are present in the system.</caption><thead><tr><th>A</th><th>Oxygen</th><th>Carbon</th><th>Hydrogen 1</th><th>Hydrogen 2</th><th>Hydrogen 3</th><th>Hydrogen 4</th></tr></thead><tbody><tr><td>Oxygen</td><td>8.7435</td><td>0.0085<br>8.5414</td><td>0.2570<br>-0.0002</td><td>-0.5225<br>-0.0064</td><td>0.1099<br>0.0490</td><td>0.0026<br>0.0953</td><td>0.0690<br>0.0042</td><td>0.1788<br>-0.0553</td><td>0.1805<br>-0.1800</td><td>-0.1810<br>-0.0085</td><td>0.1878<br>-0.1056</td><td>-0.1764<br>-0.1884</td><td>-0.1738<br>0.0070</td><td>0.2151<br>-0.0081</td><td>-0.0059<br>0.3434</td><td>-0.1885<br>0.0077</td></tr><tr><td></td><td></td><td></td><td></td><td>8.3358</td><td>0.1213</td><td>0.0028</td><td>0.0953</td><td>0.0861</td><td>0.0002</td><td>-0.1085<br>-0.1168</td><td>-0.0062<br>-0.0873</td><td>-0.1125<br>-0.1126</td><td>0.0066<br>-0.0879</td><td>-0.2949<br>-0.0094</td><td>0.0094<br>-0.0230</td><td>0.1128<br>-0.0230</td><td>0.1128<br>-0.0230</td></tr><tr><td>Carbon</td><td></td><td></td><td></td><td></td><td>4.9568</td><td>-0.0026<br>5.0359</td><td>-0.0064<br>-0.0018</td><td>0.1760<br>-0.0082</td><td>-0.1980<br>-0.0103</td><td>0.3162<br>-0.4347</td><td>-0.2255<br>-0.1000</td><td>0.1807<br>-0.1605</td><td>0.3036<br>0.2874</td><td>0.2413<br>0.1755</td><td>0.1797<br>-0.2720</td><td>-0.0297<br>-0.0003</td><td>-0.0038<br>-0.0038</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>4.9934</td><td>-0.3019<br>4.9934</td><td>-0.0116<br>-0.0018</td><td>0.1058<br>-0.0290</td><td>0.1579<br>0.0260</td><td>0.2472<br>0.0265</td><td>0.2655<br>0.0396</td><td>0.1585<br>-0.0365</td><td>0.2325<br>-0.0231</td><td>0.2872<br>-0.0231</td><td>0.0439<br>-0.0004</td><td>-0.0004<br>-0.0007</td></tr><tr><td>Hydrogen 1</td><td></td><td></td><td></td><td></td><td></td><td>0.7695<br>0.5711</td><td>0.0070<br>0.0115</td><td>0.1851<br>-0.0243</td><td>-0.0660<br>-0.0243</td><td>-0.0041<br>-0.0029</td><td>0.0115<br>-0.0041</td><td>-0.0669<br>-0.0029</td><td>0.0024<br>0.0236</td><td>0.0083<br>0.0047</td><td>0.1039<br>0.1019</td><td>0.0030<br>-0.0145</td><td>0.0540<br>-0.0145</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.9846<br>0.9846</td><td>0.0361<br>-0.0588</td><td>0.0361<br>-0.0588</td><td>0.0898<br>0.0332</td><td>0.0332<br>0.0634</td><td>0.0831<br>-0.0265</td><td>0.0265<br>-0.0015</td><td>0.0688<br>-0.0015</td><td>0.0050<br>-0.0230</td><td>0.0050<br>-0.0230</td></tr><tr><td>Hydrogen 2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.7709</td><td>0.1183<br>1.0678</td><td>-0.0810<br>-0.3065</td><td>-0.1124<br>0.0045</td><td>-0.0146<br>0.2071</td><td>-0.0292<br>-0.0140</td><td>-0.0292<br>-0.0140</td><td>0.0341<br>-0.0271</td><td>0.0050<br>-0.0271</td><td>0.0050<br>-0.0230</td></tr><tr><td>Hydrogen 3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1.0850<br>0.7787</td><td>0.0036<br>-0.1312</td><td>0.0036<br>-0.0852</td><td>0.0029<br>-0.0335</td><td>0.0029<br>-0.0581</td><td>0.0099<br>-0.0099</td><td>0.0224<br>-0.0166</td><td>0.0224<br>-0.0166</td><td>0.0224<br>-0.0166</td></tr><tr><td>Hydrogen 4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.5456<br>0.5456</td><td>0.0097<br>0.0097</td><td>0.3281<br>0.3281</td><td>0.0097<br>0.0097</td><td>0.0230<br>-0.0230</td><td>0.0230<br>-0.0230</td><td>0.0230<br>-0.0230</td><td>0.0230<br>-0.0230</td></tr></tbody></table>
the classical Hamiltonian of Eq. (D13) can be derived as its Legendre transform. Therefore, nuclear motion is classically governed by the Euler-Lagrange equation

$$ \underline{\underline{\underline{\underline{\mathcal{M}(\mathbf{R})\ddot{\mathbf{R}}}}} = -\underline{\underline{\nabla}}\underline{\underline{\mathcal{E}(\mathbf{R})}} - \frac{1}{2}\dot{\underline{\underline{\mathbf{R}}}}^T[\underline{\underline{\nabla}}\underline{\underline{\mathcal{M}(\mathbf{R})}}]\dot{\underline{\underline{\mathbf{R}}}}. \quad (G7) $$

This classical equation of motion is integrated using the Euler algorithm as described in Appendix G. If (i) we use internal coordinates, since the free motion of the c.m. can be separated as in Eq. (F14), (ii) we introduce the harmonic approximation of $E(\mathbf{R})$, and (iii) we neglect the velocity-dependent term, we obtain

$$ \ddot{\underline{\underline{\mathbf{R}}}} = -[\underline{\underline{\mathcal{M}^{-1}(\mathbf{R}_0)\underline{\underline{\mathbf{K}(\mathbf{R}_0)}}}]\underline{\underline{\mathbf{R}}}, \quad (G8) $$

with $\underline{\underline{\mathbf{K}}}$ the Hessian matrix computed from the ground-state electronic potential. The term in square brackets is evaluated at the equilibrium geometry $\mathbf{R}_0$. The diagonalization of the matrix in square brackets yields corrected $\nu + \Delta\nu$ frequencies, as $\Delta\nu$ includes the effect of electrons that follow the motion of the nuclei nonadiabatically, namely, not instantaneously.

## APPENDIX H: A MATRIX FOR $H_3O^+, NH_3, CH_4, CH_3OH$

At the equilibrium geometries reported in Appendix G 2, we have determined the A matrix for $H_3O^+, NH_3, CH_4$, and $CH_3OH$. The matrices are reported in the following tables. Table VI shows the results for $H_3O^+$, where the sum rule of Eq. (17) is verified in the caption of the table. Tables VII, VIII, and IX similarly report the A-matrices for $NH_3, CH_4$, and $CH_3OH$, respectively. The sum rules are verified in the captions as well.

The sum rule has been used to determine the number of digits of the reported values of the elements of the position-dependent mass correction matrix: We have dropped all digits for which the sum rule numerically deviates from the expected value.

[1] M. Born and R. J. Oppenheimer, Zur Quantentheorie der Molekeln, Ann. Phys. (Berlin) 389, 457 (1927) [On the Quantum Theory of Molecules, Ann. Phys. (Amsterdam) 84, 458 (1927)].

[2] P. R. Bunker and R. E. Moss, The Breakdown of the Born-Oppenheimer Approximation: The Effective Vibration-Rotation Hamiltonian for a Diatomic Molecule, Mol. Phys. 33, 417 (1977).

[3] P. R. Bunker, C. J. McLarnon, and R. E. Moss, Application of the Effective Vibration-Rotation Hamiltonian to $H_2$ and $D_2$, Mol. Phys. 33, 425 (1977).

[4] D. W. Schwenke, Beyond the Potential Energy Surface: Ab Initio Corrections to the Born-Oppenheimer Approximation for $H_2O$, J. Phys. Chem. A 105, 2352 (2001).

[5] O. L. Polyansky and J. Tennyson, Ab Initio Calculation of the Rotation-Vibration Energy Levels of $H_3^+$ and Its Isotopomers to Spectroscopic Accuracy, J. Chem. Phys. 110, 5056 (1999).

[6] A. Owens, S. N. Yurchenko, A. Yachmenev, J. Tennyson, and W. Thiel, Accurate Ab Initio Vibrational Energies of Methyl Chloride, J. Chem. Phys. 142, 244306 (2015).

[7] M. Pavanello, L. Adamowicz, A. Alijah, N. F. Zobov, I. I. Mizus, O. L. Polyansky, J. Tennyson, T. Szidarovszky, A. G. Császár, M. Berg, A. Petrignani, and A. Wolf, Precision Measurements and Computations of Transition Energies in Rotationally Cold Triatomic Hydrogen Ions up to the Midvisible Spectral Range, Phys. Rev. Lett. 108, 023002 (2012).

[8] K.-Y. Wu, Y.-H. Lien, C.-C. Liao, Y.-R. Lin, and J.-T. Shy, Measurement of the $\nu_2$ Fundamental Band of $H_3^+$, Phys. Rev. A 88, 032507 (2013).

[9] L. Wallace, P. Bernath, W. Livingston, K. Hinkle, J. Busler, B. Guo, and K. Zhang, Water on the Sun, Science 268, 1155 (1995).

[10] N. C. Handy and A. M. Lee, The Adiabatic Approximation, Chem. Phys. Lett. 252, 425 (1996).

[11] W. Kutzelnigg, The Adiabatic Approximation I. The Physical Background of the Born-Handy Ansatz, Mol. Phys. 90, 909 (1997).

[12] W. Kutzelnigg, Which Masses are Vibrating or Rotating in a Molecule?, Mol. Phys. 105, 2627 (2007).

[13] A. S. Goldhaber, Newtonian Adiabatics Unified, Phys. Rev. A 71, 062102 (2005).

[14] M. Pavanello, L. Adamowicz, A. Alijah, N. F. Zobov, I. I. Mizus, O. L. Polyansky, J. Tennyson, T. Szidarovszky, and A. G. Császár, Calibration-Quality Adiabatic Potential Energy Surfaces for $H_3^+$ and Its Isotopologues, J. Chem. Phys. 136, 184303 (2012).

[15] O. L. Polyansky, N. F. Zobov, S. Viti, J. Tennyson, P. F. Bernath, and L. Wallace, Water on the Sun: Line Assignments Based on Variational Calculations, Science 277, 346 (1997).

[16] P. R. Bunker and R. E. Moss, The Effect of the Breakdown of the Born-Oppenheimer Approximation on the Rotation-Vibration Hamiltonian of a Triatomic Molecule, J. Mol. Spectrosc. 80, 217 (1980).

[17] L. Wolniewicz, Nonadiabatic Energies of the Ground State of the Hydrogen Molecule, J. Chem. Phys. 103, 1792 (1995).

[18] K. Pachucki and J. Komasa, Nonadiabatic Corrections to Rovibrational Levels of $H_2$, J. Chem. Phys. 130, 164113 (2009).

[19] T. Furtenbacher, T. Szidarovszky, E. Mátyus, C. Fábri, and A. G. Császár, Analysis of the Rotational-Vibrational States of the Molecular Ion $H_3^+$, J. Chem. Theory Comput. 9, 5471 (2013).

[20] E. Mátyus, T. Szidarovszky, and A. G. Császár, Modelling Non-adiabatic Effects in $H_3^+$: Solution of the Rovibrational Schrödinger Equation with Motion-Dependent Masses and Mass Surfaces, J. Chem. Phys. 141, 154111 (2014).

[21] L. G. Diniz, J. R. Mohallem, A. Alijah, M. Pavanello, L. Adamowicz, O. L. Polyansky, and J. Tennyson, Vibrationally and Rotationally Nonadiabatic Calculations on $H_3^+$ Using Coordinate-Dependent Vibrational and Rotational Masses, Phys. Rev. A 88, 032506 (2013).
[22] C. Fábri, E. Mátyus, and A. G. Császár, *Numerically Constructed Internal-Coordinate Hamiltonian with Eckart Embedding and Its Application for the Inversion Tunneling of Ammonia*, Spectrochim. Acta, Part A **119**, 84 (2014).

[23] A. Scherrer, F. Agostini, D. Sebastiani, E. K. U. Gross, and R. Vuilleumier, *Nuclear Velocity Perturbation Theory for Vibrational Circular Dichroism: An Approach Based on the Exact Factorization of the Electron-Nuclear Wave Function*, J. Chem. Phys. **143**, 074106 (2015).

[24] A. Abedi, N. T. Maitra, and E. K. U. Gross, *Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function*, Phys. Rev. Lett. **105**, 123002 (2010).

[25] A. Scherrer, R. Vuilleumier, and D. Sebastiani, *Nuclear Velocity Perturbation Theory of Vibrational Circular Dichroism*, J. Chem. Theory Comput. **9**, 5305 (2013).

[26] L. A. Nafie, *Adiabatic Molecular Properties Beyond the Born-Oppenheimer Approximation. Complete Adiabatic Wave Functions and Vibrationally Induced Electronic Current Density*, J. Chem. Phys. **79**, 4950 (1983).

[27] F. G. Eich and F. Agostini, *The Adiabatic Limit of the Exact Factorization of the Electron-Nuclear Wave Function*, J. Chem. Phys. **145**, 054110 (2016).

[28] S. Baroni, P. Giannozzi, and A. Testa, *Green's-Function Approach to Linear Response in Solids*, Phys. Rev. Lett. **58**, 1861 (1987).

[29] S. Baroni, S. de Gironcoli, A. dal Corso, and P. Giannozzi, *Phonons and Related Crystal Properties from Density-Functional Perturbation Theory*, Rev. Mod. Phys. **73**, 515 (2001).

[30] X. Gonze, *Adiabatic Density-Functional Perturbation Theory*, Phys. Rev. A **52**, 1096 (1995).

[31] A. Schild, F. Agostini, and E. K. U. Gross, *A Study of the Electronic Flux Density in the Born-Oppenheimer Limit*, J. Phys. Chem. A **120**, 3316 (2016).

[32] A. Abedi, F. Agostini, Y. Suzuki, and E. K. U. Gross, *Dynamical Steps that Bridge Piecewise Adiabatic Shapes in the Exact Time-Dependent Potential Energy Surface*, Phys. Rev. Lett. **110**, 263001 (2013).

[33] F. Agostini, S. K. Min, and E. K. U. Gross, *Semiclassical Analysis of the Electron-Nuclear Coupling in Electronic Non-adiabatic Processes*, Ann. Phys. (Amsterdam) **527**, 546 (2015).

[34] J. Frenkel, *Wave Mechanics* (Clarendon, Oxford, 1934).

[35] J. L. Alonso, J. Clemente-Gallardo, P. Echeniche-Robba, and J. A. Jover-Galtier, *Comment on Correlated Electron-Nuclear Dynamics: Exact Factorization of the Molecular Wave-Function*, J. Chem. Phys. **139**, 087101 (2013).

[36] A. Abedi, N. T. Maitra, and E. K. U. Gross, *Reply to Comment on Correlated Electron-Nuclear Dynamics: Exact Factorization of the Molecular Wave-Function*, J. Chem. Phys. **139**, 087102 (2013).

[37] J. C. Tully, *Perspective on, Zur Quantentheorie der Molekeln*, Theor. Chem. Acc. **103**, 173 (2000).

[38] G. A. Hagedorn, *High Order Corrections to the Time-Dependent Born-Oppenheimer Approximation I: Smooth Potentials*, Ann. Math. **124**, 571 (1986); **126**, 219(E) (1987).

[39] G. A. Hagedorn and A. Joye, *Mathematical Analysis of Born-Oppenheimer Approximations*, in *Spectral Theory and Mathematical Physics: A Festschrift in Honor of Barry Simon's 60th Birthday, Part 1: Quantum Field Theory, Statistical Mechanics, and Nonrelativistic Quantum Systems* (American Mathematical Society, Providence, RI, 2007), Vol. 76.1, pp. 203–226.

[40] G. Rigolin, G. Ortiz, and V. H. Ponce, *Beyond the Quantum Adiabatic Approximation: Adiabatic Perturbation Theory*, Phys. Rev. A **78**, 052508 (2008).

[41] R. Requist and O. Pankratov, *Adiabatic Approximation in Time-Dependent Reduced-Density-Matrix Functional Theory*, Phys. Rev. A **81**, 042519 (2010).

[42] D. Yang and A. Rauk, *Sum Rules for Atomic Polar and Axial Tensors from Vibronic Coupling Theory*, Chem. Phys. **178**, 147 (1993).

[43] W. B. Person and J. H. Newton, *Dipole Moment Derivatives and Infrared Intensities. I. Polar Tensors*, J. Chem. Phys. **61**, 1040 (1974).

[44] P. J. Stephens, K. J. Jalkanen, R. D. Amos, P. Lazzeretti, and R. Zanasi, *Ab Initio Calculations of Atomic Polar and Axial Tensors for Hydrogen Fluoride, Water, Ammonia, and Methane*, J. Phys. Chem. **94**, 1811 (1990).

[45] D. C. Marinica, M.-P. Gaigeot, and D. Borgis, *Generating Approximate Wigner Distributions Using Gaussian Phase Packets Propagation in Imaginary Time*, Chem. Phys. Lett. **423**, 390 (2006).

[46] G. Hanna and R. Kapral, *Quantum-Classical Liouville Dynamics of Nonadiabatic Proton Transfer*, J. Chem. Phys. **122**, 244505 (2005).

[47] D. J. D. Wilson, C. E. Mohn, and T. Helgaker, *The Rotational Tensor as a Benchmark for Density-Functional Theory Calculations of Molecular Magnetic Properties*, J. Chem. Theory Comput. **1**, 877 (2005).

[48] D. W. Schwenke, *Towards Accurate Ab Initio Predictions of the Vibrational Spectrum of Methane*, Spectrochim. Acta, Part A **58**, 849 (2002).

[49] J. Crank and P. Nicolson, *A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat Conduction Type*, Proc. Cambridge Philos. Soc. **43**, 50 (1947).

[50] A. Putrino, D. Sebastiani, and M. Parrinello, *Generalized Variational Density Functional Perturbation Theory*, J. Chem. Phys. **113**, 7102 (2000).

[51] CPMD, http://www.cpmd.org, for CPMD-4.1. Copyright 2000-2017 jointly IBM Corp. and Max Planck Institute, Stuttgart.

[52] N. Troullier and J. L. Martins, *Efficient Pseudopotentials for Plane-Wave Calculations*, Phys. Rev. B **43**, 1993 (1991).

[53] A. D. Becke, *Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior*, Phys. Rev. A **38**, 3098 (1988).

[54] C. Lee, W. Yang, and R. G. Parr, *Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density*, Phys. Rev. B **37**, 785 (1988).

[55] R. A. Kendall, T. H. Dunning, and R. J. Harrison, *Electron Affinities of the First Row Atoms Revisited. Systematic Basis Sets and Wave Functions*, J. Chem. Phys. **96**, 6796 (1992).

[56] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson et al., *Gaussian 09 Revision D.01* (Gaussian Inc., Wallingford, CT, 2009).
[57] R. E. Moss, *On the Adiabatic and Non-adiabatic Corrections in the Ground Electronic State of the Hydrogen Molecular Cation*, Mol. Phys. **89**, 195 (1996).

[58] S. Fatehi and J. E. Subotnik, *Derivative Couplings with Built-in Electron-Translation Factors: Application to Benzene*, J. Phys. Chem. Lett. **3**, 2039 (2012).

[59] F. Agostini, A. Abedi, and E. K. U. Gross, *Classical Nuclear Motion Coupled to Electronic Non-adiabatic Transitions*, J. Chem. Phys. **141**, 214101 (2014).

[60] T. Watermann, A. Scherrer, and D. Sebastiani, *Linear Response Methods in Quantum Chemistry* (Springer International Publishing, New York, 2014), p. 97.