Valence Virtual Orbitals: An Unambiguous ab Initio ... - ACS Publications

Article pubs.acs.org/JPCA

Valence Virtual Orbitals: An Unambiguous ab Initio Quantification of the LUMO Concept Michael W. Schmidt,* Emily A. Hull, and Theresa L. Windus Department of Chemistry, Iowa State University, Ames, Iowa 50011, United States ABSTRACT: Many chemical concepts hinge on the notion of an orbital called the lowest unoccupied molecular orbital, or LUMO. This hypothetical orbital and the much more concrete highest occupied molecular orbital (HOMO) constitute the two “frontier orbitals”, which rationalize a great deal of chemistry. A viable LUMO candidate should have a sensible energy value, a realistic shape with amplitude on those atoms where electron attachment or reduction or excitation processes occur, and often an antibonding correspondence to one of the highest occupied MOs. Unfortunately, today’s quantum chemistry calculations do not yield useful empty orbitals. Instead, the empty canonical orbitals form a large sea of orbitals, where the interesting valence antibonds are scrambled with the basis set’s polarization and diffuse augmentations. The LUMO is thus lost within a continuum associated with a detached electron, as well as many Rydberg excited states. A suitable alternative to the canonical orbitals is proposed, namely, the valence virtual orbitals. VVOs are found by a simple algorithm based on singular value decomposition, which allows for the extraction of all valence-like orbitals from the large empty canonical orbital space. VVOs are found to be nearly independent of the working basis set. The utility of VVOs is demonstrated for construction of qualitative MO diagrams, for prediction of valence excited states, and as starting orbitals for more sophisticated calculations. This suggests that VVOs are a suitable realization of the LUMO, LUMO + 1, ... concept. VVO generation requires no expert knowledge, as the number of VVOs sought is found by counting s-block atoms as having only a valence s orbital, transition metals as having valence s and d, and main group atoms as being valence s and p elements. Closed shell, open shell, or multireference wave functions and elements up to xenon may be used in the present program.

I. INTRODUCTION I.A. Motivation and Outline. All chemists know, deep in their souls, that molecules and extended systems are composed of atoms. Furthermore, computational methods for molecules usually involve building up molecular wave functions using atom-centered basis functions (plane-wave calculations famously have difficulty building up the density peaks at atomic nuclei). Consequently, many concepts in chemistry are based upon single-particle molecular orbital (MO) theory, using linear combinations of atomic orbitals (LCAO) to expand the MOs. The valence electrons of stable molecules usually occupy bonding linear combinations of atomic orbitals (AO) and nonbonded lone pairs, often delocalized across broad regions of the molecule, while mainly avoiding the use of antibonding AO combinations. The exact electron count determines the occupancy. For example, O2 occupies all available bonding orbitals and then half-fills a π* antibond. Electron poor systems might not utilize all available bonding levels. Of particular conceptual importance are the two frontier orbitals,1,2 known as the highest occupied and lowest unoccupied MOs (HOMO and LUMO), and their energy difference, known as the “band gap”. The frontier orbitals are well-known to play a key role in explaining chemical reactivity patterns. Quantum chemistry calculations rather naturally produce the HOMO along with all lower energy occupied orbitals (HOMO − 1, ...), as the occupied canonical orbitals. The term “canonical” means those orbitals that diagonalize the Fock operator of self-consistent field (SCF) theory. However, the © XXXX American Chemical Society

lowest empty canonical orbitals often have little valence antibonding character, as will be shown below. This is well recognized among theoretical chemists, as may be seen from the number of papers seeking alternative empty orbitals (reviewed just below). The poor connection between the LUMO concept and the lowest empty canonical orbital seems to be less well understood in the wider chemical community. The purpose of this paper is to demonstrate that useful, chemically relevant, and interpretable unoccupied orbitals can be obtained from MO calculations in a simple, robust, and automated way,3 from any type of SCF or MCSCF calculation. These alternative empty MOs are called the valence virtual orbitals (VVOs).4 Accordingly, the main body of this paper deliberately avoids use of the familiar acronym LUMO. Instead, for precision, the first few empty canonical MOs are designated LCMO, LCMO + 1, LCMO + 2, ... to avoid confusion with members of the set of all valence virtual orbitals termed LVVO, LVVO + 1, etc. These two sets are compared and evaluated for their suitability as the conceptual LUMO, LUMO + 1, ... orbitals. The outline of the paper is as follows. The second half of the introduction is a review of various procedures proposed by theoretical chemists to improve upon the empty CMOs. Section II presents the equations for VVO generation in a simpler form than a decade ago,4 since the present paper is not Received: July 16, 2015 Revised: September 21, 2015

A

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

appropriate for spectroscopy or for electron correlation recovery.16 Writing the usual closed shell Fock operator with six user selectable coefficients in front of its usual kinetic energy, nuclear−electron attraction, and Coulomb and exchange operators,

concerned with the closely related atom localized orbitals known as QUAMBOs4 and QUAOs.5,6 Section III describes the characteristics and applications of VVOs. The VVOs are shown to be nearly independent of the quantitative basis set in which a calculation is performed, in their shape and in their energy. The VVO energies allow the construction of a full MO diagram including all antibonding MOs. Such MO diagrams are found to be qualitatively similar for Hartree−Fock and DFT calculations but not quantitatively identical in energy. The LVVO is found to predict the correct shape for an excited electron in the first excited states of molecules with low-lying states (for example, dyes). Frequent comparisons of VVOs to the usual empty CMOs are made along the way. VVOs are also shown to be excellent starting orbitals for more sophisticated calculations that begin to occupy antibonding orbitals, such as multiconfigurational SCF (MCSCF). The paper ends with a short summary in section IV. The valence virtual orbitals address only those aspects of chemistry that involve the valence orbitals of the atoms comprising the molecule. Quite often molecular excited states have Rydberg character, involving spatially large, nonvalence shapes. On rare occasions, such as dipole bound anions7 or solvated electrons,8 even the ground state may be nonvalence in character: all such situations are beyond the scope of this paper. I.B. Review of Schemes To Improve upon the CMOs. There is a long history of seeking virtual orbitals different from the canonical molecular orbitals that are the usual direct output of molecular orbital calculations. These are grouped into four main classes below: diagonalization of an energy operator, diagonalization of a molecular density matrix, direct imposition of atomic valence character (which includes the present work), and atomic sub-block density diagonalization. Often the goal has been to obtain better convergence of configuration interaction (CI) energies9 (e.g., Boys’ 1960 “oscillator orbitals”, based on dipole moments10). So the first two categories do not necessarily separate the virtual space into two separate subspaces, namely, valence plus nonvalence. They are nonetheless historically important and are briefly reviewed, before the two groups of methods intended to cleanly separate valence virtuals from all other virtuals. The nonvalence empty orbitals are called the external orbitals in this paper. Diagonalization of Fock-like or exchange-like operators within the canonical virtual orbital space of a converged SCF calculation leads to altered virtuals, without any change to the molecule’s occupied orbitals.11,12 The improved virtual orbitals (IVOs) are obtained by diagonalizing a Fock operator corresponding to the removal of one electron from the highest occupied orbital.13 Modified virtual orbitals (MVOs) are obtained from a Fock matrix after removal of all valence electrons, to further increase the “Coulombic attraction” of the virtual space into the “valence region” of the molecule.14 Averaged virtual orbitals (AVOs) remove one electron, but spread this hole in the electron density out evenly over all occupied valence orbitals.15 Two schemes suggested by these procedures are available in the GAMESS program, namely, removal of the user’s choice of electrons from several of the highest occupied orbitals (often 6−10, so more than 1 as for IVOs but usually not all valence electrons as for MVOs), or else to remove half of the valence electrons (one electron from every valence orbital, not just one electron as for AVOs). It should also be noted that alternative Fock operators corresponding to Hamiltonians with an unchanged number of electrons can be formulated to give virtual orbitals more

F = aT + bVNE + cJ CORE + dK CORE + eJ VALENCE + fKVALENCE

facilitates discussion of some additional energy operator schemes. Of course, the canonical virtual orbitals are obtained from the usual Fock operator (namely, a = b = 1, c = e = 2, d = f = −1 for closed shell systems). Rather than remove electron density when creating Coulomb and exchange operators, as in the various aforementioned “hole operator” schemes, Cooper and Pounder17 suggested increasing the strength of the nuclear attraction, namely, to increase the coefficient b. A number of workers18−20 have used the exchange operator because the Hamiltonian matrix element for the double excitation φi2 → φa2 is given by the exchange integral Kia (i is an occupied MO, a is an empty MO). Maximization of exchange interactions should improve singles and doubles configuration interaction (CI-SD) energies and is accomplished by diagonalization of the valence exchange operator (a = b = c = d = e = 0 with f = −1). It is most practical to use the total valence exchange operator instead of working with one occupied orbital at a time. Finally, the Korbitals proposed by Feller and Davidson21 use the closely related operator 0.06F − KVALENCE, namely, a = b = 0.06, c = e = 0.12, d = −0.06, and f = −1.06, and are one of the best virtual orbital choices for producing compact CI expansions. Anyone wishing to experiment with energy operators constructed with the six coefficients above can do so with the GAMESS program. The computational cost of all these procedures is very low, typically just a single SCF iteration’s time to build the desired energy operator, plus one diagonalization in the virtual space. Only one of several studies of the rate of CI convergence using different virtual orbitals is mentioned here,22 to illustrate that all schemes in this “energy operator” category are expected to be more effective in truncated CI calculations than use of the canonical virtual orbitals. A second major category of virtual orbital adjustment involves diagonalization of some molecular density matrix. Requiring the density diagonalization to occur only within the SCF calculation’s virtual block leaves the occupied SCF orbitals unchanged. The resulting weakly occupied natural orbitals with the larger occupation numbers will contain the important antibonding valence orbitals but will also include in−out or angular correlating shapes as well. Since the idea is to have a reasonable but economical procedure, the density matrix is usually obtained at a doubles level of theory. To name just a few, the density might be obtained using second order perturbation theory (MP2)23,24 or perhaps CI-SD.25−27 It is natural to consider more rapidly obtainable density matrices as well. In the limit of a small molecule, with at least as many valence electrons as there are valence orbitals, one might simply do a high-spin-multiplicity open shell SCF calculation, obtaining all valence orbitals, since all are occupied (e.g., H2O in a quintet state, to occupy both OH σ* orbitals). In larger molecules, one may average the density matrices obtained from several SCF calculations (ground state and a few excited states) as a means to occupy more of the valence orbitals.28 A systematic procedure for obtaining those few empty valence orbitals involved in “strong correlation” is to B

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

identifying (somewhat distorted) atom centered entities within the molecular wave function, the AIM theory is not orbital based but rather partitions the molecule’s electron density into atomic regions. In the present work, the external virtual orbital space is of little interest, after it has been successfully separated from the valence virtual space. However, selection of atomically or at least regionally localized orbitals within the external space facilitates the truncation of effort needed for recovery of dynamical electron correlation energy. Consequently the literature on this topic is now quite large, from which a very few references are given here for the interested reader. A particularly relevant paper is that of Subotnik et al., whose procedure40 first projects out the valence virtual space by a process that falls in the third category mentioned here and then localizes the remaining “hard virtuals”. Very often, orthogonality is sacrificed in order to obtain a much higher degree of localization, as pioneered in Pulay’s projected atomic orbitals.52 Two very recent attempts to localize the external space with retention of orbital orthogonality use an SVD procedure53 or a fourth power variant54,55 of the Pipek−Mezey localization. In closing, note that several of the procedures just reviewed, including ref 3, have been adapted to the study of periodic systems,56−58 even for plane wave calculations.

extract natural orbitals from the total density matrix of spinunrestricted open shell wave functions29,30 (known as the UNO-CAS method). The two diagonalization schemes (energy-type or densitytype) just discussed often produce within their “lower energy” or “greater occupancy” orbitals some additional orbitals that do not strictly speaking have valence character. Many of the most important empty orbitals will in fact be valence orbitals, such as σ* or π*, but other orbitals with just one new node (such as in−out correlating orbitals for lone pairs) will also occur. While important for improved CI convergence, perhaps, such orbitals with radial or angular correlating nodes are not necessarily very interesting for valence chemistries. Therefore, it is desirable to seek a procedure that (i) cleanly separates the canonical virtual space into valence and external subspaces, (ii) finds all valence virtual orbitals, and (iii) does so without requiring complicated algorithms or long computer times. The third group of methods targets the identification of only valence virtual orbitals. It is apparent that the objective can be achieved by projecting some prototype collection of the atomic valence orbitals onto the molecular virtual space. Projections were first used quite some time ago,31,32 albeit only after explorations of the first two diagonalization schemes were tried. Justification for the idea that valence atomic orbitals are the predominant contributors to valence molecular orbitals comes from projecting free atom SCF orbitals onto full valence type MCSCF wave functions: overlaps between the free atom AOs and their projections onto the molecular orbitals remain at or above 0.99.33 The present scheme for the projection of valence orbitals out of the virtual space of closed shell functions to yield atom localized orbitals was presented in 2004.3 This algorithm is recast in a simpler form in the present paper, as an instance of the singular value decomposition (SVD). The SVD procedure (also known in quantum chemistry as the method of corresponding orbitals) consists of basis rotations in two different orbital spaces to bring the bases into maximum coincidence with each other; namely, it is a mutual projection of each space onto the other.34−38 Many others have used similar projections of AOs onto molecular orbitals to produce the extracted polarized atomic orbitals,39,40 enveloping localized orbitals,41 intrinsic minimal atomic basis,42 molecule-adapted atomic orbitals,43 or intrinsic atomic orbitals.44 A great many interesting chemical applications of the intrinsic atomic orbitals were presented in the latter paper.44 The equivalence of the orbitals in ref 44 to those in refs 3 and 5 has been shown.45 Recent work at this university has generalized the present work’s methodology3 to all types of SCF wave functions, including multiconfigurational MCSCF,5,6 and demonstrated useful applications6 to chemical bonding and charge population analyses. Perhaps the most commonly applied procedure for generating chemically meaningful virtual orbitals is the natural bond order analysis.46−48 Its key step is diagonalization of the density matrix in atomic sub-blocks, producing atom-localized hybrid valence orbitals. Antibonding combinations of the atomic hybrids may then be generated to form a valence virtual orbital space similar to the results obtained in the present paper. Other methods in this fourth and final classification of valence orbital extraction are the effective atomic orbital method49 and the adaptive natural density partitioning (AdNDP) method.50 The famous “atoms in molecules” theory of Bader51 is not categorized here. While certainly well within the spirit of

