Orbital optimization in the perfect pairing hierarchy. Applications to full-valence calculations on linear polyacenes

TL;DR

This paper implements orbital optimization within the perfect pairing hierarchy to improve full-valence calculations on linear polyacenes, demonstrating that larger active spaces reduce apparent strong correlation effects and achieving high accuracy with novel computational methods.

Contribution

It introduces an efficient implementation of orbital optimization for the perfect pairing hierarchy, enabling large-scale full-valence calculations on polyacenes with improved correlation energy accuracy.

Findings

Orbital optimization avoids local minima and symmetry breaking issues.

Full-valence PH calculations capture over 95% of correlation energy.

Larger basis sets and active spaces weaken observed strong correlations.

Abstract

We describe the implementation of orbital optimization for the models in the perfect pairing hierarchy [Lehtola et al, J. Chem. Phys. 145, 134110 (2016)]. Orbital optimization, which is generally necessary to obtain reliable results, is pursued at perfect pairing (PP) and perfect quadruples (PQ) levels of theory for applications on linear polyacenes, which are believed to exhibit strong correlation in the {\pi} space. While local minima and {\sigma}-{\pi} symmetry breaking solutions were found for PP orbitals, no such problems were encountered for PQ orbitals. The PQ orbitals are used for single-point calculations at PP, PQ and perfect hextuples (PH) levels of theory, both only in the {\pi} subspace, as well as in the full {\sigma}{\pi} valence space. It is numerically demonstrated that the inclusion of single excitations is necessary also when optimized orbitals are used. PH is found to yield good agreement with previously published density matrix renormalization group (DMRG) data in the {\pi} space, capturing over 95% of the correlation energy. Full-valence calculations made possible by our novel, efficient code reveal that strong correlations are weaker when larger bases or active spaces are employed than in previous calculations. The largest full-valence PH calculations presented correspond to a (192e,192o) problem.

Tables

Table 1. ?tablename? 1: π 𝜋 \pi -space STO-3G correlation energy captured by the models using PQ orbitals, without and with single excitations, compared to DMRG reference data. The Hartree–Fock energy E HF superscript 𝐸 HF E^{\text{HF}} is also given.

			Without singles			With singles
Molecule	Active space	$E^{\text{HF}}$	PP	PQ	PH	PP	PQ	PH	DMRGa
2acene	(10e,10o)	-378.683899	-0.115752	-0.158155	-0.171489	-0.115752	-0.161151	-0.177136	-0.178294
3acene	(14e,14o)	-529.467421	-0.159805	-0.218652	-0.239966	-0.159805	-0.223151	-0.250444	-0.254307
4acene	(18e,18o)	-680.246082	-0.210661	-0.284669	-0.312433	-0.210661	-0.290035	-0.325936	-0.332661
5acene	(22e,22o)	-831.022093	-0.256128	-0.347849	-0.383498	-0.256128	-0.354469	-0.402041	-0.412627
6acene	(26e,26o)	-981.794501	-0.300955	-0.412119	-0.456519	-0.300955	-0.420058	-0.480374	-0.495683
8acene	(34e,34o)	-1283.332597	-0.374849	-0.531322	-0.599087	-0.374849	-0.542803	-0.639394	-0.668526
10acene	(42e,42o)	-1584.874961	-0.430551	-0.668290	-0.752407	-0.430552	-0.681619	-0.797103	-0.838580
12acene	(50e,50o)	-1886.418502	-0.488271	-0.806117	-0.909540	-0.488271	-0.822062	-0.958267	-1.007378b

Table 2. ?tablename? 2: Correlation energies in the STO-3G full valence space.

			Without singles			With singles
Molecule	Active space	Hilbert space	PP	PQ	PH	PP	PQ	PH
2acene	(48e,48o)	$1.0\times {10}^{27}$	-0.322737	-0.577926	-0.681109	-0.322737	-0.582018	-0.688946
3acene	(66e,66o)	$5.2\times {10}^{37}$	-0.440449	-0.795324	-0.945119	-0.440449	-0.801436	-0.958952
4acene	(84e,84o)	$2.8\times {10}^{48}$	-0.564961	-1.019499	-1.213242	-0.564962	-1.026913	-1.231753
5acene	(102e,102o)	$1.6\times {10}^{59}$	-0.684079	-1.240375	-1.479518	-0.684080	-1.249531	-1.504362
6acene	(120e,120o)	$9.3\times {10}^{69}$	-0.802603	-1.462114	-1.747772	-0.802604	-1.473113	-1.779061
8acene	(156e,156o)	$3.4\times {10}^{91}$	-1.023950	-1.894352	-2.279430	-1.023952	-1.910179	-2.327852
10acene	(192e,192o)	$1.3\times {10}^{113}$	-1.227908	-2.323432	-2.819420	-1.227910	-2.341611	-2.872236
12acene	(228e,228o)	$5.2\times {10}^{134}$	-1.433817	-2.750959		-1.433820	-2.772657

Table 3. ?tablename? 3: Full-valence STO-3G HONO-LUNO gaps for PP, PQ, PH.

	Without singles			With singles
2acene	1.8642	1.7288	1.7091	1.8642	1.7266	1.6998
3acene	1.8372	1.7011	1.6538	1.8372	1.7046	1.6469
4acene	1.8340	1.6852	1.6227	1.8340	1.6872	1.6045
5acene	1.8280	1.6658	1.5775	1.8280	1.6710	1.5616
6acene	1.8280	1.6535	1.5474	1.8280	1.6585	1.5211
8acene	1.8286	1.6242	1.4624	1.8286	1.6326	1.4271
10acene	1.8286	1.5095	1.3842	1.8286	1.5311	1.3216
12acene	1.8138	1.3446		1.8138	1.3905

Table 4. ?tablename? 4: π 𝜋 \pi -space HONO-LUNO gaps for PP, PQ, PH, compared with DMRG data from ref. 55 .

	Without singles			With singles
Molecule	PP	PQ	PH	PP	PQ	PH	DMRGa
2acene	1.8642	1.6990	1.6600	1.8642	1.6956	1.6433	1.6315
3acene	1.8372	1.6617	1.5770	1.8372	1.6645	1.5602	1.5303
4acene	1.8340	1.6411	1.5293	1.8340	1.6418	1.4913	1.4213
5acene	1.8280	1.6150	1.4570	1.8280	1.6190	1.4154	1.3063
6acene	1.8280	1.5992	1.4066	1.8280	1.6025	1.3409	1.1567
8acene	1.8286	1.5608	1.2559	1.8286	1.5671	1.1450	0.7885
10acene	1.8286	1.4507	1.1102	1.8286	1.4692	0.9276	0.5777
12acene	1.8138	1.2641	0.8942	1.8138	1.3063	0.7175	0.4717b

Table 5. ?tablename? 5: Correlation energies in the full cc-pVDZ valence active space. The Hartree–Fock energy E HF superscript 𝐸 HF E^{\text{HF}} is also given .

			Without singles			With singles
Molecule	Active space	$E^{\text{HF}}$	PP	PQ	PH	PP	PQ	PH
2acene	(48e,48o)	-383.383836	-0.301755	-0.533474	-0.612158	-0.301756	-0.536128	-0.617254
3acene	(66e,66o)	-536.037636	-0.412436	-0.734416	-0.848505	-0.412438	-0.738267	-0.857082
4acene	(84e,84o)	-688.687483	-0.527842	-0.940703	-1.088061	-0.527846	-0.945346	-1.099478
5acene	(102e,102o)	-841.335342	-0.639575	-1.144255	-1.325955	-0.639581	-1.149943	-1.340944

Table 6. ?tablename? 6: Full-valence cc-pVDZ HONO-LUNO gaps for PP, PQ, PH.

	Without singles			With singles
