Time-dependent optimized coupled-cluster method with doubles and perturbative triples for first principles simulation of multielectron dynamics.

We report the formulation of a new, cost-effective approximation method in the time-dependent optimized coupled-cluster (TD-OCC) framework [T. Sato et al., J. Chem. Phys. 148, 051101 (2018)] for first-principles simulations of multielectron dynamics in an intense laser field. The method, designated as TD-OCCD(T), is a time-dependent, orbital-optimized extension of the "gold-standard" CCSD(T) method in the ground-state electronic structure theory. The equations of motion for the orbital functions and the coupled-cluster amplitudes are derived based on the real-valued time-dependent variational principle using the fourth-order Lagrangian. The TD-OCCD(T) is size extensive and gauge invariant, and scales as O(N 7) with respect to the number of active orbitals N. The pilot application of the TD-OCCD(T) method to the strong-field ionization and high-order harmonic generation from a Kr atom is reported in comparison with the results of the previously developed methods, such as the time-dependent complete-active-space self-consistent field (TD-CASSCF), TD-OCC with double and triple excitations (TD-OCCDT), TD-OCC with double excitations (TD-OCCD), and the time-dependent Hartree-Fock (TDHF) methods.

1 Introduction

Recent years witnessed unprecedented progress in laser technologies, which made it possible to observe the motions of electrons at the attosecond time scale (Itatani et al. (2004); Corkum and Krausz (2007); Krausz and Ivanov (2009); Baker et al. (2006)). On the other hand, various theoretical and numerical methods have been developed for interpreting, understanding, and predicting the experiments.

The multi-configuration time-dependent Hartree-Fock (MCTDHF) method (Caillat et al. (2005); Kato and Kono (2004); Nest et al. (2005); Haxton et al. (2011); Hochstuhl and Bonitz (2011)), and the time-dependent complete-active-space self-consistent-field (TD-CASSCF) method (Sato and Ishikawa (2013); Sato et al. (2016); Sato et al. (2018a)) are the most rigorous approaches to solve time-dependent Schrödinger equation (TDSE) of many-electron systems, where the wavefunction is given by the full configuration interaction (FCI) expansion, with both CI coefficients {C (t)} and orbital functions {ψ (t)} constituting Slater determinants {Φ (t)} are propagated in time according to the time-dependent variational principle (TDVP). The TD-CASSCF method broadens the applicability of the MCTDHF method by flexibly classifying the orbital subspace into frozen-core, dynamical-core, and active. Unfortunately, the factorial computational scaling impedes large-scale applications. There are reports of various affordable size-inextensive methods (Miyagi and Madsen (2013, 2014); Haxton and McCurdy (2015); Sato and Ishikawa (2015)) developed by limiting the CI expansion of the wavefunction. Alternatively, the size-extensive coupled-cluster method, which relies on an exponential wavefunction, is a superior choice to address these problems with a polynomial cost-scaling (Kümmel (2003); Shavitt and Bartlett (2009)). We have developed an explicitly time-dependent coupled-cluster method considering optimized orthonormal orbitals within the flexibly chosen active space, called the time-dependent optimized coupled-cluster (TD-OCC) method, (Sato et al. (2018b)) including double (TD-OCCD) and double and triple excitation amplitudes (TD-OCCDT). Our method is a time-dependent formulation of the stationary optimized coupled-cluster method (Scuseria and Schaefer (1987); Sherrill et al. (1998); Krylov et al. (1998)). Kvaal (Kvaal (2012)) also developed an orbital adaptive time-dependent coupled-cluster (OATDCC) method using biorthogonal orbitals. We take note of a few reports on the time-dependent coupled-cluster methods (Huber and Klamroth (2011); Pigg et al. (2012); Nascimento and DePrince (2016)), using time-independent orbitals, and their interpretation (Pedersen and Kvaal (2019); Pedersen et al. (2021)), including the very initial attempts (Schonhammer (1978); Hoodbhoy and Negele (1978, 1979)).

The TD-OCCDT scales as O(N 8) (N= the number of active orbitals), not ideally suited for applications to larger chemical systems. Therefore, we have developed a few lower cost methods in the TD-OCC framework (Pathak et al. (2020b,c,a, 2021)). We find triple excitations are necessary, including perfect optimization of the orbitals. Therefore, we are interested in developing affordable TD-OCC methods retaining a part of the triples. The most popular coupled-cluster method that treats the triple excitation amplitudes approximately is called CCSD(T) (Raghavachari et al. (1989); Watts et al. (1993)). Bozkaya et al, (Bozkaya and Schaefer (2012)) included various symmetric and asymmetric triple excitation corrections to their optimized double (OD) method.

In this communication, we report the formulation and implementation of the CCSD(T) method in the time-dependent optimized coupled-cluster framework, TD-OCCD(T). Following our previous works (Sato et al. (2018b); Pathak et al. (2020b,c, 2021)), we exclude single excitation amplitudes but optimize the orbitals according to time-dependent variational principle (TDVP). As the first application of this method, we study electron dynamics in Kr using intense near-infrared laser fields.