II. GENERATION OF VALENCE VIRTUAL ORBITALS Two types of orbitals3,5 with predominantly atomic character can be extracted from the occupied and virtual orbital spaces of molecular orbital calculations. One set is very atomic in nature: these are termed the quasi-atomic minimal basis orbitals (QUAMBOs or QUAOs), which are localized onto individual atoms and are as close as possible to the s, p, d orbitals of the free atoms. The other set is molecular in nature: the occupied canonical orbitals of the Hartree−Fock wave function supplemented by valence virtual orbitals (VVOs), which are unoccupied and typically antibonding valence molecular orbitals. This second set (occupied plus VVOs) should, to the maximum extent possible, be formed from the core and valence orbitals of the atoms comprising the molecule. The number of orbitals in both sets is equal to the sum of core plus valence orbitals for all atoms. Both sets (QUAMBOs or occupied plus VVOs) exactly span the converged occupied orbitals of the molecular SCF and the same valence subspace of the molecule’s virtual orbitals. II.A. Separation of the VVO Space. The present paper involves only applications for the VVOs, which are hoped to provide a quantification of the frontier orbital concept. Applications for the atom-localized QUAMBO orbitals are given elsewhere.5,6 In case only VVOs are needed, the QUAMBO/VVO algorithm of Lu et al.3 can be greatly simplified, as shown below. The algorithm is also now recognized as an instance of the singular value decomposition (SVD), whose mathematics is discussed elsewhere.34−38 The algorithm’s first input is a set of ordinary canonical molecular orbitals (CMOs) from some kind of self-consistent field (SCF) calculation, such as closed shell or high- or low-spin coupled open shell SCF or even multiconfigurational SCF. The virtual subspace of the CMOs is given as a linear combination of atomic orbitals (LCAO) expansion, ϕv = ∑μ χμCμν. The user may freely choose their favorite quantitative basis set χμ for the molecular calculation. The goal is to divide this molecular virtual space into an internal (valence-like) part and its external orthogonal counterpart. C

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A The second input is a set of accurate atomic minimal basis set orbitals (AAMBS) for every atom in the molecule, denoted collectively as Ai, or as AAa when referring to a particular atom A with its various occupied orbitals a (index a includes both core and valence, but not empty atomic orbitals): AAa =

application of the rotation T to the original CMO virtual space separates it into the desired internal and external orbital spaces: ψw =

v

If the diagonalization (or SVD) routine orders λ into descending order, the first ψw are the desired internal orbital space while the remainder are the external orbital space. The number of internal orbitals created is the sum of all AAMBS orbitals minus the number of orbitals occupied in the molecular wave function. A correct program will exhibit a clear drop off in λ from near unity to around 0.1−0.2 at the boundary of these two virtual subspaces.5 As an example, the closed shell wave function for H2SO4 has 25 occupied molecular orbitals, and 2 × 1 + 9 + 4 × 5 = 31 AAMBS orbitals. When these 31 (preorthogonalized) atomic orbitals are used in the SVD process, 31 − 25 = 6 values of λ are found near unity, and their 6 columns of T generate the desired antibonding valence virtual space. The number of external orbitals depends on the user’s choice of the molecular basis set: 25 of the remaining columns of T have small but nonzero singular values, while the remainder vanish since the number of nonzero SVD diagonal elements of a rectangular matrix is equal to its smaller dimension. The final step in preparation of the valence virtual orbitals (VVOs) is a “pseudocanonicalization” to produce uniquely defined orbitals with an energy estimate. The pseudocanonicalization is accomplished by diagonalization of the relevant Fock operator,4,5,59 within its internal virtual and its external virtual blocks, separately. The pseudocanonical internal orbitals are the desired VVOs, which (a) are as valence-like as possible, due to the SVD’s maximal alignment of the CMO virtual orbital space to the AAMBS, (b) have a “pseudoeigenvalue” that represents an energy estimate for each VVO, namely, their Fock operator expectation values, (c) possess the full symmetry of the molecule (assignable to some irreducible representation of the point group), and (d) are typically delocalized across the entire molecule, just as occupied canonical orbitals are. The program creates a complete set of MOs by joining together three subspaces: the converged occupied SCF orbitals, then the pseudocanonical VVOs, and finally the pseudocanonical external molecular orbitals. Because the occupied orbitals are obtained by closed or open shell SCF, the only nonzero Fock elements are those connecting the VVO and external orbital spaces (full diagonalization of the entire virtual block of the Fock operator just regenerates the virtual CMOs of the SCF calculation). Implementation of VVO generation is straightforward. The only unusual requirements are storage of the AAMBS orbitals for all atoms and the ability to compute overlap integrals between their GTO expansions and the molecular basis set. The computational requirements are just a few matrix multiplications and diagonalizations, so VVOs may be extracted after a molecular SCF calculation converges, in far less time than a single SCF cycle. II.B. Specification of the AAMBS Atomic Orbitals. The original algorithm3 suggested using atomic orbitals expanded in the same basis set that is chosen for the molecular calculation. This is impractical, since there are a very large number of

∑ χνA A νa ν

The exact nature of the AAa orbitals is given in the following section II.B. For the moment, it is only important to know that each of these is expanded in an auxiliary basis set χAν , centered on just one atom A. The atomic orbital expansions, Aaν, are pregenerated and internally stored in the program. The original paper3 did not explicitly identify the algorithm as an SVD and therefore did not stress that the AAMBS side of the SVD must be preorthogonalized, perhaps by a symmetric orthogonalization, A*j =

∑ Ai(Sij)−1/2

Sij = ⟨Ai |Aj ⟩

i

where {Ai} is the union of all AAa. Note that each diagonal block of S is a unit matrix by virtue of the orthogonality of any atom’s AAMBS orbitals, but S contains nonzero interatomic overlaps. The ∗ emphasizes the requirement of orthogonality for the A* orbitals (i.e., ∗ is not complex conjugation). Since the SVD step below deals with the AAMBS side and CMO side as two entire function spaces, the preorthogonalization procedure chosen does not matter unless one is also interested in forming the atomically localized QUAMBOs. Note that the CMO side of the SVD below is already orthonormal. Steps i and ii of the algorithm3 represent a singular value decomposition of the overlap between the molecule’s virtual orbitals and the atoms’ preorthogonalized AAMBS space, svj* = ⟨ϕv|A*j ⟩ =

∑ ∑ Cνμ† ⟨χμ |χνA ⟩A ν*j μ

ν

Note that it must be possible to compute the rectangular overlap matrix between two different sets of Gaussian type orbitals (GTOs), shown in the center of this formula. The SVD is this generalized eigenvalue problem, T†s∗U = λ

The orthogonal transformations T and U represent rotations within the virtual orbital space and the AAMBS space, respectively, to create new functions within each space that are parallel to each other, with overlaps λ. The number of nonzero λ is at most (and usually equal to) the smaller dimension of the molecular virtual space and the AAMBS space (the latter should be smaller). Clearly, the SVD between the core and valence atomic orbitals A* and the CMO space must pick out the most atomic-like orbitals lying in the latter space. Alternatively, the SVD equation may be squared, to create an ordinary eigenvalue problem, (T†s∗U)(T†s∗U)† = T†BT = λ 2 ,

∑ ϕvTvw

where B = s∗(s∗)†

which is step ii in the original algorithm.3 Diagonalization of the symmetric matrix B means that an SVD subroutine need not be available. In the present work, which does not involve the use of atomlocalized QUAMBO orbitals, steps iii−vi in the original algorithm which generate QUAMBOs may be skipped. Instead, D

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 1. Summary of Configurations and Terms Used in the State-Averaged MCSCF Calculations Used To Prepare the Transition Metal Atomic Orbitalsa atom Sc Ti V Crb Mn Fe Co Ni Cu Zn

term d1s2 d2s2 d3s2 d4s2 d5s2 d6s2 d7s2 d8s2 d9s2 d10s2

2

D◁ 3 F◁ 4 F◁ 5 D 6 S◁ 5 D◁ 4 F◁ 3 F◁ 2 D 1 S◁

d2s1 d3s1 d4s1 d5s1 d6s1 d7s1 d8s1 d9s1 d10s1

term

ΔE

atom

4

33 19 6 22 49 20 10 1 32

Yc Zrb Nb Mob Tc Rub Rh Pd Agc Cd

F 5 F 6 D 7 S◁ 6 D 5 F 4 F 3 D 2 S◁

term d1s2 d2s2 d3s2 d4s2 d5s2 d6s2 d8s1 d9s1 d9s2 d10s2

2

D◁ 3 F◁ 4 F 5 D 6 S◁ 5 D 4 F◁ 3 D 2 D 1 S◁

d2s1 d3s1 d4s1 d5s1 d6s1 d7s1 d9s0 d10s0 d10s1

term

ΔE

4

31 14 3 31 7 21 10 19 87

F 5 F 6 D◁ 7 S◁ 6 D 5 F◁ 2 D 1 S◁ 2 S◁

An arrow indicates the lower energy L−S term, by experiment.61 The experimental energy separation ΔE between the lowest J levels of these L−S terms is given in kcal/mol. bExperimentally, there is a second L−S term in the lower energy configuration which lies between the two states that are averaged: Cr and Mo’s d5s1 5S, Zr’s d2s2 3P, and Ru’s d6s2 3F. These terms are ignored since the intent is to average configurations. cExperimentally, Y’s s2p1 2P and Ag’s d10p1 2P states lie between the two L−S states averaged but are ignored since the intent is to obtain s and d orbitals. a

configurations and terms chosen for the averaging, guided by experimental data for neutral atoms.61 The two chosen L−S terms were averaged with equal weights, and 5-fold radial degeneracy was imposed on the d orbitals. The metal calculations are technically complex and are described briefly. Since the two terms to be averaged sometimes have different spin multiplicities and usually possess complicated angular momentum couplings, a full CI determinant62 based state-averaged MCSCF program was used to perform the atomic calculations. State-specific true single configuration calculations are easily done with the occupation restricted multiple active space (ORMAS) determinant program,63 which can choose a specific occupancy such as s1dm+1. Analogous state averaged SCF-level calculations are not always possible. Even if ORMAS is used to rigorously select only two of the three possible electron configurations, a given L−S term still might not be single configurational. Consider the yttrium atom, using ORMAS to select only s2d1 and s1d2 determinants, thus ignoring s0d3 determinants. The s2d1 configuration contains only a 2D term, which is Y’s ground state. The higher s1d2 configuration contains sufficient determinants to form 4F (the configuration’s lowest L−S term), along with 2F, 4P, 2P, 2G, 2D, and 2S terms. Since 2D appears again, state averaging yields a true single configuration 4F term but a two configuration 2D term. Since even using ORMAS cannot preclude some multiconfigurational character, the metals were computed using all determinants arising from all valence electrons in a (s, d) active space, thus allowing up to three configurations to mix. So the actual yttrium AAMBS orbitals also involved determinants from s0d3 which contains (among other L−S terms) an additional 2D and 4F term, so that the Y atom’s SAMCSCF contains mixing of three 2D terms and two 4F terms. However, the orbital averaging is over the lowest 2D and 4F states, which are largely but not exactly single configurational in character. It is hoped that the valence d and s orbital sizes that result from these configurationally averaged metal atom calculations are similar to the orbitals belonging to a metal found in a molecule. The ground state of metal atoms in the gas phase is often (n − 1)dm, ns2. The decrease in relevance of this configuration after incorporation into molecules and the difficulty of dividing the valence orbitals from other orbitals for atoms found deep in the periodic table have been discussed by Schwarz.64−66

possible basis sets, although one can imagine generating each atom’s orbitals by “on the fly” atomic SCF calculations. A more practical solution is to choose some particular basis set to expand the atoms, once and for all, storing the resulting core and valence atomic orbitals. This choice leads to the aforementioned requirement to compute overlaps between the AAMBS’ Gaussians and the molecule’s Gaussians. The accurate atomic minimal basis set (AAMBS) orbitals are chosen to be nonrelativistic SCF orbitals or state-averaged multiconfigurational SCF (SA-MCSCF) orbitals, for neutral atoms. At the present time, all atoms are stored down to Xe (Z = 54). Beyond Xe, atomic orbital size changes due to scalar relativity would have to be considered. The number of AAMBS orbitals must be defined: Alkali and alkali earths are considered to be s-block, with valence ns orbitals. Main group elements are considered p-block, with valence ns and np orbitals. Transition metals are considered d-block, with valence (n − 1)d and ns. All filled (core) orbitals below these valence AOs are also part of each atom’s AAMBS. Of course, the working basis set for accurate molecular calculations should contain additional functions for certain higher energy orbitals that mix to some extent into the molecule’s occupied orbitals, for example, the np orbitals of alkali or transition metals or the nd functions that may be important in hypervalent main group compounds. At present, such functions are not part of the AAMBS. The s-block and p-block AAMBS orbitals are obtained by closed- or open-shell SCF calculations on their free atom ground states, imposing radial degeneracy on the p orbitals. It is well-known that the radial size of p orbitals is insensitive to the L−S coupling for open shell pm configurations, so the ground state term is used. Transition metals, however, require more consideration. It is well-known that molecular chemistry may involve several lowlying configurations: (n − 1)dm, ns2; (n − 1)dm+1, ns1; or (n − 1)dm+2, ns0. These atomic configurations have markedly different radial expectation values (sizes) for the d shells, since the extent of screening by the smaller (n − 1)d AO increases with its increasing occupancy.60 Accordingly, the AAMBS chosen for the metals is an average of the two lowest atomic configurations: usually but not always the gaseous state metal’s two lowest energy configurations are (n − 1)dm, ns2 and (n − 1)dm+1, ns1. The SA-MCSCF calculations averaged two Russell−Saunders terms, using the lowest term from the two lowest energy configurations. Table 1 gives details about the E

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A It remains to specify the GTO expansions used for each atom. It is desirable to use a large number of Gaussians, to effectively be at the atomic basis set limit, so that the calculations are essentially free of yet another basis set truncation choice. Accordingly very large expansions in terms of GTO primitives are used. The s- and p-block elements up to Ne (Z = 18) use even-tempered (ET) Gaussians,67 after which well-tempered basis set (WTBS) Gaussians are used, up to Xe (Z = 54).68−71 The s-block and p-block AAMBS orbitals are the following: H−He, 8s; Li−Be, 14s; B−Ne, 14s,7p; Na−Mg, 18s,9p; and Al−Ar, 18s,12p ET primitives; and K−Ca, 26s,16p; Ga−Kr, 26s,20p,14d; Rb−Sr, 28s24p22d; and In−Xe, 28s,23p,17d WTBS primitives. The metals use Sc−Ni, 26s,17p,13d; Cu−Zn, 26s,17p,14d; Y−Tc, 27s,20p,17d; and Ru−Cd, 28s,20p,17d WTBS primitives.71 All of these primitive sets approach the ground state SCF energies to within 1 mhartree or better. A general contraction of these SCF or SAMCSCF orbitals to the minimal basis set sizes defined above produced the final AAMBS orbitals. Future work could extend the VVO methodology to heavy atom all electron calculations which require scalar relativity (beyond Xe). For instance, the molecule might be treated by the infinite order two-component72 or a finite order Douglas− Kroll−Hess approximation.73−75 The inclusion of scalar relativity is well-known to result in the contraction of s orbitals, a smaller contraction of p orbitals, and expansion of d or f orbitals.76,77 It seems likely these shape changes will require storing an alternative set of AAMBS orbitals, prepared by including IOTC-level relativity in the same type of atomic calculation just described. A future extension to relativistic model core potential78,79 calculations should then be feasible: the semicore and valence orbitals retained in any MCP calculation should have the correct radial shape so that one would simply omit any core orbitals dropped by the MCP from the AAMBS side of the SVD. Extension to relativistic or nonrelativistic effective core potentials80 is problematic, since ECP orbitals are typically nodeless and therefore do not closely resemble the true radial shapes of valence atomic orbitals. II.C. Implementation. The generation of VVOs (and also QUAMBOs) for any molecule containing the elements H−Xe is implemented in the GAMESS quantum chemistry package81−83 as a user selectable option. The generation of VVOs is switched on by a single keyword. The process is entirely automatic, since no user assistance is needed to count the number of core and valence orbitals (as in the H2SO4 example given above). VVOs may be obtained for closed shell RHF, high- and lowspin open shell ROHF, and MCSCF wave functions. The number of VVOs to be found decreases as more orbitals are occupied (the sum of occupied and virtual valence molecular orbitals always equals the total number of occupied orbitals from every atom). Localization of the VVOs is also possible (see section III.C below).

