Chinese Journal of Chemical Physics  2019, Vol. 32 Issue (6): 701-707

#### The article information

Wen-yan Zhang, Fei-wu Chen

Iterative Multireference Configuration Interaction

Chinese Journal of Chemical Physics, 2019, 32(6): 701-707

http://dx.doi.org/10.1063/1674-0068/cjcp1905094

### Article history

Accepted on: July 22, 2019
Iterative Multireference Configuration Interaction
Wen-yan Zhang , Fei-wu Chen
Dated: Received on May 15, 2019; Accepted on July 22, 2019
Department of Chemistry and Chemical Engineering, School of Chemistry and Biological Engineering, University of Science and Technology Beijing, Beijing 100083, China; Beijing Key Laboratory for Science and Application of Functional Molecular and Crystalline Materials, Beijing 100083, China
Abstract: Iterative multireference configuration interaction (IMRCI) is proposed. It is exploited to compute the electronic energies of H$_2$O and CH$_2$ (singlet and triplet states) at equilibrium and non-equilibrium geometries. The potential energy curves of H$_2$O, CH$_2$ (singlet and triplet states) and N$_2$ have also been calculated with IMRCI as well as the Møller Plesset perturbation theory (MP2, MP3, and MP4), the coupled cluster method with single and double substitutions (CCSD), and CCSD with perturbative triples correction (CCSD(T)). These calculations demonstrate that IMRCI results are independent of the initial guess of configuration functions in the reference space and converge quickly to the results of the full configuration interaction. The IMRCI errors relative to the full configuration interaction results are at the order of magnitude of 10$^{-5}$ hartree within just 2-4 iterations. Further, IMRCI provides an efficient way to find on the potential energy surface the leading electron configurations which, as correct reference states, will be very helpful for the single-reference and multireference theoretical models to obtain accurate results.
Key words: Multireference configuration interaction    Configuration interaction    Full configuration interaction    Iterative multireference configuration interaction
Ⅰ. INTRODUCTION

Single reference wavefunction theories beyond Hartree-Fock approximation such as single reference perturbation theory [1-11] and coupled cluster theory [12-15] provide accurate ways to calculate the molecular ground state energy as well as the molecular properties [16-18] at the geometry near equilibrium. However, these single reference theories will break down if the molecule of interest is at quasi-degenerate excited states or far from its equilibrium geometry such as cases of the formation or breaking of chemical bonds. For these situations, the wavefunction used to describe the molecular electronic structure must have at least some multireference characteristics. Multireference configuration interaction theories (MRCI) [19-34] are developed along this line.

All-electron full configuration interaction (FCI) is size extensive and orbital-invariant [35, 36]. Its computational results are often regarded as a calibration for other wavefunction models. However, due to its large computational cost, its application is limited to small molecules. Werner and Knowles proposed internal contracted MRCI to reduce the computation cost with the accuracy loss almost negligible [21]. Feller proposed an iterative approach to estimate the FCI energies at complete basis set limit in 1993 [22]. Recently Liu and Hoffmann [26, 27] discussed possible ways to obtain the electronic correlation energy in strong correlated systems and proposed a novel iterative configuration interaction approach (iCI) whose energy converges quickly and monotonically from above to the FCI. Liu and Hoffmann [27] also discussed the relationship of iCI with the ICI (iterative configuration interaction) theory proposed by Nakatsuji et al. [37-39]. The heat-bath CI method of Umrigar and coworkers [31-33] is also similar to these approaches. There are two key stages in heat-bath CI method: generating the variational wave function and energy, and computing the perturbative energy correction [31-33].

Despite the huge progress made in multireference theories, one often ignored problem is the correct guess of initial reference states, which is in general unknown in advance except the reference state at equilibrium. A correct guess of the reference states is essential for the multireferene theoretical models such as MRCI to get results with a high accuracy. Due to this reason and also encouraged by the promising results presented by Liu and Hoffmann [26, 27], an iterative multireference configuration interaction (IMRCI) with a different updating procedure is proposed in this work. Because of huge computational cost of FCI, the purpose here is to investigate how to achieve certain target accuracy such as 10$^{-5}$ hartree with a finite computational cost. The effectiveness of IMRCI has been tested with H$_2$O, CH$_2$ (singlet and triplet states) and N$_2$. Comparison has also been made with various theoretical approaches to show their relative performances. It is also shown through these examples that IMRCI provides an efficient way to find, on the potential energy surface, the leading electron configurations which will be the correct reference states for both single reference and multireference theoretical models.