2 Methods

The second quantization representation of the Hamiltonian, including the laser field, is as follows, where represents a creation (annihilation) operator in a complete, orthonormal set of 2n bas time-dependent spin-orbitals {ψ (t)}. n bas is the number of basis functions used for expanding the spatial part of ψ , which, in the present real-space implementation, corresponds to the number of grid points, and where x = ( , σ ) represents a composite spatial-spin coordinate. h 0 is the field free one-electronic Hamiltonian and V = A(t)p in the velocity gauge, A(t) = −∫ E(t′)dt′ is the vector potential, with E(t) being the laser electric field linearly polarized along the z axis.

The complete set of 2n bas spin-orbitals (labeled with μ, ν, γ, λ) is divided into n occ occupied (o, p, q, r, s) and 2n bas − n occ virtual spin-orbitals. The coupled-cluster (or CI) wavefunction is constructed only with occupied spin-orbitals, which are time-dependent in general, and virtual spin-orbitals form the orthogonal complement of the occupied spin-orbital space. The occupied spin-orbitals are classified into n core core spin-orbitals, which are occupied in the reference Φ and kept uncorrelated, and N = n occ − n core active spin-orbitals (t, u, v, w) among which the active electrons are correlated. The active spin-orbitals are further split into those in the hole space (i, j, k, l) and the particle space (a, b, c, d), which are defined as those occupied and unoccupied, respectively, in the reference Φ. The core spin-orbitals can further be split into frozen-core space (i′′, j′′), fixed in time and the dynamical-core space (i′, j′), propagated in time (Sato and Ishikawa (2013)) (See. Figure 1 in Sato et al. (2018b) for a pictorial illustration).

FIGURE 1

Time evolution of dipole moment of Kr irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 2 × 1014 W/cm2 calculated with TDHF, TD-OCCD, TD-OCCD(T), and TD-CASSCF methods.

The real action formulation of the TDVP with orthonormal orbitals is our guiding principle, (Sato et al. (2018b)) where are (de-)excitation amplitudes, and . The stationary conditions, δS = 0, with respect to the variation of the parameters of the wavefunction ( , , and δψ ) gives us the corresponding equations of motions (EOMs), δψ is orthonormality-conserving orbital variation.

For deriving the TD-OCCD(T) method, we first construct a fourth-order Lagrangian defined in Pathak et al. (2021). We make a further approximation to the Lagrangian and write separating it into two parts, where , , , and , , and is anti-Hermitian. The double amplitudes are obtained by making of Eq. 9a stationary with respect to , , the triples by making Eq. 9b stationary with respect to , and , where p(μν) and p(μ|νγ) are the permutation operators; p(μν)A = A − A , and p(μ/νγ) = 1 − p(μν) − p(μγ).

The EOM for the orbitals can be written down in the following form Sato et al. (2016), where is the identity operator within the space spanned by the given basis, is the projector onto the occupied spin-orbital space, and where D and P are Hermitialized one- (1RDM) and two- (2RDM) particle reduced density matrices defined in Sato et al. (2018b), and is the mean-field operator (Sato and Ishikawa (2013)). The matrix element includes orbital rotations among various subspaces. Non-redundant orbital rotations are determined by . Redundant orbital rotations , , and can be arbitrary antiHermitian matrix elements. The general expressions for the RDMs are the same as in the TD-OCCDT(4) method (Pathak et al. (2021)).

3 Numerical results and discussion

Our numerical implementation has an interface with the Gaussian09 program (Frisch et al. (2009)) for checking ground state energy with the standard Gaussian basis results. We study BH molecule with double-ζ plus polarization (DZP). We have reported ground state energy computed by propagating in the imaginary time for OCCD and OCCD(T) methods in Table 1 and compared those with the optimized double and asymmetric triple excitation corrections for the orbital-optimized doubles method of Bozkaya et al., Bozkaya and Schaefer (2012). We also compare our OCCD ground state energy result with Krylov et al.,Krylov et al. (1998) within the chosen active space of six electrons correlated among the six optimized active orbitals. We obtained a perfect agreement for all available values.

TABLE 1

Comparison of the ground state energy of BH (r =2.4 bohr) molecule in DZP basis .

Method	This work	References
OCCD b	− 25.225 591 67	− 25.225 592	Bozkaya and Schaefer (2012)
OCCD(T) b	− 25.226 913 29	− 25.226 913	Bozkaya and Schaefer (2012)
OCCD c	− 25.178 285 70	− 25.178 286	Krylov et al. (1998)
OCCD(T) c	− 25.178 301 00

Gaussian09 program (Frisch et al. (2009)) is used to generate the required one-electron, two-electron, and overlap integrals, required for the imaginary time propagation of EOMs in the orthonormalized Gaussian basis. A convergence cut-off of 10–15 Hartree of energy difference is chosen in subsequent time steps.

