Cholesky Decomposition Techniques in Electronic Structure Theory

Time-Dependent Multilevel Density Functional Theory

We present a novel three-layer approach based on multilevel density functional theory (MLDFT) and polarizable molecular mechanics to simulate the electronic excitations of chemical systems embedded in an external environment within the time-dependent DFT formalism. In our method, the electronic structure of a target system, the chromophore, is determined in the field of an embedded inactive layer, which is treated as frozen. Long-range interactions are described by employing the polarizable fluctuating charge (FQ) force field. The resulting MLDFT/FQ thus accurately describes both electrostatics (and polarization) and non-electrostatic target-environment interactions. The robustness and reliability of the approach are demonstrated by comparing our results with experimental data reported for various organic molecules in solution.

Accelerating Pythonic Coupled-Cluster Implementations: A Comparison Between CPUs and GPUs

In this work, we benchmark several Python routines for time and memory requirements to identify the optimal choice of the tensor contraction operations available. We scrutinize how to accelerate the bottleneck tensor operations of Pythonic coupled-cluster implementations in the Cholesky linear algebra domain, utilizing a NVIDIA Tesla V100S PCIe 32GB (rev 1a) graphics processing unit (GPU). The NVIDIA compute unified device architecture API interacts with CuPy, an open-source library for Python, designed as a NumPy drop-in replacement for GPUs. Due to the limitations of video memory, the GPU calculations must be performed batch-wise. Timing results of some contractions containing large tensors are presented. The CuPy implementation leads to a factor of 10-16 speed-up of the bottleneck tensor contractions compared to computations on 36 central processing unit (CPU) cores. Finally, we compare example CCSD and pCCD-LCCSD calculations performed solely on CPUs to their CPU-GPU hybrid implementation, which leads to a speed-up of a factor of 3-4 compared to the CPU-only variant.

Galván I, Alavi A, Aleotti F, Aquilante F, Autschbach J, Avagliano D, Baiardi A, Bao JJ, Battaglia S, Birnoschi L, Blanco-González A, Bokarev SI, Broer R, Cacciari R, Calio PB, Carlson RK, Carvalho Couto R, Cerdán L, Chibotaru LF, Chilton NF, Church JR, Conti I, Coriani S, Cuéllar-Zuquin J, Daoud RE, Dattani N, Decleva P, de Graaf C, Delcey M, De Vico L, Dobrautz W, Dong SS, Feng R, Ferré N, Filatov(Gulak) M, Gagliardi L, Garavelli M, González L, Guan Y, Guo M, Hennefarth MR, Hermes MR, Hoyer CE, Huix-Rotllant M, Jaiswal VK, Kaiser A, Kaliakin DS, Khamesian M, King DS, Kochetov V, Krośnicki M, Kumaar AA, Larsson ED, Lehtola S, Lepetit MB, Lischka H, López Ríos P, Lundberg M, Ma D, Mai S, Marquetand P, Merritt ICD, Montorsi F, Mörchen M, Nenov A, Nguyen VHA, Nishimoto Y, Oakley MS, Olivucci M, Oppel M, Padula D, Pandharkar R, Phung QM, Plasser F, Raggi G, Rebolini E, Reiher M, Rivalta I, Roca-Sanjuán D, Romig T, Safari AA, Sánchez-Mansilla A, Sand AM, Schapiro I, Scott TR, Segarra-Martí J, Segatta F, Sergentu DC, Sharma P, Shepard R, Shu Y, Staab JK, Straatsma TP, Sørensen LK, Tenorio BNC, Truhlar DG, Ungur L, Vacher M, Veryazov V, Voß TA, Weser O, Wu D, Yang X, Yarkony D, Zhou C, Zobel JP, Lindh R. The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry

