THE JOURNAL OF CHEMICAL PHYSICS 139, 194102 (2013)
Self-interaction corrected density functional calculations of molecular Rydberg states Hildur Gudmundsdóttir,1 Yao Zhang,2 Peter M. Weber,2 and Hannes Jónsson2,3 1
Science Institute of the University of Iceland, 107 Reykjavík, Iceland Department of Chemistry, Brown University, Providence, Rhode Island 02912, USA 3 Faculty of Physical Sciences, VR-III, University of Iceland, 107 Reykjavík, Iceland 2
(Received 11 September 2013; accepted 28 October 2013; published online 15 November 2013) A method is presented for calculating the wave function and energy of Rydberg excited states of molecules. A good estimate of the Rydberg state orbital is obtained using ground state density functional theory including Perdew-Zunger self-interaction correction and an optimized effective potential. The total energy of the excited molecule is obtained using the Delta Self-Consistent Field method where an electron is removed from the highest occupied orbital and placed in the Rydberg orbital. Results are presented for the first few Rydberg states of NH3 , H2 O, H2 CO, C2 H4 , and N(CH3 )3 . The mean absolute error in the energy of the 33 molecular Rydberg states presented here is 0.18 eV. The orbitals are represented on a real space grid, avoiding the dependence on diffuse atomic basis sets. As in standard density functional theory calculations, the computational effort scales as NM2 where N is the number of orbitals and M is the number of grid points included in the calculation. Due to the slow scaling of the computational effort with system size and the high level of parallelism in the real space grid approach, the method presented here makes it possible to estimate Rydberg electron binding energy in large molecules. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829539] I. INTRODUCTION
Molecular Rydberg states, i.e., highly excited electronic states of molecules where one electron is promoted to a hydrogen-like orbital, are profoundly important for theoretical and practical reasons.1 Conceptually, they can be thought of as a positively charged molecular core that is surrounded by an electron in a hydrogen-like electronic state, typically with principal quantum number n ≥ 3. The interactions of the loosely bound electron with the electronic, vibrational, and rotational states of the ion core can be quite complex. But when considering only the electron binding energy, i.e., the energy difference between the molecule in the Rydbergexcited state and the corresponding ion state (with the same quantum numbers in the other degrees of freedom), then many of these complexities cause negligible energy changes. Experimentally, it has been found that even in large molecules the electron binding energy of a Rydberg state is a well-defined quantity that, in spectra of modest resolution, is independent of vibrational and/or rotational excitation.2 Even so, the binding energy has proven to be quite sensitive to the molecular structure.3 This has led to the recent development of spectroscopic tools that use Rydberg states for the identification of molecular structures and possibly of molecules.4–10 For many practical applications, the determination of molecular structures remains an important challenge. Examples include situations when samples are extremely dilute or where rapid throughput is of the essence. There is thus a need for additional structure determination methods such as one relying on Rydberg states. We draw attention to a parallel but much more developed spectroscopic method, nuclear magnetic resonance
0021-9606/2013/139(19)/194102/8/$30.00
(NMR), which also uses the spectroscopic determination of energy differences between molecular quantum states to determine molecular structures. While fundamentally, energy differences are measured, chemists have learned to associate these energy differences with molecular structure, turning NMR spectroscopy into a powerful tool for molecular structure determination. Although there are only relatively few spectral lines, the information content is so decisive that even complicated structures of proteins can now be elucidated. It stands to reason that an analogous approach to molecular structure determination could be based on the measurement of binding energy of electrons in Rydberg states. Experiments have already shown that the electron binding energy is remarkably sensitive to the molecular structure. But to date, the method has been applied in a comparative, rather than an absolute fashion, by identifying, but not determining, molecular structures. Consequently, the technique presently provides a fingerprint of a molecular structure and is referred to as Rydberg Fingerprint Spectroscopy (RFS). The method has proven particularly useful in a photoionization scheme when implemented in a time-resolved way, because the kinetics and dynamics of even complex molecular reactions can be followed in real time. To advance the method from a fingerprint to a structure determination tool, a technique to invert the spectra or otherwise derive the structure from the measured data is required. The aim of the work reported here is to develop a fast method for computing the energy of Rydberg states of large molecules. Rydberg states of molecules have been successfully modeled with high level wave function methods such as CASPT2, but as these methods scale poorly with system size
139, 194102-1
© 2013 AIP Publishing LLC
194102-2
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
they are only suitable for small molecules. Density functional theory (DFT)11, 12 has become a powerful tool for describing the ground state electronic structure of atoms, molecules, and solids13 and with time dependent DFT (TDDFT) it is even possible to calculate the properties of excited electronic states. Rydberg states of molecules have a very large spatial spread compared to the ground state density, which makes their simulation particularly problematic with traditional DFT. The effective potential of an electron, described by a DFT functional of the local density approximation (LDA) or generalized gradient approximation (GGA) type, goes to zero much too quickly due to self-interaction error. While this is frequently not a problem for calculating ground state energy, it means it is impossible to describe Rydberg orbitals with traditional DFT, as the effective potential is far too positive in their region of space. In the present work, we use ground state DFT with the Perdew-Zunger self-interaction correction (SIC) to calculate Rydberg orbitals of the molecules NH3 , H2 O, H2 CO, C2 H4 , and N(CH3 )3 . The SIC gives the correct long range behavior of the potential and as a result the virtual orbitals in the calculations give good estimates of the Rydberg orbitals of the molecules. The Rydberg excitation energy, or the binding energy of Rydberg electrons, is then calculated by placing an electron in the virtual Rydberg orbitals by means of the Delta Self-Consistent Field29 (SCF) method. II. SELF-INTERACTION CORRECTION
In Kohn-Sham (KS) DFT the ground state of a molecule, with an electron density ρ(r), is found by solving the eigenvalue problem 1 − ∇ 2 + vext (r) + vH (r) + vxc (r) ψi = εi ψi , (1) 2 where vext is the external potential of the nuclei, vH is the Hartree potential, accounting for the Coulomb interaction of the charge density with itself, and vxc is the exchangecorrelation (XC) potential. The occupied KS orbitals {ψ i } must span the ground state density: N
|ψi (r)|2 = ρ(r)
and possibly the gradient, of the total density. Explicit density functionals allow fast computation of the ground state energy of molecules to a good approximation. However, the Hartree energy of Eq. (3) not only contains the interaction of each electron with all the others, it also contains the electron’s interaction with itself. This self-Hartree-interaction is only partly canceled in local, explicit density functionals, leaving a spurious self-interaction error. The self-interaction in DFT gives rise to several problems, such as the delocalization of charge, the underestimation of energy barriers and the overestimation of bond energy.17–19 Another manifestation of the electron selfinteraction is the incorrect asymptotic behavior of the XC potential, which makes the description of Rydberg states particularly difficult. At large distances from the atom the correct Kohn-Sham potential has the form 1 lim vxc (r) = − , r
r→∞
(4)
whereas the potential of the explicit, local density functionals goes to zero much faster, as shown in Fig. 1. For example, the LDA potential decays exponentially. The shallow potential of the LDA and GGA functionals leads to an underestimation of the binding energy of electrons, particularly in excited orbitals. The diffuse and extended orbitals characteristic of Rydberg excitations are typically not found to correspond to bound states at this level of approximation. Various approaches have been taken to correct the longrange behavior of the XC potential in DFT in order to better describe Rydberg states. The most commonly used hybrid functionals, like B3LYP, use fractional Hartree-Fock (HF) exchange, so although they may remove some self-interaction problems they do not have the correct asymptotic potential. The Minnesota functional M06-HF is a hybrid functional designed to describe both ground and excited states using full HF exchange and fitting parameters optimized against extensive data.20 Hybrid functionals with the correct long-range potential have also been designed using the long-range correction (LC) scheme, which uses full HF exchange to account for long-range orbital-orbital interactions, dividing the electron
(2)
i
and are introduced as a means to estimate of the kinetic energy. Orbitals are also useful for the approximation of singleparticle properties, and as discussed above, we would like to use the KS orbitals to estimate the binding energy of electrons in a molecule. The usefulness of the orbital concept for describing single particles is very sensitive to the quality of the exchange-correlation potential, vxc . The electron-electron interaction energy is a sum of the Hartree energy ρ(r)ρ(r ) 1 3 (3) EH [ρ] = d rρ(r)vH (r) = d 3 rd 3 r 2 |r − r | and the exchange-correlation energy, Exc , commonly given by an explicit, local density functional such as the LDA12 or GGA.14–16 These are explicit functionals of the total density and depend only on local information about the value,
FIG. 1. Effective potential of an electron in a trimethylamine molecule (TMA) calculated along the direction of one of the N–C bonds. Dashed line: The −1/r potential centered on the central N-atom in the molecule. Solid lines: LDA (blue) and LDA-SIC (red). By applying the self-interaction correction, the correct long range form of the potential is obtained.
194102-3
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
repulsion operator into short-range and long-range parts.21 As both of these approaches require the calculation of exact exchange, they scale less favorably than DFT with system size, but they greatly improve the accuracy of TDDFT calculations of Rydberg states. A simpler, more efficient approach, is to simply give the XC potential the correct behavior in the region important for Rydberg states. The asymptotic correction (AC) to the Kohn-Sham equations gives the required −1/r behavior to the potential beyond a certain distance from the atoms, but retains the conventional DFT form in energetically important regions.22–24 The idea is that the potential is reasonably well represented by conventional DFT functionals in the region of high electron density, close to the atoms, so the AC is only applied to the potential in the outer region. The shiftand-splice AC of Casida and Salahub23 shifts the potential in the high-density region down by a constant, to account for the derivative discontinuity in the exact functional, making the orbital energy of the highest occupied molecular orbital (HOMO) roughly equal to the negative ionization potential of the molecule, and splices the DFT potential to a potential with the correct −1/r asymptotic behavior for large r. The very similar AC of Tozer and Handy22 does not change the form of the potential in the high-density region but shifts the large r potential up to match the DFT potential. The AC can give quite good estimates, for example, Tozer and Handy calculated the excitation energy of CO, N2 , H2 CO, and C2 H4 with mean absolute errors of 0.4, 0.43, 0.24, and 0.05 eV, respectively, for Rydberg excited states.22 So without actually tackling the cause of the problem, the AC estimates the Rydberg excitation energy by covering up the symptoms of an inaccurate functional. In this paper, we investigate whether the removal of orbital self-interaction makes DFT capable of describing Rydberg excited states, using the Perdew-Zunger25 SIC. The derivative discontinuity in the XC potential, as well as the correct long-range form, are obtained when the SIC is applied, making the functional ideal for calculations of Rydberg states. The SIC has previously been used with an optimized effective potential to calculate the energy of the lowest unoccupied molecular orbital (LUMO) for small molecules,26 but to the authors’ knowledge it has not been used to study Rydberg states. An orbital-by-orbital correction is made to the XC energy SIC [{ρi }] = E0 [ρ] − Exc
N
(EH [ρi ] + Exc [ρi ]).
(5)
i=1
We use a simplified notation here neglecting spin. The orbital densities, ρ i = |ψ i |2 , correspond to the occupied orbitals and E0 [ρ] is an explicit, local density functional, such as the LDA or GGA. The correction can in principle be applied to any such exchange-correlation functional, but one should bear in mind that for many-electron systems it is just an approximation. The magnitude of the many-electron self-interaction error is in general not given by the sum of the individual terms of Eq. (5).27 Recent studies of the SIC have found that the accuracy of the correction is dependent on the type of functional it is applied to, that it can greatly improve DFT estimates
of the total energy of atoms,18 and that a scaled down version of it can give good estimates of atomization energies of molecules.19 Furthermore, SIC gives orbital energies that are good estimates of electron binding energy,18 with the HOMO energy giving a good estimate of the ionization potential. The effective potential for the trimethylamine molecule (TMA) is illustrated in Fig. 1, calculated with LDA with and without SIC. It is clear that the effective potential of the LDA goes to zero much too quickly, whereas applying the SIC gives an effective potential with the correct −1/r asymptotic form. The computational cost for the SIC functional scales with system size like traditional DFT, and the functional does not introduce any empirical parameters. The purpose of the present study is to extend and test the accuracy of the SIC for calculations of Rydberg excited states of molecules. We demonstrate that DFT with the PerdewZunger self-interaction correction, applied within the KriegerLi-Iafrate28 (KLI) approximation of the optimized effective potential (OEP), can give good estimates of excited Rydberg orbitals. The total energy of the excited molecule is then obtained in a SCF29 GGA calculation, where the excited electron is placed in a fixed Rydberg state orbital. We present test calculations using this approach and compare with measurements and high level quantum chemistry calculations of the excitation energy or binding energy of the excited electron. The results show that this approach gives quite good accuracy and can be applied to large molecules. The method can, therefore, be helpful in analyzing and interpreting experimental data on Rydberg states of large molecules and can serve as a helpful tool in the analysis of RFS data. III. METHODOLOGY
The electron density can be spanned by or any set of 2 |ϕ (r)| , and thonormal orbitals {ϕj (r)} such that ρ(r) = N j j any of these sets will result in the same total energy in calculations using LDA or GGA functionals. When SIC is applied, however, the energy depends not only on the total electron density, but also on individual orbital densities. This implies that there exists a specific set of orbitals that minimizes the total energy and the energy is no longer invariant under unitary transformations of the orbitals. It is convenient to carry out the calculation using two sets of orbitals: the canonical Kohn-Sham orbitals {ψ i } that are eigenfunctions of the Kohn-Sham Hamiltonian and solutions to Eq. (1), and the energy optimal orbitals {ϕ j }, which minimize the orbital density dependent total energy, Eq. (5).34 The two basis sets are related by a unitary transformation W Wij ψi . (6) ϕj = i
An orbital density dependent (ODD) energy functional has orbital specific potentials vj (r), whose action is only defined on the corresponding occupied orbital, ϕj (r), with orbital density ρj (r) = |ϕj (r)|2 . For the SIC functional, the orbital density dependent part of the Hamiltonian is given by the ODD potential vjSIC (r) = −vH [ρj ] − vxc [ρj ].
(7)
194102-4
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
The calculation can be brought back to the computationally well-developed realm of Kohn-Sham DFT by using an OEP. This gives a local multiplicative potential that captures the effect of the SIC functional on all states, including virtual states. Otherwise, the SIC would only affect the occupied states. Here, we use the KLI28 approximation to the OEP and iteratively find the XC potential given by ⎫ ⎧ N N ⎬ ⎨
1 KLI KLI vxc (r) = ρj (r)vj (r)+ |ψi (r)|2 v¯xc,i − v¯iSIC , ⎭ ρ(r) ⎩ j =1
i=1
(8) KLI v¯xc,i
=
v¯iSIC =
KLI d 3 r|ψi (r)|2 vxc (r),
N
∗ Wij V¯jSIC k Wik ,
(9)
(10)
of Eq. (2) is replaced by the density ρ(r) =
N−1
|ψi (r)|2 + |ψex (r)|2 ,
(13)
i=1
where the Rydberg orbital is approximated by the excited orbital ψ ex ψex (r) =
M i=N
ψi |ψRyd
i
| ψi |ψRyd |2
1/2 ψi (r),
(14)
which is a linear combination of empty KS orbitals in a calculation with M orbitals. The SCF calculation gives the total energy of the excited state of the system and the excitation energy can then be calculated by subtracting the ground state energy. Sometimes one is interested in the binding energy of the excited electron, since this is what is measured in RFS, and in this case the ground state energy of the cation is subtracted, rather than the energy of the neutral molecule.
j,k=1
where the Hermitian SIC potential matrix is V¯jSIC k = 12 ϕj | vjSIC + vkSIC |ϕk . Using the KLI potential, the orbitals can be obtained from Hˆ SIC ψi = εi ψi ,
(11)
KLI = H + vxc (r) is the where the SIC Hamiltonian H 0 sum of the unitary invariant part Hˆ = − 12 ∇ 2 + vext (r) 0 0 (r), where vxc is the XC potential being cor+ vH (r) + vxc KLI , rected, and the orbital density dependent KLI potential, vxc
ˆ SIC
ˆ0
which contains the effect of the SIC. To minimize the energy with respect to the optimal orbitals, one must simultaneously solve Eq. (11) and find W. This is done by iteratively solving Eq. (11), as in a regular DFT calculation, and in every step of the self-consistent calculation the transformation W is estimated. For a given set of KS orbitals {ψ i }, the estimated W will give the energyoptimal orbitals {ϕ j } and the ODD potentials vi through Eqs. (6) and (7). The error in the energy optimal orbitals, and therefore in W, is given by the antihermitian SIC potential matrix, κ 30 κj k = ϕj | vjSIC − vkSIC |ϕk ,
(12)
which is minimized with respect to the unitary transformation W in each step of the KS iterative process. When κ ≈ 0 the energy optimal orbitals {ϕ j } and the ODD potentials vjSIC , have been found to the chosen level of accuracy, and the KLI potential (10) can be calculated and used in the Kohn-Sham Eq. (11). When self-consistency is reached κ = 0 and the Kohn-Sham equations have been solved. The eigenvalues εi turn out to give a rough approximation of the electron binding energy, and the virtual KS orbitals, {ψ i } where i > N, can be used to estimate the excited states of the molecule. DFT is a ground state theory and to calculate the energy of the Rydberg excited state, the linear-expansion SCF method29 is used with the Perdew-Burke-Ernzerhof (PBE) functional.31 Having calculated the virtual orbitals of the molecule with the SIC functional, a Rydberg-type orbital, ψRyd (r), is selected and an electron is removed from the HOMO and placed in this orbital. The ground state density
A. Implementation
The calculations were carried out with the GPAW software,32, 33 extended to include self-consistent SIC calculations.34, 35 The Rydberg orbitals were approximated using the LDA functional36 and SIC. Here, one can choose the spin multiplicity and in several cases it has been found to give better results to calculate the orbitals for the triplet state of the molecule, rather than the singlet, so the calculations presented here are performed for triplet state molecules. This is a particularly good approach when a molecule has no low lying valence excitations, and the LUMO is a Rydberg orbital, as for H2 O, NH3 , and TMA. A uniform grid with a mesh size of ∼0.15 Å was used to represent the valence electrons, with a cubic simulation cell of side length 20 Å for all molecules except TMA, where the side length was 25 Å. The core electrons were kept frozen, and to represent the core region of each atom two methods were used. In the SIC calculations, Hartwigsen-GoedeckerHutter (HGH) pseudopotentials37 were used. The SCF calculations were performed with the projector augmented wave (PAW)38 method. The PAW was not used in the first step because PAW projectors have not yet been developed for the SIC functional. All calculations were performed using the experimental geometry of the molecules.39 One advantage of using a real space grid is avoiding the need to develop special basis functions to represent the Rydberg orbital. If a Gaussian basis set were used, highly diffuse Gaussians would need to be included and a choice made where to center the diffuse basis functions. In the real space grid, no special treatment is needed for the Rydberg orbital since each point in space is equally well represented. To accurately represent the Rydberg state as a linear combination of empty KS states in the SCF calculation, as in Eq. (14), around 40 orbitals were used. Since the SIC does not increase the scaling of the computational effort with the size of the system, the method can be applied to large molecules and even clusters of molecules.40 The computational effort scales as DFT calculations with the parent functional, as NM2 for grid-based implementations of
194102-5
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
TABLE I. Comparison of eigenvalues and calculations of binding energy. The table shows eigenvalues, ε, and negative binding energy (in eV) of the first two Rydberg states of a water and an ethylene molecule in the singlet state. The eigenvalues are calculated with the LDA, PBE, and LDA-SIC functionals. The SIC + SCF method is used to calculate the energy of the excited state, which is given relative to the energy of the cation to estimate the binding energy. The experimental reference is the difference between the measured excitation energy41, 42 and ionization energy (12.62 eV for H2 O and 10.68 eV for C2 H4 ).39 ε LDA-SIC SIC + SCF Experiment
State
ε LDA
ε PBE
H2 O
3s 3py
− 0.93 0.31
− 0.94 0.28
− 6.29 − 4.13
− 5.42 − 3.81
− 5.22 − 3.52
C2 H4
3s 3py
− 0.30 0.17
− 0.33 0.13
− 4.09 − 3.32
− 3.51 − 2.94
− 3.57 − 2.88
TABLE III. Vertical excitation energy for H2 O in eV. The molecule is in the yz plane, with the electric dipole in the z direction. Columns labeled with δ give the difference between calculated and experimental values. Results of CASPT2 calculations are also shown for comparison. Excited state 3s triplet 3s singlet 3py triplet 3py singlet 3px triplet 3pz triplet 3pz singlet 3px singlet Mean absolute error a b
GGA, where N is the number of orbitals and M is the number of grid points. Similar scaling with system size would be obtained with basis set approaches. For example, the computational time for the SIC calculation on 16 Opteron cores (two quad core nodes, 2.8 GHz clock speed) was 6000 s for NH3 (the box was 20 Å long on each side, a total of 1323 grid points included, and 36 orbitals), while the calculation took 34 000 s for TMA (the box was 25 Å long on each side, a total of 1683 grid points included and 54 orbitals). The scaling of the computer time is 5.7, close to the ratio of NM2 which is 54 × 1686 /36 × 1326 = 6.4 for these two calculations. A comparison of various levels of theory for estimating the binding energy of the Rydberg electron is given in Table I. As mentioned earlier, explicit local DFT functionals will give a poor approximation of the binding energy, and Table I shows the improvement with the inclusion of SIC. An even better estimate is provided with the SCF calculations, based on Rydberg orbitals calculated with SIC. IV. RESULTS
Tables II–VI show the calculated excitation energy compared with experimental results. In the case of TMA the binding energy of the excited electron is shown, as this has been measured experimentally. Where the experimental data are available, comparison is made to both singlet and triplet state energy. For all molecules, except TMA, the results are also TABLE II. Vertical excitation energy for NH3 in eV. The rotational axis of the molecule is in the z direction. Columns labeled with δ give the difference between calculated and experimental values. Results of configuration interaction (CI) calculations are also shown for comparison. Excited state
Experimenta SIC + SCF
δ
CIb
a b
Reference 43. Reference 44.
6.39 7.93 8.26
6.15 6.39 7.62 7.98 8.02 8.25
6.14 0.00 6.37 7.86 0.05 7.88 7.87 − 0.01 8.15 0.02
7.0 7.4 8.9 9.1 9.81 9.98 10.01 10.16
7.11 7.38 8.77 8.94 9.62 9.79 9.95 9.83
− 0.05 − 0.11 0.06
CASPT2b
δ
0.11 − 0.02 − 0.13 − 0.16 − 0.19 − 0.19 − 0.06 − 0.33 0.15
7.13 7.50 9.12 9.27 9.87 9.90 9.95 10.15
0.13 0.10 0.22 0.17 0.06 − 0.08 − 0.06 − 0.01 0.10
compared with previous quantum chemical calculations. For NH3 results are compared with configuration interaction (CI) calculations, and for H2 O, H2 CO, and C2 H4 comparison is made with CASPT2 calculations. The calculated excitation energy of the 3s and 3p Rydberg states of NH3 and H2 O is shown in Tables II and III. The lowest lying vertical excitations of these molecules are Rydberg states. The agreement here is good and indicates that the approach presented can be used for interpretation and analysis of Rydberg states. Figure 2 shows the three 3p Rydberg orbitals for H2 O. These are simply unoccupied eigenstates of the SIC Hamiltonian. The 3px and 3py have a regular p-orbital shape, while the 3pz which lies in the direction of the dipole moment is quite a bit distorted. The method presented was also applied to TMA, a larger molecule which has been studied recently by photoionization spectroscopy.5 The calculated and measured values of the Rydberg electron binding energy are compared in Table IV and some of the calculated Rydberg state orbitals are shown in Fig. 3. The method was also applied to two molecules with π * valence excitations lying below the Rydberg states: H2 CO and C2 H4 . The method presented here is based on ground state DFT and is not expected to be able to account for the complexities of valence excitations. For example in C2 H4 , the π * singlet valence excitation is strongly mixed with the Rydberg TABLE IV. Binding energy of the Rydberg state electron for trimethylamine (TMA). The z axis runs through the rotational axis of molecule. The column labeled with δ gives the difference between calculated and experimental values.
δ
− 0.02
δ
Reference 41. Reference 45.
Excited state 3s triplet 3s singlet 3pxy triplet 3pxy singlet 3pz triplet 3pz singlet Mean absolute error
Experimenta SIC + SCF
3s triplet 3s singlet 3pxy triplet 3pxy singlet 3pz triplet 3pz singlet Mean absolute error a
Reference 5.
Experimenta
3.09 2.25 2.20
SIC + SCF 3.35 3.06 2.59 2.40 2.61 2.29
δ
− 0.03 0.15 0.09 0.10
194102-6
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
TABLE V. Vertical excitation energy of H2 CO in eV. The molecule is in the yz plane, with the CO bond in the z direction. Columns labeled with δ give the difference between calculated and experimental values. While there is a systematic over binding in the SIC calculation, the relative energy of the excited states is quite well reproduced. Results of CASPT2 calculations are also shown for comparison. Excited state
Experimenta
π * triplet 3.50 π * singlet 3.94 3s triplet 6.83 3s singlet 7.09 7.79 3py triplet 3pz triplet 7.96 7.97 3py singlet 8.12 3pz singlet 8.38 3px singlet Mean absolute error (Rydberg only) a b
SIC + SCF
δ
CASPT2b
δ
3.28 3.42 6.67 6.81 7.47 7.53 7.53 7.61 7.91
− 0.22 − 0.52 − 0.16 − 0.28 − 0.32 − 0.44 − 0.43 − 0.51 − 0.47 0.37
3.48 3.91
− 0.02 − 0.03
7.30
0.21
8.12 8.09 8.32
0.15 − 0.03 − 0.06 0.11
References 46 and 47. Reference 48.
3dπ state, which has the same symmetry, and this greatly increases the excitation energy. This mixing of excited states is not accounted for in the method presented here. A large error is found for this particular state, but the results are quite good for the other states, see Table VI. The Rydberg states of C2 H4 are very well described by the method, with a maximum error of 0.25 eV. Our focus here is on Rydberg states, and for these two molecules the agreement with experiment is again rather good. The errors are slightly larger for H2 CO, where the excitation energy is systematically underestimated, but the spacing between the Rydberg levels is in good agreement with experiment.
FIG. 2. Rydberg orbitals of a H2 O molecule in triplet state. Left: 3pz . Middle: 3py . Right: 3px . The 0.01 Å−3/2 iso-surface is rendered.
V. DISCUSSION
The two step procedure presented here for calculating the excitation energy or Rydberg electron binding energy of molecules using DFT and SIC has been found to give good results for several molecules, with a mean absolute error of 0.18 eV for the 33 Rydberg states calculated. This is the type of error one can expect from a DFT/GGA calculation. The average absolute error in the energy difference between adjacent
FIG. 3. Rydberg orbitals of a trimethylamine molecule, N(CH3 )3 , in triplet state. Left: 3pz . Right: 3pxy . The 0.01 Å−3/2 isosurface is rendered.
TABLE VI. Vertical excitation energy of C2 H4 in eV. Columns labeled with δ give the difference between calculated and experimental values. Results of CASPT2 calculations are also shown for comparison. Excited state
Experimenta
π * triplet 4.36 π * singlet 8.0 3s triplet 6.98 3s singlet 7.11 7.79 3pσ (3py ) triplet 7.80 3pσ (3py ) singlet 7.90 3pσ (3pz ) singlet 8.15 3pπ (3px ) triplet 3pπ (3px ) singlet 8.28 3dσ triplet 8.57 3dσ singlet 8.62 3dπ singlet 9.33 Mean absolute error (Rydberg only) a
Reference 42.
SIC + SCF
δ
CASPT2a
δ
4.43 5.56 7.02 7.12 7.59 7.69 7.73 7.94 8.03 8.40 8.48 9.13
0.07 − 2.44 0.04 0.01 − 0.20 − 0.11 − 0.17 − 0.21 − 0.25 − 0.17 − 0.14 − 0.20 0.15
4.39 8.40 7.05 7.17 7.80 7.85 7.95 8.26 8.40 8.57 8.66 9.31
0.03 0.4 0.07 0.06 0.01 0.05 0.05 0.11 0.12 0.00 0.04 − 0.02 0.05
194102-7
Gudmundsdóttir et al.
Rydberg states is 0.07 eV compared with experimental values. The results are particularly good where the lowest lying excited states have a strong Rydberg character, such as NH3 , H2 O, and TMA. By starting with a triplet state calculation as reference in the first step, one of the electrons is already placed in a Rydberg-like orbital. The results are also good for the two molecules, C2 H4 , and H2 CO, where the low lying excited states have valence-type orbitals, i.e., where the LUMO is a π * orbital. The calculations give, however, a systematic underestimate of the excitation energy, especially for H2 CO. The reason could be an overestimate of the double bond energy in DFT/GGA as discussed below. While the agreement is quite good between the calculated values and experimental values, there are several ways in which further improvements could be made. One is the treatment of the atomic core regions. The frozen core and projectors for the PAW were developed for calculations using traditional density functionals such as PBE. New PAW projectors should be developed in order to describe the inner electrons and the core region when SIC is applied. Then, PAW could be used also in the first step of the two step procedures presented here, the calculation of the Rydberg orbital. Also, the KLI method is an approximation to the OEP and the implementation of full OEP could increase the accuracy. The OEP is important for getting the SIC to apply to the virtual orbitals. The SIC procedure of Perdew and Zunger directly affects the occupied orbitals but only indirectly and insignificantly affects the virtual orbitals. Another item is the SCF calculation which was carried out here using regular DFT without SIC. While the SIC is mainly essential for obtaining a meaningful Rydberg orbital, the final estimate of the total energy in the SCF calculation could also be affected by the self-interaction correction and for consistency the SCF calculation should also be carried out with SIC. This could explain the systematic discrepancy in the excitation energy calculated for H2 CO since it involves a comparison between the energy of an oxygen double bond and a single bond in the excited state. The PBE functional is known to overestimate the bond energy in the O2 molecule by nearly 1 eV.17 Furthermore, consistency between the two steps would result in virtual orbitals in the SCF step that better represent Rydberg orbitals. Our preliminary estimates indicate that the inclusion of SIC in the SCF step would significantly affect the estimated excitation energy of H2 CO. As was shown recently, the atomization energy of molecules is not always improved by adding the selfinteraction correction. In the case of the PBE functional, it has been found that scaling the self-interaction correction by a half gives a significant improvement over normal PBE, while adding the full correction usually gives an overcorrection.17, 19 The full correction is of course necessary for the correct long-range effective potential, so this underscores the importance of developing a new XC functional to go with the selfinteraction free evaluation of the Hartree energy. With an improved self-interaction corrected functional, it will be possible to extend the method presented here to reach self-consistency between the two steps in the calculation, the estimate of the Rydberg orbital and the evaluation of the total energy of the system in the Rydberg state.
J. Chem. Phys. 139, 194102 (2013)
ACKNOWLEDGMENTS
This work was supported by the Icelandic Research Fund (H.G. and H.J.), NER network “Solar Fuel” (H.G. and H.J.), and the Division of Chemical Sciences, Geosciences, and Biosciences, the Office of Basic Energy Sciences, the U.S. Department of Energy, Grant No. DE-FG02-03ER15452 (Y.Z. and P.M.W.). We thank Simon and Peter Klüpfel for assistance and helpful discussions. Helpful discussions with J. P. Perdew are also acknowledged and assistance from S. Kümmel with the implementation of the OEP-KLI. The calculations were carried out in part on “Gardar”, the Nordic High Performance Cluster located in Iceland. 1 R.
S. Freund, in Rydberg States of Atoms and Molecules (Cambridge University, Cambridge, 1983), p. 355. 2 M. P. Minitti, J. D. Cardoza, and P. M. Weber, J. Phys. Chem. A 110, 10212 (2006). 3 N. Kuthirummal and P. M. Weber, Chem. Phys. Lett. 378, 647 (2003). 4 J. L. Gosselin and P. M. Weber, J. Phys. Chem. A 109, 4899 (2005). 5 J. D. Cardoza, F. M. Rudakov, N. Hansen, and P. M. Weber, J. Electron Spectrosc. Relat. Phenom. 165, 5 (2008). 6 F. Rudakov and P. M. Weber, Chem. Phys. Lett. 470, 187 (2009). 7 J. C. Bush, M. P. Minitti, and P. M. Weber, J. Photochem. Photobiol., A 213, 70 (2010). 8 X. Liang, M. G. Levy, S. Deb, J. D. Geiser, R. M. Stratt, and P. M. Weber, J. Mol. Struct. 978, 250 (2010). 9 S. Deb, M. P. Minitti, and P. M. Weber, J. Chem. Phys. 135, 044319 (2011). 10 F. Rudakov and P. M. Weber, J. Chem. Phys. 136, 134303 (2012). 11 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 12 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 13 W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). 14 D. C. Langreth and M. J. Mehl, Phys. Rev. B 28, 1809 (1983). 15 J. P. Perdew, Phys. Rev. Lett. 55, 1665 (1985). 16 J. P. Perdew, Phys. Rev. B 34, 7406 (1986). 17 H. Jónsson, Proc. Natl. Acad. Sci. U.S.A. 108, 944 (2011). 18 S. Klüpfel, P. Klüpfel, and H. Jónsson, Phys. Rev. A 84, 050501 (2011). 19 S. Klüpfel, P. Klüpfel, and H. Jónsson, J. Chem. Phys. 137, 124102 (2012). 20 Y. Zhao and D. Truhlar, J. Phys. Chem. A 110, 13126 (2006). 21 Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai, and K. Hirao, J. Chem. Phys. 120, 8425 (2004). 22 D. J. Tozer and N. C. Handy, J. Chem. Phys. 109, 10180 (1998). 23 M. E. Casida and D. R. Salahub, J. Chem. Phys. 113, 8918 (2000). 24 S. Hirata, C.-G. Zhan, E. Apra, T. L. Windus, and D. A. Dixon, J. Phys. Chem. A 107, 10154 (2003). 25 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 26 J. Garza, J. A. Nichols, and D. A. Dixon, J. Chem. Phys. 112, 7880 (2000). 27 A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, and G. E. Scuseria, J. Chem. Phys. 125, 194112 (2006). 28 J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 46, 5453 (1992). 29 J. Gavnholt, T. Olsen, M. Engelund, and J. Schiøtz, Phys. Rev. B 78, 075441 (2008). 30 M. R. Pederson, R. Heaton, and C. C. Lin, J. Chem. Phys. 80, 1972 (1984); 82, 2688 (1985). 31 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 32 J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B 71, 035109 (2005). 33 J. Enkovaara, C. Rostgaard, J. J. Mortensen et al., J. Phys.: Condens. Matter 22, 253202 (2010). 34 P. Klüpfel, S. Klüpfel, K. Tsemekhman, and H. Jónsson, Lect. Notes Comput. Sci. 7134, 23 (2012). 35 Á. Valdés et al., Phys. Chem. Chem. Phys. 14, 49 (2012). 36 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 37 C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev. B 58, 3641 (1998).
194102-8
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
38 P.
43 W.
39 NIST
44 R.
E. Blöchl, Phys. Rev. B 50, 17953 (1994). Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Release 15b (August 2011), see http://cccbdb.nist.gov/. 40 Calculations on Rydberg excited states of molecular clusters will be presented elsewhere. 41 A. Chutjian, R. I. Hall, and S. Trajmar, J. Chem. Phys. 63, 892 (1975). 42 L. Serrano-Andrés, M. Merchán, I. Nebot-Gil, R. Lindh, and B. O. Roos, J. Chem. Phys. 98, 3151 (1993).
R. Harshbarger, J. Chem. Phys. 54, 2504 (1971). Rianda, R. P. Frueholz, and W. A. Goddard III, Chem. Phys. 19, 131 (1977). 45 M. Rubio, L. Serrano-Andrés, and M. Merchán, J. Chem. Phys. 128, 104305 (2008). 46 D. J. Clouthier and D. A. Ramsay, Annu. Rev. Phys. Chem. 34, 31 (1983). 47 S. Taylor, D. G. Wilden, and J. Comer, Chem. Phys. 70, 291 (1982). 48 M. Merchán and B. O. Roos, Theor. Chim. Acta 92, 227 (1995).
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
Self-interaction corrected density functional calculations of molecular Rydberg states Hildur Gudmundsdóttir,1 Yao Zhang,2 Peter M. Weber,2 and Hannes Jónsson2,3 1
Science Institute of the University of Iceland, 107 Reykjavík, Iceland Department of Chemistry, Brown University, Providence, Rhode Island 02912, USA 3 Faculty of Physical Sciences, VR-III, University of Iceland, 107 Reykjavík, Iceland 2
(Received 11 September 2013; accepted 28 October 2013; published online 15 November 2013) A method is presented for calculating the wave function and energy of Rydberg excited states of molecules. A good estimate of the Rydberg state orbital is obtained using ground state density functional theory including Perdew-Zunger self-interaction correction and an optimized effective potential. The total energy of the excited molecule is obtained using the Delta Self-Consistent Field method where an electron is removed from the highest occupied orbital and placed in the Rydberg orbital. Results are presented for the first few Rydberg states of NH3 , H2 O, H2 CO, C2 H4 , and N(CH3 )3 . The mean absolute error in the energy of the 33 molecular Rydberg states presented here is 0.18 eV. The orbitals are represented on a real space grid, avoiding the dependence on diffuse atomic basis sets. As in standard density functional theory calculations, the computational effort scales as NM2 where N is the number of orbitals and M is the number of grid points included in the calculation. Due to the slow scaling of the computational effort with system size and the high level of parallelism in the real space grid approach, the method presented here makes it possible to estimate Rydberg electron binding energy in large molecules. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829539] I. INTRODUCTION
Molecular Rydberg states, i.e., highly excited electronic states of molecules where one electron is promoted to a hydrogen-like orbital, are profoundly important for theoretical and practical reasons.1 Conceptually, they can be thought of as a positively charged molecular core that is surrounded by an electron in a hydrogen-like electronic state, typically with principal quantum number n ≥ 3. The interactions of the loosely bound electron with the electronic, vibrational, and rotational states of the ion core can be quite complex. But when considering only the electron binding energy, i.e., the energy difference between the molecule in the Rydbergexcited state and the corresponding ion state (with the same quantum numbers in the other degrees of freedom), then many of these complexities cause negligible energy changes. Experimentally, it has been found that even in large molecules the electron binding energy of a Rydberg state is a well-defined quantity that, in spectra of modest resolution, is independent of vibrational and/or rotational excitation.2 Even so, the binding energy has proven to be quite sensitive to the molecular structure.3 This has led to the recent development of spectroscopic tools that use Rydberg states for the identification of molecular structures and possibly of molecules.4–10 For many practical applications, the determination of molecular structures remains an important challenge. Examples include situations when samples are extremely dilute or where rapid throughput is of the essence. There is thus a need for additional structure determination methods such as one relying on Rydberg states. We draw attention to a parallel but much more developed spectroscopic method, nuclear magnetic resonance
0021-9606/2013/139(19)/194102/8/$30.00
(NMR), which also uses the spectroscopic determination of energy differences between molecular quantum states to determine molecular structures. While fundamentally, energy differences are measured, chemists have learned to associate these energy differences with molecular structure, turning NMR spectroscopy into a powerful tool for molecular structure determination. Although there are only relatively few spectral lines, the information content is so decisive that even complicated structures of proteins can now be elucidated. It stands to reason that an analogous approach to molecular structure determination could be based on the measurement of binding energy of electrons in Rydberg states. Experiments have already shown that the electron binding energy is remarkably sensitive to the molecular structure. But to date, the method has been applied in a comparative, rather than an absolute fashion, by identifying, but not determining, molecular structures. Consequently, the technique presently provides a fingerprint of a molecular structure and is referred to as Rydberg Fingerprint Spectroscopy (RFS). The method has proven particularly useful in a photoionization scheme when implemented in a time-resolved way, because the kinetics and dynamics of even complex molecular reactions can be followed in real time. To advance the method from a fingerprint to a structure determination tool, a technique to invert the spectra or otherwise derive the structure from the measured data is required. The aim of the work reported here is to develop a fast method for computing the energy of Rydberg states of large molecules. Rydberg states of molecules have been successfully modeled with high level wave function methods such as CASPT2, but as these methods scale poorly with system size
139, 194102-1
© 2013 AIP Publishing LLC
194102-2
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
they are only suitable for small molecules. Density functional theory (DFT)11, 12 has become a powerful tool for describing the ground state electronic structure of atoms, molecules, and solids13 and with time dependent DFT (TDDFT) it is even possible to calculate the properties of excited electronic states. Rydberg states of molecules have a very large spatial spread compared to the ground state density, which makes their simulation particularly problematic with traditional DFT. The effective potential of an electron, described by a DFT functional of the local density approximation (LDA) or generalized gradient approximation (GGA) type, goes to zero much too quickly due to self-interaction error. While this is frequently not a problem for calculating ground state energy, it means it is impossible to describe Rydberg orbitals with traditional DFT, as the effective potential is far too positive in their region of space. In the present work, we use ground state DFT with the Perdew-Zunger self-interaction correction (SIC) to calculate Rydberg orbitals of the molecules NH3 , H2 O, H2 CO, C2 H4 , and N(CH3 )3 . The SIC gives the correct long range behavior of the potential and as a result the virtual orbitals in the calculations give good estimates of the Rydberg orbitals of the molecules. The Rydberg excitation energy, or the binding energy of Rydberg electrons, is then calculated by placing an electron in the virtual Rydberg orbitals by means of the Delta Self-Consistent Field29 (SCF) method. II. SELF-INTERACTION CORRECTION
In Kohn-Sham (KS) DFT the ground state of a molecule, with an electron density ρ(r), is found by solving the eigenvalue problem 1 − ∇ 2 + vext (r) + vH (r) + vxc (r) ψi = εi ψi , (1) 2 where vext is the external potential of the nuclei, vH is the Hartree potential, accounting for the Coulomb interaction of the charge density with itself, and vxc is the exchangecorrelation (XC) potential. The occupied KS orbitals {ψ i } must span the ground state density: N
|ψi (r)|2 = ρ(r)
and possibly the gradient, of the total density. Explicit density functionals allow fast computation of the ground state energy of molecules to a good approximation. However, the Hartree energy of Eq. (3) not only contains the interaction of each electron with all the others, it also contains the electron’s interaction with itself. This self-Hartree-interaction is only partly canceled in local, explicit density functionals, leaving a spurious self-interaction error. The self-interaction in DFT gives rise to several problems, such as the delocalization of charge, the underestimation of energy barriers and the overestimation of bond energy.17–19 Another manifestation of the electron selfinteraction is the incorrect asymptotic behavior of the XC potential, which makes the description of Rydberg states particularly difficult. At large distances from the atom the correct Kohn-Sham potential has the form 1 lim vxc (r) = − , r
r→∞
(4)
whereas the potential of the explicit, local density functionals goes to zero much faster, as shown in Fig. 1. For example, the LDA potential decays exponentially. The shallow potential of the LDA and GGA functionals leads to an underestimation of the binding energy of electrons, particularly in excited orbitals. The diffuse and extended orbitals characteristic of Rydberg excitations are typically not found to correspond to bound states at this level of approximation. Various approaches have been taken to correct the longrange behavior of the XC potential in DFT in order to better describe Rydberg states. The most commonly used hybrid functionals, like B3LYP, use fractional Hartree-Fock (HF) exchange, so although they may remove some self-interaction problems they do not have the correct asymptotic potential. The Minnesota functional M06-HF is a hybrid functional designed to describe both ground and excited states using full HF exchange and fitting parameters optimized against extensive data.20 Hybrid functionals with the correct long-range potential have also been designed using the long-range correction (LC) scheme, which uses full HF exchange to account for long-range orbital-orbital interactions, dividing the electron
(2)
i
and are introduced as a means to estimate of the kinetic energy. Orbitals are also useful for the approximation of singleparticle properties, and as discussed above, we would like to use the KS orbitals to estimate the binding energy of electrons in a molecule. The usefulness of the orbital concept for describing single particles is very sensitive to the quality of the exchange-correlation potential, vxc . The electron-electron interaction energy is a sum of the Hartree energy ρ(r)ρ(r ) 1 3 (3) EH [ρ] = d rρ(r)vH (r) = d 3 rd 3 r 2 |r − r | and the exchange-correlation energy, Exc , commonly given by an explicit, local density functional such as the LDA12 or GGA.14–16 These are explicit functionals of the total density and depend only on local information about the value,
FIG. 1. Effective potential of an electron in a trimethylamine molecule (TMA) calculated along the direction of one of the N–C bonds. Dashed line: The −1/r potential centered on the central N-atom in the molecule. Solid lines: LDA (blue) and LDA-SIC (red). By applying the self-interaction correction, the correct long range form of the potential is obtained.
194102-3
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
repulsion operator into short-range and long-range parts.21 As both of these approaches require the calculation of exact exchange, they scale less favorably than DFT with system size, but they greatly improve the accuracy of TDDFT calculations of Rydberg states. A simpler, more efficient approach, is to simply give the XC potential the correct behavior in the region important for Rydberg states. The asymptotic correction (AC) to the Kohn-Sham equations gives the required −1/r behavior to the potential beyond a certain distance from the atoms, but retains the conventional DFT form in energetically important regions.22–24 The idea is that the potential is reasonably well represented by conventional DFT functionals in the region of high electron density, close to the atoms, so the AC is only applied to the potential in the outer region. The shiftand-splice AC of Casida and Salahub23 shifts the potential in the high-density region down by a constant, to account for the derivative discontinuity in the exact functional, making the orbital energy of the highest occupied molecular orbital (HOMO) roughly equal to the negative ionization potential of the molecule, and splices the DFT potential to a potential with the correct −1/r asymptotic behavior for large r. The very similar AC of Tozer and Handy22 does not change the form of the potential in the high-density region but shifts the large r potential up to match the DFT potential. The AC can give quite good estimates, for example, Tozer and Handy calculated the excitation energy of CO, N2 , H2 CO, and C2 H4 with mean absolute errors of 0.4, 0.43, 0.24, and 0.05 eV, respectively, for Rydberg excited states.22 So without actually tackling the cause of the problem, the AC estimates the Rydberg excitation energy by covering up the symptoms of an inaccurate functional. In this paper, we investigate whether the removal of orbital self-interaction makes DFT capable of describing Rydberg excited states, using the Perdew-Zunger25 SIC. The derivative discontinuity in the XC potential, as well as the correct long-range form, are obtained when the SIC is applied, making the functional ideal for calculations of Rydberg states. The SIC has previously been used with an optimized effective potential to calculate the energy of the lowest unoccupied molecular orbital (LUMO) for small molecules,26 but to the authors’ knowledge it has not been used to study Rydberg states. An orbital-by-orbital correction is made to the XC energy SIC [{ρi }] = E0 [ρ] − Exc
N
(EH [ρi ] + Exc [ρi ]).
(5)
i=1
We use a simplified notation here neglecting spin. The orbital densities, ρ i = |ψ i |2 , correspond to the occupied orbitals and E0 [ρ] is an explicit, local density functional, such as the LDA or GGA. The correction can in principle be applied to any such exchange-correlation functional, but one should bear in mind that for many-electron systems it is just an approximation. The magnitude of the many-electron self-interaction error is in general not given by the sum of the individual terms of Eq. (5).27 Recent studies of the SIC have found that the accuracy of the correction is dependent on the type of functional it is applied to, that it can greatly improve DFT estimates
of the total energy of atoms,18 and that a scaled down version of it can give good estimates of atomization energies of molecules.19 Furthermore, SIC gives orbital energies that are good estimates of electron binding energy,18 with the HOMO energy giving a good estimate of the ionization potential. The effective potential for the trimethylamine molecule (TMA) is illustrated in Fig. 1, calculated with LDA with and without SIC. It is clear that the effective potential of the LDA goes to zero much too quickly, whereas applying the SIC gives an effective potential with the correct −1/r asymptotic form. The computational cost for the SIC functional scales with system size like traditional DFT, and the functional does not introduce any empirical parameters. The purpose of the present study is to extend and test the accuracy of the SIC for calculations of Rydberg excited states of molecules. We demonstrate that DFT with the PerdewZunger self-interaction correction, applied within the KriegerLi-Iafrate28 (KLI) approximation of the optimized effective potential (OEP), can give good estimates of excited Rydberg orbitals. The total energy of the excited molecule is then obtained in a SCF29 GGA calculation, where the excited electron is placed in a fixed Rydberg state orbital. We present test calculations using this approach and compare with measurements and high level quantum chemistry calculations of the excitation energy or binding energy of the excited electron. The results show that this approach gives quite good accuracy and can be applied to large molecules. The method can, therefore, be helpful in analyzing and interpreting experimental data on Rydberg states of large molecules and can serve as a helpful tool in the analysis of RFS data. III. METHODOLOGY
The electron density can be spanned by or any set of 2 |ϕ (r)| , and thonormal orbitals {ϕj (r)} such that ρ(r) = N j j any of these sets will result in the same total energy in calculations using LDA or GGA functionals. When SIC is applied, however, the energy depends not only on the total electron density, but also on individual orbital densities. This implies that there exists a specific set of orbitals that minimizes the total energy and the energy is no longer invariant under unitary transformations of the orbitals. It is convenient to carry out the calculation using two sets of orbitals: the canonical Kohn-Sham orbitals {ψ i } that are eigenfunctions of the Kohn-Sham Hamiltonian and solutions to Eq. (1), and the energy optimal orbitals {ϕ j }, which minimize the orbital density dependent total energy, Eq. (5).34 The two basis sets are related by a unitary transformation W Wij ψi . (6) ϕj = i
An orbital density dependent (ODD) energy functional has orbital specific potentials vj (r), whose action is only defined on the corresponding occupied orbital, ϕj (r), with orbital density ρj (r) = |ϕj (r)|2 . For the SIC functional, the orbital density dependent part of the Hamiltonian is given by the ODD potential vjSIC (r) = −vH [ρj ] − vxc [ρj ].
(7)
194102-4
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
The calculation can be brought back to the computationally well-developed realm of Kohn-Sham DFT by using an OEP. This gives a local multiplicative potential that captures the effect of the SIC functional on all states, including virtual states. Otherwise, the SIC would only affect the occupied states. Here, we use the KLI28 approximation to the OEP and iteratively find the XC potential given by ⎫ ⎧ N N ⎬ ⎨
1 KLI KLI vxc (r) = ρj (r)vj (r)+ |ψi (r)|2 v¯xc,i − v¯iSIC , ⎭ ρ(r) ⎩ j =1
i=1
(8) KLI v¯xc,i
=
v¯iSIC =
KLI d 3 r|ψi (r)|2 vxc (r),
N
∗ Wij V¯jSIC k Wik ,
(9)
(10)
of Eq. (2) is replaced by the density ρ(r) =
N−1
|ψi (r)|2 + |ψex (r)|2 ,
(13)
i=1
where the Rydberg orbital is approximated by the excited orbital ψ ex ψex (r) =
M i=N
ψi |ψRyd
i
| ψi |ψRyd |2
1/2 ψi (r),
(14)
which is a linear combination of empty KS orbitals in a calculation with M orbitals. The SCF calculation gives the total energy of the excited state of the system and the excitation energy can then be calculated by subtracting the ground state energy. Sometimes one is interested in the binding energy of the excited electron, since this is what is measured in RFS, and in this case the ground state energy of the cation is subtracted, rather than the energy of the neutral molecule.
j,k=1
where the Hermitian SIC potential matrix is V¯jSIC k = 12 ϕj | vjSIC + vkSIC |ϕk . Using the KLI potential, the orbitals can be obtained from Hˆ SIC ψi = εi ψi ,
(11)
KLI = H + vxc (r) is the where the SIC Hamiltonian H 0 sum of the unitary invariant part Hˆ = − 12 ∇ 2 + vext (r) 0 0 (r), where vxc is the XC potential being cor+ vH (r) + vxc KLI , rected, and the orbital density dependent KLI potential, vxc
ˆ SIC
ˆ0
which contains the effect of the SIC. To minimize the energy with respect to the optimal orbitals, one must simultaneously solve Eq. (11) and find W. This is done by iteratively solving Eq. (11), as in a regular DFT calculation, and in every step of the self-consistent calculation the transformation W is estimated. For a given set of KS orbitals {ψ i }, the estimated W will give the energyoptimal orbitals {ϕ j } and the ODD potentials vi through Eqs. (6) and (7). The error in the energy optimal orbitals, and therefore in W, is given by the antihermitian SIC potential matrix, κ 30 κj k = ϕj | vjSIC − vkSIC |ϕk ,
(12)
which is minimized with respect to the unitary transformation W in each step of the KS iterative process. When κ ≈ 0 the energy optimal orbitals {ϕ j } and the ODD potentials vjSIC , have been found to the chosen level of accuracy, and the KLI potential (10) can be calculated and used in the Kohn-Sham Eq. (11). When self-consistency is reached κ = 0 and the Kohn-Sham equations have been solved. The eigenvalues εi turn out to give a rough approximation of the electron binding energy, and the virtual KS orbitals, {ψ i } where i > N, can be used to estimate the excited states of the molecule. DFT is a ground state theory and to calculate the energy of the Rydberg excited state, the linear-expansion SCF method29 is used with the Perdew-Burke-Ernzerhof (PBE) functional.31 Having calculated the virtual orbitals of the molecule with the SIC functional, a Rydberg-type orbital, ψRyd (r), is selected and an electron is removed from the HOMO and placed in this orbital. The ground state density
A. Implementation
The calculations were carried out with the GPAW software,32, 33 extended to include self-consistent SIC calculations.34, 35 The Rydberg orbitals were approximated using the LDA functional36 and SIC. Here, one can choose the spin multiplicity and in several cases it has been found to give better results to calculate the orbitals for the triplet state of the molecule, rather than the singlet, so the calculations presented here are performed for triplet state molecules. This is a particularly good approach when a molecule has no low lying valence excitations, and the LUMO is a Rydberg orbital, as for H2 O, NH3 , and TMA. A uniform grid with a mesh size of ∼0.15 Å was used to represent the valence electrons, with a cubic simulation cell of side length 20 Å for all molecules except TMA, where the side length was 25 Å. The core electrons were kept frozen, and to represent the core region of each atom two methods were used. In the SIC calculations, Hartwigsen-GoedeckerHutter (HGH) pseudopotentials37 were used. The SCF calculations were performed with the projector augmented wave (PAW)38 method. The PAW was not used in the first step because PAW projectors have not yet been developed for the SIC functional. All calculations were performed using the experimental geometry of the molecules.39 One advantage of using a real space grid is avoiding the need to develop special basis functions to represent the Rydberg orbital. If a Gaussian basis set were used, highly diffuse Gaussians would need to be included and a choice made where to center the diffuse basis functions. In the real space grid, no special treatment is needed for the Rydberg orbital since each point in space is equally well represented. To accurately represent the Rydberg state as a linear combination of empty KS states in the SCF calculation, as in Eq. (14), around 40 orbitals were used. Since the SIC does not increase the scaling of the computational effort with the size of the system, the method can be applied to large molecules and even clusters of molecules.40 The computational effort scales as DFT calculations with the parent functional, as NM2 for grid-based implementations of
194102-5
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
TABLE I. Comparison of eigenvalues and calculations of binding energy. The table shows eigenvalues, ε, and negative binding energy (in eV) of the first two Rydberg states of a water and an ethylene molecule in the singlet state. The eigenvalues are calculated with the LDA, PBE, and LDA-SIC functionals. The SIC + SCF method is used to calculate the energy of the excited state, which is given relative to the energy of the cation to estimate the binding energy. The experimental reference is the difference between the measured excitation energy41, 42 and ionization energy (12.62 eV for H2 O and 10.68 eV for C2 H4 ).39 ε LDA-SIC SIC + SCF Experiment
State
ε LDA
ε PBE
H2 O
3s 3py
− 0.93 0.31
− 0.94 0.28
− 6.29 − 4.13
− 5.42 − 3.81
− 5.22 − 3.52
C2 H4
3s 3py
− 0.30 0.17
− 0.33 0.13
− 4.09 − 3.32
− 3.51 − 2.94
− 3.57 − 2.88
TABLE III. Vertical excitation energy for H2 O in eV. The molecule is in the yz plane, with the electric dipole in the z direction. Columns labeled with δ give the difference between calculated and experimental values. Results of CASPT2 calculations are also shown for comparison. Excited state 3s triplet 3s singlet 3py triplet 3py singlet 3px triplet 3pz triplet 3pz singlet 3px singlet Mean absolute error a b
GGA, where N is the number of orbitals and M is the number of grid points. Similar scaling with system size would be obtained with basis set approaches. For example, the computational time for the SIC calculation on 16 Opteron cores (two quad core nodes, 2.8 GHz clock speed) was 6000 s for NH3 (the box was 20 Å long on each side, a total of 1323 grid points included, and 36 orbitals), while the calculation took 34 000 s for TMA (the box was 25 Å long on each side, a total of 1683 grid points included and 54 orbitals). The scaling of the computer time is 5.7, close to the ratio of NM2 which is 54 × 1686 /36 × 1326 = 6.4 for these two calculations. A comparison of various levels of theory for estimating the binding energy of the Rydberg electron is given in Table I. As mentioned earlier, explicit local DFT functionals will give a poor approximation of the binding energy, and Table I shows the improvement with the inclusion of SIC. An even better estimate is provided with the SCF calculations, based on Rydberg orbitals calculated with SIC. IV. RESULTS
Tables II–VI show the calculated excitation energy compared with experimental results. In the case of TMA the binding energy of the excited electron is shown, as this has been measured experimentally. Where the experimental data are available, comparison is made to both singlet and triplet state energy. For all molecules, except TMA, the results are also TABLE II. Vertical excitation energy for NH3 in eV. The rotational axis of the molecule is in the z direction. Columns labeled with δ give the difference between calculated and experimental values. Results of configuration interaction (CI) calculations are also shown for comparison. Excited state
Experimenta SIC + SCF
δ
CIb
a b
Reference 43. Reference 44.
6.39 7.93 8.26
6.15 6.39 7.62 7.98 8.02 8.25
6.14 0.00 6.37 7.86 0.05 7.88 7.87 − 0.01 8.15 0.02
7.0 7.4 8.9 9.1 9.81 9.98 10.01 10.16
7.11 7.38 8.77 8.94 9.62 9.79 9.95 9.83
− 0.05 − 0.11 0.06
CASPT2b
δ
0.11 − 0.02 − 0.13 − 0.16 − 0.19 − 0.19 − 0.06 − 0.33 0.15
7.13 7.50 9.12 9.27 9.87 9.90 9.95 10.15
0.13 0.10 0.22 0.17 0.06 − 0.08 − 0.06 − 0.01 0.10
compared with previous quantum chemical calculations. For NH3 results are compared with configuration interaction (CI) calculations, and for H2 O, H2 CO, and C2 H4 comparison is made with CASPT2 calculations. The calculated excitation energy of the 3s and 3p Rydberg states of NH3 and H2 O is shown in Tables II and III. The lowest lying vertical excitations of these molecules are Rydberg states. The agreement here is good and indicates that the approach presented can be used for interpretation and analysis of Rydberg states. Figure 2 shows the three 3p Rydberg orbitals for H2 O. These are simply unoccupied eigenstates of the SIC Hamiltonian. The 3px and 3py have a regular p-orbital shape, while the 3pz which lies in the direction of the dipole moment is quite a bit distorted. The method presented was also applied to TMA, a larger molecule which has been studied recently by photoionization spectroscopy.5 The calculated and measured values of the Rydberg electron binding energy are compared in Table IV and some of the calculated Rydberg state orbitals are shown in Fig. 3. The method was also applied to two molecules with π * valence excitations lying below the Rydberg states: H2 CO and C2 H4 . The method presented here is based on ground state DFT and is not expected to be able to account for the complexities of valence excitations. For example in C2 H4 , the π * singlet valence excitation is strongly mixed with the Rydberg TABLE IV. Binding energy of the Rydberg state electron for trimethylamine (TMA). The z axis runs through the rotational axis of molecule. The column labeled with δ gives the difference between calculated and experimental values.
δ
− 0.02
δ
Reference 41. Reference 45.
Excited state 3s triplet 3s singlet 3pxy triplet 3pxy singlet 3pz triplet 3pz singlet Mean absolute error
Experimenta SIC + SCF
3s triplet 3s singlet 3pxy triplet 3pxy singlet 3pz triplet 3pz singlet Mean absolute error a
Reference 5.
Experimenta
3.09 2.25 2.20
SIC + SCF 3.35 3.06 2.59 2.40 2.61 2.29
δ
− 0.03 0.15 0.09 0.10
194102-6
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
TABLE V. Vertical excitation energy of H2 CO in eV. The molecule is in the yz plane, with the CO bond in the z direction. Columns labeled with δ give the difference between calculated and experimental values. While there is a systematic over binding in the SIC calculation, the relative energy of the excited states is quite well reproduced. Results of CASPT2 calculations are also shown for comparison. Excited state
Experimenta
π * triplet 3.50 π * singlet 3.94 3s triplet 6.83 3s singlet 7.09 7.79 3py triplet 3pz triplet 7.96 7.97 3py singlet 8.12 3pz singlet 8.38 3px singlet Mean absolute error (Rydberg only) a b
SIC + SCF
δ
CASPT2b
δ
3.28 3.42 6.67 6.81 7.47 7.53 7.53 7.61 7.91
− 0.22 − 0.52 − 0.16 − 0.28 − 0.32 − 0.44 − 0.43 − 0.51 − 0.47 0.37
3.48 3.91
− 0.02 − 0.03
7.30
0.21
8.12 8.09 8.32
0.15 − 0.03 − 0.06 0.11
References 46 and 47. Reference 48.
3dπ state, which has the same symmetry, and this greatly increases the excitation energy. This mixing of excited states is not accounted for in the method presented here. A large error is found for this particular state, but the results are quite good for the other states, see Table VI. The Rydberg states of C2 H4 are very well described by the method, with a maximum error of 0.25 eV. Our focus here is on Rydberg states, and for these two molecules the agreement with experiment is again rather good. The errors are slightly larger for H2 CO, where the excitation energy is systematically underestimated, but the spacing between the Rydberg levels is in good agreement with experiment.
FIG. 2. Rydberg orbitals of a H2 O molecule in triplet state. Left: 3pz . Middle: 3py . Right: 3px . The 0.01 Å−3/2 iso-surface is rendered.
V. DISCUSSION
The two step procedure presented here for calculating the excitation energy or Rydberg electron binding energy of molecules using DFT and SIC has been found to give good results for several molecules, with a mean absolute error of 0.18 eV for the 33 Rydberg states calculated. This is the type of error one can expect from a DFT/GGA calculation. The average absolute error in the energy difference between adjacent
FIG. 3. Rydberg orbitals of a trimethylamine molecule, N(CH3 )3 , in triplet state. Left: 3pz . Right: 3pxy . The 0.01 Å−3/2 isosurface is rendered.
TABLE VI. Vertical excitation energy of C2 H4 in eV. Columns labeled with δ give the difference between calculated and experimental values. Results of CASPT2 calculations are also shown for comparison. Excited state
Experimenta
π * triplet 4.36 π * singlet 8.0 3s triplet 6.98 3s singlet 7.11 7.79 3pσ (3py ) triplet 7.80 3pσ (3py ) singlet 7.90 3pσ (3pz ) singlet 8.15 3pπ (3px ) triplet 3pπ (3px ) singlet 8.28 3dσ triplet 8.57 3dσ singlet 8.62 3dπ singlet 9.33 Mean absolute error (Rydberg only) a
Reference 42.
SIC + SCF
δ
CASPT2a
δ
4.43 5.56 7.02 7.12 7.59 7.69 7.73 7.94 8.03 8.40 8.48 9.13
0.07 − 2.44 0.04 0.01 − 0.20 − 0.11 − 0.17 − 0.21 − 0.25 − 0.17 − 0.14 − 0.20 0.15
4.39 8.40 7.05 7.17 7.80 7.85 7.95 8.26 8.40 8.57 8.66 9.31
0.03 0.4 0.07 0.06 0.01 0.05 0.05 0.11 0.12 0.00 0.04 − 0.02 0.05
194102-7
Gudmundsdóttir et al.
Rydberg states is 0.07 eV compared with experimental values. The results are particularly good where the lowest lying excited states have a strong Rydberg character, such as NH3 , H2 O, and TMA. By starting with a triplet state calculation as reference in the first step, one of the electrons is already placed in a Rydberg-like orbital. The results are also good for the two molecules, C2 H4 , and H2 CO, where the low lying excited states have valence-type orbitals, i.e., where the LUMO is a π * orbital. The calculations give, however, a systematic underestimate of the excitation energy, especially for H2 CO. The reason could be an overestimate of the double bond energy in DFT/GGA as discussed below. While the agreement is quite good between the calculated values and experimental values, there are several ways in which further improvements could be made. One is the treatment of the atomic core regions. The frozen core and projectors for the PAW were developed for calculations using traditional density functionals such as PBE. New PAW projectors should be developed in order to describe the inner electrons and the core region when SIC is applied. Then, PAW could be used also in the first step of the two step procedures presented here, the calculation of the Rydberg orbital. Also, the KLI method is an approximation to the OEP and the implementation of full OEP could increase the accuracy. The OEP is important for getting the SIC to apply to the virtual orbitals. The SIC procedure of Perdew and Zunger directly affects the occupied orbitals but only indirectly and insignificantly affects the virtual orbitals. Another item is the SCF calculation which was carried out here using regular DFT without SIC. While the SIC is mainly essential for obtaining a meaningful Rydberg orbital, the final estimate of the total energy in the SCF calculation could also be affected by the self-interaction correction and for consistency the SCF calculation should also be carried out with SIC. This could explain the systematic discrepancy in the excitation energy calculated for H2 CO since it involves a comparison between the energy of an oxygen double bond and a single bond in the excited state. The PBE functional is known to overestimate the bond energy in the O2 molecule by nearly 1 eV.17 Furthermore, consistency between the two steps would result in virtual orbitals in the SCF step that better represent Rydberg orbitals. Our preliminary estimates indicate that the inclusion of SIC in the SCF step would significantly affect the estimated excitation energy of H2 CO. As was shown recently, the atomization energy of molecules is not always improved by adding the selfinteraction correction. In the case of the PBE functional, it has been found that scaling the self-interaction correction by a half gives a significant improvement over normal PBE, while adding the full correction usually gives an overcorrection.17, 19 The full correction is of course necessary for the correct long-range effective potential, so this underscores the importance of developing a new XC functional to go with the selfinteraction free evaluation of the Hartree energy. With an improved self-interaction corrected functional, it will be possible to extend the method presented here to reach self-consistency between the two steps in the calculation, the estimate of the Rydberg orbital and the evaluation of the total energy of the system in the Rydberg state.
J. Chem. Phys. 139, 194102 (2013)
ACKNOWLEDGMENTS
This work was supported by the Icelandic Research Fund (H.G. and H.J.), NER network “Solar Fuel” (H.G. and H.J.), and the Division of Chemical Sciences, Geosciences, and Biosciences, the Office of Basic Energy Sciences, the U.S. Department of Energy, Grant No. DE-FG02-03ER15452 (Y.Z. and P.M.W.). We thank Simon and Peter Klüpfel for assistance and helpful discussions. Helpful discussions with J. P. Perdew are also acknowledged and assistance from S. Kümmel with the implementation of the OEP-KLI. The calculations were carried out in part on “Gardar”, the Nordic High Performance Cluster located in Iceland. 1 R.
S. Freund, in Rydberg States of Atoms and Molecules (Cambridge University, Cambridge, 1983), p. 355. 2 M. P. Minitti, J. D. Cardoza, and P. M. Weber, J. Phys. Chem. A 110, 10212 (2006). 3 N. Kuthirummal and P. M. Weber, Chem. Phys. Lett. 378, 647 (2003). 4 J. L. Gosselin and P. M. Weber, J. Phys. Chem. A 109, 4899 (2005). 5 J. D. Cardoza, F. M. Rudakov, N. Hansen, and P. M. Weber, J. Electron Spectrosc. Relat. Phenom. 165, 5 (2008). 6 F. Rudakov and P. M. Weber, Chem. Phys. Lett. 470, 187 (2009). 7 J. C. Bush, M. P. Minitti, and P. M. Weber, J. Photochem. Photobiol., A 213, 70 (2010). 8 X. Liang, M. G. Levy, S. Deb, J. D. Geiser, R. M. Stratt, and P. M. Weber, J. Mol. Struct. 978, 250 (2010). 9 S. Deb, M. P. Minitti, and P. M. Weber, J. Chem. Phys. 135, 044319 (2011). 10 F. Rudakov and P. M. Weber, J. Chem. Phys. 136, 134303 (2012). 11 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 12 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 13 W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). 14 D. C. Langreth and M. J. Mehl, Phys. Rev. B 28, 1809 (1983). 15 J. P. Perdew, Phys. Rev. Lett. 55, 1665 (1985). 16 J. P. Perdew, Phys. Rev. B 34, 7406 (1986). 17 H. Jónsson, Proc. Natl. Acad. Sci. U.S.A. 108, 944 (2011). 18 S. Klüpfel, P. Klüpfel, and H. Jónsson, Phys. Rev. A 84, 050501 (2011). 19 S. Klüpfel, P. Klüpfel, and H. Jónsson, J. Chem. Phys. 137, 124102 (2012). 20 Y. Zhao and D. Truhlar, J. Phys. Chem. A 110, 13126 (2006). 21 Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai, and K. Hirao, J. Chem. Phys. 120, 8425 (2004). 22 D. J. Tozer and N. C. Handy, J. Chem. Phys. 109, 10180 (1998). 23 M. E. Casida and D. R. Salahub, J. Chem. Phys. 113, 8918 (2000). 24 S. Hirata, C.-G. Zhan, E. Apra, T. L. Windus, and D. A. Dixon, J. Phys. Chem. A 107, 10154 (2003). 25 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 26 J. Garza, J. A. Nichols, and D. A. Dixon, J. Chem. Phys. 112, 7880 (2000). 27 A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, and G. E. Scuseria, J. Chem. Phys. 125, 194112 (2006). 28 J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 46, 5453 (1992). 29 J. Gavnholt, T. Olsen, M. Engelund, and J. Schiøtz, Phys. Rev. B 78, 075441 (2008). 30 M. R. Pederson, R. Heaton, and C. C. Lin, J. Chem. Phys. 80, 1972 (1984); 82, 2688 (1985). 31 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 32 J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B 71, 035109 (2005). 33 J. Enkovaara, C. Rostgaard, J. J. Mortensen et al., J. Phys.: Condens. Matter 22, 253202 (2010). 34 P. Klüpfel, S. Klüpfel, K. Tsemekhman, and H. Jónsson, Lect. Notes Comput. Sci. 7134, 23 (2012). 35 Á. Valdés et al., Phys. Chem. Chem. Phys. 14, 49 (2012). 36 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 37 C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev. B 58, 3641 (1998).
194102-8
Gudmundsdóttir et al.
J. Chem. Phys. 139, 194102 (2013)
38 P.
43 W.
39 NIST
44 R.
E. Blöchl, Phys. Rev. B 50, 17953 (1994). Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Release 15b (August 2011), see http://cccbdb.nist.gov/. 40 Calculations on Rydberg excited states of molecular clusters will be presented elsewhere. 41 A. Chutjian, R. I. Hall, and S. Trajmar, J. Chem. Phys. 63, 892 (1975). 42 L. Serrano-Andrés, M. Merchán, I. Nebot-Gil, R. Lindh, and B. O. Roos, J. Chem. Phys. 98, 3151 (1993).
R. Harshbarger, J. Chem. Phys. 54, 2504 (1971). Rianda, R. P. Frueholz, and W. A. Goddard III, Chem. Phys. 19, 131 (1977). 45 M. Rubio, L. Serrano-Andrés, and M. Merchán, J. Chem. Phys. 128, 104305 (2008). 46 D. J. Clouthier and D. A. Ramsay, Annu. Rev. Phys. Chem. 34, 31 (1983). 47 S. Taylor, D. G. Wilden, and J. Comer, Chem. Phys. 70, 291 (1982). 48 M. Merchán and B. O. Roos, Theor. Chim. Acta 92, 227 (1995).
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