Molecule	PP	PQ	PH	PP	PQ	PH
2acene	1.9120	1.8348	1.8185	1.9120	1.8340	1.8149
3acene	1.8958	1.8198	1.7907	1.8958	1.8212	1.7871
4acene	1.8939	1.8094	1.7737	1.8939	1.8102	1.7665
5acene	1.8894	1.7984	1.7525	1.8894	1.8004	1.7454

Equations

$$⟨{\Phi }_0⟩\overline{H}{\Phi }_0=$$

$$⟨{\Phi }_i^a⟩\overline{H}{\Phi }_0=$$

$$⟨{\Phi }_{ij}^{ab}⟩\overline{H}{\Phi }_0=$$

$$⋮$$

$$E^{\text{CC}}=$$

$$\mathcal{L}=$$

$$\mathbf{C}\to$$

$$\mathbf{U}=$$

$$\mathbf{\theta }=$$

$$\frac{\partial E}{\partial {\Delta }_{tu}}=$$

$$=$$

$$\frac{{\partial }^{2}E}{\partial {\Delta }_{tu}^2}=$$

$$=$$

$$h_{pq}=$$

$$v_{pqrs}=$$

$$\mathbf{U}=$$

$${\left(\frac{\partial h_{pq}}{\partial {\theta }_{tu}}\right)}_{\mathbf{\theta }=0}=$$

$${\left(\frac{\partial v_{pqrs}}{\partial {\theta }_{tu}}\right)}_{\mathbf{\theta }=0}=$$

$$f_{pq}=$$

$${\left(\frac{\partial f_{pq}}{\partial {\theta }_{tu}}\right)}_{\mathbf{\theta }=0}={\left(\frac{\partial f_{pq}^0}{\partial {\theta }_{tu}}\right)}_{\mathbf{\theta }=0}+$$

$${\left(\frac{{\partial }^2{h}_{pq}}{\partial {\theta }_{tu}^2}\right)}_{\mathbf{\theta }=0}=$$

$${\left(\frac{{\partial }^2{h}_{pq}}{\partial {\theta }_{tu}\partial {\theta }_{ut}}\right)}_{\mathbf{\theta }=0}=$$

$$+$$

$${\left(\frac{{\partial }^2{f}_{pq}}{\partial {\theta }_{tu}^2}\right)}_{\mathbf{\theta }=0}={\left(\frac{{\partial }^2{f}_{pq}^0}{\partial {\theta }_{tu}^2}\right)}_{\mathbf{\theta }=0}+$$

$${\left(\frac{{\partial }^2{f}_{pq}}{\partial {\theta }_{tu}\partial {\theta }_{ut}}\right)}_{\mathbf{\theta }=0}={\left(\frac{{\partial }^2{f}_{pq}^0}{\partial {\theta }_{tu}\partial {\theta }_{ut}}\right)}_{\mathbf{\theta }=0}+$$

$$+$$

$${\left(\frac{{\partial }^2{v}_{pqrs}}{\partial {\theta }_{tu}^2}\right)}_{\mathbf{\theta }=0}=$$

$$+$$

$${\left(\frac{{\partial }^2{v}_{pqrs}}{\partial {\theta }_{tu}\partial {\theta }_{ut}}\right)}_{\mathbf{\theta }=0}=$$

$$+$$

$$E_0=$$

$$=$$

$$\frac{\partial E_0}{\partial {\Delta }_{ai}}=$$

Full text

\abbreviations

AO: atomic orbital, BFGS: Broyden–Fletcher–Goldfarb–Shanno, CAS-SCF: complete active space self-consistent field, CCVB: coupled cluster valence bond, CI: configuration interaction, DMRG: density matrix renormalization group, FCI: full configuration interaction, GDM: geometric direct minimization, HONO: highest occupied natural orbital, LUNO: lowest unoccupied natural orbital, MC-SCF: multiconfigurational self-consistent field, MO: molecular orbital, NO: natural orbital, NOON: natural orbital occupation number, OO-CC: orbital-optimized coupled cluster, OPDM: one particle density matrix, PP: perfect pairing, PQ: perfect quadruples, PH: perfect hextuples, PPH: perfect pairing hierarchy, TPDM: two particle density matrix, VDRM: variational reduced density matrix, VOO-CC: valence orbital-optimized coupled cluster

1 Introduction

Systems with strong correlation, in which many electronic configurations contribute significantly to the wave function, are important in many areas of chemistry such as catalysis. Unfortunately, their accurate yet efficient modeling is still an unsolved question in theoretical chemistry. The exact wave function for any system is available in theory by diagonalizing the molecular Hamiltonian in the basis of electron configurations, yielding the full configuration interaction (FCI) approach. But, FCI exhibits exponential scaling that limits its use to tiny systems [1, 2, 3, 4, 5, 6].

Restricting the number of electronic configurations allowed in the wave function yields the multiconfigurational self consistent field (MC-SCF) method [7]. However, as typically a considerable number – thousands to millions – of configurations are necessary to obtain a sufficiently converged wave function, the electronic configurations allowed in the wave function are more often defined in terms of active orbitals, as in the complete active space self-consistent field (CAS-SCF) method [8, 9, 10] or variants thereof [11, 12, 13]. CAS-SCF still maintains the exponential scaling of FCI with respect to the size of the active space, and is generally acknowledged to be limited to problems around the size of 16 electrons in 16 orbitals, denoted as (16e,16o). While ideally all valence electrons should be included in the active space, due to the limitation on the size of the feasible active space the choice of the active orbitals is generally a complicated problem, and improper choices often yield unreliable results [14]. The necessary active space can also change e.g. along a reaction path, further underlining the problems caused by the limitations of this type of an approach.

Because most of the configurations in an FCI (or CAS-SCF) wave function typically have negligible weights, it is possible to truncate the wave function without losing significant accuracy [15, 16, 17, 18]. This has been recently been exploited in various stochastic and adaptive approaches to FCI and CAS-SCF [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31], which have made quasi-FCI calculations feasible on much larger systems than before. For instance, ref. 31 describes the solution of the (28e,22o) problem in oxoMn-salen to sub-millihartree accuracy in less than a minute on 20 cores. Although these novel approaches for the solution of the FCI problem have made calculations possible on much larger systems than with conventional algorithms, due to their reliance on the FCI ansatz, they still scale exponentially in the size of the active space.

The density matrix renormalization group (DMRG) approach [32, 33, 34, 35, 36] is able to reproduce the FCI result, yet it affords polynomial scaling in linear systems. Unfortunately, for non-linear systems, like CAS-SCF or FCI, DMRG succumbs to exponential scaling as well [36]. Still, DMRG has successfully been applied to considerably larger active spaces than CAS-SCF [37, 38, 39], and is generally thought to be the best general-purpose tool for strong correlation as it can be used to treat extremely challenging problems such as open-shell transition metal compounds. In practice, DMRG is limited to ($n_{\text{el}}$e,$n_{\text{orb}}$o) problems with $n_{\text{el}}n_{\text{orb}}\leq 2000$ [36]. Furthermore, DMRG is not trivial to use due to issues with the choice of various convergence parameters. There is thus still demand for approximate methods that are able to capture the bulk of strong correlation effects while also being fast enough to be able to handle a full-valence active space. The perfect pairing hierarchy (PPH) is one possible candidate for such a method.

The methods in the PPH are obtained by truncating the equations of infinite-order coupled cluster (CC) theory [40] to be exact for a singlet state (which can be open- or closed-shell) in a ($n$e,$n$o) active space for a given number of electrons $n$ [41, 42]. The truncation of PPH bears some similarity to, but is not to be confused with the $n$CC hierarchy proposed by Bartlett and Musiał [43], where the CC diagrams are truncated for exactness for $n$ electrons. Thus, the PPH can alternatively be understood as a further truncation of $n$CC to an active space.