Ⅱ. THEORY

Multireference configuration interaction starts with a set of $p$ reference functions expressed as follows:

 $\begin{eqnarray} \{\Phi_1, \Phi_2, \Phi_3, \ldots, \Phi_p\}\nonumber \end{eqnarray}$

These reference functions may be a combination of Slater determinants with the expected electronic multiplicity or the so-called configuration state functions. They are supposed to have important contributions to the true wavefunction of an atomic or molecular system under consideration. All possible single and double excitations from these reference functions are carried out to construct the rest of the configuration functions. The wavefunction of the system is supposed to be a linear combination of these configuration functions, the combination coefficients and the corresponding electronic energy are then determined variationally. This is the so-called MRCI single and doubles. For simplicity, it is abbreviated as MRCI.

The effectiveness of MRCI depends on the choice of $p$ reference functions. Usually there are three schemes to achieve this goal: (Ⅰ) first perform a small truncated CI calculation such as the CI with single and double substitutions with respect to some leading configuration function like the electronic ground state, then choose $p$ configuration functions with the largest absolute expansion coefficients as reference functions [22, 39]; (Ⅱ) the preference functions are chosen according to their energy contribution determined by the single reference second order perturbation theory [26, 27, 41]; (Ⅲ) The $p$ reference functions are constructed with $N$ electrons distributed among $M$ active orbitals in all possible ways (the complete active space, thereafter denoted as CAS($N$, $M$)). MRCI based on the three above schemes works well for molecules at the geometry near equilibrium. However, it may fail for molecules at geometries far from equilibrium. In this case the leading configuration function might be some excited state which is in general unknown. Therefore, the schemes (Ⅰ) and (Ⅱ) may lead to wrong reference functions based on the wrong initial guess of the leading configuration function. In principle, scheme (Ⅲ) may produce the expected leading configuration function if the complete active space is large enough. However, this is not an efficient way since the dimension of the complete active space will become very large as $N$ and $M$ increase and the corresponding MRCI calculations will be very expensive.

In order to improve the efficiency and accuracy, IMRCI is proposed. The iteration procedure is expressed in three steps as follows:

(ⅰ) Choose an initial set of preference functions with scheme (Ⅰ):

 $\begin{eqnarray} \{\Phi_1^{(i)}, \Phi_2^{(i)}, \ldots, \Phi_p^{(i)}\}\nonumber \end{eqnarray}$

where $i$ denotes the $i$-th iteration.

(ⅱ) Perform an MRCI calculation and determine expansion coefficients of configuration functions variationally.

(ⅲ) Select a new set of $p$ reference functions with the largest absolute expansion coefficients. If this set of reference functions are the same as the previous ones, the iteration will be stopped. Otherwise, return to the second step and start a further iteration until the $p$ reference functions remain unchanged.

If a better accuracy is required, a relatively larger size of the reference space should be considered and the procedure starting from the step (ⅰ) to (ⅲ) has to be repeated until the self-consistency about the reference functions is reached.

Ⅲ. RESULTS AND DISCUSSION

In this preliminary test, H$_2$O, CH$_2$, and N$_2$ are selected to investigate the efficiency and convergence of IMRCI. Restricted Hartree-Fock orbitals and restricted open-shell Hartree-Fock orbitals were calculated with GAMESS [42]. MRCI, IMRCI, and FCI were performed with MELD [43] and programs in our group. MP2, MP3, MP4, CCSD, and CCSD(T) were calculated with Gaussian [44]. The 1s core orbitals of C, O, and N are frozen for all calculations presented here.

A. H$_{2}$O

Three types of H$_2$O geometries are used in IMRCI and MRCI calculations. One is the experimental geometry with the bond length 1.809 Å and bond angle 104.5$^\circ$ [45]. In order to make the wavefunction of H$_2$O have some multireference characteristics, the OH bond length in the second geometry is stretched to be 2 and 3.3 times the experimental bond length while the bond angle remains unchanged. For convenience of discussion below, these geometries are noted as geometry Ⅰ, Ⅱ, and Ⅲ, respectively. The basis set is Pople's 6-31G [46].