Molecular geometries are mostly chosen to be tightly optimized RHF/ACCT structures. In a few instances, structures found by RHF/ACCT are not as close to experiment as would be desired. Thus, DFT with the TPSS functional was used for MnO4−1 as RHF has too short MnO bonds, full valence space multiconfigurational self-consistent field (MCSCF)93 was used for dioxirane as RHF has a too short OO distance, an 8 e- in 8 orbital MCSCF was used for obenzyne as RHF has a too short CC triple bond, second order perturbation theory (MP2)94 was used for 2-norbornyl cation as RHF has too long nonclassical CC bonds, MP2 was used for P4 as RHF has too short PP distances, and finally DFT with the B3LYP functional was used for ferrocene as RHF has too long metal/ring distances. In the structures chosen, bond distances are typically within about 0.02 or 0.03 Å of available experimental data. Even when a higher level geometry is used, orbital results are mostly from closed shell RHF wave functions, although some closed shell density functional theory (DFT) results are presented, using the functionals PZ81,95 PBE,96 PBE0,97 B3LYP,98,99 wB97,100 or TPSS.101,102 The collective term DFA (density functional approximation) is used for calculations with these approximate functionals. The collective term KS (Kohn−Sham) is used to compare DFA eigenvalue data to those obtained using the exact (largely unknown) functional or some almost exact functional.103,104 For purposes of evaluating the utility of VVOs in predicting the shape of valence excited state orbitals, the natural orbitals of the relaxed density matrix specific to each excited state are obtained by time-dependent DFT105−107 (TD-DFT) using the B3LYP functional, and the Tamm−Dancoff approximation.105 The choice of closed shell examples in the rest of this paper is just for simplicity: VVOs of similar characteristics may be readily obtained for both low- and high-spin coupled open shell SCF wave functions or after MCSCF. III.B. Basis Set Dependence of the CMO and VVO Energies. Koopmans’ theorem108 asserts that the occupied orbitals of a molecule correspond roughly to ionization potentials from each occupied orbital, affording a qualitative explanation of photoelectron spectra. That is, the energies of the canonical HOMO, HOMO − 1, ... of a neutral molecule M provide information about the energies of the cations M+1. This is well accepted by the community and is borne out by data shown below: the HOMO, HOMO − 1, ... concept is concretely realized by the occupied canonical SCF (or DFA) orbitals. Textbooks discussing the Hartree−Fock equations also contain a Koopmans’ theorem interpretation of the unoccupied canonical orbitals as being related to ionic states M−1. However, again correctly, this is nearly universally accompanied by the statement that the eigenvalues of the lowest canonical molecular orbitals LCMO, LCMO + 1, ... are not even qualitatively useful for estimating electron affinities. Nonetheless, it is still common to encounter the terminology HOMO− LUMO gap, although only the occupied eigenvalues are qualitatively useful. It is easy to give a physical reason for the poor valence properties of many unoccupied CMOs. As already mentioned, a conventional interpretation for the empty orbitals is as a place to attach an additional electron, meaning they can represent (N + 1) electron states. The anions of closed shell molecules typically have only loosely bound electrons, at best. Alternatively, virtual orbitals provide approximate excited states

III. APPLICATIONS INVOLVING VALENCE VIRTUAL ORBITALS III.A. Methods. Orbital results given in this paper are from closed shell self-consistent field calculations (RHF) with the aug-cc-pVTZ basis set 84−87 (ACCT), unless otherwise indicated. Since no ACCT all-electron basis set exists for iodine and ruthenium, a similar TZ-quality Sapporo basis set88−90 was used for all atoms in IF3 and Ru(bpy)3. Orbitals are drawn with the MacMolPlt visualization program.91,92 F

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 1. Energies of RHF canonical MOs (in color) and valence virtual orbitals (in black) for five different basis sets. Eight of the panels used the same basis sets, as shown in the upper left panel. The energy scale increases going to the right, in hartree units. See text for discussion of each molecule. For ethane, only canonical MO positions are shown, and its RHF results are compared to DFT. The DFT LCMO of ethane (TPSS functional) is depicted at a reduced contour value of 0.025 bohr−3/2 appropriate for a low amplitude, spatially large MO.

canonical virtuals can only be valence antibonding orbitals and may be used as such. The quantitative comparison of empty CMOs to VVOs in this section uses four Pople-style basis sets110−112 and one essentially saturated basis set. The first basis is the very small 631G, containing only s,p functions. This is improved by adding polarization to all atoms, namely, 6-31G(d,p). The next two improvements consist of either adding more polarization or else diffuse s,p functions: 6-31G(3df,2p) or 6-31++G(d,p). Generally speaking the polarization improvement provides a larger energy lowering, but the diffuse basis’ lower CMO eigenvalues group closer to those of the final basis set, which is aug-cc-pVQZ84−86 (called ACCQ herein). This final basis is very large, containing both diffuse augmentation and numerous polarization functions up to the level of g AOs. These five basis sets are used in all panels of Figure 1, in the same order and as labeled in the upper left panel. MnO4−1 is an exception, as described below. This section focuses on the VVO energies, while their shapes are considered in the following section III.C. The first two panels in Figure 1 show only the Fock and TPSS canonical orbital spectra of ethane. Ethane has rather uninteresting valence σ* antibonding orbitals, whose VVO energies are not included in Figure 1, to concentrate on the characteristics typical of the canonical orbital eigenvalue spectrum of most molecules. Note that the occupied orbitals readily converge with respect to basis set: two of the five lines shown for C2H6 below zero are doubly degenerate, accounting for all seven valence pairs. On the other hand, it is clear that the

M*, for example, by singly exciting electrons into them (N electron states). Some of the M−1 or M* states correspond to Rydberg states, or even to a completely detached electron plus M or M+1, and the kinetic energy of any such free electron is not quantized. Whether viewed as representing M−1 or M*, such orbitals have physical meaning as scattering resonances and have been described as “being infinitely extended, with only a few orthogonality wiggles in the molecular neighborhood”.103 In other words, the canonical Fock spectrum of a neutral molecule, like the textbook H atom case, should contain a continuum of unbound electron states, which are asymptotically plane waves, starting at energy zero. Note that this is true for all neutral molecules so that the LCMO spectrum in large basis sets starts from zero, as will be demonstrated below. It should be remarked that the Gaussian orbitals used in everyday quantum chemistry applications are not particularly good plane waves! Of course, other orbitals in the virtual space do have valence character, and their extraction from this fairly uninteresting continuum, and utility in chemistry, is the point of this paper. As we will see, in some cases like C60 or indigo, there may be one (or more) empty Hartree−Fock CMOs below the continuum starting from zero, and these few discrete virtuals have valence character. Note that smaller basis sets may not contain any diffuse GTOs so that their virtual orbitals may exhibit reasonable valence character.109 In the limit of a minimal basis set (such as is used in extended Huckel theory), G

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

The small anion MnO4−1 has occupied orbitals at higher energy than typical neutrals, just as cations have deeper than usual HOMO, HOMO − 1, etc. The transition metal’s 6-31G contains s, p, and d functions. Polarization is shown in its panel by specifying these for its two atoms (Mn, O) in that order. The diffuse augmentation of the Pople-style sets is created by even-tempered ratios.114 Although there are drop-offs in the position of the anion’s LCMO eigenvalue as the basis grows, even the very large ACCQ basis set’s first canonical eigenvalue is well above zero, at +0.114 hartree. This is due to the negative charge, as a true continuum for M/e− type virtuals can only be obtained by using GTOs which are far more diffuse than normal. The smallest exponent found in the ACCQ basis for Mn is an s function with exponent 0.0128. Adding an extra p GTO shell to the central metal with the small exponent 0.001 inserts a new triply degenerate t2 LCMO at +0.034 hartree. If the extra function is instead an s GTO with the even smaller exponent 0.0001, the new LCMO is a singly degenerate a1 orbital at +0.016 hartree. As in this case, the lowest virtual CMOs of other anions seldom possess useful valence antibonding character. In contrast, the six VVOs of MnO4−1 have similar pseudoeigenvalues and shapes for all bases. The e, t2, and a1 symmetry VVOs represent antibonding combinations of metal d orbitals with the oxygens: the first degenerate pair involves antibonding interactions with O lone pairs, while the top four are MnO σ* in nature. The final column of Figure 1 returns to neutral molecules, which were selected as likely to contain a low-lying unoccupied orbital. In the singlet carbene HCCl, there should be a p orbital at carbon, perpendicular to the molecule. The first four basis sets do show such an a″ symmetry shape, of valence 2p size. However, the largest basis once again has its LCMO being very diffuse and even the incorrect a′ symmetry. In contrast, the three VVOs are very sensible. All basis sets show the same orbitals, in ascending order: the C 2p, CCl σ*, and finally the CH σ*. Dioxirane (H2CO2)115 is important in ozonolysis and has a low-lying OO σ* orbital due to ring strain, which is occupied by almost 0.1 e− in a full valence MCSCF wave function. Thus, the LCMO of the first three basis sets are this b2 symmetry OO σ*, but as usual for the two largest basis sets, the LCMO loses its valence character. In energy order, the five VVOs are the OO σ*, an antisymmetric combination of the two CO σ*, symmetric and antisymmetric CH σ*, and finally the symmetric CO σ* combination. Pyrazine116 has low-lying π* orbitals, so this moiety is often incorporated into donor/acceptor charge transfer compounds, although curiously, the isolated molecule has only a very small electron affinity. Nonetheless, small basis sets have two unique low-lying CMOs, which are its π* accepting orbitals, with a noticeable gap before reaching the remaining CMOs. These two orbitals are swamped out and lost in the final two basis sets which add diffuse orbitals in the same energy range. The first three VVOs are π* orbitals (two are much lower than the third), beneath the ten σ* antibonds, for all basis sets. For cations or for molecules with a low-lying unoccupied valence orbital and particularly when using modest basis sets, the LCMO may have a shape corresponding to the expected low-lying valence orbital. The utility of the next higher CMO, namely, the LCMO + 1, is minimal (although pyrazine has two useful empty CMOs in small bases). Thus, several decades ago when use of small basis sets was the norm, the LCMO could often be taken for the empty frontier orbital. Modern quantum

