Calculation of Low-Lying Electronic Excitations of Magnesium Monofluoride: How Well Do Coupled-Cluster Methods Work?

1. Introduction

The spectroscopy of polar diatomic molecules, such as magnesium monofluoride (MgF) has a long history, and has led in the late 1960ies to results with remarkable accuracy. Barrow and Beale [1] resolved the rotational structure of gaseous MgF and carried out a detailed analysis resulting in the determination of molecular constants for the ground and first three electronic excited states with precise values based on the two lowest vibrational states ($v=0,1$). Higher excitations were measured and reported in 1971 by Novikov and Gurvich[2]. The detailed structure of the ${\mathrm{X}}^2\sum \to {\mathrm{A}}^2\prod$ transitions including spin-orbit effects (the ${\mathrm{A}}^2\prod$ doublet is split by about 35 ${\mathrm{cm}}^{-1}$) focused on the nature of the $\Lambda$-type doubling, and it was concluded that the doublet was inverted (negative A parameter). A recent study [3] does, however, conclude that the ${\mathrm{A}}^2\prod$ doublet is normal, i.e., not inverted. Hyperfine-resolved measurements are also reported in Ref. [4].

The interest in laser cooling of MgF has resulted in a detailed study of potential energy curves and spectroscopic constants at the level of multi-reference configuration interaction (MRCI) [5,6]. These results can serve the present study as a theoretical benchmark, in addition to the experimental results [7]. The ground-state potential energy curve can be determined also from infrared spectroscopy [8] which leads to an extended Morse oscillator description of this curve. This allows us to focus directly on the calculation of electronic excitations. These studies (and the present work) ignore spin-orbit effects. Such effects can be calculated in some standard packages, e.g., using two-component density functional theory [9] as implemented in TURBOMOLE [10].

Our interest in studying the accuracy of excited-state calculations within coupled-cluster (CC) theory relates to experimental work on matrix isolated barium monofluoride (BaF) [11]. In this context earlier work on matrix isolated atomic barium [12,13] was extended to the case of BaF in argon [14,15] and neon [16]. The matrix isolation work is motivated by performing measurements of the electron electric dipole moment. While the BaF molecule is much more complicated than MgF, the use of effective core potentials makes the computational effort comparable for the two cases. The purpose of the present work is to establish expectations on accuracy with respect to methodology and basis sets in order to understand how reliable the CC methods may be in the matrix isolation environment.

The current status of describing the electronic excitations of gas-phase BaF (and its lighter homologs CaF and SrF) can be described as follows. Results from scalar relativistic CC and MRCI calculations are reported in Ref. [17]. The all-electron Fock-space CC results for the excitation energies (which include spin-orbit splittings on the order of 400-600 ${\mathrm{cm}}^{-1}$) agree with experiment at the level of about 200 ${\mathrm{cm}}^{-1}$. Better agreement (down to about 10 ${\mathrm{cm}}^{-1}$) was obtained in Ref. [18] by including QED effects and basis set extrapolation of correlation effects calculated at the level of CCSDT(Q). A similar approach based on four-component CC Dirac theory resulted in accurate ionization energies [19]. Hyperfine constants for ${\mathrm{F}}^{19}$ ground and excited-state BaF molecules were obtained in Ref. [20] and provide some confidence in the quality of the Fock-space CC wavefunctions.

To return to the topic of the present paper we note that in Ref. [21] radiative decay rates and branching fractions of the first excited ${\mathrm{A}}^2\prod$ state of ²⁴${\mathrm{Mg}}^{19}\mathrm{F}$ were reported in a comparative experimental-theoretical work. While the theory is the relativistic Fock-space CC method, comparison can be made with the present non-relativistic work due to the smallness of spin-orbit effects in MgF. The experimental data presented in Ref. [21] deviate somewhat (at the $15{\mathrm{cm}}^{-1}$ level) from the earlier work [1]. Studies which investigate equation-of-motion (EOM) coupled cluster methods have been carried out and are typically benchmarked against full configuration interaction calculations, but the goal of this work is to only compare the theoretical $T_e$ spectroscopic parameters to experimental determinations.