The developments of the open-source OpenMolcas chemistry software environment since spring 2020 are described, with a focus on novel functionalities accessible in the stable branch of the package or via interfaces with other packages. These developments span a wide range of topics in computational chemistry and are presented in thematic sections: electronic structure theory, electronic spectroscopy simulations, analytic gradients and molecular structure optimizations, ab initio molecular dynamics, and other new features. This report offers an overview of the chemical phenomena and processes OpenMolcas can address, while showing that OpenMolcas is an attractive platform for state-of-the-art atomistic computer simulations.

Static embedding with pair coupled cluster doubles based methods

Quantum embedding methods have recently been significantly developed to model large molecular structures. This work proposes a novel wave function theory in a density functional theory (WTF-in-DFT) embedding scheme based on pair-coupled cluster doubles (pCCD)-type methods. While pCCD can reliably describe strongly-correlated systems with mean-field-like computational cost, the large extent of the dynamic correlation can be accounted for by (linearized) coupled-cluster corrections on top of the pCCD wave function. Here we focus on the linearized coupled-cluster singles and doubles (LCCSD) ansatz for electronic ground states and its extension to excited states within the equation of motion (EOM) formalism. We test our EOM-pCCD-LCCSD-in-DFT approach for the vertical excitation energies of the hydrogen-bonded water-ammonia complex, micro-solvated thymine, and uranyl tetrahalides (UO2X42-, X = F, Cl, Br). Furthermore, we assess the quality of the embedding potential using an orbital entanglement and correlation analysis. The approximate embedding models successfully capture changes in the excitation energies going from bare fragments to supramolecular structures and represent a promising computational method for excited states in large molecular systems.

A Multiconfigurational Wave Function Implementation of the Frenkel Exciton Model for Molecular Aggregates

We present an implementation of the Frenkel exciton model into the OpenMolcas program package enabling calculations of collective electronic excited states of molecular aggregates based on a multiconfigurational wave function description of the individual monomers. The computational protocol avoids using diabatization schemes and, thus, supermolecule calculations. Additionally, the use of the Cholesky decomposition of the two-electron integrals entering pair interactions enhances the efficiency of the computational scheme. The application of the method is exemplified for two test systems, that is, a formaldehyde oxime and a bacteriochlorophyll-like dimer. For the sake of comparison with the dipole approximation, we restrict our considerations to situations where intermonomer exchange can be neglected. The protocol is expected to be beneficial for aggregates composed of molecules with extended π systems, unpaired electrons such as radicals or transition metal centers, where it should outperform widely used methods based on time-dependent density functional theory.

The Multistate Quantum Monte Carlo Algebraic Diagrammatic Construction Method

A multistate formulation of the recently developed quantum Monte Carlo (QMC) algebraic diagrammatic construction (ADC) method, QMCADC, is presented. QMCADC solves the Hermitian eigenvalue problem of the second-order ADC scheme for the polarization propagator stochastically by combining ADC schemes with projector quantum Monte Carlo (PQMC). It allows for massively parallel distributed computing and exploits the sparsity of the effective ADC matrix, thereby relaxing memory and processing requirements of ADC methods significantly. Here, the theory and implementation of the multistate variant of QMCADC are described, and our first proof-of-principle calculations for various molecular systems are shown. Indeed, multistate QMCADC enables sampling of an arbitrary number of low-lying excited states and can reproduce their vertical excitation energies with a marginal controllable error. The performance of multistate QMCADC is examined in terms of state-wise and overall accuracy as well as with respect to the balance in the treatments of excited states relatively to each other. The results are very promising as they show bias and imbalances among excited states to diminish as the number of sampling points increases. Furthermore, the impact of the quality of trial wave functions on the vertical excitation energies is investigated. A black-box approach for the generation of high quality trial wave functions internally is given.

Reliably assessing the electronic structure of cytochrome P450 on today's classical computers and tomorrow's quantum computers

Chemical simulation is one of the most promising applications for future quantum computers. It is thought that quantum computers may enable accurate simulation for complex molecules that are otherwise impossible to simulate classically; that is, it displays quantum advantage. To better understand quantum advantage in chemical simulation, we explore what quantum and classical resources are required to simulate a series of pharmaceutically relevant molecules. Using classical methods, we show that reliable classical simulation of these molecules requires significant resources and therefore is a promising candidate for quantum simulation. We estimate the quantum resources, both in overall simulation time and the size. The insights from this study pave the way for future quantum simulation of complex molecules.

