James E T Smith

Bastien Mussard

Adam A Holmes

Sandeep Sharma

Cheap and near exact CASSCF with large active spaces

We use the recently-developed Heat-bath Configuration Interaction (HCI) algorithm as an efficient active-space solver to perform multi-configuration self-consistent field calculations (HCISCF) with large active spaces. We give a detailed derivation of the theory and show that difficulties associated with non-variationality of the HCI procedure can be overcome by making use of the Lagrangian formulation to calculate the HCI relaxed two body reduced density matrix. HCISCF is then used to study the electronic structure of butadiene, pentacene, and Fe-porphyrin. One of the most striking results of our work is that the converged active space orbitals obtained from HCISCF are relatively insensitive to the accuracy of the HCI calculation. This allows us to obtain nearly converged CASSCF energies with an estimated error of less than 1 mHa using the orbitals obtained from the HCISCF procedure in which the integral transformation is the dominant cost. For example, an HCISCF calculation on Fe-Porphyrin model complex with an active space of (44e, 44o) took only 412 seconds per iteration on a single node containing 28 cores, out of which 185 seconds were spent in the HCI calculation and the remaining 227 seconds were mainly used for integral transformation. Finally, we also show that active-space orbitals can be optimized using HCISCF to substantially speed up the convergence of the HCI energy to the Full CI limit because HCI is not invariant to unitary transformations within the active space.

I. INTRODUCTION

Many molecular systems with interesting electronic structure exhibit strong correlation, and as a result they are not well-described by the standard single-reference quantum chemical techniques, such as density functional theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Parr | Density-functional theory of atoms and molecules[END_REF] , Møller-Plesset perturbation theory [START_REF] Møller | Note on an approximation treatment for many-electron systems[END_REF] and coupled cluster theory [START_REF] Coester | Bound states of a many-particle system[END_REF][START_REF] Čížek | On the correlation problem in atomic and molecular systems. Calculation of wavefunction components in Ursell-type expansion using quantum-field theoretical methods[END_REF][START_REF] Čížek | Coupled cluster approach[END_REF][START_REF] Purvis | A full coupled-cluster singles and doubles model: the inclusion of disconnected triples[END_REF] . These highly multireference systems include transition metal complexes, excited states of conjugated organic molecules, and systems far from equilibrium, such as those near the breaking or forming of a covalent bond.

A popular paradigm for generating a multideterminant reference for such systems is to identify the "active" subset of the electrons and orbitals that are most important for a correct qualitative description of the molecular physics, and assume that the remaining orbitals are either always occupied ("core" orbitals) or never occupied ("virtual" orbitals). Once this partitioning of orbitals is chosen, the active electrons can be fully correlated to obtain the complete active-space configuration interaction (CASCI) ground-state wavefunction and energy. In addition to correlating the active-space electrons, the orbitals can also be optimized such that the active orbitals contain the most important degrees of freedom. The traditional heuristic method for finding these degrees of freedom is to minimize the CASCI energy, resulting in the complete active-space self-consistent field [START_REF] Roos | A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach[END_REF][START_REF] Roos | The complete active space SCF method in a fock-matrix-based super-CI formulation[END_REF][START_REF] Siegbahn | The complete active space SCF (CASSCF) method in a Newton-Raphson formulation with application to the HNO molecule[END_REF][START_REF] Roos | The Complete Active Space Self-Consistent Field Method and its Applications in Electronic Structure Calculations[END_REF] (CASSCF) algorithm.

Due to the exponential scaling of an exact correlated calculation with system size, the maximum active-space size of a CASSCF is about 16 electrons in 16 orbitals, although recent developments have made it possible to go up to 20 electrons in 20 orbitals on massively parallel machines [START_REF] Vogiatzis | Pushing Configuration-Interaction to the Limit: Towards Massively Parallel MCSCF Calculations[END_REF] . Naturally, there is great interest in developing algorithms that can overcome this limit in a sys-tematically improvable way. There are many algorithms that can be used as approximate active-space solvers, including the density matrix renormalization group restricted [START_REF] Malmqvist | The restricted active space self-consistent-field method, implemented with a split graph unitary group approach[END_REF][START_REF] Celani | Multireference perturbation theory for large restricted and selected active space reference wave functions[END_REF] , and generalized [START_REF] Ma | The generalized active space concept in multiconfigurational selfconsistent field methods[END_REF] active space methods, symmetry broken coupled cluster theory [START_REF] Gomez | Singlet-paired coupled cluster theory for open shells[END_REF][START_REF] Qiu | Projected coupled cluster theory[END_REF] , reduced density matrix approaches [START_REF] Nakata | Density matrix variational theory: Application to the potential energy surfaces and strongly correlated systems[END_REF][START_REF] Mazziotti | Quantum chemistry without wave functions: two-electron reduced density matrices[END_REF][START_REF] Valdemoro | In Reduced-density-matrix mechanics with applications to many-electron atoms and molecules; Mazziotti[END_REF] , various flavors of selected configuration interaction followed by perturbation theory [START_REF] Ivanic | Identification of deadwood in configuration spaces through general direct configuration interaction[END_REF][START_REF] Huron | Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions[END_REF][START_REF] Buenker | Individualized configuration selection in CI calculations with subsequent energy extrapolation[END_REF][START_REF] Evangelista | IMPROVED ClPSl ALGORITHM[END_REF][START_REF] Harrison | Approximating full configuration interaction with selected configuration interaction and perturbation theory[END_REF][START_REF] Steiner | The efficient treatment of higher excitations in CI calculations. A comparison of exact and approximate results[END_REF][START_REF] Wenzel | Multireference Basis-Set Reduction[END_REF][START_REF] Neese | A spectroscopy oriented configuration interaction procedure[END_REF][START_REF] Abrams | Important configurations in configuration interaction and coupled-cluster wave functions[END_REF][START_REF] Bytautas | A priori identification of configurational deadwood[END_REF][START_REF] Evangelista | A driven similarity renormalization group approach to quantum many-body problems[END_REF][START_REF] Knowles | Compressive sampling in configuration interaction wavefunctions[END_REF][START_REF] Schriber | Communication: An adaptive configuration interaction approach for strongly correlated electrons with tunable accuracy[END_REF][START_REF] Liu | ICI: Iterative CI toward full CI[END_REF][START_REF] Caffarel | Using CIPSI Nodes in Diffusion Monte Carlo[END_REF][START_REF] Garniron | Hybrid stochastic-deterministic calculation of the secondorder perturbative contribution of multireference perturbation theory[END_REF] (SCI+PT), and stochastic and semistochastic methods such as Full CI quantum Monte Carlo [START_REF] Booth | Fermion Monte Carlo without fixed nodes: a game of life, death, and annihilation in Slater determinant space[END_REF][START_REF] Cleland | Communications: Survival of the fittest: accelerating convergence in full configuration-interaction quantum Monte Carlo[END_REF][START_REF] Petruzielo | Semistochastic projector monte carlo method[END_REF][START_REF] Thomas | Stochastic Multiconfigurational Self-Consistent Field Theory[END_REF] (FCIQMC).

