Cheap and near exact CASSCF with large active spaces
James E. T. Smith Adam A. Holmes Bastien Mussard Sandeep Sharma




james.e.smith@colorado.edu
Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO

Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO

Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO

sandeep.sharma@colorado.edu
Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO






	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.


Cheap and near exact CASSCF with large active spaces
    Sandeep Sharma
    December 30, 2023
====================================================





§ 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 <cit.>, Møller-Plesset <cit.> and coupled cluster theory <cit.>.
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 multi-determinant 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 <cit.> (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 <cit.>.
Naturally, there is great interest in developing algorithms that can overcome this limit in a systematically improvable way.
There are many algorithms that can be used as approximate active-space solvers,
including the density matrix renormalization group<cit.>
restricted <cit.> and generalized <cit.> active space methods, reduced density matrix approaches <cit.>,
various flavors of selected configuration interaction followed by perturbation theory<cit.> (SCI+PT), and stochastic and semistochastic methods such as Full CI quantum Monte Carlo <cit.> (FCIQMC).



Here, we use the recently-developed Heat-bath Configuration Interaction <cit.> (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<cit.>, FCIQMC-SCF<cit.>.
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. ∂ E_HCI/∂c̱≠ 0. Thus, the first-order change in energy due to change in orbital rotation parameters κ is given by


    d E_HCI/d κ = ∂ E_HCI/∂κ + ∂ E_HCI/∂c̱dc̱/d κ,

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 <ref> suggests 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.
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).<cit.> 


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 <ref> is 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 (^1A_g) and lowest triplet (^3B_2u) states, and compare the gaps to experimental<cit.> and theoretical<cit.> 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.
<cit.> 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 ^5A_g and ^3B_2g states.

The remainder of the paper is organized as follows.
In Section <ref>, 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 <ref>, we derive the working equations of the HCISCF and vHCISCF algorithms. In Section <ref>, we discuss the implementation and practical aspects of running an HCISCF calculation.
Finally, in Section <ref>, we report our HCISCF calculations for butadiene, pentacene, and Fe(P), and compare them to experiments and previous calculations on these systems.



§ 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.







 §.§ Variational Stage


In the variational stage, a set of determinants 𝒱 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 𝒱^μ and the current variational wavefunction is given by
eq:Psi0
    |Ψ_0^μ⟩=	∑_|D_i⟩∈𝒱^μ c_i^μ|D_i⟩,

where c̱^μ are the current CI coefficients.

The variational space is augmented (see Figure  <ref>) as 𝒱^μ+1=𝒱^μ+𝒞^μ(ϵ_1) where ϵ_1 is a user-defined parameter. The space of connected determinants 𝒞^μ (ϵ_1) is defined as

    𝒞^μ (ϵ_1)={|D_a⟩ |  f^μ(|D_a⟩)> ϵ_1 },

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_|D_i⟩∈𝒱^μ|H_ai c^μ_i|,

where H_ai is the Hamiltonian matrix element between |D_a⟩ and |D_i⟩.
In other words, the space 𝒞^μ(ϵ_1) is composed of all determinants |D_a⟩ for which |H_aic^μ_i|>ϵ_1 is true for at least one determinant |D_i⟩ in the current 𝒱^μ.

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
eq:CIPSI
	f_CIPSI^μ(|D_a⟩) =  | ∑_|D_i⟩∈𝒱^μH_aic^μ_i /E_0^μ-E_a|,

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 <ref>). Second, the HCI importance function 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 𝒞(ϵ_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̱^μ, ℰ_0^μ; 𝒱^μ]=
           Ψ_0^μĤ_0 Ψ_0^μ
        - ℰ_0^μ (⟨Ψ_0^μ|Ψ_0^μ⟩-1)
    	,

with respect to the coefficients c̱^μ of |Ψ_0^μ⟩ and to the Lagrange multiplier ℰ_0^μ ensuring that the new wavefunction remains normalized,
in the context of a space 𝒱^μ which is fixed at this point.




 §.§ Perturbative stage

Once a converged
variational space 𝒱,
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 = ∑_|D_i⟩, |D_j⟩∈𝒱 H_ij|D_i⟩⟨D_j| + ∑_|D_a⟩∉𝒱 H_aa|D_a⟩⟨D_a|
    ,

and perturbation V̂ = Ĥ - Ĥ_0. With this partitioning, the second-order perturbative energy correction to the variational energy is given by

    E_2=
    	∑_|D_a⟩∈𝒞1/E_0-H_aa(∑_|D_i⟩∈𝒱 H_aic_i)^2,

