Systems with weakly bound extra electrons impose great challenges to semilocal density functional approximations (DFAs), which suffer from self-interaction errors. Small ammonia clusters are one such example of weakly bound anions where the extra electron is weakly bound. We applied two self-interaction correction (SIC) schemes, viz., the well-known Perdew–Zunger and the recently developed locally scaled SIC (LSIC) with the local spin density approximation (LSDA), Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA), and the SCAN meta-GGA functionals to calculate the vertical detachment energies (VDEs) of small ammonia cluster anions (NH3)n−. Our results show that the LSIC significantly reduces the errors in calculations of VDE with LSDA and PBE-GGA functionals leading to better agreement with the reference values calculated with coupled cluster singles and doubles with perturbative triples [CCSD(T)]. Accurate prediction of VDE as an absolute of the highest occupied molecular orbital (HOMO) is challenging for DFAs. Our results show that VDEs estimated from the negative of HOMO eigenvalues with the LSIC-LSDA and Perdew–Zunger SIC-PBE are within 11 meV of the reference CCSD(T) results. The LSIC method performs consistently well for the VDE estimates, from both the total energy differences and the absolute HOMO eigenvalues.

## I. INTRODUCTION

Systems with solvated electrons have attracted much attention due to their role in many chemical and biological processes.^{1} In such systems, an extra electron is bound as in a cavity such that it polarizes the surrounding dielectric medium. Nearly two hundred years ago, Davy observed that sodium when dissolved in liquid ammonia produces a deep blue solution.^{2} Several years later, Weyl described this phenomenon in detail and attributed it to the formation of NaNH_{3} chemical complexes.^{3} Kraus studied the properties of various metal ammonia solutions and surmised that the blue color was due to solvated ammoniated electrons.^{4} Since then, a number of studies on solvated ammoninated electrons have appeared in the literature.^{5–22} Several groups^{23–26} studied ammonia cluster anions experimentally. Haberland and co-workers^{23} concluded that the smallest ammonia cluster size to bind an electron in high Rydberg states is $(NH3)35$. Sarkas and co-workers^{2} used photo-electron spectroscopy to study $(NH3)n\u2212$ cluster anions and concluded that the vertical detachment energies (VDE) for the cluster anions in the size range *n* = 41 to *n* = 1100 increase smoothly from 0.55 to 1.05 eV. By extrapolating the measured VDEs to the VDE(∞) limit, they surmised that the cluster anions exist as embryonic forms of bulk ammoniated electrons. However, Zewail and co-workers^{27} were able to produce ammonia cluster anions as small as *n* = 13.

A number of theoretical studies have been carried out on the solvated electron in the ammonia clusters.^{7–10} Path integral Monte Carlo studies by Klein and co-workers have shown that in *n* = 16 cluster the electron binding energy at 100 K is less than k_{B}T. This study also showed that in the smaller clusters the excess electron is on the surface. Sommerfeld^{11} investigated excess electrons bound to the small ammonia clusters $(NH3)n\u2212(n=2\u22124)$ using the coupled cluster method. He showed that the smallest ammonia cluster able to form a dipole-bound anion is the trimer, and its VDE is predicted to be 27 meV, somewhat smaller than that of the water dimer. Baranyi and Turi^{28} studied smaller size ammonia anions using quantum chemical and molecular dynamics methods. They employed MP2, density functional methods, and coupled cluster singles and doubles with perturbative triples [CCSD(T)] to study the excess electron binding in small clusters of size *n* = 3–8. They concluded that the excess electrons in these clusters are dipole bound with increasing VDEs as the cluster size grows. They also noted that the electron binding energy in these clusters is less than those in the water clusters. Based on the calculations on VDEs of the relaxed anionic clusters, they inferred that the clusters in the size range *n* = 12–32 weakly bind the excess electron with binding energies on the order 100 meV. In a subsequent study, Baranyi and Turi^{12} used finite temperature molecular dynamics simulations to study excess electron solvation dynamics in $(NH3)n\u2212$ ammonia clusters in the *n* = 8–32 size range. Based on their results, they concluded that the $(NH3)24\u2212$ and $(NH3)32\u2212$ can accommodate the excess electron, which is mainly localized on the surface but may be partly embedded in the cluster. Very recently, Restrepo and co-workers^{29} have studied microsolvation of electrons in ammonia clusters. They noted that the smallest cluster that can bind an electron is ammonia dimer.