Here, we use the recently-developed Heat-bath Configuration Interaction [START_REF] Holmes | Heatbath Configuration Interaction: An efficient selected CI algorithm inspired by heat-bath sampling[END_REF][START_REF] Sharma | Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory[END_REF] (HCI) algorithm, an efficient SCI+PT algorithm, as an approximate active-space solver, in order to perform efficient CASSCF-like calculations with large active spaces. We call the resulting algorithm Heat-bath Configuration Interaction Self-Consistent Field (HCISCF). Similar extensions have already been presented with other theories such as DMRG-SCF [START_REF] Ma | Second-Order Self-Consistent-Field Density-Matrix Renormalization Group[END_REF][START_REF] Ghosh | Orbital optimization in the density matrix renormalization group, with applications to polyenes and β-carotene[END_REF][START_REF] Zgid | The density matrix renormalization group self-consistent field method: Orbital optimization with the density matrix renormalization group method in the active space[END_REF] , FCIQMC-SCF [START_REF] Manni | Combining the Complete Active Space Self-Consistent Field Method and the Full Con fi guration Interaction Quantum Monte Carlo within a Super-CI Framework , with Application to Challenging Metal-Porphyrins[END_REF][START_REF] Thomas | Stochastic Multiconfigurational Self-Consistent Field Theory[END_REF] . However, unlike DMRG-SCF, formulating HCISCF is made more complicated by the fact that HCI in not a variational method since it performs second-order perturbation theory. This implies that HCI energy is not stationary with respect to the variations of its zeroth-order wavefunction parameters c ¯, i.e. ∂EHCI ∂c ¯ = 0. Thus, the first-order change in energy due to change in orbital rotation parameters κ κ κ is given by

dE HCI dκ κ κ = ∂E HCI ∂κ κ κ + ∂E HCI ∂c ¯dc d κ κ κ , (1)
where the second term does not vanish, as it does in DMRG-SCF and CASSCF. Here, and in the rest of the article, bold-face letters represent vectors. Equation 1suggests that to calculate the gradient of the HCI energy with respect to κ κ κ, one has to calculate the derivatives of the wavefunction with respect to each of the O(n) parameters in κ κ κ, where n is the number of basis functions (here we have assumed that the number of active and core orbitals is much smaller than the number of basis fuctions). In this work, we overcome this prohibitive expense by using the technique of Lagrange multipliers, which is also known as the z-vector method [START_REF] Handy | On the evaluation of analytic energy derivatives for correlated wave functions[END_REF] . The Lagrange multiplier technique replaces calculating the derivatives of the wavefunction with respect to all the parameters, with calculating only a single set of Lagrange multipliers. It is worth mentioning that similar orbital optimization over a perturbatively-corrected energy is performed in orbital-optimized Møller-Plesset perturbation theory (OO-MP2). [START_REF] Lochan | Orbital-optimized opposite-spin scaled second-order correlation: An economical method to improve the description of open-shell molecules[END_REF][START_REF] Neese | Assessment of orbital-optimized, spincomponent scaled second-order many-body perturbation theory for thermochemistry and kinetics[END_REF] A considerable simplification in the HCISCF theory arises if only the HCI variational energy is minimized, rather than the sum of the variational plus perturbative correction. In this case, the second term of Equation 1is zero because the variational HCI energy is stationary with respect to the wavefunction parameters c ¯. We have implemented such a theory as well and call it vHCISCF to distinguish it from the full HCISCF.

We apply this new algorithm to calculate the ground states for butadiene, the pentacene monomer, and the Fe(II)-Porphyrin complex, abbreviated as Fe(P), as well as the excited states for the latter two systems. For pentacene, we calculate the energies of the ground ( 1 A g ) and lowest triplet ( 3 B 2u ) states, and compare the gaps to experimental [START_REF] Burgos | Heterofission in pentacene-doped tetracene single crystals[END_REF][START_REF] Biermann | Diels-Alder Reactivity of Polycyclic Aromatic Hydrocarbons III New Experimental and Theoretical Results[END_REF][START_REF] Biermann | Diels-Alder Reactivity of Polycyclic Aromatic Hydrocarbons. 1. Acenes and Benzologs[END_REF] and theoretical [START_REF] Hachmann | The radical character of the acenes: A density matrix renormalization group study[END_REF][START_REF] Dorando | Targeted excited state algorithms[END_REF][START_REF] Hajgató | A benchmark theoretical study of the electronic ground state and of the singlet-triplet split of benzene and linear acenes[END_REF][START_REF] Zimmerman | Singlet fission in pentacene through multi-exciton quantum states[END_REF][START_REF] Yuki | Theoretical Study of the * Excited States of Oligoacenes: A Full -Valence DMRG-CASPT2 Study[END_REF][START_REF] Coto | Low-lying electronic excited states of pentacene oligomers: A comparative electronic structure study in the context of singlet fission[END_REF][START_REF] Yang | Nature of ground and electronic excited states of higher acenes[END_REF] results of other CAS-based approaches. Fe(P) is a challenging electronic-structure problem, in which theory and experiments disagree on the symmetry of the ground state wavefunction. [START_REF] Manni | Combining the Complete Active Space Self-Consistent Field Method and the Full Con fi guration Interaction Quantum Monte Carlo within a Super-CI Framework , with Application to Challenging Metal-Porphyrins[END_REF][START_REF] Hachmann | The radical character of the acenes: A density matrix renormalization group study[END_REF][START_REF] Hajgató | A benchmark theoretical study of the electronic ground state and of the singlet-triplet split of benzene and linear acenes[END_REF][START_REF] Collman | Stereochemistry, and Structure-Related Properties of alpha beta gamma delta Tetraphenylporphinatoiron(II)[END_REF][START_REF] Goff | Nuclear Magnetic Resonance Investigation of Magnetic and Electronic Properties of "Intermediate Spin" Ferrous Porphyrin Complexes[END_REF][START_REF] Kitagawa | The Resonance Raman Spectra of Intermediate-spin Ferrous Porphyrin[END_REF][START_REF] Dolphin | Synthesis and Mossbauer spectra of octaethylporphyrin ferrous complexes[END_REF][START_REF] Lang | Mössbauer effect study of the magnetic properties of S =1 ferrous tetraphenylporphyrin[END_REF][START_REF] Boyd | Paramagnetic Anisotropy, Average Magnetic Susceptibility, and Electronic Structure of Intermediate-Spin S = 1 (5,10,15,20-Tetraphenylporphyrin)iron(II)[END_REF][START_REF] Mispelter | Proton magnetic resonance characterization of the intermediate ( S =1) spin state of ferrous porphyrins[END_REF][START_REF] Strauss | Comparison of the Molecular and Electronic Structures of (2,3,7,8,12,13,17,18-Octaethylporphyrinato)iron( 11) and 4207 porphyrinato)iron( 11) Steven[END_REF][START_REF] Choe | Theoretical study of the electronic ground state of iron (II) porphine[END_REF][START_REF] Choe | Theoretical study of the electronic ground state of iron II . . . porphine . II[END_REF][START_REF] Pierloot | The CASPT2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes[END_REF][START_REF] Radon | 2 to Heme by Density Functional and Multireference ab Initio Calculations[END_REF][START_REF] Vancoillie | Multiconfigurational second-order perturbation theory restricted active space (RASPT2) studies on mononuclear first-row transition-metal systems[END_REF][START_REF] Phung | Cumulant Approximated Second-Order Perturbation Theory Based on the Density Matrix Renormalization Group for Transition Metal Complexes: A Benchmark Study[END_REF] Due to the challenging nature of this system, we investigate the effects of basis set and choice of active-space orbitals on the gap between the 5 A g and 3 B 2g states.

