On the Equivalence of Different Methods for ... - ACS Publications

Article pubs.acs.org/JCTC

On the Equivalence of Different Methods for Calculating Resonances: From Complex Gaussian Basis Set to Reflection-Free Complex Absorbing Potentials via the Smooth Exterior Scaling Transformation Anael Ben-Asher and Nimrod Moiseyev* Schulich Faculty of Chemistry, Technion-Israel Institute of Technology, Haifa 32000, Israel ABSTRACT: The work compares different methods to compute metastable states (resonances) using Gaussian functions. The subject is of current interest, as methods of computing resonances using ab initio methodologies have accumulated, and it is relevant to shed light on their differences and similarities. Since Gaussian functions are usually used in quantum chemistry, we focus on their use in the present. The illustrative numerical example is for a single particle problem. However, one can learn much from the latter on a many-particle quantum chemistry study. Our findings show that the calculations of molecular resonances by the use of complex Gaussian basis functions and by adding complex absorbing potentials to the original molecular Hamiltonian are closely equivalent one to another.

the experimental and theoretical papers published on ICD.4,5 The need for accurate complex potential energy is explained in ref 6, where the large sensitivity of the calculated cross sections due to small changes in the CPES has been shown. The potential energy surfaces are obtained by solving the electronic time-independent Schrödinger equation within the framework of the Born−Oppenheimer (BO) approximation. The BO molecular Hamiltonian is hermitian when the eigenfunctions associated with the bound and the quasi-discrete continuum eigenfunctions are expanded in terms of square integrable basis functions. The eigenvalues of the hermitian BO Hamiltonian are real. Therefore, only real potential energy surfaces which are associated with bound electronic states can be calculated using the hermitian BO Hamiltonian when Gaussians are used as a basis set. Note that the Gaussians are the most widely used basis functions in electronic structure calculations. The calculations of CPES require calculations of the eigenfunctions of the BO Hamiltonian by imposing on the asymptotes of the eigenfunctions outgoing boundary conditions. This requirement implies that the obtained solutions of the time-independent Schrödinger equation are the poles of the scattering matrix. The poles of the scattering matrix, which are obtained by imposing out going boundary conditions on the solutions of the time-independent Schrödinger equation, are bound and resonance states that decay at the asymptotes. The asymptotes of these solutions can be described as

1. INTRODUCTION AND MOTIVATION Atomic and molecular autoionization resonances are associated with complex eigenvalues of the atomic/molecular Hamiltonian systems. The real parts of the complex resonance eigenvalues provide the energy positions of the autoionization states (e.g., energies above the ionization thresholds). The imaginary parts of the complex resonance eigenvalues provide the autoionization decay rates (inverse lifetimes) of the metastable atomic/molecular systems. The molecular resonance decay rates (widths) heavily depend on the structure of the molecular system under study. In some geometrical structures, the molecules have very short lifetimes, while in other geometrical structures the same molecules have extremely long lifetime (very small widths). In some cases, even the short-lived metastable autoionization molecular states become bound states by varying the geometrical structure of the molecules. A simple example for this behavior is molecular hydrogen anion, H−2 , that at the equilibrium distance of H−2 the molecular anion in its ground electronic state has a lifetime which is equal to about one femtosecond, while beyond a critical bond length the ground state of H−2 is a bound state (at the dissociation limit, H and a stable H− are obtained, and no ionization takes place). The complex eigenvalue of the molecule as a function of its nuclear coordinates is defined as a complex-potential-energy surface (CPES), although only the real part of the complex molecular eigenvalues is associated with the energy of the molecule. The solution of the Schrödinger equation with complex potential energy surfaces (CPES) introduces a possibility to couple the electronic and the nuclear coordinates when autoionization takes place. There is a need for CPES in particular for providing predictions of new observable phenomena and for the analysis of certain experimental data. See the most recent publications of Narevicius and his coworkers on cold molecular collisions for examples.1−3 See also © XXXX American Chemical Society

Received: January 18, 2016

A

DOI: 10.1021/acs.jctc.6b00059 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

In the following, we discuss the similarity transformations, which can be regarded as the basis of the calculations of molecular autoinozation resonances and deal with the singularities in the Born−Oppenheimer Hamiltonian. 1.1. Complex Absorbing Potentials (CAPs). Modification of the original molecular Hamiltonian by adding a CAP eliminates the problem of calculating resonances that arise from the singularity in the Born−Oppenheimer Hamiltonian. In addition, it renders resonance states to the Hilbert space and enables their description using a Gaussian basis set. There is a large variety of CAPs that can be used to calculate resonances. The approach of most CAPs can be divided into two steps in the numerical calculations of the resonance decay poles of the scattering matrix. In the first step of the calculations, a CAP term is added to the Hamiltonian in order to guarantee that the asymptotes of the resonance eigenfunctions decay to zero. Let us take as an example the Zel’dovitz similarity transformation, which can be regarded as the base to the calculations of resonances by complex absorbing potentials (CAPs). The Zel’dovitz similarity transformation is given as

Ψbound/res(r1, ..., rN )→