The majority of the quantum chemical calculations on clusters are nowadays performed using the Kohn–Sham formulation of the density functional theory which requires approximation to the exchange–correlation (XC) functional. Standard approximations used for modeling the XC functionals such as various generalized gradient approximations (GGAs) or *meta*-GGAs are not accurate enough to describe weakly bound anions such as the anions of small ammonia clusters, the topic of the present study. In fact, it is well-known that these approximations cannot even describe atomic anions due to the self-interaction-error (SIE). With density functional approximations (DFAs), the eigenvalue of the highest occupied molecular orbital (HOMO) turns out to be positive for systems with excess electrons indicating an unbound electron. A more reliable indicator of excess electron binding with DFAs is the total energy. For stable anions, the total energy of anion is lower than the total energy of the neutral cluster. However, one issue with local or semilocal functionals is the excessive binding energy when the total energy differences are calculated. For these reasons, most computational studies on small ammonia cluster anions^{9,28} have used post-Hartree–Fock (HF) methods such as the Møller–Plesset perturbation method (MP2) and CCSD(T). CCSD(T) is computationally very demanding and has only been used to study small clusters, and often in these studies, structural relaxation is carried out using a faster method such as MP2. The SIE in the DFAs can be mitigated by means of mixing some percent of the Hartree–Fock exchange using various criteria. The accuracy of such hybrid approaches in predicting VDEs then can depend on the percentage of the Hartree–Fock exchange.

Alternatively, the SIE can be removed on an orbital-by-orbital basis using self-interaction-correction (SIC) methods. The most widely known of these methods is the Perdew–Zunger SIC (PZSIC) method.^{30} Recently, Zope and co-workers^{31,32} proposed a new one-electron SIC method, called the locally scaled self-interaction-correction (LSIC) method, that scales the energy densities of the orbital-wise energy corrections.^{31} The LSIC method when used with the simplest DFA, the local spin density approximation (LSDA), has been found to provide an accurate description of both the equilibrium and stretched bond properties.^{33–39} In this work, we use both the PZSIC and LSIC methods to study the weakly bound anions of ammonia clusters. The PZSIC method is used in conjunctions with three local and semilocal functionals, viz., the LSDA functional, the Perdew–Burke–Ernzerhof (PBE) GGA, and the SCAN *meta*-GGA functionals, while the LSIC is used with the LSDA and PBE-GGA functionals. This is the first application of the LSIC method to study weakly bound anions. Our results show that the LSIC predicts VDE of ammonia clusters accurately, especially when estimated from the eigenvalue of the highest occupied orbital wherein the error with respect to the benchmark CCSD(T) is only 11 meV. In general, the performance of SIC methods in predicting VDE is shown to be good. We assess their performance by comparing the calculated VDEs of ammonia clusters with the accurate CCSD(T) results.

## II. THEORETICAL METHODS AND COMPUTATIONAL DETAILS

^{30}from the LSDA functional on an orbital-by-orbital as follows:

*ρ*

_{iσ}are single orbital densities and

*σ*=

*α*or

*β*denote the spin. The self-XC energy is

^{40}Removal of SIE provides better predictions of several properties than the LSDA functional such as binding of atomic anions, chemical reaction barrier heights, polarizabilities, and ionization potentials. Unfortunately, the applications of the PZSIC do not provide universal improvements resulting in the paradoxical behavior of the PZSIC. An interested reader can find more details about the paradox of SIC in Refs. 41–43. The PZSIC corrects the tendency of DFA to overstabilize stretched bonds but also worsens the description of the uniform electron limit of the DFAs that are exact in this limit.

^{42}As a result, the PZSIC works well for the one-electron density regions such as systems with stretched bonds but fails in the many-electron regions such as systems in equilibrium. Recently, an alternative to the PZSIC method to correct for one-electron SIEs called the LSIC method was proposed.

^{31,32}The method applies PZSIC correction locally in space and uses the iso-orbital indicator to assess the size of SIC correction as follows:

*et al.*

^{31}Here, the scaling factor $z\sigma (r\u20d7)=\tau \sigma W(r\u20d7)\tau \sigma (r\u20d7)$ lies between 0 and 1 and indicates the nature of the charge density at $r\u20d7$.