The remainder of the paper is organized as follows. In Section II, we briefly review the HCI algorithm with the aim of presenting the equations that will be used in the formulation of the HCISCF and vHCISCF algorithms. In Section III, we derive the working equations of the HCISCF and vHCISCF algorithms. In Section IV, we discuss the implementation and practical aspects of running an HCISCF calculation. Finally, in Section V, we report our HCISCF calculations for butadiene, pentacene, and Fe(P), and compare them to experiments and previous calculations on these systems.

II. THE HCI ALGORITHM

Heat-bath Configuration Interaction (HCI), like other SCI+PT schemes, consists of two stages:

• a variational stage in which a set of important determinants is iteratively selected and used to compute a variational wavefunction and energy, and

• a perturbative stage in which the second-order correction to the variational energy is computed using multi-reference Epstein-Nesbet perturbation theory.

A. Variational Stage

In the variational stage, a set of determinants V is iteratively generated and their coefficients are variationally optimized to find a multi-determinantal wavefunction and energy. At a given iteration µ, the current set of determinants in noted V µ and the current variational wavefunction is given by

|Ψ µ 0 = |Di ∈V µ c µ i |D i , (2)
where c ¯µ are the current CI coefficients. The variational space is augmented (see Figure 1a) as

V µ+1 = V µ +C µ ( 1 )
where 1 is a user-defined parameter. The space of connected determinants C µ ( 1 ) is defined as

C µ ( 1 ) = |D a | f µ (|D a ) > 1 , (3)
where f µ (|D a ) is a non-negative "importance" function, which can vary between different SCI+PT schemes.

The HCI importance function is:

f µ HCI (|D a ) = max |Di ∈V µ |H ai c µ i | , (4)
where H ai is the Hamiltonian matrix element between |D a and |D i . In other words, the space C µ ( 1 ) is composed of all determinants |D a for which |H ai c µ i | > 1 is true for at least one determinant |D i in the current V µ .

By contrast, the importance function of CIPSI and its variants is inspired by perturbation theory (hence the name "CI by Perturbatively Selecting Iteratively"), and is given by

f µ CIPSI (|D a ) = |Di ∈V µ H ai c µ i E µ 0 -E a , (5)
where E a is the energy of the determinant |D a and E µ 0 is the current ground-state energy.

The advantage of the HCI importance function over the CIPSI one is twofold. First, the relevant information for the HCI importance function for a given reference (apart from its coefficient) is simply the magnitude of the connecting matrix element, |H ai |. Since the vast majority of excitations are double excitations, whose magnitudes are simple functions of only the orbitals involved (and not the specific determinants involved), all the necessary information can be sorted and stored prior to the run. As a result, the doubly-excited determinants meeting the HCI criterion can be generated without generating lists of candidates first, as would be required for CIPSI (see Figure 1b). Second, the HCI importance function FIG. 1. Outline of the HCI scheme. The dots represent the determinants in a Hilbert space; the connections between them are randomized for the purpose of this sketch. Figure (a) shows that at a given iteration µ the current variational space is V µ (encompassing the cyan dots) and is augmented by the space of connected determinant C µ ( 1) (that includes the blue dots). This sequential aggregation of determinants is done iteratively until some designated convergence threshold is met. The dotted lines represent possible excitations from every pair of occupied orbitals in a given cyan determinant |Di of a current V µ . A determinant |Da is generated and included in C µ ( 1) if the magnitude of Hia = Di| Ĥ|Da is greater than a certain i-dependent threshold,

1/ |c µ i |.
Hence, for all pairs of occupied orbitals in a given determinant (i.e. for all dotted lines) the algorithm will only browse the sorted list of magnitudes up to that threshold and generate the corresponding blue determinants |Da : no computational time is lost on generating determinants (here in green) that will not be included in C µ ( 1). can be utilized in constant time, whereas the CIPSI importance function requires performing a sum in the numerator (which may require communication across cores in a parallel run) and evaluating a diagonal element in the denominator.

After the addition of the determinants in C( 1 ), the new ground-state energy and wavefunction within the new variational space are obtained at each iteration by the Davidson procedure, which amounts to minimizing the energy functional

E[c ¯µ, E µ 0 ; V µ ] = Ψ µ 0 | Ĥ0 |Ψ µ 0 -E µ 0 ( Ψ µ 0 |Ψ µ 0 -1), ( 6
)
with respect to the coefficients c ¯µ of |Ψ µ 0 and to the Lagrange multiplier E µ 0 ensuring that the new wavefunction remains normalized, in the context of a space V µ which is fixed at this point.

B. Perturbative stage

Once a converged variational space V, variational coefficients c ¯, variational wavefunction |Ψ 0 , and corresponding variational energy E 0 have been obtained using the algorithm in the previous section, multireference perturbation theory can be performed to estimate the Full CI energy. We use the Epstein-Nesbet partitioning of the Hamiltonian by defining the zeroth-order Hamiltonian

Ĥ0 = |Di ,|Dj ∈V H ij |D i D j | + |Da / ∈V H aa |D a D a |, (7)
and perturbation V = Ĥ -Ĥ0 . With this partitioning, the second-order perturbative energy correction to the variational energy is given by

E 2 = |Da ∈C 1 E 0 -H aa   |Di ∈V H ai c i   2 , (8)
where C denotes the set of determinants that are connected to at least one determinant in V by a non-zero Hamiltonian matrix element. The vast majority of the contributions in the double sum are negligibly small, and can be discarded without significant loss of accuracy. HCI therefore approximates the perturbative energy correction as

E 2 ( 2 ) = |Da ∈C( 2) 1 E 0 -H aa   ( 2) |Di ∈V H ai c i   2 , (9)
where the symbol