Computational results of MRCI and IMRCI for H$_2$O in geometry Ⅰ are listed in Table Ⅰ. The numbers of electron configurations ($N_{\rm{EC}}$) in the reference space instead of the number of reference functions are listed in the first column. An electron configuration is a distribution of 10 electrons of H$_2$O among molecular orbitals. The number of electron configurations will not be bigger than the number of corresponding reference functions since one electron configuration can be used to construct more than one reference functions for singlet and triplet states studied here if the number of unpaired electrons in the configuration is larger than two. The electron configurations in the reference space are chosen with larger absolute expansion coefficients determined by CISD. The energy errors of MRCI relative to the FCI result ($-$76.119919 hartree) are listed in the second column. As can be seen from Table Ⅰ, the errors remain almost unchanged at the order of 10$^{-4}$ hartree as the size of electron configurations goes from 60 to 200. This demonstrates the slow convergence nature of traditional configuration interaction method. The errors of IMRCI relative to the FCI result are listed in the third column of Table Ⅰ. Because the leading configuration determined by CISD is correct for H$_2$O at equilibrium, the accuracy of IMRCI is almost the same as that of MRCI if $N_{\rm{EC}}$ is smaller than 120. As shown in the bottom of Table Ⅰ, once $N_{\rm{EC}}$ is bigger than 120, the errors of IMRCI become smaller than the errors of MRCI because some relative important electron configurations are selected into the reference space during the iterative procedures in IMCRI. These errors are at the order of magnitude of 10$^{-5}$ hartree. The iteration number of IMRCI ($N_{\rm{IMRCI}}$) is 1 to 3 as shown in the fourth column of Table Ⅰ. For the purpose of comparison, the scheme (Ⅲ) is also used to select the reference space. The numbers of electron configurations in CAS(4, 4) and CAS(6, 6) are 19 and 141, respectively. The corresponding MRCI errors ($\Delta$MRCI/CAS(4, 4) and $\Delta$MRCI/CAS(6, 6)) relative to FCI are 0.005377 and 0.004284 hartree respectively. They are even larger than that of MRCI with $N_{\rm{EC}}$=5 as shown in the first row in Table Ⅰ.

Table Ⅰ Relative energy errors of H$_2$O calculated with MRCI and IMRCI at the experimental equilibrium geometry (geometry Ⅰ). The basis set is 6-31G. FCI energy is $-$76.119919 hartree. $\Delta$MRCI/CAS(4, 4)=0.005377 hartree, $\Delta$MRCI/CAS(6, 6)=0.004284 hartree.

The results of H$_2$O in geometry Ⅱ calculated with MRCI and IMRCI are presented in Table Ⅱ. In comparison with the results in Table Ⅰ, it is easily found that the convergence of MRCI in Table Ⅱ is very slow. The energy error relative to FCI is still at the order of magnitude 10$^{-3}$ hartree although $N_{\rm{EC}}$ is as high as 100. On the other hand, the results of IMRCI converge very quickly to the FCI energy and the error relative to FCI is at the order of magnitude of 10$^{-5}$ hartree even if $N_{\rm{EC}}$ is just 80. The reason for this is that, because of the multireference nature of the wavefunction of H$_2$O in geometry Ⅱ, the less important electron configurations in geometry Ⅰ become more significant in geometry Ⅱ and are selected into the reference space in the IMCRI even if $N_{\rm{EC}}$ is still as small as 40.

Table Ⅱ Relative energy errors of H$_2$O calculated with MRCI and IMRCI in geometry Ⅱ. The OH bond is stretched to be two times of the experimental equilibrium bond length while the bond angle remains unchanged. The basis set is 6-31G. FCI energy is $-$75.879641 hartree. $\Delta$MRCI/CAS(4, 4)=0.004039 hartree, $\Delta$MRCI/CAS(6, 6)=0.003443 hartree.