An accurate assessment of how quantum computers can be used for chemical simulation, especially their potential computational advantages, provides important context on how to deploy these future devices. To perform this assessment reliably, quantum resource estimates must be coupled with classical computations attempting to answer relevant chemical questions and to define the classical algorithms simulation frontier. Herein, we explore the quantum computation and classical computation resources required to assess the electronic structure of cytochrome P450 enzymes (CYPs) and thus define a classical–quantum advantage boundary. This is accomplished by analyzing the convergence of density matrix renormalization group plus n-electron valence state perturbation theory (DMRG+NEVPT2) and coupled-cluster singles doubles with noniterative triples [CCSD(T)] calculations for spin gaps in models of the CYP catalytic cycle that indicate multireference character. The quantum resources required to perform phase estimation using qubitized quantum walks are calculated for the same systems. Compilation into the surface code provides runtime estimates to compare directly to DMRG runtimes and to evaluate potential quantum advantage. Both classical and quantum resource estimates suggest that simulation of CYP models at scales large enough to balance dynamic and multiconfigurational electron correlation has the potential to be a quantum advantage problem and emphasizes the important interplay between classical computations and quantum algorithms development for chemical simulation.

Fragment Localized Molecular Orbitals

We introduce the concept of fragment localized molecular orbitals (FLMOs), which are Hartree-Fock molecular orbitals localized in specific fragments constituting a molecular system. In physical terms, we minimize the local electronic energies of the different fragments, at the cost of maximizing the repulsion between them. To showcase the approach, we rationalize the main interactions occurring in large biological systems in terms of interactions between the fragments of the system. In particular, we study an anticancer drug intercalated within DNA and retinal in anabaena sensory rhodopsin as prototypes of molecular systems embedded in biological matrixes. Finally, the FLMOs are exploited to rationalize the formation of two oligomers, prototypes of amyloid diseases, such as Parkinson and Alzheimer.

Computation of NMR shieldings at the CASSCF level using gauge-including atomic orbitals and Cholesky decomposition

We present an implementation of coupled-perturbed complete active space self-consistent field (CP-CASSCF) theory for the computation of nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals and Cholesky decomposed two-electron integrals. The CP-CASSCF equations are solved using a direct algorithm where the magnetic Hessian matrix-vector product is expressed in terms of one-index transformed quantities. Numerical tests on systems with up to about 1300 basis functions provide information regarding both the computational efficiency and limitations of our implementation.

Efficient Adiabatic Connection Approach for Strongly Correlated Systems: Application to Singlet-Triplet Gaps of Biradicals

Strong electron correlation can be captured with multireference wave function methods, but an accurate description of the electronic structure requires accounting for the dynamic correlation, which they miss. In this work, a new approach for the correlation energy based on the adiabatic connection (AC) is proposed. The AC_n method accounts for terms up to order n in the coupling constant, and it is size-consistent and free from instabilities. It employs the multireference random phase approximation and the Cholesky decomposition technique, leading to a computational cost growing with the fifth power of the system size. Because of the dependence on only one- and two-electron reduced density matrices, AC_n is more efficient than existing ab initio multireference dynamic correlation methods. AC_n affords excellent results for singlet-triplet gaps of challenging organic biradicals. The development presented in this work opens new perspectives for accurate calculations of systems with dozens of strongly correlated electrons.

Automated Generation of Optimized Auxiliary Basis Sets for Long-Range-Corrected TDDFT Using the Cholesky Decomposition