*z*

_{σ}= 1 corresponds to one-electron density and

*z*

_{σ}= 0 corresponds to the uniform density where the gradient of density vanishes. Scaling down the self-interaction correction terms with

*z*

_{σ}thus retains the full correction for one-electron densities, making the theory exact in that limit, and eliminates the correction in the limit of a uniform density where $EXCDFA$ is exact.

^{44}Thus, all calculations herein are performed using the FLOSIC method.

^{45,46}The Fermi orbitals are first constructed as

*i*/

*j*and

*σ*are the orbital and spin indices, respectively, and

*ρ*is the electron density.

*ψ*

_{j,σ}’s are any set of orbitals that span the occupied space. In practice, canonical orbitals are used. $a\u20d7i$ is a position in space called a Fermi orbital descriptor (FOD). One such FOD is required for each Fermi orbital. The Fermi orbitals resulting from Eq. (7) are normalized but not mutually orthogonal. They are orthogonalized using the Löwdin symmetric orthogonalization scheme,

^{47,48}yielding the FLOs. The SIC energy is determined by minimizing the FOD forces

^{46,49}using the conjugate-gradient and limited memory BFGS algorithm.

^{50}The self-consistent PZSIC (with FLOs) calculations can be performed either using optimized effective potential within the Kriger–Li–Ifarate approximation

^{51}or using the Jacobi update approach.

^{52}The self-consistent LSIC calculations are performed by ignoring the variation of the scaling factor as explained in Ref. 34. This approach is equivalent to applying the local scaling correction to the potential instead of energy density. This approach can be referred to as quasi LSIC.

The VDE is defined as the energy difference between the neutral cluster and the anionic cluster evaluated at the anionic geometry. The assumption here is that the excess electron removal is a fast process without any structural relaxation. A positive VDE indicates a bound anionic system.

All calculations are performed using the FLOSIC code 0.2.^{53,54} We use the default NRLMOL basis^{55,56} with extra long-range polarization functions on the hydrogen and nitrogen atoms. The extra long-range Gaussians are chosen such that one-third$(13)$ of the (N)th exponent for hydrogen and one sixth$(16)$ of the (N)th exponent for nitrogen is the (N + 1)th exponent for the long-range Gaussians on the respective atoms. The basis sets are expanded by adding two extra contracted s- and p-type functions, one d-type diffuse functions on the nitrogen atoms and two s-, p-, and d-type functions on the hydrogen atoms. For a given atom, the extra s-, p-, and d-type functions share the same extra long-range exponents. This is done until the total energy of the ammonia anion converges. The details of the basis functions are given in the supplementary material.

## III. RESULTS AND DISCUSSION

As mentioned in the Introduction, our goal is to examine how one-electron SIC methods perform in predicting the VDEs of the ammonia clusters. For this purpose, we use the recently reported structures by Baranyi and Turi.^{28} This set of small ammonia cluster anions is shown in Fig. 1. It includes (a) six linear chain-like structures from trimer (*n* = 3) to octamer (*n* = 8), where the ammonia molecules are bound by successive hydrogen bonds and (b) two branched structures for heptamer and octamer isomers. Besides these clusters shown in Fig. 1, we also consider the $NH3\u2212$ for the purpose of comparison.

These hydrogen-bonded clusters possess a weak dipole moment, which is responsible for the extra electron binding to the cluster. In such anions, the extra electron is often located away from the molecular frame. The dipole moment model, however, has its limitations. In Table I, we present calculated dipole moments of the neutral ammonia clusters with DFA, PZSIC-DFA, and LSIC-DFA. The purpose of this comparison is to show how SIC methods affect the dipole moment of these clusters. The experimental electric dipole moment of ammonia in the gas phase is 1.47 D.^{57} The calculated dipole moments of ammonia molecules with DFAs are in good agreement with the experimental value. On the other hand, PZSIC-DFA tends to increase the dipole moment, whereas the LSIC-DFA values are between the DFA and PZSIC-DFA. This trend is similar to what was noted earlier^{39} for small molecules. The incorporation of SIC leads to a relatively deeper potential on nitrogen atoms resulting in larger charge transfer between the hydrogen donor and acceptor molecules. The use of quasi-LSIC reduces the dipole moments compared to the PZSIC since the SIC potential is scaled and, as a result, the charge transfer is also reduced. The linear chains have large dipole moments due to their structure compared to the branched structures.