The relative errors of MRCI with CAS(4, 4) and CAS(6, 6), i.e. $\Delta$MRCI/CAS(4, 4) and $\Delta$MRCI/ CAS(6, 6), are 0.004039 and 0.003443 hartree, respectively. The accuracies are not good. Because of the completeness of CAS, many less important electron configurations are included in the reference space. Thus, even the sizes of CAS(4, 4) and CAS(6, 6) are not small, the corresponding MRCI errors relative to FCI are still not satisfactory. For this reason, MRCI calculations based on the complete active space will not be discussed further although its results will be still listed in the tables below for the purpose of comparison.

The results of H$_2$O in geometry Ⅲ calculated with MRCI and IMRCI are presented in Table Ⅲ. Two basis sets, 6-31G [46] and cc-pVDZ [47] are considered for comparison. The relative errors of MRCI are still at the order of 10$^{-3}$ hartree as the same as shown in Table Ⅰ and . On the other hand, the results of IMRCI converge faster than those shown in Tables Ⅰ and . As the basis set changes from 6-31G to cc-pVDZ, the convergences of MRCI and IMRCI become slightly slow, which is consistent with the previous theoretical results reported by Feller [22]. As can be seen from Table Ⅲ, the accuracy of IMRCI/cc-pVDZ is comparable to the accuracy of Liu and Hoffmann's results with the same number of iterations although no data about the size of electron configurations were given in their original work [27].

Table Ⅲ Relative energy errors (hartree) of H$_2$O calculated with MRCI and IMRCI. The bond length is 3.3 times longer than the experimental bond length while the bond angle remains unchanged. The basis sets are 6-31G and cc-pVDZ, FCI energies for 6-31G and cc-pVDZ are $-$75.835791 and $-$75.90944379 hartree respectively. $\Delta$MRCI/6-31G and $\Delta$MRCI/cc-pVDZ for CAS(4, 4) are respectively 0.004525 and 0.007841 hartree, while $\Delta$MRCI/6-31G and $\Delta$MRCI/cc-pVDZ for CAS(6, 6) are respectively 0.003889 and 0.006461 hartree.

In comparison with the results listed in Tables Ⅰ and , the accuracy between the results of MRCI and IMRCI in Table Ⅲ is very large. It is found that the leading electron configurations in MRCI calculations are the same, i.e., 10 electrons of H$_2$O occupy the five lowest orbitals. On the other hand, the leading electron configurations in reference spaces are changed in IMRCI calculations as $N_{\rm{EC}}$ increases. They are listed in Table Ⅳ for both 6-31G and cc-pVDZ. When $N_{\rm{EC}}$ is bigger than 40, the leading electron configurations remain the same as those configurations with $N_{\rm{EC}}$=40 and therefore they are not listed in the table. In the electron configuration expressions, the number "2", "1", and "0" mean that the orbital is doubly occupied, singly occupied, and empty, respectively. As can be seen from Table Ⅳ, when $N_{\rm{EC}}$ increases to 20 and goes beyond, the leading configuration of the electronic wavefunction is changed from the equilibrium ground state configuration "2222200000" to an excited state "2221111000". This is why IMRCI results converge fast and MRCI results converge slowly. The data in Table Ⅳ also indicate that the leading electron configurations will be independent of the basis set once the size of the multireferece space is large enough such as $N_{\rm{EC}}$=20.

Table Ⅳ The two leading configurations calculated with IMRCI/6-31G and cc-pVDZ for H$_2$O. The bond length is 3.3 times longer than the experimental bond length while the bond angle remains unchanged as the one at equilibrium.The leading configuration of MRCI is 2222200000.

Since the change of the leading configuration on the potential energy surface is in general unknown in advance, the reference state of the single reference models such as MP2, MP3, CCSD and CCSD(T) on the potential energy surface is always the leading electron configuration of the ground state at equilibrium. Therefore, once the leading electron configuration on the potential energy surface was changed, there would be a corresponding jump on the potential energy surface computed with MP2, MP3, CCSD and CCSD(T). The potential energy curves of H$_2$O calculated with MP2, MP3, CCSD CCSD(T), MRCI and IMRCI are shown in FIG. 1. The basis set is 6-31G. The multireference space for MRCI and IMRCI is 30. As can be seen from the figure, due to the change of the leading electron configuration near $R/R_0$=3.3 ($R$ is the OH bond length while $R_0$ is the bond length at equilibrium), there are discontinuities on the potential energy curves of MP2, MP3, CCSD, CCSD(T), and MRCI. However, the potential energy curve of IMRCI is very smooth. FIG. 1 demonstrates the important role of the leading electron configuration on the potential energy surface. IMRCI not only computes correctly the potential energy surface but also provides an efficient way to find the leading electron configuration.

 FIG. 1 The potential energy curves of H$_2$O calculated with various theoretical models. The size of the multireference space is 30. $R$ is the OH bond length while $R_0$ is the bond length at equilibrium. The bond angle remains the same as the angle at equilibrium.