unoccupied CMOs change greatly with the basis set and do in fact start to resemble a continuum for ACCQ. By happenstance, the LCMO in all five basis sets has symmetry type a1g, but for the two final bases with diffuse functions, this orbital becomes spatially very large. This is because the plane waves with kinetic energy just above zero are being mimicked by the large (diffuse) GTOs added in the final two basis sets. Note that use of TPSS density functional theory also reaches a continuum virtual space limit for large bases, but the energy values do shift: the occupied orbitals are all raised substantially, and the LCMO for the final two bases with diffuse functions drops slightly below zero. The occurrence of some empty canonical orbitals below zero is not uncommon in DFA calculations. This paper does not explore in depth what DFA functional might be quantitatively the most accurate: the remaining panels in Figure 1 show eigenvalues and pseudoeigenvalues from the Hartree−Fock operator. Diborane’s CMO spectrum in the bottom panel of the first column of Figure 1 is similar to ethane, except the symmetry type of the LCMO now changes (as indicated) with the basis set, particularly when diffuse functions are added. Of the eight molecules shown in Figure 1, only ethane and H3O+ have the same symmetry for the LCMO for all bases. As noted above, this is just accidental for ethane, as its LCMO undergoes a big spatial expansion. The pseudoeigenvalue for each VVO is placed on top of the CMO spectrum using thick black lines. Notice that the VVO positions converge quickly with respect to the basis set improvements! Their shapes also remain nearly identical and, if visualized, are seen to be of antibonding valence character. Removal of the seven VVOs from the virtual space means that the remaining external orbitals of diborane shift around somewhat. This is not illustrated, and for the largest atomic bases, the external pseudoeigenvalues remain a nearcontinuum, just missing seven values. The middle column of Figure 1 contains ions. For H3O+, the small size means that the LCMO is a symmetric combination (a1 symmetry) of all three OH σ* antibonds. Even for the largest basis set, the LCMO looks similar to the LVVO, although the latter clearly has a more systematically converging energy. Thus, the Koopmans’ theorem argument that the virtual orbitals of M+ can be interpreted as electron attachment sites holds up, at least for the few CMO levels appearing below zero. The higher CMO levels begin to fill in a continuum above zero, as it remains true that many of the attached electron states correspond to an unbound additional electron scattering from H3O+. The two VVOs are the a1 and degenerate e symmetry combinations of the OH σ* antibonds and, as for neutrals, have similar energy and shapes no matter what basis is used. The tendency for small cations to exhibit valence character in their low-lying CMOs provided the motivation for IVO13 and MVO14 types of alternate virtual orbitals. The much larger 2-norbornyl cation (in its nonclassical geometry113) has its charge spread over three carbon atoms, so now the LCMO in the two largest basis sets is a spatially diffuse orbital (and has different symmetry than the three smaller sets). As is typical of cations, the unoccupied CMO spectrum begins below zero. The HOMO and the first two VVOs are chemically interesting and will be discussed and illustrated below in section III.C. This molecule also contains numerous CC and CH σ* antibonds, forming a dense cluster of black lines, with the highest VVO being the σ* of its shortest CC bond (at the base of the three-membered ring). H

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A chemistry applications routinely use large basis sets and clearly cannot be expected to yield meaningful CMOs, as spatially large empty CMOs commonly occur. An illustration of what is meant by “spatially large” CMOs is given for C2H6 in Figure 1, and additional examples will be presented below. In contrast, all VVOs are found to occur at nearly the same pseudoeigenvalue, no matter what basis set is used and no matter what the molecular charge is. The shapes of VVOs, too, are also nearly independent of the basis set. The set of VVOs represents all missing valence antibonding orbitals, not just the lowest one, as will be convincingly demonstrated in the next section. The basis set dependence of the lowest unoccupied canonical orbital’s eigenvalue is not well-known within the chemical community perhaps because workers may use only a single basis set for any given theoretical study. However, the results in Figure 1 are unquestionable. As the basis set grows, the LCMO of neutral and anion systems is expected to drop to zero. This has been noted before, in atoms and small hydrides,117 and clearly is a very general result, which casts profound doubt on any physical significance for the LCMO. III.C. VVO Shapes. A few illustrations are given in this section to prove the assertion that VVOs possess antibonding valence character. Additional examples are given in subsequent sections which focus on possible applications for VVOs. VVO shapes are found to be almost entirely independent of the basis set, so only ACCT orbitals are shown in the remainder of this paper. Because of the “pseudocanonicalization”, which assigns a diagonal Fock (or DFA) operator energy to each VVO, the VVOs are usually delocalized across the entire molecule. Since the number of VVOs is similar to the number of occupied valence orbitals, there is no trouble using ordinary localization algorithms to localize the VVOs. This is in contrast to the great difficulty of localizing within the external virtual space, which contains many orbitals from every atom.40,53−55 As implemented in GAMESS, any localization takes care to leave the total wave function invariant: the filled, partially filled in the case of open shell wave functions, and virtual valence orbitals are localized separately. Preserving the occupied orbital space(s) means the localized orbitals represent the same total wave function while often providing greater interpretation. If the molecule of interest does not possess any inherent delocalization (aromatic rings, three-center bonds, ...), the result of localizing the occupied orbitals will often be twocenter bonds and one-center lone pairs. In this case, the localization of the VVO space results in a set of two-center antibonds, with one antibond for every bond. For certain applications, such as using VVOs as starting orbitals for MCSCF (see below), it is convenient to localize as much as possible but without destroying orbital symmetry. This “symmetry constrained” localization is available as a user option in GAMESS. The localization of VVOs, possibly with such an imposed symmetry constraint, is available for all popular localization procedures: Edmiston/Ruedenberg self-energy,118,119 “Boys” dipole,120,121 and Pipek/Mezey population122 criterions. The Edmiston/Ruedenberg procedure is used for all localized orbital results presented here. For cases where the unoccupied frontier orbital(s) are expected to be clear-cut, the pseudocanonical VVOs themselves may already be fairly well localized. This is the case for the 2norbornyl cation, as can be seen in the top row of Figure 2. After filling all ordinary two-center bonds in the rest of the molecule, there are two electrons left over. This pair occupies the HOMO of the cation, which is a bonding orbital in the

Figure 2. This figure and all others (unless otherwise indicated) uses a contour increment appropriate to valence-sized orbitals, namely, 0.05 bohr−3/2, which is the only contour visible in the three-dimensional plots and is the increment in these two-dimensional slices. Electron occupancies (2 or 0) and some symmetry labels are shown. See text for discussion of the three molecules.

nonclassical CCC triangle involving three p orbitals, one from each carbon. The LVVO and LVVO + 1 are antibonding in the same region. Symmetry constrained localization, if it were done, would remove most of the amplitudes on more remote atoms, such as can be seen on the top hydrogen atoms in the LVVO (middle image). The hypervalent IF3 compound is shown in the center of Figure 2. Counting the 3 lone pairs at each F atom, there are a total of 16 valence orbitals in IF3, only two of which are unoccupied VVOs. Three of the canonical occupied (labeled with occupation 2) and both VVOs (occupation 0) are shown in the top row for IF3. The canonical occupied orbitals mix up the orbitals corresponding to the two center IF bond and the three center FIF bonds. Symmetry constrained localization, illustrated in the next row, clarifies the bonding. The two VVOs are of different symmetry, so symmetry constrained localization in the virtual space does not change them. However, the occupied orbitals separate into a two-center IF bond, which is paired with the a1 symmetry VVO, both shown at the right. The three orbitals on the left correspond to a typical three-center, four-electron bond:123 an axial 5p of iodine bonds to a fluorine hybrid with a phase change at I caused by its 5p (symmetry b2), followed by a nonbonding orbital with small amplitude on iodine which is a symmetric combination of the ligand hybrids; both are occupied. The final three-center orbital is the antibonding counterpart of the first, and this VVO is not occupied. Two more occupied localized orbitals at iodine are of interest but are not shown: there is a 5p perpendicular to the plane and a largely 5s lone pair in-plane. The latter orbital, in I

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A accordance with the VSEPR model,124 is responsible for the slight bending of the IFI axis toward the equatorial F. The final example in Figure 2 is tetrahedral P4,125−127 seen in the bottom row of Figure 2. The P4 is oriented with a triangular face at the bottom, and one atom at the top, so the bottom right corner is the projection of atoms above and below the plane drawn. VVOs possess the full symmetry of the molecule so that the six PP antibonds transform as two triply degenerate levels: t1 at +0.138 and t2 at +0.213 hartree for RHF/ACCT. Thus, the VVOs are completely delocalized. This is a case where full localization (without symmetry constraint) produces the most chemically relevant picture: the occupied orbitals give four equivalent lone pairs and six identical curved bond pairs. The localization of the VVOs gives six equivalent antibonds, which are also curved outward due to ring strain. Note that under unconstrained localization in singly bonded molecules, VVOs localize to one two-center antibond for every localized two-center bond, whereas lone pairs are found only in the occupied space. Because P4 possesses only lone pairs and two center bonds, its orbitals are not as interesting as the 2norbornyl+1 ion or IF3, apart from localization clearly displaying its ring strain in curved bonds and antibonds. The examples in this section show that even for a large basis set like ACCT, VVOs are found to be the same size as occupied valence orbitals, to be antibonding combinations of various valence atomic orbitals, and to provide all of the antibonds that should exist. III.D. MO Diagrams. Since the two previous sections indicate both VVO shapes and VVO energy values are chemically sensible, as well as basis set independent, it is reasonable to combine the ordinary occupied CMOs (HOMO, HOMO − 1, ...) and the unoccupied VVOs to construct MO diagrams of bonding and antibonding levels. Four examples are considered in this section. The famous molecule ferrocene128 contains sufficiently interesting bonding interactions to serve as a first example. The staggered D5d conformation found in the solid state129 is used to enable comparison to recent photoelectron spectra of the solid.130,131 Figure 3 shows a ferrocene MO diagram, constructed with B3LYP occupied orbitals as well as the VVOs that pseudocanonicalize the same DFA operator. An SCF-level RHF Fock operator produces qualitatively identical MO shapes but orders the occupied orbitals differently, whereas B3LYP gives good agreement with experimentally known occupied orbital orders130,131 (and to similar DFA calculations130,131 that were used to interpret the experimental results). RHF also gives three nearly isoenergetic doubly degenerate levels as the first three unoccupied VVOs,132 while B3LYP separates out the e1g as the ferrocene LVVO. Many diagrams for ferrocene construct this molecule from Fe2+ and two cyclopentadienyl anions (Cp−), but in accordance with section III.B above, the Fe2+ cation orbitals are very low and the Cp− anion orbitals are too high. Figure 3 suggests that using open shell B3LYP energies for neutral Fe and Cp· is a more reasonable way to construct the MO diagram for ferrocene.133 The Cp· moieties are far enough apart that formal bonding and antibonding combinations (with respect to the inversion center) of its orbitals have almost the same energy: these are slightly split at the very bottom, but their twinning often results in only thicker lines, for example, in the numerous orbitals at the top of the ferrocene diagram. Note that a successful energy analysis of ferrocene has been made using either neutral Fe/2Cp· or ionic Fe2+/2Cp−1 fragments.134

Figure 3. MO diagram of D5d ferrocene, with orbital energies from B3LYP, given in hartree units. Red levels indicate the most relevant orbitals composed of metal s and d orbitals and ring π and π* orbitals. Species names (in blue) separate occupied orbitals from the VVOs. The inset illustrates the highest red VVO: Fe 4s antibonding to ring π. Ferrocene’s relevant orbitals are its occupied CMOs a1g (−0.408), a2u (−0.359), e1g (−0.267), e1u (−0.248), a1g (−0.227), e2g (−0.198) and its empty VVOs e1g (+0.005), e2u (+0.094), e2g (+0.107), a1g (+0.199). Doubly degenerate orbitals are shown as pairs, but the 5-fold degeneracy of the metal 3d is not indicated.

The orbitals shown in red in Figure 3 form the ferrocene molecule: π levels from Cp· and the valence s/d from Fe. The high symmetry of ferrocene leads to a relatively small number of metal/ring interactions, shown as the red diagonal lines. These interactions are in perfect accord with the first explanation of the occupied orbitals by MO theory.135 That report was entirely based on group theory rather than any actual calculation in 1953! Many of the twinned π levels have different symmetry types than the metal levels in D5d and thus do not interact with the metal. Of particular interest are the symmetric combination of the lowest ring π and the metal’s 4s and 3d0, which all have a1g symmetry in D5d: the 4s makes bonding and antibonding combination with the ring, with the latter being the highest interesting VVO (shown in Figure 3’s inset). Interestingly, the metal 3d0 scarcely interacts, in spite of having the correct symmetry, and so lies at a nearly unchanged energy in the molecule. The formally symmetric combination of the middle π level of Cp· has the correct orbital symmetry to interact with Fe’s dxz and dyz pair (e1g) in a bonding and antibonding fashion: the latter is the LVVO of ferrocene. The formally antisymmetric combination of the middle two π levels has no orbital on Fe to interact with and so just creates an occupied degenerate level at slightly lower energy in ferrocene. The highest π* of Cp· just twins in ferrocene. The rest of the VVO spectrum is evidently sensible: all antibonding metal/π orbitals lie below the many twins of the Cp· σ* VVOs. A very recent XANES experiment131 interprets the pre-edge feature in the photoelectron spectrum of ferrocene as due to an e1g LUMO. The B3LYP/ACCT e1g symmetry LVVO reported here, arising from the antibonding interaction between the metal dxz/dyz and the partly filled HOMO of Cp·, agrees with this experimental inference.131 C60 is a case where the nature of the first two unoccupied orbitals is well established. Alkali doping results136 showed J

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A K3C60 is conducting while K6C60 is insulating, so clearly the first unoccupied level must be triply degenerate. Elementary group theory137 using the eleven L = 5 spherical harmonics (h-type AOs) correctly predicts a HOMO of symmetry hu (in the point group Ih), a lowest empty orbital of t1u, but the remaining t2u level found in the ungerade h spherical harmonics is not the second lowest empty level of C60. Huckel theory and also band structure calculations136 established the two lowest empty levels as first t1u and then t1g. This is confirmed by recent high level calculations138 of the electronic excitation energies, where the lowest 1Gg state arises from HOMO to t1u excitations, with somewhat higher states arising from excitations involving the t1g. The full set of VVOs below 0.4 hartree for C60 is shown in Figure 4. As expected, the LVVO and LVVO + 1 levels are the

The hybrid GGA functional PBE097 (25% Hartree−Fock exchange) just interpolates between RHF and PBE. The rangeseparated GGA named ωB97100 has an orbital spectrum much closer to Hartree−Fock than the other DFA calculations shown. Clearly the gap between HOMO and LVVO (lowest VVO) varies greatly in Figure 4. However, the order of both occupied π and unoccupied π* is unchanged except for a minor reordering of two very close lying π and two very close π*. The quantitative values of the occupied eigenvalues and unoccupied VVO pseudoeigenvalues thus depend on the nature of the DFA operator used, but their energy order and shapes do not. In particular, the LVVO and LVVO + 1 are consistently and correctly predicted to be of t1u and t1g symmetry no matter what functional is used. The valence-like t1u LCMO of C60 is also shown in Figure 4, as a blue line. Clearly, compared to HF, practical DFA calculations show a sharp drop in the LCMO energy, in contradistinction to a rise in the HOMO energy. For RHF there is no other empty CMO below zero, and there is only one more for the range separated ωB97, while the other DFA calculations have multiple empty CMO levels below zero. These shifts are discussed at the end of this section, after the remaining chemical examples. For the moment, this is expected behavior: most DFA calculations show a “DFA upshift” in their occupied orbitals, including the highest ones; this means DFA estimates of the ionization potential are usually too small. Empty DFA orbitals incur roughly the same amount of “DFA upshift” as for the occupied orbitals so that the band gap may be reasonably well predicted. Hartree−Fock levels for empty orbitals tend to be much too high, so in spite of also suffering a “DFA upshift”, the lowest empty DFA orbitals appear below those of HF. Another icosahedral symmetry molecule, namely, B12H12−2, has a 4-fold degenerate HOMO as well as a 4-fold degenerate LVVO! Conceivably, this is unique, as degeneracies of 4 or 5 can only occur in the rare icosahedral point groups. The high degeneracy of both frontier orbitals is a result anticipated in a group theoretical/simple MO calculation on B12 icosahedra performed in 1955.139 To be specific, for RHF/ACCT, the HOMO is at −0.065 (gg symmetry), while the LVVO is at a substantially higher value of +0.553 hartree (gu). Both HOMO and LVVO orbitals have amplitude only within the B12 polyhedron itself but not on the axial BH bonds (not illustrated). The LVVO shows no trace of diffuse spatial character. Note that only very large basis sets yield a HOMO energy below zero in this doubly charged anion. The final example in this subsection is an orbital symmetry diagram, showing the correlation of orbital energies along a reaction path. Figure 5 is a quantitative version of Figure 40 found in Woodward and Hoffmann’s famous book on conservation of orbital symmetry.2 The dihydrogen exchange reaction between ethane and ethylene is symmetry allowed2 but was later shown to have a substantial energy barrier.140 The original WH diagram was drawn qualitatively, using labels S and A for symmetry or antisymmetry with respect to a mirror plane perpendicular to both CC bonds. Lacking any computation, the original WH diagram simply showed three horizontal lines (constant orbital energy) in the occupied space and three more in the virtual space. Figure 5 shows that two of the three occupied orbitals involved in the hydrogen exchange actually rise in energy substantially near the transition state (see the solid lines). This is easy to understand from the elongation of the reacting CH bonds at the saddle point. Of course, the