DFA . | FLOSIC . | LSIC . | ||||||
---|---|---|---|---|---|---|---|---|

Cluster . | LDA . | PBE . | SCAN . | LDA . | PBE . | SCAN . | LDA . | PBE . |

Monomer | 1.54 | 1.49 | 1.53 | 1.69 | 1.66 | 1.68 | 1.55 | 1.66 |

Trimer | 5.09 | 4.91 | 4.95 | 5.3 | 5.25 | 5.25 | 4.89 | 4.86 |

Tetramer | 7.07 | 6.83 | 6.85 | 7.3 | 7.22 | 7.21 | 6.74 | 6.69 |

Pentamer | 9.13 | 8.81 | 8.82 | 9.38 | 9.27 | 9.27 | 8.65 | 8.59 |

Hexamer | 11.14 | 10.74 | 10.75 | 11.38 | 11.26 | 11.24 | 10.49 | 10.44 |

Heptamer | 13.2 | 12.73 | 12.72 | 13.44 | 13.32 | 13.3 | 12.43 | 12.35 |

Heptamer B | 6.33 | 6.11 | 6.15 | 6.62 | 6.5 | ⋯ | 5.99 | 5.94 |

Octamer | 15.28 | 14.74 | 14.71 | 15.53 | 15.33 | 15.35 | 14.36 | 14.21 |

Octamer B | 14.13 | 13.62 | 13.65 | 14.55 | 14.37 | ⋯ | 13.34 | 13.23 |

DFA . | FLOSIC . | LSIC . | ||||||
---|---|---|---|---|---|---|---|---|

Cluster . | LDA . | PBE . | SCAN . | LDA . | PBE . | SCAN . | LDA . | PBE . |

Monomer | 1.54 | 1.49 | 1.53 | 1.69 | 1.66 | 1.68 | 1.55 | 1.66 |

Trimer | 5.09 | 4.91 | 4.95 | 5.3 | 5.25 | 5.25 | 4.89 | 4.86 |

Tetramer | 7.07 | 6.83 | 6.85 | 7.3 | 7.22 | 7.21 | 6.74 | 6.69 |

Pentamer | 9.13 | 8.81 | 8.82 | 9.38 | 9.27 | 9.27 | 8.65 | 8.59 |

Hexamer | 11.14 | 10.74 | 10.75 | 11.38 | 11.26 | 11.24 | 10.49 | 10.44 |

Heptamer | 13.2 | 12.73 | 12.72 | 13.44 | 13.32 | 13.3 | 12.43 | 12.35 |

Heptamer B | 6.33 | 6.11 | 6.15 | 6.62 | 6.5 | ⋯ | 5.99 | 5.94 |

Octamer | 15.28 | 14.74 | 14.71 | 15.53 | 15.33 | 15.35 | 14.36 | 14.21 |

Octamer B | 14.13 | 13.62 | 13.65 | 14.55 | 14.37 | ⋯ | 13.34 | 13.23 |

The density differences between the anion and neutral states for the six linear chain molecules are shown in Fig. 2. The blue (red) color stands for the positive (negative) sign and shows the excess (deficit) density in the anion. From this plot, we see that the blue color cloud becomes more widespread toward where the extra charge is mostly located. The extra electron is located toward the positive end of the dipole. A large percentage of the charge distribution outside of the molecule indicates the diffusive nature of the extra electron in such anion systems and reaffirms the necessity of employing the additional diffusive Gaussian functions.

The VDEs for all clusters calculated using different methods employed in this work are compared with MP2 and CCSD(T) values from the literature in Table II. The CCSD(T) is chosen as the accuracy reference. Our calculations on the monomer showed that both LDA and PBE incorrectly predict positive VDE suggesting that with these functionals the monomer can bind an extra electron. This result is a manifestation of the overbinding tendencies of these functionals. On the other hand, the semilocal self-correlation free SCAN functional yields a negative VDE in agreement with the accurate quantum chemical calculations. We would like to note here that a calculated negative VDE means that the cluster anion is at higher energy with respect to the neutral cluster indicating that the extra electron is not bound in this case. However, the negative VDE itself does not have physical meaning. The discussion below in this paragraph is only to show that the LSIC results follow the CCSD(T) results closer than the PZSIC results. Both SIC-based methods (PZSIC and LSIC) yield the correct sign showing that the monomer is unable to bind an extra electron. While the application of PZSIC with LDA and PBE corrects the overbinding of VDE to −26.5 and −28.2 meV, respectively, LSIC further corrects the results for both DFA functionals to −65.4 and −44.6 meV compared to the reference CCSD(T) value of −67.4 meV. The LSIC-LSDA VDE (−65.4 meV) is in excellent agreement with CCSD(T) VDE (−67.4 meV).