B. Singlet and triplet states of CH2

Methylene is a well-known model molecule tested by various theoretical approaches [9, 28, 48]. The geometries of the singlet and triplet states of CH$_2$ are the same as those in the Taylor's work [48]. The basis set is Pople's 6-31G [46]. The computational results of the singlet and triplet states of CH$_2$ are presented in Table Ⅴ. It is well-known that the singlet and triplet states have different correlation effect. It can be seen easily from the data in Table Ⅴ that the singlet state has more mulireference characteristics than the triplet state. The convergence behaviors of MRCI and IMRCI are very similar for singlet state when $N_{\rm{EC}}$ is relatively small. However, as shown in the table, the convergence difference becomes significant as $N_{\rm{EC}}$ goes beyond 60. Errors of IMRCI results relative to the FCI energy are about one order of magnitude smaller than those of MRCI results. On the other hand, the convergence rates of MRCI and IMRCI for triplet state are almost the same as $N_{\rm{EC}}$ increases from 5 to 100, which to some extent indicates the single reference nature of the triplet state.

Table Ⅴ Relative energy errors (hartree) of the singlet and triplet states of CH$_2$ calculated with MRCI and IMRCI. FCI energies of the singlet and triplet states are $-$38.942887 and $-$38.980200 hartree, respectively.

The potential energy curves of the singlet and triplet states of CH$_2$ are presented in FIG. 2 and FIG. 3. The size of the multireference space is 30. These two figures demonstrate the similar convergence behaviors of MRCI and IMRCI as shown in Table Ⅴ. As for the singlet state, although higher accuracy results can be achieved with IMRCI if $N_{\rm{EC}}$ increases to 70 or an even larger number, yet the leading electron configuration will remain unchanged. This is the big difference between the potential energy curves of H$_2$O and CH$_2$. Nevertheless, the performances of MP2, MP3, CCSD, and CCSD(T) in the two figures are still not satisfactory.

 FIG. 2 The potential energy curves of the singlet state of CH$_2$ calculated with various theoretical methods. The size of the multireference space is 30. $R$ is the CH bond length while $R_0$ is the bond length at equilibrium. The bond angle remains the same as the angle at equilibrium.
 FIG. 3 The potential energy curves of the triplet state of CH$_2$ calculated with various theoretical methods. The size of the multireference space is 30. $R$ is the CH bond length while $R_0$ is the bond length at equilibrium. The bond angle remains the same as at equilibrium.
C. Potential energy curves of N$_2$

N$_2$ is a typical molecule with a triple bond. It is a challenging problem for theoretical approaches to obtain a correct description of the dissociation of N$-$N bond. The potential energy curves calculated with MP2, MP3, MP4, MRCI, and IMRCI are presented in FIG. 4. The basis set is cc-pVDZ [47]. Because of the failure of iterative procedure in CCSD at $R$=5.135 Bohr, the potential energy curves of CCSD and CCSD(T) are not plotted in the figure. The curves calculated with MP2 and MP4 bend downward to the horizontal $R$ axis while the MP3 curve bends in the opposite direction. Due to the discontinuity of its potential energy curve at $R$=6.635 Bohr, the performance of MRCI is also inferior to that of IMRCI.

 FIG. 4 The potential energy curve of N$_2$ calculated with various theoretical methods. $R$ is the N$-$N bond length. The size of the multireference space is 40.
Ⅳ. CONCLUSION