Figure 4. MO diagram of Ih symmetry C60. The 30 π and 30 π* are shown in red. Closed shell Hartree−Fock and several DFT functionals are used to assign the occupied orbital energies, and the valence virtual pseudoeigenvalues, in hartree units, from ACCT. The orbital ordering is essentially unchanged from the labels at left except as indicated. Variation of the HOMO and LVVO energy is discussed in the text. The blue line indicates the position of the LCMO orbital, whose shape is very similar to the LVVO in C60. In the Ih group, a, t, g, h symmetry orbitals are 1-, 3-, 4-, 5-fold degenerate, respectively.

t1u and t1g, followed by all remaining π* orbitals: all thirty of these appear before any σ* VVO occurs. Most of the highest occupied SCF orbitals are C60’s thirty π orbitals, but the lowest of these intermingle with some of the highest σ orbitals. C60’s Hartree−Fock LCMO is shown with a blue line in Figure 4. The LCMO occurs at ε = −0.0232 hartree, has the correct t1u symmetry, and is only slightly larger spatially than the LVVO which has a nearly identical pseudoeigenvalue. However, the rest of C60’s unoccupied CMO level diagram consists of numerous diffuse orbitals and resembles a continuum starting at zero energy (see ref 138 for a typical C60 CMO spectrum). In addition to ordinary closed shell Hartree−Fock, Figure 4 shows the VVO positions (and the t1u LCMO level) for various types of density functional theory. Two pure DFA functionals have nearly the same orbital energies: the GGA PBE96 and the metaGGA TPSS.101 Both have substantially raised occupied levels (compared to RHF) and substantially lowered VVOs. K

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Hartree−Fock due to the N − 1 nature of the KS potential compared to the N-electron mean field potential of Hartree− Fock. Their conclusion is the KS HOMO should be a good match to the experimental IP. Additionally, they show that molecules with low-lying states should (and do) possess one or more KS orbitals below zero energy. The lowest orbital occurs at a position above the HOMO that accurately matches the optical gap, which is probably a more precise term than the equivalent term band gap. In regards the topic of the next section, they also demonstrate that TD-DFT calculations within KS theory are expected to show excited states as almost clean occupied to empty individual single excitations. A number of these worker’s expected results are met by the SAOP functional,141 which gives near KS quality results. In practice, most DFA calculations differ from true KS calculations because the DFA exchange-correlation potential is shifted upward by a few eV.103,104 This shift in potential leads to what is called throughout this paper the “DFA upshift” in the eigenvalues, seen for C60 in Figure 4. Because the empty RHF eigenvalue (LVVO or LCMO) is much too high, even with some “DFA upshift”, the empty DFA orbitals occur lower than for RHF. However, the empty DFA orbitals are reasonably positioned relative to the occupied levels, as both “upshifts” should be similar. To be somewhat more quantitative, the following experimental facts for C60 may be considered (see Figure 4, whose scale is in hartree units). The C60 ionization potential is 7.58 eV,142 or 0.28 hartree: clearly the RHF HOMO eigenvalue lies reasonably close to the negative of this IP, while the pure DFA eigenvalues (PBE and TPSS) are “upshifted” considerably. Note that DFA calculations with inclusion of “exact exchange” through either range separation (ωB97) or hybridization (PBE0) interpolate between RHF and pure DFA. The electron affinity of C60 is measured as 2.689 eV,143 or 0.10 hartree. This value has little to do with any orbital position of neutral C60, since the −EA of a neutral corresponds most closely to the IP of its anion,103 but no data in Figure 4 are for C60−1. The band gap of C60 measured in films is 1.7 eV 144 (0.062 hartree) and in liquids 1.82 eV 138 (0.067 hartree), values that agree fairly well with the separation between the HOMO and the LVVO (or C60’s very similar valence LCMO) for both pure DFAs, an indication that both the occupied and empty orbital have the same “DFA upshift”. Clearly, research into more meaningful eigenvalues from DFA calculations is important but is not the central theme of this paper, which is concerned with the separation of the empty orbital space into valence and nonvalence orbitals. Other workers often attempt to correct what this paper terms very nonspecifically as a “DFA upshift” by addressing the selfinteraction error. This can be reduced by hybrids or rangeseparation, but there are efforts to develop new functionals with self-interaction corrections.117,145−147 Other recent papers on DFA eigenvalues include a survey148 of DFA IPs and band gaps suggesting TD-DFT corrections may improve accuracy; the utility149 of KS HOMO positions for Koopmans’-type IPs; discussion150 of the LCMO’s relationship to EAs; and an evaluation151 of range-separation on HOMO and LCMO positions. III.E. Valence Excited States. Frontier orbitals may also be expected to give insight into low-lying valence states of molecules. For example, the ground state’s unoccupied orbitals might be some indication as to the nature of the lowest valence state, upon excitation of an electron out of the occupied space into the empty space. This can be explored by comparing the

Figure 5. Orbital correlation diagram for the dihydrogen exchange reaction between ethane and ethylene. RHF/ACCT has a barrier of 82 kcal/mol. The saddle point, at path distance 0.0 bohr·amu−1/2, has D2h symmetry (imaginary frequency 2334 cm−1 is shown in the inset) . The orbital energies (given in hartree) are colored according to their symmetry type in the C2v symmetry of the reaction path: a1 = blue, a2 = red, b1 = purple, b2 = green. Solid lines represent the six reacting orbitals and correspond to the original Woodward−Hoffmann diagram. WH’s labels of S/A for symmetric/antisymmetric correspond to a1/ b2 symmetry (blue/green). All other valence orbitals are shown as dashed lines.

antibonding orbitals become somewhat less antibonding by the same cause. Neither variance in energy during the reaction is sufficient to cause an orbital crossing, so the reaction remains “allowed”. Note that all spectator orbitals (namely, the CC σ and σ* and non-reacting CH σ and σ*) are also shown in Figure 5, as dotted lines. The examples in this section were chosen to be cases where the nature of the lowest unoccupied orbital was already wellestablished. VVOs matched those expectations and also provided sufficiently accurate predictions for all remaining valence antibonding levels to enable the construction of “semiquantitative” full MO diagrams. Before leaving the topic of MO diagrams, a few remarks about the accuracy of both HOMO and LVVO eigenvalues are in order. Of course, the quality of the VVO pseudoeigenvalues cannot be expected to be any better than the usual Koopmans’ theorem level of accuracy for occupied orbitals. For example, the HOMO eigenvalue only approximates an experimental IP value. In addition, the C60 example in this section shows both occupied and VVO energy values depend strongly on the DFA functional. But one can certainly hope that choosing any particular functional would predict trends in HOMO/LVVO positions within a series of related molecules. Readers interested in the prospects for more quantitative VVO pseudoeigenvalues may consult some of the recent literature. Of particular note is recent work by Baerends and coworkers103,104 on the physical significance of unoccupied true Kohn−Sham orbitals (KS, as opposed to conventional DFA) in relation to electron affinities and the optical gap. They expect the true behavior of KS orbitals to differ fundamentally from L

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 2. Excitation Energies to the Lowest Singlet and Triplet States (in eV) and Oscillator Strengths f, Computed at the TDDFT/B3LYP/ACCT Level at RHF/ACCT Geometries Except As Indicateda no. of empty valence CMOs

no. of VVOS with ε < 0

TD-DFT/B3LYP states

molecule

charge

RHF

B3LYP

RHF

B3LYP

E(S1)

f(S1)

Ru(bpy)3b S7 blueOLEDc S9 p-litmusd S3 indigoe luciferinf arsenicing 3TCh S2 thymine S2 q-litmusd S8 b-litmusd orangei S4 cAMPj pentobarbk benzylic 8l

+2

8

15

27

57

0

0

8

0

15

0

0

3

0

8

0 0 0 0

1 0 0 0

3 3 4 1

0 0 0 0

6 7 6 4

0

0

1

0

3

−1

0

1

0

0

−1 −1

0 0

1 1

0 0

0 0

−1 −1 −1

0 0 0

2 0 2

0 0 0

0 0 0

2.58 2.99 2.67 4.10 4.19 4.73 2.62 4.25 3.83 4.72 4.92 4.93 5.33 2.37 3.63 1.51 2.38 3.31 3.69 3.75 1.26

0.000 0.125 0.000 0.046 0.006 0.039 0.395 0.366 0.033 0.005 0.133 0.000 0.149 0.007 0.210 0.000 0.001 0.996 0.000 0.001 0.002

E(T1)

3.91 1.44 3.18 3.40 3.73

2.02 1.49 2.05

0.84

a

Where S1 has no appreciable intensity, the lowest state with some intensity at this level of theory is also listed. The number of unoccupied ground state canonical molecular orbitals with apparent valence character is given, as well as the number of ground state valence virtual orbitals with pseudoeigenvalues less than zero, at either RHF or B3LYP levels. bGeometry optimization by MP2/SPK-DZP. Results in table increase the basis to SPK-TZP (no diffuse). cThe blue OLED O2S(C6H4-DMAC-DPS)2 of Zhang, Q.; Li, B.; Huang, S.; Nomura, H.; Tanaka, H.; Adachi, C. Nat. Photonics 2014, 8, 326−332. Geometry optimization was by B3LYP/ACCD, and because of the large size of this system, the excited state calculations also use ACCD rather than ACCT. dProtonated, quinoidal, and benzoquinodal forms of phenolphthalein. See Kunimoto, K.-K.; Sugiura, H.; Kato, T.; Senda, H.; Kuwae, A.; Hanai, K. Spectrochim. Acta 2001, A57, 265−271 for practical reasons the ACCD basis was used for RHF geometry optimization. eIndigo dye, see refs 153−155. fFirefly luciferin, from White, E. H.; Capra, F.; McElroy, W. D. J. Am. Chem. Soc. 1961, 83, 2402−2403. g Stereoisomer (S)-arsenicin A, see refs 157−159. hLamivudine, aka 3TC, CAS no. 134678-17-4. iMethyl orange dye, CAS no. 547-58-0. jCyclic adenosine monophosphate, CAS no. 60-92-4. kPentobarbital (barbituate), CAS no. 76-74-4. lAnion number 8 of Perrotta, R. R.; Winter, A. H.; Coldren, W. H.; Falvey, D. E. J. Am. Chem. Soc. 2011, 133, 15553−15558). RHF/ACCD geometry was used.

the neutral and anionic molecules of Table 2, B3LYP yields one or more empty CMOs with apparent valence character. To understand the commonplace B3LYP prediction of valence LCMOs in Table 2, recall the discussion at the end of the previous section, regarding the expected behavior of KS calculations and thus many DFA calculations too. Figure 6 illustrates these remarks for the neutral thymine molecule.152 RHF and the 100% Hartree−Fock exchange PZ8195 functional clearly have diffuse LCMO orbitals, which bear no resemblance whatsoever to the TD-DFT natural orbital which accepts the excited electron in the S1 state (shown at lower right). As one might expect, the ground state’s lowest valence virtual orbital (LVVO) does predict nicely the shape of the excited state’s orbital (the LVVO for RHF, PZ81, and B3LYP all look the same, so only the last is illustrated, bottom row center). The B3LYP LCMO takes on valence character, so like the LVVO, the B3LYP LCMO also correctly predicts the S1 state’s orbital shape. The apparent correctness of the thymine B3LYP LCMO (lower left) is expected behavior for KS calculations in molecules with a low-lying excited state103,104 and also occurs in most DFA applications. This result for thymine is quite typical for most other molecules in Table 2. Note that valence character for empty DFA canonical orbitals does not occur in all molecules: see for example the spatially large TPSS-level LCMO for C2H6 in Figure 1. Of course, ethane lacks low-lying valence excited states.

natural orbitals (both hole and particle) found for each excited state to the occupied CMOs and to either unoccupied CMOs or VVOs. Because the chemical community has adopted time-dependent density functional theory with the B3LYP functional as the most common treatment of excited states, this section presents mainly TD-DFT results (within the Tamm/Dancoff approximation105) using the B3LYP functional and the ACCT basis set. One set of results is presented for the equivalent wave function theory level: TD-HF (which is also known as singly excited configuration interaction, or CIS). The molecules are chosen to have low-lying valence-type excited states, mainly of π* character, and are fairly large in size. Some were designed to have low-lying valence triplet states. One cation was considered, along with several neutrals and anions. Table 2 summarizes the molecules considered: only pentobarbital anion lacks a S1 state with valence character. In keeping with results presented above, one expects the lowest empty Hartree−Fock canonical MOs in a basis set containing diffuse functions to be spatially large except possibly for cations or if there is an exceptionally low-lying valence state (see C60). The second column of Table 2 shows that to be the case for all but one of its neutral and anionic systems: RHF does indeed predict a spatially large nonvalence lowest canonical molecular orbital (LCMO), assessed qualitatively by visualizing them,91 except for indigo. However, in all but one of M

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 6. Thymine orbitals. The LCMOs for RHF or for the pure correlation DFT functional PZ81 shown at the top are clearly similar. The B3LYP calculation’s LCMO + 1 at upper right resembles these two, but the B3LYP LCMO shown at lower left has valence character. This valence LCMO closely resembles the LVVO at bottom center, and both empty B3LYP orbitals are good predictors of the shape of the natural orbital to which the electron is excited in the first singlet state, shown at lower right. The three valence orbitals at the bottom are drawn at 0.05 bohr−3/2, while the three diffuse orbitals at the top are shown at the reduced amplitude value of 0.025 bohr−3/2.

Figure 7. Ground state and first excited state orbitals of indigo, treated by RHF/ACCT and CIS/ACCT, respectively. The first empty CMO closely resembles the first VVO. Higher canonical orbitals are spatially large, as shown for LCMO + 1, while the first six VVOs are various π* shapes like the lowest two shown here. The natural orbitals of the lowest singlet state of indigo show S1 is well described as a HOMO to LVVO (or LCMO) transition. Orbital energy values are in hartree, and occupation numbers are in electron units. The amplitude shown for the spatially large LCMO + 1 is reduced to 0.01 bohr−3/2, while all valence orbitals are drawn at 0.05 bohr−3/2.