The layout of the paper is as follows. We begin in Section 2 with a short summary of how to use CC methods for electronic excitations. In Section 3 we present our results for excitation energies of MgF using three methods and two basis set families, and compare them with experimentally determined values, and with the relativistic Fock-space CC results of Ref. [20]. We draw some conclusions in Section 4.

2. Theory

Electronic excitation energies for molecules can be calculated in many quantum chemistry packages with a selection of methods. For the purpose of the present work the accuracy of time-dependent density functional theory methods is not sufficient, since we would like to reach percent-level agreement with measured values or better. This requires the inclusion of electron correlation with a systematic treatment as provided by configuration interaction or coupled cluster methodology.

A fast method in this realm is the CC2 method coupled with response theory [22]. It is, e.g., implemented in TURBOMOLE [10], alongside with a faster version based on Møller-Plesset perturbation theory, namely the ADC(2) method [23]. These allow for fast calculations based on the resolution-of-identity methods, and efficient basis sets of triple- and quadruple-zeta quality [24], (e.g, def2-QZVPPD). Transition moments and other excited-state first-order properties are calculated quickly even for larger molecules [25]. We will refer to this method as CC2-RPA. It is mostly restricted to not use the full symmetry, i.e., it is running with C1 symmetry by default in TURBOMOLE. One can use the properties of the excited states, such as transition moments to get some idea of which excited states have been generated.

If one requires higher accuracy, one has available a general excited-state method at the level of CC theory including single and double excitations, named EOM-CCSD (EOM stands for equation of motion) [26,27]. Since the class of polar molecules to which MgF belongs is dominated by single excitations, it should provide reasonable accuracy. One can ask for a number of excitations for each irreducible representation. We used Psi4 [28] to perform the calculations. The CCEOM package in Psi4 also allows to perform higher-accuracy calculations at the level of EOM-CC3 theory. These are demanding calculations, but are considered the best there is as a general-purpose CC method for excited states [29]. The EOM-CC3 method (as implemented in Psi4) first calculates the CC3 ground state, then performs an EOM-CCSD for multiple roots (which in our work has the trend of being higher in energy than the usual EOM-CCSD roots), and then one has to pick one root at a time to calculate the EOM-CC3 value. Thus, the method is more cumbersome to use. Faster implementations of EOM-CC3 than generally available have been developed [30].

In addition to the mentioned methods one can also attempt so-called $\Delta$-CC methods [31]. One complements a standard ground-state CC calculation with a run where the SCF calculation is carried our for an excited state obtained by switching two orbitals (e.g., highest occupied and lowest unoccupied) to seed an SCF calculation for a Slater determinant which represents the excited state. This step may or may not work, as the SCF calculation may collapse back to the ground state. When it works, the resulting orbitals can be used in a correlated calculation at the CCSD or CCSD(T) level or higher. For the MgF molecule we were not able to carry the $\Delta$-CC approach, since the $\Delta$$-SCF$ step failed, but it is interesting to note that in cases where the method works the excitation energies from $\Delta$-CCSD and EOM-CCSD in general do not agree with each other [32]. The agreement may become better at the level of the CC3 model. Such results for the barium monofluoride (BaF) molecule using a pseudopotential for the inner electrons of barium for the ${\mathrm{A}}^{\prime }$$\Delta$, $A\prod$, and $B\sum$ states will be reported in a separate publication.

In terms of basis sets we used the def2-family at triple- and quadruple-zeta level and made a guess at the complete basis set (CBS) limit by using an extrapolation of the form

$$E_{\mathrm{cbs}}=\frac{3^3{E}_{\mathrm{TZ}}-4^3{E}_{\mathrm{QZ}}}{3^3-4^3}\approx 1.73E_{\mathrm{QZ}}-0.73E_{\mathrm{TZ}},$$

(1)

where $E_{\mathrm{TZ}}$ and $E_{\mathrm{QZ}}$ are the energies obtained with def2-TZVPPD and def2-QZVPPD respectively.

We also used the correlation-consistent augmented basis sets aug-cc-PVnZ with $\mathrm{n}=3,4,5$ [33,34]. From the calculated energies $E_{\mathrm{n}}$ the CBS limit can be obtained using