Range-separated hybrid functionals making use of a smooth separation of the Coulomb operator in terms of the error function and its complement have proven to be a valuable tool for improving Kohn-Sham density functional theory (DFT) calculations. This holds in particular for obtaining accurate excitation energies from linear-response time-dependent DFT. Evaluating the long-range exchange contributions represents one of the most time-consuming tasks in such calculations. Prefitted auxiliary basis sets can be employed to speed up this step. Here, we present a way to generate auxiliary basis sets optimized to fit the long-range exchange contributions only, contrary to the common optimization strategies on the basis of the full Coulomb operator. For this purpose, we use the atomic Cholesky decomposition technique. The basis sets are generated on-the-fly using the specific range-separation parameter defined in the exchange-correlation functional. We obtain excitation energies and oscillator strengths which are of similar or better accuracy than those obtained with conventional resolution-of-the-identity auxiliary basis sets while drastically reducing the number of auxiliary functions required. This is demonstrated for the QUESTDB#5 benchmark set. In addition, we outline the benefits of this approach in sequences of calculations employing varying range-separation parameters, as is the case in the optimally tuned range-separation strategy. Finally, we illustrate the efficiency of this approach for real-world examples, namely, a chlorophyll tetramer from photosystem II and a carotenoid-porphyrin-C₆₀ triad.

Toward the Minimal Floating Operation Count Cholesky Decomposition of Electron Repulsion Integrals

As quantum chemistry calculations deal with molecular systems of increasing size, the memory requirement to store electron-repulsion integrals (ERIs) greatly outpaces the physical memory available in computing hardware. The Cholesky decomposition of ERIs provides a convenient yet accurate technique to reduce the storage requirement of integrals. Recent developments of a two-step algorithm have drastically reduced the memory operation (MOP) count, leaving the floating operation (FLOP) count as the last frontier of cost reduction in the Cholesky ERI algorithm. In this report, we introduce a dynamic integral tracking, reusing, and compression/elimination protocol embedded in the two-step Cholesky ERI method. Benchmark studies suggest that this technique becomes particularly advantageous when the basis set consists of many computationally expensive high-angular-momentum basis functions. With this dynamic-ERI improvement, the Cholesky ERI approach proves to be a highly efficient algorithm with minimal FLOP and MOP count.

Multilevel Density Functional Theory

Following recent developments in multilevel embedding methods, we introduce a novel density matrix-based multilevel approach within the framework of density functional theory (DFT). In this multilevel DFT, the system is partitioned in an active and an inactive fragment, and all interactions are retained between the two parts. The decomposition of the total system is performed upon the density matrix. The orthogonality between the two parts is maintained by solving the Kohn–Sham equations in the MO basis for the active part only, while keeping the inactive density matrix frozen. This results in the reduction of computational cost. We outline the theory and implementation and discuss the differences and similarities with state-of-the-art DFT embedding methods. We present applications to aqueous solutions of methyloxirane and glycidol.

eT 1.0: An open source electronic structure program with emphasis on coupled cluster and multilevel methods

The e^T program is an open source electronic structure package with emphasis on coupled cluster and multilevel methods. It includes efficient spin adapted implementations of ground and excited singlet states, as well as equation of motion oscillator strengths, for CCS, CC2, CCSD, and CC3. Furthermore, e^T provides unique capabilities such as multilevel Hartree-Fock and multilevel CC2, real-time propagation for CCS and CCSD, and efficient CC3 oscillator strengths. With a coupled cluster code based on an efficient Cholesky decomposition algorithm for the electronic repulsion integrals, e^T has similar advantages as codes using density fitting, but with strict error control. Here, we present the main features of the program and demonstrate its performance through example calculations. Because of its availability, performance, and unique capabilities, we expect e^T to become a valuable resource to the electronic structure community.

Curing basis set overcompleteness with pivoted Cholesky decompositions