The only exception in Table 2 to finding a valence LCMO in B3LYP ground states is the relatively small pentobarbital anion, which has only three nonconjugated CO π occupied orbitals, in addition to being negatively charged. The spatially large LCMO for B3LYP-level pentobarbital does successfully predict the shape of the natural orbital of its S1 electronic state, which is found to possess Rydberg character. Two of the remaining examples of valence states from Table 2 will now be discussed in more detail. Note that natural orbitals of each specific excited state are generated to convert the possibly numerous excitation amplitudes for computations with CMOs into a single orbital from which the electron is removed and another orbital which accepts it. As the next example, consider the blue dye molecule indigo,153 which has very low-lying excited states.154,155 In fact, this means that even for RHF the LCMO has valence character, as shown in Figure 7. By TD-HF (aka CIS) with the ACCT basis set, the S1 state is found to lie at 3.77 eV, not surprisingly less accurate than B3LYP’s 2.62 eV compared to experiment’s 2.28 eV.156 Note that the LVVO is almost identical to the LCMO. The Hartree−Fock LCMO and LVVO are both excellent predictors of the shape of the upper orbital in the S1 state, which differs only slightly from them, at the atom indicated by the red arrow in Figure 7. Not surprisingly, the HOMO of the ground state is also very similar to the orbital from which the electron is excited in S1. Note that the LCMO + 1 lies at an almost identical energy to the LCMO and is now a spatially large orbital. No illustrations are provided in Figure 7 for states above S1, but a short description follows. There are many densely spaced and spatially very large orbitals above the LCMO + 1, whereas the first six VVOs are various kinds of π* orbitals, followed by the first valence σ*. The S2 excited state (4.51 eV by TD-HF/ACCT) involves four natural orbitals, with electron occupancies 1.788, 1.265 (this is the HOMO − 1), 0.735 (this resembles the LVVO), and 0.212 (this resembles the LVVO + 1). The S3 state (4.63 eV by TD-HF/ACCT) is n → π* type, with occupations 1.814, 1.197, 0.802, and 0.183 where once again the more important accepting orbital resembles the Hartree−Fock LCMO/LVVO.

A more interesting example is arsenicin, As4O3(CH2)3, the first molecule containing more than one As to be isolated from a living organism (a sponge found in New Caledonia).157−159 The molecule can be considered to be derived from adamantane, by substituting its CH groups by As and half its CH2 groups by O. It is also related to the laboratory compound As4O6. The experimental spectrum shows four peaks in the visible region, albeit “in the absence of an obvious chromophore”.158 However, the observed transitions were explained by TD-DFT calculations, as due to excitation into σ* orbitals. The B3LYP functional (and others) were found to give a good match to experimental observations. Results in the present work are obtained by TD-DFT using the B3LYP functional and the larger ACCT basis set. The LVVO and LVVO + 1 of arsenicin are illustrated in Figure 8, along with the HOMO and HOMO − 1. The four peaks in the experimental spectrum (according to computed intensities) are identified as S1 (3.82 eV), S4 (4.66 eV), S10 (5.21 eV), and S13 (5.56 eV). S1 is a clear-cut example of a HOMO to LVVO excitation, while the next peak in the spectrum involves four orbitals but with the primary character being HOMO − 1 to LVVO + 1. This state also slightly depopulates the HOMO by placing about 0.2 e− in an orbital resembling the LVVO (these two are not shown to keep the figure simple). Like indigo, this is a molecule where canonical virtuals also do well: the B3LYP LCMO and LCMO + 1 (not illustrated) both closely resemble the LVVO and LVVO + 1. In general, the remaining molecules in Table 2 have both the LCMO and LVVO from B3LYP calculations appearing to be nearly identical, as was explicitly shown for thymine (in Figure 6). Also, as for all three examples shown, the two fractionally occupied natural orbitals of the S1 states are quite well predicted by the shapes of the HOMO and LCMO/LVVO. Typically, the S2 state’s natural orbitals correspond to the excitation from HOMO − 1 to the LCMO/LVVO. Higher excitation energies (S3 on up) do utilize the LVVO + 1 as an upper orbital but frequently involve sufficient energy to N

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

the VVO orbital shapes: if one freezes the occupied orbitals at their SCF shapes but optimizes the three weakly occupied orbitals for HNO, the CAS-SCF energy is −129.9753 hartree. Thus, the missing 10.5% of the MCSCF energy is partitioned as 8.0% due to relaxation (optimization) of the three VVOs, with 2.5% due to relaxation of the six occupied SCF orbitals of HNO. This is reasonable, as the occupied orbitals of the SCF may variationally optimize their usage of the diffuse or polarization parts of the working basis set during the SCF step. In contrast, the singular value decomposition selects VVOs based on their closeness to occupied atomic orbitals, and thus the VVOs have little contribution from the diffuse/ polarization basis set extensions. Two additional examples, similar to those given elsewhere,5 are diborane and dioxirane. Full valence CI within a valence orbital space containing all occupied MOs and all VVOs recovers 86.3% and 87.3% of the optimized full valence CASSCF energies for diborane (14 e− in 12ϕ) or dioxirane (18 e− in 14ϕ), respectively. It is reasonable to consider the union of occupied SCF orbitals with VVOs, followed by a full CI in those orbitals, as “a poor man’s MCSCF”. In cases where the molecule is too large to permit full valence MCSCF calculations, VVOs can still guide the calculation in two ways. For o-benzyne160,161 a reasonable active space is obvious on chemical grounds: the three π and three π* orbitals perpendicular to the ring plus the in-plane π and π*. Here, the symmetry constrained localization, applied separately in the occupied and VVO spaces, generates excellent starting orbitals. The initial MCSCF iteration using these localized starting orbitals already recovers 83% of the correlation energy contained in this eight electron, eight orbital active space. Of course there are also cases where the active space is much more difficult to anticipate from chemical intuition. One can simply generate and visualize VVOs to consider the low energy VVOs as active orbital candidates. If the active space is still not obvious, one can perform a limited CI calculation, such as singles and doubles only, f rom all occupied orbitals into all VVOs, which is feasible for rather large molecules. The resulting natural orbitals of this CI-SD will reveal those orbitals whose occupation deviates significantly from 2.0 or 0.0: these fractionally occupied orbitals constitute an automatically generated set of “active orbitals” for the system’s multiconfigurational wave function. One can repeat this process at several geometries typical of the PES being examined to ensure an entire set of important active orbitals is identified. Of course, full valence MCSCF is not the only kind of wave function that involves the occupation of antibonding levels. The GVB perfect pairing (PP) wave function is another example. This wave function contains pairs of orbitals (geminals) in which a bonding orbital is doubly excited into its antibonding counterpart. This is precisely the type of orbital pairs found when localizing the occupied and the VVO spaces separately. Consequently, the glycine example proposed as a GVB-PP(10) convergence test162 is readily convergent from VVOs, after preparing ten perfectly paired starting geminals by matching each of the ten localized closed shell bonding orbitals with its corresponding ten antibonding orbitals, obtained by localizing the VVOs. The RHF energy of glycine is −282.8373 hartree, the initial iteration with the ten starting geminals is −282.9475 hartree, and convergence to the optimal GVB-PP(10) energy of −283.0194 hartree is smooth. The initial geminals generated from VVOs thus provided some 60% of the GVB-PP(10) correlation energy.

Figure 8. Comparison of the ground state’s B3LYP/ACCT occupied CMOs and empty VVOs in As4O3(CH2)3 to the natural orbitals of the TDDFT 1B S1 and 1A S4 valence excited states. Orbital energies are in hartree for the former, while the latter are labeled by their electron occupation numbers. The view is along the C2 axis of the molecule.

scramble its shape with the LVVO and to have mixed occupancies. Typically, the higher states have four natural orbitals with mixed occupancies and resemble combinations of the ground state’s HOMO − 1, HOMO, LVVO, and LVVO + 1. Arsenicin is a notable counterexample, as its S4 state can be qualitatively considered as a simple HOMO − 1 to LVVO + 1 single excitation. Clearly the chemical community’s oft-made choice of B3LYP TD-DFT calculations for large molecule excited state predictions is related to the circumstance that the ground state’s LCMO will often have valence character, when the S1 state is a valence state, whereas RHF often does not. The LVVO from either Hartree−Fock or DFA calculations is capable of predicting the lowest valence state, even in cases like pentobarbital anion where S1 is not a valence state, but higher states are. III.F. VVOs as MCSCF Starting Orbitals. In addition to the conceptual value of VVOs, discussed above, the VVOs are also found to have practical use in designing and carrying out multireference calculations. As shown elsewhere,5 the occupied SCF orbitals plus the VVOs form a very good approximation to the full valence space MCSCF wave function. Performing a full CI calculation within the valence orbital space recovers around 80−90% of the correlation energy found after CAS-SCF orbital optimization. This is impressive, considering the O(N3) effort in preparing VVOs, compared to the O(N5) integral transformation work in running a MCSCF program. As an example, the SCF energy of HNO is −129.8516 hartree, the full valence CI using its closed shell VVOs is −129.9652 hartree, and a fully relaxed CASSCF is −129.9784 hartree (energies are from the aug-cc-pVQZ basis set). The VVOs thus provide 89.5% of the near-degeneracy correlation recovered by the CAS-SCF. Most of the defect is in O

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

■

IV. SUMMARY This paper has illustrated the previously known result (at least among theoretical chemists) that for most neutral molecule SCF calculations, the lowest canonical MO changes shape with an increasing basis set, becoming spatially large and moving to an energy near zero. The LCMO therefore does not correspond well to the LUMO concept in chemistry. The circumstances in which the LCMO is a better approximation to a valence antibonding level are as follows: if the system has a positive charge or if there are low-lying excited states or if the theoretical computation involves small basis sets or if a DFA functional is used. The converse of that list is the circumstances for which the LCMO is expected to be less useful. In contrast, the lowest VVO does converge with respect to basis set, to a reasonable energy pseudoeigenvalue, independent of the molecular charge, or the energy needed to access its first valence excited state. Consequently the VVOs allow the construction of MO diagrams including all valence orbitals, in qualitatively reasonable positions. The LVVO is also successful in predicting the nature of the first excited valence state, as illustrated for a number of molecules which have a valence state at low excitation energy. Finally, VVOs provide excellent starting orbitals for multireference computations. This paper has not described the uses for orbitals localized onto specific atoms, which are easily obtained by allowing a localization to mix together the occupied CMOs and the VVO spaces freely. The use of oriented atom localized orbitals in bonding analyses (bond orders and kinetic bond orders) and in charge population analyses can be found elsewhere.5,6

■

Article

REFERENCES