DFA . | PZSIC . | LSIC . | Reference . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Cluster . | LDA . | PBE . | SCAN . | LDA . | PBE . | SCAN . | LDA . | PBE . | MP2 . | CCSD(T) . |

Monomer | 52.2 | 58.4 | −83.8 | −26.5 | −28.2 | −47.6 | −65.4 | −44.6 | −73.6 | −67.4 |

Trimer | 212.4 | 196.4 | 84.5 | 95.0 | 73.2 | 37.3 | 5.8 | 48.2 | 3.2 | 21.7 |

Tetramer | 276.7 | 254.4 | 120.3 | 150.7 | 117.6 | 91.5 | 47.4 | 94.3 | 37.7 | 62.3 |

Pentamer | 327.3 | 299.5 | 158.4 | 190.7 | 160.8 | 129.3 | 73.6 | 132.9 | 67.8 | 97.3 |

Hexamer | 369 | 339.1 | 185.5 | 227.3 | 192.8 | 159.6 | 105.3 | 152.3 | 91.4 | 124.0 |

Heptamer | 400.5 | 368 | 210.7 | 241.3 | 197.3 | 121.6 | 92.0 | 184 | 111.0 | 147.0 |

Heptamer B | 331.8 | 288.7 | 144.7 | 172.5 | 129.8 | ⋯ | 11.1 | 93.6 | 35.1 | 70.3 |

Octamer | 429.6 | 388.9 | 270.6 | 282.7 | 251.6 | 206.3 | 130.7 | 191.5 | 127.0 | 164.0 |

Octamer B | 564.4 | 483.1 | 337.8 | 410.3 | 331.4 | ⋯ | 218.2 | 288.3 | 193.3 | 250.5 |

MAE | 232.7 | 200.8 | 65.8 | 97.1 | 61.8 | 28.6 | 28.3 | 30.1 | 30.8 |

DFA . | PZSIC . | LSIC . | Reference . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Cluster . | LDA . | PBE . | SCAN . | LDA . | PBE . | SCAN . | LDA . | PBE . | MP2 . | CCSD(T) . |

Monomer | 52.2 | 58.4 | −83.8 | −26.5 | −28.2 | −47.6 | −65.4 | −44.6 | −73.6 | −67.4 |

Trimer | 212.4 | 196.4 | 84.5 | 95.0 | 73.2 | 37.3 | 5.8 | 48.2 | 3.2 | 21.7 |

Tetramer | 276.7 | 254.4 | 120.3 | 150.7 | 117.6 | 91.5 | 47.4 | 94.3 | 37.7 | 62.3 |

Pentamer | 327.3 | 299.5 | 158.4 | 190.7 | 160.8 | 129.3 | 73.6 | 132.9 | 67.8 | 97.3 |

Hexamer | 369 | 339.1 | 185.5 | 227.3 | 192.8 | 159.6 | 105.3 | 152.3 | 91.4 | 124.0 |

Heptamer | 400.5 | 368 | 210.7 | 241.3 | 197.3 | 121.6 | 92.0 | 184 | 111.0 | 147.0 |

Heptamer B | 331.8 | 288.7 | 144.7 | 172.5 | 129.8 | ⋯ | 11.1 | 93.6 | 35.1 | 70.3 |

Octamer | 429.6 | 388.9 | 270.6 | 282.7 | 251.6 | 206.3 | 130.7 | 191.5 | 127.0 | 164.0 |

Octamer B | 564.4 | 483.1 | 337.8 | 410.3 | 331.4 | ⋯ | 218.2 | 288.3 | 193.3 | 250.5 |

MAE | 232.7 | 200.8 | 65.8 | 97.1 | 61.8 | 28.6 | 28.3 | 30.1 | 30.8 |