The simplest model in the PPH, perfect pairing [44, 45, 46, 47, 48, 49] (PP), is exact for a $n=2$ singlet, whereas the perfect quadruples [41] (PQ) and perfect hextuples [42] (PH) models are exact for $n=4$ and $n=6$, respectively. The PP, PQ, and PH models form a hierarchy of approximations to CAS-SCF. With our recent, efficient implementation of the PQ and PH models [50], this hierarchy can be used to perform approximate CAS-SCF calculations on systems of unforeseen size, as will be seen later on in the manuscript.

While CC theory is invariant to orbital rotations in the occupied-occupied and virtual-virtual blocks, the results of truncated CC approaches generally depend on the set of occupied orbitals used in the calculation. In orbital-optimized CC (OO-CC) theory [51] occupied-virtual (ov) rotations are performed to minimize the CC energy. Furthermore, as core orbitals are chemically inactive, an active space can be introduced, as is also done in CAS-SCF, yielding the valence orbital-optimized CC (VOO-CC) method [52]. The introduction of an active space necessitates optimization of rotations between active and inactive occupied orbitals, as well as active and inactive virtual orbitals, in addition to the ov rotations in OO-CC. In the PPH, mixings within the active occupied, as well as active virtual orbital spaces must further be considered, because the truncation is based on the orbitals.

In the present manuscript, we describe how orbital optimization can be efficiently implemented for the models within the PPH, and present their applications to full-valence calculations on linear polyacenes. Ever since preliminary calculations predicted linear polyacenes to have singlet diradical ground states [53], both linear and cyclacenes have been studied intensively using a variety of methods [54]. Famously, DMRG calculations have been used to show that the larger linear polyacenes exhibit strong correlation for their $\pi$ electrons [55]. Other approaches used have included density functional theory [56, 57, 58, 59, 60, 61], multiconfiguration pair-density functional theory [62], projected Hartree–Fock theory [63], the random phase approximation [64], configuration interaction [60] (CI), adaptive CI [27], GW theory [65], variational two-electron reduced density matrix (VRDM) theory [66, 67, 68, 69], Møller–Plesset perturbation theory [70, 71, 72, 73, 61], spin-flip methods [74, 75], CAS-SCF [53, 71, 72, 59, 76] as well as restricted-active space self-consistent field theory [77], valence bond [78] and CC valence bond (CCVB) theory [79, 80], CC theory [70, 71, 72, 80] and multireference averaged quadratic CC theory [81, 82], an algebraic diagrammatic construction scheme [83], as well as PPH methods with approximate orbitals [50].

We chose these systems for our study, because the correlation exhibited is too strong for successful description with conventional CC approaches, such as CC with full single and double excitations (CCSD) [71, 70, 80], and because the CAS-SCF, VRDM, and DMRG studies have been restricted to the $\pi$ space only, as full-valence calculations would otherwise be too costly for these methods for anything beyond the smallest acenes. As far as we are aware, the only full-valence calculations on large polyacenes employing a method potentially capable of accurate treatment of strong correlation that have been published are the CCVB calculations performed contemporaneously in our group [79, 80]. Thanks to our recent implementation [50], even the larger polyacenes can be treated accurately in full-valence active spaces using PPH methods, yielding an independent way to check the correctness of the recent CCVB results.

As the result of our application, we find that orbital optimization can be routinely applied to large systems, and that it significantly improves the accuracy of the PPH models. We show that even when optimal orbitals are used, it is necessary to include single excitations as the optimal orbitals are not Brueckner orbitals. We show that PH successfully describes strong correlation in the STO-3G $\pi$ space compared to literature DMRG reference values [55]. In agreement with the recent CCVB results of ref. 80, and in contrast to the $\pi$ active space results obtained using e.g. DMRG [55] and VDRM [66, 68, 69], we show the strong correlation is quenched when all valence orbitals are included in the active space. Furthermore, going beyond the STO-3G basis for the first time with a high-level full-valence method shows that improvement of the basis set leads to even more quenching of the strong correlation for the same active space.

In the following, $i$, $j$ and $k$ denote occupied spin-orbitals, $a$ and $b$ denote unoccupied spin-orbitals, and $p$, $q$, $r$, $s$, $t$ and $u$ denote general spin-orbitals.

2 Theory

Analytic gradients and orbital-optimized CC theory have been discussed extensively elsewhere [51, 52, 84], but a brief overview of the optimization procedure is given below for the sake of completeness. The CC equations are given by

[TABLE]

where $\overline{H}$ is the effective Hamiltonian, the correlation energy $E_{c}$ is given by (LABEL:CC-E), and (LABEL:CC-S,CC-D) determine the single and double excitation ampitudes, respectively. Regardless of the truncation of the amplitudes, the CC energy is determined solely by the single and double excitations as

[TABLE]

where $E_{0}$ is the Hartree–Fock reference energy, $f_{i}^{a}$ is a Fock matrix element and $\Braket{ij}{ab}$ is a two-electron integral, and $t_{i}^{a}$ and $t_{ij}^{ab}$ are the single and double excitation amplitudes, respectively.

The orbitals enter the problem through the matrix elements $f_{i}^{a}$ and $\Braket{ij}{ab}$. However, the orbital rotation gradient cannot be calculated from (LABEL:E-CC), because the amplitudes $t$ depend on the matrix elements in the first order. To be able to optimize the orbitals, one must demand the simultaneous satisfaction of the CC equations through Lagrangian multipliers $\lambda$ [85]

[TABLE]

that are known as de-excitation amplitudes, in analogy to the excitation amplitudes $t$.

Demanding that the variation of the CC Lagrangian ((LABEL:L)) vanish both with respect to the $t$ and the $\lambda$ amplitudes, equations are obtained for the $\lambda$ and the $t$ amplitudes, respectively. Having solved for the values of $t$ and $\lambda$, the terms in (LABEL:L) can be regrouped into the standard form

[TABLE]

where $\gamma_{pq}$ and $\Gamma_{pqrs}$ are then identified as the one-particle (OPDM) and two-particle density matrices (TPDM), respectively. Unlike the CC energy ((LABEL:E-CC)), the CC Lagrangian ((LABEL:L,_E-DM)) is correct to the first order in the matrix elements. Thus, $\partial\mathcal{L}/\partial\theta_{pq}$ is a proper gradient for a rotation of the orbitals through an angle $\boldsymbol{\theta}$

[TABLE]

The in-detail equations for the gradient are given in the Appendix.

3 Computational Details

The geometric direct minimization (GDM) approach [86], which uses a limited-memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) approximate Hessian [87], is used to optimize the orbitals [88, 49]. No parallel transport of previous gradients in the BFGS is used in our implementation, as we are viewing the optimization problem on the Stiefel manifold where the efficient implementation of parallel transport is an unsolved problem [89]. To speed up the convergence of the optimization, also the diagonal Hessian $\partial^{2}\mathcal{L}/\partial\theta_{pq}^{2}$ is calculated (see the Appendix) and used to initialize the BFGS solver.

Once a search direction $\boldsymbol{\theta}$ has been chosen with GDM, a line search is performed to minimize the Lagrangian ((LABEL:E-DM)); that is, $\mathcal{L}(\exp(\ell\boldsymbol{\theta}))$ is minimized with respect to the step size $\ell$. Because $\mathcal{L}$ is a quasiperiodic function along $\ell$ [90, 91], the period of its fastest oscillation $\lambda_{0}$ can be estimated from the largest absolute eigenvalue $|\omega|_{\text{max}}$ of $\boldsymbol{\theta}$ as $\lambda_{0}=\pi/2|\omega|_{\text{max}}$. This expression can be derived when the Lagrangian has fourth-order dependence on the orbitals [92], as is the case here due to the two-electron integrals. Relying on earlier experience [92], we chose to perform the line search using parabolic interpolation with a trial step size of $\lambda_{0}/5$. Once a optimal value $\ell$ has been found, the orbitals are updated $\boldsymbol{C}\to\boldsymbol{C}\exp(\ell\boldsymbol{\theta})$ and the new reference energy $E_{0}$ and density matrices $\gamma$ and $\Gamma$ are evaluated. Then, the new gradient and diagonal Hessian are computed, followed by either a new line search, or the determination that successful convergence has been achieved.

