Article pubs.acs.org/JPCA
Modeling the Photoelectron Spectra of MoNbO2− Accounting for Spin Contamination in Density Functional Theory Lee M. Thompson and Hrant P. Hratchian* Chemistry and Chemical Biology, University of California, Merced, California 95343, United States
Downloaded by UNIV OF MANITOBA on September 7, 2015 | http://pubs.acs.org Publication Date (Web): August 3, 2015 | doi: 10.1021/acs.jpca.5b04625
S Supporting Information *
ABSTRACT: Spin contamination in density functional studies has been identified as a cause of discrepancies between theoretical and experimental spectra of metal oxide clusters such as MoNbO2. We perform calculations to simulate the photoelectron spectra of the MoNbO2 anion using broken-symmetry density functional theory incorporating recently developed approximate projection methods. These calculations are able to account for the presence of contaminating spin states at single-reference computational cost. Results using these new tools demonstrate the significant effect of spin-contamination on geometries and force constants and show that the related errors in simulated spectra may be largely overcome by using an approximate projection model.
■
INTRODUCTION Transition metal oxide surfaces are widely used as heterogeneous catalysts for atom-transfer reactions in applications including oxidation, hydrogenation, and desulfurization.1 A number of studies have shown that catalytic activity in these materials is localized along edge extended defects and at surface point defects, in particular, those resulting from oxygen vacancies.2 Development of new materials and improvement of existing catalysts require a fundamental understanding of possible defect types and their roles in enhancing or diminishing reactivity. Direct examination of these factors using only experimental methods is challenging for two reasons. First, the relatively low abundance of defects compared to regions of perfect lattice structure results in dilution effects. Second, under typical reaction conditions the surface is undergoing a number of structural changes, including defect formation, surface restructuring, subsurface oxygen diffusion, reagent adsorption and desorption, and reduction and oxidation reactions.3−5 As an alternative, many experimentalists have turned to small metal oxide cluster molecules as prototype systems of defects in bulk surfaces. Justification for this model comes from the local nature of defect sites, which exhibit similar incomplete valencies and nontraditional oxidation states to oxygen-deficient metal oxide clusters.6 To explore the types of electronic and molecular structures of metal oxide clusters, tools such as photoelectron spectroscopy are used for experimental characterization.7−16 The resulting spectra provide information on the electronic and vibrational states of the molecule, though interpretation can be complicated by several factors. Some of these difficulties arise from the specific class of molecules being studied. For example, there are a number of possible oxidation states all lying close in energy, and for each oxidation state there are a number of possible atomic arrangements that can all potentially exist as the global minimum energy structure. Other issues with interpretation are general to the spectroscopic method. In © 2015 American Chemical Society
particular, spectral congestion arising from multiple electronic transitions with numerous active vibrational modes generally scales with cluster size. This renders it increasingly difficult to distinguish bands arising from overtone and combination modes from those due to fundamental vibrations. Furthermore, as better resolution can be obtained using low electron kinetic energies, it is necessary to combine several spectra taken with different photon energies to span an energy range.17,18 This procedure can result in information loss in relative peak intensities. For these reasons, theoretical electronic structure methods have become a vital part of the toolkit for determining structure−function relationships in metal oxide catalysis.19−22 Electronic structure calculations allow prediction of the likely minimum energy structure(s), vibrational and electronic states of these structures, and prediction of the photoelectron spectra to compare with experiment. However, the electronic strain inherent in some of these systems can be problematic for single-reference methodologies, including almost all methods based on widely used density functional theory (DFT), and result in calculations that suffer from severe spin contamination errors. As this error is not uniform across different regions of the potential energy surface (PES), its presence can result in erroneous or questionable predictions of geometries, electronic state ordering, and vibrational modes.23−25 The use of multireference methods, which can resolve these shortcomings, can be ruled out as unfeasible for structures containing several transition metal centers due to computational cost and complexity, especially in studies that require second and higher order derivatives of the energy to obtain simulated spectra.26,27 A recent study of MoNbO2 found serious issues with the use of conventional hybrid density functionals in determining lowest energy structures and comparing computed and Received: May 14, 2015 Revised: June 25, 2015 Published: July 1, 2015 8744
DOI: 10.1021/acs.jpca.5b04625 J. Phys. Chem. A 2015, 119, 8744−8751
Article
The Journal of Physical Chemistry A Table 1. Lowest Energy Anionic Structures Computed Using B3PW91 without Correction for Spin Contamination structure
Downloaded by UNIV OF MANITOBA on September 7, 2015 | http://pubs.acs.org Publication Date (Web): August 3, 2015 | doi: 10.1021/acs.jpca.5b04625
002 002 011 011 011 011 011 020 110
electronic
nuclear
state
symmetry
⟨Ŝ2⟩
energy (eV)
R(Mo−Nb) (Å)
Cs C2v C1 C1 C1 C1 C1 C2v Cs
12.093 6.786 3.744 12.151 12.242 6.847 7.226 2.294 0.185
0.000 0.091 0.101 0.152 0.183 0.194 0.205 0.379 0.441
3.005 2.758 2.560 2.938 2.871 2.675 2.786 2.193 2.153
7
(1) (2) (1) (2)
A′ 5 B1 3 A 7 A 7 A 5 A 5 A 3 A2 1 A′
relative
theoretical photoelectron spectra.28 Close-lying structures with varying degrees of spin-contamination prevented a categorical assignment of the lowest energy structure and spin state, and cast doubt on assignments based on comparison of theoretical and experimental photoelectron spectra. In this study we apply our recently developed derivative methods for approximate projection (AP) models to theoretically determine structures, state ordering, and photoelectron spectra of MoNbO2−.29−35 We show that this method overcomes spin contamination effects and is able to reconcile problematic structures for which theory and experiment significantly disagree. Using the AP method, we reassign the structure and spin of the initial state and the associated transition responsible for the experimental photoelectron spectrum of MoNbO2−.
to systems for which even analytical projection schemes are too expensive. The AP method applies the BS solution to a Heisenberg Hamiltonian framework to resolve the system into substituent spin-pure states and gives rise to the energy expression LS LS HS EAP = αE BS + (1 − α)Eexact
where 2
α=
2
2
HS ⟨S ̂ ⟩exact − ⟨S ̂ ⟩LS BS
(3)
EHS exact
In eqs 2 and 3, is the energy of the BS solution, is the energy of the contaminating state which is taken to be spin2 ̂2 HS pure, ⟨Ŝ2⟩LS BS and ⟨S ⟩exact are the corresponding S expectation values, and Sz,LS is the z component of the total spin in the lowspin state. Recently developed analytical gradients and second derivatives have enabled efficient calculation of geometries, energies, and properties within an AP framework.34,35 Calculations were performed with a local development version of Gaussian37 using the UB3PW91 correlation functional with SDD relativistic pseudopotentials on the metal centers, an augmented SDD basis set for metal valence electrons, and aug-cc-pVTz for electrons on the oxygen atoms, referred to as SDDPlusTZ.28,38 The B3PW91 exchangecorrelation functional has been shown to be quite good for studies involving transition metal compounds, including systems with BS ground states.39,40 For each of the possible anionic spin states (singlet, triplet, quintet, and septet), and for each of the possible neutral spin states, (doublet, quartet, sextet, and octet), all possible structural isomers were considered. These structures were generated systematically by labeling each oxygen depending on whether it was at the Mo terminus (X), bridging the metal centers (Y), or at the Nb terminus (Z). The structure nomenclature XYZ indicates the number of oxygen atoms at each location relative to the metal centers. For example, 011 has one bridging oxygen and one terminal oxygen on Nb. Using this systematic approach, six structural isomers were identified. Conformers of each structural isomer were determined by considering the most symmetric point group and all possible subgroups, generating 232 initial geometries. Optimization was initially performed using the unprojected BS solution on each of the 232 structures. The stability of the initial and final wave functions were tested to ensure that they corresponded to the lowest energy electronic symmetry for each spin state. From these optimized structures, AP calculations were performed with the same functional and basis, removing the next highest multiplicity spin state
METHODS Broken symmetry (BS) approaches provide a convenient method for exploring open-shell and spin-polarized systems at single reference cost. An issue with BS approaches is that it is not possible to describe symmetry-adapted configurations using a single determinant, so that the resulting solutions to the wave function may not necessarily be an eigenfunction of the total spin-squared operator S2. As S2 and the non-relativistic Hamiltonian operator Ĥ commute, they necessarily require that any solution to the Hamiltonian must also be an eigenfunction of S2, and solutions that do not uphold this commutation relationship may suffer from spin contamination. The degree of spin contamination is given by the deviation of ⟨S2⟩ from that of the pure spin state, and in the case of DFT is complicated by the two-electron nature of the S2 operator. It is well-established that within the Kohn−Sham framework, treatment of the single determinant in a manner consistent with Hartree−Fock (HF) wave functions is an acceptable approximation.36 Therefore, ⟨S2⟩ can be given by
∑ |Sij ̅|2 ij
HS ⟨S ̂ ⟩exact − Sz ,LS(Sz ,LS + 1)
ELS BS
■
⟨S2⟩ = Sz(Sz + 1) + nβ −
(2)
(1)
where nβ and Si j are the number of β electrons and the matrix of α−β spatial overlaps of occupied molecular orbitals (MOs). When the self-consistent field (SCF) density breaks symmetry, ∑i j|Si j|2 < nβ, and it can be shown that the SCF energy of the spin-contaminated state is higher than that of the spin-pure state. In such cases the BS state can be shown to result from a mixing of higher-spin states with the pure-spin state of interest, resulting in spin contamination. While there are several possible methods that can be applied to resolve the issue of spin contamination, the Yamaguchi AP method can do so at mean-field cost and so is ideal for applying 8745
DOI: 10.1021/acs.jpca.5b04625 J. Phys. Chem. A 2015, 119, 8744−8751
Article
The Journal of Physical Chemistry A
is the presence of spin contamination, which is frequently measured by the deviation of ⟨S2̂ ⟩ from the expectation value of a pure spin state eq 1. Table 1 shows that the degree of spin contamination is not systematic across all members of this set, indicating that this source of error cannot likely be corrected with simple postprocessing approaches, such as numerical scaling of relative energies or frequencies. Instead, an explicit formalism addressing such effects is required. Further analysis of the ⟨S2̂ ⟩ values in Table 1 indicates that the high-spin states are the least contaminated, as would be expected given that the highest spin state can always be described by a single determinant to be spin pure. Furthermore, there is a correlation between the structure and the extent of spin contamination. The 002 structures (7A′ and 5B1) are effectively spin pure (
Modeling the Photoelectron Spectra of MoNbO2− Accounting for Spin Contamination in Density Functional Theory Lee M. Thompson and Hrant P. Hratchian* Chemistry and Chemical Biology, University of California, Merced, California 95343, United States
Downloaded by UNIV OF MANITOBA on September 7, 2015 | http://pubs.acs.org Publication Date (Web): August 3, 2015 | doi: 10.1021/acs.jpca.5b04625
S Supporting Information *
ABSTRACT: Spin contamination in density functional studies has been identified as a cause of discrepancies between theoretical and experimental spectra of metal oxide clusters such as MoNbO2. We perform calculations to simulate the photoelectron spectra of the MoNbO2 anion using broken-symmetry density functional theory incorporating recently developed approximate projection methods. These calculations are able to account for the presence of contaminating spin states at single-reference computational cost. Results using these new tools demonstrate the significant effect of spin-contamination on geometries and force constants and show that the related errors in simulated spectra may be largely overcome by using an approximate projection model.
■
INTRODUCTION Transition metal oxide surfaces are widely used as heterogeneous catalysts for atom-transfer reactions in applications including oxidation, hydrogenation, and desulfurization.1 A number of studies have shown that catalytic activity in these materials is localized along edge extended defects and at surface point defects, in particular, those resulting from oxygen vacancies.2 Development of new materials and improvement of existing catalysts require a fundamental understanding of possible defect types and their roles in enhancing or diminishing reactivity. Direct examination of these factors using only experimental methods is challenging for two reasons. First, the relatively low abundance of defects compared to regions of perfect lattice structure results in dilution effects. Second, under typical reaction conditions the surface is undergoing a number of structural changes, including defect formation, surface restructuring, subsurface oxygen diffusion, reagent adsorption and desorption, and reduction and oxidation reactions.3−5 As an alternative, many experimentalists have turned to small metal oxide cluster molecules as prototype systems of defects in bulk surfaces. Justification for this model comes from the local nature of defect sites, which exhibit similar incomplete valencies and nontraditional oxidation states to oxygen-deficient metal oxide clusters.6 To explore the types of electronic and molecular structures of metal oxide clusters, tools such as photoelectron spectroscopy are used for experimental characterization.7−16 The resulting spectra provide information on the electronic and vibrational states of the molecule, though interpretation can be complicated by several factors. Some of these difficulties arise from the specific class of molecules being studied. For example, there are a number of possible oxidation states all lying close in energy, and for each oxidation state there are a number of possible atomic arrangements that can all potentially exist as the global minimum energy structure. Other issues with interpretation are general to the spectroscopic method. In © 2015 American Chemical Society
particular, spectral congestion arising from multiple electronic transitions with numerous active vibrational modes generally scales with cluster size. This renders it increasingly difficult to distinguish bands arising from overtone and combination modes from those due to fundamental vibrations. Furthermore, as better resolution can be obtained using low electron kinetic energies, it is necessary to combine several spectra taken with different photon energies to span an energy range.17,18 This procedure can result in information loss in relative peak intensities. For these reasons, theoretical electronic structure methods have become a vital part of the toolkit for determining structure−function relationships in metal oxide catalysis.19−22 Electronic structure calculations allow prediction of the likely minimum energy structure(s), vibrational and electronic states of these structures, and prediction of the photoelectron spectra to compare with experiment. However, the electronic strain inherent in some of these systems can be problematic for single-reference methodologies, including almost all methods based on widely used density functional theory (DFT), and result in calculations that suffer from severe spin contamination errors. As this error is not uniform across different regions of the potential energy surface (PES), its presence can result in erroneous or questionable predictions of geometries, electronic state ordering, and vibrational modes.23−25 The use of multireference methods, which can resolve these shortcomings, can be ruled out as unfeasible for structures containing several transition metal centers due to computational cost and complexity, especially in studies that require second and higher order derivatives of the energy to obtain simulated spectra.26,27 A recent study of MoNbO2 found serious issues with the use of conventional hybrid density functionals in determining lowest energy structures and comparing computed and Received: May 14, 2015 Revised: June 25, 2015 Published: July 1, 2015 8744
DOI: 10.1021/acs.jpca.5b04625 J. Phys. Chem. A 2015, 119, 8744−8751
Article
The Journal of Physical Chemistry A Table 1. Lowest Energy Anionic Structures Computed Using B3PW91 without Correction for Spin Contamination structure
Downloaded by UNIV OF MANITOBA on September 7, 2015 | http://pubs.acs.org Publication Date (Web): August 3, 2015 | doi: 10.1021/acs.jpca.5b04625
002 002 011 011 011 011 011 020 110
electronic
nuclear
state
symmetry
⟨Ŝ2⟩
energy (eV)
R(Mo−Nb) (Å)
Cs C2v C1 C1 C1 C1 C1 C2v Cs
12.093 6.786 3.744 12.151 12.242 6.847 7.226 2.294 0.185
0.000 0.091 0.101 0.152 0.183 0.194 0.205 0.379 0.441
3.005 2.758 2.560 2.938 2.871 2.675 2.786 2.193 2.153
7
(1) (2) (1) (2)
A′ 5 B1 3 A 7 A 7 A 5 A 5 A 3 A2 1 A′
relative
theoretical photoelectron spectra.28 Close-lying structures with varying degrees of spin-contamination prevented a categorical assignment of the lowest energy structure and spin state, and cast doubt on assignments based on comparison of theoretical and experimental photoelectron spectra. In this study we apply our recently developed derivative methods for approximate projection (AP) models to theoretically determine structures, state ordering, and photoelectron spectra of MoNbO2−.29−35 We show that this method overcomes spin contamination effects and is able to reconcile problematic structures for which theory and experiment significantly disagree. Using the AP method, we reassign the structure and spin of the initial state and the associated transition responsible for the experimental photoelectron spectrum of MoNbO2−.
to systems for which even analytical projection schemes are too expensive. The AP method applies the BS solution to a Heisenberg Hamiltonian framework to resolve the system into substituent spin-pure states and gives rise to the energy expression LS LS HS EAP = αE BS + (1 − α)Eexact
where 2
α=
2
2
HS ⟨S ̂ ⟩exact − ⟨S ̂ ⟩LS BS
(3)
EHS exact
In eqs 2 and 3, is the energy of the BS solution, is the energy of the contaminating state which is taken to be spin2 ̂2 HS pure, ⟨Ŝ2⟩LS BS and ⟨S ⟩exact are the corresponding S expectation values, and Sz,LS is the z component of the total spin in the lowspin state. Recently developed analytical gradients and second derivatives have enabled efficient calculation of geometries, energies, and properties within an AP framework.34,35 Calculations were performed with a local development version of Gaussian37 using the UB3PW91 correlation functional with SDD relativistic pseudopotentials on the metal centers, an augmented SDD basis set for metal valence electrons, and aug-cc-pVTz for electrons on the oxygen atoms, referred to as SDDPlusTZ.28,38 The B3PW91 exchangecorrelation functional has been shown to be quite good for studies involving transition metal compounds, including systems with BS ground states.39,40 For each of the possible anionic spin states (singlet, triplet, quintet, and septet), and for each of the possible neutral spin states, (doublet, quartet, sextet, and octet), all possible structural isomers were considered. These structures were generated systematically by labeling each oxygen depending on whether it was at the Mo terminus (X), bridging the metal centers (Y), or at the Nb terminus (Z). The structure nomenclature XYZ indicates the number of oxygen atoms at each location relative to the metal centers. For example, 011 has one bridging oxygen and one terminal oxygen on Nb. Using this systematic approach, six structural isomers were identified. Conformers of each structural isomer were determined by considering the most symmetric point group and all possible subgroups, generating 232 initial geometries. Optimization was initially performed using the unprojected BS solution on each of the 232 structures. The stability of the initial and final wave functions were tested to ensure that they corresponded to the lowest energy electronic symmetry for each spin state. From these optimized structures, AP calculations were performed with the same functional and basis, removing the next highest multiplicity spin state
METHODS Broken symmetry (BS) approaches provide a convenient method for exploring open-shell and spin-polarized systems at single reference cost. An issue with BS approaches is that it is not possible to describe symmetry-adapted configurations using a single determinant, so that the resulting solutions to the wave function may not necessarily be an eigenfunction of the total spin-squared operator S2. As S2 and the non-relativistic Hamiltonian operator Ĥ commute, they necessarily require that any solution to the Hamiltonian must also be an eigenfunction of S2, and solutions that do not uphold this commutation relationship may suffer from spin contamination. The degree of spin contamination is given by the deviation of ⟨S2⟩ from that of the pure spin state, and in the case of DFT is complicated by the two-electron nature of the S2 operator. It is well-established that within the Kohn−Sham framework, treatment of the single determinant in a manner consistent with Hartree−Fock (HF) wave functions is an acceptable approximation.36 Therefore, ⟨S2⟩ can be given by
∑ |Sij ̅|2 ij
HS ⟨S ̂ ⟩exact − Sz ,LS(Sz ,LS + 1)
ELS BS
■
⟨S2⟩ = Sz(Sz + 1) + nβ −
(2)
(1)
where nβ and Si j are the number of β electrons and the matrix of α−β spatial overlaps of occupied molecular orbitals (MOs). When the self-consistent field (SCF) density breaks symmetry, ∑i j|Si j|2 < nβ, and it can be shown that the SCF energy of the spin-contaminated state is higher than that of the spin-pure state. In such cases the BS state can be shown to result from a mixing of higher-spin states with the pure-spin state of interest, resulting in spin contamination. While there are several possible methods that can be applied to resolve the issue of spin contamination, the Yamaguchi AP method can do so at mean-field cost and so is ideal for applying 8745
DOI: 10.1021/acs.jpca.5b04625 J. Phys. Chem. A 2015, 119, 8744−8751
Article
The Journal of Physical Chemistry A
is the presence of spin contamination, which is frequently measured by the deviation of ⟨S2̂ ⟩ from the expectation value of a pure spin state eq 1. Table 1 shows that the degree of spin contamination is not systematic across all members of this set, indicating that this source of error cannot likely be corrected with simple postprocessing approaches, such as numerical scaling of relative energies or frequencies. Instead, an explicit formalism addressing such effects is required. Further analysis of the ⟨S2̂ ⟩ values in Table 1 indicates that the high-spin states are the least contaminated, as would be expected given that the highest spin state can always be described by a single determinant to be spin pure. Furthermore, there is a correlation between the structure and the extent of spin contamination. The 002 structures (7A′ and 5B1) are effectively spin pure (