The tendency of LDA and PBE to overbind an electron can also be seen in the rest of the systems in this study. The main results of the benchmark for the linear structures from Fig. 1(a) are summarized in Fig. 3. To facilitate the discussion of the results, we will refer to the VDE obtained from the total energy differences of neutral and anionic clusters as ΔSCF VDE. Figure 3 includes the deviation from the reference value $ECCSD(T)VDE$ of (a) ΔSCF VDE, (b) VDE estimated from the absolute of the HOMO eigenvalue, and (c) the mean average error (MAE) of both. For convenience, here, we introduce two shorthand notations. The deviations of VDE and HOMO energies are defined as $\Delta EKVDE=EKVDE\u2212ECCSD(T)VDE$ and $\Delta \epsilon KHOMO=\u2212\epsilon KHOMO+ECCSD(T)VDE$, respectively, where $EKVDE$ is the VDE energy calculated with the method *K* specified in the legend section in Fig. 3. VDE values for the reference method MP2 are from the total energy differences of anion and neutral clusters.

Even though the LDA and PBE predict the correct sign, the ΔSCF VDE values with these functionals are generally significantly larger than the more accurate quantum chemical methods as can be seen from the mean absolute errors for these functionals. With the DFAs, the error in VDEs is reduced when one moves up on the Perdew–Schmidt Jacob’s ladder from LDA to GGA to *meta*-GGA. The errors of the results for SCAN are much smaller and do not scale strongly with the system size as LDA and GGA do. By explicitly removing the self-interaction in an orbital-by-orbital fashion, a substantial reduction of error for DFAs at different rungs of the Jacob’s ladder can be seen from the results for DFA-PZSIC methods. The PZSIC-DFA MAE values range down from 97.1 meV for PZSIC-LDA to 28.6 meV for PZSIC-SCAN. The PZSIC results show that the removal of SIE corrects the overbinding tendency of the DFAs.

The best agreements are observed for the LSIC methods, which scale down the SIC energy density at a given position using a local scaling factor. This is possibly because LSIC is well-defined for the LDA functional with no gauge dependency.^{58} Our earlier experience shows that LSIC also often performs well with PBE functional^{33} even though the formal gauge dependency occurs. We applied LSIC with both LDA and PBE functionals in this study. The ΔSCF VDE values with LSIC are in excellent agreement with the reference values in all cases with MAE values of 28.3 and 30.1 meV for LSIC-LDA and LSIC-PBE, respectively. These results are comparable with, and in most cases better than, MP2 which has a MAE of 30.8 meV.

In exact DFT, the highest occupied eigenvalue equals the negative of the ionization potential.^{59–62} Thus, for an anionic cluster, the negative of the eigenvalue of the HOMO is the VDE of electron. We also estimate the VDE of ammonia anions from the HOMO eigenvalues, which are presented in Table III for the linear systems in Fig. 1(a). With all three semilocal DFAs, the HOMO eigenvalues fail in predicting the binding of the excess electron. Not only do they deviate by a large amount, they do not even yield the correct sign, which indicates that the HOMOs are unbound orbitals in the complete basis set limit. The positive eigenvalue is a result of SIE in these functionals, which makes the asymptotic potential seen by the electron too shallow. On the other hand, the results of HOMO energies for both groups of SIC-based methods compare well with the reference values with correct sign. The MAE for the HOMO eigenvalues with PZSIC-PBE and PZSIC-SCAN show excellent agreement with the reference values. Figure 3(b) shows the HOMO energy deviation from the reference values $ECCSD(T)VDE$ for PZSIC-DFA, LSIC-DFA, and $EMP2VDE$. It shows that LSIC-LDA performs the best; PZSIC-PBE comes second, both of which are significantly better than MP2 in most cases. This improvement is due to the removal of self-interaction which yields improved asymptotic description of the potential resulting in negative (bound) eigenvalues for the extra electron in all SIC methods.

DFA . | PZSIC . | LSIC . | Reference . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Cluster . | LDA . | PBE . | SCAN . | LDA . | PBE . | SCAN . | LDA . | PBE . | MP2 . | CCSD(T) . |

Monomer | −492.1 | −448.8 | −467.9 | −42.1 | −50.9 | −50.3 | −48.3 | −59.7 | −73.6 | −67.4 |