The description of weakly bound electronic states is especially difficult with atomic orbital basis sets. The diffuse atomic basis functions that are necessary to describe the extended electronic state generate significant linear dependencies in the molecular basis set, which may make the electronic structure calculations ill-convergent. We propose a method where the overcomplete molecular basis set is pruned by a pivoted Cholesky decomposition of the overlap matrix, yielding an optimal low-rank approximation that is numerically stable, the pivot indices determining a reduced basis set that is complete enough to describe all the basis functions in the original overcomplete basis. The method can be implemented either by a simple modification to the usual canonical orthogonalization procedure, which hides the excess functions and yields fewer efficiency benefits, or by generating custom basis sets for all the atoms in the system, yielding significant cost reductions in electronic structure calculations. The pruned basis sets from the latter choice allow accurate calculations to be performed at a lower cost even at the self-consistent field level, as illustrated on a solvated (H₂O)₂₄ ^- anion. Our results indicate that the Cholesky procedure allows one to perform calculations with accuracies close to standard augmented basis sets with cost savings which increase with the size of the basis set, ranging from 9% fewer functions in single-ζ basis sets to 28% fewer functions in triple-ζ basis sets.

Theoretical Investigation of Excited-State Intramolecular Proton Transfer and Photoisomerization of 2-(Iminomethyl)phenol

Time-dependent density functional theory and high-level ab initio calculations were performed to investigate the excited-state intramolecular proton transfer (ESIPT) and subsequent isomerization of 2-(iminomethyl)phenol (IMP). According to the results of the correlated theoretical methods, ESIPT is a barrierless process; subsequently, the isomerization (rotation of the torsion angle) of IMP also readily occurs. Fictitious intermediates are found due to an insufficient theoretical level. The molecular structure of the conical intersection (CI) during the isomerization process is optimized, and its branching plane is characterized. Both the gradient difference vector and the derivative coupling vector are significantly correlated to the C═O and H₂N-C antiparallel stretching coordinates, and the dynamic electron correlation effect is crucial to optimize the molecular structure of the real CI of IMP. The relaxation pathway from the CI in the S₀ state was examined; the dominant pathway proceeds to the trans-keto form of IMP. However, if the C═O and H₂N-C antiparallel stretching mode is sufficiently populated, then the reaction proceeds to the cis-keto form of IMP and can eventually recover to the cis-enol form.

Low-scaling first-order properties within second-order Møller-Plesset perturbation theory using Cholesky decomposed density matrices

An efficient implementation of energy gradients and of hyperfine coupling constants in second-order Møller-Plesset perturbation theory (MP2) is presented based on our fully atomic orbital (AO)-based approach. For the latter, an unrestricted AO-based MP2 formulation is introduced. A reduction in the dependency of the computational efficiency on the size of the basis set is achieved by a Cholesky decomposition and the prefactor is reduced by the resolution-of-the-identity approximation. Significant integral contributions are selected based on distance-including integral estimates (denoted as QQR-screening) and its reliability as a fully controlled screening procedure is demonstrated. The rate-determining steps are shown via model computations to scale cubically in the computation of energy gradients and quadratically in the case of hyperfine coupling constants. Furthermore, a significant speed-up of the computational time with respect to the canonical formulation is demonstrated.

A full-pivoting algorithm for the Cholesky decomposition of two-electron repulsion and spin-orbit coupling integrals

A significant reduction in the computational effort for the evaluation of the electronic repulsion integrals (ERI) in ab initio quantum chemistry calculations is obtained by using Cholesky decomposition (CD), a numerical procedure that can remove the zero or small eigenvalues of the ERI positive (semi)definite matrix, while avoiding the calculation of the entire matrix. Conversely, due to its antisymmetric character, CD cannot be directly applied to the matrix representation of the spatial part of the two-electron spin-orbit coupling (2e-SOC) integrals. Here, we present a computational strategy to achieve a Cholesky representation of the spatial part of the 2e-SOC integrals, and propose a new efficient CD algorithm for both ERI and 2e-SOC integrals. The proposed algorithm differs from previous CD implementations by the extensive use of a full-pivoting design, which allows a univocal definition of the Cholesky basis, once the CD δ threshold is made explicit. We show that 2δ is the upper limit for the errors affecting the reconstructed 2e-SOC integrals. The proposed strategy was implemented in the ab initio program Computational Emulator of Rare Earth Systems (CERES), and tested for computational performance on both the ERI and 2e-SOC integrals evaluation. © 2017 Wiley Periodicals, Inc.