The line search minimization of (LABEL:E-DM) has two obvious approaches. First, one may allow the density matrices $\gamma$ and $\Gamma$ to relax at every point along the line, which appears the obvious solution. However, one could also consider (LABEL:E-DM) as a model energy functional by assuming the density matrices $\gamma$ and $\Gamma$ to be fixed. After some experimentation, we chose to pursue the latter approach, as in addition to the slight benefit of eliminating the need to solve new amplitudes during the line search, it also makes the search direction and preconditioning in GDM exact. Note that even when the density matrices are frozen within the line search, every optimization step must result in a decrease of the true CC energy, as the true, relaxed energy will always lie below the energy predicted by the model functional.

The gradient and diagonal Hessian equations are autogenerated in a similar fashion to the CC equations, following the machinery detailed in ref. 50. The density matrices and integrals are represented as dense subtensors in the molecular orbital (MO) basis, which is sparse in the PPH. Also the gradient and the diagonal Hessian are generated in the paired MO basis, as this enabled re-use of the code developed in ref. 50. However, the Fock matrix response term in the gradient (the two-electron integral terms in (LABEL:f-grad) in the Appendix) is excluded from the automatic evaluation, as it would result in the calculation of 3-index molecular integrals for PP and 4-index molecular orbitals for PQ, which would dominate the computational cost of the models. Instead, the response term is evaluated with hand-written code in the atomic orbital (AO) basis, allowing both for better scaling and for taking advantage of AO integral sparsity, reducing the cost of evaluating the term to the same level as Hartree–Fock theory.

The methods have been implemented in a general-purpose library, which may be made publicly available in the future. The input for a calculation is composed of the matrix elements, that is the MO basis Fock and B matrices, which are read in from the disk. The B matrices can be precomputed using either Cholesky [93] or resolution-of-the-identity [94] techniques. After the amplitudes have been converged, the computer generated implementation outputs the energy, the OPDM and the TPDM, as well as the orbital gradient and the diagonal Hessian as matrices if their calculation has been requested. This minimal interface guarantees straightforward interfacing possibilities to various programs.

The results of the present manuscript have been obtained by interfacing the library with the Erkale program [95, 96]. Hartree–Fock occupied orbitals, localized [91] with the generalized Pipek–Mezey method [97] with Becke charges and then paired with corresponding virtual orbitals using the Sano procedure [98], have been used to initialize the orbital optimization by GDM. Cholesky decompositions [93, 99] are used for all integral computations, with a $10^{-10}$ screening threshold for the two-electron integral calculations, and a $10^{-6}$ threshold for the Cholesky procedure. The AO Cholesky decomposition is run only once for a given molecule, at the beginning of the Hartree–Fock calculation, constituting a trivial fraction of the total runtime of the models.

Because single excitations are analogous to orbital rotations [100] and thus imitate orbital relaxation effects within the active space, the orbital optimization is run without single excitations as is done in OO-CC and VOO-CC. A gradient convergence threshold of $\left\|\partial\mathcal{L}/\partial\theta_{pq}\right\|\leq 10^{-5}$ is used for the PPH orbital optimization. Having converged the orbitals with PP or PQ, single-point calculations are run with PP, PQ, and PH, and density matrices are calculated. Natural orbitals (NOs) and NO occupation numbers (NOONs) are obtained by diagonalizing the OPDM.

4 Results

While the automated procedures would allow us to formulate an orbital optimization procedure for all the models in the PPH, in the present work orbital optimization is undertaken only with PP and PQ, which have been shown to exhibit excellent scaling with system size [50]. Orbital optimization with PH is not considered, because the truncation of the intermediates that is performed in PH requires consistent definitions of all the intermediates in the $t$, $\lambda$ as well as the density matrix equations, but this is not guaranteed in our code at present. We estimate the effect of the inconsistent truncation on the NOONs to be negligible (changes in the occupation numbers of the order of 0.001).

To facilitate comparison with DMRG reference values, the molecular geometries for the polyacenes were adopted from ref. 55. With the molecule placed on the $xy$ plane, $\sigma$ and $\pi$ orbital symmetries were assigned by determining the parity of the orbital $\phi$ upon reflection about the plane $\phi(x,y,-z)=\pm\phi(x,y,z)$. Like the original Pipek–Mezey method [101], the generalized Pipek–Mezey method [97] used in the present work properly separates $\sigma$ and $\pi$ orbitals, but lacks severe mathematical shortcomings of the original formulation by Pipek and Mezey that is based on Mulliken charges [102].

Orbital optimizations restricted to the $\pi$ and $\sigma$ active spaces, respectively, were attempted in the STO-3G basis set [103] by deleting the virtual orbitals of the opposite class. Because STO-3G is a minimal basis set, there are no inactive virtual orbitals and the optimization problem only needs to be solved within the active space. For both the $\sigma$- and the $\pi$-space, the PP optimizations converged within a few dozen iterations. For PQ, the $\sigma$-space orbital optimizations converged in a few hundred iterations. However, the convergence of the PQ $\pi$-space calculations was poor: while for the smaller acenes convergence was achieved in several hundred iterations, for 6acene the orbital optimization was still not converged at 500 iterations, and for even larger acenes the $\lambda$ amplitude equations diverged during the orbital optimization. (These problems can be understood in light of the results discussed later on in the manuscript: the restriction to the $\pi$ space exaggerates correlation effects, which then make orbital optimization hard.) For this reason, we will not consider orbital optimization restricted to the $\sigma$ or $\pi$ space further in the manuscript.

For the full-valence STO-3G orbital optimization with PP, solutions that break $\sigma$-$\pi$ symmetry by having some orbitals of mixed character were found for multiple molecules. In contrast, we did not encounter any broken symmetry solutions with full-valence STO-3G PQ optimized orbitals. We have thoroughly tested this in the case of 2acene. In analogy to our recent study in self-interaction corrected density functional theory [92], 100 different initial guesses were considered for the initial occupied orbitals, after which the orbitals were localized with the generalized Pipek–Mezey method and corresponding virtuals determined with the Sano algorithm. Orbital optimization initialized with these guesses was performed using both PP and PQ. For PP, a quarter of the calculations converged on a symmetry broken solution, while three quarters found the $\sigma$-$\pi$ preserving solution which is higher in energy. For PQ, all calculations converged to the same solution that preserves $\sigma$-$\pi$ symmetry. Calculations were also performed without the initial generalized Pipek–Mezey / Sano localization procedure, in which case all PP calculations converged onto the symmetry broken solution, and PQ again converged to the $\sigma$-$\pi$ preserving solution.