Trimer | −462.6 | −437.6 | −394.9 | 63.8 | 36.4 | 39.3 | 38.7 | 15.8 | 3.2 | 21.7 |

Tetramer | −425.6 | −401.7 | −493.4 | 111.8 | 72.5 | 80.2 | 75.1 | 46.1 | 37.7 | 62.3 |

Pentamer | −404.2 | −376.4 | −488.5 | 146.4 | 106.6 | 121.2 | 106.1 | 72.1 | 67.8 | 97.3 |

Hexamer | −372.1 | −346.4 | −451.6 | 177.2 | 130.8 | 142.8 | 129.4 | 92.8 | 91.4 | 124 |

Heptamer | −361.2 | −331.2 | −451.3 | 200 | 156.7 | 164.3 | 150.2 | 110.7 | 111 | 147 |

Heptamer B | −400.4 | −376.2 | −468.2 | 128.1 | 79.9 | ⋯ | 94.9 | 66.01 | 35.1 | 70.3 |

Octamer | −438.7 | −421.9 | −257 | 223.2 | 171.4 | 184.6 | 170.5 | 127.5 | 127 | 164 |

Octamer B | −493.7 | −328.3 | −575.1 | 340.3 | 253.1 | ⋯ | 250.6 | 179.7 | 193.3 | 250.5 |

MAE | 524.5 | 482.0 | 546.4 | 53.2 | 9.6 | 19.0 | 10.8 | 26.0 | 30.8 |

DFA . | PZSIC . | LSIC . | Reference . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Cluster . | LDA . | PBE . | SCAN . | LDA . | PBE . | SCAN . | LDA . | PBE . | MP2 . | CCSD(T) . |

Monomer | −492.1 | −448.8 | −467.9 | −42.1 | −50.9 | −50.3 | −48.3 | −59.7 | −73.6 | −67.4 |

Trimer | −462.6 | −437.6 | −394.9 | 63.8 | 36.4 | 39.3 | 38.7 | 15.8 | 3.2 | 21.7 |

Tetramer | −425.6 | −401.7 | −493.4 | 111.8 | 72.5 | 80.2 | 75.1 | 46.1 | 37.7 | 62.3 |

Pentamer | −404.2 | −376.4 | −488.5 | 146.4 | 106.6 | 121.2 | 106.1 | 72.1 | 67.8 | 97.3 |

Hexamer | −372.1 | −346.4 | −451.6 | 177.2 | 130.8 | 142.8 | 129.4 | 92.8 | 91.4 | 124 |

Heptamer | −361.2 | −331.2 | −451.3 | 200 | 156.7 | 164.3 | 150.2 | 110.7 | 111 | 147 |

Heptamer B | −400.4 | −376.2 | −468.2 | 128.1 | 79.9 | ⋯ | 94.9 | 66.01 | 35.1 | 70.3 |

Octamer | −438.7 | −421.9 | −257 | 223.2 | 171.4 | 184.6 | 170.5 | 127.5 | 127 | 164 |

Octamer B | −493.7 | −328.3 | −575.1 | 340.3 | 253.1 | ⋯ | 250.6 | 179.7 | 193.3 | 250.5 |

MAE | 524.5 | 482.0 | 546.4 | 53.2 | 9.6 | 19.0 | 10.8 | 26.0 | 30.8 |

A side-by-side comparison of MAE of both ΔSCF VDEs and the VDE estimates from the absolute HOMO eigenvalues is shown in Fig. 3(c). This comparison highlights that with the same SIC functional, the VDE estimates from the HOMO eigenvalues are closer to the reference values than the VDEs. LSIC-LDA consistently performs the best in both VDE estimates from the absolute HOMO eigenvalues and ΔSCF VDEs. Both LSIC-PBE and PZSIC-SCAN offer reasonably accurate predictions that are comparable to MP2. On the other hand, PZSIC-PBE although performing well for the VDE estimates from the absolute HOMO eigenvalues does rather poorly for the ΔSCF VDEs.