$$E_{\mathrm{cbs}}\approx 1.583E_5-0.597E_4+0.004E_3$$

(2)

which follows from Equation (2) in Ref. [35].

For the MgF molecule the known bond lengths for the ground and excited states in question are very close ($R_0\approx 1.75$Å). Thus, we restrict our calculations here to a single Mg-F separation, since we are interested in comparing values of $T_e$, i.e., for excitation energies the difference in the minimal energies between the potential energy curves. Detailed calculations of such curves at the level of four-component (relativistic) Fock space CC theory are given, e.g., in Figure 1 of Ref. [21] and can be used for justification of this shortcut. For other molecules, such as, e.g., CaF, SrF and BaF one has to calculate the potential energy curves, i.e., energies as a function of R, e.g., as given in Ref. [17] in order to determine accurate adiabatic excitation energies. From these curves one can then determine the spectroscopic parameters ${\omega }_e$ and ${\omega }_e{x}_e$, which follow from the nuclear vibrational excitation energies in accord with

$$E_v=T_e+{\omega }_e(v+\frac{1}{2})+{\omega }_e{x}_e{(v+\frac{1}{2})}^2.$$

(3)

In the present work, however, the focus is entirely on the determination of accurate values of $T_e$ for a few excited states relative to the ${\mathrm{X}}^2\sum$ ground state.

A final general comment for this section can be made: each of the theory models and for each of the chosen basis sets the ground-state is calculated at some level of precision. All CC calculations start from a common self-consistent field result, which for our calculations is based on spin-unrestricted Hartree-Fock theory. The changes in ground-state energy with basis set at this level are not too significant. The differences between these values from these calculations are appreciable, and they arise mainly for two reasons: (i) the correlation contributions included in CC2, CCSD, and CC3 are different, and (ii) the truncated basis set introduces a significant error when computing correlation effects. Excitation energies are calculated relative to ground states which are numerically quite different. Nevertheless, they can be compared directly to each other when looking for excited-state $T_e$ values.

3. Results

3.1. Mg(${\mathrm{X}}^2\sum \to {\mathrm{A}}^2\prod$)

The first excitation is to the $A^2{\prod }_J$ states. The present non-relativistic calculation can be compared to the experimental analysis of Ref. [1] which lists not $T_e$ but the transition energies to the lowest two vibrational states $v=0,1$. The $J=\left(\frac{1}{2}\right),\left(\frac{3}{2}\right)$ doublet has a separation of $\approx 35-36{\mathrm{cm}}^{-1}$ which will be of the order of the accuracy goal for calculating this excitation. Figure 1 shows this doublet as a pair of solid lines, and the relativistic calculation from Ref. [21] as a pair of dashed lines at the right. This calculation gives results which are higher by about $60{\mathrm{cm}}^{-1}$, but it should be noted that Ref. [21] includes an analysis of the data of Ref. [7] which results in higher values for the excitation energies than the original data of Barrow and Beale [1] by about $17{\mathrm{cm}}^{-1}$. Such a difference exceeds the accuracy of the original experimental analysis, and can be related to the use of the effective Hamiltonian involved in the analysis of the data. An example of such a discrepancy was reported in Ref. [36]. For our purposes one may then argue that the best nonrelativistic result would fall not in the middle between the two black lines, but closer to the shown ${\mathrm{A}}^2{\prod }_{3/2}$ energy.

On the left in Figure 1 we observe that the extrapolated CC2-RPA results underestimate the transition energy: the extrapolated result using the def2 basis sets (blue square) disagrees by about $120{\mathrm{cm}}^{-1}$. The aug-cc basis set calculations fare a bit better, reducing the discrepancy to about $80{\mathrm{cm}}^{-1}$. The open symbols show that the results from the EOM-CCSD method overestimate the excitation energy. The relationship between the two basis set families remains similar to the results from the CC2-RPA method. Finally, the EOM-CC3 results (solid circles) turn out to be lower than EOM-CCSD, with the aug-cc basis set results leading to an extrapolated value that falls within $25{\mathrm{cm}}^{-1}$ of the results from Ref. [1], or Ref. [7] as re-analyzed in Ref. [21]. The EOM-CC3 method with the aug-cc basis set family can be considered a great success in reaching better than per mille relative accuracy for the transition energy. It would be interesting to find out whether the relativistic Fock-space CC method of Ref. [21] can be pushed to higher order to reach this level of precision on a system which involves 21 electrons.