( 2) denotes a "screened sum" in which terms smaller in magnitude than a user-defined parameter 2 are discarded. Hence, the important terms involve determinants in the set C( 2 ) that are connected to at least one determinant in V by a Hamiltonian matrix element larger in magnitude than 2 / |c i |. In the limit of 2 → 0, the exact perturbation correction is recovered. In order to obtain a good approximation to the perturbative correction 2 must be much smaller than 1 , but a nearly exact approximation to the perturbative correction can be obtained at a much reduced cost by choosing a small but non-zero 2 parameter.

The naive evaluation of Eq. 9 requires one to simultaneously store the entire set of determinants C( 2 ) in order to combine contributions to the sum before they are squared. Even with the use of a non-zero parameter 2 , the number of determinants in C( 2) can be extremely large, resulting in a memory bottleneck.

The second-order energy correction can alternatively be estimated stochastically with the same accuracy, circumventing the memory bottleneck at the cost of the intrusion of an unbiased stochastic error. We refer the reader to Ref. 65 for details on this unbiased sampling procedure, which uses a sampled subset of the determinants in the zeroth-order wavefunction.

The stochastic error can be reduced by implementing a semistochatic algorithm, in which a deterministic perturbative calculation is performed using a large parameter d 2 to avoid the memory bottleneck, and the error is corrected stochastically using a tight parameter 2 , as follows (note the superscripts D and S for expressions that are evaluated deterministically or stochastically, respectively):

E 2 ( 2 ) = E D 2 ( d 2 ) + E S 2 ( 2 ) -E S 2 ( d 2 ) . ( 10
)
The key point is that E S 2 ( 2 ) and E S 2 ( d 2 ) are evaluated using the same set of sampled determinants; consequently, their errors are highly correlated and the difference between the two calculations (shown in the square brackets in Eq 10) has a much reduced stochastic noise.

For the formulation of the self-consistent procedure, it is useful to point out that in place of Equation 8, the second-order correction to the energy can equivalently be obtained by minimizing the Hylleraas functional

H[d ¯; c ¯, V, C( 2 )] = Ψ 2 1 | Ĥ0 -E 0 |Ψ 2 1 + 2 Ψ 2 1 | V |Ψ 0 , (11)
with respect to the coefficients d ¯of the first-order wavefunction

|Ψ 2 1 = a∈C( 2) d a |D a , (12)
with c ¯, V and C( 2 ) held fixed. At its minimum, the Hylleraas functional gives the optimal values of d ¯, and the value of the functional is equal to the second-order correction evaluated using Eq 8. This will be used in Section III to derive formulas for HCISCF calculations.

It should be mentioned at this point that E 2 ( 2 ) is not a strict upper limit to E 2 , since the introduction of the non-zero parameter 2 not only truncates the size of the space C( 2), but also changes the perturbation V by ignoring small matrix elements. The Hylleraas functional formulation of perturbation theory seems to show that merely truncating the size of C( 2 ) would provide a variational upper bound to the second-order energy i.e. E 2 ( 2 ) ≥ E 2 ; however, because the perturbation V is simultaneously changed, the strict variationality is destroyed. Instead, given 2 and the zeroth-order wavefunction (c ¯and V), the inequality

H[d ¯; c ¯, V, C( 2 )] ≥ E 2 ( 2 )
holds, where the equality holds at its minimum.

III. HCI SELF-CONSISTENT FIELD

The HCISCF algorithm presented here is designed as a CASSCF-like procedure, in which the CASCI is replaced by an HCI calculation in the active space. However, an important distinction between HCISCF and CASSCF is that unlike CASCI, the unconverged active-space HCI energy is in general not invariant to active-active rotations. For instance, using natural orbitals as opposed to canonical orbitals can result in a drastically improved convergence of HCI to the CASCI or Full CI limit. This flexibility will be utilized to optimize the active-space orbitals in order to accelerate convergence.

Since HCI uses both a variational step and secondorder perturbation theory, we have the choice of optimizing the total HCI energy or just the variational HCI energy with respect to the orbital coefficients. In this section we describe both these procedures and call them HCISCF and vHCISCF respectively.

A. HCISCF

The HCI energy is given by the sum of the zeroth-order and second-order energies calculated by minimizing respectively the energy functional in Eq 6 with respect to c ¯and E µ 0 , and the Hylleraas functional in Eq 11 with respect to c ¯, E 0 , and d ¯. Formally, the minimization of these functionals is performed by setting to zero their partial derivatives with respect to the parameters c ¯, E 0 and d ¯. Thus the most natural way of deriving the HCISCF procedure is by setting to zero the derivative of the HCI functional with respect to the parameter κ κ κ of the orbital coefficients. The HCI functional reads:

E HCI [κ κ κ, c ¯, E 0 , d ¯] = E[κ κ κ, c ¯, E 0 ] + H[κ κ κ, c ¯, d ¯], (13)
where the dependence on the parameters κ κ κ is shown and the dependence on V and C is dropped. Its derivative with respect to κ κ κ reads:

dE HCI dκ κ κ = ∂E HCI ∂κ κ κ + ∂E HCI ∂c ¯dc d κ κ κ + ∂E HCI ∂E 0 dE 0 dκ κ κ + ∂E HCI ∂d ¯dd d κ κ κ . ( 14
)
As mentioned in the introduction, this formulation requires the calculation of the derivatives of the HCI pa-rameters with respect to κ κ κ, which makes it an impractical approach. Instead, we introduce the Lagrangian

L[κ κ κ, c ¯, E 0 , d ¯, λ λ λ c , λ λ λ d ] = E[κ κ κ, c ¯, E 0 ] + H[κ κ κ, c ¯, d ¯] + λ λ λ † c ∂E ∂c ¯ † + λ λ λ † d ∂H ∂d ¯ † , (15)
which is a function of the variables κ κ κ, c ¯, and d ¯, and the set of Lagrange multipliers E 0 , λ λ λ c , and λ λ λ d . If the partial derivatives of the Lagrangian with respect to the Lagrange multipliers are set to zero, one obtains the governing equations of HCI from which the optimal values of c ¯and d ¯can be recovered. Note that at these optimal values of c ¯and d ¯, the Lagrangian is by construction equal to the HCI energy functional. Furthermore, the Lagrange multipliers λ λ λ c and λ λ λ d can be used to impose the stationarity of the Lagrangian with respect to c ¯and d ¯. Thus, the Lagrangian can be made stationary with respect to all its parameters, allowing the calculation of the derivative of the HCI energy with respect to κ κ κ as

dE HCI dκ κ κ = dL dκ κ κ = ∂L ∂κ κ κ , (16)
leading to the desired minimization of the HCI energy with respect to κ κ κ under the constraint that ∂E ∂c ¯= 0 and

∂H ∂d ¯= 0.
Let us first derive the equations to obtain the Lagrange multipliers:

∂L ∂d ¯= ∂H ∂d ¯+ λ λ λ † d ∂ 2 H ∂d ¯∂d ¯ † = 0 (17) ∂L ∂c ¯= ∂E ∂c ¯+ ∂H ∂c ¯+ λ λ λ † c ∂ 2 E ∂c ¯∂c ¯ † + λ λ λ d ∂ 2 H ∂c ¯∂d ¯ † = 0 (18)
Here, we have assumed that all parameters are real numbers, although extension to complex numbers is straightforward without any additional complications. From Equation 17, we can infer that λ λ λ d = 0 because ∂H ∂d ¯= 0. Equation 18 simplifies to