It is not hard to understand why PP would have more local solutions than PQ. In PP orbital optimization (i.e. no single excitations), the only excitations are intra-pair ones: the double excitation $t_{i\bar{i}}^{i^{\star}\bar{i}^{\star}}$ from the occupied alpha and beta orbital $i$ and $\bar{i}$, respectively, to the corresponding paired virtuals $i^{\star}$ and $\bar{i}^{\star}$, respectively. Furthermore, each amplitude $t_{i\bar{i}}^{i^{\star}\bar{i}^{\star}}$ can be solved analytically, independent of the other pairs [49]. In contrast, in PQ orbital optimization, several inter-pair double excitations $t_{j\bar{i}}^{i^{\star}\bar{i}^{\star}}$, $t_{i\bar{j}}^{i^{\star}\bar{i}^{\star}}$, $t_{i\bar{i}}^{j^{\star}\bar{i}^{\star}}$, $t_{i\bar{i}}^{i^{\star}\bar{j}^{\star}}$, $t_{i\bar{i}}^{j^{\star}\bar{j}^{\star}}$, $t_{i\bar{j}}^{i^{\star}\bar{j}^{\star}}$, $t_{j\bar{i}}^{j^{\star}\bar{i}^{\star}}$ are added (as well as same-spin double excitations), as well as the triple excitations $t_{ij\bar{j}}^{i^{\star}j^{\star}\bar{j}^{\star}}$, $t_{ij\bar{j}}^{j^{\small\star}i^{\star}\bar{i}^{\star}}$, $t_{ji\bar{i}}^{i^{\star}j^{\star}\bar{j}^{\star}}$, and $t_{ji\bar{i}}^{j^{\star}i^{\star}\bar{i}^{\small\star}}$ and quadruple excitations $t_{ij\bar{i}\bar{j}}^{i^{\star}j^{\star}\bar{i}^{\star}\bar{j}^{\star}}$. Thus, PP describes single Kekulé structures, whereas PQ adds in resonance terms which couple the structures, eliminating the local solutions. Because the PQ orbitals thus appear well defined, we chose to use the full-valence optimized PQ orbitals for all the single-point calculations in our study.

The scaling of the STO-3G calculations on one of the authors’ (S.L.) workstation, using a single Intel i7-4770 processor core, is shown in 1. Already in the STO-3G basis, where there are no inactive virtual orbitals, the integral calculations and the build of the orbital gradient of the Lagrangian are close to being the rate determining steps for the orbital optimization. Because everything else is built only within the active space and thus does not scale with the size of the basis, the computation of the integrals and the orbital gradient contraction will become even more dominant in larger bases: PP has $O(N_{\text{act}})$ matrix elements and $O(N_{\text{act}}N_{\text{bf}})$ orbital gradient matrix elements, while PQ has $O(N_{\text{act}}^{2})$ matrix elements and $O(N_{\text{act}}^{2}N_{\text{bf}})$ orbital gradient matrix elements, where $N_{\text{act}}$ and $N_{\text{bf}}$ are the number of active orbitals and basis functions, respectively, when the Fock response terms are handled otherwise as discussed above in the Computational Detail section. An optimized AO implementation has been reported for PP [49]; similar techniques might prove beneficial for PQ as well. While PQ scales asymptotically as $N_{\text{act}}^{4}$ [41, 50], PQ is still cubic scaling for the 12acene calculation with $N_{\text{act}}=228$, suggesting that enormous calculations will be possible with a fully optimized implementation.

4.1 STO-3G basis, $\pi$ space

To establish the accuracy of the orbital-optimized models, we begin in the STO-3G basis [103] where accurate DMRG data is available in the $\pi$ space [55]. Calculations are performed using the $\pi$ orbitals extracted from the full-valence optimized PQ orbitals. Because $\sigma$-$\pi$ correlation is expected to be much less important than $\pi$-$\pi$ correlation, the $\pi$ orbitals extracted from the full-valence PQ optimization should not be exceedingly different from the ones that would be optimal for a $\pi$-space-only PQ calculation, as all the $\pi$ orbitals of STO-3G are anyhow included in the full-valence $\pi$ active space.

The correlation energy captured by PP, PQ, and PH, along with the Hartree–Fock energy are shown in 1. Because PQ optimized orbitals are used for all models, results are given both with and without single excitations included in the single-point calculation, as these represent relaxations in the active occupied – active virtual block. (Note however, that the singles do not relax the other orbital degrees of freedom in PPH, namely, the active occupied – active occupied and active virtual – active virtual rotations.) One would expect this active space orbital relaxation effect to be minimal. However, it has also been argued that OO-CC theories that rely on projective approaches fail to reproduce FCI results [104], because the singles amplitudes at the optimized orbitals are not guaranteed to vanish, as they do for true Brueckner orbitals [105]. Thus, the results obtained including single excitations should be more accurate than the ones obtained without singles. But, due to the non-variational character of projected CC theory, the inclusion of the single excitations may also result in a significant underestimation of the total energy. The importance of the inclusion of single excitations is clear from the results in 1. The energies estimated with single excitations are noticeably lower than the ones estimated without single excitations, and are in better agreement with DMRG.

The agreement of PH with DMRG is excellent. Without singles, PH captures >90% of the correlation energy for all systems, and with singles, PH captures >95% of the correlation energy for all acenes, even though the $\pi$ electrons have a significant degree of delocality.

Comparing the numbers of 1 with the analogous data using approximate orbitals (corresponding to the initial guess used in the present work) in ref. 50 shows remarkable improvement in the PQ results due to the orbital optimization carried out in the present work. While PQ orbitals are suboptimal for PP, also the PP energies estimated with PQ orbitals are significantly lower than those estimated with approximate orbitals.

PQ orbitals are also suboptimal for PH. However, the comparison with the energies of ref. 50 suggests that PQ orbitals yield lower energies for 2acene–6acene, while for 8acene–12acene the energies estimated from the approximate orbitals are below the present estimates obtained with the PQ orbitals. Because PH orbital optimization has not been attempted in the present work due to technical issues, it remains to be seen whether this phenomenon is caused by non-variational effects in PH, or the more delocalized nature of PH optimal orbitals compared to PQ optimal orbitals as speculated in ref. 50.

The $\pi$-space NOONs reproduced by PP, PQ, and PH on PQ orbitals are shown in 2. Like the earlier results with approximate orbitals [50], the NOONs reproduced by PH with PQ orbitals are in qualitative agreement with the DMRG results of ref. 55. However, important differences are noticeable. Compared to the earlier approximate orbital results [50], the signature of strong correlation is apparent with optimized orbitals already in PQ. Also, the overcorrelation of PH for 12acene has been removed, implying it was caused by nonvariational behavior originating from a poor orbital guess – this may be like the difference between nonvariational CCSD energies for bond-stretching with restricted vs. unrestricted orbitals. The remaining differences between the PH and DMRG results can be attributed to two factors: suboptimal orbitals (PQ instead of PH) and the remaining truncation error in PH.

[FIGURE:]

4.2 STO-3G basis, full-valence calculations

Having discussed the accuracy of the models in the $\pi$ space, we are in a good position to continue onto the full valence space. The excellent accuracy in the $\pi$ space is noteworthy, because as discussed in ref. 50, the perfect pairing hierarchy can be seen as a local correlation model. As the models are able to successfully treat the delocalized $\pi$ electrons in the polyacenes, they will be even better in treating the additional correlation due to the $\sigma$ electrons, which are localized.

The results of our full-valence STO-3G calculations are shown in 2 for the correlation energies, and 3 for the NOONs. As is evident from the plot, the strong correlation that is clear in the $\pi$ space calculation is quenched in the full-valence calculation. \Tabrefhonogap-fv-sto3g shows the gap between the highest occupied NO (HONO) and the lowest unoccupied NO (LUNO) in the full-valence calculation, which is significantly larger than the one in the $\pi$-space only calculation in 4.

The single excitations are again found to play a significant role. Contrary to the previous case dealing with the subset of $\pi$ orbitals, here the orbitals have been optimized for the same problem. Thus, if the orbitals had a quasi-Brueckner character, there would be little difference between the results obtained with and without single excitations. In contrast, the results in 2 clearly show that this is not the case: single excitations are necessary even when the orbitals have been optimized.

[FIGURE:]

4.3 cc-pVDZ basis

Full-valence orbital optimizations were performed up to 5acene in the cc-pVDZ basis set [106]. The resulting correlation energies, NOONs, and HONO-LUNO gaps are shown in 5, 4, and 6, respectively. The correlation energies are smaller in cc-pVDZ than in STO-3G, and the HONO-LUNO gaps are larger. While the strong correlation of the $\pi$ electrons was already quenched in the full valence space of STO-3G, the signature of strong correlation is even more strongly quenched in the cc-pVDZ basis, even though the size of the used active space is the same in STO-3G and cc-pVDZ. As a result of the smaller amount of correlation, the PP, PQ and PH results are closer to each other in cc-pVDZ than in STO-3G. It is thus clear that the limited flexibility of the STO-3G basis set results in an exaggeration of the correlation effects, and the question about the polyradicaloid character of the polyacenes may still be unresolved.