Communication: Almost error-free resolution-of-the-identity correlation methods by null space removal of the particle-hole interactions

We present a method to improve upon the resolution-of-the-identity (RI) for correlation methods. While RI is known to allow for drastic speedups, it relies on a cancellation of errors. Our method eliminates the errors introduced by RI which are known to be problematic for absolute energies. In this way, independence of the error compensation assumption for relative energies is also achieved. The proposed method is based on the idea of starting with an oversized RI basis and projecting out all of its unphysical parts. The approach can be easily implemented into existing RI codes and results in an overhead of about 30%, while effectively removing the RI error. In passing, this process alleviates the problem that for many frequently employed basis sets no optimized RI basis sets have been constructed. In this paper, the theory is presented and results are discussed exemplarily for the random phase approximation and Møller-Plesset perturbation theory.

An integral-factorized implementation of the driven similarity renormalization group second-order multireference perturbation theory

We report an efficient implementation of a second-order multireference perturbation theory based on the driven similarity renormalization group (DSRG-MRPT2) [C. Li and F. A. Evangelista, J. Chem. Theory Comput. 11, 2097 (2015)]. Our implementation employs factorized two-electron integrals to avoid storage of large four-index intermediates. It also exploits the block structure of the reference density matrices to reduce the computational cost to that of second-order Møller-Plesset perturbation theory. Our new DSRG-MRPT2 implementation is benchmarked on ten naphthyne isomers using basis sets up to quintuple-ζ quality. We find that the singlet-triplet splittings (ΔST) of the naphthyne isomers strongly depend on the equilibrium structures. For a consistent set of geometries, the ΔST values predicted by the DSRG-MRPT2 are in good agreements with those computed by the reduced multireference coupled cluster theory with singles, doubles, and perturbative triples.

Multireference Perturbation Theory with Cholesky Decomposition for the Density Matrix Renormalization Group

We present a second-order N-electron valence state perturbation theory (NEVPT2) based on a density matrix renormalization group (DMRG) reference wave function that exploits a Cholesky decomposition of the two-electron repulsion integrals (CD-DMRG-NEVPT2). With a parameter-free multireference perturbation theory approach at hand, the latter allows us to efficiently describe static and dynamic correlation in large molecular systems. We demonstrate the applicability of CD-DMRG-NEVPT2 for spin-state energetics of spin-crossover complexes involving calculations with more than 1000 atomic basis functions. We first assess, in a study of a heme model, the accuracy of the strongly and partially contracted variant of CD-DMRG-NEVPT2 before embarking on resolving a controversy about the spin ground state of a cobalt tropocoronand complex.

Ab Initio Investigation of the Ground States of F2P(S)N, F2PNS, and F2PSN

A recent spectroscopic experiment identified difluorothiophosphoryl nitrene (F₂P(S)N) and found that it showed rich photochemistry. However, a discrepancy between the experimental results and the quantum chemical calculations was reported. Thus, high-level ab initio calculations using the coupled cluster singles and doubles with perturbative triples and second-order multiconfigurational perturbation theory were performed to elucidate this inconsistency. The discrepancy arose due to the failure to consider the triplet state of difluoro(thionitroso)phosphine (F₂PNS). In this work, we identify that the global minimum of the system is the triplet state of F₂PNS, which allows us to explain the inconsistency between the experimental and theoretical results. All calculated results give consistent results with the recent experimental results.

The multilevel CC3 coupled cluster model

We present an efficient implementation of the closed shell multilevel coupled cluster method where coupled cluster singles and doubles (CCSD) is used for the inactive orbital space and CCSD with perturbative triples (CC3) is employed for the smaller active orbital space. Using Cholesky orbitals, the active space can be spatially localized and the computational cost is greatly reduced compared to full CC3 while retaining the accuracy of CC3 excitation energies. For the small organic molecules considered we achieve up to two orders of magnitude reduction in the computational requirements.