We also considered additional geometries as depicted in Fig. 1(b) of the heptamer and octamer clusters to examine the VDE dependence on the geometries. These structures are referred to as branched systems. The results of both ΔSCF VDE and VDE from the absolute HOMO eigenvalues for the branched clusters are also included in Tables II and III labeled as heptamer B and octamer B, respectively. The trends seen in the linear clusters that FLOSIC-DFA improves over DFAs in predicting the excess electron binding also holds for the branched structures too. The LSIC-PBE consistently underestimates the VDEs for all clusters including the branched isomers, whereas LSIC-LDA consistently overestimates. There is no notable difference in the errors of VDEs as well as for the HOMO eigenvalues for linear and branched clusters.

After the completion of present work, another study by Restrepo and co-workers on the microsolvation of electrons in ammonia clusters was reported.^{29} This work shows that the ammonia dimer can also bind an extra electron. We, therefore, also considered the dimer geometry that binds an electron from this work and computed vertical detachment energy using the PZSIC-LSDA and LSIC-LSDA. The VDE from the total energy differences and as negative of the HOMO eigenvalue in the PZSIC-LSDA are 50.3 and 21.6 meV compared to CCSD(T) value of 2 meV by Restrepo and co-workers. Thus, PZSIC overestimates the VDE which is consistent with the above-mentioned results for other cluster anions. Likewise, the LSIC-LSDA significantly improves over the PZSIC-LSDA with VDE from the total energy differences and as negative of the HOMO eigenvalue of 5.3 and 5.9 meV compared to the CCSD(T) value of 2 meV. This further strengthens the primary conclusion of the work that the LSIC method even with the simplest density functional such as LSDA provides good estimates of VDE of ammonia cluster anions.

## IV. CONCLUSIONS

We have investigated the role of SIEs on the VDEs of ammonia clusters for three rungs of non-empirical semilocal functionals, viz., LDA, PBE-GGA, and SCAN *meta*-GGA using two one-electron self-interaction methods, PZSIC and LSIC. The LSIC calculations are carried out with LDA and PBE functionals only due to stronger gauge dependency in the case of SCAN *meta*-GGA. The results show that the VDEs of the excess electron in ammonia clusters is not reliably predicted with the semilocal density functionals. All three functionals exhibit overbinding tendencies though it is reduced for the SCAN functional. As expected, uncorrected DFAs perform very poorly for the VDE estimates from the HOMO eigenvalues. The PZSIC method improves the accuracy of the ΔSCF VDEs as well as the VDE estimates as absolute of HOMO eigenvalues. The most notable feature of this work is that the LSIC leads to a substantial improvement over the PZSIC in predicting the VDE from the total energy differences for the LDA and PBE functionals. The MAE of the calculated VDE values with PZSIC-SCAN, LSIC-LDA, and LSIC-PBE are comparable to those with the MP2 method. The prediction of VDE of weakly bound anions especially from the eigenvalue of the HOMO is a challenging test for density functional approximations even with the Hartree–Fock component. In this regard, it is remarkable that the LSIC-LDA (MAE = 10.8 meV) and the PZSIC-PBE (MAE = 9.6 meV) provide excellent estimates of the electron removal energy as absolute of the HOMO eigenvalue. The absolute of the HOMO eigenvalues of PZSIC-PBE, PZSIC-SCAN, and LSIC-LDA provide more accurate estimates of the VDE energies than the ΔSCF VDEs of the corresponding functionals. These results combined with earlier results on water cluster anions^{33} show that the removal of self-interaction errors using one-electron self-interaction theories can provide accurate VDEs (better than the MP2 method) for hydrogen-bonded anionic systems.

## SUPPLEMENTARY MATERIAL

The supplementary material contains details of the basis sets and the optimized Fermi orbital descriptor positions.

## ACKNOWLEDGMENTS

This work was supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, as part of the Computational Chemical Sciences Program under Award No. DE-SC0018331. Support for computational time at the Texas Advanced Computing Center and at NERSC is acknowledged. We thank Mr. Prakash Mishra for assistance with calculations on dimer anions during the revision version of the manuscript.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Peter Ufondu**: Data curation (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Visualization (equal); Writing – original draft (equal). **Po-Hao Chang**: Investigation (equal); Software (equal); Writing – original draft (equal). **Tunna Baruah**: Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Rajendra R. Zope**: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

## REFERENCES

_{3})

_{n = 41–1100}

_{3}

*Small Particles and Inorganic Clusters*

*Advances in Atomic, Molecular, and Optical Physics*

*Z*atoms?