MRCI is a well-established and reliable configuration interaction approach applicable for the electronic ground state as well as the excited states. It is especially useful for molecules at equilibrium although its slow convergence to the FCI results is often observed. The MRCI accuracy depends strongly on the initial guess of the configuration functions in the reference space. A wrong initial guess which often occurs for molecules far from its equilibrium geometry will lead to very low MRCI accuracy. Similar to iCI [26, 27] and ICI [37-39], preliminary test results presented here demonstrated that IMRCI can improve the accuracy of MRCI for molecules both at equilibrium and non-equlibrium geometries with a few iterations due to the fact that the most important electronic configurations can be selected into the reference space through the iterative procedure presented here. It is often thought that the problems of the single reference theoretical models on the potential energy surfaces can be solved with the multireference theoretical methods. However, this is not true since the multireference models are strongly dependent on the guess of initial multireference wavefunctions which, on the other hand, is in general unknown in advance. As we discussed for H$_2$O molecule, the leading electron configuration will not remain the same and will change along the potential energy surface. IMRCI provides an efficient way to find the leading electron configurations on the potential energy surface. This information about the reference states will be very helpful for multirefence perturbation theory and multireference coupled cluster theory to obtain results with a higher accuracy.

Ⅴ. ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China (No.21473008 and No.21873011).