3.2. Mg(${\mathrm{X}}^2\sum \to {\mathrm{B}}^2\sum$)

Excitations to higher electronic states can be expected to come out with less accuracy, particularly since they are in the same symmetry group. Figure 2 shows that this is clearly the case: the ordinate spans a wider range. All methods approach the experimental excitation energy from above, and exhibit qualitatively a similar trend when comparing the def2 and aug-cc basis sets. The differences are larger (on the order of $300{\mathrm{cm}}^{-1}$), and the aug-cc basis set calculations clearly do better. Remarkably the EOM-CC3 result for the aug-cc basis set falls right onto the experimental value when extrapolated to the CBS limit. The Fock-space CC results of Ref. [21] are more than $400{\mathrm{cm}}^{-1}$ too high.

3.3. Mg(${\mathrm{X}}^2\sum \to {\mathrm{C}}^2\sum$)

The next-higher state was considered by Barrow and Beale [1] to be a Rydberg state. The energy scale in Figure 3 is expanded to cover $2000{\mathrm{cm}}^{-1}$. Interestingly, this range is required mostly to show the difference between results from the two basis set families. The EOM-CC3 method gives an excellent result in the CBS limit. In order to assess whether this is a lucky coincidence, we report in Table 1 the actual EOM-CC3 data for the aug-cc basis sets. These data allow one to observe some interesting trends, particularly as to how the extrapolation works.

For the excitation to the ${\mathrm{A}}^2\prod$ the progression of numbers approaches the CBS from below in a systematic way. The extrapolation to a value which is about $100{\mathrm{cm}}^{-1}$ above the $\mathrm{n}=5$ answer is based on the fact that the quadruple- and quintuple-zeta basis result differ by almost $200{\mathrm{cm}}^{-1}$.

For the ${\mathrm{B}}^2\sum$ excitation the progression of $\mathrm{n}=3,4,5$ results is not monotonic. It follows, however, from Equation (2) that the estimated CBS limit result is governed by the increase in excitation energy between the $\mathrm{n}=4,5$ results. The variation of the excitation energy with basis set quality is at the $200{\mathrm{cm}}^{-1}$ level. This is an interesting observation in the context of the relativistic calculation of Ref. [21] which overestimates experiment by about twice this amount. It leads to the conclusion that the EOM-CC3 method includes the important correlation contributions almost independent of the value of n when using the aug-cc-pVnZ basis. Given that the def2 basis result was less successful (cf. Figure 2), one can argue that the augmentation of the basis set is important to achieve the reported accuracy.

For the case of ${\mathrm{C}}^2\sum$ excitation the sequence of excitation energies for $\mathrm{n}=3,4,5$ is systematic and again approaches the extrapolated result from below. The data for all three excited states presented in Table 1 support the notion that the triple-zeta (n=3) calculations can be used to obtain a quick overview of the situation accurate at the few-hundred ${\mathrm{cm}}^{-1}$ level.

4. Conclusions

Electronically excited states of magnesium fluoride, a basic polar diatomic molecule were calculated at few levels of coupled-cluster theory. The very fast CC2-RPA method as implemented in TURBOMOLE yields good first estimates of the excitation spectrum, particularly in combination with the aug-cc basis set family.

The EOM-CCSD method which is a turnkey method to obtain the spectrum tends to overestimate the excitation energies. While the def2 basis works reasonably well for the lower two excitations, it leads to overestimated energies for the ${\mathrm{C}}^2\sum$ state both in the CC2-RPA and the EOM-CCSD methods.

The EOM-CC3 method requires more effort (both operationally and in terms of computing time), but represents the most accurate of these toolsets. For molecules where the $\Delta$$-CC$ methodology works, i.e., higher roots than for the ground state can be found $\Delta$$-CCSD\left(T\right)$ computations should represent an economical alternative to it.