[FIGURE:]

5 Summary and Discussion

We have presented the efficient implementation of orbital optimization for the models in the perfect pairing hierarchy, and applied the new implementation to calculations on linear polyacenes. Our results show that orbital optimization is feasible for PQ for huge system sizes – up to (228e,228o) in the present work – and that the optimization significantly improves the quality of the results. The cost of orbital optimization in PP and PQ is dominated by the evaluation of matrix elements, the cost of which can be reduced by various optimization techniques, such as evaluation in the atomic orbital basis.

While it was found that PP may produce broken symmetry solutions that lie energetically below the $\sigma$-$\pi$ symmetry preserving solution, only a single PQ solution – which furthermore preserves $\sigma$-$\pi$ symmetry – was found for every molecule in the polyacene series. For this reason, full-valence optimized PQ orbitals were used for all the single-point calculations in the present work.

As has been argued in the literature [104], even when optimized orbitals are used, it is necessary to include single excitations in the treatment, because the Brueckner condition [105] may not be satisfied in projected CC theories which are non-variational. We have demonstrated that this is true also within the PPH: the inclusion of single excitations leads to significant energy lowerings in full-valence PQ calculations using optimized orbitals. In principle this energy lowering may also be due to non-variational effects, but we did not witness any breaking of variationality in the present work.

We have shown that PH used on top of PQ orbitals successfully captures the strong correlation in the $\pi$ space of the polyacenes in STO-3G, capturing over 95% of the correlation energy relative to DMRG. While PH is already successful in the delocalized $\pi$ space, because it can be understood as a three-pair local correlation model [50], it should perform even better for the additional description of the localized $\sigma$ electrons. We have shown that once the $\sigma$ electrons are included to complete the full valence active space, the strong correlation of the $\pi$ electrons of the polyacenes is quenched in the STO-3G calculations. A similar result has been obtained contemporaneously in our group through coupled cluster valence bond calculations [80]. The reason why the results of ref. 80 agree with ours can be understood by analysis of the results of the PH calculations. In the method of ref. 80, excitations are included to the approximate quadruples level. In our PH calculations, the role of quintuple and hextuple excitations (which CCVB lacks) is small, indicating they describe dynamic correlation effects that have negligible importance on the natural orbital occupation numbers.

We have further presented the first accurate ab initio calculations on polyacenes that go beyond the minimal STO-3G basis set. Our calculations in the cc-pVDZ basis set, employing the same active spaces as in STO-3G, show that the adoption of a non-minimal basis set further reduces the signature of strong correlation in the full valence space. Thus, the polyradicaloid character of the smaller acenes claimed in works by multiple authors may then be purely an effect of an insufficient active space and basis set, highlighting the need for further calculations, as the onset of polyradicaloid behavior in the acene series may be delayed when all valence electrons are correlated.

Orbital optimization for strong correlation problems is often non-trivial due to the existence of broad minima. For this reason, many approaches that guarantee second-order convergence have been developed for MC-SCF [107, 108, 10, 109, 110, 111, 112, 113, 114, 115, 116] and DMRG [117, 118, 119] calculations. While the Lagrangian formulation ((LABEL:E-DM)) is not correct to second order in the orbitals, it has been recently used for a quadratic convergence OO-CC implementation [84]. A similar approach should be feasible for the PPH as well – possibly supplemented by the use of the proxy function approach (freezing the density matrices between orbital update steps) which has been used in the present work. Although the cost per iteration as well as storage requirements in such an approach would be significantly higher than in the present one based on GDM, a second-order approach may yield converged orbitals in considerably fewer iterations and could be pursued in future work.

Funding

This work was supported by the Director, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division of the U.S. Department of Energy, under Contract No. DE-AC02-05CH11231.

Appendix

Here, we give the equations necessary to implement the rotation gradients and diagonal Hessians of the integrals. The rotation matrix $\boldsymbol{U}$ that is used to rotate the orbitals $\boldsymbol{C}$

[TABLE]

is parametrized as

[TABLE]

where $\boldsymbol{\theta}$ is an anti-hermitian matrix of rotation parameters. Due to the antihermicity, we can parametrize the rotations with either only the upper or only the lower triangular part of $\boldsymbol{\theta}$. Using the former choice we write

[TABLE]

or more symbolically

[TABLE]

The derivative with respect to the irreducible parameters $\boldsymbol{\Delta}$ is obtained as

[TABLE]

and the diagonal Hessian as

[TABLE]

The starting point to obtain the rotation gradients and diagonal Hessians is to write integrals rotated by an angle $\boldsymbol{\theta}$ as

[TABLE]

where the capital indices $P$, $Q$, $R$ and $S$ refer to the reference set of orbitals, $h_{pq}$ are one-electron integrals, and the shorthand notation $v_{pqrs}=\Braket{pq}{rs}$ is used for the two-electron integrals. Then, by substituting the Taylor expansion

[TABLE]

the gradient with respect to the rotation parameters $\boldsymbol{\theta}$ can be calculated, and evaluated at the reference orbitals ($\boldsymbol{\theta}=\boldsymbol{0}$) to give the orbital rotation gradient

[TABLE]

where $\delta_{pq}$ is the Kronecker delta symbol. Similarly, the gradient for the two-electron integral is

[TABLE]

The Fock matrix is given by

[TABLE]

where the sum over $k$ runs over occupied spin-orbitals, and thus the gradient from the Fock matrix elements also gathers a two-electron response term

[TABLE]

Here, $(\partial f_{pq}^{0}/\partial_{tu})_{\boldsymbol{\theta}=\boldsymbol{0}}$ denotes the gradient of the “frozen” Fock matrix, which is given by (LABEL:h-grad). Similar expressions to \eqrangerefh-gradf-grad have been presented in ref. 51. Using the same methodology, the diagonal Hessian components can be found out to be

[TABLE]

These equations are sufficient to calculate the CC contribution to the orbital gradient and diagonal Hessian in (LABEL:E-DM). The contributions from the reference energy

[TABLE]

are known to be[120]

[TABLE]

?refname?

• [1]

P.J. Knowles and N.C. Handy, J. Chem. Phys. 91 (4), 2396 (1989).

• [2]

P.J. Knowles and N.C. Handy, Comput. Phys. Commun. 54 (1), 75–83 (1989).

• [3]

E. Rossi, G.L. Bendazzoli, S. Evangelisti and D. Maynau, Chem. Phys. Lett. 310 (5-6), 530–536 (1999).

• [4]

L. Thøgersen, J. Olsen, D. Yeager, P. Jørgensen, P. Sałek and T. Helgaker, J. Chem. Phys. 121 (1), 16–27 (2004).

• [5]

Zhengting Gan and R. Harrison, in ACM/IEEE SC 2005 Conf. (, , 2005), pp. 22–22.

• [6]

Z. Gan, D.J. Grant, R.J. Harrison and D.A. Dixon, J. Chem. Phys. 125 (12), 124311 (2006).

• [7]

J. Hinze, J. Chem. Phys. 59 (12), 6424 (1973).

• [8]

P. Siegbahn, A. Heiberg, B. Roos and B. Levy, Phys. Scr. 21 (3-4), 323–327 (1980).

• [9]

B.O. Roos, Int. J. Quantum Chem. 18 (S14), 175–189 (1980).

• [10]

P.E.M. Siegbahn, J. Almlöf, A. Heiberg and B.O. Roos, J. Chem. Phys. 74 (4), 2384–2396 (1981).

• [11]

J. Olsen, B.O. Roos, P. Jørgensen and H.J.A. Jensen, J. Chem. Phys. 89 (4), 2185 (1988).

• [12]

T. Fleig, J. Olsen and C.M. Marian, J. Chem. Phys. 114 (11), 4775–4790 (2001).

• [13]

D. Ma, G. Li Manni and L. Gagliardi, J. Chem. Phys. 135 (4), 044128 (2011).

• [14]

V. Veryazov, P.Å. Malmqvist and B.O. Roos, Int. J. Quantum Chem. 111 (13), 3329–3338 (2011).

• [15]

J. Ivanic and K. Ruedenberg, Theor. Chem. Acc. 106 (5), 339–351 (2001).

• [16]

J. Ivanic and K. Ruedenberg, Theor. Chem. Accounts Theory, Comput. Model. (Theoretica Chim. Acta) 107 (4), 220–228 (2002).

• [17]

L. Bytautas, J. Ivanic and K. Ruedenberg, J. Chem. Phys. 119 (16), 8217 (2003).

• [18]

L. Bytautas and K. Ruedenberg, Chem. Phys. 356 (1-3), 64–75 (2009).

• [19]

G.H. Booth, A.J.W. Thom and A. Alavi, J. Chem. Phys. 131 (5), 054106 (2009).

• [20]

G.H. Booth and A. Alavi, J. Chem. Phys. 132 (17), 174104 (2010).

• [21]

G.H. Booth, D. Cleland, A.J.W. Thom and A. Alavi, J. Chem. Phys. 135 (8), 084104 (2011).

• [22]

J.J. Shepherd, G.H. Booth and A. Alavi, J. Chem. Phys. 136 (24), 244101 (2012).

• [23]

J.J. Shepherd, G. Booth, A. Grüneis and A. Alavi, Phys. Rev. B 85 (8), 081103 (2012).

• [24]

J.J. Shepherd, A. Grüneis, G.H. Booth, G. Kresse and A. Alavi, Phys. Rev. B 86 (3), 035111 (2012).

• [25]

C. Daday, S. Smart, G.H. Booth, A. Alavi and C. Filippi, J. Chem. Theory Comput. 8 (11), 4441–4451 (2012).

• [26]

R.E. Thomas, Q. Sun, A. Alavi and G.H. Booth, J. Chem. Theory Comput. 11 (11), 5316–5325 (2015).

• [27]

J.B. Schriber and F.A. Evangelista, J. Chem. Phys. 144 (16), 161106 (2016).

• [28]

N.M. Tubman, J. Lee, T.Y. Takeshita, M. Head-Gordon and K.B. Whaley, J. Chem. Phys. 145 (4), 044112 (2016).

• [29]

T. Zhang and F.A. Evangelista, J. Chem. Theory Comput. 12 (9), 4326–4337 (2016).

• [30]

A.A. Holmes, H.J. Changlani and C.J. Umrigar, J. Chem. Theory Comput. 12 (4), 1561–1571 (2016).

• [31]

S. Sharma, A.A. Holmes, G. Jeanmairet, A. Alavi and C.J. Umrigar, J. Chem. Theory Comput. 13 (4), 1595–1604 (2017).

• [32]

S.R. White, Phys. Rev. Lett. 69 (19), 2863–2866 (1992).

• [33]

G.K.L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465–81 (2011).

• [34]

K.H. Marti and M. Reiher, Phys. Chem. Chem. Phys. 13 (15), 6750 (2011).

• [35]

S. Wouters and D. Van Neck, Eur. Phys. J. D 68 (9), 272 (2014).

• [36]

R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang and G.K.L. Chan, J. Chem. Phys. 142 (3), 034102 (2015).

• [37]

Y. Kurashige and T. Yanai, J. Chem. Phys. 130 (23), 234114 (2009).

• [38]

Y. Kurashige, G.K.L. Chan and T. Yanai, Nat. Chem. 5 (8), 660–666 (2013).

• [39]

S. Wouters, T. Bogaerts, P. Van Der Voort, V. Van Speybroeck and D. Van Neck, J. Chem. Phys. 140 (24), 241103 (2014).

• [40]

J. Čížek, J. Chem. Phys. 45 (11), 4256 (1966).

• [41]

J.A. Parkhill, K. Lawler and M. Head-Gordon, J. Chem. Phys. 130 (8), 084101 (2009).

• [42]

J.A. Parkhill and M. Head-Gordon, J. Chem. Phys. 133 (2), 024103 (2010).

• [43]

R.J. Bartlett and M. Musiał, J. Chem. Phys. 125 (20), 204105 (2006).

• [44]

A.C. Hurley, J. Lennard-Jones and J.A. Pople, Proc. R. Soc. A Math. Phys. Eng. Sci. 220 (1143), 446–455 (1953).

• [45]

W.J. Hunt, P.J. Hay and W.A. Goddard III, J. Chem. Phys. 57 (2), 738 (1972).

• [46]

I.I. Ukrainskii, Theor. Math. Phys. 32 (3), 816–822 (1977).

• [47]

W.A. Goddard and L.B. Harding, Annu. Rev. Phys. Chem. 29 (1), 363–396 (1978).

• [48]

J. Cullen, Chem. Phys. 202 (2-3), 217–229 (1996).

• [49]

G.J.O. Beran, B. Austin, A. Sodt and M. Head-Gordon, J. Phys. Chem. A 109 (40), 9183–92 (2005).

• [50]

S. Lehtola, J. Parkhill and M. Head-Gordon, J. Chem. Phys. 145 (13), 134110 (2016).

• [51]

C.D. Sherrill, A.I. Krylov, E.F.C. Byrd and M. Head-Gordon, J. Chem. Phys. 109 (11), 4171 (1998).

• [52]

A.I. Krylov, C.D. Sherrill, E.F.C. Byrd and M. Head-Gordon, J. Chem. Phys. 109 (24), 10669 (1998).

• [53]

M. Bendikov, H.M. Duong, K. Starkey, K.N. Houk, E.A. Carter and F. Wudl, J. Am. Chem. Soc. 126 (24), 7416–7417 (2004).

• [54]

H.F. Bettinger, Pure Appl. Chem. 82 (4), 905–915 (2010).

• [55]

J. Hachmann, J.J. Dorando, M. Avilés and G.K.L. Chan, J. Chem. Phys. 127 (13), 134309 (2007).

• [56]

D.e. Jiang and S. Dai, J. Phys. Chem. A 112 (2), 332–335 (2008).

• [57]

D.H. Ess, E.R. Johnson, X. Hu and W. Yang, J. Phys. Chem. A 115 (1), 76–83 (2011).

• [58]

J.D. Chai, J. Chem. Phys. 136 (15), 154104 (2012).

• [59]

A.E. Torres, P. Guadarrama and S. Fomine, J. Mol. Model. 20 (5), 2208 (2014).

• [60]

H.F. Bettinger, C. Tönshoff, M. Doerr and E. Sanchez-Garcia, J. Chem. Theory Comput. 12 (1), 305–312 (2016).

• [61]

S. Rayne and K. Forest, Can. J. Chem. 94 (3), 251–258 (2016).

• [62]

S. Ghosh, C.J. Cramer, D.G. Truhlar and L. Gagliardi, Chem. Sci. 8 (4), 2741–2750 (2017).

• [63]

P. Rivero, C.A. Jiménez-Hoyos and G.E. Scuseria, J. Phys. Chem. B 117 (42), 12750–12758 (2013).

• [64]

Y. Yang, E.R. Davidson and W. Yang, Proc. Natl. Acad. Sci. 113 (35), E5098–E5107 (2016).

• [65]

J. Wilhelm, M. Del Ben and J. Hutter, J. Chem. Theory Comput. 12 (8), 3623–3635 (2016).

• [66]

