Article pubs.acs.org/JCTC
Resonance Energies and Lifetimes from the Analytic Continuation of the Coupling Constant Method: Robust Algorithms and a Critical Analysis Thomas Sommerfeld,*,† Joshua B. Melugin,† Prakash Hamal,† and Masahiro Ehara*,‡ †
Department of Chemistry and Physics, Southeastern Louisiana University, SLU 10878, Hammond, Louisiana 70402, United States Institute for Molecular Science, Research Center for Computational Science, Myodai-ji, Okazaki 444-8585, Japan
‡
ABSTRACT: The energy of a metastable state can be computed by adding an artificial stabilizing potential to the Hamiltonian, increasing the stabilization until the metastable state is turned into a bound one, and then further increasing the stabilization until enough bound-state data have been collected so that these can be extrapolated back to vanishing stabilization. The lifetime of the metastable state can be obtained from the same data, but only if the extrapolation is performed by analytic continuation. This extrapolation method is called analytic continuation of the coupling constant (ACCC). Here we introduce preconditioning schemes for two of the three established extrapolation algorithms and critically compare results from all three extrapolation schemes in a variety of situations: As examples for resonance states serve the π* temporary anions of ethylene and formaldehyde as well as a model potential, which provides a case where input data with full numeric precision are available. In the data collection step, three different stabilizing potentials are employed, a Coulomb potential, a short-range Coulomb potential, and a soft-box Voronoi potential. Effects of different orders of the extrapolating Padé approximant are investigated, and last, the energy range of input data for the extrapolation is studied. Moreover, all ACCC results are compared to resonance parameters that have been independently obtained with the same theoretical method, but with a different continuum approachcomplex scaling for the model and complex absorbing potentials for the temporary anions.
1. INTRODUCTION An electronic resonance is an electronic state of an (N + 1) -electron system that lies energetically above the N-electron system, and therefore in the scattering continuum of the N-electron state. Examples are temporary anions such as N2− or the benzene anion, small dianions in the gas phase such as CO32− or SO42−, core-ionized atoms and molecules decaying via the Auger process, and the analogous Auger-like decay of innervalence ionized clusters. These states are characterized by their so-called resonance position, Er, that is, the energy above the N-electron state, and by their width, Γ, which playscum grano salisthe role of a first-order decay constant. In other words, there is a typical decay time, τ = ℏ/Γ, which is inversely proportional to the width.1,2 Characterizing a resonance with standardor at least more or less standardquantum chemistry methods is a challenge. If a standard basis set is used, the continuum is so poorly described that any description of the decay is suppressed, and the resonance is effectively trapped in a basis-set box. Using more and more diffuse functions gives a better and better discretization of the continuum; however, then the resonance mixes with the continuum states, and from a straightforward calculation one obtains at best a rough idea of the resonance position but at worst the lowest energy of the discretized continuum. © XXXX American Chemical Society
Almost from the beginning of computational chemistry, various so-called L2-methods have been put forward to address this challenge. All L2-methods have one thing in common: The original Hamiltonian is parametrized, and Er and Γ are extracted from a study of the behavior of its eigenvalues as functions of the parameter. Examples for these methods are complex scaling and complex basis functions,3−5 complex absorbing potentials,6,7 the stabilization method,8,9 and extrapolation methods1,10 (see ref 11 for recently identified limitations of the stabilization method). In most of these methods, several electronic states are needed for the eigenvalue analysisif not explicitly, then at least to identify the resonance stateand therefore these continuum methods areas a rulecombined with an electronic-structure method that yields several electronic states. Examples are configuration interaction,12 symmetryadapted cluster-configuration interaction and equation-of-motion coupled-cluster methods,13−18 or Green’s functions methods.19−21 That is, when doing calculations for resonances a continuummethod/electronic-structure-method/basis-set combination is used, and it is this “three-dimensional space” that makes it effectively impossible to directly compare different combinations used by different research groups in a meaningful way. Received: December 19, 2016 Published: April 20, 2017 A
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
computed for increasing strengths of an artificial stabilizing potential, and both the ab initio calculations of the electron binding energy and the specific stabilizing potentials are described. In the second step, the obtained data are backextrapolated to vanishing stabilization via analytic continuation, and our specific data analysis algorithms are outlined below. 2.1. Ab Initio Methods. First the geometries of the neutral molecules C2H4 and CH2O were optimized using the PW6B95 density functional30 and Dunning’s correlation-consistent triple-ζ (cc-pVTZ) basis set31 as implemented in the ORCA package.32 After that, both continuum methods used for the molecular resonances, ACCC and the complex absorbing potential (CAP) method, involve the computation of electron attachment energies. Yet, the specific requirements of both methods are substantially different. For example, in ACCC, in principle, only the lowest attachment energy is needed, while a CAP calculation typically needs a substantial number. In contrast, the electron structure calculation has to be repeated many times for ACCC, while a single solution is sufficient for the projected CAP formalism. Since our goal is a direct comparison of ACCC and CAP resonance parameters, the symmetry-adapted cluster−configuration interaction (SAC−CI) method,33,34 an electronic-structure method that works well for both continuum approaches, ACCC and CAP, was employed. The SAC−CI method is a high-level method for computing electron attachment energies. It is essentially equivalent to equation-of-motion coupled-cluster methods for electron attachment,35−38 that is, it is size-consistent, by construction balanced, and with respect to the closed-shell neutral, the anion states are described by taking into consideration all one-particle and two-particle−one-hole configurations explicitly and all three-particle−two-hole configurations implicitly. Moreover, SAC−CI is a direct method; that is, after the correlated electronic ground state of the neutral molecule has been computed, several attachment states ((N + 1)-electron states) are computed directly without any reference to self-consistentfield wave functions for open-shell anions, and the method naturally yields more than one state. Here the investigated molecules are small, so we do not use any of the further approximations available in SAC−CI (variational SAC−CI, configuration selection through perturbation criteria); however, use is made of the direct algorithm where σ-vectors are directly calculated including all the product (nonlinear) terms.39 The SAC−CI calculations were done with the GAUSSIAN09 suite of programs, Revision B.01.40 The CAP SAC−CI combination (CAP/SAC−CI) has been extensively studied in the literature.16,41−44 Specifically, the potential used as a CAP is a cutoff harmonic potential in the distance to the nearest atom, rnearest (cf. eq 5, that is, a smoothed Voronoi-shaped CAP, which is zero in the vicinity of the molecule and then grows quadratically with the distance to the nearest atom.41,42 For the CAP calculations the cc-pVTZ valence basis set was augmented with a (2s5p2d) set of diffuse functions on C and O (even-scaled exponents starting from the smallest exponents for the respective angular momentum in the original basis set; even-scaling factor of 2.0 for s and d functions, and 1.5 for p functions) and a (2s2p) set on H (even-scaling factor of 2). For an ACCC calculation the SAC−CI calculation has to be repeated for many strengths of the artificial stabilizing potential, that is, for many values of λ (see subsection 2.2). In addition to the advantage of needing fewer electron attachment energies,
Here, the main focus is extrapolation methods in the form of the analytic continuation of the coupling constant (ACCC).10,22,23 As in all L2-methods, the ACCC method consists of two steps: In the first step the parameter in a parametrized Hamiltonian is changed, and data in the form of (parameter, energy) pairs are collected. In the ACCC method an artificial stabilizing potential, V, is added to the Hamiltonian, H(λ) = H + λV, and the parameter λ is increased until the resonance is turned into a bound state. Then λ is further increased until enough data are collected so that the boundstate energies can be back-extrapolated to λ = 0. Since only bound-state data are used for the back-extrapolation, standard quantum chemistry ground-state methods can be employed and all wave functions are obviously square integrable (L2). While it is not strictly necessary to compute more than one state in the ACCC scheme, it is often helpful to compute two, because the position and shape of the avoided crossing of these two states in an energy vs λ plot can be vital for the data analysis algorithms. The second step is then the back-extrapolation to λ = 0. A straightforward extrapolation of the real energies E(λ) just yields approximations of Er, but no lifetime. Lifetimes can only be obtained if the extrapolating function is allowed to become complex when the resonance state crosses the threshold from the continuum to the bound states, and this is achieved by the ACCC method.10,22,23 Currently, there are three variants of ACCC,24−27 the original method, an “inverse” method, and a “regularized” method. All three use variable transformations and perform the back-extrapolation by fitting the transformed data to Padé approximants, and all three are in detail described in subsection 2.3. While each variant in itself has been thoroughly tested, what is so far lacking are studies comparing the three different ACCC schemes with each other andas directly as possibleto alternative methods. One study compares normal ACCC to CAPbased methods for exactly the same electronic-structure method and basis set, and in addition considers one additional stabilizing potential.28 Moreover, after this manuscript had been submitted, a study that compares the inverse and regularized ACCC methods and considers three different stabilizing potentials appeared.29 Since one strong emphasis of ref 29 is basis-set and method convergence, but two aspects of this study are not considered in ref 29, ref 29 and this work complement each other well. The main goal of this study is to directly compare the three established variants of the ACCC method with each other and with independently obtained resonance parameters that are computed using the same electronic structure method but a different approach of addressing the continuum. For this purpose, we introduce first two preconditioned algorithms for the data analysis steps of the ACCC and inverse ACCC variants that lead to more robust computational protocols. In addition, two other aspects of the ACCC method will be examined, the sensitivity of the ACCC results with respect to the range of data used for the back-extrapolation and the type of artificial stabilizing potential V added to the physical Hamiltonian. The guinea pigs for the investigations are three resonances: that of a model system as well as two well-known temporary anion resonances, the 2B2g resonance of the ethylene and the 2B1 resonance of the formaldehyde anions.
2. COMPUTATIONAL METHODS This section discusses details of both steps of the ACCC method. In the first step, the electron binding energy must be B
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
where the R⃗ α are the positions of the nuclei. This definition of r is combined with a one-dimensional soft-box potential
another benefit of ACCC is that it requires a much smaller augmentation with diffuse basis functions than CAP calculations, even though the standard augmentation with just a single diffuse function is typically insufficient.28 Therefore, in all ACCC calculations the cc-pVTZ set was augmented with a (2s2p1d) set on C and O (even-scaling factor 3) and a (1s1p) set on H (even-scaling factors of 3 for all angular momenta). 2.2. Stabilizing Potentials. The first step of the ACCC method is the data collection step. An artificial stabilizing potential V is added to the molecular Hamiltonian H H (λ ) = H + λ V
VV (r ) = s(r ) − 1
where s(r) is a soft-step function ⎧ 0: r≤a ⎪ 3 s(r ) = ⎨ ρ (10 − ρ(15 − 6ρ)): a < r < b ⎪ ⎩1: r≥b
(7)
with
(1)
where λ is the strength parameter or in this context the so-called coupling constant. Switching on the potential V must transform the resonance into a bound state; that is, it must bind the excess electron to the molecule. The only condition V must in principle fulfill is that it should be short-ranged in the usual scattering theory sense.2 Nevertheless, owing to ease of implementation, most ACCC applications for temporary anions have successfully used a long-range Coulomb potential; noteworthy exceptions are refs 28 and 29. Here we compare three different stabilizing potentials: a Coulomb potential, a short-range Coulomb potential, and a Voronoi-shaped soft-box potential. The Coulomb and Voronoi-box potentials are identical to those used in ref 28. The particular Coulomb potential used here can be thought of as replacing each physical nuclear charge, Zα, by a scaled nuclear charge45 Zα(λ) = (1 + λ)Zα
(6)
ρ=
r−a b−a
(8)
being a scaled distance, so that s(r) increases smoothly from 0 to 1 within the r range [a, b]. This one-dimensional potential is shown in Figure 1.
(2)
In terms of a stabilizing potential VC this can be written as VC =
∑− α
Zα |r − R α |
(3) Figure 1. “Radial” part of the soft-box Voronoi potential VV defined in eq 6.
where Rα are the nuclear positions. The short-range Coulomb potential retains the short-range part of VC but cuts off the long-range part by multiplication with a Gaussian function ⎛ |r − R | ⎞ Zα α ⎟ exp⎜ − 2 |r − R α | r ⎝ ⎠ 0
In contrast to VC and VSRC, the Gaussian integrals of VV must be evaluated numerically, and for this purpose an angular Lebedev grid46 with 3470 points, and a radial grid of 500 equidistant points (Euler−Maclaurin quadrature with an upper limit for the radial integration of 5.5 times the half-width of the product Gaussian). is placed at the center of the respective product Gaussian (cf. ref 41). All three potentials are compared for the formaldehyde molecule in Figure 2. 2.3. Preconditioned ACCC and Inverse ACCC Algorithms. The second step of the ACCC method is backextrapolation to λ = 0. Yetas mentioned abovestraightforward extrapolations of the bound energies will only yield the resonance positionto compute the complex resonance energy it is necessary to extrapolate by “analytic continuation” so that the energy can become complex once it crosses the threshold. In the ACCC method this goal is achieved by a variable transformation: Instead of E(λ), k(y) is extrapolated, where k = −2E is a momentum-like variable and y = λ − λ 0 , where λ0 is the position where the relevant state crosses the threshold (E(λ0) = 0). Note that y is real for λ > λ0 but imaginary for λ < λ0 and that in contrast to ref 24 we use the normal orbitalenergy-like quantum chemistry convention: positive energies mean that an electron is unbound; negative energies mean that it is bound. The extrapolation itself is achieved by fitting the kj(yj) data points, which are all at real y, by a Padé approximant, and
2
VSRC =
∑− α
(4)
Here, r0 is a parameter that determines the cutoff length. In the text we will focus on results obtained with r0 = 5 bohrs; for significantly larger values of r0, results obtained with VSRC and VC become effectively identical due to the finite basis set; for considerably shorter r0 (2 or 1 bohr) the strength parameter λ of VSRC has to be increased so drastically to create a bound state that the extrapolation procedure becomes far more unstable. However, there is clearly optimization potential in this parameter. The third stabilizing potential is a soft, finite-depth Voronoishaped box. In other words, this potential consists of a “radial” soft-box combined with an “angular” Voronoi surface around the molecule, that is, a surface that wraps around a molecule similar to a van der Waals surface. We will refer to this potential as a soft-box Voronoi or just a Voronoi potential. Of course, a truly radial coordinate is well-defined in atoms only, and the Voronoi surface is obtained by substituting it with the distance to the nearest nucleus r = rnearest( r ⃗) = min(| r ⃗ − R⃗ α|) α
(5) C
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Figure 2. Cuts through VC (eq 3), VSRC (eq 4), and VV (eq 6) for formaldehyde. For VC and VSRC the color scale runs from −27 hartrees (white) to 0 (dark red) in steps of 3 hartrees; for VV the color scale runs from −1 hartree (white) to 0 (red) in steps of 0.2 hartree.
First, a good guess for λ0 is obtained by fitting data for the lowest state of the Hamiltonian, E(λ), to a suitable function and locating the point where the function moves through threshold. The details of this procedure are unimportant; a low-order polynomial interpolation would probably do the trick. However, we tend to use a generalized Padé approximant47 and solve the resulting system of linear equations by singular-value decomposition (SVD) so that more input-data points than output Padé parameters can be used (see discussion in the next paragraph for more detail). Moreover, we use energies from both sides of the threshold, which may at first sight seem not entirely sound. Note, however, that in practical calculations SAC−CI electron attachment energies are smooth functions of λ (e.g., Figure 3) that cross the threshold without any ado. It is
evaluating the fit at λ = 0, which yields the complex resonance 1 energy − 2 k( −λ 0)2 = Eres = Er − i Γ/2.10,22,23 In the following we will briefly describe the ACCC algorithm as given in ref 24 as well as our preconditioned algorithm. Both start out from a data set (λj, Ej) of L points where all Ej < 0. The differences lie in the procedure for how the optimal Padé approximant is found. In the original ACCC algorithm24 the basic assumption is that λ0 is only approximately known and must be determined by the quality of the fit of the Padé approximant itself. In an outer loop λ0 is changed in some way, and in an inner loop the quantity χ2 =
1 L
∑ j
1 PN (y) − kj σj 2 Q N (y)
(9)
is minimized, where P(y)/Q(y) is the Nth-order Padé approximant currently fit to the data N
P(y) =
∑ pyi i i=0
N
and
Q (y ) = 1 +
∑ qiyi i=1
(10)
and the 1/σj2 are fitting weights, which can be interpreted as relative errors σj of each data point. The inner minimization step is a nonlinear optimization for which the authors of ref 24 recommend the Powell or Nelder−Mead simplex method, and the Padé approximant obtained by this means is referred to as a Padé-III approximant.10 Regardless of whether the outer loop simply steps through a certain λ0 range, or a more sophisticated minimization procedure is used, the λ0 value that minimizes χ2 is to be found and used. The 2010 algorithm works well for Padé order N = 1; however, for N > 1 the χ2 functions turn out to be extremely rugged and show a large amount of local minima. Even for a fixed λ0 there are typically several minima of similar depth. As all nonlinear optimizations, the 2010 algorithm can be considerably stabilized by using a good initial guess: Select 2N + 1 data points from the original data set, find the Padé approximant that exactly fits this subset by solving the standard Padé system of equations (Padé-II approximant), and use the parameters of the Padé-II approximant as a guess for the optimization of the Padé-III approximant. Nonetheless, if one plots χ2 as a function of λ0as a rulea jagged curve with strong oscillations and finite jumps is obtained for N > 1. Here we suggest an alternative algorithm that optimizes the Padé parameters and λ0 simultaneously. At first glance that may look like a step in the wrong direction, because it adds further nonlinear character to an inherently unstable algorithm making it even more unstable. Howeveras a rulethe resulting algorithm becomes very robust if it is provided with a good guess that is obtained as follows.
Figure 3. Three lowest states of the formaldehyde anion stabilized with a soft Voronoi box (cf. eq 6). For λ = 0 all states are unbound, but the continuum is discrete due to the Gaussian basis set. If λ is increased the “resonance” state is stabilized more strongly than the discretized continuum, consequently undergoes avoided crossings with the discretized-continuum states, and finally becomes a bound state.
clearly imperative to choose a stabilizing potential so that this behavior of the energies is ensured. This procedure is only needed once, say for ACCC with the low-order Padé approximant [2,2]. After that λ0 remains practically unchanged, and can be used as a guess for all susequent ACCC calculations with higher Padé orders. The next agenda item is the selection of the input data to be used in the ACCC algorithm. In this communication we will select all bound-state energies between an upper and a lower limit, and this range will be systematically expanded by decreasing the lower limit. The upper range should be chosen so that it is sufficiently far from the lowest avoided crossing between the lowest two adiabatic states (Figure 3) and may D
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
(see, e.g., ref 48) with numeric gradients for solving the nonlinear optimization problem of eq 9 simultaneously with finding the optimal λ0. The BFGS algorithm is a quasi Newton method that should find the local minimum downhill from the current position; it has been shown to have good performance even for nonsmooth optimizations, and findsin our experiencethe optimal local χ2 in the overwhelming number of cases quickly and efficiently. In other words, while there are exceptionsthis is still a highly nonlinear fitting procedure after allthe resulting algorithm provides a much improved ACCC algorithm in terms of robustness and in terms of drastically reduced user interaction. The same strategy can be applied to the inverse ACCC (IACCC) method, where instead of k(y) the inverse function λ(k) is fit with a Padé approximant. The position of the threshold, λ0, is not required in the IACCC method, so the SVD can be used directly to obtain a guess for the Padé parameters. The nonlinear fitting procedure for IACCC is almost analogous to eq 9
Figure 4. Radial part of the model potential defined in eq 14.
well be at substantially negative energies.28 However, with the basis set used here, the pseudocontinuum is so purely described that all bound states up to threshold can be included without approaching the avoided crossing too closely. Second, a good guess for the Padé-III approximant is obtained by a variant of the computation of finding Padé-II parameters. Instead of using just a subset of 2N + 1 data points that a Padé-II approximant will match exactly, one can simply use all data points. This leads to an overcomplete set of linear equations for the Padé parameters kj = PN (y) − kjQ N (y)
χ2 =
1 L
∑ j
1 PN′ (k) − λj σj 2 Q N′ (k)
(12)
the only difference is that the correct threshold behavior of λ(k) is typically built into the Padé approximant; that is, P′N(k) and Q′N(k) do not have linear terms PN′ (k) = p0 + p2 k 2 + p3 k3 + ...
and
Q N′ (k) = 1 + q2k 2 + q3k3 + ...
(13)
We find that for IACCC the linear fit obtained from SVD even though it may formally not be equivalent to a true Padé-III approximantis very close to a local minimum of the nonlinear procedure from eq 12 in the sense that the two χ2 values often agree with each other to two significant figures. We emphasize again that solving the SVD is effectively for free, while the nonlinear minimization procedure becomes time intensive for higher Padé orders. Unfortunately, SVD cannot be applied to the so-called regularized analytic continuation (RAC) method,26,27 because
(11)
however, the SVD procedure handles that easily and returns the best fit to this linear problem, which is, owing to its linearity, easy to compute. It is obviously an improvement on a Padé-II approximant, butunfortunatelynot quite the solution of eq 9. Nevertheless, it makes an excellent guess for the nonlinear optimization algorithm. Third, after these two preconditioning steps, we use the Broyden−Fletcher−Goldfarb−Shanno (BFGS) algorithm
Figure 5. Resonance parameters, Er and Γ, of the model potential described in the text (eq 14) in dependence on the range of input data selected for the back-extrapolation. The units are “reduced units” in reference to the CS resonance energy, ECS; that is, the input range on the x-axis is [0:E/ECS r ] CS and the computed resonance position and width on the y-axis are Er/ECS r and Γ/Γ , respectively. For ACCC and IACCC, blue, green, and yellow code Padé [4,4], [5,5], and [6,6] approximants; for RAC, blue and green code [3,1] and [4,2] approximants. E
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
of 10−6 or 10−5 eV. Note that this error is not random but systematic and due to convergence thresholds in the integral, self-consistent-field, and correlation calculations. Moreover, the model will help to introduce the reduced units we use for comparison of the different ACCC variants. The model itself is a simple radial potential
Table 1. CAP Results for the 2B2g and 2B1 States of the Ethylene and Formaldehyde Temporary Anions C2H4− CH2O−
uncorrected corrected uncorrected corrected
Er (eV)
Γ (eV)
2.111 2.023 1.241 1.175
0.506 0.458 0.387 0.382
2
VM(r ) = (ar 2 − b)e−cr +
the RAC Padé approximant is written in a product form that cannot be brought into a linear form regarding the fitting parameters. Instead, one must deal with a nonlinear fit from the start, which we found to be in itself only stable for the least flexible [2,1] approximant. For all higher approximants we tried[3,1], [3,2], and [4,2]we had to fall back to a twostep fitting procedure, where first α and β describing the resonance pole (cf. ref 27) are held fixed at their optimal values from the [2,1] approximant while the new parameters are fitted, and only then are all parameters allowed to change freely.
l(l + 1) 2r 2
(14)
consisting of a damped harmonic term that creates a potential well and barrier structure at the origin, and an angular momentum term (Figure 4). We use atomic units, a particle with mass of m = 1, angular momentum l = 1, and the three parameters of the potential are a = 0.1, b = 1.2, and c = 0.1. All calculations for this model are performed using a discrete variable representation based on 600 sine functions leading to an equidistant-grid representation in the r range [0,40]. Employing this grid representation and the complex scaling method yields a CS resonance energy of ECS = r = 1.7587 eV and a width of Γ 0.397 eV, a value that can be considered as a practically exact reference value for comparison with different flavors of ACCC provided the same grid is used. To obtain input data for the ACCC, IACCC, and RAC backextrapolations, a soft-box potential in the spirit of the Voronoi stabilization for molecules was used (eq 6 with a = 2 and b = 10). Then the three extrapolation algorithms were applied several times with larger and larger sets of input data. The smallest set used for the extrapolation consisted of all bound data in the
3. RESULTS 3.1. Model. We start with results for a simple model system. The primary purpose of the model is to investigate the performance of the different Padé extrapolations given precise input data. The model itself is solved with a discrete variable representation,49 and the computed energies have full numeric precision. In contrast, electron binding energies computed with ab initio packages typically have errors in the range
Figure 6. C2H4− resonance position, Er, in dependence on the range of input data selected for the back-extrapolation. The units are “reduced units” CAP in reference to the CAP resonance position, ECAP r , that is, the input range on the x-axis is 0 to E/Er , and the computed resonance position on the . For ACCC and IACCC, blue, green, and yellow code Padé [4,4], [5,5], and [6,6] approximants; for RAC, blue and green code y-axis is Er/ECAP r [3,1] and [4,2] approximants. F
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Figure 7. C2H4− resonance width, Γ, in dependence on the range of input data selected for the back-extrapolation. The units are “reduced units” in CAP . reference to the CAP resonance energy, ECAP; that is, the input range on the x-axis is 0 to E/ECAP r , and the computed width on the y-axis is Γ/Γ For ACCC and IACCC, blue, green, and yellow code Padé [4,4], [5,5], and [6,6] approximants; for RAC, blue and green code [3,1] and [4,2] approximants.
range from [0:−3ECS r ] (50 energies), and this set was extended point by point to the range [0:−4.2ECS r ] (65 energies). The results for both Er and Γ are shown in Figure 5 in reduced units; that is, instead of showing the absolute energies in electronvolts, the energies relative to the complex-scaling result, CS Er/ECS r and Γ/Γ , are plotted vs the input-data range. The resonance positions computed with preconditioned ACCC and IACCC are both independent of the input-data range and show excellent agreement with the CS result; that is, in Figure 5 the reduced position is very close to 1.0 regardless of the Padé order and the input-data range. We emphasize that despite the need of fitting λ0 the preconditioned ACCC algorithm is very robust. It quickly finds an optimal λ0−Padé parameter combination, and apart from Padé [1,1], all optimized λ0 values are practically identical; that is, the local minima in the (λ0, Padé parameter) space form one family, and the ACCC algorithm follows that family to higher and higher Padé orders. The RAC algorithm, on the other hand, gives somewhat less convincing agreement with CS; however, it still shows good consistency with reduced Er values in the range between 0.95 and 1.03. For all three algorithms, the agreement of the resonance width with the CS result is worse than the agreement of the resonant positions. For preconditioned ACCC and IACCC this means that the agreement is not excellent but merely good with reduced resonance widths in the range of 0.96−1.04. However, for RAC the computed widths are too small by roughly a factor of 4. (Note that in Figure 5, in various panels many data points of the lower order approximants cannot be seen, because on the scale of the figures the data points of the higher order
approximants are plotted right on top of them. The same will be true for all similar subsequent figures.) 3.2. π* Resonances of C2H4 and CH2O. In this section the π* temporary anion resonances of ethylene (H2CCH2) and formaldehyde (H2CO) are studied. Our focus is the performance of the three different ACCC variants with different stabilizing potentials. For this purpose it is important to use common sets of input data if extrapolation variants are compared for the same stabilizing potential, as well as a similar number and range of input data if results for different stabilizing potentials are compared with each other. Moreover, it is critical to have an independent basis of comparison that can be related as directly as possible to the ACCC results, in other words, a theoretical result that has been obtained at identical geometries, with an identical electronic-structure method and with an identical valence basis set. Here the CAP/SAC−CI method as described in section 2.1 is used for this purpose. Our current CAP results are collected in Table 1; for a recent comprehensive compilation of previous experimental and theoretical results for these two resonances, see ref 18. The CAP results in Table 1 agree well with other CAP values;16,18 however, even for closely related calculations, such as SAC−CI and EOM-CCSD with cc-pVTZ valence basis sets, there are, of course, small deviations, since different geometries, different CAPs, and different augmentations of the cc-pVTZ basis set have been used. It should also be mentioned that the 2 B2g resonance of the ethylene anion has been thoroughly investigated using the IACCC method,25 and the result closest to the geometry considered here is Er = 1.918 eV and Γ = 0.511 eV. G
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Figure 8. CH2O− resonance position, Er, in dependence on the range of input data selected for the back-extrapolation. The units are “reduced units” CAP in reference to the CAP resonance position, ECAP r ; that is, the input range on the x-axis is 0 to E/Er , and the computed resonance position on the CAP y-axis is Er/Er . For ACCC and IACCC, blue, green, and yellow code Padé [4,4], [5,5], and [6,6] approximants; for RAC, blue and green code [3,1] and [4,2] approximants.
These values are in good agreement with our findings; however, they cannot be directly compared to our CAP/SAC−CI results, since different CH bond lengths and bond angles, different valence basis sets, and different augmentations with diffuse functions, as well as different electron correlation methods, have been used.25 Turning to SAC−CI based ACCC calculations, the results for the resonance position and width of the 2B2g state of H2CCH−2 are collected in Figures 6 and 7 and the corresponding results for the 2B1 state of H2CO− can be found in Figures 8 and 9. Each of these figures collects results for three different extrapolation methods and three different stabilization potentials, and in the figures these nine combinations are abbreviated, for example, as C-ACCC for Coulomb-ACCC, SRC-RAC for short-range Coulomb RAC, or V-IACCC for Voronoi IACCC. For each of the nine combinations the input-data range is CAP varied starting from [0,−3ECAP r ] to about [0,−5.5Er ], and the results are shown in reduced units, where the respective corrected CAP resonance energy serves as a reference value. In contrast to the model discussed in the last subsection, the SAC−CI electron binding energies serving as input data for the different extrapolation methods have an absolute precision of only about 10−6 eV effectively due to the convergence of the SAC−CI iteration itself. Nonetheless, the agreement for the resonance positions (Figures 6 and 8) is still good for most combinations. For ethylene, in particular the Coulomb IACCC combination shows excellent agreement with the CAP position; for formaldehyde, it is the short-range Coulomb ACCC and IACCC combinations, the Voronoi ACCC and
RAC combinations, and Coulomb IACCC for Padé orders beyond [4,4]. RAC Coulomb and short-range Coulomb combinations show the worst agreement (see scales in Figures 6 and 8) as well as large jumps in their input range dependence seemingly indicative of the inherently nonlinear fitting procedure underlying the RAC method. Yet, for some reason these instabilities can be avoided if a Voronoi potential is used as a stabilizing potential, which renders the RAC extrapolation very stable and astonishingly independent of the Padé order. Unfortunately, the agreement of the computed widthswith the CAP result and with each otheris less than satisfactory (see Figures 7 and 9). For ethylene, only Coulomb and Voronoi stabilization potentials combined with ACCC and IACCC get good agreement with the CAP width, but not for all Padé orders and not over the whole input-data range. An interesting case is Coulomb RAC, which agrees well with the CAP result, but only for Padé [3,1] and only if the input-data range is larger than [0,−4ECAP r ]. Other combinations yield widths in the right order of magnitude, with all short-range Coulomb widths being significantly too large and the Voronoi RAC widths being significantly too small. Even more unfortunate, very few of these trends carry over to formaldehyde (Figure 9). For example, for formaldehyde, the short-range Coulomb ACCC extrapolation gives remarkable good agreement with the CAP width, whereas the Coulomb ACCC and IACCC methods yield far too small values. Only the good agreement of Coulomb RAC Padé [3,1] and the too low widths in Voronoi RAC seem consistent with the ethylene results. In other words, as typical for virtually all methods for resonances, “widths are hard” in the H
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Figure 9. CH2O− resonance width, Γ, in dependence on the range of input data selected for the back-extrapolation. The units are “reduced units” in CAP . reference to the CAP resonance energy, ECAP; that is, the input range on the x-axis is 0 to E/ECAP r , and the computed width on the y-axis is Γ/Γ For ACCC and IACCC, blue, green, and yellow code Padé [4,4], [5,5], and [6,6] approximants; for RAC, blue and green code [3,1] and [4,2] approximants.
resonance/ACCC-variant/stabilizing potential, the effect of the range of bound-state input data is investigated. All ranges start at zero, and the lower limit was extended from about −3 times the respective resonance position to between 4−6 times this value. Last, but not least, the ACCC findings are compared to resonance parameters computed with a different continuum method without changing any other aspects of the calculation. For the model, complex scaling is employed for this purpose; for ethylene and formaldehyde the CAP method is used. The model allows us to study the ACCC method under ideal conditions; that is, first, the input data for the extrapolation have full numerical precision, and second, the “exact” solution of the model problem is known from the complex scaling calculations. Under these conditions all variants show excellent to good agreement for the resonance position; ACCC and IACCC yield good widths, and only the RAC width is significantly too small. The ab initio data for the molecular resonances show far more mixed results. One clear trend is that resonance positions can probably be expected to be reliable, while larger deviations should be expected for the widths. To get a fair first estimate of the fluctuations out of a single data set of bound-state data, one should, on the one hand, investigate the dependency on the input-data range, the Padé order (cf. ref 25), and the extrapolation method itself, a task that may seem tedious but is computationally inexpensive. In a more thorough, but also considerably more costly, study one would also consider different stabilizing potentials, because there seems to be no clear trend whichif anyis superior to the others, and at least for the
sense that they are far more sensitive to the particular details of the calculation: The results depend strongly on the artificial stabilizing potential used to obtain the input data, the different varieties of the ACCC method yield widths that are only consistent to within an order of magnitude or so, and the results depend far more strongly on the input-data range than the resonance positions.
4. CONCLUSIONS In this study the performance of the three variants of the ACCC method, the original ACCC (fitting k( λ − λ 0 )24), IACCC (fitting λ(k)25), and RAC (fitting λ(k) with tailor-made Padé approximants26,27), is examined. Based upon an initial linear fit with the help of singular-value decomposition, preconditioned algorithms for the nonlinear fitting steps of ACCC and IACCC are introduced, and these methods allow us to easily and robustly converge high-order Padé approximants. Using the preconditioned algorithms for ACCC and IACCC as well as the established RAC approach, the three ACCC variants are compared for a model potential and for the π* temporary anions of ethylene and formaldehyde. For the model, bound data are obtained with only a single artificial stabilizing potential, whereas for the molecular resonances three different stabilizing potentials are compared: a Coulomb potential equivalent to scaling all nuclear charges of the respective molecule, a short-range version of this Coulomb potential, which damps the long-range Coulomb tails with a Gaussian function, and a soft-box Voronoi potential that wraps around the molecule like a van der Waals surface. For all combinations, I
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Journal of Chemical Theory and Computation
■
three stabilizing potentials considered here, the answer is system dependent. In this regard there is a clear lack of data in the literature; more systematic studies where everything apart from the stabilizing potential is kept fixed are needed. Just to keep everything in perspective, similar critical study would be highly beneficial for other continuum methods, too. For example, one could image to examine the stabilization method: One would compare stabilization graphs created with its exponent scaling incarnation8,9 and with explicit hard-50 and soft-box potentials. The analysis would then be done with three different approaches: Taylor−Mandelshtam phase-shift fitting,51 analytic continuation of different avoid crossings,52 and the Löwdin method.53 Since the latter two analysis methods yield (at least) one result per avoided crossing, one could create a large number of results, and there would clearly be a certain scatter for the energies and widths; however, this scatter is unfortunately unknown, because such a systematic study has not been done. For CAPs the situation is slightly better, since there are some systematic studies that vary the cutoff radius and others that also compare different mathematical forms of CAPs (see, e.g., refs 18 and 41), but, for example, the question of whether the complex energy should be stabilized or resonance position and width should be stabilized separately has not been settled (cf. refs 6 and 54), and the same is true for the question of whether the so-called first-order correction is really improving things (see in particular the disscussion in ref 6). It is probably fair to say that the scatter between different CAP variants can be as large as 0.2 eV (mostly due to to the CAP shape and the correction) but it can be potentially much larger in the width for small cutoff radii.41 The ACCC method clearly has a lot of potential, in particular when considering resonances of intermediate-size systems. For the small molecules considered here the CAP method is by almost 2 orders of magnitude more efficient than the ACCC approach, because owing to the projected technique only a single SAC−CI calculation is required. However, ACCC has two clear advantages: First, the valence basis set needs to be augmented with a much smaller set of diffuse functions, alleviating the near-linear dependency problems the CAP method faces for intermediate-size molecules, and second, the number of SAC−CI eigenvalues needed for the subsequent data analysis are much smaller, so that ACCC must become more efficient or even the only viable approach for larger molecules.
■
Article
REFERENCES
(1) Jordan, K. D.; Voora, V. K.; Simons, J. Negative Electron Affinities from Conventional Electronic Structure Methods. Theor. Chem. Acc. 2014, 133, 1445. (2) Kukulin, V. I.; Krasnopolsky, V. M.; Horácě k, J. Theory of Resonances; Kluwer Academic, Dordrecht, The Netherlands, 1989; pp 88−154. (3) McCurdy, C.; Rescigno, R. Extension of the Method of Complex Basis Functions to Molecular Resonances. Phys. Rev. Lett. 1978, 41, 1364. (4) Reinhardt, W. P. Complex coordinates in the theory of atomic and molecular structure and dynamics. Annu. Rev. Phys. Chem. 1982, 33, 223−255. (5) Samanta, K.; Yeager, D. L. Complex multiconfigurational selfconsistent field based methods to investigate electron-atom/molecule stattering resonances. Advances in Chemical Physics 2012, 150, 103− 142. (6) Riss, U. V.; Meyer, H.-D. Calculation of resonance energies and widths using the complex absorbing potential method. J. Phys. B: At., Mol. Opt. Phys. 1993, 26, 4503. (7) Santra, R.; Cederbaum, L. S. Non-Hermitian electronic theory and application to clusters. Phys. Rep. 2002, 368, 1−117. (8) Hazi, A. U.; Taylor, H. S. Stabilization Method of Calculating Resonance Energies: Model Problem. Phys. Rev. A: At., Mol., Opt. Phys. 1970, 1, 1109. (9) Hazi, A. U. A purely L2 method for calculating resonance widths. J. Phys. B: At. Mol. Phys. 1978, 11, L259. (10) Kukulin, V. I.; Krasnopolsky, V. M.; Horácě k, J. Theory of Resonances; Kluwer Academic: Dordrecht, The Netherlands, 1989; pp 219−270. (11) Falcetta, M. F.; Fair, M. C.; Tharnish, E. M.; Williams, L. M.; Hayes, N. J.; Jordan, K. D. Ab initio calculation of the cross sections for electron impact vibrational excitation of CO via the 2II shape resonance. J. Chem. Phys. 2016, 144, 104303. (12) Sommerfeld, T.; Riss, U. V.; Meyer, H.-D.; Cederbaum, L. S.; Engels, B.; Suter, H. U. Temporary Anions - Calculation of Energy and Lifetime by Absorbing Potentials: the N−2 2Πg resonance. J. Phys. B: At., Mol. Opt. Phys. 1998, 31, 4107. (13) Sajeev, Y.; Mishra, M. K.; Vaval, N.; Pal, S. Fock space multireference coupled cluster calculations based on an underlying bivariational self-consistent field on Auger and shape resonances. J. Chem. Phys. 2004, 120, 67. (14) Sajeev, Y.; Santra, R.; Pal, S. Analytically continued Fock space multireference coupled-cluster theory: Application to the 2Πg shape resonance in e-N2 scattering. J. Chem. Phys. 2005, 122, 234320. (15) Sajeev, Y.; Pal, S. A general formalism of the Fock space multireference coupled cluster method for investigating molecular electronic resonances. Mol. Phys. 2005, 103, 2267−2275. (16) Ehara, M.; Sommerfeld, T. CAP/SAC-CI method for calculating resonance states of metastable anions. Chem. Phys. Lett. 2012, 537, 107. (17) Ghosh, A.; Vaval, N.; Pal, S. Equation-of-motion coupled-cluster method for the study of shape resonance. J. Chem. Phys. 2012, 136, 234110. (18) Zuev, D.; Jagau, T.-C.; Bravaya, K. B.; Epifanovsky, E.; Shao, Y.; Sundstrom, E.; Head-Gordon, M.; Krylov, A. I. Complex absorbing potentials within the EOM-CC family of methods: Theory implementaion and benchmarks. J. Chem. Phys. 2014, 141, 024102. (19) Santra, R.; Cederbaum, L. S. Complex absorbing potentials in the framework of electron propagator theory I. General formalism. J. Chem. Phys. 2002, 117, 5511. (20) Feuerbacher, S.; Sommerfeld, T.; Santra, R.; Cederbaum, L. S. Complex absorbing potentials in the framework of electron propagator theory. II. Application to temporary anions. J. Chem. Phys. 2003, 118, 6188. (21) Tsednee, T.; Liang, L.; Yeager, D. L. Electron-atom resonances: The complex-scaled multiconfigurational spin-tensor electron propagator method for the 2P Be− shape resonance problem. Phys. Rev. A: At., Mol., Opt. Phys. 2015, 91, 022506.
AUTHOR INFORMATION
Corresponding Authors
*(T.S.) E-mail:
[email protected]. *(M.E.) E-mail:
[email protected]. ORCID
Thomas Sommerfeld: 0000-0001-8105-5414 Masahiro Ehara: 0000-0002-2185-0077 Funding
T.S. acknowledges support by the National Science Foundation under Grant CHE-1565495. M.E. acknowledges the financial support from a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (JSPS Grants 16H04104 and 16H06511). Notes
The authors declare no competing financial interest. J
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation (22) Kukulin, V. I.; Krasnopol'sky, V. M. Description of few-body systems via analytic continuation in the coupling constant. J. Phys. A: Math. Gen. 1977, 10, L33. (23) Krasnopolsky, V. M.; Kukulin, V. I. Theory of resonance states based on analytic continuation in the coupling constant. Phys. Lett. A 1978, 69, 251−254. (24) Horácě k, J.; Mach, P.; Urban, J. Calculation of S-matrix poles by means of analytic continuation in the coupling constant: Application to the 2Πg state of N2−. Phys. Rev. A: At., Mol., Opt. Phys. 2010, 82, 032713. (25) Horácě k, J.; Paidarová, I.; Č urík, R. Determination of the resonance energy and width of the 2B2g shape resonance of ethylene with the Method of analytic continuation of the coupling constant. J. Phys. Chem. A 2014, 118, 6536. (26) Horácě k, J.; Paidarová, I.; Č urík, R. On a simple way to calculate electronic resonances for polyatomic molecules. J. Chem. Phys. 2015, 143, 184102. (27) Č urík, R.; Paidarová, I.; Horácě k, J. The 2Πg shape resonance of acetylene anion: An investigation with the RAC method. Eur. Phys. J. D 2016, 70, 146. (28) Sommerfeld, T.; Ehara, M. Short-range stabilizing potential for computing energies and lifetimes of temporary anions with extrapolation methods. J. Chem. Phys. 2015, 142, 034105. (29) White, A. F.; Head-Gordon, M.; McCurdy, C. W. Stabilizing potentials in bound state analytic continuation methods for electronic reonances in polyatomic molecules. J. Chem. Phys. 2017, 146, 044112. (30) Kossmann, S.; Kirchner, B.; Neese, F. Performance of modern density functional theory for the prediction of hyperfine structure: meta-GGA and double hybrid functionals. Mol. Phys. 2007, 105, 2049−2071. (31) Dunning, T. H., Jr. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007. (32) Neese, F. with Becker, U.; Bykov, D.; Ganyushin, D.; Hansen, A.; Izsak, R.; Liakos, D. G.; Kollmar, C.; Kossmann, S.; Pantazis, D. A.; Petrenko, T.; Reimann, C.; Riplinger, C.; Roemelt, M.; Sandhöfer, B.; Shapiro, I.; Sivalingam, K.; Wezisla, B. Orca, a package of programs. For the current version see https://orcaforum.cec.mpg.de/ (page last accessed January 2016). (33) Nakatsuji, H. Cluster Expansion of the Wavefunction. Electron Correlations in Ground and Excited States by SAC (SymmetryAdapted-Cluster) and SAC-CI Theories. Chem. Phys. Lett. 1979, 67, 329. (34) Nakatsuji, H. Cluster Expansion of the Wavefunction. Calculation of Electron Correlations in Ground and Excited States by SAC and SAC-CI Theories. Chem. Phys. Lett. 1979, 67, 334. (35) Nooijen, M.; Bartlett, R. J. Equation of motion coupled cluster method for electron attachment. J. Chem. Phys. 1995, 102, 3629. (36) Stanton, J. F.; Bartlett, R. J. The equation-of-motion, coupledcluster method. A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited-state properties. J. Chem. Phys. 1993, 98, 7029−7039. (37) Mertins, F.; Schirmer, J.; Tarantelli, A. Algebraic propagotor approaches and intermediate-state representations. II. The equationof-motion methods for N, N ± 1, and N ± 2 electrons. Phys. Rev. A: At., Mol., Opt. Phys. 1996, 53, 2153. (38) Ehara, M.; Ishida, M.; Toyota, K.; Nakatsuji, H. SAC-CI General-R Method: Theory and Applications to the Multi-Electron Processes. Reviews in Modern Quantum Chemistry; World Scientific: Singapore, 2002; pp 293−319, DOI: 10.1142/9789812775702_0011. (39) Fukuda, R.; Nakatsuji, H. Formulation and implementation of direct algorithm for the symmetry-adapted cluster and symmetryadapted cluster-configulation interaction method. J. Chem. Phys. 2008, 128, 094105. (40) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima,
T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09, Revision B.01; Gaussian: Wallingford, CT, USA, 2009. (41) Sommerfeld, T.; Ehara, M. Complex Absorbing Potentials with Voronoi Isosurfaces Wrapping Perfectly around Molecules. J. Chem. Theory Comput. 2015, 11, 4627. (42) Ehara, M.; Fukuda, R.; Sommerfeld, T. Projected CAP/SAC-CI method with smooth Voronoi potential for calculating resonance states. J. Comput. Chem. 2016, 37, 242. (43) Kanazawa, Y.; Ehara, M.; Sommerfeld, T. Low-lying π* Resonances of Standard and Rare DNA or RNA Bases Studied by the Projected CAP/SAC-CI Method. J. Phys. Chem. A 2016, 120, 1545. (44) Ehara, M.; Kanazawa, Y.; Sommerfeld, T. Low-lying π* Resonances Associated with Cyano Groups: A CAP/SAC-CI Study. Chem. Phys. 2017, 482, 169. (45) Nestmann, B. M.; Peyerimhoff, S. D. Calculation of the discrete component of resonance states in negative ions by variation of nuclear charges. J. Phys. B: At. Mol. Phys. 1985, 18, 615. (46) Lebedev, V. I.; Laikov, D. N. A quadrature formula for the sphere of the 131st algebraic order of accuracy. Dokl. Math. 1999, 59, 477. (47) Jordan, K. D. Construction of potential energy curves in avoided crossing situations. Chem. Phys. 1975, 9, 199. (48) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes, 3rd ed.; Cambridge University Press: New York, 2007; pp 487−55. (49) Light, J. C. Discrete Variable Representation in Quantum Dynamics. In Time-Dependent Quantum Molecular Dynamics; Broeckhove, J., Lathouwers, L., Eds.; Springer: New York, 1992; p 185, DOI: 10.1007/978-1-4899-2326-4. (50) Maier, C. H.; Cederbaum, L. S.; Domcke, W. A spherical-box approach to resonances. J. Phys. B: At. Mol. Phys. 1980, 13, L119. (51) Mandelshtam, V. A.; Ravuri, T. R.; Taylor, H. S. Calculation of the density of resonance states using the stabilization method. Phys. Rev. Lett. 1993, 70, 1932. (52) McCurdy, C. W.; McNutt, J. F. On the possibility of analytically continuing stabilization graphs to determine positions and widths accurately. Chem. Phys. Lett. 1983, 94, 306. (53) Löwdin, P.-O. Approximate calculation of lifetimes of resonance states in the continuum from real stabilisation graphs. Int. J. Quantum Chem. 1985, 27, 495. (54) Jagau, T.-C.; Zuev, D.; Bravaya, K. B.; Epifanovsky, E.; Krylov, A. I. A fresh look a resonances and complex absorbing potentials: Density matrix-based approach. J. Phys. Chem. Lett. 2014, 5, 310.
K
DOI: 10.1021/acs.jctc.6b01228 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX