Extent of Spin Contamination Errors in DFT/Plane-wave Calculation of Surfaces: A Case of Au Atom Aggregation on a MgO Surface

Abstract

The aggregation of Au atoms onto a Au dimer (Au₂) on a MgO (001) surface was calculated by restricted (spin-un-polarized) and unrestricted (spin-polarized) density functional theory calculations with a plane-wave basis and the approximate spin projection (AP) method. The unrestricted calculations included spin contamination errors of 0.0–0.1 eV, and the errors were removed using the AP method. The potential energy curves for the aggregation reaction estimated by the restricted and unrestricted calculations were different owing to the estimation of the open-shell structure by the unrestricted calculations. These results show the importance of the open-shell structure and correction of the spin contamination error for the calculation of small-cluster-aggregations and molecule dimerization on surfaces.

1. Introduction

Dimerization and its reverse process, homolytic cleavage, are basic chemical reactions, which are important for investigating heterogeneous reactions on surfaces [1,2,3,4,5,6,7,8,9,10,11,12]. For instance, dissociation of hydrogen on metal surfaces is an elementary reaction for hydrogenation by heterogeneous metal catalysts, and dimerization of nitrogen atoms and nitric monoxide (NO) molecules are regarded as elemental reactions for NO_x elimination from exhaust gases. In addition, both dissociation and dimerization of metal clusters on their support surfaces are crucial for clarifying the degradation mechanism of nanosized metal catalysts, especially gold catalysts [10,13,14,15,16,17,18].

We cannot estimate the potential energy curves for the dimerization of atoms (or/and radical molecules) and hemolytic cleavage of molecules by first-principle calculations with spin-unpolarized (spin-restricted) single reference schemes [19]. When the monomer has unpaired electrons, the reaction systems have an open-shell structure; however, spin-unpolarized calculations cannot estimate the open-shell structure, and the estimated values include artificial unstabilization. The artificial error is called static correlation. To avoiding this error, we perform multi-reference calculations or spin-polarized (spin-unrestricted) calculations [19,20,21,22,23,24,25,26,27].