(1) Fukui, K.; Yonezawa, T.; Shingu, H. A Molecular Orbital Theory of Reactivity in Aromatic Hydrocarbons. J. Chem. Phys. 1952, 20, 722− 725. (2) Woodward, R. B.; Hoffmann, R. The Conservation of Orbital Symmetry; Verlag Chemie GmbH: Weinheim, Germany, 1971. (3) Lu, W. C.; Wang, C. Z.; Schmidt, M. W.; Bytautas, L.; Ho, K. M.; Ruedenberg, K. Molecule Intrinsic Minimal Basis Sets. I. Exact Resolution of ab initio Optimized Molecular Orbitals in Terms of Deformed Atomic Minimal-basis Orbitals. J. Chem. Phys. 2004, 120, 2629−2637. (4) Lu, W. C.; Wang, C. Z.; Schmidt, M. W.; Bytautas, L.; Ho, K. M.; Ruedenberg, K. Molecule Intrinsic Minimal Basis Sets. II. Bonding Analyses for Si4H6 and Si2 to Si10. J. Chem. Phys. 2004, 120, 2638− 2651. (5) West, A. C.; Schmidt, M. W.; Gordon, M. S.; Ruedenberg, K. A Comprehensive Analysis of Molecule Intrinsic Quasi-atomic, Bonding, and Correlating Orbitals. I. Hartree-Fock Wave Functions. J. Chem. Phys. 2013, 139, 234107. (6) West, A. C.; Schmidt, M. W.; Gordon, M. S.; Ruedenberg, K. A Comprehensive Analysis in Terms of Molecule Intrinsic Quasi-atomic Orbitals. II. Strongly Correlated MCSCF Wave Functions. J. Phys. Chem. A, submitted.201510.1021/acs.jpca.5b03399 (7) Gutowski, M.; Jordan, K. D.; Skurski, P. Electronic Structure of Dipole-bound Anions. J. Phys. Chem. A 1998, 102, 2624−2633. (8) Jordan, K. D.; Johnson, M. A. Downsizing the Hydrated Electron’s Lair. Science 2010, 329, 42−43. (9) Illas, F.; Merchan, M.; Pelissier, M.; Malrieu, J.-P. Inexpensive Determination of Valence Virtual MOs for CI Calculations. Chem. Phys. 1986, 107, 361−380. (10) Foster, J. M.; Boys, S. F. Canonical Configuration Interaction Procedure. Rev. Mod. Phys. 1960, 32, 300−302. (11) Huzinaga, S.; Arnau, C. Virtual Orbitals in Hartree-Fock Theory. J. Chem. Phys. 1971, 54, 1948−1951. (12) Beebe, N. H. F. Modification of Virtual Orbitals. Int. J. Quantum Chem. 1979, 15, 589−600. (13) Hunt, W. J.; Goddard, W. A. Excited States of H2O using Improved Virtual Orbitals. Chem. Phys. Lett. 1969, 3, 414−418. (14) Bauschlicher, C. W. The Construction of Modified Virtual Orbitals (MVOs) which are Suited for Configuration Interaction Calculations. J. Chem. Phys. 1980, 72, 880−885. (15) Mogensen, B. J.; Rettrup, S. Average Virtual Orbitals in Configuration Interaction Studies with Application to the low-lying Singlet States of the Carbon Monoxide and Acetone Molecules. Int. J. Quantum Chem. 1992, 44, 1045−1056. (16) Davidson, E. R. Selection of the Proper Canonical RoothaanHartree-Fock Orbitals for Particular Applications. I. Theory. J. Chem. Phys. 1972, 57, 1999−2005. (17) Cooper, I. L.; Pounder, C. N. M. A Simple Procedure for the Improvement of CI Convergence in Molecular Systems. J. Chem. Phys. 1979, 71, 957−960. (18) Das, G.; Wahl, A. C. New Techniques for the Computation of Multiconfiguration Self-Consistent Field Wavefunctions. J. Chem. Phys. 1972, 56, 1769−1775. (19) Whitten, J. L. Remarks on the Description of Excited Electronic States by Configuration Interaction Theory and a Study of the 1(π→ π*) State of H2CO. J. Chem. Phys. 1972, 56, 5458−5466. (20) Luken, W. L.; Seiders, B. A. B. Interaction-optimized Virtual Orbitals. I. External Double Excitations. Chem. Phys. 1985, 92, 235− 246. (21) Feller, D.; Davidson, E. R. An Approximation to Frozen Natural Orbitals through the Use of the Hartree-Fock Exchange Potential. J. Chem. Phys. 1981, 74, 3977−3979. (22) Wasilewski, J. Modified Virtual Orbitals in Limited CI Calculations. Int. J. Quantum Chem. 1991, 39, 649−656. (23) Jensen, H. J. A.; Jorgensen, P.; Agren, H.; Olsen, J. Second-order M?ller-Plesset Perturbation Theory as a Configuration and Orbital Generator in Multiconfiguration Self-consistent Field Calculations. J. Chem. Phys. 1988, 88, 3834−3838.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Telephone: 1-515-2949796. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Prof. Klaus Ruedenberg for his decades-long conferral to them of chemical bonding theory, as well as for recognizing the VVO process is an instance of the singular value decomposition and for the need to preorthogonalize. The authors thank Dr. Aaron West for his assistance in generating VVOs after open shell and multiconfigurational SCF calculations. The authors thank Prof. Mariusz Klobukowski for providing full details about the well-tempered basis sets. M.W.S. thanks Prof. Keiko Takano and her students for their most gracious hospitality during his visit to Ochanomizu University in the fall of 2007, when the closed shell QUAMBO/VVO program was finished. M.W.S. thanks Prof. Hiromi Nakai for mentioning the curious case of 4-fold degeneracy of both frontier orbitals in B12H12−2. The authors thank Klaus Ruedenberg, Aaron West, and Joe Brom for their comments on the manuscript. E.A.H. and T.L.W.’s contribution was supported by the National Science Foundation under Grant No. OISE-0730114 for the Partnerships in International Research and Education (PIRE). M.W.S.’s contribution was supported by the National Science Foundation under Grant CHE-1147446. P

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (24) Caballol, R.; Malrieu, J.-P. Improved Non-valence Virtual Orbitals for CI Calculations. Chem. Phys. 1990, 140, 7−18. (25) Bender, C. F.; Davidson, E. R. A Natural Orbital Based Energy Calculation for Helium Hydride and Lithium Hydride. J. Phys. Chem. 1966, 70, 2675−2685. (26) Jafri, J. A.; Whitten, J. L. Iterative Natural Orbitals for Configuration Interaction Using Perturbation Theory. Theoret. Chim. Acta 1977, 44, 305−313. (27) Abrams, M. L.; Sherrill, C. D. Natural Orbitals as Substitutes for Optimized Orbitals in Complete Active Space Wavefunctions. Chem. Phys. Lett. 2004, 395, 227−232. (28) Krausbeck, F.; Mendive-Tapia, D.; Thom, A. J. W.; Bearpark, M. J. Choosing RASSCF Orbital Active Spaces for Multiple Electronic States. Comput. Theor. Chem. 2014, 1040−1041, 14−19. (29) Yamaguchi, K. The Electronic Structure of Biradicals in the Unrestricted Hartree-Fock Appoximation. Chem. Phys. Lett. 1975, 33, 330−335. (30) Pulay, P.; Hamilton, T. P. UHF Natural Orbitals for Defining and Starting MC-SCF Calculations. J. Chem. Phys. 1988, 88, 4926− 4933. (31) Iwata, S. Valence Type Vacant Orbitals for Configuration Interaction Calculations. Chem. Phys. Lett. 1981, 83, 134−138. (32) Chambaud, G.; Gerard-Ain, M.; Kassab, E.; Levy, B.; Pernot, P. Valence-Bond Calculations with Polarized Atomic Orbitals. Chem. Phys. 1984, 90, 271−289. (33) Ruedenberg, K.; Schmidt, M. W.; Gilbert, M. M. Are Atoms Intrinsic to Molecular Electronic Wavefunctions? II. Analysis of FORS Orbitals. Chem. Phys. 1982, 71, 51−64. (34) Stewart, G. W. On the Early History of the Singular Value Decomposition. SIAM Rev. 1993, 35, 551−566. (35) See the appendix of Amos, A. T.; Hall, G. G. Single Determinant Wave Functions. Proc. R. Soc. London, Ser. A 1961, 263, 483−493. (36) King, H. F.; Stanton, R. E.; Kim, H.; Wyatt, R. E.; Parr, R. G. Corresponding Orbitals and the Nonorthogonality Problem in Molecular Quantum Mechanics. J. Chem. Phys. 1967, 47, 1936−1941. (37) See appendix A of Sax, A. F. Localization of Molecular Orbitals on Fragments. J. Comput. Chem. 2012, 33, 1495−1510. (38) See the appendix of ref 5. (39) Lee, M. S.; Head-Gordon, M. Extracting Polarized Atomic Orbitals from Molecular Orbital Calculations. Int. J. Quantum Chem. 2000, 76, 169−184. (40) Subotnik, J. E.; Dutoi, A. D.; Head-Gordon, M. Fast Localized Orthonormal Virtual Orbitals which Depend Smoothly on Nuclear Coordinates. J. Chem. Phys. 2005, 123, 114108. (41) Auer, A. A.; Noojien, M. Dynamically Screened Local Correlation Method Using Enveloping Localized Orbitals. J. Chem. Phys. 2006, 125, 024104. (42) Laikov, D. M. Intrinsic Minimal Atomic Basis Representation of Molecular Electronic Wavefunctions. Int. J. Quantum Chem. 2011, 111, 2851−2867. (43) Szczepanik, D.; Mrozek, J. Minimal Set of Molecule-adapted Atomic Orbitals from Maximum Overlap Criterion. J. Math. Chem. 2013, 51, 2687−2698. (44) Knizia, G. Intrinsic Atomic Orbitals: an Unbiased Bridge between Quantum Theory and Chemical Concepts. J. Chem. Theory Comput. 2013, 9, 4834−4843. (45) Janowski, T. Near Equivalence of Intrinsic Atomic Orbitals and Quasiatomic Orbitals. J. Chem. Theory Comput. 2014, 10, 3085−3091. (46) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Intermolecular Interactions from a Natural Bond Orbital, Donor-Acceptor Viewpoint. Chem. Rev. 1988, 88, 899−926. (47) Nemukhin, A. V.; Weinhold, F. Natural Bond Orbitals in Multiconfigurational Expansions: Local Treatment of Electron Correlation in Molecules. J. Chem. Phys. 1992, 97, 1095−1108. (48) Weinhold, F.; Landis, C. Valency and Bonding; Cambridge University Press: Cambridge, U.K., 2005. (49) Mayer, I. Atomic Orbitals from Molecular Wave Functions: the Effective Minimal Basis. J. Phys. Chem. 1996, 100, 6249−6257.

(50) Zubarev, D. Yu.; Boldyrev, A. I. Developing Paradigms of Chemical Bonding: Adaptive Natural Density Partitioning. Phys. Chem. Chem. Phys. 2008, 10, 5207−5217. (51) Bader, R. F. W. Atoms in Molecules, a Quantum Theory; Clarendon Press: Oxford, U.K., 1994. (52) Pulay, P. Localizability of Dynamic Electron Correlation. Chem. Phys. Lett. 1983, 100, 151−154. (53) See section VI of ref 5 and also West, A. C. Weighted Orthogonalization of Atomic Orbitals: a Stable Alternative to the Carlson-Keller Method. Comput. Theor. Chem. 2014, 1045, 73−77. (54) Jansik, B.; Host, S.; Kristensen, K.; Jorgensen, P. Local Orbitals by Minimizing Powers of the Orbital Variance. J. Chem. Phys. 2011, 134, 194104. (55) Hoyvik, I.-M.; Jansik, B.; Jorgensen, P. Pipek-Mezey Localization of Occupied and Virtual Orbitals. J. Comput. Chem. 2013, 34, 1456−1462. (56) Lu, W. C.; Wang, C. Z.; Ruedenberg, K.; Ho, K. M. Transferability of the Slater-Koster Tight-Binding Scheme from a Environment-Dependent Minimal-Basis Perspective. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 205123. (57) Dunnington, B. D.; Schmidt, J. R. Generalization of Natural Bond Orbital Analysis to Periodic Systems: Applications to Solids and Surfaces via Plane-wave Density Functional Theory. J. Chem. Theory Comput. 2012, 8, 1902−1911. (58) Galeev, T. R.; Dunnington, B. D.; Schmidt, J. R.; Boldyrev, A. I. Solid State Adaptive Natural Density Partitioning: a Tool for Deciphering Multicenter Bonding in Periodic Systems. Phys. Chem. Chem. Phys. 2013, 15, 5022−5029. (59) Glaesemann, K. R.; Schmidt, M. W. On the Ordering of Orbital Energies in High-spin ROHF. J. Phys. Chem. A 2010, 114, 8772−8777. (60) Rappe, A. K.; Smedley, T. A.; Goddard, W. A. Flexible d Basis Sets for Scandium through Copper. J. Phys. Chem. 1981, 85, 2607− 2611. (61) Moore, C. E. Atomic Energy Levels; National Bureau of Standards: Washington DC, 1949 and 1952; Volumes I and II. (62) Ivanic, J.; Ruedenberg, K. Identification of Deadwood in Configuration Spaces through General Direct Configuration Interaction. Theor. Chem. Acc. 2001, 106, 339−351. (63) Ivanic, J. Direct Configuration Interaction and Multiconfigurational Self-Consistent-Field method for Multiple Active Spaces with Variable Occupations. I. Method and II. Application to oxoMn(salen) and N2O4. J. Chem. Phys. 2003, 119, 9364−9376 and 9377−9385. (64) Wang, S.-G.; Schwarz, W. H. E. Icon of Chemistry: the Periodic System of Chemical Elements in the New Century. Angew. Chem., Int. Ed. 2009, 48, 3404−3415. (65) Schwarz, W. H. E.; Rich, R. L. Theoretical Basis and Correct Explanation of the Periodic System: Review and Update. J. Chem. Educ. 2010, 87, 435−443. (66) Schwarz, W. H. E. The Full Story of the Electron Configurations of the Transition Elements. J. Chem. Educ. 2010, 87, 444−448. (67) Schmidt, M. W.; Ruedenberg, K. Effective Convergence to Complete Orbital Bases and to the Atomic Hartree-Fock Limit though Systematic Sequences of Gaussian Primitives. J. Chem. Phys. 1979, 71, 3951−3962. (68) Huzinaga, S.; Miguel, B. A Comparison of the Geometrical Sequence Formula and the Well-tempered Formulas for Generating GTO Basis Orbital Exponents. Chem. Phys. Lett. 1990, 175, 289−291. (69) Huzinaga, S.; Klobukowski, M. Well-tempered Gaussian Basis Sets for the Calculation of Matrix Hartree-Fock Wavefunctions. Chem. Phys. Lett. 1993, 212, 260−264. (70) Huzinaga, S.; Miguel, B.; Klobukowski, M. Well-Tempered Basis Set parameters (He-Rn); 2011, unpublished (private communication from M.K). (71) The WTBS parameters do not expand the valence 5s orbital of Pd well, since they were optimized for the 4d10 ground state. The WTBS-style basis used here for Pd is the same size as Rh and Ag, which bracket Pd in the periodic table and are the average of those two elements’ WTBS parameters: Pd uses α = 0.016 382 347 5, β = 1.897 947 95, γ = 1.612 041 55, and δ = 7.110 160 45. Q

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

(96) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868; Phys. Rev. Lett. 1997, 78, 1396. (97) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: the PBE0Model. J. Chem. Phys. 1999, 110, 6158−6170. (98) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5642. (99) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio Calculation of Vibrational Absorption and Circular Dichroism Spectra using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (100) Becke, A. D. Density-functional Thermochemistry. V. Systematic Optimization of Exchange-correlation Functionals. J. Chem. Phys. 1997, 107, 8554−8560. (101) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical MetaGeneralized Gradient Approximations Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. (102) Perdew, J. P.; Tao, T.; Staroverov, V. N.; Scuseria, G. E. MetaGeneralized Gradient Approximation: Explanation of a Realistic Nonempirical Density Functional. J. Chem. Phys. 2004, 120, 6898− 6911. (103) Baerends, E. J.; Gritsenko, O. V.; van Meer, R. The KohnSham Gap, the Fundamental Gap and the Optical Gap: the Physical Meaning of Occupied and Virtual Kohn-Sham Orbital Energies. Phys. Chem. Chem. Phys. 2013, 15, 16408−16425. (104) van Meer, R.; Gritsenko, O. V.; Baerends, E. J. Physical Meaning of Virtual Kohn-Sham Orbitals and Orbital Energies: an Ideal Basis for the Description of Molecular Excitations. J. Chem. Theory Comput. 2014, 10, 4432−4441. (105) Hirata, S.; Head-Gordon, M. Time-dependent Density Functional Theory within the Tamm-Dancoff Approximation. Chem. Phys. Lett. 1999, 314, 291−299. (106) Dreuw, A.; Head-Gordon, M. Single-Reference ab initio Methods for the Calculation of Excited State of Large Molecules. Chem. Rev. 2005, 105, 4009−4037. (107) Elliott, P.; Furche, F.; Burke, K. Excited States from TimeDependent Density Functional Theory. Rev. Comput. Chem. 2009, 26, 91−166. (108) Koopmans, T. Uber die Zuordnung von Wellenfunktionen und Eigenwerten zu den Einzelnen Elektronen eines Atoms. Physica 1934, 1, 104−113. (109) Heinrich, N.; Koch, W.; Frenking, G. On the Use of Koopmans’ Theorem to Estimate Negative Electron Affinities. Chem. Phys. Lett. 1986, 124, 20−25. (110) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theoret. Chim. Acta 1973, 28, 213−222. (111) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent Molecular Orbital Methods 25. Supplementary Functions for Gaussian Basis Sets. J. Chem. Phys. 1984, 80, 3265−3269. (112) Rassolov, V. A.; Pople, J. A.; Ratner, M. A.; Windus, T. L. 631G* Basis Sets for Atoms K through Zn. J. Chem. Phys. 1998, 109, 1223−1229. (113) Scholz, F.; Himmel, D.; Heinemann, F. W.; Schleyer, P. v. R.; Meyer, K.; Krossing, I. Crystal Structure Determination of the Nonclassical 2-Norbornyl Cation. Science 2013, 341, 62−64. (114) The ratio of the two smallest exponents present in the 631G(f,d) basis is used to extrapolate to diffuse augmentation exponents sp = 0.014 and d = 0.128 for Mn, which were combined with the standard O diffuse sp = 0.0845 GTO shell. (115) Suenram, R. D.; Lovas, F. J. Dioxirane. Its Synthesis, Microwave Spectrum, Structure, and Dipole Moment. J. Am. Chem. Soc. 1978, 100, 5117−5122. (116) Song, J. K.; Lee, N. K.; Kim, S. K. Photoelectron Spectrum of Pyrazine Anion Clusters. J. Chem. Phys. 2002, 117, 1589−1594. (117) Garza, J.; Nichols, J. A.; Dixon, D. A. The Optimized Effective Potential and the Self-interaction Correction in Density Functional