∑ d j(Φj(r1, ..., rN − 1) j

·exp( +i pj ·rN /ℏ)

(1)

where ( is the antisymmetrizer operator, and N is the number of electrons in the atomic or molecular system under study. For bound states, all components of pj get purely imaginary values, pj/ℏ = +ikbound (the wave vectors kbound get real and j j positive values) and exp( +ipj ·rN /ℏ) = exp(−kbound ·rN ) j

Consequently, as r →∞

∑ dj(Φj(r1, ..., rN − 1)·

Ψbound(r1, ..., rN ) →

j

·rN ) → 0 exp( −kbound j

(2)

The above equation shows that bound states are square integrable functions. The autoionization resonances are associated with eigenfunctions that are not in the Hilbert space and exponentially res res diverge since pres j /ℏ = kj gets complex values where Im[pj ] < res res res res 0 and consequently kj = Re(kj ) + iIm(kj ) = k r,j−ikres i, j . Therefore, Ψres(r1, ..., rN ) →

̂ Ψres(r1, ..., rN ) → limλ → 0+ ∑ d j(Φj(r1, ..., rN − 1) · SZel’dovich j

e

·e

+i k res j ·rN

(4)

→0

Although [ŜZel′dovich Ψres(r1,...,rN)] is a square integrable function, it shows a local exponential divergence behavior in intermediate region in the spatial space. It can be proved (see eq 5.6 in ref 7) that the use of the Zel’dovich similarity transformation is equivalent to the inclusion of an additional term in the original physical Hamiltonian, which is given by

∑ dj(Φj(r1, ..., rN − 1)·exp(+i k rres,j ·rN)· j

exp( +k ires , j · rN ) → ∞

−λ r 2N

(3)

Zel’dovich

̂ −i lim λVCAP

In order to be able to calculate the autoionization resonances by using Gaussians as a basis set, within the BO approximation, we should make a similarity transformation, which bring the resonance exponential diverged functions back to the Hilbert space. For the different types of similarity transformations, Ŝ, which

|λ|→ 0

where Zel’dovich

̂ VCAP

̂ res(r1, ..., rN ) → 0 SΨ

=

ℏ me

N

∑ [ripi ̂

+ pi ̂ ri]

i=1

and me is the electronic mass. In the second step of the calculations, the artificial effect of the CAP on the solutions of the time-independent Schrödinger equation is removed by taking the limit of the strength CAP parameter to zero. This approach can be extended for a large variety of CAPs that are introduced far away from the interaction region, while the strength CAP parameter is reduced as much as possible (due to the use of finite number of basis functions, it is hard to take the limit of λ → 0). See, for example, two methods that discuss the removal of the CAPs after they were used for the calculations of molecular complex resonances.8−10 Although the two methods are similar in their spirit, they are different in the computational algorithms that are used for removing the artificial reflections, which are introduced by the CAPs. For the use of CAP in the calculations of molecular autoionization resonances, see refs 11 and 12. The main purpose of CAPs is that the introduction of CAPs into the atomic/molecular Hamiltonian enables one to use computational algorithms that originally developed for bound states to calculate resonance positions and widths. However, as explained here, CAPs provide accurate resonance positions and widths (associated with the real and imaginary parts of the complex eigenvalues of the CAP augmented Hamiltonian) only in the limit when the CAP strength goes to zero. See, for

where r→∞ see Chapter 5 in ref 7. It should be stressed here that contrary to atomic calculations, when calculating molecules, the similarity transformations noted above should take into account the singularity in the Born−Oppenheimer Hamiltonian. Let us explain this important point in some more detail. The electron−nucleus attractive potential terms are inversely proportional to the distance of any one of the electrons from the nuclei, where within the Born−Oppenhiemer approximation they are held fixed in space. Therefore, it is clear that there are singularities in the electron−nucleus attractive potential terms whenever |ri −Rj| = 0, where ri stands for the ith electronic coordinates in 3D space, and Rj stands for the jth nucleus coordinates. Because of these singularities, the electron−nucleus attractive potential terms are analytical functions of the electronic coordinates only inside a sphere, which has a radius of Rj. Therefore, a uniform complex scaling (CS), i.e., analytical dilation of the coordinates, ri → ri exp(i θ), that suppresses the exponential divergence of the resonance functions is not applicable in calculation of molecules. B


Article


It should be stressed here that two conditions should be satisfied in order to obtain the resonances by the RF-CAPs even when the CAP strength parameter λ = 1. Condition I: The CAP is introduced in the noninteracting region where the electron−nuclei potential terms are almost vanished. Condition II: The electronic repulsion between an electron located out of the interaction region with all other electrons located in the interaction region is negligible. These conditions are required to avoid the singularity of BO Hamiltonian, as discussed above, and keep the Hamiltonian unchanged. For more explanation of these two conditions, see ref 15. In our work, however, we show that excellent results for the resonance positions and widths can be obtained when these two conditions are not satisfied. This can be achieved when the dependence of the resonance complex eigenvalues on λ is removed by taking the limit of the complex eigenvalues by Padé approximant to λ → 0. 1.3. Uniform Complex Basis Functions (CBFs). The CAPs and RF-CAPs presented above enable the calculations of molecular autoionization resonances. They avoid the singularity in the Born−Oppenheimer Hamiltonian by imposing absorbing boundary conditions in the region in space where the nondilation electron−nucleus attractive potential terms vanish. Alternatively, it is possible to use complex basis functions which are associated with uniform complex scaling of all Gaussian basis functions in the calculation of molecular resonances. In ref 22, it was shown that although the molecular Hamiltonian within the framework of the Born−Oppenheimer approximation is not dilation analytic (as explained above), the Hamiltonian matrix elements as obtained by using the Gaussian basis sets (as in the standard electron structure codes) are analytical functions of the scaling parameter. Therefore, instead of carrying out the uniform complex scaling ri → ri η, where η = αexp(i θ), one can carry out analytical continuation of the Hamiltonian matrix elements, where the Gaussian exponential parameters are scaled by 1/η2 = (α−1exp(−iθ))2. Namely, complex basis functions (CBFs) are used in the variational calculations. Here, all Gaussian basis functions were uniformly scaled by the complex factor 1/η. The molecular autoionization resonances are associated with the stationary solutions for which the first order derivatives of the complex eigenvalues of the non-hermitian symmetric Hamiltonian matrix are vanished. See ref 22 for the first use of this type of CBFs for the calculations of the shape type resonance of H−2 and of the Feshbach autoionization of doubly excited hydrogen molecule. In our calculations presented below, we show that indeed using CBFs is equivalent to using a contour integration in the complex plane, where all the grid points are lying on a straight line, which is rotated into the complex coordinate plane by the angle θ. 1.4. Partial CBFs: Scaling Only the Diffuse Basis Functions. As explained in the previous subsection, one can avoid the singularity in the Born−Oppenheimer Hamiltonian by carrying out the analytical continuation of the Hamiltonian matrix elements rather than the operator.22 Originally, the analytical continuation of the Hamilotnian matrix elements were carried out by uniformly scaling all Gaussian basis functions that were used in the calculations of the molecular Hamiltonian matrix elements. However, as explained previously in ref 19, there are an infinite number of analytical continuations that can be carried out for calculating resonances. All of them scale the asymptotes of the resonance states by η = αexp(iθ). This is equivalent to the method which has previously

example, the energies and widths (inverse lifetimes) of metastable molecular ion states as calculated by introducing CAPs into equation of motion electron attachment coupled cluster (EOM-EA-CC) within electronic structure packages.10,13,14 1.2. RF-CAPs: Reflection-Free Complex Absorbing Potentials. Unlike the CAPs described above that provide accurate resonance positions and widths in the limit as the CAP strength parameter goes to zero, the reflection-free CAP (RFCAP) provides accurate results also when the λ ≠ 0 and even when the RF-CAP does not vanish in the interaction region. The RF-CAPs are universal and are calculated from a complex contour of integration when the Hamiltonian matrix elements are calculated.15 The imaginary part of the contour of integration F(r) is as close as one wishes to zero in the interaction region (not necessarily covering the entire interaction region) and outside of the integration region F(r) → r exp(i θ) . When the contour of integration F(r) is a smooth (analytical) function of the electronic coordinates, it is associated with the smooth exterior scaling (SES) transformation.15 We can summarize it by saying that SES stands for the transformation that enables the derivation of the RF-CAPs, such that ̂ Ψres(r1, ..., rN ) → SSES

∑ dj(Φj(r1, ..., rN − 1)·e+i(e

iθ res k j )·rN

→0

j

(5)

when θ is larger than the phase angle of However, the contour of integration F(r) can be a piecewise potential that gets a zero value in the interaction region (to avoid the singularity points in the Coulombic potential terms in the molecular Hamiltonian) and beyond a critical point in the coordinate space F(r) = r exp(i θ) . In such a case, instead of the smooth exterior scaling transformation, one can use the exterior complex scaling (ECS) transforation ŜECS to obtain the converged asymptotes of the resonance wave functions.16−18 See ref 19 for describing all different type of complex scaling transformations which are associated with a large family of complex contour of integrations in the calculation of the Hamiltonian matrix elements that make the resonance solutions be square integrable. As it was proved in ref 15, the SES transformation is equivalent to the inclusion of the RF-CAP that has the universal functional behavior given by kres j .

̂ ‐CAP = λ[V0(r) + V1(r)∂r + V2(r)∂ 2r ] VRF

where the RF-CAP strength parameter, λ, is equal to unity (and not λ → 0 as in the standard CAP described in the previous subsection). Here, V0,V1, and V2 are complex potentials which vanish at the interaction region where the imaginary part of F(r) is vanished (not in the mathematical sense but in all significant digits which depend on the accuracy of the numerical calculations). These complex potentials depends mainly on three characteristic parameters of the complex contour of integration F(r): (1) the parameter that controls the region in space where the contour F(r) departs from the real axis and penetrates into the complex coordinate plane, (2) the parameter that controls the “rate” of the penetration of the contour into the complex plane, and (3) the parameter of the rotational angle θ that makes the asymptotes of the resonance (metastable) states to vanish as shown in eq 5. In refs 20 and 21, an RF-CAP was obtained from a SES contour to calculate atomic and molecular autoionization Feshbach type resonances. C


Article

Journal of Chemical Theory and Computation been suggested by Rescigno and McCurdy.23 The idea is only to scale the most diffused basis functions by a complex factor. The scaling of the basis functions by a complex factor is equivalent to the rotation of the electronic coordinates of the diffused basis functions to the complex coordinate plane. Therefore, the use of CBFs is somehow similar to the use of ECS/SES contour of integrations.15−18,24 Recently, by taking the CBFs approach, shape-type resonances of molecular anions were calculated.25 The advantages of calculating resonances by only complex scaling the most diffused basis functions were shown for Feshbach autoionization atomic resonances.26 This work compares the different methods of computing metastable states (resonances) using Gaussian functions. The subject is of current interest as methods of computing resonances using ab initio methodologies have accumulated, and it is relevant to shed light on their differences and similarities. Since Gaussian functions are usually used in quantum chemistry, their use in the present work is particularly relevant. It is clear from the above introduction that the different methods that were developed for the calculations of molecular resonances (and thereby CPES) are closely related one to another. However, it is not absolutely clear how the CBFs method is related to the calculations of resonances by CAPs. Since these two methods are the most promising approaches for calculating resonances and CPES for large molecular systems, we focus in this work on the study of the relationship between these methods. In this work, we show under which conditions the use of CBFs for calculating molecular autoionization resonances is equivalent to the use of CAPs for calculating the molecular resonances. The paper is organized as follows. In the next section, we show how the complex SES contour of integrations can be computed from complex Gaussian basis functions. The most diffuse Gaussians are complex scaled to keep the complex contour of integration close to the real axis in the interaction region. In the third section of the paper, we show how by using complex SES contours that were obtained from the use of partial CBFs reflection-free CAPs (RF-CAPs) are derived. In the fourth section, an illustrative numerical example is given where we show how the resonances are calculated by the use of complex Gaussian basis set and how they are calculated by RFCAPs that were derived from the nonuniform CBFs. The illustrative numerical example is for a single-particle problem. However, one can learn much from the latter with a manyparticle quantum chemistry study. Concluding remarks are given in the last section emphasizing possible future use of our findings in the calculations of CPES that are needed, for example, in the study of cold molecular collisions where ionization and dissociation take place simultaneously.

When the Gaussian basis sets are used, the matrix X, which presents the operator x̂, is obtained. The corresponded transformed matrix is given by

? = S−1/2XS−1/2

Similarly, @ and A are calculated. For the sake of simplicity and without loss of generality, we concentrate here on the derivation of the complex contour of integration from the xdependent complex Gaussian basis functions in single particle problem. In a very similar way, the other complex contours in the y and z coordinates can be calculated. The contour of integration in the spatial x space is associated with the grid points {xj}j=1,2,.. that are associated with the eigenvalues of ?

?T(jx) = xj T(jx)

(8)

The above calculation is also the first step in the discrete variable representation (DVR) method. In the DVR method, the Hamiltonian matrix is given by [T(x)]t HT(x) = [T(x)]t KT(x) + V diag

(9)

where Vdiag is a diagonal matrix where the jth component on the diagonal is V(xj) and T(x) are the values of the different j orthonormal basis functions at the grid point xj multiplied by the weight factor of each basis function. The diffused Gaussian basis functions {xn exp(−αjx2) }n=0,1,2···,j=0,1,...,jmax, where α0>α1>α2... can be analytical, continued into the complex plane such that for αj < αth, αj → αj exp(−i2 θ) . Consequently, {xj}j=1,2,.. are grid points on the complex contour of integration. For the even-tempered Gaussians where αj = α0ϵj−1, the complex integration contours are presented in Figures 1 and 2 for θ = 0.25. The contour in Figure 1 was

Figure 1. Discrete complex integration contour obtained by CBFs is depicted. The basis set used is constructed of even-tempered Gaussians, where the diffused functions are analytically continued into the complex plane. As one can see the complex integration contour behaves like SES. The contour was obtained for jmax= 41, α0 = 1000, ϵ0 = 1.4125, αth = 0.1,n = 0,1, and θ = 0.25.

2. FROM COMPLEX BASIS FUNCTIONS (CBF) TO SMOOTH EXTERIOR COMPLEX SCALING (SES) CONTOURS 2.1. From CBF to Grid Points on Complex Contours of Integration. When Gaussian basis sets are used, the Hamiltonian matrix elements are given by H = K +V, where K and V are, respectively, the kinetic Hamiltonian and potential Hamiltonian matrices. The overlapping matrix is S. The spectrum of the molecular Hamiltonian is obtained by the diagonalization of the transformed Hamiltonian matrix / = S−1/2HS−1/2 = 2 + =

(7)

obtained for jmax= 41, α0 = 1000, ϵ0 = 1.4125, αth = 0.1, andn = 0,1. The contours in Figure 2 were obtained for jmax= 21, α0 = 1000, ϵ0 = 1.4125, and αth = 0.1, with a different value of n for each contour. As one can see, different sets of CBF generate different complex integration contours, which resembles the SES contour behavior as presented in ref 15. It should be noted that when αth = α0, all the basis functions are scaled, and the obtained contour is like uniform complex scaling, as shown in Figure 3.

(6) D


Article


Figure 2. Each discrete contour depicts the eigenvalues of X (in atomic units) for complex diffused Gaussians basis functions with different values of n. In the legend, one can see the match between contour and the Gaussians basis we used. All the contours were obtained for jmax = 21, α0 = 1000, ϵ0 = 1.4125, αth = 0.1, and θ = 0.25.

Figure 4. (a) The local rotation (in radians) of every grid point (eigenvalues of X in atomic units) for complex Gaussian basis into the complex plane is depicted for different values of αth, where n = 0,1, jmax = 41, α0 = 1000, and ϵ0 = 1.4125. (b) Zooms into the regions where Im(X) ≃ 0 (2X0) for each value of αth. As clearly shown, 2X0 (unscaled region) is larger as αth is smaller.

Figure 3. Depicted discrete contour was obtained for jmax= 41, α0 = 1000, ϵ0 = 1.4125,n = 0,1, and θ = 0.25, while αth = α0, such that all the basis functions are scaled. As shown, the complex scaling of all basis functions is equal to uniform complex scaling.

/(xj ′, y

, z ),(xj , yk , zq) k′ q′

The local rotation of every grid point into the complex coordinate plane is defined by Imxj θjlocal = arctan Rexj (10)

= [2 (x)]xj ′, xj δy

, y δzq ′ , zq

k′ k

+ [2 (z)]zq′, zq δxj ′, xjδy

+ [2 (y)]y

, y δxj ′ , xjδzq ′ , zq

k′ k

,y

k′ k

+ V (xj , yk , zq)δxj ′, xjδy

, y δzq ′ , zq

k′ k

(11)

When CBFs are used to generate a SES-like contour as described above such that the potential matrix elements remain sufficiently close to the real axis, then the imaginary part of the Hamiltonian matrix elements is given by

In Figure 4(a), we show the rotation angle of every grid point xj into the complex plane. The grid points that remain close to the real axis determine the region in space where V(xj) = V(Rexj). In the calculations of molecular resonances, we would like for this to region include the nuclei positions that stay unscaled since the singularity of the BO Hamiltonian is discussed above. As is shown in Figure 4, the unscaled region can be controlled by varying the value of αth. Figure 4(b) [zoom in of Figure 4(a)] demonstrates that we can determine the region − x0 ≤ x ≤ + x0, where the imaginary part of the SES contour obtained from CBF calculations remains close to the real axis, as close as we wish. For example, in Figure 4(b), a limit of 8 × 103 was chosen to determine the desired region. Similarly, the grid points on the complex integration contours {yk}k=1,2,.. and {zq}q=1,2,.. and their corresponding eigenvector matrices ; (y) and ; (z) can be calculated. The one-particle Hamiltonian matrix elements for the 3D case are given by

Im[/(xj ′, y

]

, z ),(xj , yk , zq) k′ q′

= Im[2 (x)]xj ′, xj δy

, y δzq ′ , zq

k′ k

+ Im[2 (y)]y

, y δxj ′ , xjδzq ′ , zq

k′ k

(z )

+ Im[2 ]zq′, zq δxj ′, xjδy

,y

k′ k

(12)

In the next subsection, we show how analytical expressions for the complex contours are obtained. These analytical expressions are used in the next section for the calculations of the RF-CAP. It should be noted that as mentioned above the potential matrix elements are not absolutely real and only remain close to the real axis because the complex contour are smooth. Therefore, as is discussed in the Introduction, molecular resonance calculations using the RF-CAP obtained by this contour, when λ = 1, will obtain less accurate results. E


Article

Journal of Chemical Theory and Computation 2.2. From Grid Points on the Complex Contour of Integration to Analytical Expressions. There are many different methods that can be used to fit the numerical complex grid point xj, which were obtained from the CBF, to the analytical formula. Here, we focus on two different approaches. The first one is the fit of the numerical complex grid point xj to eq 18 in ref 15. In this approach, there are three parameters that are optimized to provide the best fit of the analytical formula to the numerical complex grid point xj. One of these parameters is θ that determines the rotation of the contour to the complex plane when |x| > x0. The parameter x0 determines the region − x0 ≤ x ≤ + x0, where the contour remains as close as possible to the real axis. The third fitting parameter is λ, which determines the “rate” that the contour is detached from the real axis and gets to the complex plane with the rotational angle θ. Exterior scaling transformation16−18 is obtained when λ →∞. See Figure 5 for the results which were obtained by the fitting of eq 18 in ref 15 to the grid points obtained as described

formula for the contour in the complex coordinate plane can play a major role in our derivation of the RF-CAP as is described in the next session. Let us explain more on every one of the fitting methods we have used here. The first method is fitting the grid points on the complex integration contour to analytical expression. We have chosen eq 18 in ref 15 for more convenient reasons. The RF-CAPs are derived for any type of SES transformation. The analytical function we have borrowed from ref 15 for the fitting procedure is given by ⎛ cosh(λ(ξ − x0)) ⎞ 1 F(ξ) = ξ + (eiθ − 1)⎜ξ + ln ⎟ 2λ cosh(λ(ξ + x0)) ⎠ ⎝ (13)

In order to find the values of the parameters in the above expression, we fit Re(xj) + iIm (xj) = F(ξj), when ξj = |Re(xj) + iIm(xj)|. This expression for ξ is an approximation of the unscaled xj because inside the region where the imaginary value of xj is very small, ξ is approximately equal Re(x), and when ξ →∞, ξ →ξ exp(i θ) . The results presented in Figure 5 were obtained for λ = 0.7046, x0 = 4.367, and θ = 0.2682. For the fitting procedure using the Padé approximation, we have used a simplified version of Padé. The Schlessinger truncated continued fraction has the form27 y1 y(x) ≡ CM(x) = a1(x − x1) 1+ a2(x − x 2) 1 + ⋮a

M (x − xM )

(14)

where the aj coefficients are chosen such that CM(xj) = yj , i = 1, 2, ..., M

(15)

In the same way as the first fitting, in the Padé approximation fitting, we use ξ = |Re(x) + iIm(x)| and use Padé to get the two sets of coefficients such that CM(ξ) = Re(x) + iIm(x) . Once the aM coefficients are determined, we perform an analytical continuation into the complex plane by evaluating CM(ξ). Thus, we used the Padé approximation method to transform real grid points to complex grid points. Convergence of the extrapolated function, CM(x), with respect to M is routinely checked, and the difference between CM(x) and CM−1(x) is reported as the Padé error. The results presented in Figure 6 are obtained for M = 28.

Figure 5. Complex contour obtained by fitting the grid points (obtained from the diagonalization of the operator x̂ presented by CBFs) to the analytical expression given in eq 13. The parameters that yields the best fit are λ = 0.7046, x0 = 4.367, and θ = 0.2682.

above and are presented in Figure 1. In the second fitting procedure, we applied the Padé approximation using the Schlessinger method.27 The results of this fitting are presented in Figure 6. The fact that the above fittings provide an analytical

3. FROM SES COMPLEX CONTOURS OBTAINED BY USING CBFS TO REFLECTION-FREE CAPS The analytical expressions derived in the previous section, which provide complex integration contours based on CBF calculations, are denoted by Fq(x). The contour obtained by eq 18 in ref 15 is associated with q = 1, and Fq=2(x) is associated with the use of Padé as a fitting procedure. As one can see from the derivation of the RF-CAP given in ref 15, the first-, second-, and third-order derivatives of F(x) with respect of x are used for the construction of the RF-CAP from Fq(x). Using these derivatives, three complex potentials are obtained, V0(x), V1(x), andV2(x), by using eqs 9, 10, and 16 in ref 15. The RF-CAP is given by ̂ ‐CAP = V0(x) + V1(x) VRF

Figure 6. Complex contour obtained by fitting the grid points (obtained form the diagonalization of the operator x̂ presented by CBFs) to the analytical expression given in eq 14, i.e., using the Padé approximation (with M = 28).

d d2 + V2(x) 2 dx dx

(16)

The Hamiltonian matrix elements when using the RF-CAP are given by F


Article

Journal of Chemical Theory and Computation ̂ ‐CAP + V c H = K + VRF

(17)

c

when V = V(F(x)) is the complex potential Hamiltonian matrix, which is calculated using the complex contour for numeric integration, and the complex kinetic Hamiltonian matrix is equal to K + V̂ RF‑CAP. By choosing V = Vc, the potential Hamiltonian matrix can be analytically calculated without the need to scale the potential energy term in the Hamiltonian. The requirement of V = Vc implies that one should use the ECS or to SES if and only if the ECS and SES contours penetrate into the complex coordinate plane only at the noninteracting regions. As a matter of fact, it is hard to satisfy this condition when long-range potentials (as the Coulombic potentials) are used. In such cases, the use of an unscaled potential term in the Hamiltonian is an approximation that introduces inaccuracies in the calculations of the resonance positions and widths. The degree of the inaccuracies depends on the selection of the shape of the SES (or ECS) contour that is used in the construction of the RF-CAP. Here, we show that the SES contour, which is obtained from the use of the complex Gaussian basis, does not satisfy the desired condition of V = Vc even in the case where short-range potentials are used. See Figure 7 for V0(x), V1(x), and V2(x), which were obtained from formula F1, and see Figure 8 for V0(x), V1(x), and V2(x), which were obtained from the second fitting procedure F2. As one can see, although every one of these complex local potential terms is generated by different formula, the behavior of each potential is similar since they all describe the same contour. Nonetheless, there are some differences between the complex potentials that are generated from F1 and from F2. The reasons for these differences are the inaccuracies in the fitting, especially in F1, presented in Figure 5, where we impose a special formula on the contour. Therefore, we will use the RF-CAP, which was derived from F2 , to illustrate the RFCAP calculations.

4. ILLUSTRATIVE NUMERICAL EXAMPLE The 1D model Hamiltonian that we use to illustrate the CBF calculation and RF-CAP calculations is given by 2

Ĥ =

Px̂ + V (x ) 2

(18)

where ⎞ ⎛ x2 2 V (x) = ⎜ − 0.8⎟e−0.1x ⎠ ⎝2

(19)

This symmetrical double barrier potential has been used before as a test case for new theories and computational algorithms for calculating resonances. The poles of the scattering matrix in this case consist of one bound state and several narrow and infinite large number of broad and overlapping resonances. We focus here on the narrowest resonances. The exact values for the resonance positions and widths are known from the literature. However, we have revaluated them by using a large number of particles in a box basis function to diagonalize the complexscaled Hamiltonian. Converged results to any desired accuracy for the resonance positions and widths were obtained by increasing the number of the particles in the box basis functions and by increasing the size of the symmetric box. Two methods are used here to calculate the resonances positions and widths. In the first method, CBF, which consists of complex eventempered Gaussian functions, are used for the calculations of

Figure 7. V0(x), V1(x), and V2(x) obtained from the complex contour in Figure 5, which results from fitting the CBF contour to eq 13 (eq 18 in ref 15).

the resonance positions and widths. In Figure 1, we represent the complex contour that is obtained when the Gaussian CBFs were used to represent the position operator (see a detailed explanation above). In the second method, we use the RF-CAP, which has been derived from the complex integration contour given in Figure 1. One may expect that the two methods shall yield the same results since as we have shown above they are related to one another by a similarity transformation. However, it is clear that the RF-CAP-based solutions are much less accurate than the CBF-based solution mainly since the G


Article


the non-hermitian sector of the domain of the Hamiltonian, while the Padé approach enables us to get rid of the artificial effects introduced by the CAP due to a series of approximations as mentioned above (one of the approximations is using the CAP without taking into consideration that the RF-CAP derived from the Gaussian CBFs is pronounced when the potential is not exponentially small). 4.1. Complex Basis Functions (CBF)-Based Calculations. The basis set used in this method consists of diffuse Gaussians {xn exp(−αjx2)}n=0,1,2···,j=0,1,...,jmax, where α0 > α1 > α2... are analytically continued into the complex plane such that for αj < αth, αj → αj exp(−i2 θ), and the c-product is used for the calculations of the Hamiltonian matrix elements. The basis set used in the calculations is even-tempered Gaussians αj = α0ϵj−1, when jmax= 41, α0 = 1000, ϵ0 = 1.4125, and n = 0,1. The scaling parameter, θ = 0.25, and αth = 0.1 are held fixed. Since each one of the resonances becomes a stationary solution in the complex variational space with different value of θ, it is necessary, in principle, to evaluate an integration counter for each state separately. However, it is interesting to see that even when we use only a single general integration counter with a fixed value of θ (without carrying out θ and αth trajectory calculations) it is possible to obtain the entire complex spectrum of the model Hamiltonian. Since our objective is not to accurately calculate the resonance energies using the CBF approach but rather to generate a RF-CAP, it is sufficient and desirable to use one complex integration contour for calculating all states. In addition, note that for a complete basis set the resonance positions, Er−Eth, and widths, Γ, should get the exact values for any value of θ for which tan θ > Γ/(2(Er−Eth)) and regardless of the value of αth. The threshold energy in our case is vanished, Eth = 0. Figure 9 presents the spectrum of the eigenvalues of the Hamiltonian given in eq 18, which are associated with the bound, resonances, and rotating continuum. As shown from the

Figure 9. The bound state and the resonances calculated by the use of complex Gaussian basis functions (CBF) are shown adjacent to the exact values of these states. The exponent parameters of the diffuse even-tempered Gaussians were scaled by exp(−2iθ), where θ = 0.25 rad. The diffuse functions are defined as those for which α > αth where αth = 0.1. Note that each one of the resonance stationary solutions requires a different value of θ (obtained by carrying a θ trajectory). However, it is interesting to see that even when we use only a fixed value of θ it is possible to obtain the entire complex spectrum of the model Hamiltonian. Since our objective is not to accurately calculate the resonance energies using the CBF approach but rather to generate a RF-CAP, it is sufficient and desirable to use one complex integration contour for calculating all states.

Figure 8. V0(x), V1(x), and V2(x) obtained from the complex contour in Figure 6, which results from fitting the CBF contour using the Padé approximation.

amplitude of the RF-CAP derived from the CBFs is not equal to zero in the interaction region, where V(x) is pronounced (i.e., does not get exponentially small values), and therefore, the use of real potential in the calculations rather than complexscaled potential introduces numerical errors. This problem will be eliminated later by removing the CAP by using the Padé approach. The idea is that CAPs enable us to “penetrate” into H


Article


1.4125, and n = 0,1. The spectrum of our model Hamiltonian, which was calculated numerically by introducing the RF-CAP as derived from the CBFs, is presented in Figure 11. As we

results presented in Figure 9, the continuum eigenvalues are rotated into the complex energy plane, and the bound state and the resonances remain stable at the complex plane. By using different values of αth or θ, one can control the number and the numerical accuracy of the obtained resonances. In Figure 10,

Figure 11. The bound and resonances energies obtained by the RFCAP approach (using eq 20 and VRF‑CAP from Figure 8), as well as the exact values.

expected, the resonances and the bound state obtained are not as accurate as those obtained from the use of complex Gaussian basis functions. The reason for this is that the SES contour, which has been obtained from the representation of the position operator x̂ by complex Gaussian basis functions, is not equal to zero in the interaction region, meaning V ≠ Vc (see for example the shape of the CAP, which was used in ref 10 in the calculations of shape-type resonances of molecular anions where CAP gets nonzero values in the region where the Coulombic potential is almost vanished.) In order to find more accurate values, we removed the RFCAP. A similar approach was taken in ref 10. However, our way of removing the CAP is different. Rather than using perturbation theory, we use the special analytical properties of the Padé approximate. The Padé approximate enables one to carry out analytical continuation of the function under study far away from the available complex grid points. Let us describe the method we have used that was adopted from the approach described in ref 8. First, we found the bound state and the resonances which were obtained for the modified Hamiltonian, which is given by

Figure 10. (a) The θ trajectory of the second (narrowest) resonance state. (b) Zooms into the cusp region of (a).

the narrowest resonance is presented as a function of θ. One can see that the value of θ = 0.25 is far from the optimal angle for which a stationary solution in the complex variational space dE is obtained. Namely, at the angle for which dθ = 0. Since the optimal angle gets different values for different resonances, we have chosen to represent all results for the different resonance positions and widths for the same rotational angle of θ = 0.25. 4.2. RF-CAP-Based Calculations. In this method, the Hamiltonian is given by Ĥ =

2 ⎞ ⎛ x2 2 Px̂ ̂ ‐CAP + ⎜ − 0.8⎟e−0.1x + VRF 2 ⎠ ⎝2

Ĥ (λ) =

2 ⎛ x2 ⎞ 2 Px̂ ̂ ‐CAP + ⎜ − 0.8⎟e−0.1x + λVRF 2 ⎝2 ⎠

(21)

for different discrete values of λj = Δ × (j + j0), where j = 1,...M. The values for λ are chosen for each resonance or bound state due to the fact that for the chosen λ there are constant differences between the resonance values (or the bound state values) which are obtained by consecutive λ. The parameters of λ for each resonance or bound state are bound state − Δ = 0.05, j0 = 0.1, M = 9; the first resonance − Δ = 0.1, j0 = 1.4, M = 11; the second resonance − Δ = 0.1, j0 = 1.8, M = 11; the third resonance − Δ = 0.1, j0 = 3.4, M = 11. Then, we used the Padé method, as described above, to get expressions for the bound state and the resonances values such that y(x) as eq 14 is the complex energy as a function of λ. Namely, yj = Ej, where Ej is the resonance complex eigenvalue that was obtained when the CAP strength parameter is equal to λ = λj where j = 1,...M.

(20)

where V̂ RF‑CAP is as it is given in eq 16 with the V0(x), V1(x), and V2(x) local universal potentials (i.e., do not depend on the problem under study), which are shown in Figure 8. We should re-emphasize here that these universal potentials depend only on our chosen complex-scaled even-tempered Gaussians as explained above. The basis functions set we used is eventempered diffuse Gaussians {xn exp(−αjx2) }n=0,1,2···,j=0,1,...,jmax, where α0 > α1 > α2..., αj = α0ϵj−1, jmax= 41, α0 = 1000, ϵ0 = I


Article

Journal of Chemical Theory and Computation When the a1M coefficients in eq 14 are determined, we can substitute λ = 0 and get the accurate values for the resonance positions and widths. See Figure 12 for an example of the

one used within the smooth exterior scaling (SES) transformation.15,19 It should be mentioned that these two propositions to solve the problem of the nonanalyticity of the Born−Openheimer Hamiltonian where shown by Simon to be based on rigorous mathematical grounds.24 Yet, it is the first time that this claim has been supported by numerical results. Moreover, our numerical calculations show that by using a complex Gaussian basis set a SES contour is obtained, FCBF(x). Subsequently, we show how one can obtain analytical expressions for the SES contour of integration from the grid of points obtained from the CBF calculations. Finally, we show how by using this SES contour of integration we can derive RFCAP. The main message of the work is that, in principle, the CBF, SES, and RF-CAPs methods for calculating resonances are all equivalent to one another. Nevertheless, as explained in this paper, it is expected to obtain more accurate results by the use of CBF rather than RF-CAPs due to the fact that the potential used is approximated V(x) (see ref 15 for the need to modify the potential when the CAP does not get exponential small values in the interacting region). The reason for this is that the SES contour does not remain on the real axis in the interaction region, although it is very close to the real axis. Consequently, the use of unscaled potential terms in Hamiltonian when using RF-CAP obtained by this contour is not justified since V(F(x)) ≠ V(x).15 However, it simplifies the numerical calculations significantly. It is important to note that the use of RF-CAPs or other type of CAPs is simpler in the electronic structure calculations than the use of CBFs, which requires more severe modifications in the codes that are used for the calculations of molecular potential energy surfaces.28 Furthermore, our calculations show that a very promising approach in calculations of resonances is to use CAPs and then remove their artificial effect (see refs 8 and 10 for different approaches that remove CAPs). A comparison between the CBF and removing-CAP methods is out of the scope of this work; however, it is hard to believe that one method will be found superior to the other one in all cases.

Figure 12. Padé removal of the artificial effects of the RF-CAP potential in the calculation of the second (narrowest) resonance energy. The energies given by the RF-CAP appraoch (using eq 21 and VRF‑CAP from Figure 8) for different values of λ in black. The RF-CAP energies in green are used within Padé in order to eliminate the RFCAP and provide the final result at λ = 0.

results for the position and width of the second narrowest resonance by removing RF-CAP by the use of Padé. It is quite amazing to see that for λ = 0 in the Padé approximate an accurate resonance position and width were obtained in spite of the fact that direct diagonalization of the Hamiltonian when λ = 0 gives real eigenvalues only (i.e., zero resonance widths). A summary of the results for the resonance positions and widths as obtained from the different approaches (which all were proved to be equivalent if no numerical approximations are taken) is presented in Table 1. The different approaches which were taken are the CBF method and the RF-CAP method, before and after removing the CAP by Padé.

■

AUTHOR INFORMATION

Corresponding Author

5. CONCLUDING REMARKS AND FUTURE APPLICATIONS FOR CALCULATIONS OF COMPLEX MOLECULAR POTENTIAL ENERGY SURFACES This paper demonstrates how different methods for calculating molecular resonance are related one to another. We have shown how the use of CBF, as first proposed by McCurdy and Rescigno23 (and closely related to Moiseyev and Corcoran approach22), provides a complex integration contour for the Hamiltonian matrix elements. This contour is comparable to the

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the I-Core: the Israeli Excellence Center “Circle of Light”, by the Israel Science Foundation Grants No. 298/11 and No. 1530/15. A.B.A.

Table 1. Comparison of the Bound and the Lowest Three Resonance Energies of the 1D Model Hamiltonian Obtained by Different Computational Methodsa first resonance

bound state exact values SES-CBFb RF-CAP(based CBF) removed-RF-CAP(Padé)

second resonance

third resonance

ReE

ImE

ReE

ImE

ReE

ImE

ReE

ImE

−0.298 −0.298 −0.298 −0.298

0 −1.112 × 10−14 −0.00609 −1.05 × 10−13

0.621 0.621 0.6204 0.621

−5.827 × 10−5 −5.834 × 10−5 −0.0102 −5.844 × 10−5

1.327 1.327 1.322 1.327

−0.01545 −0.01528 −0.0304 −0.01529

1.7846 1.7968 1.8 1.7860

−0.1738 −0.1645 −0.1712 −0.1733

a

All calculations use a Gaussian basis set. The computational methods are SES-CBF that uses of CBFs which form a SES contour of integration; RFCAP that uses a RF-CAP based on the SES contour of integration which was obtained by the CBFs calculation. Finally, we present results which remove the effect of this RF-CAP using the Padé approximation. bThese values are obtained with αth = 0.1 and θ = 0.25 and not with the optimal scaling parameters J


Article

Journal of Chemical Theory and Computation acknowledges Dr. Ido Gilari, Mr. Idan Haritan, Dr. Debarati Bhattacharya, Dr. Arik Landau, and Ms. Hanna Martiskainen for most helpful discussions and assistance in entering a new field of research.

■

REFERENCES

(1) Henson, A. B.; Gersten, S.; Shagam, Y.; Narevicius, J.; Narevicius, E. Science 2012, 338, 234−238. (2) Lavert-Ofir, E.; Shagam, Y.; Henson, A. B.; Gersten, S.; Kłos, J.; Ż uchowski, P. S.; Narevicius, J.; Narevicius, E. Nat. Chem. 2014, 6, 332−335. (3) Shagam, Y.; Klein, A.; Skomorowski, W.; Yun, R.; Averbukh, V.; Koch, C. P.; Narevicius, E. Nat. Chem. 2015, 7, 921. (4) Trinter, F.; Schöffler, M.; Kim, H.-K.; Sturm, F.; Cole, K.; Neumann, N.; Vredenborg, A.; Williams, J.; Bocharova, I.; Guillemin, R.; Simon, M.; Belkacem, A.; Landers, A. L.; Weber, T.; SchmidtBöcking, H.; Dörner, R.; Jahnke, T. Nature 2014, 505, 664−666. (5) Gokhberg, K.; Kolorenč, P.; Kuleff, A. I.; Cederbaum, L. S. Nature 2014, 505, 661−663. (6) Scheit, S.; Averbukh, V.; Meyer, H.; Moiseyev, N.; Santra, R.; Sommerfeld, T.; Zobeley, J.; Cederbaum, L. J. Chem. Phys. 2004, 121, 8393. (7) Moiseyev, N. Non-Hermitian Quantum Mechanics; Cambridge University Press; Cambridge, 2011. (8) Lefebvre, R.; Sindelka, M.; Moiseyev, N. Phys. Rev. A: At., Mol., Opt. Phys. 2005, 72, 052704. (9) Sajeev, Y.; Vysotskiy, V.; Cederbaum, L.; Moiseyev, N. J. Chem. Phys. 2009, 131, 211102. (10) Jagau, T.-C.; Zuev, D.; Bravaya, K. B.; Epifanovsky, E.; Krylov, A. I. J. Phys. Chem. Lett. 2013, 5, 310−315. (11) Santra, R.; Cederbaum, L. S. Phys. Rep. 2002, 368, 1−117. (12) Muga, J. G.; Palao, J.; Navarro, B.; Egusquiza, I. Phys. Rep. 2004, 395, 357−426. (13) Jagau, T.-C.; Krylov, A. I. J. Phys. Chem. Lett. 2014, 5, 3078− 3085. (14) Zuev, D.; Jagau, T.-C.; Bravaya, K. B.; Epifanovsky, E.; Shao, Y.; Sundstrom, E.; Head-Gordon, M.; Krylov, A. I. J. Chem. Phys. 2014, 141, 024102. (15) Moiseyev, N. J. Phys. B: At., Mol. Opt. Phys. 1998, 31, 1431. (16) Simon, B. Phys. Lett. A 1979, 71, 211−214. (17) Nicolaides, C. A.; Beck, D. R. Phys. Lett. A 1978, 65, 11−12. (18) Gyarmati, B.; Vertse, T. Nucl. Phys. A 1971, 160, 523−528. (19) Moiseyev, N.; Hirschfelder, J. J. Chem. Phys. 1988, 88, 1063− 1065. (20) Sajeev, Y.; Sindelka, M.; Moiseyev, N. Chem. Phys. 2006, 329, 307−312. (21) Sajeev, Y.; Moiseyev, N. J. Chem. Phys. 2007, 127, 034105. (22) Moiseyev, N.; Corcoran, C. Phys. Rev. A: At., Mol., Opt. Phys. 1979, 20, 814−817. (23) McCurdy, C., Jr; Rescigno, T. Phys. Rev. Lett. 1978, 41, 1364. (24) Morgan, J.; Simon, B. J. Phys. B: At. Mol. Phys. 1981, 14, L167. (25) White, A. F.; Head-Gordon, M.; McCurdy, C. W. J. Chem. Phys. 2015, 142, 054103−1−11. (26) Landau, A.; Haritan, I.; Kaprálová-Ž ďánská, P. R.; Moiseyev, N. Mol. Phys. 2015, 113, 3141−3146. (27) Schlessinger, L. Phys. Rev. 1968, 167, 1411. (28) White, A. F.; McCurdy, C. W.; Head-Gordon, M. J. Chem. Phys. 2015, 143, 074103.

K


On the Equivalence of Different Methods for ... - ACS Publications

Recommend Documents