∂H ∂c ¯+ λ λ λ † c ∂ 2 E ∂c ¯∂c ¯ † = 0, (19)
because λ λ λ d = 0 and ∂E ∂c ¯= 0. This equation can be solved to evaluate λ λ λ c (in the analytic gradient theory, this corresponds to the z-vector equation).

Finally, the gradient of the HCI energy with respect to κ κ κ is given by

∂L ∂κ κ κ = c ¯ † ∂H 0 ∂κ κ κ c ¯+ d ¯ † ∂H 0 ∂κ κ κ d ¯+ 2d ¯ † ∂V ∂κ κ κ c ¯+ 2λ λ λ † c ∂H 0 ∂κ κ κ c = ijkl ∂H 0,ijkl ∂κ κ κ Γ c,c ijkl + ∂H 0,ijkl ∂κ κ κ Γ d,d ijkl + 2 ∂V ijkl ∂κ κ κ Γ d,c ijkl + 2 ∂H 0,ijkl ∂κ κ κ Γ λc,c ijkl , (20)
where the partial derivatives of the two-body integrals (H 0,ijkl and V ijkl ) with respect to κ κ κ are contracted with the transition two-body reduced density matrices (Γ Γ Γ; see Computational Details) between two states shown as superscript (the one-body terms are not shown here to avoid proliferation of terms). A careful look at the equation reveals that the first term is the reduced density matrix of the variational wavefunction and the second two terms are the unrelaxed reduced density matrices of the perturbative correction, while the last term arises due to the change in the second-order energy with the relaxation of the zeroth-order wavefunction as the orbitals are optimized. The partial derivatives of the two-body in-tegrals with respect to κ κ κ are calculated using the usual techniques, which are described in detail in Ref. 98.

B. vHCISCF

When only the variational energy E HCI is optimized instead of the total HCI energy, one can calculate the FIG. 2.

The basic scheme for the HCISCF procedure is outlined here. The HCISCF module of PySCF is used to interface the Dice program with PySCF. Each iteration consists of a single Dice run, which returns the energy and the reduced density matrices to PySCF. These are in turn used to update the orbitals and the active space two electron integrals, which are passed on to Dice.

gradient with respect to κ κ κ as

dE HCI dκ κ κ = ∂E ∂κ κ κ + ∂E ∂c ¯dc d κ κ κ = ∂E ∂κ κ κ = ijkl ∂H 0,ijkl ∂κ κ κ Γ c,c ijkl . (21)
Thus, to calculate the energy gradient, one does not need to evaluate the Lagrange multipliers, and the simple variational two-body reduced density matrix Γ c,c ijkl is sufficient.

IV. COMPUTATIONAL DETAILS

We briefly describe the calculation of the reduced density matrices encountered in Eq. (20). The density matrices Γ c,c ijkl and Γ λc,c ijkl are strictly limited to states that contain determinants from the variational space V. During the variational calculation, all single and double excitations between pairs of determinants in the variational wavefunction are generated in order to evaluate the Hamiltonian. With this data available, the reduced density matrices can be calculated in a straightforward manner by looping over such connections and accumulating the contributions. The reduced density matrix Γ d,d ijkl contains the first-order wavefunction as both the bra and the ket; however, because Ĥ0 only contains the diagonal elements in the space of the connections C( 2 ), the reduced density matrix simply reads:

Γ d,d ijkl = |Da ∈C( 2) d 2 a D a |a † i a † j a k a l |D a . (22)
Finally, for the evaluation of Γ d,c ijkl , all the determinants |D i in V that are connected by a Hamiltonian matrix element |H ia | > 2 to a determinant |D a are stored in a list. These connections are used at the end of the perturbative calculations to evaluate Γ d,c ijkl as

Γ d,c ijkl = |Di ∈V |Da ∈C( 2) |Hia|> 2 c i d a D i |a † i a † j a k a l |D a . ( 23
)
Thus out of all the different reduced density matrices only the evaluation of Γ d,c

ijkl adds a non-trivial memory cost over the HCI calculation, because not only do we have to store all the determinants |D a in the connected space C( 2 ), but for each of these determinants we also have to store a list of variational determinants that are connected to them. This memory bottleneck can again be overcome with the use of semistochastic perturbation theory, but in the paper we have not done so and have limited ourselves to using a value of 2 for which the deterministic calculations can be performed. Although using the stochastic perturbation theory poses no challenge, we will see in the results section that the optimized orbitals converge relatively rapidly even when only a vH-CISCF with a loose 1 is performed.

All calculations in this work were performed using our HCISCF module in the PySCF software package, which interfaces PySCF with the Dice program. Dice is used to calculate the HCI energy and reduced density matrices. At each iteration of the HCISCF procedure, PySCF updates the orbital coefficients by calculating the energy gradient using the reduced density matrices obtained from the previous iteration. The updated orbitals are used to calculate the active-space Hamiltonian with which Dice calculates the HCI energy and reduced density matrices. This procedure is illustrated in Fig 2 and is carried out until convergence. After convergence of the orbitals and CI coefficients, a final HCI calculation is performed where we use a smaller 1 and 2 to obtain near full configuration interaction (FCI) energy in the optimized active space.

V. RESULTS AND DISCUSSION

We perform benchmark calculations on three different systems: butadiene, pentacene, and Fe-porphyrin. In addition to getting system-specific information, these calculations are meant to provide heuristics for running HCISCF calculations. We will investigate the following aspects:

• How tightly do we need to converge the HCI energies during the HCISCF procedure to obtain CASSCF-quality active-space orbitals?

• Do vHCISCF and HCISCF converge to similar or substantially different active-space orbitals?

• How much energy relaxation can be obtained just by optimizing the active-space orbitals, while keeping the active space itself fixed? Such a calculation where only the active-active orbital rotations are allowed during SCF will be called aHCISCF.

In the result section we use the acronyms vHCI and SHCI respectively to indicate the energies of the variational step and calculations where semistochastic perturbation theory was used. All SHCI calculations performed here have stochastic noise of less than 0.05 mHa.

A. Butadiene

The HCI calculations on butadiene were performed with the same ANO-L-pVDZ basis set and geometry as the one used in Ref. 99. All electrons except the 1s orbitals were fully correlated to give an active space of (82o, 22e).

