The effect of the Perdew-Zunger self-interaction correction to density functionals on the energetics of small molecules

Self-consistent calculations using the Perdew-Zunger self-interaction correction (PZ-SIC) to local density and gradient dependent energy functionals are presented for the binding energy and equilibrium geometry of small molecules as well as energy barriers of reactions. The effect of the correction is to reduce binding energy and bond lengths and increase activation energy barriers when bond breaking is involved. The accuracy of the corrected functionals varies strongly, the correction to the binding energy being too weak for the local density approximation but too strong for the gradient dependent functionals considered. For the Perdew, Burke, and Ernzerhof (PBE) functional, a scaling of the PZ-SIC by one half gives improved results on average for both binding energy and bond lengths. The PZ-SIC does not necessarily give more accurate total energy, but it can result in a better cancellation of errors. An essential aspect of these calculations is the use of complex orbitals. A restriction to real orbitals leads to less accurate results as was recently shown for atoms [S. Kl\"upfel, P. Kl\"upfel, and H. J\'onsson, Phys. Rev. A {84}, 050501 (2011)]. The molecular geometry of radicals can be strongly affected by PZ-SIC. An incorrect, non-linear structure of the \ce{C2H} radical predicted by PBE is corrected by PZ-SIC. The \ce{CH3} radical is correctly predicted to be planar when complex orbitals are used, while it is non-planar when the PZ-SIC calculation is restricted to real orbitals.


I. INTRODUCTION
Density functional theory (DFT) 1,2 has become a powerful tool for physicists and chemists to describe the electronic structure of atoms, molecules, and solids. While being exact in theory, practical applications of DFT rely on approximations of the exchange-correlation (xc) functional. The local spin density approximation (LSD) 2 is in most cases too crude an approximation for the study of molecular systems. Functionals based on the generalized gradient approximation (GGA) 3-8 allow for a more accurate description of the inhomogeneous electron densities of molecules but also turn out not to be accurate enough in many cases. The accuracy of GGA functionals for applications of chemical interest can be further improved by admixture of a fraction of exact exchange in the form of hybrid functionals 9 . The B3LYP 10-14 hybrid functional has become widely used in molecular studies. More recently a number of new functionals have been developed to reproduce certain chemical or molecular properties to high accuracy, while being less accurate for others 15 .
In the spin-unrestricted Kohn-Sham (KS) formalism 2,16 the energy of a system of N electrons in an external potential v ext (r) is defined through the spin-densities ρ ↑ + ρ ↓ = ρ by where V ext is the electrostatic interaction energy of the density with the external field and E H is the Hartree energy, the classical Coulomb repulsion of the charge density with itself. T s is the kinetic energy of the noninteracting reference system, constructed from the set of KS-orbitals {ϕ i } to produce the same density as the exact wave function. The exchange-correlation energy, E xc , contains all remaining contributions to the exact energy. The xc-energy can be split into the sum of exchange, E x , and correlation energy, E c . For any one-electron density, ρ 1 , the two conditions and are fulfilled by the exact functional. The second condition can be satisfied by a semi-local form of the correlation functional as the one proposed by Lee, Yang, and Parr (LYP) 13 . It is, however, not possible to formulate an exchange functional that can, for any possible ρ 1 , compensate the non-local Hartree energy merely from local information of the density. For approximate semilocal functionals condition (2), or both conditions, (2) and (3), are violated. This gives rise to a self-interaction error (SIE), Perdew and Zunger proposed a self-interaction correction scheme (SIC) in which the SIE of the individual orbitals, as defined by Eq. (4), is subtracted from the total energy 17 . The Perdew-Zunger self-interaction correction (PZ-SIC) energy functional, depends not only on the total spin-densities, but also on the orbital densities, ρ N = (ρ 1 , . . . , ρ N ) with ρ i = |ϕ i | 2 , and can in principle be applied to any approximate xcfunctional.
For the exact functional, the correction term vanishes for any many-electron system, as can be seen from Eqs. (2) and (3). For approximate functionals, the PZ-SIC is accurate for any one-electron density but for manyelectron systems it is in general only an approximate correction, as the magnitude of the many-electron selfinteraction error does not have to be the sum of the individual SIE terms of Eq. (4) 18 . While this orbital based estimate cannot be expected to eliminate all selfinteraction for many-electron systems, it may improve the accuracy of approximate functionals. The purpose of the present study is to test the accuracy of PZ-SIC for a few commonly used functionals when applied to molecules.
Errors in the energy due to approximation of the xcfunctional stem from both an incomplete cancellation of the electron self-interaction and an inaccurate description of the inter-electronic interaction. The accuracy of the PZ-SIC does, therefore, depend on the functional approximation it is used with. From earlier studies, it was concluded that PZ-SIC often overcorrects errors in calculated observables of many-electron systems such as equilibrium bond lengths and atomization energy. A scaled down modification of the PZ-SIC functional has been proposed in the form of (6) and for the Perdew, Burke, and Ernzerhof (PBE) functional a factor α in the range of 0.4 to 0.5 could improve the description of some observables [19][20][21] . More elaborate scaling schemes have also been proposed 22,23 .
A recent study of atoms with PZ-SIC showed that the total energy can be lowered significantly by allowing the orbitals to be complex functions 24 . This improved significantly the results when PZ-SIC was applied to the PBE functional, while a restriction to real orbitals led to much larger errors than those of the uncorrected functional 25 . According to the variational principle, the addition of an imaginary component to the orbitals can only lower the total energy. The effect on equilibrium geometry or energy differences such as atomization energy or energy barriers is not monotonous in a similar way. To our knowledge, previously published fully variational studies of the energetics of molecules using stationary PZ-SIC have exclusively been based on real orbitals. Only recently complex orbitals were used in a study focusing on bond-lengths of molecules within an approximate Kohn-Sham interpretation of the PZ-SIC 26 .
We present results on the ground state geometry and atomization energy of a set of 17 molecules, the equilibrium structure of two 'problematic' radicals and the energy barrier of four reactions. We studied the effect of SIC using three different functionals: the local spin density approximation (using Slater exchange [27][28][29] and the Perdew-Wang parameterization of correlation 30 , SPW92), the two generalized gradient approximations of Perdew, Burke, and Ernzerhof 31 , and Becke's exchange functional 12 and correlation of Lee, Yang, and Parr 13 (BLYP). The results of self-consistent calculations using PZ-SIC (SIC) as well as the scaled down modification with a factor of one half (SIC/2) are presented. For comparison, calculations using the two hybrid functionals, PBE0 32 and B3LYP 10 , were also carried out, as well as less accurate SIC calculations using real orbitals.

II. COMPUTATIONAL METHOD
The energy minimum with respect to variation of the orbitals under the constraint of orthonormality is described by the two sets of equations 33,34 , where the orbital-specific Hamiltonians are defined bŷ The Lagrange multipliers can be determined by projection of Eq. (7) as λ ji = ϕ j |Ĥ i |ϕ i . In contrast to semilocal functionals, the matrix of Lagrange multipliers is not Hermitian for any unitary transformation among the occupied orbitals. The minimum energy is determined both by the space spanned by the set of orbitals and by the linear combination of the orbitals.
The effective potential will not be the same for all orbitals and this places SIC outside the domain of Kohn-Sham DFT. It can be treated as a true Kohn-Sham functional by means of the optimized effective potential method (OEP) 35,36 , but in many applications the functional is treated in the generalized Kohn-Sham framework, i.e., the energy is minimized with respect to variation of the orbitals, resulting in different potentials for each one of them.
Analytical gradients of the energy of SIC functionals have been derived 37 and can be used in a direct minimization of the energy. The efficiency can be improved by an additional step in the iterative minimization. Before the orbitals are altered according to the energy gradient, the unitary transformation that minimizes the SIC energy terms is found. By such a 'unitary optimization' the convergence rate can be greatly improved 38,39 . An efficient algorithm for the unitary optimization has recently been developed 40 .
The calculations were carried out with the Gaussiantype orbital based program Quantice 41 . As analytical gradients of the atomic positions are not available in the program, the molecular structure has been optimized manually. A sequence of internuclear distances and angles was sampled and the minimum energy configuration found by cubic interpolation. The equilibrium structure was confirmed by comparison of the interpolated energy with the calculated energy at the interpolated geometry. This manual scheme limits the size of the molecules and the number of structural degrees of freedom that can be optimized in the current version of this software. In the geometry optimization, all analogous bonds in a molecule were constrained to have the same length to reduce the computational effort compared to a completely unconstrained structural relaxation.
For all calculations, atom-centered grids of 75 radial shells 42 of a 302-point Lebedev-Laikov grid 43 were combined to form a multicenter integration grid 44 . The Cartesian representation of correlation-consistent polarized valence-triple-zeta 45,46 (cc-pVTZ) basis sets were used for most calculations. For the reaction barriers involving only hydrogen, a quadruple-zeta basis set 46 (cc-pVQZ) was used. The convergence criterion of the electronic minimization was 10 −8 Ha 2 for the norm of the energy gradient. All calculations were using unrestricted orbitals, starting from a random initialization.

III. PZ-SIC ORBITALS
The orbital density dependence of the PZ-SIC energy expression results in a set of well defined 'optimal orbitals,' defined by a unitary transformation of the canonical orbitals. For atomic systems, these orbitals have often been found to resemble sp n hybrid orbitals 47 . When the orbitals are allowed to be complex functions, the optimal orbitals still resemble hybrid orbitals but can have significantly different shape and orientation than hybrids of real orbitals. For neon, e.g., real orbitals produce a set of sp 3 orbitals with tetrahedral orientation, but complex orbitals produce a set with tetragonal orientation without nodal surfaces 24 . For molecules, the optimal orbitals often take forms that are consistent with 'chemical intuition,' i.e. they can be interpreted as lone pairs, single, or multiple bonds. Also, orbitals consistent with the more 'exotic' three-center or banana bonds can form. Also here, the shape of the complex orbitals is often rather different from that of the real ones, but the interpretation above can still be retained in many cases. Figure 1 shows the optimized valence orbitals of one spin channel for N 2 , obtained from PBE+SIC using complex (a) and real (b) orbitals. The real orbitals are of two kinds: three spatially degenerate orbitals that add up to a triple bond, and two lone pairs. The triple bond is built up from two different kinds of orbitals. One orbital has the charac- ter of a sigma-bond with rotational symmetry about the molecular axis. The other two orbitals are degenerate and symmetric about a plane going through the molecular axis. The complex lone pairs do not differ much qualitatively from the real orbitals. The real orbitals are rather localized and in staggered orientation, but the complex orbitals are more delocalized and to a larger extent share the same space.
To test the accuracy of the various density functional approximations, the predicted atomization energy, i.e., the difference in the total energy of the atoms constituting a molecule and the total energy of the molecule, was calculated and compared with experimental estimates corrected for zero-point energy 48 or, in the case of H 2 , with an accurate theoretical result 49 . The mean error (ME) and mean absolute error (MAE) for each functional approximation is shown in Figure 2. For comparison, results of SIC and SIC/2 calculations restricted to real orbitals are also shown. The numerical values are listed in Table I, except for the results obtained with real orbitals. The LSD energy shows the largest deviation of all functionals, strongly over binding the molecules. The errors are reduced by applying SIC (see LSD+SIC), but this only partially eliminates the errors. The GGA functionals, PBE and BLYP, reduce the errors of LSD significantly, but still predict most molecules to be too stable. The binding energy is reduced by SIC for both functionals, but the correction is too large. For PBE+SIC, the MAE is actually slightly increased and is doubled for BLYP+SIC. The mean deviation is greatly reduced by applying SIC scaled by one-half, PBE+SIC/2, while the MAE is just slightly smaller than for the uncorrected functional. BLYP+SIC/2 gives smaller errors than BLYP+SIC but still predicts the binding energy to be too low. Calculations using SIC that are restricted to real orbitals predict lower atomization energy on average.
Vydrov et al. studied the effect of PZ-SIC on the heat of formation using several functionals using calculations restricted to real orbitals 25 . They concluded that PZ-SIC only improves the results for LSD, while larger deviations are found when the correction is used with GGA functionals. As shown here, it is important to allow the orbitals to be complex functions. This reduces the over correction, but does not eliminate it. Better agreement with the reference data is obtained by scaling the SIC. No fitting of the scaling factor was carried out, but a factor of one-half chosen to illustrate the trend. The scaled SIC Optimized orbitals corresponding to the single bond in F 2 calculated with (a) PBE+SIC/2 and (b) PBE+SIC. Both spin-up and spin-down orbitals are shown. For PBE+SIC/2, the total density is not spin polarized. For PBE+SIC, the orbitals are localized to some extent on one of the atoms and the electron density is spin-polarized.
used with the PBE functional gives a smaller mean error than the PBE0 hybrid functional, however, the MAE indicates that some molecules are over bound while others are too unstable. A systematic under binding is found for molecules containing hydrogen, except for CH 2 and C 2 H 2 .
An extreme case of the under binding obtained from the BLYP+SIC functional is F 2 , which is predicted to be unstable. Also, PBE+SIC gives severe underestimation of the binding energy and elongation of the bond as shown in Table II. The bond energy is overestimated by all uncorrected functionals. For LSD, SIC/2 and SIC reduce the bond energy, the latter giving a value closer to experiment. For PBE, the correction greatly reduces the binding energy, resulting in an underestimation already for SIC/2 and predicting a very weakly bound molecule for SIC. This effect is even more pronounced for BLYP, where SIC predicts the molecule not to be bound at all. Usually, SIC/2 gives results that are intermediate of those obtained by SIC and by the uncorrected functional. This is not the case, however, for F 2 in the GGA functionals. The bond length is reduced by SIC/2, but increased when full SIC is applied.
The optimal orbitals obtained using PBE+SIC/2 and PBE+SIC are qualitatively different, as shown in Figure 3. The two orbitals (one each for spin-up and spindown) corresponding to the single bond are similar for PBE+SIC/2. The slight difference in shape does not result in spin polarization of the total density, as it is compensated by the density of the lone pairs. For PBE+SIC, the orbitals are distorted and localize one on each of the two nuclei. The total density is spin polarized and the electronic structure can be interpreted as an intermediate state towards two separated fluorine atoms. The effect of SIC on the molecule and the single atoms is not balanced. The total energy of the molecule is predicted to be too high relative to that of the atoms, resulting in a weak bond. For BLYP+SIC the 'correction' is unbalanced to such an extent that the energy of F 2 is above that of the atoms for all nuclear separations. Building the bond in a diatomic molecule can be described as a delocalization of atomic hybrid orbitals over both nuclei, accompanied by changes in the shape of the orbitals to maintain orthogonality as well as relaxation of the orbitals not participating in the bond. The SIC energy of an orbital varies strongly depending on how localized it is. The valence atomic orbitals become more compact with increasing atomic number within a row of the periodic table. At the same time, the magnitude by which SIC reduces the atomization energy, increases for all functionals going from N 2 to O 2 to F 2 . Moving from the second to the third row of the periodic table, the valence orbitals become more delocalized, and the changes in atomization energy due to SIC are smaller as can be seen by comparing P 2 with N 2 . Preliminary results of the binding energy of larger molecules at unrelaxed geometry reveal similar trends. These trends indicate that the effect of SIC on atomization energy is more pronounced if localized atomic orbitals participate in the bonding. However, without taking into account the changes in the orbital shape and the rearrangement of non-bonding orbitals, such a simplified interpretation is insufficient to explain all the observed effects and a more detailed study is required.
Calculations of observables from energy differences usually are more accurate than the individual total energy values because of partial cancellation of errors. One source of errors is the limited basis set used in the cal-culations. This error can in theory be eliminated by systematically increasing the size of the basis set until a complete basis is reached, or in a more practical way, until the calculated energy differences do not change significantly. Even if a complete basis set is used, an error remains from the approximate energy functional. A functional with large errors in the predicted total energy must be seen as too crude an approximation of the exact functional. This, however, does not mean that it can not be a useful functional for practical calculations. If the  functional has a rather constant error per electron for all systems, the energy difference of systems of the same number of electrons still can be predicted accurately. Figure 4 correlates the errors in total energy of the molecules to that of the constituent atoms. Each circle depicts one molecule, the x-axis shows the error per electron of the molecule, the y-axis that of the atoms. The diagonal dashed line corresponds to a perfect cancellation of errors, for points above or below the line, atomization energy is overestimated or underestimated, respectively. The diameter of the circles indicates the number of electrons in each system. The molecules containing hydrogen are indicated by circles filled in grey.
In this representation, the effect of SIC and the admixture of exact exchange can be studied in more detail than from the atomization energy alone. The total energy of both the atoms and the molecules is overestimated by PBE, around 0.2 eV per electron for most molecules. The largest deviation is found for P 2 , the largest system of the test set. With the exception of H 2 (the smallest point), the errors of the molecules and atoms are similar, but some spread around the diagonal is observed. Half SIC and full SIC reduce the absolute errors for most systems with most points moving closer to the origin. However, the magnitude of the correction is in many cases different for atoms than molecules. Going from PBE to PBE+SIC/2 to PBE+SIC, the points are shifted down and the vertical spread is reduced, which corresponds to smaller and more similar errors per electron for the atoms. At the same time, however, the horizontal spread increases, indicating a less systematic correction of the errors in the total energy of the molecules. This can be observed in particular for the molecules containing hydrogen, indicated as grey points. When PBE is used, the errors for the molecules are more systematic than the errors for the atoms. The opposite trend is observed when PBE+SIC is used. For PBE0, the absolute magnitude of the errors becomes less systematic for both atoms and molecules, the spread is increased along the diagonal. At the same time, the spread perpendicular to the diagonal is reduced, improving the overall cancellation of errors.
For BLYP the total errors are much smaller than for PBE, and the cancellation is slightly better. BLYP+SIC/2 reduces the total errors for many systems, but underestimates the atomization energy. For SIC, the points move further below the diagonal but also spread both vertically and horizontally, indicating an unsystematic effect of SIC on the total energy of both atoms and molecules. B3LYP actually increases the magnitude of the errors over BLYP, predicting too low energy for all the atoms and molecules. The errors are, however, well balanced and cancellation of errors results in the superior performance of B3LYP with respect to atomization energy, as shown in Figure 2. This improved cancellation of errors can to some extent be explained by the origin of this hybrid functional. B3LYP is based on the B3PW91 functional, for which the parameters had been fit to reference data composed of the total energy of ten atoms and 106 energy differences 11 . Such a fitting procedure places less weight on the accuracy of total energy than on cancellation of errors. The same parameters are used in the B3LYP functional that, despite its different functional components, predicts energy differences often very accurately while atomic total energy is predicted to be too low.

V. EQUILIBRIUM GEOMETRY
The equilibrium geometry of a molecule is found as the minimum of the potential energy surface, determined from the total energy of the molecule at different geometries. Again, a constant error at all geometries will still allow for the prediction of the correct equilibrium structure, while varying errors can result in quantitatively or even qualitatively wrong structures. Figure 5 shows the mean and mean absolute error for the equilibrium bond lengths of the set of molecules and functionals studied, excluding F 2 , as it was shown to be not stable in BLYP+SIC. The numerical values are listed in Table III. The uncorrected functionals overestimate the bond lengths on average slightly by ∼1pm, the GGA functionals generally overestimate the bond length, whereas some molecules are predicted to have too short bonds by LSD. PZ-SIC results in a strong overcorrection. For LSD+SIC all bonds are predicted to be much too short. For PBE+SIC and BLYP+SIC, the overcorrection is smaller. Here, all bond lengths except that of Li 2 (and F 2 ) are predicted to be too short.
The scaled SIC overcorrects LSD but gives on average improved results for the GGA functionals. Still, as in the case of atomization energy, the mean absolute error shows large fluctuations. A restriction to real orbitals has less effect on the bond length than on the atomization energy but gives slightly worse results except for BLYP+SIC/2 and BLYP+SIC. The hybrid functionals give the highest accuracy. Figure 6 shows the deviation of the equilibrium bond angles from experimental geometries for the non-linear molecules. For H 2 O and NH 3 , the angle from LSD and the hybrid functionals are in very good agreement with experiment, while the GGA functionals predict angles 0.5 • -1.0 • too small. SIC/2 and SIC predict larger angles for all functionals, with a monotonous increase from SIC/2 to SIC. For the full correction, the angle in water is close to that of a regular tetrahedral structure, indicated by a dotted line, in ammonia it even exceeds this angle, in particular for BLYP+SIC. The localized nature of the optimal orbitals can motivate an interpretation along the lines of valence shell electron pair repulsion (VSEPR) theory. The increase in bond angles indicates a relatively stronger interaction between bonding electron pairs compared to the repulsion between a lone pair and a bond pair. In methylene, the bond angle is predicted to be too large by the LSD functional but is quite accurately predicted by the GGA functionals. The angle is reduced in LSD+SIC, to good agreement with experiment, while it increases for the GGA+SIC functionals. In all cases, the angles predicted by using real orbitals with SIC are lower than when complex orbitals are used, and except for methylene closer to experiment. In contrast to most equilibrium bond lengths, restriction to real orbitals has a strong effect on the equilibrium bond angles. The hybrid functionals here also give results that are closer to experimental results.

VI. STRUCTURE OF MOLECULAR RADICALS
Gräfenstein et al. found that the equilibrium structure of the CH 3 radical predicted by BLYP+SIC is non-planar 51 , in disagreement with both experimental observation 50 and ab initio calculations 51 . We determined the ground state geometry, restricted to C 3v symmetry, and found that all the uncorrected functionals predict the proper planar structure. Results obtained from applying SIC are listed in Table IV as the energy difference between the equilibrium structure and the planar structure of lowest energy, and as the out-of-plane angle, i.e., the angle enclosed by a C-H bond and the plane defined by the hydrogen atoms. Using real orbitals, the non-planar geometry is preferred in both SIC and SIC/2 for all the functionals. Both the angle and energy difference increase from SIC/2 to SIC and are larger for the GGA functionals than for LSD. Using complex orbitals, however, the self-interaction corrected GGA functionals predict the geometry to be qualitatively correct, while it merely reduces the angle and energy difference for LSD, retaining the incorrect pyramidal structure. The destabilization of the planar structure by using SIC with real orbitals can be understood from the hybridization of the optimized orbitals. Figure 7 depicts the complex and real optimized valence orbitals of the spin majority for the planar equilibrium structure predicted by PBE+SIC. The complex orbitals corresponding to C-H σ-bonds lie in the plane. The orbital of the unpaired electron is delocalized symmetrically over both sides of the plane with an increased density between two of the bonding orbitals. The shape of these two orbitals differs slightly from the shape of the third one. The total and spin density have, despite the reduced symmetry of the optimal orbitals, the full symmetry of the molecule.
The shape of the real orbitals is quite different. The orbital of the unpaired electron takes the form of a real sp 3 hybrid. As the orbitals have to be orthogonal, the binding orbitals are forced into an unnatural shape, 'bend-ing' out of plane between carbon and hydrogen. While an sp 2 configuration would seem more favorable, this would force the unpaired orbital to be an unhybridized p-orbital, which is higher in energy. The SIC energy contribution of the unpaired orbital is lowered to such an extent by the hybridization, that it compensates for the higher SIC energy terms of the 'banana-bonds.' The total energy can be lowered further, by moving the hydrogen atoms out of the plane, which results in a geometry in TABLE IV. Equilibrium structure of the CH 3 radical. The 'out of plane' angle, α, in degrees and energy difference between planar and pyramidal structure, ∆E, in meV is shown for the various SIC functionals for complex (c.) and real (r.) orbitals. The uncorrected functionals all predict the correct planar ground state. better agreement with the sp 3 hybridization on the carbon atom. For LSD+SIC (not shown), the unpaired optimal orbital appears as an intermediate between the real and complex orbital shown for PBE+SIC. The orbital is partially delocalized but is not symmetric with respect to the plane. For atoms, it was found that extending the variational space to complex orbitals lowers the total energy in GGA+SIC more than in LSD+SIC 24 . This can play a role in the equilibrium geometry. The SIC energy of the unpaired and binding orbitals will be affected differently by the additional complex degrees of freedom and a comparison between GGA+SIC and LSD+SIC might give insight into the origin of the insufficient correction found for LSD+SIC. However, as the bond lengths are different for real and complex orbitals, all terms of the energy functional change, making such an analysis more difficult.
Recently, Oyeyemi et al. found that PBE predicts an incorrect equilibrium structure of the ethynyl radical, C 2 H 52 . The structure is not linear, but rather bent by ∼166 • . This was attributed to an over delocalization of the electron density due to the self-interaction error of approximate functionals. By including exact exchange in the form of hybrid functionals, the correct linear ground state geometry is obtained. Figure 8 shows the energy of the bent radical relative to the energy of the linear structure. For all angles, bond lengths have been optimized. For PBE the energy drops slightly when the molecule starts to bend, with an optimal angle in good agreement with the previous study 52 . For PBE+SIC and PBE+SIC/2, the linear geometry is favored. The selfinteraction correction corrects the Hartree self-energy, as does exact exchange. In cases where inaccuracy of approximate functionals stems from the spurious selfrepulsion of the orbitals, (scaled) PZ-SIC corrects in a way that is analogous to hybrid functionals which include scaled exact exchange.

VII. REACTION BARRIERS
The energy barrier for four reactions has been calculated to study the effect of SIC on the activation energy. The energy difference between the reactant minimum and the lowest saddle point of the potential energy surface is calculated. The reactions are listed in Table V   teraction (CI) 53 . The barrier heights are calculated with respect to the reactants, using the accurate total energy of the hydrogen molecule 49 and the exact energy of the hydrogen atom as reference values. For HFH, only the collinear, symmetric saddle point was considered, as a reference for this barrier is available 54 , and the computational effort was too large to find the non-collinear saddle point which is slightly lower in energy 55 . The inversion barrier height of NH 3 had been calculated using the coupled cluster method, CCSD(T) 56 . Figure 9 shows the deviation of the calculated saddle points from the reference values. For H 4 and H 3 , the barrier height predicted by the uncorrected functionals is always too low, and is increased by SIC/2 and SIC. For both reactions, SIC improves the energy barrier but it is still underestimated for H 4 . The bond lengths at the saddle point decrease for SIC/2 and further for SIC, as is also found for most molecules. Compared to the bond lengths predicted by the uncorrected functionals, this increases the deviation for H 4 and overcorrects the slightly too large bonds for H 3 .
The HFH barrier is underestimated by the uncorrected functionals and increases with SIC, monotonously from SIC/2 to SIC, the latter giving an overcorrection. For the GGA functionals the bond lengths do not change monotonously but behave similar to the F 2 bond length discussed above. The bond length decreases for SIC/2 and agrees better with the reference value, but then increases for PBE+SIC, while hardly changing for BLYP+SIC. For the functionals studied here, the GGA+SIC/2 barriers give the best results, being more accurate than the hybrid functionals. For this reaction, in contrast to the hydrogen barriers, different results are obtained when using real orbitals. Here, the bond length is larger and the barrier height is slightly lower than for complex orbitals. These three barriers describe bond breaking situations and are underestimated by the uncorrected functionals. The ammonia inversion is qualitatively different, as bonds merely rearrange and their lengths at equilibrium and saddle point do not differ significantly. The calculated bond lengths in the planar configuration are elongated by less than 2 pm compared to the equilibrium structure 56 . The inversion barrier calculated with the uncorrected functionals is in good agreement with the reference, but the bonds are predicted to be too long. The bond length decreases and the barrier increases with SIC for LSD, but for the GGA functionals the barrier is lowered and becomes underestimated. Using real orbitals, on the other hand, the barrier increases strongly. This qualitatively different effect of SIC can be explained by the structure of the 'lone pair,' as was done for the methyl radical. For the planar molecule, the orbitals of both spins look very similar to those in Figure 7. For the complex orbitals, spin up and spin down orbitals are rotated with respect to each other by 120 • . In the real case, the spin-up and spin-down orbitals are symmetric with respect to reflection through the plane. As for CH 3 , this configuration of the real orbitals is higher in energy than that of the complex orbitals, which explains the large difference in the predicted energy barrier.
As for the atomization energy, the errors in the calculated reaction barriers can be analyzed in greater detail by comparison of the errors in total energy of the species involved. Figure 10 shows the error per electron of the total energy of the saddle point structure (S), the separated reactants (R) and that of the separated atoms (A), calculated using PBE, PBE+SIC/2, and PBE+SIC. Points on the dashed line correspond to a perfect description of the total energy, the vertical difference between neighboring points corresponds to deviations in the atomization energy (A-R) and energy barrier (R-S). For the H 3 systems, the PBE energy of the hydrogen atoms is quite accurate and the saddle point energy is only slightly underestimated. The underestimation of both the barrier height and atomization energy originates mainly from an inaccurate description of the hydrogen molecule. When SIC is applied, the energy of the atom is lowered insignificantly and the energy of the molecule is raised slightly, making both total and atomization energy of the molecule less accurate. However, PBE+SIC introduces a large error at the saddle point, which cancels the error in the molecule. For SIC/2, the errors introduced are less well balanced and the barrier is still underestimated.
The H 4 barrier is severely underestimated by the PBE functional stemming from equally large errors of opposite sign in the reactants and saddle point structure. The energy of the saddle point is predicted accurately by PBE+SIC. Each optimized orbital is to great extent localized at one of the hydrogen atoms and neighboring orbitals are of opposite spin. In such a configuration, most of the exchange-correlation energy can be expected to be self-interaction energy, in which case PZ-SIC usually performs well. In this case the barrier is improved by SIC, but as the energy of the hydrogen molecule is too high, a still significant error in the barrier remains.
For the atomic reference of HFH and NH 3 , both the error in the total energy and the change when applying SIC are dominated by the heavy atoms. The reduction of the error by SIC is smaller for the HF molecule than for the atoms, thus making the binding energy less accurate compared to PBE. The error in the saddle point energy is increased and exceeds that of the reactants. Here, half SIC gives a better cancellation of error than both PBE and PBE+SIC. For the NH 3 molecules, the effect is smaller and the barrier height changes only slightly.
For the H 3 reaction, we have calculated the energy along a path connecting the saddle point and the reactants. For a number of distances d H-H 2 , the bond length of the H 2 fragment was relaxed. In Figure 11, the relaxed energy of the LSD and PBE functionals is shown for several separations of the fragments. The leftmost points correspond to the bond length and energy barrier of the H 3 saddle point. Increasing values of d H-H 2 describe the separation into a hydrogen atom and a H 2 molecule. In LSD, the combined system has lower energy than the separated fragments for all separations, the energy barrier is negative and H 3 is predicted to be a stable molecule. LSD+SIC/2 and LSD+SIC increase the H 3 energy and predict a reaction barrier. The height is however still underestimated and at intermediate distances, a system more stable than the reactants is found. In PBE, no intermediate configuration is more stable than the reactants, but the energy difference is globally underestimated. PBE+SIC agrees well with the CI results while the PBE0 hybrid does not describe the system as well. However, as shown in Figure 10 the improved description by the PBE+SIC functional results only from a better cancellation of errors and not from a better description of the system.

VIII. SUMMARY AND CONCLUSION
The results presented here on the energetics of small molecules provide insight into both the strengths and the shortcomings of the Perdew-Zunger self-interaction correction. The qualitatively wrong equilibrium structure of C 2 H predicted by PBE is corrected by SIC. In the case of the CH 3 radical, SIC had been found to predict an incorrect geometry, but in the present study this was shown to be an artifact of a restriction to real orbitals. This molecule and the H 3 potential energy surface emphasize the importance of using SIC with GGA functionals to obtain more accurate results, as had already been shown in a study of the total energy of atoms 24 . However, not all GGA functionals are suited for the application of SIC. The F 2 molecule is found to be unstable when SIC is applied to the BLYP functional and the errors in to-tal energy are large and unsystematic for the rest of the molecules in the test set as well as the isolated atoms.
Analysis of the total energy revealed that for PBE+SIC, the errors per electron in the atomic systems are reduced and show less fluctuations than that of the molecules, resulting in an unsystematic effect on the atomization energy. The Perdew-Zunger SIC does not in general result in an overcorrection when applied to PBE, as the predicted total energy is still too high. However, an overcorrection in calculated atomization energy is usually observed which can be improved by scaling the SIC. The simple scaling scheme of using a constant factor of one-half for all orbitals when applied with PBE reduces the mean error in atomization energy to less than that of the PBE0 hybrid functional, but still gives significant absolute errors. For systems with only a few electrons, this half-SIC approximation does not perform better than full SIC in calculations of energy barriers. More flexible scaling 22,23 where full SIC is retained for one-electron systems or isolated orbitals could work better for such systems. These scaling schemes have so far only been used in combination with real orbitals and their performances would need to be reassessed using complex orbitals.