Combining the Complete Active Space Self-Consistent Field Method

Article pubs.acs.org/JCTC

Combining the Complete Active Space Self-Consistent Field Method and the Full Configuration Interaction Quantum Monte Carlo within a Super-CI Framework, with Application to Challenging MetalPorphyrins Giovanni Li Manni,* Simon D. Smart, and Ali Alavi* Max-Planck-Institut für Festkörperforschung, Heisenbergstraße 1, 70569 Stuttgart, Germany S Supporting Information *

ABSTRACT: A novel stochastic Complete Active Space Self-Consistent Field (CASSCF) method has been developed and implemented in the Molcas software package. A two-step procedure is used, in which the CAS configuration interaction secular equations are solved stochastically with the Full Configuration Interaction Quantum Monte Carlo (FCIQMC) approach, while orbital rotations are performed using an approximated form of the Super-CI method. This new method does not suffer from the strong combinatorial limitations of standard MCSCF implementations using direct schemes and can handle active spaces well in excess of those accessible to traditional CASSCF approaches. The density matrix formulation of the Super-CI method makes this step independent of the size of the CI expansion, depending exclusively on one- and two-body density matrices with indices restricted to the relatively small number of active orbitals. No sigma vectors need to be stored in memory for the FCIQMC eigensolvera substantial gain in comparison to implementations using the Davidson method, which require three or more vectors of the size of the CI expansion. Further, no orbital Hessian is computed, circumventing limitations on basis set expansions. Like the parent FCIQMC method, the present technique is scalable on massively parallel architectures. We present in this report the method and its application to the free-base porphyrin, Mg(II) porphyrin, and Fe(II) porphyrin. In the present study, active spaces up to 32 electrons and 29 orbitals in orbital expansions containing up to 916 contracted functions are treated with modest computational resources. Results are quite promising even without accounting for the correlation outside the active space. The systems here presented clearly demonstrate that large CASSCF calculations are possible via FCIQMC-CASSCF without limitations on basis set size.

1. INTRODUCTION Most chemical systems of practical interest feature electronic structures dominated by more than one electronic configuration. These are often referred to as strongly correlated systems. Transition metal complexes, playing crucial roles in the active sites of natural enzymes and being core components of human-made catalysts, represent an important subset of such entangled systems.1,2 Near degeneracy effects are generally invoked to explain the complex nature of such systems. A multitude of valence orbitals sharing the same position in the energy spectrum poses a real challenge for ab initio methods based on a single-configurational representation. Among single-configurational methods, the Hartree−Fock approximation3,4 (in the form of the Roothaan equations5,6) and density functional theory (in its Kohn−Sham formulation,7 KS-DFT) are widely used to generate molecular orbitals and provide qualitative, and in some cases quantitative, means to rationalize chemical processes which are otherwise hard to explain. In general, these methods are quite accurate for closed shell systems and to date are the only methods that can be routinely applied to large systems. However, being single reference methods by construction, they are incapable of © 2016 American Chemical Society

representing multiconfigurational systems where orbitals and electrons undergo significant rearrangements, such as in bond formation or weakening. Other methods attempt to improve on Hartree−Fock solutions by using the HF wave function (or other single determinantal wave functions) as reference for an a posteriori correction, such as in the coupled-cluster (CC) approach. These approaches, relying on the exactness of the reference wave function, would identically fail. In fact, for any single determinantal wave function (say the Hartree−Fock wave function), the orbitals are optimized under the field of a single electronic configuration, making the method a poor choice as an orbital generator for electronically entangled systems. Multiconfigurational self-consistent field (MCSCF) theories have been developed since the early years of quantum chemistry as a natural extension of the HF method for strongly correlated systems. In these approaches, orbitals are variationally optimized simultaneously with the configuration interaction expansion coefficients. As orbitals are optimized under the field Received: December 15, 2015 Published: January 25, 2016 1245

DOI: 10.1021/acs.jctc.5b01190 J. Chem. Theory Comput. 2016, 12, 1245−1258

Article

Journal of Chemical Theory and Computation

the SplitGAS method,17 in which active subspaces of GAS type define two different CI expansions, a principal space (P), containing ∼109 SDs, and an extended space (Q), containing a much larger expansion. By means of Löwdin’s partitioning technique, the size of the secular problem in (P+Q) is reduced to the dimensionality of the principal space. With SplitGAS, active spaces containing 24 electrons and 48 active orbitals have been reported for the Cr2 ground state potential energy curve. The success of the above methods relies on the truncation of the complete active space according to some chemical knowledge of the system under analysis, as documented in the literature. There exist cases where this truncation is not straightforward, and wrong choices may lead to erroneous results.18 Correlation effects not included at the MCSCF level of theory are usually recovered a posteriori, such as by using a second-order perturbation theory (CASPT2 or RASPT2)19−21 with optimized CASSCF or RASSCF wave functions as the reference. The multiconfiguration pair-density functional theory (MC-PDFT)22−25 represents a valid and cheap alternative to CASPT2, but still relies on the qualitative correctness of the reference wave function, of CAS or GAS type, and its limit lies once again on the preceding multiconfigurational step. It is important also to mention the success of the Density Matrix Renormalization Group (DMRG) approach in tackling strong correlation effects, including some DMRG-SCF implementations.26−33 In this report, we present a novel CASSCF approach that retains the flexibility and simplicity of CAS wave functions, while circumventing the exponential scaling wall encountered by standard CASSCF implementations. The standard Davidson CI eigensolver is replaced by the stochastic Full Configuration Interaction Quantum Monte Carlo (FCIQMC) approach,34−42 and it is named FCIQMC-CASSCF. Analogous to other twostep MCSCF implementations, a macro-iteration in the FCIQMC-CASSCF calculation can be divided into two parts. In the first step, given the current set of molecular orbitals, the CI-coefficients are optimized using FCIQMC, and the one- and two-body reduced density matrices are accumulated. Following this, in the second stage, orbital optimization is performed using an approximated form of the Super-CI method.8−11 To the best of our knowledge, only one other MCSCF scheme43 using FCIQMC as its CI eigensolver has been developed to date. The two methods differ in the orbital optimization step with the method developed by Thomas et al. using a Newton− Raphson approach44,45 rather than the approximated Super-CI scheme reported here. The method presented here has been developed in a locally modified version of Molcas 8.46 We have tested the FCIQMC-CASSCF method on the freebase porphyrin, H2P (also known as porphin in the old literature), and two metalloporphyrins, MgP and FeP, containing the Mg(II) and Fe(II) ions, respectively, at the center of the aromatic macrocycle. For H2P and MgP, the Q absorption bands have been investigated, while for the FeP compound, six low-lying electronic states have been investigated aiming to identify the ground state. We turned our attention to this class of compounds for several reasons. Natural porphyrins (suitably modified) provide the key material for photosynthesis (chlorophyll), respiration (hemoglobin), and electron-transport (cytochromes) and therefore play crucial roles in life. Synthetic porphyrins are also studied for their potential applications in solar cells,47,48 as catalysts in the water splitting process49−51 and even as molecular switches.52 Given their fundamental role, electronic spectra of porphyrins have

generated by multiconfigurational many-electron wave functions, they are not optimal for any specific configuration (the HF determinant for instance). Rather, they represent the best orbitals for the set of configurations considered. Multiconfigurational methods did not gain popularity until the complete active space (CAS) wave function8−11 was developed, as there are practical difficulties in “picking” the right configurations in earlier methods. In the CAS approach, these difficulties have been partially mitigated by the introduction of the concept of active space. The orbital expansion is split into three parts: inactive (also referred to as core orbitals), active, and virtual orbitals. The inactive and virtual orbitals are doubly occupied and empty, respectively, in all configurations of the MC expansion. By contrast the CI expansion is complete across all possible (spatial- and spin-symmetry allowed) distributions of the remaining electrons in the active orbitals. Therefore, the expectation value of the occupation numbers of the corresponding natural orbitals may vary continuously between 0 and 2. The CAS wave function has several practical (on the user side) and computational (on the developer side) advantages. Most importantly, the problem of choosing “dominant configurations” is changed into one of choosing “active orbitals”, which is more easily accessible to the chemistry community. Of course, these advantages come with a cost; the size of the CI expansion grows exponentially with the number of electrons and orbitals in the active space. This exponential scaling has represented the main limitation of the method for many years, limiting it to being able to correlate at most 18 electrons and 18 active orbitals. A CAS(18,18) for a singlet spin state and without space symmetry constraints contains ∼2 × 109 Slater determinants. Multiple arrays, each of ∼20 GB, must be allocated in memory to accommodate expansion coefficients, sigma vectors, update vectors, and other data related to the ordering of the configurations. Specifically, Molcas would require about 100 GB of memory to perform such a calculation. A CAS(20,20) for a singlet spin state would already contain ∼34 × 109 Slater determinants, corresponding to ∼270 GB per array. This type of calculation would require computers with ∼2 TB of memory, which are not accessible on single node computers to our knowledge. Developments in shared/ distributed memory would probably allow us to exceed this limit in the near future. Working in the space of configurational state functions (CSFs) reduces the size of the CI vectors (this can be efficiently achieved using Shavitt’s Graphical Unitary Group Approach, GUGA),12,13 but the memory requirements remain high, with a CAS(20,20) requiring about 50 GB per array. There have been several attempts to tackle this “exponential wall” in MCSCF calculations. In the Restricted Active Space (RAS) approach,14,15 for instance, the active space is divided in three subsets: RAS1, RAS2, and RAS3. Restrictions on the electron occupations are introduced for the RAS1 and RAS3 subspaces, such that in RAS1 a maximum number of holes are allowed and in RAS3 a maximum number of electrons are allowed. The remaining active electrons are distributed in all possible ways among the RAS2 orbitals. The Generalized Active Space (GAS)16 approach is a generalization of the RAS concept. Rather than three spaces, an arbitrary number of active subspaces may be chosen. For each GAS space, cumulative minimum and maximum occupation numbers are chosen. This has enabled the use of larger active spaces, for instance, the GAS-5(20,32) for the Gd dimer. It is also worth mentioning 1246