In this section, calculations were done using either Hartree-Fock canonical orbitals or optimized orbitals obtained from an aHCISCF calculation (where only the active-active rotations are allowed) with a relatively loose 1 = 3 × 10 -4 Ha. HCISCF energies typically are quadratically convergent because PySCF is able to perform pseudo-second-order optimization by estimating the Hessian, but the rate of convergence of aHCISCF calculations become substantially worse because of the strong coupling between the CI coefficients and orbital rotation parameters. The cost of performing the orbital optimization is, however, more than made up for by the improved convergence of the HCI energies with the optimized orbitals. Table I shows that with approximately the same number of determinants in the variational space, the vHCI energy is more than 24 mHa lower when optimized active space orbitals are used as opposed to the Hartree-Fock canonical orbitals. The effect is smaller for the full HCI energy, where the relaxation is only of 1.4 mHa, but one can observe that the rate of convergence of the full HCI energies with the number of determinants in the variational space still substantially improves when optimized orbitals are used. As a result, the HCI energies calculated using the optimized orbitals can be accurately extrapolated to the FCI limit.

As in a previous paper 100 , we estimate the Full CI limit by extrapolating the SHCI energy to the E 2 → 0 limit. In this work, we have improved the procedure in the previous paper to obtain better extrapolated values. We first perform a single SHCI calculation using the smallest 1 allowed by the available computational resources to obtain the vHCI and SHCI energies. For butadiene an 1 of 10 -5 Ha was used and this resulted in 56276315 determinants in the variational space. In the next step 33% of the least important determinants by their coefficient in the variational wavefunction were discarded and a SHCI FIG. 3. Extrapolation of the SHCI total energy to the FCI limit for the ground state of butadiene with an (82o, 22e) active space, using optimized aHCISCF orbitals. The dotted green line shows the DMRG energy calculated with a bond dimension of M=6000 and is believed to be converged to better than 1 mHa accuracy.

calculation was performed with this smaller variational space to obtain higher vHCI and SHCI energies. This procedure was iterated 3 more times to obtain a set of SHCI and vHCI energies. We then plot the total SHCI energy versus the PT correction (E 2 ) and perform a linear extrapolation to E 2 → 0, as shown in Figure 3. This procedure was inspired by the one used in DMRG calculations, where a large M (number of retained renormalized states) calculation is performed, a few subsequent sweeps are carried out with progressively smaller values of M and the energies and discarded weights obtained from these calculations are used to perform a linear extrapolation to a zero discarded weight.

B. Pentacene

Linear acenes, such as tetracene and pentacene, are promising candidates for singlet fission application because the gap to the lowest lying triplet state (T0) is roughly half of that of the first singlet excited state (S1). They also show great promise as organic semiconductors due to their large carrier mobility and low production cost. With recent applications to light emitting diodes, photovoltaic cells, and field effect transistors, these systems have been the subject of many theoretical and experimental studies. [START_REF] Burgos | Heterofission in pentacene-doped tetracene single crystals[END_REF][START_REF] Biermann | Diels-Alder Reactivity of Polycyclic Aromatic Hydrocarbons III New Experimental and Theoretical Results[END_REF][START_REF] Biermann | Diels-Alder Reactivity of Polycyclic Aromatic Hydrocarbons. 1. Acenes and Benzologs[END_REF][START_REF] Hachmann | The radical character of the acenes: A density matrix renormalization group study[END_REF][START_REF] Dorando | Targeted excited state algorithms[END_REF][START_REF] Hajgató | A benchmark theoretical study of the electronic ground state and of the singlet-triplet split of benzene and linear acenes[END_REF][START_REF] Zimmerman | Singlet fission in pentacene through multi-exciton quantum states[END_REF][START_REF] Yuki | Theoretical Study of the * Excited States of Oligoacenes: A Full -Valence DMRG-CASPT2 Study[END_REF][START_REF] Coto | Low-lying electronic excited states of pentacene oligomers: A comparative electronic structure study in the context of singlet fission[END_REF][START_REF] Yang | Nature of ground and electronic excited states of higher acenes[END_REF] The ground and excited state of acenes become progressively more multireference as their chain length is increased, and multireference methods are necessary to obtain an accurate description of their gap.

We calculate the energies of the pentacene 1 A g ground state and 3 B 2u lowest triplet state [START_REF] Hajgató | A benchmark theoretical study of the electronic ground state and of the singlet-triplet split of benzene and linear acenes[END_REF][START_REF] Zimmerman | Singlet fission in pentacene through multi-exciton quantum states[END_REF][START_REF] Yuki | Theoretical Study of the * Excited States of Oligoacenes: A Full -Valence DMRG-CASPT2 Study[END_REF][START_REF] Yang | Nature of ground and electronic excited states of higher acenes[END_REF] in the cc-pVDZ basis set on the optimized singlet and triplet geometries. TABLE I. Results of HCI calculations on butadiene with an active space of (82o, 22e) with the ANO-L-pVDZ basis set. Two sets of HCI calculations were performed: calculations using the canonical Hartree-Fock orbitals ("Canonical") and calculations using the aHCISCF optimized orbitals ("Optimized"). Nvar, vHCI and SHCI respectively show the number of determinants included in the variational space, the variational energy (+155.0 Ha) and the final SHCI energy (+155.0 Ha) with stochastic error bar of 0.04 mHa. The results of the calculations using the optimized orbitals were used to extrapolate the SHCI energy (see text for more details). The extrapolated SHCI energy is shown along with results obtained using other methods. The PT calculation in HCISCF was performed with 2 = 10 -5 Ha. We also report extrapolated SHCI energies (see text) using the optimized active space orbitals obtained from the SCF calculations. It is interesting to note that although the HCISCF energies are far from converged, the SHCI energies agreement to within 0.1 mHa indicates that the active space orbitals are most likely converged. This allows us to calculate both vertical and relaxed (well-to-well) 1 A g → 3 B 2u excitation energies. The active space chosen for all calculations is composed of the 11 π and 11 π * orbitals. We begin by assessing the effect of the accuracy of the HCI calculation on the active-space orbitals obtained after an HCISCF calculation. In Table II, we show vH-CISCF and HCISCF energies obtained using different 1 cutoffs. We also show the energies obtained by performing a final "tight" HCI calculation on the three different optimized active space orbitals. It is interesting to note that even though the HCISCF energies are clearly un-converged, the final HCI energies agree to within 0.1 mHa in these three cases. This indicates that the optimized active-space orbitals obtained from the three SCF calculations are virtually identical.

Table III shows our gap results for pentacene. The HCI vertical gap for the singlet geometry are slightly higher than that of Kurashige et al. [START_REF] Yuki | Theoretical Study of the * Excited States of Oligoacenes: A Full -Valence DMRG-CASPT2 Study[END_REF] ; however, the difference of 1.5 kcal/mol can be due to the different basis set used, or to small differences in the choice of initial orbitals and subsequent convergence of the DMRG-SCF and HCISCF calculations. Our vertical excitation energy calculated using the triplet geometry agrees reasonably well with the experimental excitation energy of 19.8 kcal/mol reported in Ref. 74 (a direct comparison with the Kurashige et al. paper is not pertinent, as they only report CASPT2 calculations). However, it is worth noting that this energy is significantly different from both the vertical excitation energies calculated using the singlet geometry and the well-to-well excitation energies.