(72) Barysz, M.; Sadlej, A. J. Infinite-Order Two-Component Theory for Relativistic Quantum Chemistry. J. Chem. Phys. 2002, 116, 2696− 2704. (73) Douglas, M.; Kroll, N. M. Quantum Electrodynamical Corrections to the Fine Structure of Helium. Ann. Phys. 1974, 82, 89−155. (74) Hess, B. A. Relativistic Electronic-structure Calculations Employing a Two-component No-pair Formalism with External-field Projection Operators. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 33, 3742−3748. (75) Nakajima, T.; Yanai, T.; Hirao, K. Relativistic Electronic Structure Theory. J. Comput. Chem. 2002, 23, 847−860. (76) Pyykko, P. Relativistic Effects in Structural Chemistry. Chem. Rev. 1988, 88, 563−594. (77) Pyykko, P. Relativistic Effects in Chemistry: More Common than You Thought. Annu. Rev. Phys. Chem. 2012, 63, 45−64. (78) Huzinaga, S. 1994 Polanyi Award Lecture Concept of Active Electrons in Chemistry. Can. J. Chem. 1995, 73, 619−628. (79) Klobukowski, M.; Huzinaga, S.; Sakai, Y. Model Core Potentials: Theory and Applications. In Computational Chemistry: Reviews of Current Trends; Leszczynski, J., Ed.; World Scientific: Singapore, 1999; Vol. 3, pp 49−74. (80) Krauss, M.; Stevens, W. J. Effective Potentials in Molecular Quantum Chemistry. Annu. Rev. Phys. Chem. 1984, 35, 357−385. (81) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (82) Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: GAMESS a Decade Later. In Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E.. Eds.; Elsevier: Amsterdam, 2005; pp 1167−1189. (83) GAMESS. http://www.msg.chem.iastate.edu/GAMESS/ GAMESS.html. (84) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron Affinities of the First-row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (85) Dunning, T. H.; Peterson, K. A.; Wilson, A. K. Gaussian Basis Sets for Use in Correlated Molecular Calculations. X. The Atoms Aluminum through Argon Revisited. J. Chem. Phys. 2001, 114, 9244− 9253. (86) Balabanov, N. B.; Peterson, K. A. Systematically Converging Basis Sets for Transition Metals. I. All-electron Correlation Consistent Basis sets for the 3d Elements Sc-Zn. J. Chem. Phys. 2005, 123, 064107. (87) Peterson Research Group. http://tyr0.chem.wsu.edu/~kipeters/ basis-bib.html. (88) Noro, T.; Sekiya, M.; Koga, T. Contracted Polarization Functions for the Atoms Helium through Neon. Theor. Chem. Acc. 1997, 98, 25−32. (89) Koga, T.; Yamamoto, S.; Shimazaki, T.; Tatewaki, H. Contracted Gaussian-type Basis Functions Revisited. Theor. Chem. Acc. 2002, 108, 41−45. (90) Segmented Gaussian Basis Set. http://setani.sci.hokudai.ac.jp/ sapporo/Welcome.do. (91) Bode, B. M.; Gordon, M. S. MacMolPlt: a Graphical User Interface for GAMESS. J. Mol. Graphics Modell. 1998, 16, 133−138. (92) http://brettbode.github.io/wxmacmolplt/. (93) Schmidt, M. W.; Gordon, M. S. The Construction and Interpretation of MCSCF Wavefunctions. Annu. Rev. Phys. Chem. 1998, 49, 233−266. (94) Aikens, C. M.; Webb, S. P.; Bell, R. L.; Fletcher, G. D.; Schmidt, M. W.; Gordon, M. S. A Derivation of the Frozen-orbital Unrestricted Open Shell and Restricted Closed Shell MP2 Analytic Gradient Expressions. Theor. Chem. Acc. 2003, 110, 233−253. (95) Perdew, J. P.; Zunger, A. Self-interaction Correction to Densityfunctional Approximations for Many-electron Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048−5079. R

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Theory: Application to Molecules. J. Chem. Phys. 2000, 112, 7880− 7890. (118) Edmiston, C.; Ruedenberg, K. Localized Atomic and Molecular Orbitals. Rev. Mod. Phys. 1963, 35, 457−465. (119) Raffenetti, R. C.; Ruedenberg, K.; Janssen, C. L.; Schaefer, H. F. Efficient Use of Jacobi Rotations for Orbital Optimization and Localization. Theoret. Chim. Acta 1993, 86, 149−165. (120) Boys, S. F. Localized Orbitals and Localized Adjustment Functions. In Quantum Science of Atoms, Molecules, and Solids; Lowdin, P.-O., Ed.; Academic Press: New York, 1966; pp 253−262. (121) The equations of ref 120 have never been implemented, as “Boys” localization programs actually follow the suggestion of a similar criterion given in ref 118. (122) Pipek, J.; Mezey, P. Z. A Fast Intrinsic Localization Procedure Applicable for ab initio and Semiempirical Linear Combination of Atomic Orbital Wavefunctions. J. Chem. Phys. 1989, 90, 4916−4926. (123) Musher, J. I. The Chemistry of Hypervalent Molecules. Angew. Chem., Int. Ed. Engl. 1969, 8, 54−68. (124) Gillespie, R. J.; Nyholm, R. S. Inorganic Stereochemistry. Q. Rev., Chem. Soc. 1957, 11, 339−380. (125) Schmidt, M. W.; Gordon, M. S. On the Observability of Cubic P8. Inorg. Chem. 1985, 24, 4503−4506. (126) Raghavachari, K.; Haddon, R. C.; Binkley, J. S. Small Elemental Clusters: Theoretical Study of P, P2, P4, and P8. Chem. Phys. Lett. 1985, 122, 219−223. (127) Jerabek, P.; Frenking, G. Comparative Bonding Analysis of N2 and P2 versus Tetrahedral N4 and P4. Theor. Chem. Acc. 2014, 133, 1447−1456. (128) Werner, H. At least 60 Years of Ferrocene: the Discovery and Rediscovery of the Sandwich Complexes. Angew. Chem., Int. Ed. 2012, 51, 6052−6058. (129) Haaland, A.; Nilsson, J. E. The Determination of Barriers to Internal Rotation by Means of Electron Diffraction. Ferrocene and Ruthenocene. Acta Chem. Scand. 1968, 22, 2653−2670. (130) Lancaster, K. M.; Finkelstein, K. D.; DeBeer, S. Kβ X-ray Emission Spectroscopy Offers Unique Chemical Bonding Insights: Revisiting the Electronic Structure of Ferrocene. Inorg. Chem. 2011, 50, 6767−6774. (131) Atkins, A. J.; Bauer, M.; Jacob, C. R. The Chemical Sensitivity of X-ray Spectroscopy: High Energy Resolution XANES versus X-ray Emission Spectroscopy of Substituted Ferrocenes. Phys. Chem. Chem. Phys. 2013, 15, 8095−8105. (132) For ferrocene, the closed shell SCF occupied active CMO energies are a1g = −0.569, a1g = −0.508, a2u = −0.500, e2g = −0.427, e1u = −0.340, e1g = −0.338 with empty VVO energies e2u = +0.298, e2g = +0.302, e1g = +0.309, a1g = +0.366 hartree. (133) For technical reasons, the DFA calculations on Fe and Cp· used an open shell DFA program that is unable to rigorously enforce orbital degeneracy: this has been introduced into Figure 3 by averaging over the nearly degenerate orbitals in Fe and Cp·. Additionally, open shell energies are not uniquely defined (ref 59 and Plakhutin, B.; Davidson, E. R. Canonical Form of the Hartree-Fock Orbitals in Open-shell Systems. J. Chem. Phys. 2014, 140, 014102). When tested, this procedure gave eigenvalues very close to those from the open shell SCF level program that can enforce orbital degeneracy. No such concerns apply to the energy levels for ferrocene, which is a closed shell system containing only filled or empty degenerate orbitals. (134) Rayon, V. M.; Frenking, G. Bis(benzene)chromium is a δbonded Molecule and Ferrocene is a π-bonded Molecule. Organometallics 2003, 22, 3304−3308. (135) Dunitz, J. D.; Orgel, L. E. Bis-cyclopentadienyl Iron: a Molecular Sandwich. Nature 1953, 171, 121−122. (136) Haddon, R. C. Electronic Structure, Conductivity, and Superconductivity of Alkali Metal Doped C60. Acc. Chem. Res. 1992, 25, 127−133. (137) Rioux, F. Quantum Mechanics, Group Theory, and C60. J. Chem. Educ. 1994, 71, 464−465.

(138) Fukuda, R.; Ehara, M. Electronic Excitations of C60 Fullerene Calculated using the ab initio Cluster Expansion Method. J. Chem. Phys. 2012, 137, 134304. (139) Longuet-Higgins, H. C.; de V. Roberts, M. The Electronic Structure of an Icosahedron of Boron Atoms. Proc. R. Soc. London, Ser. A 1955, A230, 110−119. (140) Feller, D. F.; Schmidt, M. W.; Ruedenberg, K. Concerted Dihydrogen Exchange between Ethane and Ethylene. SCF and FORS Calculations of the Barrier. J. Am. Chem. Soc. 1982, 104, 960−967. (141) Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. Approximation of the Exchange-Correlation Kohn-Sham Potential with a Statistical Average of Different Orbital Model Potentials. Chem. Phys. Lett. 1999, 302, 199−207. (142) de Vries, J.; Steger, H.; Kamke, B.; Menzel, C.; Weisser, B.; Kamke, W.; Hertel, I. V. Single-photon Ionization of C60 and C70 Fullerene with Synchrotron Radiation: Determination of the Ionization Potential of C60. Chem. Phys. Lett. 1992, 188, 159−162. (143) Wang, X.-B.; Ding, C.-F.; Wang, L.-S. High Resolution Photoelectron Spectroscopy of C60−. J. Chem. Phys. 1999, 110, 8217− 8220. (144) Zhou, W.-Y.; Xie, S.-S.; Qian, S.-F.; Wang, G.; Qian, L.-X. Photothermal Deflection Specta of Solid C60. J. Phys.: Condens. Matter 1996, 8, 5793−5800. (145) Garza, J.; Nichols, J. A.; Dixon, D. A. The Hartree product and the Description of Local and Global Quanties in Atomic Systems: A Study with Kohn-Sham Theory. J. Chem. Phys. 2000, 112, 1150−1157. (146) Garza, J.; Nichols, J. A.; Dixon, D. A. The Role of the Localmultiplicative Kohn-Sham Potential on the Description of Occupied and Unoccupied Orbitals. J. Chem. Phys. 2000, 113, 6029−6034. (147) Pederson, M. R.; Ruzsinszky, A.; Perdew, J. R. Communication: Self-Interaction Correction with Unitary Invariance in Density Functional Theory. J. Chem. Phys. 2014, 140, 121103. (148) Zhang, G.; Musgrave, C. B. Comparison of DFT Methods for Molecular Orbital Eigenvalue Calculations. J. Phys. Chem. A 2007, 111, 1554−1561. (149) Chong, D. P.; Gritsenko, O. V.; Baerends, E. J. Interpretation of the Kohn-Sham Orbital Energies as Approximate Vertical Ionization Potentials. J. Chem. Phys. 2002, 116, 1760−1771. (150) Yang, W.; Cohen, A. J.; Mori-Sanchez, P. Derivative Discontinuity, Bandgap, and Lowest Unoccupied Molecular Orbital in Density Functional Theory. J. Chem. Phys. 2012, 136, 204111. (151) Tsuneda, T.; Song, J.-W.; Suzuki, S.; Hirao, K. On Koopmans’ Theorem in Density Functional Theory. J. Chem. Phys. 2010, 133, 174101. (152) For thymine, the B3LYP LCMO and LVVO both appear slightly below zero, and the B3LYP HOMO energy is elevated (akin to DFA results given for C60). The HOMO energies with the ACCT basis are RHF = −0.356, PZ81 = −0.406, B3LYP = −0.253. Unoccupied orbital energies are contained in Figure 5. (153) von Baeyer, A. Zur Geschichte der Indigo-Synthese. Ber. Dtsch. Chem. Ges. 1900, 33, LI−LXX. (154) Jacquemin, D.; Preat, J.; Wathelet, V.; Perpete, E. A. Substitution and Chemical Environment Effects on the Absorption Spectrum of Indigo. J. Chem. Phys. 2006, 124, 074104. (155) Ngan, V. T.; Gopakumar, G.; Hue, T. T.; Nguyen, M. T. The Triplet State of Indigo: Electronic Structure Calculations. Chem. Phys. Lett. 2007, 449, 11−17. (156) Luttke, W.; Hermann, H.; Klessinger, M. Theoretically and Experimentally Determined Properties of the Fundamental Indigo Chromophore. Angew. Chem., Int. Ed. Engl. 1966, 5, 598−599. (157) Lu, D.; Rae, A. D.; Salem, G.; Weir, M. L.; Willis, A. C.; Wild, S. B. Arsenicin A, a Natural Polyarsenical: Synthesis and Crustal Structure. Organometallics 2010, 29, 32−33. (158) Arulmozhiraja, S.; Coote, M. L.; Lu, D.; Salem, G.; Wild, S. B. Origin of the Unusual Ultraviolet Absorption of Arsenicin A. J. Phys. Chem. A 2011, 115, 4530−4534. (159) Lu, D.; Coote, M. L.; Ho, J.; Kilah, N. L.; Lin, C.-Y.; Salem, G.; Weir, M. L.; Willis, A. C.; Wild, S. B. Resolution and Improved Synthesis of (±)-Arsenicin A: a Natural Adamantane-type TetraarsenS

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A ical Possessing Strong Anti-acute Promelocytic Leukemia Cell Line Activity. Organometallics 2012, 31, 1808−1816. (160) Wentrup, C. The Benzyne Story. Aust. J. Chem. 2010, 63, 979− 986. (161) Groner, P.; Kukolich, S. G. Equilibrium Structure of Gas Phase o-Benzyne. J. Mol. Struct. 2006, 780−781, 178−181. (162) Muller, R. P.; Langlois, J.-M.; Ringnalda, M. N.; Friesner, R. A.; Goddard, W. A. A Generalized Direct Inversion in the Iterative Subspace Approach for Generalized Valence Bond Wave Functions. J. Chem. Phys. 1994, 100, 1226−1235.

T

DOI: 10.1021/acs.jpca.5b06893 J. Phys. Chem. A XXXX, XXX, XXX−XXX