Article


molecular orbital rotation parameters. The expansion coefficients Ci0 of eq 1 may be introduced as a matrix exponential of an anti-Hermitial matrix

been extensively studied, both experimentally and theoretically. Because of the size of these molecules, most computational investigations have been carried out by means of density functional theory (DFT),53−55 although there are always concerns about the reliability of DFT in handling multireference systems such as the iron porphyrins in which multiple spin states lie very close together in energy. The free-base porphyrin, with 26 valence electrons distributed in 24 near degenerate π valence orbitals, already poses challenges for ab initio simulations. Recently some correlated ab initio methods have been employed, including coupled cluser (CC),56,57 Diffusion Monte Carlo (DMC),58 second-order perturbation theory (CASPT2),59,60 symmetry-adapted cluster-configuration interaction (SAC−CI),61,62 multireference configuration interaction (MRCI),63,64 and even some hybrid approaches, such as the DFT/MRCI of Grimme.65 This list could easily be extended. These calculations have not been made without a number of assumptions to simplify the computational costs. Our contribution in this field is important as it represents the first fully variational calculation on this class of macrocycles which correlates the entire (π − π*) system (plus the iron valence orbitals and their electrons for the ferrous porphyrin) without further assumptions, simplifications, or truncation of the CI expansion. The manuscript is organized as follows: In Section 2, we describe the formulation of the method and highlight some relevant algorithmic details. In Section 3, we discuss the computational details for the systems investigated, with a focus on basis set, geometry, and active space. In Section 4, our results are analyzed and discussed. Concluding remarks are offered in Section 5.

Ci0 = [exp( −S)]i0

or equivalently in an operator representation, e , with Ŝ defined as Ŝ =

T̂ =

(6)

i>j

with the excitation operators defined as Epq = ap†αaqα + ap†β aqβ

and (apσ) is the usual creation (annihilation) operator for a molecular orbital p with spin σ. Variations of the MCSCF states can thus be written as ̂

̂

0′⟩ = eT e S 0⟩

pq

̂ ̂ ̂ ̂ E(T, S) = ⟨0 e−Se−T Hê T e S 0⟩

1 ⟨0 [[Ĥ , 2 1 + ⟨0 [[Ĥ , 2

1 2

∑ gpqrsdpqrs pqrs

⎛ ∂ 2E ⎞ ⎟⎟ p = 0 n ⎝ ∂pm ∂pn ⎠

T̂ ], S]̂ 0⟩ + ... (9)

(2)

Hessian and gradient components, respectively. Given the extraordinary size of the Hessian matrix, this equation is never solved exactly in real calculations. In the unfolded two-step (2) Newton−Raphson method, the coupling terms E(2) oc and Eco of the Hessian are neglected, giving the uncoupled equations

∑ ⎜⎜

0

1 ̂ S]̂ 0⟩ T̂ ], T̂ ] 0⟩ + ⟨0 [[Ĥ , S], 2

The first term gives the energy, and the next two terms represent the gradient with respect to orbitals and CI coefficients, respectively. The following terms are Hessian terms, which can be divided into three components: orbital− orbital part (oo), configuration−configuration part (cc), and orbital−configuration parts (oc and co). By separating orbital components from CI componentswhere E(2) and E(1) are the

where indices p, q, r, and s refer to generic spatial orbitals. The information about the MOs is encoded entirely within the oneand two-electron integrals, h and g, and the wave function in the one- and two-body reduced density matrices, D and d. In MCSCF, this energy expression is minimized with respect to both CI coefficients and molecular orbitals, which represent the variational parameters, pm. Truncating a Taylor expansion of the energy to second order in terms of the variational parameters and setting its derivative to zero leads to a system of linear Newton equations

m,n

(8)

Substituting eq 8 into eq 2, a Baker−Campbell−Hausdorf f (BCH) expansion to second order is obtained for the MCSCF energy

expanded in an orthonormal determinantal basis of size m, the energy may be obtained as the expectation value of the spinfree, non-relativistic Hamiltonian

∑ hpqDpq +

(7)

a†pσ

(1)

i

m

∑ Tij(Eiĵ − Ejî )

= ⟨0 Ĥ 0⟩ + ⟨0 [Ĥ , T̂ ] 0⟩ + ⟨0 [Ĥ , S]̂ 0⟩

∑ Ci0 i⟩

⎛ ∂E ⎞ ⎟⎟ + 1 ∂ p 2 ⎝ m ⎠0

(5)

+

E = ⟨0 Ĥ 0⟩ = hnuc +

− 0⟩⟨K )

States |K⟩ are the orthogonal complement states to the target state |0⟩. The connection between the parameters SK0 and the expansion coefficients Ci0 becomes obvious. Similarly, an exponential form of the molecular orbital unitary tranŝ formation can be introduced, eT, with T̂ defined by

m

∑ ⎜⎜