C. Iron (II) Porphyrin

Fe(II) porphyrin (Fe(P)) are the active centers of several important biological proteins such as hemoglobin, myoglobin and catalase. Experimental work suggests that the ground state of Fe(P) is a triplet state belonging to either the 3 A 2g or the 3 E g irreducible representation of the D 4h point group [START_REF] Collman | Stereochemistry, and Structure-Related Properties of alpha beta gamma delta Tetraphenylporphinatoiron(II)[END_REF][START_REF] Goff | Nuclear Magnetic Resonance Investigation of Magnetic and Electronic Properties of "Intermediate Spin" Ferrous Porphyrin Complexes[END_REF][START_REF] Kitagawa | The Resonance Raman Spectra of Intermediate-spin Ferrous Porphyrin[END_REF][START_REF] Dolphin | Synthesis and Mossbauer spectra of octaethylporphyrin ferrous complexes[END_REF][START_REF] Lang | Mössbauer effect study of the magnetic properties of S =1 ferrous tetraphenylporphyrin[END_REF][START_REF] Boyd | Paramagnetic Anisotropy, Average Magnetic Susceptibility, and Electronic Structure of Intermediate-Spin S = 1 (5,10,15,20-Tetraphenylporphyrin)iron(II)[END_REF][START_REF] Mispelter | Proton magnetic resonance characterization of the intermediate ( S =1) spin state of ferrous porphyrins[END_REF][START_REF] Strauss | Comparison of the Molecular and Electronic Structures of (2,3,7,8,12,13,17,18-Octaethylporphyrinato)iron( 11) and 4207 porphyrinato)iron( 11) Steven[END_REF] ; however, a majority of theoretical studies have predicted a quintet 5 A g ground state [START_REF] Manni | Combining the Complete Active Space Self-Consistent Field Method and the Full Con fi guration Interaction Quantum Monte Carlo within a Super-CI Framework , with Application to Challenging Metal-Porphyrins[END_REF][START_REF] Hachmann | The radical character of the acenes: A density matrix renormalization group study[END_REF][START_REF] Hajgató | A benchmark theoretical study of the electronic ground state and of the singlet-triplet split of benzene and linear acenes[END_REF][START_REF] Choe | Theoretical study of the electronic ground state of iron (II) porphine[END_REF][START_REF] Choe | Theoretical study of the electronic ground state of iron II . . . porphine . II[END_REF][START_REF] Pierloot | The CASPT2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes[END_REF][START_REF] Radon | 2 to Heme by Density Functional and Multireference ab Initio Calculations[END_REF][START_REF] Vancoillie | Multiconfigurational second-order perturbation theory restricted active space (RASPT2) studies on mononuclear first-row transition-metal systems[END_REF][START_REF] Phung | Cumulant Approximated Second-Order Perturbation Theory Based on the Density Matrix Renormalization Group for Transition Metal Complexes: A Benchmark Study[END_REF] instead. It is important to note that here and in previous theoretical studies the calculations were performed on a Fe(P) model cluster. Although the model cluster is not identical to the one studied experimentally, their electronic structures are assumed to be very similar. This assumption needs further examination but for this work we will not pursue this further.

The theoretically predicted quintet ground state may be an artifact of the calculation protocol, such as the size of the active space, the method used to calculate the dynamical correlation or the basis set. Here we will explore the effect of the active space and basis set size on the calculated quintet-triplet gap, while we will leave the exploration of the effect of the dynamical correlation method for future work.

The largest CASSCF-like (FCIQMC-CASSCF) calculations performed to date on Fe(P) model cluster have been those by Li Manni et al. [START_REF] Manni | Combining the Complete Active Space Self-Consistent Field Method and the Full Con fi guration Interaction Quantum Monte Carlo within a Super-CI Framework , with Application to Challenging Metal-Porphyrins[END_REF] which fully correlated 32 electrons in a space of 29 orbitals including the 20 C 2 pz , 4 N 2 pz , and all 5 Fe 3 d orbitals. Results from RASSCF and RASPT2 calculations [START_REF] Vancoillie | Multiconfigurational second-order perturbation theory restricted active space (RASPT2) studies on mononuclear first-row transition-metal systems[END_REF] have suggested that using a second set of Fe d orbitals can ease the electron repulsion in occupied d orbitals and lower the triplet state energy. This effect was also observed in a DMRG calculation [START_REF] Olivares-Amaya | The ab-initio density matrix renormalization group in practice[END_REF] where the ground state triplet and quintet calculations were performed using a large active space (44o, 44e) obtained from a 5 A g Hartree-Fock calculation. This large active space included the 29 orbitals of Li Manni in addition to 15 additional orbitals including the 5 Fe 4 d , 1 Fe 4 px , 1 Fe 4 py , 3 N 2 px , and 3 N 2 py orbitals. However, picking these orbitals from the results of a canonical Hartree-Fock calculation is non-trivial. Unlike in the original work, here we perform orbital optimization to minimize the effect of the original orbital choice. We perform HCISCF calculations using two different active spaces, the (29o, 32e) active space of Li Manni and the (44o, 44e) active space used by Olivares-Amaya in the DMRG calculations. Two different HCISCF calculations with (29o, 32e) active space were performed, one with the cc-pVDZ basis set and another with the cc-pVTZ basis set. All calculations were performed using the optimized triplet structure from Ref. 101, which is also the structure used in the DMRG study.

We begin by studying the effect of the accuracy of the HCI calculations in HCISCF on the final optimized active space obtained. Similar to pentacene, three different HCISCF calculations, two of which were vHCISCF with different values of the 1 threshold and the third was a full HCISCF calculation in which the effect of the PT correction was included. The results summarized in Table IV show that the E SHCI energies calculated with different HCISCF orbitals agree with each other to better than 0.2 mHa. This suggests that the active space obtained by the different HCISCF calculations are nearly identical and the difference in HCISCF energies in the second column is largely due to the different HCI tolerance. These results suggest that a cheap variational HCI calculation can be used in the HCISCF procedure with little loss in accuracy of the final orbitals obtained. This observation is in agreement with the pentacene results. TABLE IV. vHCISCF and HCISCF calculations performed on Fe(P) using various thresholds of 1. The PT calculation in HCISCF was performed with 2 = 10 -5 Ha. We also report extrapolated SHCI energies (see text) using the optimized active space orbitals obtained from the SCF calculations. It is interesting to note that (in agreement with pentacene results) although the HCISCF energies are far from converged, the SHCI energies agreement to within 0.1 mHa indicates that the active space orbitals are most likely converged. Table V summarizes the results of our calculations performed with different active spaces and basis set. For each calculation, we show the vHCISCF energy and the SHCI energy extrapolated to the FCI limit calculated using the optimized active space. The extrapolation procedure is identical to the one used for the butadiene and pentacene calculations, with the largest calculation performed with 1 = 10 -5 Ha. The vHCISCF energies are themselves not important because they are quite far from convergence; however, based on the results shown in Table IV we expect the active space to be converged.