where 𝒞 denotes the set of determinants that are connected to at least one determinant in 𝒱 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)=
        ∑_|D_a⟩∈𝒞(ϵ_2)1/E_0-H_aa(∑_|D_i⟩∈𝒱^(ϵ_2) H_aic_i)^2,

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
𝒞(ϵ_2) that are connected to at least one determinant in 𝒱 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. <ref> requires one to simultaneously store the entire set of determinants 𝒞(ϵ_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 𝒞(ϵ_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.  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 ϵ_2^d 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_2^D(ϵ_2^d) + [ E_2^S(ϵ_2) -  E_2^S(ϵ_2^d)].

The key point is that E_2^S(ϵ_2) and E_2^S(ϵ_2^d) 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 <ref>) 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 <ref>, the second-order correction to the energy can equivalently be obtained by minimizing the Hylleraas functional

    H[ḏ; c̱, 𝒱, 𝒞(ϵ_2)]
          = Ψ_1^ϵ_2Ĥ_0-E_0 Ψ_1^ϵ_2
       +2Ψ_1^ϵ_2V̂Ψ_0
    ,

with respect to the coefficients ḏ of the first-order wavefunction

    |Ψ_1^ϵ_2⟩=∑_a ∈𝒞(ϵ_2) d_a |D_a⟩
    ,

with c̱, 𝒱 and 𝒞(ϵ_2) held fixed.

At its minimum, the Hylleraas functional gives the optimal values of ḏ, and the value of the functional is equal to the second-order correction evaluated using Eq <ref>. This will be used in Section <ref> 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 𝒞(ϵ_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 𝒞(ϵ_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 𝒱), the inequality H[ḏ; c̱, 𝒱, 𝒞(ϵ_2)] ≥ E_2(ϵ_2) holds, where the equality holds at its minimum.




§ 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 second-order 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.



 §.§ 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 <ref> with respect to c̱ and ℰ_0^μ, and the Hylleraas functional in Eq <ref> with respect to c̱, ℰ_0, and ḏ.
Formally, the minimization of these functionals is performed by setting to zero their partial
derivatives with respect to the parameters c̱, ℰ_0 and ḏ.
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̱, ℰ_0, ḏ]=E[κ, c̱, ℰ_0] + H[κ, c̱, ḏ]
    ,

where the dependence on the parameters κ is shown and the dependence on 𝒱 and 𝒞 is dropped.
Its derivative with respect to κ reads:

    d E_HCI/d κ    = ∂ E_HCI/∂κ
    	+ ∂ E_HCI/∂c̱dc̱/d κ
         
    	+ ∂ E_HCI/∂ℰ_0dℰ_0/d κ
    	+ ∂ E_HCI/∂ḏdḏ/d κ.

As mentioned in the introduction, this formulation requires the calculation of the derivatives of the HCI parameters with respect to κ, which makes it an impractical approach.

Instead, we introduce the Lagrangian

	
    ℒ[κ, c̱, ℰ_0, ḏ, λ_c, λ_d]
    		= E[κ, c̱, ℰ_0]
    		+ H[κ, c̱, ḏ]
    		+ λ_c^†∂ E/∂c̱^†
    		+ λ_d^†∂ H/∂ḏ^†,


which is a function of the variables κ, c̱, and ḏ, and the set of Lagrange multipliers ℰ_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 ḏ can be recovered.
Note that at these optimal values of c̱ and ḏ, 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 ḏ. 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


    d E_HCI/d κ   = d ℒ/d κ = ∂ℒ/∂κ
    ,

leading to the desired minimization of the HCI energy with respect to κ under the constraint that ∂ E/∂c̱=0 and ∂ H/∂ḏ=0.

Let us first derive the equations to obtain the Lagrange multipliers:

    ∂ℒ/∂ḏ    = ∂ H/∂ḏ + λ_d^†∂^2 H/∂ḏ∂ḏ^†    =0 
    ∂ℒ/∂c̱    = ∂ E/∂c̱ + ∂ H/∂c̱ + λ_c^†∂^2 E/∂c̱∂c̱^†  + λ_d∂^2 H/∂c̱∂ḏ^†    =0

Here, we have assumed that all parameters are real numbers, although extension to complex numbers is straightforward without any additional complications. From Equation <ref>, we can infer that λ_d=0 because  ∂ H/∂ḏ=0. Equation <ref> simplifies to

    ∂ H/∂c̱ + λ_c^†∂^2 E/∂c̱∂c̱^†    =0
    ,

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


    ∂ℒ/∂κ   = c̱^†∂ H_0/∂κc̱
        +  ḏ^†∂ H_0/∂κḏ
        + 2ḏ^†∂ V/∂κc̱
        + 2λ_c^†∂ H_0/∂κc̱
       = ∑_ijkl∂ H_0,ijkl/∂κΓ_ijkl^c,c
        +  ∂ H_0,ijkl/∂κΓ_ijkl^d,d
        + 2∂ V_ijkl/∂κΓ_ijkl^d,c
        + 2∂ H_0,ijkl/∂κΓ_ijkl^λ_c,c,


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 integrals with respect to κ are calculated using the usual techniques, which are described in detail in Ref..






 §.§ vHCISCF


When only the variational energy E_ HCI is optimized instead of the total HCI energy, one can calculate the gradient with respect to κ as

    d E_HCI/d κ   =  ∂ E/∂κ + ∂ E/∂c̱dc̱/d κ
       =∂ E/∂κ
       = ∑_ijkl∂ H_0,ijkl/∂κΓ_ijkl^c,c.

Thus, to calculate the energy gradient, one does not need to evaluate the Lagrange multipliers, and the simple variational two-body reduced density matrix Γ_ijkl^c,c is sufficient.




§ COMPUTATIONAL DETAILS



We briefly describe the calculation of the reduced density matrices encountered in Eq. (<ref>). The density matrices Γ_ijkl^c,c and Γ_ijkl^λ_c,c are strictly limited to states that contain determinants from the variational space 𝒱. 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 Γ_ijkl^d,d 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 𝒞(ϵ_2), the reduced density matrix
simply reads:


    Γ_ijkl^d,d = ∑_|D_a⟩∈𝒞(ϵ_2) d_a^2 D_aa_i^† a_j^† a_k a_lD_a.

Finally, for the evaluation of Γ_ijkl^d,c, all the determinants |D_i⟩ in 𝒱 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 Γ_ijkl^d,c as

    Γ_ijkl^d,c = ∑_|D_i⟩∈𝒱
    |D_a⟩∈𝒞(ϵ_2)
    |H_ia| > ϵ_2 c_i d_aD_ia_i^† a_j^† a_k a_lD_a.



Thus out of all the different reduced density matrices only the evaluation of Γ_ijkl^d,c 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 𝒞(ϵ_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 vHCISCF 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 <ref> 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.






§ 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.



 §.§ Butadiene

The HCI calculations on butadiene were performed with the same ANO-L-pVDZ basis set and geometry as the one used in Ref. . 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 <ref> 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 <cit.>, we extrapolate the HCI energy to the Full CI limit by extrapolating to where the perturbative correction goes to zero. In this work, we do this by first performing a single HCI calculation in which the smallest ϵ_1 allowed by the available computational resources is used, which for butadiene is approximately 10^-5 Ha. Some subsequent single step HCI calculations are performed in which a fraction of the least important determinants are discarded (in this work: 0.33). After a few such iterations, we plot the total HCI energy versus the PT correction and perform a linear extrapolation to E_2→ 0, as shown in Figure <ref>. 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.









 §.§ 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.<cit.>
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 ^1A_g ground state and ^3B_2u lowest triplet state<cit.> in the cc-pVDZ basis set on the optimized singlet and triplet geometries.
This allows us to calculate both vertical and relaxed (well-to-well) ^1A_g→^3B_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 <ref>, we show vHCISCF 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 <ref> shows our gap results for pentacene. The HCI vertical gap for the singlet geometry are slightly higher than that of Kurashige <cit.>; 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. 
(a direct comparison with the Kurashige 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.






 §.§ 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 ^3A_2g or the ^3E_g irreducible representation of the D_4h point group <cit.>; however, a majority of theoretical studies have predicted a quintet ^5A_g ground state<cit.> 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. <cit.> 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<cit.> 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<cit.> where the ground state triplet and quintet calculations were performed using a large active space (44o, 44e) obtained from a ^5A_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. , 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 <ref>, show that the active space orbitals converge much more quickly than the HCISCF energies do. This observation is in agreement with the pentacene results.





Table <ref> summarizes the results of our calculations performed with different active spaces and basis set. For each calculation, we show the vHCISCF energy and the converged HCI energy extrapolated to the FCI limit. 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 <ref> 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 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 ^5A_g and the ^3B_1g states, a Hartree-Fock calculation was performed to target the ^5A_1g state, after which the initial guess of the active space was chosen. Unlike in Ref. , at the first iteration of the HCISCF calculations we observe that the ^5A_1g is still the lower energy state. Further, the CASCI energies of the ^5A_1g and the ^3B_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. , indicating that our initial active-space orbitals are more appropriate. Although at the first iteration the  ^5A_1g  is lower in energy than the  ^3B_1g state, we observe that after HCISCF convergence this ordering is reversed and we obtain the results shown in Table <ref>.

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.



§ 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)<cit.>. MRLCC is formulated as a perturbation theory and uses the Fink's partitioning<cit.> 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.





§ 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.


biochem