∑ SK 0( K ⟩⟨0 K ≠0

2. METHOD Given a normalized CI wave function 0⟩ =

(4) −Ŝ

⎧ E(2)S = −E(1) ⎪ cc c ⎨ ⎪ E(2)T = −E(1) ⎩ oo o

(3)

The variational parameters, pm, are introduced in the form of unitary transformations of the wave function within the space spanned by both the configuration expansion coefficients (also referred to as the CI vector or the CI coefficients) and the

(11)

The first of these equations is promptly recognized as the CI secular problem. In FCIQMC-CASSCF, this step is performed 1247


Article

Journal of Chemical Theory and Computation using FCIQMC as discussed in Section 2.2. The orbital rotation parameters are obtained by solving the second equation. E(2) oo can be written in terms of one- and two-electron integrals in a MO basis and one- and two-body reduced density matrices. Without going into the details, it is important to mention that the orbital Hessian elements contain integrals of the type (ab| pq) and (ap|bq), where indices a and b run over all orbitals and p and q over the inactive and active orbitals. Ignoring spacial symmetry, the size of this matrix is (nin)2, where ni is the number of active and inactive orbitals, and n is the total number of orbitals. For example, if a calculation is performed with ni = 90 and n = 900 (these numbers reflect the size of the calculations reported here), ∼50 GB of memory would be required to store the full Hessian. The operation count for the transformation from AO to MO basis is proportional to (nin)4 and may also be large. 2.1. Super-CI for Orbital Rotation. The Super-CI approach was first developed by Grein66 and Ruedenberg.67 The optimal orbitals are found by solving the Super-CI secular equations, defined by the wave function at the point of expansion, |0⟩, and all possible singly excited states ̂ − Eqp ̂ ) 0⟩ p → q⟩ = (Epq

is adopted for the matrix elements coupling singly excited states, ⟨p → q |Ĥ |r → s⟩, where the indices p, q, r, s, and t run over the entire orbital space. Simplifications arise from the partitioning of the orbital expansion into different subspaces (inactive, active, and virtual space). No approximations are made for the terms coupling the reference wave function |0⟩ with itself or with the singly excited states, ⟨0|H|p → q⟩, as no third-order density matrices are required for those elements. The motivation behind the one-particle model Hamiltonian lies in the fact that in MCSCF, although one searches for the point in energy where the gradient vanishes, the precise path taken to that point is less crucial. Equation 15 shares some similarities with the Fock matrix in HF theory. Due to the invariance with respect to inactive−inactive, active−active, and virtual−virtual rotations, only interspace single excitations must be addressed. One important computational aspect arising from the choice of the one-electron model Hamiltonian operator is that most of the terms are easily evaluated directly from a list of atomic, one, and two electron integrals (similarly to the closed-shell HF method). A small additional effort is required due to the evaluation of molecular two-electron integrals of the type (qv| xy), where q runs over the entire orbital space, while v, x, and y are limited to the active space. The resultant array is of size (n· n3a) (na is the number of active orbitals), in comparison to the much larger Hessian matrix, (nin)2, required for the Newton− Raphson approach. One- and two-reduced density matrices are also required, but as their indices run only over the active space, they do not represent a limiting factor. The operation count for the AO to MO transformation is reduced with respect to the Newton approach, being of the order of (nan)4. Although the asymptotic convergence of the Super-CI approach is only first order and theoretically slows down in proximity to the stationary point, in Molcas, the first-order convergence has been enhanced using a quasi-Newton update approach.15 2.2. FCIQMC as CI Eigensolver. In FCIQMC-CASSCF, the FCIQMC method developed by Alavi et al.34−36 is used to solve the CI secular problem. FCIQMC is a projector technique that propagates the CI vector in imaginary time using a stochastic algorithm, which in the long-time limit leads to an accurate sampling of the ground-state CI wave function. The starting point for this method is the imaginary-time Schrödinger equation

(12)

Here, the effect of the exitation operators (Ê pq − Ê qp) on each configuration of |0⟩ is that of creating two new configurations in which the spin−orbital of index q is replaced by p and vice versa. In the particular case of CAS wave functions, p and q are orbitals belonging to different subspaces, as rotation between two orbitals in the same subspace are redundant. Therefore, only inactive−active, active−virtual, and inactive−virtual rotations are considered when generating the |p → q⟩ states. The Super-CI wave function, |SCI⟩, is therefore defined as SCI⟩ = 0⟩ +

∑ χpq p → q⟩ (13)

p>q

The method involves solving the corresponding secular equations to obtain the Super-CI linear coefficients, χpq. These coefficients are subsequently used in matrix exponential form to define a unitary transformation of the orbitals appearing in |0⟩ (or some approximation to it68). Far from the stationary points, the wave function is improved by adding all single excitation contributions to |0⟩. Iteration after iteration, |0⟩ will approach the wave function at the stationary point. Coefficients χpq will decrease as the SCF optimization proceeds, and at convergence, all coefficients will vanish, implying that no more rotations are needed and |0⟩ will reveal the variationally optimized wave function. If solved exactly, the Super-CI secular problem would be more costly than the Newton−Raphson method as it implies the evaluation of third-order density matrix elements, Qvxyztu, with indices in the active space.11 To avoid this, an approximated Super-CI approach is employed in the method presented in this paper, in which an effective one-electron Hamiltonian Ĥ eff =

∑ f pq Epq̂

∂Ψ = −(Ĥ − ES)Ψ ∂τ

where ES is an in principle arbitrary energy offset (referred to as the “shift”). Starting from an initial wave function, Ψ(0), the long-time limit becomes lim Ψ(τ ) → Ψ0⟩⟨Ψ0 Ψ(0)⟩e−τ(E0 − ES)

τ →∞

with the Fock matrix elements f pq =

∑ Dpr hrq + ∑ dprstgqrst r

r ,s,t

(17)

where |Ψ0⟩ and E0 are the normalized ground-state wave function and energy, respectively, of the system. If ES = E0, then we obtain a solution Ψ(τ) that is proportional to the exact ground state |Ψ0⟩, with a finite and nonvanishing norm. In the context of FCIQMC, ES plays the role of population control and also serves a useful measure of the total energy of the system, independent of the projected or variational measures to be discussed later on. The initial wave function can be any wave function with a non-zero overlay with |Ψ0⟩, and in practice, the Hartree−Fock wave function is often used. In some cases, multideterminantal wave functions may be used as reference wave functions to improve the efficiency of the

(14)

p,q

(16)

(15) 1248


Article

Journal of Chemical Theory and Computation method.38,69 Similarly, by choosing initial wave functions in

determinant becomes equal to the FCI coefficient of that determinant:

different symmetry sectors, the appropriate ground-states

Ci ∝ ⟨Ni⟩τ

solutions in those symmetry sectors can be obtained.

This is achieved via a sequence of stochastic steps which update the walker distribution each iteration such that a population dynamics simulates eq 18. Although the wave function is never explicitly represented, it is stochastically sampled such that accurate energy and other estimates are obtained. The method has some points in common with Dif f usion Monte Carlo (DMC).71,72 The FCIQMC dynamics is divided into three parts. A spawning step in which child particles are created from their parents into various locations in the Hilbert space, a death step in which parent walkers are stochastically removed from the simulation, and an annihilation step in which surviving parent walkers and newly spawned walkers are removed from the simulation if they are present on the same determinant in the Hilbert space but with opposite signs. In the spawning event, particles are stochastically created from the parent walkers with weight proportional to |Hij|, while in the death step, parent particles are removed with magnitude proportional to Hii − ES. A simplified algorithm for each FCIQMC iteration is offered in Figure 1.

Short-time integration of eq 16 with a time step, Δτ, provides an iterable working equationIn terms of the amplitudes Ci(τ), eq 18 becomes Ci(τ + Δτ ) = Ci(τ ) − Δτ(Hii − ES)Ci(τ ) − Δτ ∑ HijCj(τ ) (19)

j≠i

A nonstochastic scheme for the above equation would require the storage of several arrays of size equal to the full set of Ci(τ) coefficients. This rapidly becomes impractical and faces the same exponential scaling problem as direct CI schemes. Instead, the concept of walkers (or particles) is introduced in FCIQMC. A walker is a Kronecker delta function that resides on a specific determinant; thus, if the α-th walker resides on |iα⟩, the corresponding Kronecker delta is δi,iα. Each walker also has a sign, sα = ± 1, associated with it. Therefore, given a set of Nw walkers, {i1, i2, ..., iNw}, the signed number of walkers, Ni, which reside on |i⟩ is given by Ni =

∑ sαδi ,i

α

(21)

(20)

α 39

In practice, a hash algorithm is used to communicate and address all walkers that reside on a particular determinant to the same processor and location in memory (i.e., entry in a hash table). The above sum can therefore be restricted to the elements of that specific entry of the hash table (which is typically only a few walkers, much smaller than the entire set of walkers), and Ni is stored for each occupied determinant. Furthermore, whenever a new walker is created (via the spawning step of the FCIQMC algorithm, as explained below), it is communicated to its designated processor and entry in the hash table. If the determinant of the new walker already exists in the hash table (i.e., it is currently occupied), the new walker adds or subtracts to the population on that determinant, thereby achieving the annihilation step of the algorithm. Otherwise, a new entry in the hash table is created for the new walker. Thus, the key point of the FCIQMC algorithm is that instead of allocating memory for accommodating the entire CI vector, only determinants populated by walkers are stored. The socalled “deadwood” (negligible configurations) are either rarely populated or not populated at all during the dynamics and therefore do not contribute significantly to the memory requirements or to the computational cost as unoccupied determinants need not be enumerated in an iteration. As a technical detail, the implementation of FCIQMC used in this paper permits Ni to have continuous floating point magnitudes above the value of 1, which allows for a finer representation of the walker populations and reduces stochastic and systematic errors of the simulations.69,70 The purpose of FCIQMC is to generate many distributions of walkers, such that in the large Nw and long-time limit the expectation value of the number of walkers on each

Figure 1. FCIQMC algorithm showing spawning, death, and annihilation steps for each iteration.

When averaged over a long period of imaginary time, the walker population on each determinant, Ni, become proportional to the CI expansion coefficient, Ci. The averages associated with each determinant are, however, never explicitly obtained as the memory requirements to store these values would pose the exponential bottleneck that the method attempts to avoid. Instead, derived properties are accumulated once the walker distribution has reached a dynamic equilibrium. The most straightforward of these values is the projected energy, where the numerator and denominator of the expression EP = 1249

⟨∑i Ni⟨D0 Ĥ Di⟩⟩τ ⟨N0⟩τ

(22) DOI: 10.1021/acs.jctc.5b01190 J. Chem. Theory Comput. 2016, 12, 1245−1258

Article


to as the α band. Although Gouterman’s model can qualitatively describe the Q bands, it cannot provide a satisfactory quantitative interpretation due to the limited correlation recovered using only the four-orbitals. Gouterman’s orbitals are part of a more extensive π−π* system in which near-degeneracy effects are not negligible, suggesting that the entire π−π* system must be correlated to obtain quantitative accuracy. In the literature, it has been argued that even the σ orbitals are required.63,77 Metal-porphyrins are very interesting as biological systems make abundant use of them. Among these compounds, the Mg(II) porphyrins and Fe(II) porphyrins are the most common and occur in abundance in nature. While Mg(II) porphyrins are closed shell diamagnetic compounds, the Fe(II) porphyrins may appear in their high-spin (quintet), low-spin (singlet), and intermediate-spin (triplet) configurations depending on the coordinating ligands. To the best of our knowledge, while many six-coordinated Fe(II) porphyrins have been unambiguously characterized in terms of spin multiplicity and electronic configuration,78 a definitive assignment of the ground state of the four-coordinated system is still missing. Experimental studies on the Fe(II) tetraphenilporphyrin (FeTTP) compound agree that the ground state is a triplet but disagree in the details of the electronic configuration. A 3A2g ground state with configuration (dxy)2(dxz, dyz)2(dz2)2 was suggested by Mössbauer,79,80 magnetic,81 and H NMR82,83 measurements of the Fe(II) tetraphenilporphyrin (FeTPP). A 3 Eg state was instead suggested by Raman spectroscopy,84 with configuration (dxy)2(dxz, dyz)3(dz2)1 for the iron octaethylporphyrin (FeOEP). It is worth noting that the functionalization of the macrocycle may affect the relative ordering of the first lowlying states. A large number of theoretical studies have been reported for ferrous porphyrin.85−92 Most of them have identified a quintet ground state, 5A1g, suggesting either a systematic error in the theoretical means, related to the limited size of active space, or that the ground state for the unfunctionalized macrocycle may indeed be the high-spin manifold. Hirao et al.91 performed a CASSCF/MRMP2 investigation of the Fe(II) porphyrin system. The active space for their CASSCF calculations consisted of the five 3d orbitals and 6 π-type orbitals of the porphyrin ring and their electrons, CAS(10,11). Later, Hirao and Lindh90 enlarged the active space by adding two more πtype orbitals and four electrons, CAS(14,13). Both calculations resulted in a 5A1g ground state, with leading (dz2)2(dx2−y2)1(dxy)1(dxz)1(dyz)1 configuration. To further investigate the inconsistency between the MRMP2 results and experimental findings, Pierloot et al.92 performed CASPT2 and RASPT2 calculations, choosing three different active spaces. In CAS(8,11), six electrons from the iron and two from the σ Fe− N orbital were correlated. In order to reduce the electron− electron repulsion within the 3d orbitals, a set of five 4d orbitals was also added to the active space (the so-called “double-shell”). 5 A1g was again found to be the ground state, but the effect of the “double-shell” was to make the triplet 3A2g state nearly degenerate to the ground state, separated by only 0.09 eV. The same authors92 enlarged the active space by means of the RASPT2 method, up to a RAS(34,35), by also including the πorbitals, with only single and double excitations allowed from RAS1 and into RAS3. 5A1g was still found to be the ground state for the Fe(II) compound with the first triplet state (3A2g) lying 0.2 eV above. In their third choice of active space, the semi-core (3s3p) orbitals were included, whereas the π system was not,

are accumulated and averaged separately. (D0 is the reference (most populated) determinant, and ⟨N0⟩τ is the time-averaged number of walkers on D0. The sum in the numerator can be rigorously limited to the single and double excitations of D0, as all other matrix elements ⟨D0|Ĥ |Di⟩ vanish). Similarly, the oneand two-electron reduced density matrices are accumulated using the replica method37 to evaluate the CASSCF energy as for eq 2 and the matrices required for the Super-CI step. The target total walker population, Nw, represents an optimization parameter in addition to basis set and active space. In general, the larger the number of walkers is, the smaller is the systematic error associated with the initiator approximation. As discussed later in this manuscript, even for large systems and large active spaces, the variance of the energy and other properties such as the elements of the reduced density matrices quickly reduces with the number of walkers. The walker population is increased until negligible variation in the probed quantities (in particular the energy and the unitary transformations of the CASSCF procedure) is observed.

3. COMPUTATIONAL DETAILS 3.1. Porphyrins as Model System. The ground state of H2P is the singlet 1Ag state. All free-base porphyrins show a similar absorption spectrum with two absorption bands of moderate intensity in the visible region, namely, the Qx (1.98− 2.02 eV) and Qy (2.33−2.42 eV) bands, and a much stronger B band (also known as Soret band) in the near-UV region (3.13− 3.33 eV). The Qx and Qy bands are presumed to be singlet− singlet excitations associated with the 1Ag → 1B3u and 1Ag → 1 B2u transitions.73−76 The splitting of the Qx and Qy bands is ascribed to the fact that the x and y axes are not equivalent in the presence of the pyrrolic hydrogen atoms, as explained by Gouterman’s four-orbital model. In this model, each band comprises two singly excited configurations from the two highest occupied orbitals (HOMOs), 5b1u and 2au, to the lowest unoccupied molecular orbitals (LUMOs), 4b2g and 4b3g (Figure 2). When the two pyrrolic hydrogens are removed to form a cation or a metalloporphyrin, the asymmetry along the x and y axes is removed, the molecule shifts from D2h to D4h symmetry, and the frontier orbitals become degenerate. As a result, the Qx and Qy bands merge into a single Q-band, sometimes referred

Figure 2. HOMOs and LUMOs of porphyrin in Gouterman’s fourorbital model. X and Y excitation labels refer to the Qx and Qy absorption bands. 1250


Article


the CI space within the active space have been made, we believe that our results are not biased by fortuitous error cancellation. No σ orbitals have been included in the active space, and no dynamical correlation outside the active space has been accounted for. These two aspects may be regarded as the weak points of the simulations presented here and will be further addressed in a future study. In order to compare FCIQMC-CASSCF with deterministic CASSCF implementations, calculations have been performed with a much smaller active space, consisting of 14 active electrons and 14 active orbitals of π−π* character. For comparison purposes, RASSCF calculations have been performed in which 26 electrons and 24 active orbitals have been partially correlated. In the RAS(26,24) active space, 11 orbitals are in the RAS1 space, four orbitals in RAS2 (Gouterman’s orbitals), and nine in RAS3. Only single and double excitations out of RAS1 and into RAS3 have been allowed. Second-order perturbation theory using the optimized RASSCF wave function as reference has also been employed on this system. For the RASPT2 calculation, the first 24 low energy σ orbitals have been frozen. All virtual orbitals have been included in the PT2 approach. For the metal-porphyrins, a highly symmetric D4h idealized geometry was used (see SI.2). This geometry was derived from the FBP with small changes to enforce D4h symmetry. Calculations were carried out in the D2h point group (the largest point group Molcas can handle to date). In this point group, the π-orbitals belong to the b1u, b2g, b3g, and au irreducible representations, with the z axis perpendicular to the molecular plane. Table 2 shows the relationship between

leading to a (16,15) active space. Using this choice, the triplet 3 A2g was found to be the ground state with the quintet state only 0.04 eV above. It is worth noting that this result holds only in the adiabatic limit. When vertical excitations are evaluated at the 5A1g geometry, the quintet once again becomes the ground state with the triplets about 0.3 eV above. As larger active spaces are accessible using the FCIQMC-CASSCF approach, the π orbitals as well as the valence 3d orbitals of the iron center have been explicitly correlated, allowing for a deeper view of the interaction between the aromatic macrocycle and the metal center. 3.2. Geometry, Basis Set, and Active Space. The atomic coordinates for the free-base porphyrin have been taken from Nagashima et al.64 (coordinates given in SI.1) and were used for both the ground and excited states. The molecule was placed on the xy plane, the x axis containing the pyrrolic hydrogen atoms and constrained to the D2h point group, such that the π − π* orbital system belongs to the b1u, b2g, b3g, and au irreducible representations. We used generally contracted atomic natural orbitals (ANO) basis sets,93,94 obtained from the C,N(14s9p4d3f2g)/H(8s4p3d1f) primitive functions. Two contraction schemes were employed. In the smaller one, the primitive functions were contracted to C,N(3s2p1d)/H(2s1p), giving a basis set of split-valence double-ζ plus polarization quality, which will be referred to as ANO-VDZP. This basis set choice led to a total of 406 basis functions for the free-base porphyrin. A larger contraction scheme has also been adopted for the free-base porphyrin, containing C,N(4s3p2d1f)/H(3s2p1d) contracted functions. This basis set is of triple-ζ quality, providing 916 basis functions. It will be referred to as ANO-VTZP. As discussed by Gwaltney and Bartlett95 a f lexible basis set for this system should include polarization and diffuse functions, and our choice is certainly consistent with their statement. The Hartree−Fock SCF orbitals consist of 81 occupied and 325 vacant molecular orbitals in the ANO-VDZP basis set (835 vacant in the ANO-VTZP basis). For the CASSCF calculations, an active space was chosen with 26 electrons distributed among the 24 π valence orbitals (see Table 1 for details), for which the short notation “CAS(26,24)” will be used.

Table 2. Relationship between the D2h and D4h Point Groups

Inactive Active Secondary ANO-VDZP Secondary ANO-VTZP

B3u

B2u

B1g

B1u

B2g

B3g

Au

20 0

17 0

17 0

14 0

0 7

0 6

0 6

0 5

59

58

56

55

23

22

21

20

143

140

136

133

72

69

67

64

D4h

3d orbitals

Ag B3u, B2u B1g B1u B2g, B3g Au

A1g, B1g Eu A2g, B2g A2u, B2u Eg A1u,B1u

3dz2, 3dx2−y2 3dxy 3dxz, 3dyz

the irreducible representations of D2h and D4h, together with the distribution of the 3d orbitals among the irreducible representations. This mapping is required to be able to assign the term symbol to the iron(II) porphyrin states obtained in D2h. No geometry optimization was performed. Both metal atoms are in oxidation state (II). The ANO-VDZP basis set was also used for the metal atoms. Mg{17s12p6d2f2g} and Fe{21s15p10d6f4g2h} primitive functions were contracted to Mg{4s3p1d} and Fe{5s4p2d1f}. This resulted in a total of 414 and 430 basis functions for the Mg(II) and Fe(II) compounds, respectively. For the Mg(II) porphyrin, the effect of the metal on the Qband has been analyzed. Only the B3u component will be discussed here, as identical results were obtained for the B2u component, as the Qx and Qy bands in H2P are degenerate for metal-porphyrins. The same (26,24) active space was chosen as was used in H2P, as the Mg orbitals are not expected to have a multiconfigurational nature and do not mix with the π orbitals of the macrocycle. For the Fe(II) porphyrin, our main focus was the characterization of the ground state. A CAS(32,29) was chosen for this system, where the five 3d orbitals from the iron atom

Table 1. Distribution of Molecular Orbitals among Inactive, Active, and Secondary Spaces for Each Irreducible Representation Ag

D2h

The wave function has been optimized for the ground state, Ag, and the low-lying 1B3u and 1B2u excited states responsible for the Q bands. No orbitals have been frozen or deleted at the CASSCF level of theory. Canonical Hartree−Fock orbitals have been used as the starting orbitals for the FCIQMC-CASSCF calculations. To our knowledge, the choice of both basis set and active space represents the largest one- and many-body basis sets ever used for correlated methods applied to this compound. Further, since no assumptions or truncations of 1

1251


Article


between the two approaches at each point is of the order of 10−5 Eh. Only a small improvement is observed when the number of walkers is increased to 5 × 105, indicating that sufficient walkers are already present to represent the structure of the wave function in the Hilbert space for this model system. The small size of the deviations from deterministic-CASSCF demonstrate the robustness of the FCIQMC-CASSCF approach and its ability to provide results to chemical accuracy. 4.2. Q Bands for Free-Base Porphyrin. The CASSCF optimization curves for the CAS(26,24) of H2P in its ground state, 1Ag, and the two 1B3u and 1B2u excited states are given in Figure 4.

and their six electrons have been added to the CAS(26,24) of the free-base porphyrin. A detailed analysis of the effects related to the extended one-particle basis set, semi-core, and doubleshell orbitals goes beyond the scope of the present paper, which is aiming to show that larger CASSCF calculations can be performed, and will be the topic of a later study. Within FCIQMC, the initiator formulation of the method has been used35,36 with a threshold value of na = 3.0 together with the semi-stochastic method38,69 using a deterministic subspace consisting of |+| = 10000 most populated determinants. For all calculations, a single determinantal wave function has been used as reference. The calculations were run in replica mode37 in order to sample the reduced density matrices. At the early stages of the CASSCF prodecure, each replica consisted of 2 × 105 walkers (for the Fe(II) porphyrin 5 × 105 walkers have been used). After CASSCF convergence was obtained, more refined solutions were obtained by increasing the target number of walkers first to 5 × 106 and then 10 × 106 in consecutive steps. For the iron(II) states, up to 20 × 106 have been used. The entire walker dynamics occurred in 2 × 105 iterations, which was sufficient for all simulations to equilibrate (see also SI.3 and SI.4 for more details). Reduced density matrices were sampled after the first 1 × 105 iterations. The time-step, Δτ, was found via an automatic search procedure39 for each simulation and took typical values in the range 5 − 10 × 10−4 a.u. For the largest active spaces we studied here (32 electrons in 29 orbitals, for the FeP), the CPU costs were approximately as follows: with 2 × 105 walkers, each FCIQMC iteration took 4 × 10−3 seconds on 160 processors (8 nodes), increasing linearly with the number of walkers to 0.2 s with 10 × 106 walkers. A typical FCIQMC simulation, therefore, took roughly 1000 s for each orbital rotation step of the CASSCF procedure at 5 × 105 walkers and 20,000 s at 10 × 106 walkers.

4. RESULTS 4.1. Stochastic- and Deterministic-CASSCF. Figure 3 presents the CASSCF optimization of the (14,14) active space for the 1Ag ground state of the free-base porphyrin, comparing the deterministic and the stochastic pathways. The difference

Figure 4. Optimization of the (26,24) active space for the free-base porphyrin within the ANO-VDZP basis set. The three subplots refer, in order, to the 1Ag, 1B3u, and 1B2u states. Red squares, blue circles, and green triangles refer to population dynamics with 2 × 105, 5 × 106, and 10 × 106 walkers, respectively.

Hartree−Fock canonical orbitals in the ANO-VDZP basis set were used as the starting orbitals for the CASSCF optimization, and 2 × 105 walkers were employed for the initial dynamics. After about seven macroiterations, the CASSCF energies are converged to within 10−4 Eh. Starting from the final natural orbitals, a new CASSCF optimization was performed, increasing the number of walkers to 5 × 106. Small orbital rotations were observed, and after five more iterations, CASSCF energies were converged to within 10−5 Eh. We found it quite surprising that a small walker distribution is able to generate a convenient averaged field to allow for an effective orbital optimization. Unless further analysis of this approach disprove this behavior, these authors suggest to always perform a rather fast and inexpensive FCIQMC-CASSCF calculation with a relatively small number of walkers (in the range 105− 106) and use the optimized natural orbitals as starting orbitals in subsequent FCIQMC-CASSCF calculations with enlarged walker population. This procedure guarantees a fast initial orbital optimization and a limited number of macroiterations at the high-population regime. Enlarging the walker population has the effect of improving the sampling of the Hilbert space,

Figure 3. Optimization of the (14,14) active space for the free-base porphyrin within the ANO-VDZP basis set. Differences between the deterministic and the stochastic schemes are given in the lower subplot; 200 and 500 K labels refer to a population dynamic with 2 × 105 and 5 × 105 walkers, respectively. 1252


Article


The RAS(26,24)SCF and RASPT2 excitation energies are also reported in Table 3. It is interesting to note that the truncation of the CI expansion due to the RAS constraints (only single and double excitations from RAS1 and to RAS3) has a noticeable impact on the excitation energies, such that the Qx band is overestimated by about 0.5 eV and the Qy band by about 0.1 eV relative to FCIQMC-CASSCF. The RASPT2 approach appears more sensitive for the Qy band than for the Qx band. In fact, RASPT2 does not improve the excitation energy for the 1Ag → 1B3u transition over the RASSCF results. On the contrary, the position of the Qy band at RASPT2 level is about 0.5 eV closer to the experimental value when compared to the RASSCF result. This finding reinforces our previous conclusion that dynamical correlation is more important for the Qy band than for the Qx band. One open question still remains to these authors: why would two qualitatively similar states behave so differently in terms of correlation effects? The optimized active natural orbitals for the 1Ag state at CAS(26,24) level of theory are depicted in Figure 5. The natural orbitals for the 1B3u and 1B2u states have similar shapes. From the occupation numbers of the natural orbitals, the highly multiconfigurational nature of this compound in its ground state and its excited states is evident. The ground state is dominated by the (b1u)10(b2g)6(b3g)6(au)4 configuration with a weight of 34.4%. The other relevant configurations are mainly double excitations from the dominant one, namely, the (b1u)10(b2g)8(b3g)6(au)2 (5.0%) and (b1u)8(b2g)6(b3g)8(au)4 (4.9%) configurations. The 1B3u state is dominated by (b1u)9(b2g)7(b3g)6(au)4 (15.4%) and (b1u)10(b2g)6(b3g)7(au)3 (14.1%), with the following configurations accounting for less the 2% each. The 1 B 2 u state is dominated by (b1u)10(b2g)7(b3g)6(au)3 (18.9%) and (b1u)9(b2g)6(b3g)7(au)4 (12.7%). The extraordinarily low weight of the leading determinants is a reflection of the multiconfigurational nature of this system. These findings are in agreement with Gouterman’s four-orbital model. In fact, the 1B3u state arises from single electron excitations from the 5b1u orbital to the 4b2g orbital and from the 2au orbital to the 4b3g orbital, while the 1B2u state arises from single electron excitations from orbital 5b1u to orbital 4b3g and from orbital 2au to orbital 4b2g (Figure 2). It is worth noting that the band gap estimated by FCIQMCCASSCF for the 1Ag → 1B2u transition is at about the value experimentally referred to the By band. It is interesting to note that also the occupation number of the natural orbitals are close to the ones reported by Ohno et al.63 This might also be an indication that the FCIQMC-CASSCF wave function was optimized to the second 1B2u state (the one responsible for the Soret band). We tried various optimization schemes to verify such possibility. However, all our attempts led to the same wave function. 4.3. Q-Band of Mg(II) Porphyrin. Already at the Hartree− Fock level of theory, it is clear that the interaction between the magnesium orbitals and the π−π* system of the aromatic macrocycle is small. The doubly occupied orbitals of the magnesium lie on the low side of the inactive space. In the virtual manifold, some MOs have contributions from the Mg atomic orbitals, but they are very high in energy and not mixed with the orbitals of the π−π* system. As a consequence, the π orbitals remain structurally similar to those of the free-base porphyrin, with differing orbital and transition energies. Table 4 collects excitation energies corresponding to the 1A1g → 1Eu transition (in D4h nomenclature) for the Mg(II)

and consequently, more of the correlation energy was recovered, with total electronic energies being lowered by about 10 mEh for the FBP compound. When the walker population was further increased to 10 × 106, total CASSCF energies only changed by 1 mEh. We consider this walker population sufficient to provide a reliable walker population dynamics for this system. A further increase in the total number of walkers to 20 × 106 only impacts the total energies by less than a mEh, below the required accuracy. Excitation energies for each choice of walker distribution are reported in Table 3. The impact of the total number of walkers Table 3. Excitation Energies (eV) Corresponding to 1Ag → 1 B3u and 1Ag → 1B2u Transitions method ANO-VDZP 200 K ANO-VDZP 5M ANO-VDZP 10M ANO-VTZP 5M ANO-VTZP 10M ANO-VDZP RASSCF ANO-VDZP RASPT2 exp.96

Ag → 1B3u (Qx)

1

2.19 2.05 2.02 2.08 2.06 2.55 2.54 2.02

1

Ag → 1B2u (Qy) 3.29 3.25 3.24 3.27 3.24 3.38 2.81 2.39

is smaller on the relative energies than the absolute energies. Although enlarging from 2 × 105 (200 K) to 5 × 106 (5 M) walkers has a noticeable effect, when the total population is increased to 10 × 106 (10 M) walkers, deviations of only 0.03 and 0.01 eV are observed for the transitions. The basis set also plays a minor role in the evaluation of the absorption bands. The 1Ag → 1B3u excitation energy is in very good agreement with the experimental value. The most surprising finding is that the 1Ag → 1B2u transition, that has been assigned to the Qy band in the literature, is located at ∼3.3 eV, blue-shifted by 0.9 eV with respect to the experimental value of 2.4 eV. The calculated value is sensitive to neither the basis set nor the walker population. We associate this discrepancy to two potential factors: (1) size of the active space and (2) geometry. In the literature, it is widely argued that σ-orbitals play a major role in providing quantitative accuracy for this system, while the active space of our choice includes only the valence π−π* system. When considering only the π−π* transitions, Ohno et al.63 obtained excitations energies around 0.2 eV larger for the Qx band than those obtained using FCIQMC-CASSCF. The Qy band is also overestimated in their approach when only the π orbitals are correlated. When they correlated σ orbitals both bands were well reproduced. In this respect, our finding confirms their conclusion: correlation outside the π-system needs to be addressed in order to have a complete view on correlation effects for this system. Experimental observations of H2P cooled in supersonic expansion, observed by fluorescence excitation, show a complex, structured Qy band with “dramatic” line broadening, an order of magnitude larger than observed for the Qx band.97,98 This broadening has been attributed to an adiabatic vibronic interstate coupling between the 1B3u and 1B2u states. Our excitation energy was estimated at fixed geometry and did not account for interstate coupling. Therefore, it can also be argued that relaxation effects and coupling of states may be invoked to explain the deviation of the simulated Qy band with respect to the experimental one. 1253


Article


Figure 5. Active natural orbitals for the CAS(26,24) of the free-base porphyrin. The occupation numbers of the natural orbitals for the 1Ag state (no brackets), 1B3u state (square brackets), and 1B2u state (round brackets) are given.

Table 4. Excitation Energies (eV) Corresponding to 1A1g → 1 Eu Transition method (26,24) 200 K (this work) (26,24) 5 M (this work) (26,24) 10 M (this work) SAC−CI99 MRPT100 CAS(15,18)SCF101 CAS(15,18)PT2101 exp.102

correction (PT2) underestimated the excitation energy by 0.42 eV. The improvement of our approach over the CAS(15,18)SCF and the PT2 results is notable, deviating from the experimental value by only about 0.3 eV. The SAC−CI approach proposed by Hasegawa is in better agreement with the experimental value, underestimating the excitation energy by only 0.2 eV. Similarly to the free-base porphyrin, the ground state of the Mg(II) porphyrin is dominated by the Hartree−Fock configuration, (b1u)10(b2g )6(b3g) 6(au) 4, which contributes 39.4% (using D2h nomenclature). Other relevant configurations (weight >4%) consist of pairs of single excitations from the reference configuration. This again shows the highly multiconfigurational nature of this class of compounds due to the extended aromatic macrocycle. The 1Eu excited state shows even more enhanced multiconfigurational character, in the same way as the excited states of the metal-free counterpart. By comparison, the CASSCF results from ref 101 show much less multiconfigurational character for the Mg(II) porphyrin compound, probably due to the limited size of the active space. Rubio et al. also discussed the somewhat “positive” effect of removing the polarization

A1g → 1Eu

1

2.69 2.52 2.50 2.01 2.00 3.06 1.78 2.20

porphyrin obtained using FCIQMC-CASSCF on a (26,24) active space (same as for H2P). Experimental values and some other theoretical values available in the literature are given for comparison. Rubio et al. show there is a non-negligible dependence of the excitation energy on the size of the active space (see Table 4 in ref 101). Their largest active space, consisting of 15 electrons in 18 active orbitals, gave an excitation energy of 3.06 eV at the CASSCF level of theory. On the other hand, the perturbative 1254


Article


calculations the separation was fixed at 2.05 Å. A similar length was used in their work (2.03−2.07 Å). By contrast, the experimental Fe−N bond length in FeTPP obtained by X-ray diffraction is 1.97 Å. Geometry optimization using a CAS(32,29) active space to consider the relaxation of the aromatic cycle would be an interesting further study. Functionalization of the macrocycle, as for example, in TPP or OEP, may directly (by modifying the electron density) or indirectly (for example by changing the Fe−N bond length) exchange the relative order of states. In this respect, it is worth noting that another fourcoordinate ferrous porphyrin, iron(II) octamethyltetrabenzporphyrin (FeOTBP), is surprisingly different from FeTPP and FeOEP. Experimental data have supported an 5A1g ground state for this compound.103 The Fe−N bond length for this compound has been estimated to be 2.03 Å.87 This finding is in line with Lindh’s suggestion that longer Fe−N bond lengths favor high-spin multiplets and supports also our findings. Spin− orbit coupling effects may also be invoked to discuss the relative energy of quasi-degenerate states. This aspect has been already discussed by Pierloot et al.92 Finally, although the entire set of valence orbitals has been correlated using FCIQMCCASSCF, important correlations (semi-core and double-shell effects) are still missing. Pierloot et al.92 instead accounted for the correlation outside the active space via the CASPT2 approach. These aspects have not been addressed in the present analysis.

functions from the basis set (see Table 5 in ref 101). They argued that the lack of polarization functions in the basis set fortuitously improves agreement with experiment. We did not investigate the effects of the truncation of the one-electron basis any further. By looking at the Mulliken charges, Rubio et al. concluded that there exists back-donation from the macrocycle to the metal site and a net charge of +1.10 on the Mg atom instead of the formal +2 charge. Our optimized CAS(26,24) calculation gives a +1.45 Mulliken charge on the central metal atom, confirming to some extent the back-donation. 4.4. Fe(II) Porphyrin Case. Using FCIQMC-CASSCF with a large (32,29) active space, the quintet 5A1g state is confirmed as the ground state for the ferrous porphyrin, with the dominant configuration, (3dz2)2(3dx2−y2)1(3dyz)1(3dxz)1(3dxy)1, contributing 37.3%. Surprisingly, the other relevant configurations (with weight >5%) are charge-transfer configurations in which electrons from the metal center are excited into the π orbitals. The dominant configuration agrees with the prediction of Lindh et al., but no significant contributions from the chargetransfer configurations were observed in their calculations. They concluded that the iron d-electron space and the porphyrin π-electron space are well separated, whereas our findings show the opposite. The FCIQMC-CASSCF results show a close interaction between the aromatic ring and the metal center via non-negligible contributions from chargetransfer configurations. The CASSCF convergence for the iron(II) porphyrin is shown in Figure 6. It is noticeable that during the first five

5. CONCLUSIONS A new FCIQMC-CASSCF method has been presented as a stable and reliable tool to tackle strongly correlated systems with a multitude of entangled electrons and orbitals. Application to the free-base porphyrin, Mg(II) porphyrin, and Fe(II) porphyrin demonstrates the capability of the method. Active spaces containing up to 32 active electrons and 29 orbitals are easily accessible by using the FCIQMC method as the CI eigensolver. One-body bases containing up to 916 contracted functions have been employed without difficulty. By not requiring the explicit evaluation and storage of orbital Hessian elements, the Super-CI approach has a significant advantage over second-order schemes. The comparison of the Newton−Raphson and Super-CI approach for large active spaces is of particular interest. For free-base porphyrin, good agreement is obtained for the Qx absorption band, while a considerable discrepancy is observed for the Qy band. The authors do not fully understand this result and presume that two factors may play an important role in this large deviation: (1) missing correlation energy and (2) adiabatic vibronic interstate coupling. Satisfactory results have been obtained for the Mg(II) porphyrin. For the most challenging system presented, the Fe(II) porphyrin, 5A1g has been confirmed as the ground state. This is in agreement with the Lindh and Pierloot results in the diabatic limit, but further investigations are required as the effects of the double shell and semi-core orbitals need to be addressed. This would require a CAS(40,42), certainly accessible by the method presented here. Geometric distortions and functionalization of the aromatic ring may play important roles in the relative stability of the low-lying quintet and triplet states. Finally, it is worth mentioning that the coupling of this approach to a posteriori corrections aiming at the missing correlation will be the main focus of our upcoming research.

Figure 6. Convergence path for the Fe(II) porphyrin in a (32,29) active space.

iterations the CASSCF energy decreases by about 0.75 Eh (2 eV). Only small deviations are observed during later iterations. This first stage was relatively cheap as only 5 × 105 walkers were used, and when the walker population was increased, only negligible orbital rotations were observed. The lower energy is purely due to the more quantitatively effective sampling of the accessible configurations in the Hilbert Space. Active natural orbitals and their occupation numbers for the 5A1g ground state are shown in Figure 7. It has been observed87 that small perturbations might change the relative ordering of the states. Lindh et al. discussed the effect of the Fe−N bond length on the position of the 3dx2−y2 orbital in the energy spectrum90 and concluded that a larger Fe−N separation might favor the high-spin states. For our 1255


Article


Figure 7. Active natural orbitals for the CAS(32,29) of the iron(II) porphyrin in its 5A1g ground state.

■

ASSOCIATED CONTENT

evaluation of one- and two-electron integrals in the parallel framework, and Armin Schuhmacher for keeping the MPI-FKF cluster up and running at all times even when overloaded.

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b01190. Contains cartesian coordinates (Å) for the free-base porphyrin (Listing 1) and metal-porphyrins (Listing 2), Molcas input example (Listing 3), and list of relevant keywords used within the FCIQMC program for the calculations here presented (Listing 4). (PDF)

■

■

REFERENCES

(1) Loew, G. H.; Harris, D. L. Chem. Rev. 2000, 100, 407−420. (2) Jensen, K. P.; Ryde, U. J. Biol. Chem. 2004, 279, 14561−14569. (3) Hartree, D. R. Math. Proc. Cambridge Philos. Soc. 1928, 24, 89− 110. (4) Fock, V. Z. Phys. 1930, 61, 126−148. (5) Roothaan, C. C. J. Rev. Mod. Phys. 1951, 23, 69−89. (6) Roothaan, C. C. J. Rev. Mod. Phys. 1960, 32, 179−185. (7) Kohn, W.; Sham, L. Phys. Rev. 1965, 140, A1133−A1138. (8) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. M. Chem. Phys. 1980, 48, 157−173. (9) Roos, B. O. Int. J. Quantum Chem. 1980, 18, 175−189. (10) Siegbahn, P. E. M.; Almlöf, J.; Heiberg, A.; Roos, B. O. J. Chem. Phys. 1981, 74, 2384−2396. (11) Siegbahn, P.; Heiberg, A.; Roos, B. O.; Levy, B. Phys. Scr. 1980, 21, 323−327. (12) Shavitt, I. Int. J. Quantum Chem. 1977, 12 (S11), 131−148. (13) Shavitt, I. Int. J. Quantum Chem. 1978, 14 (S12), 5−32. (14) Olsen, J.; Roos, B. O.; Jorgensen, P.; Jensen, H. J. A. J. Chem. Phys. 1988, 89, 2185−2192.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (G.L.M.). *E-mail: [email protected] (A.A.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge Valera Veryazov and Victor Vysotskiy for making parallel libraries available in Molcas, facilitating the integration of the NECI program with the Molcas modules, Roland Lindh for making available the 1256


Article

Journal of Chemical Theory and Computation (15) Malmqvist, P.-Å; Rendell, A.; Roos, B. O. J. Phys. Chem. 1990, 94, 5477−5482. (16) Ma, D.; Li Manni, G.; Gagliardi, L. J. Chem. Phys. 2011, 135, 044128. (17) Li Manni, G.; Ma, D.; Aquilante, F.; Olsen, J.; Gagliardi, L. J. Chem. Theory Comput. 2013, 9, 3375−3384. (18) Vogiatzis, K. D.; Li Manni, G.; Stoneburner, S. J.; Ma, D.; Gagliardi, L. J. Chem. Theory Comput. 2015, 11, 3010−3021. (19) Andersson, K.; Malmqvist, P.-Å; Roos, B. O.; Sadlej, A. J.; Wolinski, K. J. Phys. Chem. 1990, 94, 5483−5488. (20) Andersson, K.; Malmqvist, P.-Å; Roos, B. O. J. Chem. Phys. 1992, 96, 1218−1226. (21) Malmqvist, P.-Å; Pierloot, K.; Shahi, A. R. M.; Cramer, C. J.; Gagliardi, L. J. Chem. Phys. 2008, 128, 204109. (22) Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen, J.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. 2014, 10, 3669−3680. (23) Carlson, R. K.; Li Manni, G.; Sonnenberger, A. L.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. 2015, 11, 82−90. (24) Carlson, R. K.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. 2015, 11, 4077−4085. (25) Ghosh, S.; Sonnenberger, A. L.; Hoyer, C. H.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. 2015, 11, 3643−3649. (26) Keller, S.; Dolfi, M.; Troyer, M.; Reiher, M. J. Chem. Phys. 2015, 143, 244118. (27) Kurashige, Y.; Chan, G. K.-L.; Yanai, T. Nat. Chem. 2013, 5, 660−666. (28) Chan, G. K.-L.; Head-Gordon, M. J. Chem. Phys. 2002, 116, 4462−4476. (29) Schollwöck, U. Rev. Mod. Phys. 2005, 77, 259−315. (30) Marti, K. H.; Ondík, I. M.; Moritz, G.; Reiher, M. J. Chem. Phys. 2008, 128, 014104. (31) Hedegård, E. D.; Knecht, S.; Kielberg, J. S.; Jensen, H. J. A.; Reiher, M. J. Chem. Phys. 2015, 142, 224108. (32) Hu, W.; Chan, G. K.-L. J. Chem. Theory Comput. 2015, 11, 3000−3009. (33) Sharma, S.; Sivalingam, K.; Neese, F.; Chan, G. K.-L. Nat. Chem. 2014, 6, 927−933. (34) Booth, G. H.; Thom, A. J. W.; Alavi, A. J. Chem. Phys. 2009, 131, 054106. (35) Cleland, D.; Booth, G. H.; Alavi, A. J. Chem. Phys. 2010, 132, 041103. (36) Cleland, D.; Booth, G. H.; Alavi, A. J. Chem. Phys. 2011, 134, 024112. (37) Overy, C.; Booth, G. H.; Blunt, N. S.; Shepherd, J. J.; Cleland, D.; Alavi, A. J. Chem. Phys. 2014, 141, 244117. (38) Blunt, N. S.; Smart, S. D.; Kersten, J. A.-F.; Spencer, J. S.; Booth, G. H.; Alavi, A. J. Chem. Phys. 2015, 142, 184107. (39) Booth, G. H.; Smart, S. D.; Alavi, A. Mol. Phys. 2014, 112, 1855−1869. (40) Booth, G. H.; Cleland, D. M.; Thom, A. J. W.; Alavi, A. J. Chem. Phys. 2011, 135, 084104. (41) Blunt, N. S.; Smart, S. D.; Booth, G. H.; Alavi, A. J. Chem. Phys. 2015, 143.13411710.1063/1.4932595 (42) Blunt, N. S.; Alavi, A.; Booth, G. H. Phys. Rev. Lett. 2015, 115, 050603. (43) Thomas, R. E.; Sun, Q.; Alavi, A.; Booth, G. H. J. Chem. Theory Comput. 2015, 11, 5316−5325. (44) Levy, B. Int. J. Quantum Chem. 1970, 4, 297−313. (45) Hinze, J. J. Chem. Phys. 1973, 59, 6424−6432. (46) Aquilante, F.; et al. J. Comput. Chem. 2016, 37, 506−541. (47) Walter, M. G.; Rudine, A. B.; Wamser, C. C. J. Porphyrins Phthalocyanines 2010, 14, 759−792. (48) Yella, A.; Lee, H.-W.; Tsao, H. N.; Yi, C.; Chandiran, A. K.; Nazeeruddin, M.; Diau, E. W.-G.; Yeh, C.-Y.; Zakeeruddin, S. M.; Grätzel, M. Science 2011, 334, 629−634. (49) Swierk, J. R.; Méndez-Hernández, D. D.; McCool, N. S.; Liddell, P.; Terazono, Y.; Pahk, I.; Tomlin, J. J.; Oster, N. V.; Moore, T. A.; Moore, A. L.; Gust, D.; Mallouk, T. E. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 1681−1686.

(50) De Wael, K.; Adriaens, A. Talanta 2008, 74, 1562−1567. (51) Dogutan, D. K.; McGuire, R.; Nocera, D. G. J. Am. Chem. Soc. 2011, 133, 9178−9180. (52) Feringa, B. L. Molecular Switches; Wiley Online Library, 2001; Vol. 42. (53) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Chem. Rev. 2005, 105, 2279−2328. (54) Shaik, S.; Hirao, H.; Kumar, D. Acc. Chem. Res. 2007, 40, 532− 542. (55) Siegbahn, P. E. M.; Borowski, T. Acc. Chem. Res. 2006, 39, 729− 738. (56) Nooijen, M.; Bartlett, R. J. J. Chem. Phys. 1997, 106, 6449−6455. (57) Kowalski, K.; Krishnamoorthy, S.; Villa, O.; Hammond, J. R.; Govind, N. J. Chem. Phys. 2010, 132, 154103. (58) Aspuru-Guzik, A.; El Akramine, O.; Grossman, J. C.; Lester, W. A. J. J. Chem. Phys. 2004, 120, 3049−3050. (59) Serrano-Andrés, L.; Merchán, M.; Rubio, M.; Roos, B. O. Chem. Phys. Lett. 1998, 295, 195−203. (60) Merchán, M.; Ortí, E.; Roos, B. O. Chem. Phys. Lett. 1994, 226, 27−36. (61) Tokita, Y.; Hasegawa, J.; Nakatsuji, H. J. Phys. Chem. A 1998, 102, 1843−1849. (62) Nakatsuji, H.; Hasegawa, J.; Hada, M. J. Chem. Phys. 1996, 104, 2321−2329. (63) Yamamoto, Y.; Noro, T.; Ohno, K. Int. J. Quantum Chem. 1992, 42, 1563−1575. (64) Nagashima, U.; Takada, T.; Ohno, K. J. Chem. Phys. 1986, 85, 4524−4529. (65) Parusel, A. B. J.; Grimme, S. J. Porphyrins Phthalocyanines 2001, 5, 225−232. (66) Banerjee, A.; Grein, F. J. Chem. Phys. 1977, 66, 1054−1062. (67) Ruedenberg, K.; Cheung, L. M.; Elbert, S. T. Int. J. Quantum Chem. 1979, 16, 1069−1101. (68) Banerjee, A.; Grein, F. Int. J. Quantum Chem. 1976, 10, 123− 134. (69) Petruzielo, F. R.; Holmes, A. A.; Changlani, H. J.; Nightingale, M. P.; Umrigar, C. J. Phys. Rev. Lett. 2012, 109, 230201. (70) Overy, C.; Booth, G. H.; Blunt, N. S.; Shepherd, J. J.; Cleland, D. M.; Alavi, A. J. Chem. Phys. 2014, 141, 244117. (71) Anderson, J. B. J. Chem. Phys. 1975, 63, 1499−1503. (72) Anderson, J. B. J. Chem. Phys. 1976, 65, 4121−4127. (73) Gouterman, M. J. Mol. Spectrosc. 1961, 6, 138−163. (74) Falk, J. E. Porphyrins and Metalloporphyrins; Elsevier: Amsterdam, 1964. (75) Smith, K. M. Porphyrins and Metalloporphyrins; Elsevier: Amsterdam, 1975. (76) Weiss, C.; Kobayashi, H.; Gouterman, M. J. Mol. Spectrosc. 1965, 16, 415−450. (77) Merchán, M.; Ortí, E.; Roos, B. O. Chem. Phys. Lett. 1994, 221, 136−144. (78) Scheidt, W. R.; Reed, C. A. Chem. Rev. 1981, 81, 543−555. (79) Collman, J. P.; Hoard, J. L.; Kim, N.; Lang, G.; Reed, C. A. J. Am. Chem. Soc. 1975, 97, 2676−2681. (80) Lang, G.; Spartalian, K.; Reed, C. A.; Collman, J. P. J. Chem. Phys. 1978, 69, 5424−5427. (81) Boyd, P. D. W.; Buckingham, A. D.; McMeeking, R. F.; Mitra, S. Inorg. Chem. 1979, 18, 3585−3591. (82) Goff, H.; La Mar, G. N.; Reed, C. A. J. Am. Chem. Soc. 1977, 99, 3641−3647. (83) Mispelter, J.; Momenteau, M.; Lhoste, J. M. J. Chem. Phys. 1980, 72, 1003−1012. (84) Kitagawa, T.; Teraoka, J. Chem. Phys. Lett. 1979, 63, 443−445. (85) Matsuzawa, N.; Ata, M.; Dixon, D. A. J. Phys. Chem. 1995, 99, 7698−7706. (86) Kozlowski, P. M.; Spiro, T. G.; Bérces, A.; Zgierski, M. Z. J. Phys. Chem. B 1998, 102, 2603−2608. (87) Liao, M.-S.; Scheiner, S. J. Chem. Phys. 2002, 116, 3635−3645. (88) Smith, D. M. A.; Dupuis, M.; Straatsma, T. P. Mol. Phys. 2005, 103, 273−278. 1257


Article

Journal of Chemical Theory and Computation (89) Chen, H.; Lai, W.; Shaik, S. J. Phys. Chem. B 2011, 115, 1727− 1742. (90) Choe, Y.-K.; Nakajima, T.; Hirao, K.; Lindh, R. J. Chem. Phys. 1999, 111, 3837−3844. (91) Choe, Y.-K.; Hashimoto, T.; Nakano, H.; Hirao, K. Chem. Phys. Lett. 1998, 295, 380−388. (92) Vancoillie, S.; Zhao, H.; Tran, V. T.; Hendrickx, M. F. A.; Pierloot, K. J. Chem. Theory Comput. 2011, 7, 3961−3977. (93) Widmark, P.-O.; Malmqvist, P.-Å; Roos, B. O. Theor. Chem. Acc. 1990, 77, 291−306. (94) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2004, 108, 2851−2858. (95) Gwaltney, S. R.; Bartlett, R. J. J. Chem. Phys. 1998, 108, 6790− 6798. (96) Edwards, L.; Dolphin, D. H.; Gouterman, M.; Adler, A. D. J. Mol. Spectrosc. 1971, 38, 16−32. (97) Even, U.; Jortner, J. J. Chem. Phys. 1982, 77, 4391−4399. (98) Arabei, S.; Kuzmitsky, V.; Solovyov, K. Opt. Spectrosc. 2007, 102, 692−704. (99) Hasegawa, J.; Hada, M.; Nonoguchi, M.; Nakatsuji, H. Chem. Phys. Lett. 1996, 250, 159−164. (100) Hashimoto, T.; Choe, Y.-K.; Nakano, H.; Hirao, K. J. Phys. Chem. A 1999, 103, 1894−1904. (101) Rubio, M.; Roos, B. O.; Serrano-Andrés, L.; Merchán, M. J. Chem. Phys. 1999, 110, 7202−7209. (102) Starukhin, A.; Shulga, A.; Waluk, J. Chem. Phys. Lett. 1997, 272, 405−411. (103) Sams, J. R.; Tsin, T. B. Chem. Phys. Lett. 1974, 25, 599−601.

1258


Combining the Complete Active Space Self-Consistent Field Method

Recommend Documents