Time-dependent coupled-cluster method with time-varying orbital functions, called time-dependent optimized coupled-cluster (TD-OCC) method, is formulated for multielectron dynamics in an intense laser field. We have successfully derived the equations of motion for CC amplitudes and orthonormal orbital functions based on the real action functional, and implemented the method including double excitations (TD-OCCD) and double and triple excitations (TD-OCCDT) within the optimized active orbitals. The present method is size extensive and gauge invariant, a polynomial cost-scaling alternative to the time-dependent multiconfiguration self-consistent-field method. The first application of the TD-OCC method of intense-laser driven correlated electron dynamics in Ar atom is reported.

The strong-field physics and ultrafast science are rapidly progressing toward a world-changing goal to directly measure and control the electron motion in matters.1–5 Reliable theoretical and computational methods are indispensable to understand and predict strong-field and ultrafast phenomena, which frequently involve both bound excitations and ionizations, where dynamical electron correlation plays a key role.

One of the most advanced methods to describe such electron dynamics is the multiconfiguration time-dependent Hartree-Fock (MCTDHF) method6–13 and more generally the time-dependent multiconfiguration self-consistent-field (TD-MCSCF) method.14,15 In TD-MCSCF, the electronic wavefunction is given by the configuration-interaction (CI) expansion, ΨCI(t) = ICI(tI(t), and both CI coefficients {CI(t)} and occupied spin-orbital functions {ψp(t)} constituting the Slater determinants {ΦI(t)} are propagated in time. (Hereafter spin-orbital functions are simply referred to as orbitals. See, e.g, Refs. 16–19 for TDCI methods using fixed orbitals and Ref. 15 for a broad review of ab initio wavefunction-based methods for multielectron dynamics.) Though powerful, the full CI-based MCTDHF method suffers from the factorial scaling of the computational cost with respect to the number of electrons, which restricts its applicability to the systems of moderate size.

The time-dependent complete-active-space self-consistent-field (TD-CASSCF) method20 significantly reduces the cost by flexibly classifying occupied orbitals into frozen-core, dynamical-core, and active orbitals as defined in Fig. 1. More approximate and thus less demanding models have been actively investigated,21–24 which rely on the truncated CI expansion within the active orbitals given, e.g, as

ΨCI=(C0+CiaÊia+CijabÊijab+CijkabcÊijkabc+)Φ,
(1)

where Êijkabc=ĉaĉbĉcĉkĉjĉi, with cμ and cμ being creation and annihilation operators, respectively, for ψμ. (Einstein summation convention is implied throughout for orbital indices.) The operator Êijkabc excites electron(s) from orbital(s) i, j, k, … occupied in the reference determinant Φ (hole) to those a, b, c, … not occupied in the reference (particle). These methods achieve a polynomial, instead of factorial, cost scaling, and state-of-the-art real-space implementations23,25–28 have proved their great utility; however, truncated-CI-based methods share a general drawback of not being size extensive.29 

FIG. 1.

Illustration of the orbital subspacing (for a spin-degenerate case). The up and down arrows represent electrons, and horizontal lines represent spatial orbitals, classified into frozen-core (doubly occupied and fixed in time) and dynamical-core (doubly occupied but propagated in time) orbitals and active orbitals (correlated and propagated). The active orbital space is further split into the hole and particle subspaces. The index convention for orbitals of each group is given in the parentheses.

FIG. 1.

Illustration of the orbital subspacing (for a spin-degenerate case). The up and down arrows represent electrons, and horizontal lines represent spatial orbitals, classified into frozen-core (doubly occupied and fixed in time) and dynamical-core (doubly occupied but propagated in time) orbitals and active orbitals (correlated and propagated). The active orbital space is further split into the hole and particle subspaces. The index convention for orbitals of each group is given in the parentheses.

Close modal

Therefore, one naturally seeks an alternative to Eq. (1), where the truncated CI is replaced with the size-extensive coupled-cluster (CC) expansion29–32 

ΨCC=exp(τiaÊia+τijabÊijab+τijkabcÊijkabc+)Φ,
(2)

where both excitation amplitudes {τijkabc(t)} and orbitals are propagated in time, just as the stationary orbital-optimized CC (OCC) method33–36 optimizes both excitation amplitudes and orbitals to minimize the CC energy functional. This idea was recently put forward by Kvaal,37 who, based on Arponen’s seminal work,38 developed a time-dependent coupled-cluster method using time-varying biorthonormal orbitals, called orbital-adapted time-dependent coupled-cluster (OATDCC) method, and applied it to the collision problem in a one-dimensional model Hamiltonian. This method, however, has a difficulty in obtaining the ground state of the system, required as an initial state of simulations, and has not yet been applied to real three-dimensional or time-dependent Hamiltonian. We also notice recent investigations on the time-dependent coupled-cluster methods using fixed orbitals.39,40 Their primary focus, however, has been placed on retrieving spectral information rather than the time-domain electron dynamics itself.

In this communication, we report our successful formulation and implementation of the time-dependent coupled-cluster method using time-varying orthonormal orbitals, called time-dependent optimized coupled-cluster (TD-OCC) method, and present its first application to electron dynamics of Ar atom subject to an intense laser field. The Hartree atomic units are used throughout unless otherwise noted.

We are interested in the dynamics of a many-electron system interacting with external electromagnetic fields governed by the Hamiltonian written in the second quantization as

Ĥ(t)=hνμ(t)Êνμ+12gνδμγÊνδμγ,
(3)
hνμ(t)=dxψμ*[h0+Vext(t)]ψν,
(4)
gνδμγ=dx1dx2ψμ*(x1)ψγ*(x2)r121ψν(x1)ψδ(x2),
(5)

where [x = (r, σ) labels spatial-spin coordinate.] h0 and Vext are the field-free one-electron Hamiltonian and the interaction of an electron and external fields, given, e.g., in the electric dipole approximation as Vext = E(t) · r in the length gauge or Vext = −iA(tr in the velocity gauge, with E(t) and A(t) = −tdtE(t′) being the electric field and vector potential, respectively.

Note that we include in Eq. (3) not only a given number No of occupied orbitals {ψp} but also Nvvirtual orbitals {ψα} which are orthogonal to {ψp}. The number N = No + Nv of all orbitals {ψμ} is equal to the number of basis functions to expand orbitals (Fig. 1), and in the basis-set limit (N), Ĥ of Eq. (3) is equivalent to the first-quantized Hamiltonian. The index convention given in Fig. 1 is used for orbitals.

Following the previous work,29,37,41 we begin with the coupled-cluster Lagrangian L,

L=ΨL(Ĥit)ΨR=Φ(1+Λ^)eT^(Ĥit)eT^Φ,
(6)
T^=τiaÊia+τijabÊijab+τijkabcÊijkabc+=τIAÊIA=T^1+T^2+T^3+,
(7)
Λ^=λaiÊai+λabijÊabij+λabcijkÊabcijk+=λAIÊAI=Λ^1+Λ^2+Λ^3+,
(8)

where Φ=iiiĉiĉiĉi is a reference determinant. For the notational brevity, we use composite indices I = ijk… and A = abc… to denote general-rank excitation (de-excitation) operators and amplitudes as ÊIA (ÊAI) and τIA (λAI). Note that we use the orthonormal orbitals {ψμ} and ⟨Φ| = |Φ⟩.

As the guiding principle, we adopt the manifestly real action functional,

S=Rt0t1Ldt=12t0t1L+L*dt,
(9)

and require the action to be stationary

2δS=δΨLĤΨR+ΨL|Ĥ|δΨRiδΨL|Ψ̇R+iΨ̇L|δΨR+c.c=0,
(10)

with respect to the variation of left- and right-state wavefunctions. This real-action formulation using orthonormal orbitals (also adopted in Ref. 41 for a gauge-invariant response theory) distinguishes our method from the OATDCC method,37 which is based on the bivariational principle with the complex-analytic action functional S = ∫dtL using biorthonormal orbitals.

The time derivative and variation of the left- and right-state wavefunctions can be written as20,22,24,42

Ψ̇R=(tT^+X^)eT^Φ,
(11a)
δΨR=(δT^+Δ^)eT^Φ,
(11b)
Ψ̇L=Φ[(1+Λ^)eT^(tT^+X^)+tΛ^eT^],
(11c)
δΨL=Φ[(1+Λ^)eT^(δT^+Δ^)+δΛ^eT^],
(11d)

where tT^=τ̇IAÊIA, δT^=δτIAÊIA, tΛ^=λ̇AIÊAI, δΛ^=δλAIÊAI, X^=XνμÊνμ, Δ^=ΔνμÊνμ, Xνμ=ψμ|ψ̇ν, and Δνμ=ψμ|δψν. The matrices X and Δ are anti-Hermitian to enforce the orthonormality of orbitals ψμ(t)|ψν(t)=δνμ at any time.

We insert Eqs. (11) into Eq. (10) and group terms for each kind of variation to obtain

2δS=δτIAΦ(1+Λ^)eT^[ĤiX^,ÊIA]eT^Φ+iλ̇AI+δλAIΦIAeT^(ĤiX^)eT^Φiτ̇IA+c.c+ΔνμΦ(1+Λ^)eT^[Ĥi(tT^)iX^,Êνμ]eT^Φ+iΦ(tΛ^)eT^ÊνμeT^ΦΦ(1+Λ^)eT^[Ĥi(tT^)iX^,Êμν]eT^Φ*+iΦ(tΛ^)eT^ÊμνeT^Φ*.
(12)

The equation of motion (EOM) for τIA is derived from δS/δλAI(t)=0 and that for λAI from δS/δτIA(t)=0 as

iτ̇IA=ΦIAeT^(ĤiX^)eT^Φ,
(13)
iλ̇AI=Φ(1+Λ^)eT^[ĤiX^,ÊIA]eT^Φ.
(14)

The EOMs for orbitals are to be obtained from δS/δΔνμ(t)=0. After straightforward yet tedious algebraic manipulations, we found that (i) orbital pairs within the same orbital space (with no further subclassification in Fig. 1) are redundant as well known for the stationary single-reference CC methods, i.e., Xji,Xji,Xji,Xba can be arbitrary anti-Hermitian matrix elements, and (ii) the expressions of all the other Xνμ terms except the hole-particle contributions Xia (and Xai=Xia*) are formally identical to those for the TD-CASSCF method.20 Thus we follow the discussion in Ref. 20 and write the final expression for the orbital EOM, with the hole-particle terms {Xia} left unspecified until the next section,

iψṗ=(1P^)F^pψp+iψqXpq,
(15)
F^pψp=ĥψp+ŴsrψqPorqs(D1)po,
(16)
Dqp=12(ρqp+ρpq*),Pqspr=12(ρqspr+ρprqs*),
(17)
ρqp=Φ(1+Λ^)eT^ÊpqeT^Φ,
(18)
ρqspr=Φ(1+Λ^)eT^ÊprqseT^Φ,
(19)

where P^=pψpψp, Ŵsr is the mean-field operator,20 and {Xqp} except {Xia} are given in Ref. 20. (Note that the Hermitian matrix RiX was used as the working variable in Ref. 20. See also Ref. 25 for the velocity gauge treatment of frozen-core orbitals.) We note the natural appearance of the Hermitialized reduced density matrices (RDMs) D and P is the direct consequence of relying on the real action functional with orthonormal orbitals.

In order to derive a working equation for the hole-particle rotations {Xia}, and thus complete the derivation of EOMs, one has to specify a particular ansatz of the CC Lagrangian, i.e., which terms to be included in Eqs. (7) and (8). We, therefore, define the TD-OCC method by the ansatz

T^=T^2+T^3+,Λ^=Λ^2+Λ^3+,
(20)

that is, with single amplitudes T^1 and Λ^1 omitted. The same form has been adopted in the ground-state OCC methods33–35 and OATDCC method. We have also tried to retain single amplitudes (T^=T^1+T^2+, Λ^=Λ^1+Λ^2+) and apply the time-dependent variational principle. However, the resultant EOMs were stable neither in real time nor in imaginary time propagation. A similar difficulty has been reported for the stationary OCC method,33 which is in essence due to the similar physical roles played by single excitations and hole-particle rotations.

For the ansatz of Eq. (20), we derive the system of equations to be solved for Xia [by coupling Eqs. (13) and (14) and δS/δΔai=0], as

  i(δbaDijDbaδij)Xjb+12AibajXjb*=Bia,
(21)
Aibaj=Φ(1+Λ^)eT^[ÊIA,Êbj]eT^ΦΦIAeT^ÊaieT^Φ  Φ(1+Λ^)eT^[ÊIA,Êai]eT^ΦΦIAeT^ÊbjeT^Φ,
(22)
Bia=FjaDijDbaFbi*  12Φ(1+Λ^)eT^[ÊIA,Ĥ]eT^ΦΦIAeT^ÊaieT^Φ  +12Φ(1+Λ^)eT^[ÊIA,Êai]eT^ΦΦIAeT^ĤeT^Φ.
(23)

A special simplification arises when only double amplitudes are included where the cluster operator eT^ consists only of even-order (double, quadruple, …) excitations, in which case Eq. (22) and the last two terms of Eq. (23) vanish, and X^ needs not be included in Eqs. (13) and (14).

In summary, the TD-OCC EOMs consist of Eq. (13) for T^2,T^3,, Eq. (14) for Λ^2,Λ^3,, and Eq. (15) for orbitals, with {Xia} obtained by solving Eq. (21).

The present TD-OCC method is, as the name suggests, a direct time-dependent extension of the stationary OCC method.33–35 

We have implemented the TD-OCC method with double excitations (TD-OCCD) and double and triple excitations (TD-OCCDT) for atom-centered Gaussian basis functions and the atomic Hamiltonian with finite-element discrete variable representation (FEDVR) basis.25,28 These codes share the same implementation for the basis-independent procedures [Eqs. (13), (14), (17)–(19), and (21)–(23)].

It turns out, in contrast to the OATDCC method, that the imaginary-time relaxation is feasible for the TD-OCC method to obtain the ground state, which is quite convenient when using grids or locally supported basis like FEDVR to expand orbitals. This is confirmed in Table I, which compares the total energies of the ground state of BH molecule with double-ζ plus polarization (DZP) basis44 obtained by imaginary time propagation with results of corresponding stationary methods.35 The OCC and CASSCF methods correlate all six electrons among six active orbitals. The agreement of OCCD energies up to the reported (six) decimal places in Ref. 35 clearly demonstrates our correct implementation and the feasibility of the imaginary-time method.

TABLE I.

Total energies (Hartree) of BH molecule (the bond length 2.4 Å) with the DZP basis set.

This workaReference 35 
HF −25.124 742 010 −25.124 742 
OCCD −25.178 285 704 −25.178 286 
OCCDT −25.178 312 565  
CASSCF −25.178 334 889 −25.178 335 
This workaReference 35 
HF −25.124 742 010 −25.124 742 
OCCD −25.178 285 704 −25.178 286 
OCCDT −25.178 312 565  
CASSCF −25.178 334 889 −25.178 335 
a

The overlap, one-electron, and two-electron repulsion integrals over Gaussian basis functions are generated using Gaussian09 program (Ref. 43), and used to propagate EOMs in imaginary time in the orthonormalized Gaussian basis, with a convergence threshold of 10−15 Hartree of energy difference in subsequent time steps.

Finally, we report our initial application of the TD-OCC method to Ar atom in a strong laser field. We used the spherical FEDVR basis given as the product Ylm(θ, ϕ)fk(r) of spherical harmonics Ylm with llmax and the FEDVR basis fk(r), which divides the range of radial coordinate 0 < r < Rmax into finite elements of length 4 (with several finer elements near the nuclei), each supporting 23 local DVR functions. We first obtained the ground states of TDHF, TD-OCC, and TD-CASSCF methods by the imaginary time relaxation. For the latter two methods, the Ne core was kept frozen at HF orbitals and eight valence electrons were correlated among 13 optimized active orbitals, with their initial characters being 3s3p3d4p4s. Starting from the ground states, we then simulated the electron dynamics of Ar atom subject to a linearly polarized (in the z direction) laser pulse with a wavelength of 800 nm, a peak intensity of 6 × 1014 W/cm2, and a foot-to-foot pulse duration of three optical cycles with a sin2-shaped envelope, within the dipole approximation in the velocity gauge. The EOMs are propagated using the fourth-order exponential Runge-Kutta method.45 The Ne core was frozen in all the methods and 13 active orbitals were propagated in TD-OCC and TD-CASSCF methods. The basis parameters were calibrated to achieve convergence with Rmax = 240 (with the mask absorbing boundary25) and lmax = 63.

Figure 2 shows high-harmonic generation (HHG) spectra calculated as the modulus squared of the Fourier transform of the expectation value of the dipole acceleration ⟨d2z^/dt2⟩(t) which, in turn, is obtained with the modified Ehrenfest expression25 using RDMs of Eq. (17). All methods predicted qualitatively similar spectra, with a minimum at 53 eV (≈34th order, indicated with an arrow) close to the Cooper minimum of photoionization spectrum,46 which was experimentally observed.47 The TDHF method, however, failed to reproduce fine structures of the TD-CASSCF spectrum especially at higher plateau as seen in the inset of Fig. 2. The description is largely improved by TD-OCCD, and the TD-OCCDT spectrum well reproduces the TD-CASSCF one in virtually all details.

FIG. 2.

The HHG spectra of Ar exposed to a laser pulse with a wavelength of 800 nm and an intensity of 6 × 1014 W/cm2. Comparison of the results of TDHF, TD-OCCD, TD-OCCDT, and TD-CASSCF methods. The inset shows a close-up of the spectra from 50th to 80th harmonic order.

FIG. 2.

The HHG spectra of Ar exposed to a laser pulse with a wavelength of 800 nm and an intensity of 6 × 1014 W/cm2. Comparison of the results of TDHF, TD-OCCD, TD-OCCDT, and TD-CASSCF methods. The inset shows a close-up of the spectra from 50th to 80th harmonic order.

Close modal

We also compare the probabilities of finding one [Fig. 3(a)] and two [Fig. 3(b)] electron(s) outside a sphere of radius R0 = 20, which measure the single and double ionization probabilities, respectively. These probabilities are much more sensitive to the description of correlated motions of electrons than HHG, which hinders a correct description with TDHF. The rapid convergence of TD-OCC results to the TD-CASSCF ones for such nonperturbatively nonlinear processes strongly promises the TD-OCC method to be a vital tool to investigate correlated high-field phenomena.

FIG. 3.

The probabilities, as a function of time, of finding one (a) and two (b) electrons outside a sphere of radius R0 = 20 a.u. Comparison of the results of TDHF, TD-OCCD, TD-OCCDT, and TD-CASSCF methods.

FIG. 3.

The probabilities, as a function of time, of finding one (a) and two (b) electrons outside a sphere of radius R0 = 20 a.u. Comparison of the results of TDHF, TD-OCCD, TD-OCCDT, and TD-CASSCF methods.

Close modal

We have successfully formulated a new time-dependent coupled-cluster method called the TD-OCC method and implemented the first two, and the most important, series of approximations, TD-OCCD and TD-OCCDT. The present method is size extensive and gauge invariant, a polynomial cost-scaling alternative to the TD-MCSCF method. It would open new possibilities of high-accuracy first-principle investigations of multielectron dynamics in ever-unreachable large target systems. The rigorous derivation, details of the implementation, as well as other ansatz for TD-CC theories and comparison thereof, will be presented in a forthcoming article.

This research was supported in part by a Grant-in-Aid for Scientific Research (Grants Nos. 25286064, 26390076, 26600111, 16H03881, and 17K05070) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan and also by the Photon Frontier Network Program of MEXT. This research was also partially supported by the Center of Innovation Program from the Japan Science and Technology Agency, JST, and by CREST (Grant No. JPMJCR15N1), JST. Y.O. gratefully acknowledges support from the Graduate School of Engineering, the University of Tokyo, Doctoral Student Special Incentives Program (SEUT Fellowship).

1.
M.
Protopapas
,
C. H.
Keitel
, and
P. L.
Knight
,
Rep. Prog. Phys.
60
,
389
(
1997
).
2.
P.
Agostini
and
L. F.
DiMauro
,
Rep. Prog. Phys.
67
,
813
(
2004
).
3.
F.
Krausz
and
M.
Ivanov
,
Rev. Mod. Phys.
81
,
163
(
2009
).
4.
L.
Gallmann
,
C.
Cirelli
, and
U.
Keller
,
Annu. Rev. Phys. Chem.
63
,
447
(
2013
).
5.
M.
Nisoli
,
P.
Decleva
,
F.
Calegari
,
A.
Palacios
, and
F.
Martín
,
Chem. Rev.
117
,
10760
(
2017
).
6.
J.
Zanghellini
,
M.
Kitzler
,
C.
Fabian
,
T.
Brabec
, and
A.
Scrinzi
,
Laser Phys.
13
,
1064
(
2003
).
7.
T.
Kato
and
H.
Kono
,
Chem. Phys. Lett.
392
,
533
(
2004
).
8.
J.
Caillat
,
J.
Zanghellini
,
M.
Kitzler
,
O.
Koch
,
W.
Kreuzer
, and
A.
Scrinzi
,
Phys. Rev. A
71
,
012712
(
2005
).
9.
M.
Nest
,
T.
Klamroth
, and
P.
Saalfrank
,
J. Chem. Phys.
122
,
124102
(
2005
).
10.
T.
Kato
and
H.
Kono
,
J. Chem. Phys.
128
,
184102
(
2008
).
11.
D.
Hochstuhl
and
M.
Bonitz
,
J. Chem. Phys.
134
,
084106
(
2011
).
12.
D. J.
Haxton
,
K. V.
Lawler
, and
C. W.
McCurdy
,
Phys. Rev. A
86
,
013406
(
2012
).
13.
D. J.
Haxton
and
C. W.
McCurdy
,
Phys. Rev. A
90
,
053426
(
2014
).
14.
T. T.
Nguyen-Dang
,
M.
Peters
,
S.-M.
Wang
,
E.
Sinelnikov
, and
F.
Dion
,
J. Chem. Phys.
127
,
174107
(
2007
).
15.
K. L.
Ishikawa
and
T.
Sato
,
IEEE J. Sel. Top. Quantum Electron.
21
,
8700916
(
2015
).
16.
N.
Rohringer
,
A.
Gordon
, and
R.
Santra
,
Phys. Rev. A
74
,
043420
(
2006
).
17.
S.
Pabst
,
Eur. Phys. J. Spec. Top.
221
,
1
(
2013
).
18.
D.
Hochstuhl
and
M.
Bonitz
,
Phys. Rev. A
86
,
053424
(
2012
).
19.
S.
Bauch
,
L. K.
Sørensen
, and
L. B.
Madsen
,
Phys. Rev. A
90
,
062508
(
2014
).
20.
T.
Sato
and
K. L.
Ishikawa
,
Phys. Rev. A
88
,
023402
(
2013
).
21.
H.
Miyagi
and
L. B.
Madsen
,
Phys. Rev. A
87
,
062511
(
2013
).
22.
H.
Miyagi
and
L. B.
Madsen
,
Phys. Rev. A
89
,
063416
(
2014
).
23.
D. J.
Haxton
and
C. W.
McCurdy
,
Phys. Rev. A
91
,
012509
(
2015
).
24.
T.
Sato
and
K. L.
Ishikawa
,
Phys. Rev. A
91
,
023417
(
2015
).
25.
T.
Sato
,
K. L.
Ishikawa
,
I.
Březinová
,
F.
Lackner
,
S.
Nagele
, and
J.
Burgdörfer
,
Phys. Rev. A
94
,
023405
(
2016
).
26.
R.
Sawada
,
T.
Sato
, and
K. L.
Ishikawa
,
Phys. Rev. A
93
,
023434
(
2016
).
27.
J. J.
Omiste
,
W.
Li
, and
L. B.
Madsen
,
Phys. Rev. A
95
,
053422
(
2017
).
28.
Y.
Orimo
,
T.
Sato
,
A.
Scrinzi
, and
K. L.
Ishikawa
, “
Implementation of infinite-range exterior complex scaling to the time-dependent complete-active-space self-consistent-field method
” (unpublished).
29.
T.
Helgaker
,
P.
Jørgensen
, and
J.
Olsen
,
Molecular Electronic-Structure Theory
(
Wiley
,
2002
).
30.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry
(
Dover
,
Mineola
,
1996
).
31.
H. G.
Kümmel
,
Int. J. Mod. Phys. B
17
,
5311
(
2003
).
32.
I.
Shavitt
and
R. J.
Bartlett
,
Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory
(
Cambridge University Press
,
2009
).
33.
G. E.
Scuseria
and
H. F.
Schaefer
, III
,
Chem. Phys. Lett.
142
,
354
(
1987
).
34.
C. D.
Sherrill
,
A. I.
Krylov
,
E. F. C.
Byrd
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
4171
(
1998
).
35.
A. I.
Krylov
,
C. D.
Sherrill
,
E. F. C.
Byrd
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
10669
(
1998
).
36.
A.
Köhn
and
J.
Olsen
,
J. Chem. Phys.
122
,
084116
(
2005
).
37.
S.
Kvaal
,
J. Chem. Phys.
136
,
194109
(
2012
).
39.
C.
Huber
and
T.
Klamroth
,
J. Chem. Phys.
134
,
054113
(
2011
).
40.
D. R.
Nascimento
and
A. E.
DePrince
, III
,
J. Chem. Theory Comput.
12
,
5834
(
2016
).
41.
T. B.
Pedersen
,
H.
Koch
, and
C.
Hättig
,
J. Chem. Phys.
110
,
8318
(
1999
).
42.
R. P.
Miranda
,
A. J.
Fisher
,
L.
Stella
, and
A. P.
Horsfield
,
J. Chem. Phys.
134
,
244101
(
2011
).
43.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
 et al, gaussian 09, Revision A.2,
Gaussian, Inc.
,
Wallingford
,
CT
,
2009
.
44.
R. J.
Harrison
and
N. C.
Handy
,
Chem. Phys. Lett.
95
,
386
(
1983
).
45.
M.
Hochbruck
and
A.
Ostermann
,
Acta Numer.
19
,
209
(
2010
).
46.
J. W.
Cooper
,
Phys. Rev.
128
,
681
(
1962
).
47.
H. J.
Wörner
,
H.
Niikura
,
J. B.
Bertrand
,
P. B.
Corkum
, and
D.
Villeneuve
,
Phys. Rev. Lett.
102
,
103901
(
2009
).