G. Gidofalvi and D.A. Mazziotti, J. Chem. Phys. 129 (13), 134108 (2008).

• [67]

K. Pelzer, L. Greenman, G. Gidofalvi and D.A. Mazziotti, J. Phys. Chem. A 115 (22), 5632–5640 (2011).

• [68]

J. Fosso-Tande, D.R. Nascimento and A.E. DePrince, Mol. Phys. 114 (3-4), 423–430 (2015).

• [69]

J. Fosso-Tande, T.S. Nguyen, G. Gidofalvi and A.E. DePrince, J. Chem. Theory Comput. 12 (5), 2260–2271 (2016).

• [70]

M.S. Deleuze, L. Claes, E.S. Kryachko and J.P. François, J. Chem. Phys. 119 (6), 3106–3119 (2003).

• [71]

B. Hajgató, D. Szieberth, P. Geerlings, F. De Proft and M.S. Deleuze, J. Chem. Phys. 131 (22), 224321 (2009).

• [72]

B. Hajgató, M. Huzak and M.S. Deleuze, J. Phys. Chem. A 115 (33), 9282–9293 (2011).

• [73]

S.R. Yost and M. Head-Gordon, J. Chem. Phys. 145 (5), 054105 (2016).

• [74]

D. Casanova and M. Head-Gordon, Phys. Chem. Chem. Phys. 11 (42), 9779–90 (2009).

• [75]

C.U. Ibeji and D. Ghosh, Phys. Chem. Chem. Phys. 17 (15), 9849–9856 (2015).

• [76]

S. Battaglia, N. Faginas-Lago, D. Andrae, S. Evangelisti and T. Leininger, J. Phys. Chem. A 121 (19), 3746–3756 (2017).

• [77]

F. Aiga, J. Phys. Chem. A 116 (1), 663–669 (2012).

• [78]

Z. Qu, D. Zhang, C. Liu and Y. Jiang, J. Phys. Chem. A 113 (27), 7909–7914 (2009).

• [79]

D.W. Small, K.V. Lawler and M. Head-Gordon, J. Chem. Theory Comput. 10 (5), 2027–2040 (2014).

• [80]

J. Lee, D.W. Small, E. Epifanovsky and M. Head-Gordon, J. Chem. Theory Comput. 13 (2), 602–615 (2017).

• [81]

F. Plasser, H. Pašalić, M.H. Gerzabek, F. Libisch, R. Reiter, J. Burgdörfer, T. Müller, R. Shepard and H. Lischka, Angew. Chemie Int. Ed. 52 (9), 2581–2584 (2013).

• [82]

S. Horn, F. Plasser, T. Müller, F. Libisch, J. Burgdörfer and H. Lischka, Theor. Chem. Acc. 133 (8), 1511 (2014).

• [83]

S. Knippenberg, J.H. Starcke, M. Wormit and A. Dreuw, Mol. Phys. 108 (19-20), 2801–2813 (2010).

• [84]

U. Bozkaya, J.M. Turney, Y. Yamaguchi, H.F. Schaefer and C.D. Sherrill, J. Chem. Phys. 135 (10), 104103 (2011).

• [85]

H. Koch, H.J.A. Jensen, P. Jørgensen, T. Helgaker, G.E. Scuseria and H.F. Schaefer, J. Chem. Phys. 92 (8), 4924 (1990).

• [86]

T. Van Voorhis and M. Head-Gordon, Mol. Phys. 100 (11), 1713–1721 (2002).

• [87]

J. Nocedal and S. Wright, Numerical Optimization Springer Series in Operations Research and Financial Engineering (Springer New York, New York, 1999).

• [88]

T. Van Voorhis and M. Head-Gordon, J. Chem. Phys. 117 (20), 9190–9201 (2002).

• [89]

A. Edelman, T.A. Arias and S.T. Smith, SIAM J. Matrix Anal. Appl. 20 (2), 303–353 (1998).

• [90]

T. Abrudan, J. Eriksson and V. Koivunen, Signal Processing 89 (9), 1704–1714 (2009).

• [91]

S. Lehtola and H. Jónsson, J. Chem. Theory Comput. 9 (12), 5365–5372 (2013).

• [92]

S. Lehtola, M. Head-Gordon and H. Jónsson, J. Chem. Theory Comput. 12 (7), 3195–3207 (2016).

• [93]

N.H.F. Beebe and J. Linderberg, Int. J. Quant. Chem. 12, 683–705 (1977).

• [94]

O. Vahtras, J. Almlöf and M. Feyereisen, Chem. Phys. Lett. 213 (5-6), 514–518 (1993).

• [95]

J. Lehtola, M. Hakala, A. Sakko and K. Hämäläinen, J. Comput. Chem. 33 (18), 1572–1585 (2012).

• [96]

S. Lehtola, ERKALE - HF/DFT from Hel 2016. $<$https://github.com/susilehtola/erkale$>$.

• [97]

S. Lehtola and H. Jónsson, J. Chem. Theory Comput. 10 (2), 642–649 (2014).

• [98]

T. Sano, J. Mol. Struct. THEOCHEM 528 (1-3), 177–191 (2000).

• [99]

H. Koch, A. Sánchez de Merás and T.B. Pedersen, J. Chem. Phys. 118 (21), 9481 (2003).

• [100]

D.J. Thouless, Nucl. Phys. 21, 225–232 (1960).

• [101]

J. Pipek and P.G. Mezey, J. Chem. Phys. 90 (9), 4916 (1989).

• [102]

R.S. Mulliken, J. Chem. Phys. 23 (10), 1833 (1955).

• [103]

W.J. Hehre, R.F. Stewart and J.A. Pople, J. Chem. Phys. 51 (6), 2657 (1969).

• [104]

A. Köhn and J. Olsen, J. Chem. Phys. 122 (8), 084116 (2005).

• [105]

K.A. Brueckner, Phys. Rev. 100 (1), 36–45 (1955).

• [106]

T.H. Dunning, J. Chem. Phys. 90 (2), 1007 (1989).

• [107]

B.H. Lengsfield, J. Chem. Phys. 73 (1), 382–390 (1980).

• [108]

B.H. Lengsfield and B. Liu, J. Chem. Phys. 75 (1), 478–480 (1981).

• [109]

H.J. Werner, J. Chem. Phys. 74 (10), 5794 (1981).

• [110]

D.L. Yeager, D. Lynch, J. Nichols, P. Jørgensen and J. Olsen, J. Phys. Chem. 86 (12), 2140–2153 (1982).

• [111]

P. Jørgensen, P. Swanstrøm and D.L. Yeager, J. Chem. Phys. 78 (1), 347–356 (1983).

• [112]

H.J.A. Jensen and P. Jørgensen, J. Chem. Phys. 80 (3), 1204–1214 (1984).

• [113]

H.J. Werner and P.J. Knowles, J. Chem. Phys. 82 (11), 5053–5063 (1985).

• [114]

P. Knowles, K. Somasundram, N. Handy and K. Hirao, Chem. Phys. Lett. 113 (1), 8–12 (1985).

• [115]

H.J.A. Jensen, P. Jørgensen and H. Ågren, J. Chem. Phys. 87 (1), 451–466 (1987).

• [116]

M.R. Hoffmann, C. David Sherrill, M.L. Leininger and H.F. Schaefer III, Chem. Phys. Lett. 355 (1-2), 183–192 (2002).

• [117]

D. Ghosh, J. Hachmann, T. Yanai and G.K.L. Chan, J. Chem. Phys. 128 (14), 144117 (2008).

• [118]

Y. Ma, S. Knecht, S. Keller and M. Reiher, arXiv Prepr. 1611.05972 (2016).

• [119]

Q. Sun, J. Yang and G.K.L. Chan, Chem. Phys. Lett. (2017).

• [120]

M. Head-Gordon and J.A. Pople, J. Phys. Chem. 92 (11), 3063–3069 (1988).