Fragment-Based Restricted Active Space Configuration Interaction with Second-order corrections Embedded in Periodic Hartree-Fock Wavefunction
We present a computational scheme for restricted-active-space configuration interaction (RASCI) calculations combined with second-order perturbation theory (RASCI-PT2) on a fragment of a periodic system embedded in the periodic Hartree-Fock (HF) wavefunction. This method allows one to calculate the electronic structure of localized strongly-correlated features in crystals and surfaces. The scheme was implemented via an interface between the Cryscor and Q-Chem codes. To evaluate the performance of the embedding method, we explored dissociation of fluorine atom from a lithium fluoride surface and partially fluorinated graphane layer. The results show that RASCI 1 and RASCI-PT2 embedded in periodic HF are able to produce well-behaved potential energy surfaces and accurate dissociation energies.
Introduction<!>Theoretical Background Embedding Workflow<!>RASCI and RASCI-PT2 with hole and particle approximation<!>Crystal geometry<!>Crystal geometry Basis set<!>LiF<!>HF-graphane<!>Fragment<!>Potential energy curves<!>Convergence with fragment and supercell sizes<!>Conclusions
<p>Theoretical modeling plays an important role in the development of new materials. For solid state calculations, Kohn-Sham (KS) density functional theory (DFT) is presently the most popular electronic structure method due to a good balance between accuracy and computational cost. However, standard KS-DFT within currently available functional approximations is known to fail for multireference systems 1,2 e.g., transition metal complexes, open-shell defects, systems far from equilibrium, etc. In this situation, good accuracy can be reached with the so-called multireference methods, which are, however, computationally very expensive and hence applications to periodic systems are limited.</p><p>In solids, multireference features have often local character, e.g. a local excitation, a vacancy, a catalytic reaction on a surface, etc. One possible approach in this case is to build a finite cluster by carving out the region of interest from the bulk. In this way one can employ molecular methods implemented in existing molecular quantum-chemical packages.</p><p>Unfortunately, validity of the finite cluster model strongly depends on the type of the system.</p><p>It is in principle applicable to either purely covalent solids, where the dangling bonds can be terminated with hydrogen atoms, or purely ionic ones, where no bonds need to be cut.</p><p>However, for systems with a complex binding pattern or strongly disperse frontier bands the cluster approach is often unreliable. Besides, the long-range electrostatic interactions, present in solids, are missing in isolated (non-embedded) clusters, which can have a profound effect on the results.</p><p>Embedding approaches can be seen as an extension of the cluster model to study local effects in crystals. Several embedding methods have been developed. [3][4][5][6][7][8][9][10][11][12][13][14] They can be categorized into four classes: wavefunction partitioning, density-based, high-level-cluster embedded in low-level-periodic model and quantum/classical force field approaches. These embedding classes can also be combined in hybrid schemes. As an example of density-based embedding one can mention the WFT-in-DFT embedding models, pioneered by Carter and co-workers. 6,15,16 Since then different flavors of DFT embedding have been developed and applied to molecular 9,11,[17][18][19][20][21][22][23][24][25][26] and periodic systems. 15,16,27,28 Wavefunction partitioning and the related WFT-in-WFT embedding scheme can be seen as an alternative to the DFT embedding. This approach is presently a very active field of research. 14,17,[29][30][31][32][33][34] One of those methods is fragment-based direct-local-ring coupled-cluster doubles (d-LrCCD) embedded in the periodic Hartree-Fock wavefunction, developed by some of us. 35 Although d-LrCCD itself does not guarantee more accurate results than MP2, it is much more stable for crystals with small band gap. For such systems it becomes instrumental to build hierarchical schemes that include higher-order corrections, evaluated on finite clusters (i.e., a hybrid embedding scheme, see above). 36,37 Very recently the periodic HF embedding was employed in conjunction with the Full CI quantum Monte Carlo (FCIQMC) technique. 38 In this work we continue exploring the fragment approach with the periodic HF embedding. We present an interface for the second order perturbation theory based on the restricted-active-space configuration interaction (RASCI): RASCI-PT2 embedded in the periodic Hartree-Fock wavefunction. There are several flavors of the RASCI method; in this work we employ RASCI with one-hole and one-particle approximation. This variant of RASCI is an approximation to CASSCF (within the same active space RAS2), as the orbital relaxation is not performed in RASCI and is only partially taken into account by the single excitations from RAS1 and to RAS3. However, due to its computational simplicity RASCI allows for faster calculations and somewhat bigger active spaces than CASSCF.</p><p>The purpose of this paper is to present the embedded RASCI-PT2 approach and to assess its performance on a particular example of F atom dissociation from a (001) lithium fluoride (LiF) surface and from a partially fluorinated graphane (HF-graphane). In these systems the dynamic correlation is substantial due to the presence of the fluorine atom, and the amount of the static correlation can be tuned by increasing the C−F or Li−F bond distance up to dissociation. In the HF-graphane case, RASCI-PT2 is compared against the distinguishable cluster 39 (DC) approximated coupled cluster with singles doubles and triples (DC-CCSDT) method, 40 which is expected to deliver accurate results for fluorine dissociation.</p><!><p>The embedding scheme used here is related to that of Masur and coworkers, 35 where it was used in conjunction with local ring-CCD. In contrast to that method, here we use local orbitals for the fragment definition only, and once the fragment is defined the orbitals are canonicalized within the fragment. The detailed specification of the current model is to be presented elsewhere 38 and will be not reiterated here. We provide only the general concept, which can be summarized as follows:</p><p>• The orbital space of the fragment is defined by a group of localized occupied and virtual orbitals that are spatially close to the site of interest. This allows for a seamless specification of the fragment and its embedding into the HF mean-field from the rest of the crystal. The occupied orbitals -Wannier functions (WFs) -are mutually orthogonal and correspond to the periodic HF solution. The virtual space is spanned by the so called projected atomic orbitals (PAOs), i.e. AOs projected from the occupied space of the periodic HF solution. In this way the orbital space of the fragment is by construction orthogonal to the HF occupied space of the embedding environment.</p><p>• Once the occupied and virtual spaces of the fragment are defined, the Fock matrix block, corresponding to the fragment orbitals, is cut out from the periodic Fock matrix and diagonalized. The resulting eigenfunctions are the canonical occupied and virtual orbitals of the fragment. These canonical orbitals are then used as the basis for the in-fragment post HF treatment. The active spaces for the RASCI and RASCI-PT2 procedures are chosen on the basis of the eigenvalues of these Fock matrices (i.e. orbital energies).</p><p>• The matrix elements of the fragment's Hamiltonian in the basis of the fragment's canonical orbitals are calculated using the periodic Cryscor code, where the AO support extends into the environment region. The procedure consists of evaluation of the electron repulsion integrals by means of the periodic local density fitting, 41,42 construction of the fragment's core Hamiltonian h frag from the periodic Fock matrix by subtracting the Coulomb and exchange contributions belonging to the fragment. This procedure thus implicitly includes the Coulomb and exchange field of the environment in h frag .</p><p>• The formal nuclear repulsion energy of the fragment E frag nuc is defined such that the fragment's HF energy, computed using h frag and electron-electron repulsion integrals in the basis of the fragments occupied orbitals, coincides with the periodic HF energy. This choice of E frag nuc implies that the whole fragment energy (HF energy plus the correlation correction) can also be interpreted as the energy per cell, where the fragment approximation restricts the correlated orbital space to the fragment orbitals and neglects the correlated couplings beyond the fragment.</p><p>The workflow is outlined in Figure 1. The implementation of this scheme involves three quantum chemistry packages, i.e., Crystal17 , 43 Cryscor , 44 and Q-Chem 45 The final goal is to run the RASCI-PT2 method using a set of orbitals from the converged periodic HF calculation as a basis set, and the Hamiltonian that includes the periodic mean field.</p><p>In the initial stage the periodic restricted closed-shell HF wavefunction is calculated by Crystal17. WFs, which are then constructed by localization of converged HF Bloch functions, 46 as well as the Fock matrix and the relevant auxiliary quantities, like the basis set, geometry, etc, are passed to Cryscor. 44 The localized virtual orbitals, i.e. PAOs, are computed on the Cryscor side, 47 as well as other quantities such as the electron repulsion integrals (ERIs), the fragment's core Hamiltonian h frag and E frag nuc . The ERIs, h frag and E frag nuc together with the number of electrons (for neutral fragments: the half of number of the fragment's WFs) form the essential data needed for the fragment RASCI-PT2 treatment.</p><p>They are passed to the molecular code Q-Chem, which proceeds with the RASCI and RASCI-PT2 calculation.</p><p>Technically, Q-Chem receives the one-h frag pq and two-electron part (pq|rs) of the Hamiltonian in the basis of the fragments orbitals p, q, r and s in the format equivalent to that of Ref. 48. These orbitals take the role of the basis functions for the fragment, and h frag pq and</p><p>(pq|rs) substitute the corresponding AO-based quantities of Q-Chem. Since all the transformations from the AO basis and E frag nuc have already been evaluated on the periodic side, the atoms and the corresponding basis functions themselves are irrelevant for the Q-Chem part. The basis-orbital overlap matrix is replaced with an identity matrix as the fragment canonical orbitals are orthonormal. Q-Chem starts with a HF calculation, which, in case of a singlet state, converges in one iteration, as the fragment's orbitals by construction obey the Brillouin condition for the fragment's Fock matrix. At this point, RASCI and the subsequent second-order perturbative correction are ready to go. The RASCI and RASCI-PT2 approaches of the Q-Chem program can be found in Refs. 49-54. We specify it in some detail in the next subsection.</p><!><p>Generally in the RASCI method, 55 the active orbitals of fragment are divided into three subspaces, i.e. RAS1, RAS2, and RAS3. The RAS2 subspace is similar to the complete active space of the CASSCF or CASCI method and allows all possible configurations of its</p><!><p>Basis set</p><!><p>Fragment seed atoms RASCI parameters is the periodic HF energy; F periodic is the periodic Fock matrix (in the AO basis); φ W F are the occupied WFs; (pq|rs) are the two-electron integrals in the basis of fragment canonical orbitals denoted by the indices p, q, r and s; h frag pq is the matrix of the fragment's one-electron operator, that includes the HF mean field from the embedding; E frag nuc is the formal nuclear energy to equalize the fragment and periodic HF energies; N el -number of electrons in the fragment; N orb -number of orbitals in the fragment; E frag RASCI−PT2 is the target RASCI-PT2 energy of the embedded fragment.</p><p>n electrons in m orbitals. In contrast, in the RAS1 and RAS3 subspaces the occupations are restricted to a maximum number of holes and electrons, respectively. Following the formalism of Ref. 54 we write the excitation operator as R = r0 + rh + rp + rhp + r2h + r2p + ...,</p><p>where r0 generates the full determinantal space within RAS2, while further operators extend this space by certain classes of determinants. The superscripts h and p represent, respectively, hole and electron excitations, such that, for example, rh generate configurations with one hole in RAS1, rp -with one electron in RAS3, etc.</p><p>In this work we employ the hole and particle approximation within RASCI, which truncates the operator (1) to the first three terms, or in other words allows only one hole or one electron in RAS1 or RAS3, respectively. Since there are only single excitations from RAS1 or to RAS3, the size extensivity issues are not as severe as in RASCI (or RASSCF)</p><p>with two holes and electrons in RAS1 and RAS3 spaces or as in MRCI. Furthermore, with the hole and particle approximation it usually becomes unproblematic to extend RAS1 and RAS3 to the full occupied and virtual spaces, respectively. In this form, RASCI can be seen as an approximation to CASSCF, where the single excitations to some extent mimic orbital relaxations.</p><p>The RASCI method with the hole and particle approximation is not able to capture most of the dynamical correlation. With a standard active space RASCI only recovers less than 20 % of the CCSD(T) correlation energy. 54 Therefore for quantitative accuracy RASCI is to be followed by a perturbative treatment. We employ the state-specific multi-reference perturbation theory of Ref. 54 as implemented in the Q-Chem code.</p><p>The method employs the following zeroth order Hamiltonian:</p><p>where |0 and E 0 are, respectively, the Hamiltonian eigenfunction and eigenvalue for the state of interest in the space generated by r0 . In other words, |0 and E 0 correspond to the CASCI solution in the RAS2 space. The configurations |k form the orthogonal complement to |0 and (together with |0 ) span the space generated by the operators of eq. ( 1) that generate at most two holes in RAS1 and two electrons in RAS3. In other words the RASCI-PT2</p><p>wavefunction includes up to the doubly excited configurations with respect to the determinants of RAS2. For the energies E k we use the generalized Davidson-Kapuy partitioning, 56 which represents them as differences of the diagonal elements of a generalized Fock matrix, see equations ( 6) and ( 8) of Ref. 54.</p><p>Since the zeroth-order Hamiltonian corresponds to the CASCI solution in RAS2 rather than the actual RASCI solution, the determinants with a hole in RAS1 or a particle in RAS3 occur among the |k states. In order to avoid double counting, the contributions from these configurations must be excluded from the second order energy correction, which then takes the form:</p><p>where H and P denote these hole and particle spaces, respectively.</p><!><p>As a first test we consider dissociation of a fluorine atom from a surface of an ionic LiF crystal. Dissociation of a LiF molecule is a standard test for a multireference treatment to study the ionic-neutral crossing. To test our embedding method we consider the periodic analog of this system. To model the (001) surface of LiF we used a three-layer slab with a supercell containing 8 LiF units per layer. The model is shown in Figure 2. The structure was relaxed at B3LYP/DZVP level and is explicitly given in Supporting Information (see Figure S2). For the embedded RASCI and RASCI-PT2 calculations we considered a fragment with one fluorine atom surrounded by five lithium atoms, as also shown in Figure 2. the LiF molecule the CASSCF prediction for the crossing can noticeably deviate from that of CASPT2. However, the CASSCF results depend strongly on the choice of the active space, and for an active space comparable to that used here, the differences between CASSCF and CASPT2 crossing positions are also quite small. 57 We finally note, that compared to RASCI, RASCI-PT2 predicts lower excitation energies at dissociation, but higher excitation energies near the minimum.</p><!><p>As a second example we consider covalent-bond breaking. A single bond dissociation is a standard test to assess an ability of a method to accurately describe a multireference case.</p><p>When a simple (e.g. hydrogen-hydrogen) covalent bond is stretched, the orbital energies of the occupied bonding HF orbital σ and the virtual anti-bonding HF orbital (σ * ) are getting close. Therefore, in order to describe the ground (singlet) state of such a system even qualitatively correctly, at least two configurations are needed in the wavefunction expansion.</p><p>In this work we explore a complicated case: dissociation of the C-F bond in a partially fluorinated graphane. Firstly, this is a periodic system and the simple bonding-antibonding orbital picture is not directly applicable, as the canonical orbitals are delocalized Bloch waves. Secondly, the ground state of an isolated fluorine atom is triply degenerate, so at dissociation several orbitals, not just two, become close in the energies. Finally, since the number of electrons in the vicinity of the bond is large (fluorine alone has six lone-pair valence electrons), the dynamic correlation becomes absolutely essential for quantitative accuracy.</p><p>Specification of the model parameters The model is illustrated in Figure 4. The structure of periodic partially fluorinated graphane was optimized at the B3LYP/6-21G* level. The geometry can be found in the Supporting Information in the form of Crystal17 input (see Figure S1). We considered several fragments with progressively increasing size: from 8 to 56 electrons. Figure 4 gives a graphical representation of the fragments by the nuclei denoted with larger spheres. The orbital space of each fragment includes PAOs, centered on these atoms, and WFs, which are either bonding or lone-pair WFs corresponding to these atoms. A detailed specification of the fragments is given in Table 1.</p><p>Table 1: Specification of the fragments. N el is the number of electrons in the fragment (twice the number of WFs), N orb is the number of orbitals in the fragment (both occupied and virtual). We also list the active spaces RAS2 that we employ for different fragments. The RASCI and RASCI-PT2 calculations were performed according to the above de-scribed scheme that involves three codes: Crystal17, Cryscor and Q-Chem. The initial periodic HF calculation employed the 6-21G* basis set. Throughout the periodic stages of the calculations (HF, orbital localization, construction of PAOs) the 8x8 k-grids were used.</p><!><p>For the fluorine dissociation, only the C-F distance was varied, the rest of the structure was kept frozen.</p><p>For the reference values we performed DC-CCSDT calculations for the first three fragments using an interface to the Molpro 58 and GeCCo 59 codes. For the smallest fragment, we also did MP2, CCSD and CCSDTQ, the first two with Molpro and the latter with the MRCC code, 60 interfaced with Molpro. In all these calculations the necessary quantities for high-level methods were also generated by Cryscor.</p><p>Similar to other multireference methods, RASCI is not black-box, as one needs to choose the active space. In our case this concerns RAS2, as RAS1 and RAS3 are then defined automatically to include the complete closed-shell and virtual spaces, respectively. RASCI is expected to be more sensitive to the choice of the active space than CASSCF, since it does not reoptimize the HF orbitals, which might be far from optimal in the multireference regime. In this work we test two types of active spaces, type 1: with plenty of electrons, but with just one empty orbital to excite to; and type 2: a standard active space with a moderate number of electrons in approximately the same number of orbitals. The exact specification of the active spaces for each fragment is also given in Table 1.</p><!><p>We start with the smallest fragment I, consisting of just two atoms, F and C, corresponding to the dissociating bond. The potential energy curves for this fragment are presented in Figure 5. As expected the standard single reference methods (here HF, MP2, CCSD) do not dissociate the fluorine atom even qualitatively correctly. The coupled cluster method with up to quadruple excitations (CCSDTQ) is able to dissociate single bonds and, moreover, accurately captures the dynamic correlation. It could be used as a reference, but it is very computationally expensive and is feasible only for the smallest fragment. Distinguishable cluster method with up to triples excitations DC-CCSDT is known to perform well for fluorine dissociation, 40 and indeed it deviates from CCSDTQ only in the sub-millihartree region along the whole curve (see supplementary information). We thus use DC-CCSDT as the reference for RASCI-PT2. The multireference treatments RASCI and RASCI-PT2 provide a qualitatively correct dissociation for this fragment regardless of the active space choice. RASCI is clearly missing a big share of dynamic correlation, but RASCI-PT2 recovers it virtually completely for both active space choices (and even slightly overcorrelates). A similar picture is observed for the fragment II, presented in Figure 6. However, one can notice here that the RASCI-PT2 curves are already slightly above the DC-CCSDT one. In this fragment the deviations of the RASCI-PT2 total energies from those of DC-CCSDT, which are shown in Figure S3, are around 2-20 millihartree, which is higher than in fragment I. One of the reasons for that is a general loss of size extensivity in a RASCI (or CASCI/CASSCF) energy if the active space is not expanded proportionally to the system size. We note, however, that in RASCI-PT2 this effect is not as pronounced as in RASCI itself. Besides, the the errors in the dissociation energies (Figure S4), which are much more informative than the total energies, are similar to those of fragment I. As is shown in Figure 7, for the fragment III the RASCI-PT2 potential energy curve, corresponding to the "standard" active space (10,12), exhibits some non-physical features.</p><p>We were not able to get rid of these artifacts by varying the active space by an orbital or two. By visualizing the active orbitals (shown in Figures S7 and S8) we found both types of active space to be large enough to include all the relevant orbitals, i.e. sigma bonding and antibonding, p x and p y of the F atom. Therefore, the active space is adequate and absence of orbital rotations in RASCI cannot be the origin of this issue.</p><p>To investigate this problem further, we considered a molecule C 4 H 9 F as finite cluster prototype of Fragment III. The molecule was cut out from the periodic structure with the dangling bonds terminated by hydrogen atoms. The results are shown in Figure S6. Similar to the fragment calculations, RASCI (10,12)PT2 manifests here non-physical artifacts.</p><p>Switching to the spin-flip version of the RASCI-PT2 method, [61][62][63] which uses the triplet state as the reference, does not solve the problem. At the same time, both CASPT2 (10,9)   and CASCI(10,9)PT2 (the latter using the same RHF orbitals as in RASCI(10,9)PT2), calculated with Molpro, provide physically correct curves. This suggests that the problems in the RASCI (10,12)PT2 calculations are either caused by the RASCI-PT2 method itself or its implementation in Q-Chem.</p><p>Importantly, the active space of type 1 - (22,12) -does not cause such problems. Indeed, the RASCI (22,12)PT2 curve is smooth and parallel to the three available DC-CCSDT data points (for the fragment of this size DC-CCSDT in its present implementation is already quite expensive). As expected the deviations of the RASCI (22,12)PT2 total energies from DC-CCSDT grow with fragment size and are noticeably bigger than for fragment II (see Figure S3). But the error in the dissociation energies (Figure S4) remains in the millihartree region and is insensitive to the fragment size.</p><!><p>For chemistry, the total energies themselves are of limited interest, as any chemical process is governed by energy differences. In this study the relevant energy difference is the dissociation energy, which we define here as the difference between the total energies at R C−F = 3.38 Å and R C−F = 1.38 Å. The results for the dissociation energy computed on different fragments are given in Figure 8. The DC-CCSDT method is again used as a reference.</p><p>RASCI, despite the qualitatively correct description, quantitatively grossly underestimates the dissociation energy, especially so with the small active space (type 1). The The active spaces of type 1 and 2 are marked with one and two stars, respectively. The illustration and specification of the fragment size can be found in Figure 4 and Table 1, respectively.</p><p>RASCI-PT2 dissociation energies are much more accurate. They agree with the reference values within a few kcal/mol, which is a standard uncertainty expected from a second-order perturbative treatment. Furthermore, the RASCI-PT2 results are less sensitive to the choice of the active space than RASCI, until it runs into the unphysical regime in fragment III with the active space of type 2. RASCI-PT2 with the active space of type 1, on the other hand, behaves robustly and can be used to explore the convergence with fragment size. Fragments IV and V are too big for DC-CCSDT in its present implementation, but RASCI-PT2 demonstrates that after fragment III the dissociation energy is essentially converged, at least within a few kcal/mol uncertainty.</p><p>Finally we investigate the dependence of the results on the supercell chosen. The supercell size is an important parameter of the model. Indeed a defect, in the case of dissociated fluorine atom, is periodically repeated with the period defined by the supercell. The fragment is defined to cover the defect in the reference cell, but the periodic replicas of the defect can still influence the fragment via the embedding. In our model the embedding remains at the uncorrelated HF level, which can be insufficient for an adequate description of the defect replicas. 64 One way to cope with this problem would be to increase the embedding level or to couple the correlated replicas of the fragments. 65,66 Yet, for extended periodic systems this is difficult. A much simpler procedure would be to "dilute" the defect in the crystal by increasing the supercell size until convergence is reached.</p><p>In this work in addition to the initial unit cell, we consider two supercells: 3x3 and 5x5 (both treated with a single Γ k-point in the Brillouin zone). The RASCI and RASCI-PT2 dissociation energies for fragments I-III and type 1 active space are summarized in the Figure 9 (and Table S6). There is dependence on the supercell size and with the choice of the initial unit cell, the influence of the replicas of the defect is non-negligible. However, for RASCI-PT2 the difference between 3x3 and 5x5 dissociation energies is already not greater than 2 kcal/mol, which can be considered converged. The convergence of RASCI is somewhat slower, but for the fragment III it also seems to be converged already with the 3x3 supercell.</p><!><p>In this work we presented an embedding scheme for investigation of local multireference features inside crystalline systems. The method allows for RASCI-PT2 treatment of fragments embedded in the periodic Hartree Fock wavefunction. Fragments are defined by groups of localized occupied and virtual orbitals, corresponding to the periodic Hartree-Fock solution.</p><p>These orbitals serve as a basis set for the fragment-based multireference treatment. The embedding field, consisting of the Coulomb potential from the crystalline nuclei and the HF mean-field Coulomb and exchange potentials from the electrons outside the fragment, is effectively included in the one-electron Hamiltonian of the fragment. An embedded RASCI-PT2 treatment, based on RASCI with hole and particle approximation, captures both dynamic and static correlation in the region surrounding the site of interest and can deliver quantitatively accurate results at a relatively low for a multireference method computational cost.</p><p>This makes it possible to apply this technique to rather extended fragments and probe the convergence of the results with the fragment size.</p><p>The method has been tested on a fluorine atom dissociation from a surface of LiF crystal and a partially fluorinated graphane. In the former case the qualitative picture of avoided crossing between the ionic and neutral state has been reproduced. For the latter system, a simple choice of the active space that includes most of valence electrons but just one virtual orbital to excite to is sufficient to reach a good accuracy for the dissociation energy with a deviation from DC-CCSDT of just a few kcal/mol. The convergence of the dissociation energy with fragment size was achieved at a 14-atom fragment (fragment III) and 3x3 supercell.</p><p>To estimate the computational cost of this technique: a single point embedded RASCI-PT2 calculation for the fragment III with the (22,12) active space on a single thread Intel of Xeon E5-2640 v4 @ 2.40GHz took about 6.5 minutes. We believe that this method will be a reliable and computationally affordable tool to study multireference problems in solids.</p>