Thus, multi-reference or spin-polarized schemes are adopted for the first-principle calculations of the potential energy curves for dimerization. The multi-reference calculations are generally more accurate than spin-polarized calculations because the methods account for the multifigurational states, but the computational costs are very high [20,21]. In contrast, the computational costs of spin-polarized calculations are not very high, so these methods can be applied to surface reaction systems, but the methods include an artificial error, such as spin contamination error [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. However, it was shown that relevant corrections result in the same potential energy curves as those by multireference calculations even when spin-polarized calculations are used [31,32,33,34].

For first-principle calculations of heterogeneous catalysts, large slab surface models, large metal cluster models, or both are needed. The models consist of hundreds to thousands of atoms, and multireference calculations cannot be adopted due to the high computational costs. Therefore, we adopt the spin-polarized single reference scheme for calculating the adsorptions of molecules that have open-shell structures. Generally, the calculations are based on the density functional theory (DFT) [37], the core electrons are treated by pseudo potentials [38,39,40], and the wave functions are expanded by plane-wave basis sets (DFT/plane-wave method).

The effects of spin contamination error and static correlation on dimerizations (homolytic cleavages) in the gas phase are well known, and the errors can be corrected. By contrast, the errors that are included in DFT/plane-wave calculations for dimerizations (hemolytic cleavages) on surfaces have not been corrected, even though it is possible that the errors affect the results of DFT/plane-wave calculations. This is due to the difficulty of the corrections.

Calculating dimerizations on a surface by DFT/plane-wave, large supercells are necessary for ignoring artificial interactions between a given adsorbate and the adsorbate in the neighboring cell. Hence, it was difficult to perform even spin-polarized calculations owing to the high computational costs. However, with the availability of high-performance of computers, high-cost spin-polarized DFT/plane-wave calculations have become feasible; i.e., nowadays, we can estimate the effects of static correlation on the DFT/plane-wave calculations for surface reactions.

On the other hand, there was no scheme that can estimate the errors in DFT/plane-wave calculations, for correcting spin contamination errors. However, it was shown that the approximate spin projection (AP) scheme [28,41] for the correction of spin contamination errors can be applied to the DFT/plane-wave method (AP-DFT/plane-wave) [35,36]. By using AP-DFT/plane-wave method, we can investigate the effects of spin contamination on surface reactions. Recently, it was clarified that (1) the spin contamination error affects the interaction energy between Au atoms dissociating on a MgO (001) surface [35], (2) the rate-determining step for NO dimerization on a ZrO₂/Cu catalyst model is changed when the error is corrected [36], and (3) the error occurs in NO dimerization on a TiO₂/Ag catalyst model even when the NO-NO distance is outside the region where the spin contamination error does not occur in the gas phase [36].

Thus, it has become possible to investigate the effects of spin contamination and static correlation on DFT/plane-wave calculations for describing the potential energy surfaces of surface reactions, although detailed investigations are yet to be performed. The present study is the first to investigate the aforementioned effects. The investigated surface reaction is the aggregation of Au atoms on a MgO (001) surface (Au₂ dissociation on the MgO (001) surface) shown in Equation (1):

$$(\mathrm{Au}+\mathrm{Au})/\mathrm{MgO}\rightleftharpoons {\mathrm{Au}}_2/\mathrm{MgO}$$

(1)

The Au/MgO system is a typical gold catalyst system [14,16,35,42,43,44], and reaction (1) is related to the sintering of Au clusters on the MgO surface, which pertains to a crucial degradation of the gold catalyst [13,14,15,16,17,18], and its suppression. The potential energy curves of the closed singlet (LS: low spin), triplet (HS: high spin), and open singlet (BS: broken symmetry, which is a common name of open-shell structures estimated by spin-polarized calculations) states were estimated, and the spin contamination error in the open singlet state was estimated by the AP scheme (AP: BS state with AP correction). In addition, the corresponding gas-phase reaction (Equation (2)) was calculated in the same manner as the surface reaction (Equation (1)):

$$\mathrm{Au}+\mathrm{Au}\rightleftharpoons {\mathrm{Au}}_2$$

(2)

Using the calculation results, we investigated the effects of the errors on reaction (1) and the effects of molecule/surface interactions on the magnitude of the errors.

2. Computational Procedure

2.1. Method

The wave functions were calculated by the DFT scheme [37] with the PW91 exchange-correlation functional [45]. Spin-polarized DFT was adopted for calculating the open-shell structures (HS and BS), and spin-unpolarized DFT was adopted for the closed-shell structure (LS). For the interpolation of the correlation part, the formula proposed by Vosko, Wilk, and Nusair was used [46]. The core electrons were treated by the ultra-soft-pseudo-potential method [38], and the numbers of valence electrons of Mg, O, and Au were 2, 6, and 11, respectively. The energy cut-off values were 400 eV (wave function) and 2400 eV (augmented charge). The only Γ-point was calculated; the energy dependence on k-point sampling has already been investigated in our prior work [35]. The DFT calculations were performed by the Vienna Ab initio Simulation Package (VASP) [47,48,49,50]. For population analysis, the Bader charge analysis [51,52,53,54] was used. The Visualization for Electronic and STructural Analysis program (VESTA) [55] was used for visualization of the geometry and spin density distribution.

2.2. AP Scheme

The spin contamination errors in the energies of the open singlet states were estimated by the AP scheme [28,41]. In the AP method, the spin contamination error is corrected via a projection to the Heisenberg Hamiltonian:

$${\widehat{\mathbf{H}}}_{\mathrm{HB}}=-2J_{\mathrm{ab}}{\widehat{\mathbf{S}}}_{\mathrm{a}}{\widehat{\mathbf{S}}}_{\mathrm{b}}$$

(3)

${\widehat{\mathbf{S}}}_{\mathrm{x}}$ is a spin operator on site x, and $J_{\mathrm{ab}}$ is the effective magnetic exchange interaction between spins. The projection on the Heisenberg Hamiltonian indicates that Equations (4) and (5) are established:

$$J_{\mathrm{ab}}=\frac{E_{\mathrm{HB}}^{\mathrm{LS}}-E_{\mathrm{HB}}^{\mathrm{HS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{HB}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{HB}}^{\mathrm{LS}}}=\frac{E_{\mathrm{SA}}^{\mathrm{LS}}-E_{\mathrm{SA}}^{\mathrm{HS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{LS}}}$$

(4)

$$J_{\mathrm{ab}}=\frac{E_{\mathrm{HB}}^{\mathrm{LS}}-E_{\mathrm{HB}}^{\mathrm{HS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{HB}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{HB}}^{\mathrm{LS}}}=\frac{E_{\mathrm{U}}^{\mathrm{LS}}-E_{\mathrm{U}}^{\mathrm{HS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{LS}}}$$

(5)

Here, E and $\langle {\widehat{\mathbf{S}}}^2\rangle$ are the expectation values of the Hamiltonian $\langle \mathsf{\Psi }\left|\widehat{\mathbf{H}}\right|\mathsf{\Psi }\rangle$ and the square of the total spin operator $\langle \mathsf{\Psi }|{\widehat{\mathbf{S}}}^2|\mathsf{\Psi }\rangle$, respectively. The superscripts LS and HS indicate the values of the low spin state and high spin state, respectively. The subscripts (HB, SA, exact, and U) show the method for estimating the expectation values; namely, HB is the value of the Heisenberg model Hamiltonian, SA is estimated by highly accurate spin-adapted calculations such as multireference calculations, exact means the exact value (the eigenvalue of ${\widehat{\mathbf{S}}}^2$), and U indicates the value estimated by spin-polarized (unrestricted) calculations.

Because the spin contamination error is artificial contamination from the higher spin state [19], the HS state is not affected by the error. Therefore, the HS values of U and SA are approximately the same, and we obtain:

$$E_{\mathrm{SA}}^{\mathrm{LS}}=\frac{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{LS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{LS}}}{E}_{\mathrm{U}}^{\mathrm{LS}}+\frac{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{LS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{LS}}}{E}_{\mathrm{U}}^{\mathrm{HS}}$$

(6)

The spin contamination error in energy is the difference between $E_{\mathrm{SA}}^{\mathrm{LS}}$ and $E_{\mathrm{U}}^{\mathrm{LS}}$:

$$\mathrm{SCE}=\frac{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{LS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{BS}}}{{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{exact}}^{\mathrm{HS}}-{\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}^{\mathrm{BS}}}(E_{\mathrm{U}}^{\mathrm{HS}}-E_{\mathrm{U}}^{\mathrm{BS}})$$

(7)

Here, LS is rewritten as BS, because the low-spin state including the spin contamination error has a broken-symmetry electronic structure. The ${\langle {\widehat{\mathbf{S}}}^2\rangle }_{\mathrm{U}}$ values were estimated by Wang’s scheme [56]:

2.3. Model

A 4 × 4 MgO (001) slab, as shown in Figure 1a, was used for the MgO (001) surface. The lattice constant was optimized by using the MgO bulk structure. The number of atomic layers in the surface model was 4, and the bottom layer was fixed to mimic the MgO bulk structure. The dependencies of the Au atom adsorption energy and O atom desorption energy on the number of atomic layers were checked. Figure 1b shows a schematic view of the calculated reaction, such as aggregation of Au atoms (Au₂ dissociation) on the MgO (001) surface. The diffusion direction of the Au atom is [010], and the position of the diffusing Au atom (²Au) is represented by fractional coordinate r. For the structure shown on the left side of Figure 1b, r = 0.500, and for that on the right side, r = 0.000. By varying r from 0.500 to 0.000 in 0.025 increments, the potential energy of reaction (1) was estimated. The geometry of the MgO slab (except for the bottom layer) and the [001] coordination of the Au atoms were optimized at each r value. Figure 1c–f show the optimized typical structures and their r values.

2.4. Definitions of Energies

The energy baseline of this work is the non-adsorbed state, such as Au + Au + MgO. When using this baseline definition, the energy of the calculated ¹Au²Au/MgO system is the same as the adsorption energy shown in Equation (8):

$$E_{\mathrm{ads}}=E\left({}_{}{}^1\mathrm{A}\mathrm{u}{}_{}{}^2\mathrm{A}\mathrm{u}/\mathrm{MgO}\right)-E\left(\mathrm{MgO}\right)-2E\left(\mathrm{Au}\right)$$

(8)

E(¹Au²Au/MgO) is the total energy of the Au₂ (or ¹Au and ²Au atoms)-adsorbed MgO model, E(MgO) is the total energy of the optimized 4 × 4 MgO (001) slab, and E(Au) is the total energy of the Au atom in the 30 × 30 × 30 Å supercell that is sufficient volume to calculate Au atom energy. By calculating E_ads at each r value, the potential energy is estimated.

E_ads can be divided into four energies; E_int(¹Au/²Au) is the interaction energy between ¹Au (non-diffused Au atom) and ²Au (diffused Au atom) defined by Equation (9); E_int(¹Au/MgO) is the interaction energy between ¹Au and MgO defined by Equation (10); E_int(²Au/MgO) is the interaction energy between ²Au and MgO defined by Equation (11); and E_dis(MgO) is the distortion energy of MgO due to the adsorption of Au atoms defined by Equation (12):

$$E_{\mathrm{int}}({}_{}{}^1\mathrm{A}\mathrm{u}/{}_{}{}^2\mathrm{A}\mathrm{u})=E\left({}_{}{}^1\mathrm{A}\mathrm{u}{}_{}{}^2\mathrm{A}\mathrm{u}/\mathrm{MgO}\right)-E{\left({}_{}{}^1\mathrm{A}\mathrm{u}/\mathrm{MgO}\right)}_{\mathrm{fix}}-E{\left({}_{}{}^2\mathrm{A}\mathrm{u}/\mathrm{MgO}\right)}_{\mathrm{fix}}+E{\left(\mathrm{MgO}\right)}_{\mathrm{fix}}$$

(9)

$$E_{\mathrm{int}}({}_{}{}^1\mathrm{A}\mathrm{u}/\mathrm{MgO})=E{\left({}_{}{}^1\mathrm{A}\mathrm{u}/\mathrm{MgO}\right)}_{\mathrm{fix}}-E{\left(\mathrm{MgO}\right)}_{\mathrm{fix}}-E\left(\mathrm{Au}\right)$$

(10)

$$E_{\mathrm{int}}({}_{}{}^2\mathrm{A}\mathrm{u}/\mathrm{MgO})=E{\left({}_{}{}^2\mathrm{A}\mathrm{u}/\mathrm{MgO}\right)}_{\mathrm{fix}}-E{\left(\mathrm{MgO}\right)}_{\mathrm{fix}}-E\left(\mathrm{Au}\right)$$

(11)

$$E_{\mathrm{dis}}\left(\mathrm{MgO}\right)=E{\left(\mathrm{MgO}\right)}_{\mathrm{fix}}-E\left(\mathrm{MgO}\right)$$

(12)

E(X) is the total energy of X, and the subscript fix in Equations (9)–(12) indicates that the total energy is calculated using the fixed structure which is in the optimized ¹Au²Au/MgO system.

The previous work [35] argued that the percentage of spin contamination error in E_int(¹Au/²Au) is large. Therefore, the divided energies (E_int(¹Au/²Au), E_int(¹Au/MgO), E_int(²Au/MgO), and E_dis(MgO)) were also investigated. In addition, variations in the effective exchange integral J_ab (Equation (5)) and the energy difference BS and HS, ΔE (Equation (13)), were investigated; these values reflect the stability of the BS state:

$$\mathsf{\Delta }E=\left|E\left(\mathrm{HS}\right)-E\left(\mathrm{LS}\right)\right|$$

(13)

3. Results and Discussion

All the optimized structures are included in the Supplementary Materials. Figure 2 shows the distance data for each r. d(¹Au-²Au) is the distance between ¹Au and ²Au, d(¹Au-¹O) is the distance between ¹Au and ¹O, and d(²Au-n.n.) is the distance between ²Au and its nearest neighbor atom. The nearest neighbor atoms are summarized in the Supplementary Materials (Table S1).

Figure 3 depicts the calculated potential energy curves for reaction (1); the horizontal axis in Figure 3a,b is r and d(¹Au-²Au), respectively. The calculated E_ads reflects the total energy; therefore, the stability of the electronic states and adsorption structure can be judged from the value of E_ads.

The dependencies of E_int(¹Au/²Au), E_int(¹Au/MgO), E_int(²Au/MgO), and E_dis(MgO) on r are shown in Figure 4. For E_int(¹Au/²Au), the values of the LS (closed singlet state), HS (triplet state), BS (open singlet state without correction of the spin contamination error), and AP (open singlet state with correction of the spin contamination error by the AP scheme) were estimated. For other energies such as E_int(¹Au/MgO), E_int(²Au/MgO), and E_dis(MgO), only the LS, HS, and BS values were estimated. For estimating these energies, the energy of an open singlet structure, which is a state including the spin contamination error, was not used; therefore, these energies do not include the spin contamination error, and the AP values are identical to the BS value. Figure 5 shows the percentage of the four energies in E_ads.

From the results shown in Figure 2a, it is confirmed that Au atoms aggregates to Au₂ as r decreases. We classified the calculation results with the nearest neighbor atom of ²Au, and made discussions. Section 3.1 shows the results for the interaction of ²Au with ³O (r = 0.500–0.450, typical structure is r = 0.500 shown in Figure 1c). Section 3.2 shows the results for the interaction of ²Au with ²Mg (r = 0.425–0.325, the typical structure is r = 0.375 shown in Figure 1d). Section 3.3 shows the results for the interaction of ²Au with ²O (r = 0.300–0.200, the typical structure is r = 0.250 shown in Figure 1e). Section 3.4 shows the results for the interaction of ²Au with ¹Au or ¹Mg (r = 0.175–0.000, the typical structure is r = 0.000 shown in Figure 1f).

3.1. Dissociated Au₂ Adsorption onto MgO: r = 0.500–0.450

In this r region, r = 0.500–0.450, the Au atoms interact with separate O²⁻ ions. d(¹Au-²Au) is long (Figure 2a), and the E_int(¹Au/²Au) values of the open-shell states (HS and BS) are zero (Figure 4a,b). However, the E_int(¹Au/²Au) values of the closed-shell state (LS) are not zero, and the value for r = 0.500 (d(¹Au-²Au) = 8.4984 Å) is 0.40 eV. It is difficult to assume that a repulsion of 0.40 eV will experimentally occur; therefore, it is reasonable to assume that the unstabilization is an artificial error, such as static correlation. The spin-unpolarized DFT calculations cannot estimate the potential energy curves of dissociation of the covalent bonds in the gas phase due to the static correlation, and the artificial error occurs as Au₂ dissociation in the gas phase [35]. The dominant interaction between the Au atoms in Au₂ is the covalent interaction between the Au 6s orbitals. The highest occupied molecular orbital (HOMO) is a bonding orbital, and the lowest unoccupied molecular orbital (LUMO) is an anti-bonding orbital. The frontier orbitals do not change when Au₂ is adsorbed onto metal oxides such as MgO [35,43], Al₂O₃ [18,57], and TiO₂ [58,59,60,61]. Hence, the HOMO and LUMO of Au₂ can be described by using the 6s orbital of ¹Au (${\phi }_1$) and 6s orbital of ²Au (${\phi }_2$):

$${\phi }_{HOMO}=\mathrm{A}\left({\phi }_1+{\phi }_2\right)$$

(14)

$${\phi }_{LUMO}=\mathrm{B}\left({\phi }_1-{\phi }_2\right)$$

(15)

A and B are linear combination constants. Then, the frontier orbitals of the calculated ¹Au²Au/MgO systems can be simplified:

$${\phi }_{Au-Au}=\mathrm{C}{\phi }_1+\mathrm{D}{\phi }_2$$

(16)

In the gas phase, the two Au atoms are identical: |C| = |D|. In the adsorbate/surface system, the two Au atoms are not identical, and then $|\mathrm{C}|\ne \left|\mathrm{D}\right|$. The spin-unpolarized calculations use the same orbitals for electrons with different spins. Using the scheme, the occupied spin orbitals of Au₂ systems with the 6s orbitals are:

$$\{\begin{array}{c}{\chi }_1=\alpha {\phi }_{Au-Au}\\ {\chi }_2=\beta {\phi }_{Au-Au}\end{array}$$

(17)

and the Slater determinant is:

$$\begin{array}{c}\frac{1}{\sqrt{2}}\left|\begin{array}{cc}{\chi }_1\left(r_1,{\omega }_1\right)& {\chi }_2\left(r_1,{\omega }_1\right)\\ {\chi }_1\left(r_2,{\omega }_2\right)& {\chi }_2\left(r_2,{\omega }_2\right)\end{array}\right|\\ =\frac{\mathrm{CD}}{\sqrt{2}}\left|\begin{array}{cc}\alpha \left({\omega }_1\right){\phi }_1\left(r_1\right)& \beta \left({\omega }_1\right){\phi }_2\left(r_1\right)\\ \alpha \left({\omega }_2\right){\phi }_1\left(r_2\right)& \beta \left({\omega }_2\right){\phi }_2\left(r_2\right)\end{array}\right|+\frac{\mathrm{CD}}{\sqrt{2}}\left|\begin{array}{cc}\alpha \left({\omega }_1\right){\phi }_2\left(r_1\right)& \beta \left({\omega }_1\right){\phi }_1\left(r_1\right)\\ \alpha \left({\omega }_2\right){\phi }_2\left(r_2\right)& \beta \left({\omega }_2\right){\phi }_1\left(r_2\right)\end{array}\right|\\ +\frac{\mathrm{CD}}{\sqrt{2}}\left|\begin{array}{cc}\alpha \left({\omega }_1\right){\phi }_1\left(r_1\right)& \beta \left({\omega }_1\right){\phi }_1\left(r_1\right)\\ \alpha \left({\omega }_2\right){\phi }_1\left(r_2\right)& \beta \left({\omega }_2\right){\phi }_1\left(r_2\right)\end{array}\right|+\frac{\mathrm{CD}}{\sqrt{2}}\left|\begin{array}{cc}\alpha \left({\omega }_1\right){\phi }_2\left(r_1\right)& \beta \left({\omega }_1\right){\phi }_2\left(r_1\right)\\ \alpha \left({\omega }_2\right){\phi }_2\left(r_2\right)& \beta \left({\omega }_2\right){\phi }_2\left(r_2\right)\end{array}\right|\end{array}$$

(18)

Here, $r_x$ and ${\omega }_x$ are the geometry coordination and spin coordination of electron x. The third and fourth terms of Equation (18) are ion terms, whose schematic view is shown in Figure 6. These terms cause static correlation because they are artificially included and unstable in dissociated structures. Even if $|\mathrm{C}|\ne \left|\mathrm{D}\right|$, e.g., in an adsorbed system, the terms do not disappear and affect the calculated energy. In other words, even if we consider a surface reaction system, when the dissociated adsorbates have localized spins, the DFT/plane-wave calculations of the aggregation reactions and dissociation reactions include static correlation, and the calculation fails. This is similar to the failure of spin-unpolarized calculations of the dissociations of diatomic molecules composed of two elements. Thus, the calculation results (Figure 3 and Figure 4) explicitly show the importance of investigating the open-shell structure by DFT/plane-wave calculations of the surface reaction.

Figure 7 summarizes the calculation results for the gas-phase reaction, reaction (2). The calculation results are obtained by removing MgO from the optimized ¹Au²Au/MgO system and calculating the fixed structures. The results of E_int(¹Au/²Au) for the surface reaction (Figure 4a) and those for the gas-phase reaction (Figure 7a) clearly differ; the former are constant, while the latter are not constant. The reason for this trend can also be inferred from Equation (18). Because of the electrostatic interactions with the MgO surface, the energy difference between the ion (third and fourth) and radical (first and second) terms fluctuates.

3.2. Dissociating Au₂ Adsorption onto MgO (1): r = 0.425–0.325

In this r region, the dissociating Au (²Au) interacts with ²Mg. The potential energy surface (Figure 3) shows that the transition state of reaction (1) is included in this r region. The structures in which the ²Au atom is adsorbed around the ²Mg atom are the least stable. However, the heights of the reaction (1) (aggregation of Au atoms) of LS, HS, BS, and AP are different: 0.17, 0.45, 0.40 and 0.32 eV, respectively (

Table 1. Activation energies of reaction (1).

State	Aggregation of Au Atoms	Au2 Dissociation
LS	0.17 eV (0.53 eV 1)	2.47 eV
HS	0.45 eV	0.00 eV (2.39 eV 2)
BS	0.40 eV	2.34 eV
AP	0.32 eV	2.26 eV

¹ Estimated using the BS value as baseline; ² Estimated using the LS value as baseline.

).

The LS states are the least stable among calculated states. Since the LS state at r = 0.500 (the dissociated structure) is also unstable, the activation barrier of the aggregation reaction is estimated to be low (0.17 eV). Assuming that the energy of the open-shell structure is adopted as the baseline (the energy for the dissociated adsorption structure: r = 0.500), the activation barrier is 0.53 eV. The instability is due to the static correlation. However, the effects on the surface reaction at r = 0.325–0.425 (Figure 4a) are smaller than those of the gas-phase reaction (Figure 7a). This is because the ion terms in Equation (18) for the surface system is not as unstable as those of the gas phase. There are electrostatic interactions such as ${}_{}{}^1\mathrm{A}{\mathrm{u}}^{\mathsf{\delta }+}\cdots {}_{}{}^1\mathrm{O}{}^{2-}$ and ${}_{}{}^2\mathrm{A}{\mathrm{u}}^{\delta -}\cdots {}_{}{}^2\mathrm{M}{\mathrm{g}}^{2+}$ in the ¹Au²Au/MgO system. The results of Bader charge analysis (Figure 8a,b) show the charge polarization such as ${\mathrm{Au}}^{\mathsf{\delta }+}\cdots {\mathrm{Au}}^{\delta -}$. The Bader charges also show that the LS state overestimates the polarization as compared with the BS state. The overestimation causes a large distortion of the MgO slab (E_dis(MgO)) and the LS state becomes more unstable.

The second least stable state is the HS state. In the HS state, there is no artificial unstabilization at r = 0.500; therefore, the activation barrier of the aggregation reaction is the highest (0.45 eV). This value is higher than that of AP by 0.1 eV, and it is obvious that the HS calculation misestimates the activation barrier of reaction (1). In addition, due to the lack of ${\mathrm{Au}}^{\mathsf{\delta }+}\cdots {\mathrm{Au}}^{\delta -}$ polarization, the absolute values of E_int(¹Au/²Au) and E_dis(MgO) are smaller than those for the other states.

The BS state is the ground state in the region r = 0.325–0.425, similar to the other r regions. Two interactions stabilize this state: one is the electrostatic interaction between the Au atoms and MgO (similar to that in the LS state at r = 0.325–0.425) shown in Figure 8, and the other is a spin delocalization to the O sites shown in Figure 9 and Figure 10. When the Au atom interacts with the MgO (001) surface, the spin density on the Au atom is delocalized to the interacted O, and Au-O becomes a spin site [35,62,63]. Figure 9 and Figure 10 show that the spin density on the Au atom delocalizes only to the O sites but not to the Mg site; the adsorption of the Au atoms affects the electronic structure of only the O atoms in MgO (001). Thus, the radical spin still remains and localized at the Au-O sites when the Au atoms are adsorbed; therefore, the spin contamination occurs, and the total energy of the BS states is affected by the error.

At this r region, the spin contamination error should be investigated; the results are shown in Figure 11, with related values such as $\langle {\widehat{\mathit{S}}}^2\rangle$, J, and $\mathsf{\Delta }E$. The values of SCE, $\langle {\widehat{\mathit{S}}}^2\rangle$, and J of reactions (1) and (2) are summarized in

Table 2. The values of SCE, $\langle {\widehat{\mathit{S}}}^2\rangle$, and J of BS states for reactions (1) (surface reaction) and (2) (gas-phase reaction).

r	SCE/eV		$\langle {\widehat{\mathit{S}}}^2\rangle$		J/eV		d(1Au-2Au)/Å
	Surface	Gas	Surface	Gas	Surface	Gas
0.000	0.000	0.000	0.0000	0.0000	−1.208 1	−0.894 1	2.52917
0.025	0.000	0.000	0.0000	0.0000	−1.193 1	−0.894 1	2.53112
0.050	0.000	0.000	0.0000	0.0000	−1.148 1	−0.893 1	2.54011
0.075	0.000	0.000	0.0000	0.0000	−1.073 1	−0.891 1	2.55653
0.100	0.000	0.000	0.0000	0.0000	−0.966 1	−0.891 1	2.57930
0.125	0.000	0.000	0.0000	0.0000	−0.823 1	−0.893 1	2.61054
0.150	0.000	0.000	0.0000	0.0000	−0.600 1	−0.880 1	2.66684
0.175	0.000	0.000	0.0000	0.0000	−0.407 1	−0.626 1	2.99674
0.200	0.000	0.000	0.0000	0.0001	−0.207	−0.399 1	3.40169
0.225	−0.092	0.000	0.7747	0.0001	−0.195	−0.207	3.82430
0.250	−0.049	−0.095	0.9567	0.7240	−0.093	−0.083	4.24920
0.275	−0.020	−0.054	1.0168	0.9215	−0.037	−0.031	4.67412
0.300	−0.009	−0.023	1.0394	0.9915	−0.017	−0.012	5.09924
0.325	−0.009	−0.009	1.0394	1.0178	−0.015	−0.004	5.52656
0.350	−0.075	−0.003	0.8399	1.0283	−0.088	−0.001	5.96645
0.375	−0.084	−0.001	0.7947	1.0325	−0.101	0.000	6.39464
0.400	−0.066	0.000	0.8554	1.0344	−0.075	0.000	6.81369
0.425	−0.001	0.000	1.0616	1.0353	−0.001	0.000	7.22546
0.450	0.000	0.000	1.0619	1.0358	−0.001	0.000	7.64864
0.475	0.000	0.000	1.0605	1.0360	0.000	0.000	8.07348
0.500	0.000	0.000	1.0599	1.0360	0.000	0.000	8.49840

¹ Estimated by using the LS value.

.

As shown in Figure 11a and Table 2, the spin contamination errors in this r region are large, and the maximum is 0.084 eV. This result is surprising because the spin contamination error does not occur at the Au-Au distances in the gas phase (Figure 7b and Table 2). The spin contamination error affects the calculated energies of BS states of the surface system even at r = 0.325–0.425 because the two Au atoms are not identical due to the interactions with the MgO surface, and a clear energy difference is generated between HS and BS. This mechanism is confirmed by the results of J and $\mathsf{\Delta }E$, summarized in Table 2.

Figure 3 shows that the spin contamination error also influences the estimation of the activation barrier and the structure of the transition state. The activation barriers without AP correction are 0.40 eV (aggregation reaction of Au atoms) and 2.34 eV (Au₂ dissociation reaction), while those with AP correction are 0.32 eV and 2.26 eV, respectively. In the aggregation reaction, the percentage of the error is large (20%). In addition, the transition states of BS and AP are different. The former is the structure at r = 0.375 (the ²Au atom is at the on-top site of ²Mg), while the latter is the structure at r = 0.400, although the energy difference between the structures at r = 0.375 and 0.400 of AP is small (0.01 eV). Considering the asymmetry of the structure of the surface adsorption system, it is not surprising that the position of the ²Au atom in the transition state shifts from the on-top site of ²Mg. In fact, even before correction, the energies of BS at r = 0.350 and 0.400 are different, and the energy at r = 0.400 is more unstable.

Note that the correction of the spin contamination error in this study was performed only for energy. The corrections to the first derivative (force) and the second derivative (frequency) of energy have not been performed yet, and the activation energies and transition state structures in this study are roughly estimated. Nevertheless, the results of the present work are sufficient to show the importance of the investigation and correction of artificial errors such as spin contamination and static correlation in spin-polarized DFT/plane-wave calculations for surface reactions, because it is obvious that the errors affect the calculated results of the potential energy curves for the aggregation of Au atoms (Au₂ dissociation) on the MgO (001) surface.

3.3. Dissociating Au₂ Adsorption onto MgO (2): r = 0.300–0.200

In the r region of this section (r = 0.300–0.200), when the Au-Au bond is being created, the ground state is BS and there is an energy difference between HS and BS [35], as confirmed by the potential energy curves shown in Figure 3. The ²Au diffusing on MgO interacts with ²O, and there are two spin states (¹Au-¹O and ²Au-²O) in the BS states structures of r = 0.200–0.300 (Figure 9).

At r = 0.200, because the covalent bond between the Au atoms remains, the LS state without the spin is more stable than the HS state with spins at the Au-O sites. The BS state does not have any large localized spin, and the electronic state is almost the same as that of LS state, as shown by the results of the $\langle {\widehat{\mathit{S}}}^2\rangle$ value (Figure 11a), spin density distribution (Figure 9m), and local magnetic moment (Figure 10).

From r = 0.225, the results of the LS, HS, and BS states clearly differ. At r = 0.225–0.300, the ground state is the BS state, and the HS state artificially affects the BS state (spin contamination). The maximum and minimum SCE values of the surface reaction are 0.092 eV (r = 0.225) and 0.009 eV (r = 0.300), respectively. On the other hand, the maximum, minimum, and second minimum SCE values of the gas-phase reaction are 0.0095 eV (r = 0.250), 0.000 eV (r = 0.225), and 0.023 eV (r = 0.300), respectively.

These results indicate that the maximum SCE value of the gas phase is larger than that of the surface reaction. This can be explained by the mechanism reported in our previous study [35]. In the dissociated Au₂ on MgO, the covalent interaction between the Au atoms is weakened by the interaction between the adsorbed Au and the adsorbing O^2-. Then, the open singlet state of the surface system becomes less stable than that of the gas-phase system; the energy difference between the HS and BS states decreases, and the SCE value becomes small. SCE(surf.) is smaller than SCE(gas) at r = 0.250–0.300 because of this destabilization of the BS state by the Au/MgO interaction.

However, at r = 0.225, the value of SCE(surf.) is larger than that of SCE(gas) (SCE(gas) is zero), which is explained by the mechanism reported in our other previous study [36]. The Au-Au interaction is weakened by the Au-MgO interaction, which leads to destabilization of not only the BS state but also the LS state. In the surface system, even when the Au-Au distance is short so that there are covalent interactions in the gas phase, there is no clear covalent interaction between the Au atoms, and the LS state of the surface system becomes less stable than the BS state of the surface reaction system. The destabilization of the LS state is confirmed from the results shown in Figure 3a. Because of this destabilization, the BS state becomes the ground state, the spin is localized at the Au-O sites, as shown by the spin density distribution in Figure 9l, and spin contamination occurs.

From the results described in this Section and Section 3.2, it is concluded that in the surface reaction systems, the region where spin contamination affects the total energies becomes wider than that in the gas-phase reaction systems.

On the other hand, Figure 3a and Figure 4a show that the LS state is the least stable among the calculated states at r = 0.250–0.300 due to the static correlation. The artificial unstabilization shown in equation (18) occurs even in the structures where the Au-Au bond is creating (is breaking) at the surface.

The spin contamination and static correlation affect the calculation energies from r = 0.225. Figure 5 shows that the percentages of the Au/Au interaction (E_int(¹Au/²Au)) and the Au/MgO interaction (E_int(¹Au/MgO) and E_int(²Au/MgO)) change from r = 0.225, i.e., from r = 0.225, the Au/MgO interactions become larger than the Au-Au interaction. This result suggests that the spin contamination and static correlation affect the calculation results when the covalent interaction within the adsorbate weakens and the interaction between the dissociated fragments and the surface becomes dominant. This will be used as a criterion for judging whether the effect of artificial errors on surface reactions should be considered.

3.4. Non-Dissociated Au₂ Adsorption onto MgO: r = 0.175–0.000

In this r region, the distance between ¹Au and ²Au is so short that the ¹Au-²Au covalent interaction is dominant. When there is a clear covalent interaction between ¹Au and ²Au, the open singlet state (BS state) does not exist and converges to the same electronic structure as the closed singlet state (LS state) [35]. Therefore, in this r region, the results of LS and BS identical, and the spin contamination error does not affect the results.

The d(¹Au-²Au) values of HS are longer than those of the other states in the region r = 0.000–0.150. This is because the Au atoms do not interact covalently in the r region, as confirmed by the E_int(¹Au/²Au) results of HS shown in Figure 4a.

The d(¹Au-MgO) values of HS are also longer than those of the other states because of the insignificant interaction between the Au atoms of HS. When a Au cluster is adsorbed onto metal oxides, it acquires a slight positive charge and interacts with an O²⁻ ion electrostatically [16,17,18,35,43,44,58,59,60,61,64,65,66,67]. Furthermore, the charge on the Au cluster consisting of two or more Au atoms is polarized, and the Au atoms closer to the surface are more positively charged [17,44,60,65,67]. In other words, when there are Au-Au covalent interactions, the Au cluster can interact with MgO more strongly than the Au atom. In the HS state, which lacks Au-Au covalent interactions, Au₂ is not polarized and d(¹Au-MgO) becomes longer, as confirmed by the Bader atomic charges shown in Figure 8.

4. Conclusions

In order to clarify the effects of spin contamination error and static correlation, the aggregation reaction of Au atoms (Au₂ dissociation) on the MgO (001) surface have been investigated in detail. The calculation results show that the singlet state calculated by spin-polarized DFT is the ground state throughout the reaction, and that the calculated energy is affected by the spin contamination error. The investigated reaction can be divided into four steps:

(1)

Two isolated Au atoms are adsorbed onto the on-top of O site in MgO (001), which is the initial state of the aggregation of Au atoms

(2)

One Au atom interacts with Mg and the other Au atom interacts with O, and there is a slight electrostatic interaction between the Au atoms, which includes the transition state of the Au atom aggregation

(3)

Generation of Au-Au covalent interaction

(4)

Non-dissociated adsorption state of Au₂ on MgO, which is the final structure of the Au atom aggregation.

The spin-unpolarized calculation that gives only the closed-shell structure (LS) can calculate step (4) accurately, but the calculation fails at other steps due to the static correlation. Using the results of the HS states, we can estimate the energy of step (1). However, the HS state calculation misestimates the Au-Au interaction and cannot treat steps (2–4). The BS calculation results, which are the ground state of all the steps, are influenced by the spin contamination error, except for steps (1) and (4). Therefore, the BS calculations overestimate the activation barrier and underestimate the Au-Au interaction energy; for treating all the steps, AP calculation is needed. These results explicitly show the importance of the open-shell structure calculations (spin-polarized DFT/plane-wave) with a spin contamination error correction for investigating the aggregation of small metal clusters on metal oxide surfaces.

It is shown that “the maximum value of the spin contamination error in surface reactions is smaller than that in gas-phase reactions, whereas the area causing spin contamination in the surface reaction is wider than that in the gas-phase reaction.” In addition, it is predicted that “the effects of spin contamination and static correlation emerge when the dissociating adsorbate/surface interaction becomes larger than the covalent interaction in the adsorbate.” These will serve as a useful guide for correcting the errors in spin-polarized DFT/plane-wave calculations for surface reactions.

The spin contamination in force and frequency has not been corrected. For an accurate estimation of activation energies and transition states, corrections such as the AP-opt scheme [29,30] are indispensable. Developing a new method (AP-opt-DFT/plane-wave) that enables corrections is one of the important issues to be resolved in the future.