With the smaller (29o, 32e) active space we observe that the quintet state is lower in energy than the triplet state by more than 16 kcal/mol. This result remains virtually unchanged as we go from the cc-pVDZ basis set to the cc-pVTZ basis set. However, when the active space is enlarged to (44o, 44e) we see a switching of the energy ordering and find that the triplet is the ground state. It is worth mentioning that identifying the additional 5 Fe 4 d , 1 Fe 4 px , 1 Fe 4 py , 3 N 2 px , and 3 N 2 py orbitals from the canonical Hartree-Fock orbitals is not trivial and so TABLE V. Calculated energies for the Fe(porphyrin) with active spaces of (29o, 32e) and (44o, 44e) and cc-pVDZ and cc-pVTZ basis sets. An initial vHCISCF calculation was performed with 1 = 10 -4 Ha. The optimized active space was used to perform a more accurate SHCI calculation which was extrapolated to near FCI energies with an estimated error shown in the table. The error bar is calculated as 25% of the difference between the extrapolated energy and the most accurate SHCI energy. N d (vHCISCF) is the number of determinants included in the variational space during vHCISCF calculations and N d (SHCI) is the largest number of determinants in the variational space of the SHCI calculation before performing the extrapolation. Eex is the vertical excitation, and TOO and TCI are the wall time in seconds needed for the orbital optimization and HCI calculation steps during a single vHCISCF iteration, using a single node with two 14-core 2.4 GHz Intel Broadwell processors and a combined memory of 128 GB RAM.

Basis

Sym. instead of trying to pick the orbitals by visual inspection we have chosen to include the appropriate number of orbitals (six A g , three B 3u , three B 2u , four B 1g , seven B 1u , eight B 2g , eight B 3g and five A u orbitals) from each irreducible representation in the active space. For both the 5 A g and the 3 B 1g states, a Hartree-Fock calculation was performed to target the 5 A 1g state, after which the initial guess of the active space was chosen. Unlike in Ref. 99, at the first iteration of the HCISCF calculations we observe that the 5 A 1g is still the lower energy state. Further, the CASCI energies of the 5 A 1g and the 3 B 1g states were -2244.94423 and -2244.90676 Hartree when calculated using just the variational HCI with a small 1 = 10 -4 . Although these energies are approximate upper bound of the true energies, they are still lower than the nearly converged DMRG energies reported in Ref. 99, indicating that our initial active-space orbitals are more appropriate. Although at the first iteration the 5 A 1g is lower in energy than the 3 B 1g state, we observe that after HCISCF convergence this ordering is reversed and we obtain the results shown in Table V.

EvHCISCF

Our results strongly suggest that past theoretical results disagreed with experiments because an insufficiently large active space was used, or an inaccurate method for including dynamical correlation was used. The active space suggested by chemical intuition would not include the high lying virtual orbitals such as the 5 Fe 4 d , since it should be possible to capture the energy relaxation due to these orbitals with a dynamical correlation method. Such methods are currently being developed in our group and we plan to use these methods with the smaller active space to see if the correct spin ordering can be obtained.

VI. CONCLUSIONS

The results presented in this work can be used to draw the following three conclusions. First, by using the Lagrangian formulation we can calculate the relaxed reduced density matrices which allow us to straightforwardly use the CASSCF program in PySCF to perform HCISCF. Second, the converged active space orbitals obtained from HCISCF are relatively insensitive to the accuracy of the HCI calculation and consequently loose 1 and 2 thresholds can be used. Third, for large active spaces where getting converged HCI energies become difficult, an initial aHCISCF calculation can be performed to optimizes the active space orbitals while keeping the active space itself fixed. These optimized orbitals can vastly improve the convergence of the HCI calculations to the Full CI or CASCI limit, resulting in great speedups.

Here we have exclusively focused on the development of the HCISCF method as a cheap and accurate approximation to CASSCF. To get quantitatively accurate results it is essential to include dynamical correlation effects by allowing excitations outside of the active space. We have recently worked on developing a particularly accurate method for calculating the dynamical correlation called the multireference linearized coupled cluster theory (MRLCC) [102][103][104] . MRLCC is formulated as a perturbation theory and uses the Fink's partitioning 105,106 of the Hamiltonian. We are currently working on combining the internally-contracted MRLCC with the HCISCF calculation, which will be the focus of a forthcoming paper.

  FIG.1. Outline of the HCI scheme. The dots represent the determinants in a Hilbert space; the connections between them are randomized for the purpose of this sketch. Figure(a) shows that at a given iteration µ the current variational space is V µ (encompassing the cyan dots) and is augmented by the space of connected determinant C µ ( 1) (that includes the blue dots). This sequential aggregation of determinants is done iteratively until some designated convergence threshold is met. Figure(b)shows how HCI takes advantage of its importance function to add new determinants into C µ ( 1). The dotted lines represent possible excitations from every pair of occupied orbitals in a given cyan determinant |Di of a current V µ . A determinant |Da is generated and included in C µ ( 1) if the magnitude of Hia = Di| Ĥ|Da is greater than a certain i-dependent threshold,

TABLE III .

 III Pentacene ground (singlet) state and lowest triplet state energies calculated at the singlet and triplet geometries reported by Kurashige et al.[START_REF] Yuki | Theoretical Study of the * Excited States of Oligoacenes: A Full -Valence DMRG-CASPT2 Study[END_REF] . An initial vHCISCF calculation was performed with an 1 of 8.5×10-5 Ha. The optimized active space obtained after this calculation was used to perform more accurate HCI calculations which were extrapolated to obtain near FCI energies with an estimated error shown in the table. Eex and E ref are respectively the vertical excitation calculated in this work and by Kurashige et al[START_REF] Yuki | Theoretical Study of the * Excited States of Oligoacenes: A Full -Valence DMRG-CASPT2 Study[END_REF] . TOO and TCI, are the wall time in seconds needed for the HCI calculation and orbital optimization step during a single vHCISCF iteration, using a single node with two 14-core 2.4 GHz Intel Broadwell processors and a combined memory of 128 GB RAM.

	Sym. EvHCISCF	ESHCI	Eex	E ref	TOO	TCI
	(Ha)	(Ha)	(kcal/mol) (kcal/mol)	(sec) (sec)
		Singlet Geometry
	1 Ag 3 B2u -841.5457 -841.5936	-841.6174(6) -841.5720(8)	28.5	27.0	50 70	33 24
		Triplet Geometry
	1 Ag 3 B2u -841.5556 -841.5823	-841.6050(7) -841.5751(9)	18.6	-	57 57	26 31

ACKNOWLEDGEMENTS

This work was supported through the startup package of the University of Colorado, Boulder. We would also like to thank Cyrus Umrigar for carefully reading the manuscript and suggesting several improvements.