Six electrons correlated within the full basis set.

Six electrons correlated within the six optimized active orbitals.

We have used a spherical-finite-element-discrete-variable representation (FEDVR) basis for representing orbital functions, Sato et al. (2016); Orimo et al. (2018) where Y and f (r) are spherical harmonics and the normalized radial-FEDVR basis function, respectively. The expansion of the spherical harmonics continued up to the maximum angular momentum L max, and the radial FEDVR basis supports the range of radial coordinate 0 ≤ r ≤ R max, with cos1/4 mask function used as an absorbing boundary for avoiding unphysical reflection from the wall of the simulation box. We have used l max = 72, and the FEDVR basis supporting the radial coordinate 0 < r < 300 using 78 finite elements each containing 25 DVR functions. The absorbing boundary is switched on at r = 180 in all our simulations. The Fourth-order exponential Runge-Kutta method (Hochbruck and Ostermann (2010)) is used to propagate the EOMs with 20000 time steps for each optical cycle. We run the simulations for a further 6,000 time steps after the end of the pulse. In all correlation calculations, eight electrons of 4s4p orbitals are considered as active and correlated among thirteen active orbitals. We report simulation results computed using a three-cycle laser pulse with a central wavelength of 800 nm having intensity 2 × 1014 W/cm2 and a period of T = 2π/ω 0 ∼ 2.67 fs.

We report the time evolution of dipole moment of Kr in Figure 1 and in Figure 2 single electron ionization probability. Time-dependent dipole moment is evaluated as a trace using 1RDMs. For the single electron ionization probability, we computed the probability of finding an electron outside a sphere of a radius of 20 a.u. using RDMs defined in Refs. 19; 20; 37. We compare the results of TD-CASSCF, TD-OCCD(T), TD-OCCD, and TDHF methods.

FIGURE 2

Time evolution of single ionization probability of Kr irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 2 × 1014 W/cm2 calculated with TDHF, TD-OCCD, TD-OCCD(T), and TD-CASSCF methods.

We observe a substantial underestimation (both in Figure 1, and Figure 2) by the TDHF method due to the lack of correlation treatment. All correlation methods perform according to their ability to treat electron correlation. We also computed results using the TD-OCCDT method but not reported here since those results are not identifiable from the TD-CASSCF results within the graphical resolution.

Next, we report high-harmonic generation in Figure 3. It is calculated by squaring the modulus I(ω) = |a(ω)|2 of the Fourier transform of the expectation value of the dipole acceleration with a modified Ehrenfest expression (Sato et al. (2016)). In panel (c) of Figure 3, we plot the absolute relative deviation (δ(ω), of the spectral amplitude a(ω) from the TD-CASSCF value for each method. All methods qualitatively predict similar HHG spectra with TDHF underestimates the spectral intensity. The relative deviation of results from TD-CASSCF ones follows the general trend TDHF TD-OCCD TD-OCCD(T) TD-OCCDT, the same as what we observe for the time-dependent dipole moment and single ionization probability. We also simulated results with lower and higher intensity. However, the trend remains the same.

FIGURE 3

The HHG spectra (A,B) and the relative deviation (C) of the spectral amplitude from the TD-CASSCF spectrum from Kr irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 2 × 1014 W/cm2 with various methods.

Finally, we make a tally of computational costs for all the methods considered in this article. All simulations performed using an Intel(R) Xeon(R) Gold 6,230 central processing unit (CPU) with 40 processors with a clock speed of 2.10 GHz, and report total simulations time in Table 2. Further, we report a reduction in the computational cost for various TD-OCC methods relative to the TD-CASSCF. We see a massive 63% cost reduction for the TD-OCCD(T) method, which is larger than for the TD-OCCDT method (58%), and a minimal increase from the TD-OCCD method.

TABLE 2

Comparison of the total simulation time (in min) spent for TD-CASSCF, TD-OCCDT, TDCCD(T), and TD-OCCD methods.

Method	Time (min)	Cost reduction (%)
TD-CASSCF	47303	…
TD-OCCDT	19697	58
TD-OCCD(T)	17504	63
TD-OCCD	17494	63

Time spent for the simulation of Kr atom for 66000 time steps (0 ≤ t ≤ 3.3T) of a real-time simulation (I 0 = 2 × 1014 W/cm2 and λ = 800 nm), using an Intel(R) Xeon(R) Gold 6230 CPU with 40 processors having a clock speed of 2.10 GHz.

4 Concluding remarks

We have reported the formulation and implementation of the TD-OCCD(T) method. As the first application, we employed this method to study laser-driven dynamics in Kr exposed to an intense near-infrared laser pulse. We observe a 63% cost reduction in comparison to the TD-CASSCF method without losing much accuracy. Therefore, we conclude that TD-OCCD(T) method will certainly be beneficial in exploring highly accurate ab initio simulations of electron dynamics in larger chemical systems.