Non-adiabatic mass correction to the rovibrational states of molecules. Numerical application for the H + 2 molecular ion
aa r X i v : . [ phy s i c s . c h e m - ph ] N ov Non-adiabatic mass correction to the rovibrational states ofmolecules. Numerical application for the H +2 molecular ion Edit M´atyus Institute of Chemistry, E¨otv¨os Lor´and University,P´azm´any P´eter s´et´any 1/A, Budapest, H-1117, Hungary (Dated: November 6, 2018)
Abstract
General transformation expressions of the second-order non-adiabatic Hamiltonian of the atomicnuclei, including the kinetic-energy correction terms, are derived upon the change from laboratory-fixed Cartesian coordinates to general curvilinear coordinate systems commonly used in rovibra-tional computations. The kinetic-energy or so-called “mass-correction” tensor elements are com-puted with the stochastic variational method and floating explicitly correlated Gaussian functionsfor the H +2 molecular ion in its ground electronic state. (Further numerical applications for the He +2 molecular ion are presented in the forthcoming paper, Paper II.) The general, curvilinearnon-adiabatic kinetic energy operator expressions are used in the examples and non-adiabatic rovi-brational energies and corrections are determined by solving the rovibrational Schr¨odinger equationincluding the diagonal Born–Oppenheimer as well as the mass-tensor corrections. . INTRODUCTION The mass (or better: kinetic-energy) correction terms in the second-order non-adiabaticrovibrational Hamiltonian has been re-discovered in several very different contexts since theadvent of quantum mechanics. In contrast to the diagonal Born–Oppenheimer (DBOC) cor-rection, the mass-correction terms have been rarely included in rovibrational computations,apart from the common arguments of using atomic instead of nuclear masses. The practiceof using atomic masses is confirmed by a posteriori by obtaining a better agreement withthe experimentally observed rovibrational transitions (with or without accounting for anyrelativistic and radiative effects in the computations).The earliest account of the kinetic-energy correction terms which we mention here isfrom 1964 by Fisk and Kirtman [1] who wrote about a velocity-dependent effective potentialenergy surface (PES), which they obtained obtained using van Vleck’s perturbation theoryapproach. Bunker and Moss derived effective non-adiabatic rovibrational Hamiltonians us-ing van Vleck’s perturbation theory for di- [2] and triatomic [3] molecules and used theirformalism for diatomics [4]. Later Schwenke elaborated on these correction terms and com-puted them in curvilinear coordinates for diatomic molecules [5] as well as for the triatomicH O molecule [6].Pachucki and Komasa introduced non-adiabatic perturbation theory [7] and arrived atthe same mass-correction functions for diatomic molecules, which they successfully used for aseries of systems [8, 9], including the H molecule in its ground electronic state [8]. Kutzelnigg[10] and Jaquet and Kutzelnigg [11] derived and computed mass-correction functions bystarting out from a careful consideration of the separation of the center of mass and the totalmass of the molecule in an adiabatic theory. There were also more empirical proposals forthe “rotational” and “vibrational” mass correction functions of di- and triatomic molecules,which connected the mass correction to the electron density assignable to the nuclei [12].In an elegant and formal series of work, Teufel and co-workers introduced adiabatic pertur-bation theory in quantum dynamics and used it for the Born–Oppenheimer (BO) separationof the electronic and nuclear degrees of freedom, the coupling of which can be characterizedby the square root of the electron-to-nucleus mass ratio. The authors identified the almostinvariant subspace for the electronic problem and using this subspace they derived an ef-fective Hamiltonian for the quantum mechanical motion of the atomic nuclei corresponding2o an isolated electronic state (separated with gaps from all other electronic states). Ex-panding this effective Hamiltonian in terms of increasing orders of the coupling parameter,corrections are obtained to the (zeroth order) BO Hamiltonian of the atomic nuclei. Thesecond-order Hamiltonian contains the DBOC and mass-correction tensor obtained also inother perturbative procedures. Teufel and co-workers mention the derivation of Weigert andLittlejohn from 1993 [13] using Weyl calculus for formally diagonalizing the multicomponentwave equation, an approach which resulted a similar effective Hamiltonian for the nuclearmotion.Most recently, independent of earlier work, the mass correction tensor [14] was derivedfrom exact factorization through nuclear velocity perturbation theory [15], which successfullyprovides a theoretical framework also for vibrational circular dichroism.The relation of the mass-correction tensor to the computation of magnetic properties hasbeen observed (for diatomics) also by Bunker and Moss [2]. Herman and Asgharian [16] andHerman and Ogilvie [17] pointed out its connection to the electronic contribution to therotational and vibrational (paramagnetic) g -factors (see also Ref. [18]). For example, the g -factors computed using a linear response method with full-CI/aug-cc-pVTZ by Sauer etal. [19] for HeH + were found to be in a few % agreement with computations of Pachucki andKomasa [9] using the optimized explicitly correlated Gaussian functions (ECGs). Ogilvieand co-workers [20] (as well as Pachuchki and Komasa [8]) computed the vibrational androtational g -factors and adiabatic corrections for the hydrogen molecule. The vibrational g -factor was computed for the bending mode of HCN in Ref. [21]Obviously, there has been a substantial theoretical and computational progress over thepast half a century concerning the theory and computational applications of the complete second-order rovibrational Hamiltonian, which includes the mass-correction tensor. We findit fascinating that the same quantity, i.e., the same mass-correction term, appear in a va-riety of contexts and essentially the same quantity has been obtained starting from variousdirections and by using very different (perturbational) formalisms. It would be interestingto explore the various aspects and theoretical connections between the different derivations.After recognizing all these earlier, somewhat parallel developments, it is interesting tonotice that the mass-correction tensor is not routinely included in the nowadays numerically“exact” rovibrational computations, in which however the PES (almost routinely) includesthe diagonal Born–Oppenheimer correction, which is also a second-order term in the non-3diabatic Hamiltonian (see for example, Ref. [10, 22]). The accuracy of the present-dayhigh-resolution (precision) spectroscopy measurements (see for example Ref. [23] discussedin Paper II [24]) implies that one has to go beyond the empirical non-adiabatic correctionsin which the nuclear masses are arbitrarily replaced with some effective (constant, usuallythe atomic) masses (which are thought to account for some of the mass-correction tensoreffects). In particular, there are at least two important “families” of small effects to be ac-counted for in the computations: “non-adiabatic corrections” and “relativistic and radiativecorrections”—both represent challenges for a rigorous theoretical description.We wonder why the computation and the use of the mass-correction tensor did notbecome routine in rovibrational studies. We think and agree with the authors of Ref. [14]that for a widespread applicability it would be important to compute the mass-correctiontensor in Cartesian coordinates (Schwenke in Ref. [6] also mentions this direction as possiblefurther development for his curvilinear derivation). In particular, the widespread and generalcomputation of the DBOC has become possible by its formulation and computation in simpleCartesian coordinates following Handy and co-workers [25–27] (instead of using the system-and coordinate-dependent form). The laboratory-frame Cartesian coordinate expression ofthe DBOC was later confirmed in a stringent numerical test by Cencek and Kutzelnigg [28]and was explained in formal terms by Kutzelnigg [29]. The DBOC is a scalar quantity, whilethe mass-correction factor is a tensor. As it was also discussed by Schwenke [6], the mass-correction tensor is an inherently more complicated mathematical object. Nevertheless, itcan be computed using Cartesian coordinates with a selected frame for the nuclei from theelectronic energies and wave-functions as it was demonstrated in Ref. [14] at the equilibriumstructures of the H , H O, and CH OH molecules.In what follows, we derive the rovibrational Hamiltonian including the mass-correctionterm starting from laboratory-frame Cartesian coordinates to a general curvilinear coordi-nate system, in the spirit of the numerical-kinetic energy operator approach used in theGENIUSH protocol [30]. For the electronic structure computations, (due to the lack of anywidely available electronic structure method) we employ our in-house developed computerprogram, QUANTEN (QUANTum mechanical treatment of Electrons and atomic Nuclei),which uses the variational method and explicitly correlated Gaussian functions to solvethe Schr¨odinger equation. If all charges belong to the quantum system we solve a pre-Born–Oppenheimer problem [31–33]. For the present work, we extended QUANTEN for the case of4xed external charges to solve the electronic Schr¨odinger equation. In the second part of thepaper, we explain how the necessary adiabatic and non-adiabatic quantities are computed.Applications for a variety of poly-atomic and poly-electronic molecules will become pos-sible when the Cartesian mass-correction tensor can be computed with an efficient electronicstructure package over a broad range of molecular configurations. The implementation ofthe mass-correction tensor computed in Cartesian coordinates in a general curvilinear rovi-brational program ( e.g., the GENIUSH program) should be straightforward based on theexpressions derived in the first part of the article. To demonstrate the applicability of theexpressions and to show numerical examples we make the calculations explicit and givenumerical results for the homonuclear diatomic H +2 molecular ion and further applicationsfollow for He +2 in Paper II [24]. 5 I. THE SECOND-ORDER NON-ADIABATIC HAMILTONIAN
In this section, we re-iterate the second-order non-adiabatic Hamiltonian derived by Panati,Spohn and Teufel [22] and adjust some of the notation to that used in the derivation from ex-act factorization [14]. Both derivations have been carried out in laboratory-frame Cartesiancoordinates, and these general, N -atomic expressions provide the most convenient start-ing point for our work. The second-order non-adiabatic Hamiltonian of the atomic nucleicorresponding to a single electronic state (in Hartree atomic units ~ = m e = a = 1) is:ˆ H (2) = − N X i =1 X a m n ∂ R ia + N X i,j =1 X a,b m n ∂ R ia (cid:20) m n A ia,jb (cid:21) ∂ R jb + V + N X i =1 X a m n U ia (1)where ∂ R ia = ∂/∂ R ia is short notation for the partial derivative with respect to the laboratory-frame Cartesian coordinates, R ia ( a = X, Y, Z ). The first and thrid term is the BO kineticenergy operator and the potential energy surface (PES), respectively. The latter equals tothe electronic energy, E (el) , the eigenvalue of the electronic Schr¨odinger equation:ˆ H (el) ψ = E (el) ψ (2)with ˆ H (el) = − N el X i =1 ∆ x i + N el X i =1 N el X j>i | x i − x j | − N el X i =1 N X k =1 Z k | x i − R k | + N X k =1 N X l>k Z k Z l | R k − R l | (3)for the nuclear configuration R k ( k = 1 , , . . . , N ) ( Z k ( k = 1 , , . . . , N ) label the electriccharge of the nuclei).The correction terms, multiplied by the second power of p /m n are U ia = 12 D ∂ R ia ψ | (1 − ˆ P ) ∂ R ia ψ E el , (4)which gives rise to the well-known diagonal Born–Oppenheimer correction (DBOC) to thePES, U = P Ni =1 P a U ia /m n , while the correction tensor to the kinetic energy is A ia,jb ( R ) = 2 D ∂ R ia ψ | ( ˆ H (el) − E (el) ) − (1 − ˆ P ) ∂ R jb ψ E el , (5)6here ˆ P = | ψ ih ψ | is a projector to the ψ electronic eigenstate. We usually consider theground electronic state, but in principle E (el) and ψ can correspond to any isolated, elec-tronically excited state which is separated with a gap from rest of the electronic states[22].Note that we have introduced a factor of 2 in the expression of A to synchronize thenotation with Ref. [14]. Panati, Spohn, and Teufel derive the formalism for identical nuclearmasses, m n . The equations can be generalized to different nuclear masses by assuming that R ia are mass-scaled coordinates, i.e., R ia = p m i /m n R ′ ia .For later use, we define the second-order non-adiabatic kinetic energy operator as thesum of the terms containing the differential operators in ˆ H (2) asˆ T (2) = − N X i =1 X a m n ∂ R ia + N X i,j =1 X a,b m n ∂ R ia (cid:20) m n A ia,jb (cid:21) ∂ R jb = − N X i,j =1 X a,b m n ∂ R ia (cid:20) δ ia,jb − m n A ia,jb (cid:21) ∂ R jb = − N X i,j =1 X a,b m n ∂ R ia M ia,jb ∂ R jb , (6)where the elements of the effective mass tensor have been defined as M ia,jb = δ ia,jb − m n A ia,jb , (7)and we refer to A ia,jb as (elements of) the mass-correction tensor.For later convenience, a “condensed-index” labeling is introduced for the vector and tensorquantities as: R ia → R I and A ia,jb → A IJ and M ia,jb → M IJ (8)with ( ia ) → I = 3( i −
1) + a , ( jb ) → J = 3( j −
1) + b , and a, b ∈ { X ) , Y ) , Z ) } . Theexpanded ( ia ) and the condensed ( I ) indices will be used in an interchangeable manner.7 II. COORDINATES AND TRANSFORMATIONS TO DESCRIBE THE QUAN-TUM MECHANICAL MOTION OF THE ATOMIC NUCLEI
We may write the BO kinetic energy operator of the atomic nuclei in a compact form withconstant masses as ˆ T (0) = − m n div grad , (9)while the second-order non-adiabatic kinetic energy operator in a similarly compact form isˆ T (2) = − m n div M grad , (10)where M is the coordinate-dependent matrix defined in Eq. (7). At this starting point, allquantities are in laboratory-frame Cartesian coordinates (LFCC).Rovibrational computations can be efficiently performed (see for example Refs. [34, 35]),if the laboratory-fixed Cartesian coordinates are replaced with a physically motivated (curvi-linear) coordinate set, ξ . This physically motivated set includes a set of internal coordinatesthat are well suited to describe the internal motions (vibrations), three orientation angleswhich describe the orientation of the body-fixed frame with respect to the laboratory frame(rotations), and three coordinates which describe the translation of the center of mass (trans-lations). In the BO framework, Eq. (9), the mass corresponding to the overall translation ofthe nuclei is the sum of the nuclear masses. The effective mass corresponding to the overalltranslation in the non-adiabatic Hamiltonian, Eq. (10), will be discussed in Section V. A. General curvilinear coordinates
In a BO (or at least, constant-mass computations), one has to re-write ˆ T (0) according tothe R , R , . . . , R N ⇒ ξ , ξ , . . . , ξ N coordinate transformation:ˆ T (0) = − m n div grad ⇒ ˆ T (0) ξ = − m n div ξ grad ξ (11)8here div ξ and grad ξ are the divergence and the gradient in the new coordinates (no sub-script means plain laboratory-fixed Cartesian coordinates). General expressions for the curvi-linear form, ˆ T (0) ξ are routinely used in variational rovibrational computations [30, 36, 37].Concerning the transformation of ˆ T (2) to curvilinear coordinates, one has to consider thetransformationˆ T (2) = − m n div M grad ⇒ ˆ T (2) ξ = − m n div ξ M ( ξ ) grad ξ (12)together with the transformation, M → M ( ξ ) .In all earlier rovibrational computations which included mass-correction terms, tailor-made ˆ T (2) non-adiabatic kinetic energy operators have been derived corresponding to specific ξ choices ( i.e., diatomic molecules, or triatomic Radau coordinates [6]), and the resultingcoordinate-dependent mass coefficients have been computed from electronic-structure theory.Our aim in the present paper is to derive a general curvilinear expression for the non-adiabatic kinetic energy operator starting from the laboratory-frame expressions. It will pro-vide us with general formulae not only for the transformation of the differential operators butalso for the transformation of the mass tensor. Having all these transformation expressionsat hand the mass-correction tensor computed in plain Cartesian coordinates by electronicstructure theory will be straightforwardly applicable in rovibrational computations. Hence,we would arrive in some sence to a generalization of Handy’s method of DBOC for themass-correction tensor. We may say that we arrived at a general expression in curvilinearcoordinates if all operators are expressed in terms of the metric tensor and/or the Jacobi ten-sor, which are fundamental mathematical objects for the R , R , . . . , R N ⇒ ξ , ξ , . . . , ξ N coordinate change (and are routinely evaluated in the GENIUSH program over a grid toconstruct the curvilinear kinetic-energy operator terms during the computation).In curvilinear coordinates, ξ n ( n = 1 , , . . . , D , here D = 3 N ) with the covariant metrictensor g µν = δ ij ∂R i ∂ξ µ ∂R j ∂ξ ν , the divergence of an F vector field (using Einstein’s summationnotation) is div ξ F = ˜ g − / ∂ µ ˜ g / F µ (13)9here ∂ µ = ∂/∂ξ µ and ˜ g = det g µν . The gradient of a function φ in the new coordinates isgrad ξ φ = ∂ µ φ = g µν ∂ ν φ, (14)which includes the g µν contravariant metric tensor which is the inverse of the covariantmetric tensor, g αβ g βγ = δ γα . a. BO kinetic energy operator in curvilinear coordinates Using this notation, we canre-write the differential operator in the kinetic energy operator with constant masses asdiv grad φ = ˜ g − / ∂ µ ˜ g / g µν ∂ ν φ (15)with the normalization condition Z φ ∗ φ ˜ g / d ξ . . . d ξ D = 1 , (16)and thus the corresponding volume element isd V = ˜ g / d ξ . . . d ξ D . (17)Following Podolsky’s work [38] we introduce wave functions normalized according to Z χ ∗ χ d ξ . . . d ξ D , (18)and thus we can re-write the Schr¨odinger equation and the kinetic energy operator into amore symmetric form by first inserting φ = ˜ g − / χ and then multiplying with ˜ g / from theleft: 2( E − V ) φ = ˜ g − / ∂ µ ˜ g / g µν ∂ ν φ (19)2( E − V ) ˜ g − / χ = ˜ g − / ∂ µ ˜ g / g µν ∂ ν ˜ g − / χ (20)2( E − V ) χ = ˜ g − / ∂ µ ˜ g / g µν ∂ ν ˜ g − / χ. (21)10he last equation, Eq. (21), has been referred to as the Podolsky form of the kinetic energyoperator in curvilinear coordinates:ˆ T (BO)Pod = 12 m n ˜ g − / ∂ µ ˜ g / g µν ∂ ν g − / , (22)which is the form of the kinetic energy of the atomic nuclei implemented in the GENIUSHprogram [30, 39] and used in several rovibrational computations, e.g., in Refs. [40, 41]. b. Second-order non-adiabatic kinetic energy operator in curvilinear coordinates In asimilar spirit, we re-write the “div M grad” operator to curvilinear coordinates, ξ µ ( µ =1 , , . . . , D ), with the g µν metric and J iα Jacobian tensors as:div M grad φ = ˜ g − / ∂ µ ˜ g / M µν ∂ ν φ = ˜ g − / ∂ µ ˜ g / g µα M αβ g βν ∂ ν φ = ˜ g − / ∂ µ ˜ g / g µα J iα M ij J jβ g βν ∂ ν φ. (23)Note that the M αβ element of the M tensor corresponds to the α th, β th components incurvilinear coordinates, which is related to the i th, j th Cartesian components by M αβ = ∂R i ∂ξ α M ij ∂R j ∂ξ β , (24)and the J Jacobian tensor collects the derivatives of the Cartesian coordinates with respectto the curvilinear coordinates: J iα = ∂R i ∂ξ α .For later convenience, we introduce the short notationdiv M grad φ = ˜ g − / ∂ µ ˜ g / ˜ G µν ∂ ν φ. (25)with the effective ˜ G µν matrix which includes not only the mass-weighted metric tensor butalso the second-order non-adiabatic corrections to the kinetic energy operator (compare withEqs. (5) and (7) and note that the condensed-index labeling defined in Eq. (8) is used in11his section): ˜ G µν = g µα J iα M ij J jβ g βν . (26)It is important to recognize that elements of the mass-correction tensor are computed fromthe electronic wave function with the atomic nuclei positioned in a certain body-fixed (BF)frame [6]. We denote the mass matrix (mass-correction matrix, see Eq. (7)) correspondingto this BF frame with ¯ M ij ( ¯ A ij ). In order to obtain the mass (correction) tensor in thelaboratory-fixed (LF) frame, we have to account for the rotation from the BF frame to theLF frame. Let us represent this rotation with an O matrix, and thus the relation betweenthe LF mass matrix (mass-correction matrix), M ( A ), and the BF mass matrix, ¯ M ( ¯ A ),is M ij = (cid:0) O ¯ M O T (cid:1) ij (27)with A ij = (cid:0) O ¯ A O T (cid:1) ij (28)and ˜ G µν = g µα J iα (cid:0) O ¯ M O T (cid:1) ij J jβ g βν (29)= g µα (cid:0) J T O ¯ M O T J (cid:1) αβ g βν , (30)which completes the expression for the effective ˜ G µν matrix including the ¯ A mass correc-tion tensor computed in electronic structure theory with a certain embedding (BF) of theatomic nuclei. Besides the mass-correction matrix values, the effective ˜ G tensor contains themetric tensor and the Jacobi matrix elements—mathematical objects defined by a generalcoordinate transformation, and the orientation matrix which defines the BF frame used tocompute the mass-correction matrix elements.12hereby, the complete non-adiabatic kinetic energy operator expression in general curvi-linear coordinates is ˆ T (2) φ = − m n div M grad φ = − m n ˜ g − / ∂ µ ˜ g / ˜ G µν ∂ ν φ (31)with the volume element d V = ˜ g / d ξ . . . d ξ D . . The effective ˜ G matrix in general curvilinearcoordinates is calculated from the ¯ A ij body-fixed mass-correction tensor elements as:˜ G µν = g µα ( S ¯ M S T ) αβ g βν = 1 m n g µα ( SIS T ) αβ g βν − m g µα ( S ¯ A S T ) αβ g βν , (32)where S = J T O . (33) J , ˜ g / , g µν is the Jacobi tensor, the Jacobi determinant, and the contravariant metric tensorof the coordinate transformation, respectively. Furthermore, O is the rotation matrix whichtransforms ¯ A obtained from electronic structure theory in a selected body-fixed frame (whichmay be different from the body-fixed frame corresponding to the new ξ coordinates of therovibrational Hamiltonian) to the laboratory-fixed frame.The general expression, Eq. (31), can be re-written to the Podolsky form, following thesame reasoning as for constant masses in Eqs. (15)–( ?? ), by absorbing the ˜ g / factor in theoperator: ˆ T (2)Pod = − m n ˜ g − / ∂ µ ˜ g / ˜ G µν ∂ ν ˜ g − / , (34)for which the wave functions are normalized with the volume element:d V Pod = d ξ . . . d ξ D . (35)13 V. COMPUTATION OF THE SECOND-ORDER CORRECTION TERMS USINGAN EXPLICITLY CORRELATED GAUSSIAN BASIS SET
We solve the electronic Schr¨odinger equation, Eq. (2), and to compute the mass-correctioncorrection tensor with the QUANTEN program using floating explicitly correlated Gaussianfunctions (fECG) and the stochastic variational method [42] with regular refinements asimplemented in Refs. [31, 32]. Since the electronic problem with fixed nuclei does not havethe full O(3) rotation-inversion symmetry, which the pre-Born–Oppenheimer problem has,the integral expressions are considerably simpler and were obtained from the integrals derivedfor the generator function (see for example the Supplementary Material of Ref. [31]): g ( r ; A , s ) = exp (cid:20) − r T ( A ⊗ I ) r + s T r (cid:21) , (36)which is related to a floating explicitly correlated Gaussian functions (fECG) centered at R ∈ R N el as f ( r ; A , R ) = exp (cid:20) −
12 ( r − R ) T ( A ⊗ I )( r − R ) (cid:21) = exp (cid:20) − R T ( A ⊗ I ) R (cid:21) × exp (cid:20) − r T ( A ⊗ I ) r + R T ( A ⊗ I ) r (cid:21) , (37)The first term in the product is a constant with respect to the integration for the electroniccoordinates, and thus can be accounted for by simple multiplication.The non-linear parameters of the fECG basis functions ( R centers and A exponents) aregenerated in a stochastic variational optimization procedure and they are regularly refinedusing Powell’s method [43] to minimize the electronic energy, E (similarly to the pre-BOapproach [31, 32]). We have used the full point-group symmetry for the energy minimization,and have exploited the idempotency of the symmetry projector, so its explicit effect had to becalculated only for the ket functions. We have also exploited the convenient transformationproperties of fECGs under the effect of symmetry operators, which can be translated to thetransformation of the parameterization of an fECG [31, 33].In order to generate a potential energy curve for a diatomic molecule (and other correctionquantities), we have carried out a full basis-set optimization (energy minimization) at onlya few nuclear configurations. Then, we used the idea of Cencek and Kutzelnigg [28] for14escaling the centers upon making a small displacement of the nuclear configurations (alsoused by Pavanello and Adamowicz for the computation of triatomic molecules [44]). Bymaking sufficiently small nuclear displacements, one obtains a very good starting (non-linear)parameter set for the fECGs, which can be refined at a moderate computational cost. Forthe numerical examples shown later (Section VI and Paper II) we used a 0.1 bohr step sizefor the intermolecular distance and carried out 1-2 full refinement cycle was sufficient tomaintain the accuracy of the electronic energies along the potential energy curve. A. Diagonal Born–Oppenheimer correction
Since the pragmatic approach of Handy and co-workers [25–27], it is possible to computethe diagonal Born–Oppenheimer correction (DBOC) using the simple and general expressionin laboratory-fixed Cartesian coordinates: U ia ( R ) = 12 (cid:28) ∂ψ∂R ia (cid:12)(cid:12)(cid:12) ∂ψ∂R ia (cid:29) el (38)We compute the wave function derivatives numerically using the rescaling idea of Cencekand Kutzelnigg [28] (since the displacements were very small, on the order of 10 − bohr, thefull refinement of the basis set was not necessary): ∂ψ∂R ia ( R ) ≈ ψ ( R + ∆ ia ) − ψ ( R − ∆ ia )∆ ia (39)with the ∆ ia displacement vector (of the atomic nuclei), which labels a ∆ ia displacementalong the ia degree of freedom. Then, Eq. (38) becomes U ia ( R ) ≈ (1 − S ( ± ) ia ( R ))(∆ ia ) , (40)where we exploited that the electronic wave function is real and normalized and S ( ± ) ia ( R ) = h ψ ( R + 12 ∆ ia ) | ψ ( R − ∆ ia ) i el . (41)15 . Mass-correction tensor The mass-correction tensor, Eq. (5), for an isolated electronic state (
E, ψ )is written inCartesian coordinates as (the bar over A indicates that the electronic-structure computationsare carried out for a selected embedding of the nuclei):¯ A ia,jb ( R ) = 2 h ∂ R ia ψ | ( ˆ H (el) − E (el) ) − (1 − ˆ P ) ∂ R jb ψ i el . (42)In order to calculate the effect of the resolvent, we introduce an auxiliary basis set { f n , n =1 , , . . . , N aux } : ¯ A ia,jb ( R ) = 2 N aux X n =1 N aux X m =1 h ∂ R ia ψ | f n i el (cid:0) F − (cid:1) nm h f m | ∂ R jb ψ i el (43)with ( F ) nm = h f n | ( ˆ H (el) − E (el) )(1 − ˆ P ) | f m i el . (44)Representation of the resolvent by the direct summation over excited electronic states isimpractical, because it would require to compute (and converge) a very large number ofelectronic states to tightly converge the resolvent (and the mass-correction tensor). Instead,we ensure the convergence of the mass matrix elements by enlarging the auxiliary basis setsimilarly to Refs. [8, 14]. As for the DBOC, the wave function derivatives are computed byfinite differences using the fECG rescaling idea [28]: h ∂ R ia ψ | f n i el ≈ ia (cid:16) T (+) ia,n − T ( − ) ia,n (cid:17) , (45)where we have introduced the short notation T ( ± ) ia,n = h ψ ( R ± ∆ ia ) | f n i el = N ( ± ) ia X k =1 ( c ( ± ) ia ) k h ( g ( ± ) ia ) k | f n i el . (46)16n Eq. (46), ( g ( ± ) ia ) k is the k th fECG basis function obtained for the R ± ∆ ia nuclear geometryand ( c ( ± ) ia ) k is the corresponding linear combination coefficient of the electronic state ψ ( r ; R ± ∆ ia ) resulting from the variational solution of the electronic Schr¨odinger equation in the( g ( ± ) ia ) k basis set ( k = 1 , , . . . , N ( ± ) ia ). a. Evaluation of the resolvent Instead of directly constructing the F − matrix, it iscomputationally more accurate and stable (see also Ref. [14]) to consider( F x j b ) m = h f m | ∂ R jb ψ i el (47)solve this system of linear equations for x j b , and then evaluate Eq. (43) as the sum of vectorproducts. Furthermore, we consider F ( E (el) + ε ) instead of F ( E (el) ) with a small ε real orimaginary value in order to avoid numerical instabilities due to having either explicitly orimplicitly the inverse of a singular matrix (note that E is an eigenvalue of ˆ H in Eq (44)). Theoptimal value of ε is determined in a series of computations by maximizing both numericalstability and accuracy. b. Auxiliary basis set In numerical applications, it is important to ensure that Eq. (43)is converged with respect to the enlargement of the auxiliary basis set. A reasonable startingpoint is to use the basis set optimized for the solution of the electronic Schr¨odinger equation.This choice often gives an accurate estimate (at least within a 1-2 % of the exact value), butusually additional functions are necessary to obtain really accurate mass-correction terms.Pachucki and Komasa optimized the auxiliary basis set by using the variational property ofthe rotational and vibrational mass-correction terms arising in the curvilinear expression ofthe mass-correction tensor for diatomic molecules [8]. In Ref. [14] a plane-wave expansionwas used both for the electronic-state representation and also for the auxiliary basis set.We do not have any direct and general optimization strategy to build an optimal andconverged auxiliary basis set of fECGs for general polyatomic and polyelectronic molecules.Instead, we simply enlarge the auxiliary basis set until convergence was achieved by usingthe electronic basis sets optimized at neighbohring nuclear configurations. As a practicalhandle, we have noticed that if the resolution of identity was poorly represented by theauxiliary basis functions, the resolvent was also inaccurate. The resolution of identity in a17on-orthogonal set of functions isˆ I = N aux X i =1 N aux X j =1 | f i i · h f i | ˆ I | f j i − · h f j | = N aux X i =1 N aux X j =1 | f i i · S − ij · h f j | (48)and we evaluate the derivative overlaps also by inserting the auxiliary basis set: h ∂ k Ψ | ∂ k Ψ i el ≈ h ∂ k Ψ | f i i el · S − ij · h f j | ∂ k Ψ i el with S ij = h f i | f j i el . (49)The accuracy of Eq. (49) is not a sufficient condition for the convergence of Eq. (43), it ismerely a useful indicator for the reliability of the results. c. Frame transformation It is important to notice that the mass-tensor elements areevaluated with a certain frame definition of the nuclear geometry. This choice of the body-fixed frame in the electronic-structure computations is described by the O rotation matrixwhich connects this body-fixed frame and the laboratory-fixed frame. The body-fixed frameused to compute the mass-correction tensor may be different from the body-fixed frameused in the rovibrational computations. The general curvilinear form of the second-ordernon-adiabatic kinetic energy operator, Eqs. (25)–(33), includes the corresponding framedefinitions. 18 . APPLICATIONS FOR HOMONUCLEAR DIATOMICSA. Evaluation of the general curvilinear kinetic energy expressions for homonu-clear diatomic molecules The quantum mechanical motion of atomic nuclei in di- and poly-atomic molecules isefficiently described if the laboratory-frame (LF) Cartesian coordinates, R , . . . , R N , arereplaced with internal coordinates, orientation angles, and the Cartesian coordinates of thenuclear center of mass (NCM), R NCM .The general curvilinear kinetic-energy operator expressions developed in Section III canbe directly implemented in the GENIUSH program and the derived effective mass-matrixexpressions can be evaluated at grid points similarly to the mass-weighted metric tensor asit has been implemented for the constant mass case in Ref. [30]. In this section, we shallgo through the calculations with the coordinate-dependent mass-correction tensor step bystep and derive the expressions explicitly in order to better understand the formalism andto highlight the formal and numerical properties of the derived expressions for the simplecase of homonuclear diatomic molecules. a. Preparatory calculations with spherical polar coordinates
The transformation fromthe ( r , r , r ) “flat coordinates” to spherical polar coordinates, ρ ∈ [0 , ∞ ), θ ∈ ( − π, π ), and φ ∈ [0 , π ), is defined by r r r = ρ sin θ cos φρ sin θ sin φρ cos θ . (50)the corresponding Jacobian matrix is [57]: J sp = sin θ cos φ ρ cos θ cos φ − ρ sin θ sin φ sin θ sin φ ρ cos θ sin φ ρ sin θ cos φ cos θ − ρ sin θ . . (51)19he covariant metric tensor, the elements of which are calculated as g µν = δ ij ∂r i ∂ξ µ ∂r i ∂ξ ν , is: g µν = . .. ρ .. . ρ sin θ , (52)and it is inverted to obtain, the contravariant metric tensor: g µν = . .. ρ .. . ρ sin θ . (53)The Jacobi determinant for this coordinate transformation reads as:˜ g − / = [det( g µν )] − / = det( J sp ) = ρ sin θ. (54)The electronic structure computations and the evaluation of the mass-correction tensor iscarried out in a frame selected for the atomic nuclei. For homonuclear diatomic moleculesit is a natural choice to position the nuclei symmetrically with respect to the origin alongthe z axis. The rotation matrix from this “ z -axis embedding” to the laboratory-fixed frame(expressed with the θ and φ spherical angles) is: O sp = cos φ − sin φ . sin φ cos φ .. . · cos θ . sin θ. . − sin θ . cos θ = cos θ cos φ − sin φ sin θ cos φ cos θ sin φ cos φ sin θ sin φ − sin θ . cos θ . (55)20hus, the S transformation matrix defined in Eq. (33) for these curvilinear coordinates andbody-fixed frame is: S sp = J Tsp O sp = . . r . .. r sin θ . . (56)Note also the relationship between the covariant metric tensor and the S sp transformationmatrix: g − S sp S Tsp = S sp S Tsp g − = I . (57) B. Transformation to spherical polar and center-of-mass Cartesian coordinates
For diatomic molecules, the six laboratory-fixed Cartesian coordinates are replacedwith three NCM Cartesian coordinates and the three spherical polar coordinates. So, the( R , . . . , R ) → ( ξ , . . . , ξ ) transformation from the R i “flat coordinates” to the ξ µ newcoordinates is R R R R R R = − r + R − r + R − r + R r + R r + R r + R = − ρ sin θ cos φ + R − ρ sin θ sin φ + R − ρ cos θ + R ρ sin θ cos φ + R ρ sin θ sin φ + R ρ cos θ + R (58)where the r internuclear displacement vector is written in terms of spherical polar coordi-nates and R NCM collects the center of mass coordinates of the atomic nuclei (calculated with21he nuclear masses). The Jacobian matrix of this transformation is J = ′′ ∂R i ∂ξ µ ′′ = . . − J sp . .. . . . J sp . .. . , (59)where J sp is the Jacobian matrix corresponding to the spherical polar coordinates given inEq. (51). The covariant metric tensor is g µν = ′′ δ ij ∂R i ∂ξ µ ∂R i ∂ξ ν ′′ = . . . . .. ρ . . . .. . ρ sin θ . . .. . . . .. . . . .. . . . . , (60)and the Jacobi determinant reads as ˜ g / = ρ sin θ. (61)The rotation of the 6-dimensional position vectors expressed in the z -axis embedding, whichis also the body-fixed frame which we use in electronic structure theory to compute themass-correction tensor, and in the laboratory frame is the direct product of the rotationmatrix given in Eq. (55) with the 2 × O = O sp ⊗ I , (62)22nd thus, the S transformation matrix defined in Eq. (33) for this 6-dimensional coordinatetransformation is S = J T O = − S sp 12 S sp O sp O sp , (63)where S sp was given in Eq. (56). C. Calculation of the BO and non-adiabatic kinetic energy operators in curvilinearcoordinates
For a start, let us consider (and reproduce) the kinetic energy operator for the constant-mass case, i.e., for a constant, diagonal mass matrix (1 /m n I with the I × T (0) = − m n div I grad= − m n ˜ g − / ∂ µ ˜ g / g µα ( SI S T ) αβ g βν ∂ ν = − m n ˜ g − / ∂ µ ˜ g / g µν ∂ ν = − m n ρ sin θ (cid:16) ∂ ρ , ∂ θ , ∂ φ , ∂ R , ∂ R , ∂ R (cid:17) ρ sin θ . . . . .. ρ . . . .. . ρ sin θ ) . . .. . . . .. . . . .. . . . . ∂ ρ ∂ θ ∂ φ ∂ R ∂ R ∂ R = −
12 1 ρ sin θ ( ∂ ρ , ∂ θ , ∂ φ ) ρ sin θ m n . .. m n ρ .. . m n ρ sin θ ∂ ρ ∂ θ ∂ φ + ˆ T (0)NCM , (64)23here ∂ µ = ∂/∂ξ µ and the kinetic energy operator of the nuclear center of mass was intro-duced as: ˆ T (0)NCM = −
12 12 m n (cid:16) ∂ R + ∂ R + ∂ R (cid:17) . (65)It is important to remember that ˆ T (0)NCM is separable from the rovibrational kinetic energyoperator.In Eq. (64) we have reproduced the well-known diatomic rotation-vibration kinetic energyexpression, which is the first term including the inverse of the reduced mass of the twoidentical nuclei, 1 /µ = 1 /m n + 1 /m n = 2 /m n , and the translational kinetic energy of thenuclear center of mass, which is the second term including the inverse of the total nuclearmass, 1 /M = 1 / (2 m n ).Next, let us calculate the curvilinear kinetic energy operator including the mass-correctiontensor for the present example. In addition to the general expressions, to better highlighthow the formalism is used, we also give the numerical results for the H +2 molecular ion for ρ = 3 bohr internuclear separation.The general form of the mass-correction tensor for a homonuclear diatomic molecule inthe z -axis embedding can be deduced from simple symmetry arguments ( A, B, a, b ∈ R + ,further conditions apply, vide infra ):¯ A ij = A . . a . .. A . . a .. . B . . ba . . A . .. a . . A .. . b . . B { ρ ∈ R + , z -axis embedding } , (66)24hich takes the numerical values¯ A ij = . . . . . .. . . . . .. . . . . . . . . . . .. . . . . .. . . . . . { for ρ = 3 bohr } . (67)The curvilinear expression for A αβ is obtained by using the S transformation matrix givenin Eq. (63) of the present coordinate transformation: A αβ = ( S ¯ A S T ) αβ = B − b . . . . .. A − a ρ . . . .. . A − a ρ sin θ . . .. . . B + b ) . .. . . . A + a ) .. . . . . A + a ) { ρ ∈ R + , ( ρ, θ, φ ) } , (68)which equals to = . . . . . .. . ρ . . . .. . . ρ sin θ . . .. . . . . .. . . . . .. . . . . . { for ρ = 3 bohr } . (69)25n the next step, we multiply A αβ with the contravariant metric tensor, and obtain thefollowing simple expression A µ β = g µα A αβ = g µα ( S ¯ A S T ) αβ = B − b . . . . .. A − a . . . .. . A − a . . .. . . B + b . .. . . . A + a .. . . . . A + a { ρ ∈ R + , ( ρ, θ, φ ) } (70)= . . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . . { for ρ = 3 bohr } . (71)26urthermore, by multiplying this matrix with g βν from the right, we get A µν = A µ β g βν (72)= ρ sin θ B − b ) . . . . .. A − a ) ρ . . . .. . A − a ) ρ sin θ . . .. . . B + b . .. . . . A + a .. . . . . A + a { ρ ∈ R + , ( ρ, θ, φ ) } (73)= ρ sin θ . . . . . .. . ρ . . . .. . . ρ sin θ . . .. . . . . .. . . . . .. . . . . . { ρ = 3 , ( ρ, θ, φ ) } , (74)and we write down the effective ˜ G µν matrix defined in Eq. (32) as:˜ G µν = 1 m n g µα ( SI S T ) αβ g βν − m n2 g µα ( S ¯ A S T ) αβ g βν (75)= (cid:20) δ µ β − m n A µ β (cid:21) m n g βν (76)= m n (cid:20) − A ρ ρm n (cid:21) . . . . .. m n (cid:20) − A θ θm n (cid:21) ρ . . . .. . m n " − A φ φm n ρ θ . . .. . . m n " − A R R m n . .. . . . m n " − A R R m n .. . . . . m n " − A R R m n = m n h − B − bm n i . . . . .. m n h − A − am n i ρ . . . .. . m n h − A − am n i ρ θ . . .. . . m n h − B + bm n i . .. . . . m n h − A + am n i .. . . . . m n h − A + am n i (77)for ρ ∈ R + (spherical polar coordinates, z -axis embedding).27y inserting the mass-correction values for ρ = 3 bohr we have= m n h − . m n i . . . . .. m n h − . m n i ρ . . . .. . m n h − . m n i ρ sin θ . . .. . .. . . m n h − m n i I . . . for H +2 with ρ = 3 bohr (spherical polar coordinates, z -axis embedding). (78)Using this effective ˜ G µν matrix and exploiting that its off-diagonal elements are all zero,we write the second-order nonadiabatic kinetic energy operator for a diatomic molecule (forcomparison see the constant-mass operator in Eq. (64)–(65)) asˆ T (2) = −
12 ˜ g − / ∂ µ ˜ g / ˜ G µν ∂ ν = −
12 1 ρ sin θ ( ∂ ρ , ∂ θ , ∂ φ ) ρ sin θ ˜ G ρρ . .. ˜ G θθ .. . ˜ G φφ ∂ ρ ∂ θ ∂ φ + ˆ T (2)NCM . (79)The second-order non-adiabatic kinetic energy operator for the nuclear center of mass motionis ˆ T (2)NCM = −
12 ˜ G NCM (cid:16) ∂ R NCM , + ∂ R NCM , + ∂ R NCM , (cid:17) (80)and ˜ G NCM corresponds to the diagonal elements of the lower right block of the effectivetensor given in Eq. (77).The resulting expressions have a couple of important and interesting properties. By look-ing at the numerical values computed for H +2 (at ρ = 3 and at other internuclear distancevalues) we may observe that the non-zero, diagonal elements in the NCM block can be re-written into a physically meaningful form using the 1 − x ≈ (1 + x ) − approximation, valid28or small x values:12 m n (cid:20) − m n (cid:21) ≈ m n [1 + 1 / (2 m n )] = 12 m n + 1 = 12 ˜ M tot (81)where ˜ M tot equals the total mass of the H +2 molecular ion, two protons and one electron, inatomic units ( m el = 1). Using the same approximation, we re-write the diagonal elements ofthe “rovibrational” block as2 m n (cid:20) − δ ˜ mm n (cid:21) ≈ m n (1 + δ ˜ m/m n ) = 2 m n + δ ˜ m = 2˜ m n (82)in which ˜ m n can be interpreted as an effective nuclear mass. The “effective vibrational mass”,corresponding to the ρ degree of freedom, is thus2 m n (cid:20) − A ρ ρ m n (cid:21) ≈ m n + A ρ ρ = 2˜ m vib , and ˜ m vib = m n + A ρ ρ (83) { ρ ∈ R + bohr, spherical polar coordinates } ,which has the value for ρ = 3 bohr:2 m n (cid:20) − . m n (cid:21) ≈ m n + 0 .
424 = 2˜ m vib and ˜ m vib = m n + 0 .
424 (84) { ρ = 3 bohr, spherical polar coordinates } .The effective rotational mass, which is the same for the θ and φ degrees of freedom ( A Ω Ω = A θ θ = A φ φ ), is 2 m n (cid:20) − A Ω Ω m n (cid:21) ≈ m n + A Ω Ω = 2˜ m rot , and ˜ m rot = m n + A Ω Ω (85) { ρ ∈ R + , spherical polar coordinates } m n (cid:20) − . m n (cid:21) ≈ m n + 0 .
144 = 2˜ m rot and ˜ m rot = m n + 0 .
144 (86) { ρ = 3 bohr, spherical polar coordinates } .There are two important, general properties of the second-order non-adiabatic kineticenergy operator, ˆ T (2) —which manifest themselves also in this example—, in relation withthe transformation of the effective mass matrix to translationally invariant (TI) coordinatesand the Cartesian coordinates of the nuclear center of mass (NCM). A formal derivationfor both properties was given in Ref. [14]. First, the effective mass matrix is always blockdiagonal in a TI-NCM representation and there is not any coupling between the blockcorresponding to translationally invariant coordinates and the block of the nuclear-center-of-mass Cartesian coordinates. Hence, the translational kinetic energy can always be separatedfrom the rovibrational kinetic energy operator in ˆ T (2) (which is indeed a very importantproperty). Second, the translational kinetic energy term in ˆ T (2) , can be rearranged (withinthe 1 − /m n ≈ (1 + m n ) − approximation) to a form in which the mass associated to thetranslational degrees of freedom is the total mass of the molecule (nuclei plus electrons).Hence, in ˆ T (2) not only the rovibrational but also the translational kinetic energy gains acorrection term, which increases the total nuclear mass with the mass of the electrons.In short, the overall translation remains exactly separable from the internal (rotational-vibrational) degrees of freedom in the second-order kinetic-energy operator, but has aneffective mass, which is equal to the total mass of the molecule. It is interesting to notethat these general properties provided the starting point for the non-adiabatic theory forsecond-order “mass” effects of Kutzelnigg [10]. This theory was later used by Jaquet andKutzelnigg for diatomics [11], and by Jaquet and Khoma for modeling non-adiabatic effectsin H +3 [45, 46]. 30 I. COMPUTATION OF THE MASS CORRECTION FUNCTIONS AND NON-ADIABATIC CORRECTIONS TO THE ROVIBRATIONAL ENERGIES
We consider the second-order non-adiabatic Hamiltonian with the effective ˜ G µν tensor in-cluding the mass-correction terms, see Eqs. (31)–(32), (cid:18) − m n ˜ g − / ∂ µ ˜ g / ˜ G µν ∂ ν + U + V (cid:19) φ = E (n) φ (87)and U = P Ni =1 P a U ia /m n is the diagonal Born–Oppenheimer correction to the BO potentialenergy, V .For diatomic molecules, there is not any coupling term between the rotational and vi-brational degrees of freedom neither in the BO, ˆ T (0)rv = ˆ T (0) − ˆ T (0)NCM in Eq. (64), nor in thesecond-order kinetic energy operator, ˆ T (2)rv = ˆ T (2) − ˆ T (2)NCM in (79). Hence, the angular partof ˆ T (2)rv can be integrated with the Y JM ( θ, φ ) spherical harmonic functions (similarly to thestandard solution of diatomics with ˆ T (0)rv ), and we are left with the numerical solution of theradial equation: (cid:18) −
12 1 ρ ∂∂ρ ρ m n (cid:20) − A ρ ρ m n (cid:21) ∂∂ρ + J ( J + 1) ρ m n (cid:20) − A Ω Ω m n (cid:21) + U ( ρ ) + V ( ρ ) (cid:19) f J ( ρ ) = E (n) J f J ( ρ ) (88)with the volume element ρ d ρ . Instead of solving Eq. (88), we proceed similarly to Pachuckiand Komasa [8] and use the operator identity1 ρ ∂∂ρ ρ X ( ρ ) ∂∂ρ f ( ρ ) = 1 ρ ∂∂ρ X ( ρ ) ∂∂ρ ρf ( ρ ) − ρ ∂X∂ρ f ( ρ ) (89)to obtain (cid:18) − ∂∂ρ m n (cid:20) − A ρ ρ m n (cid:21) ∂∂ρ − ρ m ∂ A ρ ρ ∂ρ + J ( J + 1) ρ m n (cid:20) − A Ω Ω m n (cid:21) + U ( ρ ) + V ( ρ ) (cid:19) φ J ( ρ ) = E (n) J φ J ( ρ ) (90)with the volume element d ρ . To obtain rovibrational (in fact, rovibronic) energies and wavefunctions, we solve Eq. (90) using the discrete variable representation (DVR) and associated31aguerre polynomials, L ( α ) n with α = 2 for the radial (vibrational) degree of freedom. TheDVR points are scaled to an [ R min , R max ] interval. The n number of DVR points and func-tions, as well as R min and R max are determined as convergence parameters, and their typicalvalue is around n = 300 − R min = 0 . R max = 30 −
100 bohr.32 . Mass-correction curves and non-adiabatic corrections to the rovibrational en-ergies of H +2 in its ground electronic state, ˜ X Σ + g Using the variational method with fECG basis functions and the computational proceduredescribed in Section IV we have computed the mass-correction functions for the H +2 molecularion in its ˜ X Σ +g (ground) electronic state over the interval of the ρ internuclear distance. Theresulting rotational and vibrational mass-correction curves computed in the present workfor H +2 are visualized in Figure 1 (the computed values are deposited in the SupplementaryMaterial). The figure also shows the adiabatic potential energy curve (including DBOC) toallow a visual comparison of the improtant features of the mass-correction and the potentialenergy curve.In Table I we present the non-adiabatic corrections obtained with these mass-correctioncurves in comparison with Moss’s results [47, 48]. (The non-adiabatic correction to a rovibra-tional state is defined as the difference of the eigenvalues obtained with the BO Hamiltonianwith the potential energy curve including the DBOC minus the corresponding eigenvaluesobtained using the full non-adiabatic Hamiltonian.) Excellent agreement is observed over allbound and long-lived resonance states. (Note that the disagreement for the ( v, J ) = (0 , ,
4) states remains unexplained and might be due to a typographical error in theearlier references.)As to the convergence of the presented results, the electronic state, ψ , used to computethe mass-correction curves, had an energy eigenvalue E converged to better than 0.1 µE h .The auxiliary basis set (constructed from functions of Σ +g and Π +u symmetries) used tobuild the F matrix, Eq. (44), was compiled using the parameters optimized for the ground-state energy. The most stable and accurate mass-correction values were obtained by usinga complex ε value ( ε = i · − ) for computing the resolvent in Eq. (47). This setup wassufficient to converge the non-adiabatic correction energies accurate to at least < .
005 cm − .In order to assess the necessary accuracy of the electronic energy and wave function toobtain accurate non-adiabatic rovibrational energies, we have observed that the mass cor-rection functions (and hence, non-adiabatic rovibrational energies) were robust with respectto inaccuracies in the electronic state: a 10–20 µE h error of the electronic energy intro-duced a less than 0.5 % error in the mass-correction values, which had an almost negligible( < .
005 cm − ) effect on the non-adiabatic correction energies. At the same time, it was33ound to be important that a sufficiently large auxiliary basis set is used in particular forcomputing the vibrational mass-correction function at large internuclear separations.Using the computed, rigorous mass-correction curves, it will be interesting to comparethe (second-order) non-adiabatic rovibrational bound and resonance states with their three-particle variational counterparts recently computed by Korobov [49], and the recently mea-sured shape resonances of the hydrogen molecular ion and its deuterated isotopologue[50, 51]. This will offer a unique opportunity to directly test the numerical accuracy ofsecond-order non-adiabatic theory with respect to the numerically exact few-particle varia-tional results. Such a direct comparision has recently become available in the literature forthe rotational states of the hydrogen molecule [52].34 V + U [ E h ] ρ [bohr]PES+DBOC δ ˜ m [ m e ] ρ [bohr] δ ˜ m rot δ ˜ m vib FIG. 1: H +2 molecular ion in its ˜ X Σ +g electronic state (ground state): mass correction functionsto the vibrational, δ ˜ m vib = A ρ ρ , and the rotational, δ ˜ m rot = A θ θ = A φ φ degrees of freedom whenthe kinetic energy operator is expressed in spherical polar coordinates, Eq. (70). For the definitionof the effective rotational and vibrational mass in diatomic molecules see Eqs. (83)–(86). The thinline takes the constant value 1 / . m e , which corresponds to an equal share of the singleelectron between the two atomic nuclei. A B L E I : N o n - a d i a b a t i cc o rr ec t i o n s f o r t h e r o v i b r a t i o n a l s t a t e s o f H + : c o m p a r i s o n ( d e v i a t i o n ) f r o m M o ss ’ r e s u l t s ( c o m pu t e d w i t h m p / m e = . nd E h = . c m − ) [ , ]. v / J − . − .
002 0.059 a − . − .
001 0 . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . . − . − . − . − .
002 0 . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . b − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . • − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . • − . − . − . − . − . − . − . − . − . − . − . − . − .
001 0 . − . • − . − . − . − . − . − . − . − . − . − . − . − .
001 0 . • − . − . − . − . − . − . − . − . − .
001 0 . • − .
001 0 . − . − . − . − .
001 0 . •
18 0 .
000 0 .
000 0 .
000 0 . •
19 0 .
000 0 . v / J
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 410 − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . • − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − .
003 0 . • − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . e • − . − . − . − . − . − . − . − . − . − . − . − . − . − .
002 0 . d • − . − . − . − . − . − . − . − . − . − . − . − .
002 0 . • − . − . − . − . − . − . − . − . − . − .
002 0 . • − . − . − . − . − . − . − . − . − . • − . − . − . − . − . − . − .
002 0 . c • − . − . − . − . − .
001 0 . • − . − . − . − . • − . − . • • a,b The stability and estimated accuracy of our correction terms is much better than the order of magnitude of the deviation. c,d,e
The accuracy of these results has been labelled in the present work as well as Ref. [47]. The symbol • is used for states which have been predictedby Moss without giving their energy and non-adiabatic correction explicitly in Refs. [47, 48]. II. SUMMARY AND CONCLUSIONS
General curvilinear expressions have been derived for the second-order non-adiabatic ki-netic energy operator. The derivations have been carried out using the general Jacobianand metric tensors of the coordinate transformation, thereby they are in direct connectionwith the numerical-kinetic energy operator approach (with constant masses) also used inthe GENIUSH protocol [30]. While for the constant-mass case one has to transform the“div grad” operator to curvilinear coordinates, we had to consider the “div M grad” op-erator with the M coordinate-dependent tensor quantity within the operator. As a result,should a general mass-correction tensor surfaces become available for polyatomic molecules,their implementation in the (polyatomic) rovibrational GENIUSH program has been madestraightforward.At the moment, we are not aware of any widely available electronic structure packagewhich we could use to compute the mass-correction tensor over a wide range of nuclear con-figurations of polyatomic molecules (note however that Ref. [14] points into a very promis-ing direction). Hence, we have modified (simplified) our in-house preBO code originallydeveloped for the computation of isolated few-particle quantum systems [31–33, 53, 54]to few-electron systems being in interaction with external charges (nuclei) and use float-ing explicitly-correlated Gaussian functions for the spatial part of the basis set. We re-port the first numerical results using this computational setup and present the rigorousmass-correction values for the H +2 molecular ion in its ground electronic state over the[0 . , +2 within a few 0.001 cm − for the bound and long-lived resonance states. In a forthcomingpaper [24] further applications of the developed formalism and methodology are presentedfor the He +2 molecular ion in its ground electronic state for which a highly accurate poten-tial energy and DBOC curve is available [55]. It is observed that full account of the rigorousmass correction functions substantially reduce the deviation of experimental [23, 56] andcomputed rotational(-vibrational) energy intervals. Supplementary Material
The Supplementary Material contains non-adiabatic mass correction values for H +2 ( ˜ X Σ +g ).37 cknowledgment Financial support of the Swiss National Science Foundation through a PROMYS Grant (no.IZ11Z0 166525) is gratefully acknowledged. The author is thankful to Krzysztof Pachuckifor discussions about non-adiabatic perturbation theory, and to Federica Agostini for theexplanation of the mass-correction tensor derivation starting from exact factorization.38
1] G. A. Fisk and B. Kirtman, J. Chem. Phys. , 3516 (1964).[2] P. R. Bunker and R. E. Moss, Mol. Phys. , 417 (1977).[3] P. R. Bunker and R. E. Moss, J. Mol. Spectrosc. , 217 (1980).[4] P. R. Bunker, C. J. McLarnon, and R. E. Moss, Mol. Phys. , 417 (1977).[5] D. W. Schwenke, J. Chem. Phys. , 1693 (2001).[6] D. W. Schwenke, J. Phys. Chem. A , 2352 (2001).[7] K. Pachucki and J. Komasa, J. Chem. Phys. , 034102 (2008).[8] K. Pachucki and J. Komasa, J. Chem. Phys. , 164113 (2009).[9] K. Pachucki and J. Komasa, J. Chem. Phys. , 204314 (2012).[10] W. Kutzelnigg, Mol. Phys. , 2627 (2007).[11] R. Jaquet and W. Kutzelnigg, Chem. Phys. , 69 (2008).[12] L. G. Diniz, J. R. Mohallem, A. Alijah, M. Pavanello, L. Adamowicz, O. L. Polyansky, andJ. Tennyson, Phys. Rev. A , 032506 (2013).[13] S. Weigert and R. G. Littlejohn, Phys. Rev. A , 3506 (1993).[14] A. Scherrer, F. Agostini, D. Sebastiani, E. K. U. Gross, and R. Vuilleumier, Phys. Rev. X ,031035 (2017).[15] A. Scherrer, F. Agostini, D. Sebastiani, E. K. U. Gross, and R. Vuilleumier, J. Chem. Phys. , 074106 (2015).[16] R. M. Herman and A. Asgharian, J. Mol. Spectrosc. , 305 (1966).[17] R. M. Herman and J. F. Ogilvie, Adv. Chem. Phys. , 187 (1998).[18] J. F. Ogilvie, The Vibrational and Rotational Spectrometry of Diatomic Molecules (AcademicPress, 1998).[19] S. P. A. Sauer, H. J. A. Jensen, and J. F. Ogilvie, Adv. Quant. Chem. , 319 (2005).[20] K. L. Bak, S. P. A. Sauer, J. Oddershede, and J. F. Ogilvie, Phys. Chem. Chem. Phys. ,1747 (2005).[21] P. A. Braun, T. K. Rebane, and K. Ruud, Chem. Phys. , 299 (1996).[22] G. Panati, H. Spohn, and S. Teufel, ESAIM: Mathematical Modelling and Numerical Analysis , 297 (2007).[23] L. Semeria, P. Jansen, and F. Merkt, J. Chem. Phys. , 204301 (2016).
24] E. M´atyus, Computation of the non-adiabatic mass-correction functions and rovibrationalenergy levels for the He +2 molecular ion in its ground electronic state, J. Chem. Phys. accepted(2018).[25] N. C. Handy and A. M. Lee, Chem. Phys. Lett. , 425 (1996).[26] A. G. Ioannou, R. D. Amos, and N. C. Handy, Chem. Phys. Lett. , 52 (1996).[27] N. C. Handy, Y. Yamaguchi, and H. F. Schaefer III, J. Chem. Phys. , 4481 (1986).[28] W. Cencek and W. Kutzelnigg, Chem. Phys. Lett. , 383 (1997).[29] W. Kutzelnigg, Mol. Phys. , 909 (1997).[30] E. M´atyus, G. Czak´o, and A. G. Cs´asz´ar, J. Chem. Phys. , 134112 (2009).[31] E. M´atyus and M. Reiher, J. Chem. Phys. , 024104 (2012).[32] E. M´atyus, J. Phys. Chem. A , 7195 (2013).[33] E. M´atyus, Mol. Phys. p. arXiv:1801.05885 (2018).[34] J. Tennyson, J. Chem. Phys. , 120901 (2016).[35] T. Carrington Jr., J. Chem. Phys. , 120902 (2017).[36] D. Luckhaus, J. Chem. Phys. , 1329 (2000).[37] D. Lauvergnat and A. Nauts, J. Chem. Phys. , 8560 (2002).[38] B. Podolsky, Phys. Rev. , 812 (1928).[39] C. F´abri, E. M´atyus, and A. G. Cs´asz´ar, J. Chem. Phys. , 074105 (2011).[40] J. Sarka, A. G. Cs´asz´ar, S. C. Althorpe, D. J. Wales, and E. M´atyus, Phys. Chem. Chem.Phys. , 22816 (2016).[41] J. Sarka, A. G. Cs´asz´ar, and E. M´atyus, Phys. Chem. Chem. Phys. , 15335 (2017).[42] Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum-Mechanical Few-BodyProblems , 184303 (2012).[45] R. Jaquet and M. Khoma, J. Phys. Chem. A , 7016 (2017).[46] R. Jaquet and M. V. Khoma, Mol. Phys. p. doi.org/10.1080/00268976.2018.1464225 (2018).[47] R. E. Moss, Mol. Phys. , 1541 (1993).
48] R. E. Moss, Mol. Phys. , 195 (1996).[49] V. I. Korobov, Mol. Phys. , 93 (2018).[50] M. Beyer and F. Merkt, Phys. Rev. Lett. , 093001 (2016).[51] M. Beyer and F. Merkt, J. Mol. Spectrosc. , 147 (2016).[52] K. Pachucki and J. Komasa, Phys. Chem. Chem. Phys. , 247 (2018).[53] E. M´atyus, J. Hutter, U. M¨uller-Herold, and M. Reiher, Phys. Rev. A , 052512 (2011).[54] E. M´atyus, J. Hutter, U. M¨uller-Herold, and M. Reiher, J. Chem. Phys. , 204302 (2011).[55] W.-C. Tung, M. Pavanello, and L. Adamowicz, J. Chem. Phys. , 104309 (2012).[56] A. Carrington, C. H. Pyne, and P. J. Knowles, J. Chem. Phys. , 5979 (1995).[57] In matrices, we shall use “.” for “0” to enhance readability in matrix expressions., 5979 (1995).[57] In matrices, we shall use “.” for “0” to enhance readability in matrix expressions.