Reference
 [1] C. Møller, and M. S. Plesset, Phys. Rev. 46, 618(1934). DOI:10.1103/PhysRev.46.618 [2] C. W. Murray, and E. R. Davidson, Chem. Phys. Lett. 187, 451(1991). DOI:10.1016/0009-2614(91)80281-2 [3] R. D. Amos, J. S. Andrews, N. C. Handy, and P. J. Knowles, Chem. Phys. Lett. 185, 256(1991). DOI:10.1016/S0009-2614(91)85057-4 [4] P. J. Knowles, J. S. Andrews, R. D. Amos, N. C. Handy, and J. A. Pople, Chem. Phys. Lett. 186, 130(1991). DOI:10.1016/S0009-2614(91)85118-G [5] W. J. Lauderdale, J. F. Stanton, J. Gauss, J. D. Watts, and R. J. Bartlett, Chem. Phys. Lett. 187, 21(1991). DOI:10.1016/0009-2614(91)90478-R [6] D. Jayatilaka, and T. J. Lee, Chem. Phys. Lett. 199, 211(1992). DOI:10.1016/0009-2614(92)80108-N [7] T. J. Lee, and D. Jayatilaka, Chem. Phys. Lett. 201, 1(1993). DOI:10.1016/0009-2614(93)85024-I [8] P. M. Kozlowski, and E. R. Davidson, Chem. Phys. Lett. 226, 440(1994). DOI:10.1016/0009-2614(94)00763-2 [9] F. Chen, J. Chem. Theory. Comput. 5, 931(2009). DOI:10.1021/ct800546g [10] F. Chen, M. Wei, and W. Liu, China Chem. Sci. 54, 446(2011). DOI:10.1007/s11426-010-4199-1 [11] S. R. Yost, and M. Head-Gordon, J. Chem. Phys. 145, 054105(2016). DOI:10.1063/1.4959794 [12] J. ${\rm{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over C} }}$í${\rm{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over z} }}$ek, J. Chem. Phys. 45, 4256(1966). DOI:10.1063/1.1727484 [13] G. D. Purvis, and R. J. Bartlett, J. Chem. Phys. 76, 1910(1982). DOI:10.1063/1.443164 [14] K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon and, Chem. Phys. Lett. 157, 479(1989). DOI:10.1016/S0009-2614(89)87395-6 [15] D. I. Lyakh, M. Musiaz, V. F. Lotrich, and R. J. Bartlett, Chem. Rev. 112, 182(2012). DOI:10.1021/cr2001417 [16] T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen, and K. Ruud, Chem. Rev. 112, 543(2012). DOI:10.1021/cr2002239 [17] J. Cao, M. Wei, and F. Chen, Acta Phys. Chim. Sin. 32, 1639(2016). [18] Y. C. Chou, Comput. Theor. Chem. 1111, 50(2017). DOI:10.1016/j.comptc.2017.04.003 [19] I. Shavitt, The Method of Configuration Interaction, H. F. Schaefer Eds., Methods of Electronic Structure Theory, Modern Theoretical Chemistry, Boston: Springer, Vol. 3, 189-275 (1997). [20] B. Brooks, and H. F. Schaefer, J. Chem. Phys. 70, 5092(1970). [21] H. J. Werner, and P. J. Knowles, J. Chem. Phys. 89, 5803(1988). DOI:10.1063/1.455556 [22] D. Feller, J. Chem. Phys. 98, 7059(1993). DOI:10.1063/1.464749 [23] C. D. Sherrill, and H. F. Schaefer, Adv. Quantum. Chem. 34, 143(1999). DOI:10.1016/S0065-3276(08)60532-8 [24] P. G. Szalay, T. Müller, G. Gidofalvi, H. Lischka, and R. Shepard, Chem. Rev. 112, 108(2012). DOI:10.1021/cr200137a [25] X. Chen, Z. Chen, and W. Wu, J. Chem. Phys. 141, 194113(2014). DOI:10.1063/1.4901729 [26] W. Liu, and M. R. Hoffmann, Theor. Chem. Acc. 133, 1481(2014). DOI:10.1007/s00214-014-1481-x [27] W. Liu, and M. R. Hoffmann, J. Chem. Theory. Comput. 12, 1169(2016). DOI:10.1021/acs.jctc.5b01099 [28] J. Cao, and F. Chen, Acta Phys. Chim. Sin. 33, 1130(2017). [29] J. B. Schriber, and F. A. Evangelista, J. Chem. Phys. 144, 161106(2016). DOI:10.1063/1.4948308 [30] J. B. Schriber, and F. A. Evangelista, J. Chem. Theory. Comput. 13, 5354(2017). DOI:10.1021/acs.jctc.7b00725 [31] A. A. Holmes, N. M. Tubman, and C. J. Umrigar, J. Chem. Theory Comput. 12, 3674(2016). DOI:10.1021/acs.jctc.6b00407 [32] S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, J. Chem. Theory. Comput. 13, 1595(2017). DOI:10.1021/acs.jctc.6b01028 [33] A. D. Chien, A. A. Holmes, M. Otten, C. J. Umrigar, S. Sharma, and P. M. Zimmerman, J. Phys. Chem. A 122, 2714(2018). DOI:10.1021/acs.jpca.8b01554 [34] H. Lischka, D. Nachtigallová, A. J. A. Aquino, P. G. Szalay, F. Plasser, F. B. C. Machado, and M. Barbatti, Chem. Rev. 118, 7293(2018). DOI:10.1021/acs.chemrev.8b00244 [35] Z. Rolik, A. Szabados, and P. R. Surján, J. Chem. Phys. 128, 144101(2008). DOI:10.1063/1.2839304 [36] W. A. Vigor, J. S. Spencer, M. J. Bearpark, and A. J. W. Thom, J. Chem. Phys. 144, 094110(2016). DOI:10.1063/1.4943113 [37] H. Nakatsuji, and M. Ehara, J. Chem. Phys. 122, 194108(2005). DOI:10.1063/1.1898207 [38] H. Nakatsuji, J. Chem. Phys. 113, 2949(2000). DOI:10.1063/1.1287275 [39] H. Nakatsuji, and M. Ehara, J. Chem. Phys. 117, 9(2002). DOI:10.1063/1.1487830 [40] J. Ivanic, and K. Ruedenberg, Theor. Chem. Acc. 106, 339(2001). DOI:10.1007/s002140100285 [41] S. Evangelisti, J. P. Daudey, and J. P. Malrieu, Chem. Phys. 75, 91(1983). DOI:10.1016/0301-0104(83)85011-3 [42] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Natsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347(1993). DOI:10.1002/jcc.540141112 [43] L. E. McMurchie, S. T. Elbert, S. R. Langhoff, and E. R. Davidson, MELD, Modifications by D. Feller and D. C. Rawlings, http://php.indiana.edu/~davidson/meld.htm. [44] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V.Barone, B. Mennucci, G. A.Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E.Peralta, F.Ogliaro, M.Bearpark, J. J. Heyd, E. Brothers, K. N.Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C.Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. Cross, V. Bakken, C.Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian 09, Revision D.01, Wallingford CT: Gaussian Inc., (2009). [45] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, New York: Macmillan, (1982). [46] W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. 56, 2257(1972). DOI:10.1063/1.1677527 [47] T. H. Dunning Jr., J. Chem. Phys. 90, 1007(1989). DOI:10.1063/1.456153 [48] C. W. Bauschlisher, and R. P. Taylor, J. Chem. Phys. 85, 6510(1986). DOI:10.1063/1.451431