THE JOURNAL OF CHEMICAL PHYSICS 142, 134112 (2015)
The variational subspace valence bond method Graham D. Fletcher Argonne National Laboratory, 9700 South Cass Ave., Lemont, Illinois 60439, USA
(Received 19 January 2015; accepted 23 March 2015; published online 7 April 2015) The variational subspace valence bond (VSVB) method based on overlapping orbitals is introduced. VSVB provides variational support against collapse for the optimization of overlapping linear combinations of atomic orbitals (OLCAOs) using modified orbital expansions, without recourse to orthogonalization. OLCAO have the advantage of being naturally localized, chemically intuitive (to individually model bonds and lone pairs, for example), and transferrable between different molecular systems. Such features are exploited to avoid key computational bottlenecks. Since the OLCAO can be doubly occupied, VSVB can access very large problems, and calculations on systems with several hundred atoms are presented. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4916743]
I. INTRODUCTION
In molecular electronic structure theory, Hartree-Fock (HF)1 provides the starting point for many ab initio methods. One purpose of HF is to provide a set of molecular orbitals for further analysis. A key property of molecular orbitals is their mutual orthogonality which has two major advantages: (1) it obviates the need to evaluate N! terms arising from a fermionic wave function, where N is the number of electrons and (2) it provides orbitals with variational support, preventing them from coalescing on the lowest energy solution to the Fock equations, yielding the shell structure of chemistry. However, orthogonality has several disadvantages that ultimately hinder applications of mainstream quantum chemistry to large problems. The tendency for molecular orbitals to be delocalized over the entire molecular system incurs high-order data structures for intermediate quantities that give rise to bottlenecks in storage (and/or communication for distributed, parallel processing), adversely impacting both efficient concurrency and the cost as a function of problem size. Molecular orbitals must be generated simultaneously, with further consequences for concurrency, and obtaining their initial guess requires a non-trivial energy calculation, often of low parallel scalability. Finally, molecular orbitals demand multiple spatial configurations in order to model excited states and regions of the nuclear potential energy surface away from the equilibrium geometry (consider, for example, the need for “anti-bonding” orbitals when breaking chemical bonds). Nevertheless, interest remains high in scaling quantumbased methods up to systems with hundreds to thousands of atoms, reflecting the impact of quantum effects well into the mesoscale, where ab initio methods obviate the need to anticipate exactly where and when the quantum effects will be important. In this aim, fragmentation approaches achieve both cost and concurrency scalability by subdividing a system into manageable parts. The fragment molecular orbital (FMO) method2 is, arguably, the most advanced of such approaches. However, the success of FMO is critically dependent upon the availability of a chemically intuitive fragmentation. The definition of fragment boundaries across chemical bonds can 0021-9606/2015/142(13)/134112/8/$30.00
be a source of instability3 in which orthogonality plays a central role, limiting this situation to small basis sets. The above issues cannot be entirely avoided so long as orthogonality is part of the model. Not surprisingly, then, the situation can be largely obviated by using overlapping orbitals. It is well-known from spin coupled valence bond (SCVB) theory that the wave function for a single product of overlapping orbitals can model many classes of chemical phenomena (see the work of Cooper and co-workers,4 and also van Lenthe and co-workers5). Overlapping orbitals also turn out to be “naturally” localized. For example, the minimal definition of a bonding orbital involves basis functions on two atomic centers. Thus, localization offers computational advantages because the orbital definitions rapidly become independent of the system size as the system grows. Localization, along with advances in computer hardware, ought to bring many more problems within the scope of methods based purely on overlapping orbitals, especially considering that determinants can be used to compute N! terms at a cost of order-N 3. However, the number of determinants remains an obstacle and is linked to the issue of variational support. The use of overlapping orbitals was first explored in the valence bond (VB) method in which determinant-based “structures” (being spin eigenfunctions) are built from atomic orbitals (AOs) and an eigenproblem in the structure basis is solved. VB tends to need many structures because the AO, individually, represent the chemistry rather poorly. Nevertheless, VB has retained a key role in chemical interpretation.6 SCVB solves the aforementioned problem using a single structure of overlapping linear combinations of AO (OLCAO), since the latter can directly model chemical objects such as bonds and lone pairs. The primary computational task of SCVB is to optimize the OLCAO, whose chemical forms then arise naturally by minimizing the energy. However, the number of determinants needed to spincouple N electrons tends to increase exponentially with N. The attempt to reduce the number of determinants by doubly occupying the OLCAO fails due to orbital coalescence (see, for example, Ref. 7, where it is termed “self-annihilation”). Furthermore, SCVB is reliant upon the antisymmetrizer for
142, 134112-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-2
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
variational support and, in general, will encounter coalescence for problems larger than eight or nine electrons, regardless of double-occupancy.8 The strategy of embedding SCVB within a complete active space self-consistent field (CASSCF) wave function ensures that the OLCAO are supported against collapse by a lower manifold of orthogonal orbitals.9 However, the computational limitations of SCVB become, essentially, those of CASSCF. Elsewhere, orbital overlap might be a component of a method, even if it is not the foundation.10 The effective fragment potential method of Jensen and co-workers is based on molecular orbitals, where the overlap between different fragments is approximated in order to model exchangerepulsion effects.11 In FMO, the overlap between orbitals of different fragments is purely incidental.2 On the other hand, the variational support of OLCAO, without recourse to orthogonalization, could access many computational advantages and extend their use far beyond what is currently possible. Such an approach is described in Secs. II–VIII. Notethat the present work significantly extends and generalizes earlier work12–14 and clarifies its connection with HF.
II. THEORY AND METHOD
For reasons that will become clear later, the current method will be called the variational subspace valence bond (VSVB) method. While VSVB shares many conceptual ingredients with SCVB,4,12 it can give very different results, as will be shown in Secs. III–VIII. In its most general form, the VSVB wave function is multi-structural, i ′ = Ci ΨVSVB . (1) ΨVSVB i
However, the focus of the present work is on the single configuration form and, hereafter, the superscript is dropped. The main departure from traditional SCVB is that OLCAO are allowed to be doubly occupied, with zero overall spin arising from them. Thus, the single-configuration VSVB wave function for N electrons is a product of “spin coupled” and “doubly occupied” wave functions, ˆ N ΞN ΨVSVB = AΦ , SC DOCC
(2)
where Aˆ is the antisymmetrizer, NSC is the number of spin coupled electrons, and NDOCC is the number of electrons occupying doubly occupied orbitals, with N = NSC + NDOCC .
(3)
The “SCVB” wave function involves a product of NSC orbitals coupled to an overall spin, S, NSC Φ NSC = *. Φi (i)+/ ΘSNSC. (4) , i=1 The spin coupling in (4) may be generally expanded over a basis of spin eigenfunctions ΘSNSC = d i ΘS,i (5) NSC i
(see Ref. 15, for example). As with SCVB, multiple spin couplings are an important component of VSVB although, for simplicity, Secs. III–VIII give examples with single spin couplings. The “doubly occupied” wave function is given by Ξ NDOCC =
NSC+N DOCC /2
¯ i ( j), Φi ( j − 1)Φ
(6)
i=NSC+1
where j = NSC + 2i, an overbar indicates β-spin, and its absence indicates α-spin. A useful feature of VSVB is that either (4) or (6) is allowed to be null, that is, NSC, NDOCC ≥ 0.
(7)
The orbitals in both (4) and (6) have LCAO form Φi =
Mi
Cµi η µ ,
(8)
µ ∈{η}
where η µ are the AO and Ci are their weights. Each LCAO is expanded over a subset of the molecular basis functions, {η}, in a manner described below, where Mi is the size of the subset. The OLCAO have the property
Φi Φ j , 0,i , j. (9) The molecular total energy may be expressed in terms of integrals and densities as follows: ΨVSVB Hˆ ΨVSVB EVSVB = ⟨ΨVSVB |ΨVSVB ⟩ = h i j Di j + gijkl Dijkl + Vnuc, (10) i, j
i, j,k,l
where h, g are the 1, 2-electron integrals over OLCAO weighted by densities, D, of matching order, and Vnuc is the nuclear repulsion energy. The D may be expressed in terms of the determinants of matrices containing the overlap integrals (9). The weights in (8) are determined variationally, subject only to the requirement of normalization. From a suitable guess at Ci , the ith OLCAO can be improved by setting the appropriate first derivatives of (10) to zero and solving H(i)Ci = S(i)Ci Ei ,
(11)
where the elements of the Hamiltonian matrix are (i) ˆ 1(1)Φ2(2) · · · η µ (i) · · · Φ N (N)ΘSN H µν = AΦ × Hˆ Φ1(1)Φ2(2) · · · η ν (i) · · · Φ N (N)ΘSN . (12) Hˆ is the Born-Oppenheimer molecular Hamiltonian. Analogous expressions apply for the unit operator. The order of (11) is Mi . The improved orbital corresponds to the lowest root of (11) for the ground state. Higher roots can be used to excite selected electrons because the solutions of (11) correspond to orthogonal N-electron states.13,14 Analogous expressions apply for the optimization of spin coupling weights (5). The simple first-order optimization method outlined in Eqs. (11) and (12) is intended primarily for illustrative purposes. In practice, the experience from SCVB4 shows that
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-3
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
second-order methods can significantly reduce the number of iterations needed to achieve convergence because the wave function parameters often interact strongly (for example, the OLCAO of a spin coupled pair). If S = NSC/2, (2) can be written as a single determinant. Alternatively, if two spatially distinct orbitals are to be singletcoupled, a spin function such as 1 Θ20 = √ (α(1) β(2) − β(1)α(2)) 2
(13)
may be used, needing two determinants. For n such pairs, 2n determinants would be needed. Given that the majority of chemical processes rearrange few electrons, n need not typically be large, provided OLCAO are used. However, the optimization of doubly occupied OLCAO remains to be addressed for which the variational properties of HF must be examined. The case of closed-shell, or Restricted HF (RHF), provides the simplest starting point. The RHF energy, ERHF, requires that the variational subspace of each occupied orbital contains a unique degree of freedom not found in the others. Thus, for a system with M basis functions and n˜ = NDOCC occupied orbitals, the total number of degrees of freedom for the ith orbital, Mi , is M − n˜ + 1. The latter condition is implicitly satisfied by orthogonality. In the context of OLCAO, recognizing that the weights of the normalized functions of the molecular basis set constitute equivalent units of variational freedom, the orbital expansions in (8) can be chosen to satisfy the same criterion without assuming orthogonality.14 Compared to RHF, the Mi may be preserved by explicitly removing n˜ − 1 degrees of freedom from each OLCAO expansion (8) while simultaneously restoring n˜ − 1 degrees of freedom by allowing overlap. Since the structure of the degrees of freedom is preserved from RHF, the model quality is exactly that of RHF. In so far as ERHF is obtained, the OLCAO are variationally supported. When Mi = M − n˜ + 1, the method will be referred to as variational subspace HartreeFock (VSHF), with EVSHF = ERHF. Note that terms such as “degrees of freedom,” “AO,” and “basis function” are, in the present work, synonymous. However, the distinctions can be useful if, for example, a degree of freedom itself corresponds to a LCAO, as may be the case when forming a symmetryadapted basis set. It is necessary to decide which degrees of freedom best represent the intended OLCAO types (core, lone pair, bonding, etc.), and declare them to be the unique degrees of freedom for each OLCAO. Where there are equivalent centers, a symmetry-adapted basis might be convenient in this respect, though this is not a requirement of VSHF. The latter choices define the variational subspaces and, therefore, the OLCAO expansions (8). In VSHF, n˜ = NDOCC, though, in general, n˜ counts the distinct variational subspaces. When solving (11), collapse is prevented by the non-coincident nature of the variational subspaces. During optimization, the energy penalty associated with OLCAO abandoning their unique degrees of freedom discourages coalescence. Thus, the use of variational subspaces is analogous to orthogonality. As ever, numerical success depends on the suitability of the molecular basis set and the wave function guess. However, the
TABLE I. Expansions of two doubly occupied orbitals to model the 1 S ground state of Be over a basis of five functions. The guess and optimum weights are shown. Variational subspaces AO
Φ1
s1 s2 s3 s4 s5
• • • •
Guess
Φ2
Φ1 1.0 0.0 0.0 0.0
• • • •
Optimized
Φ2
0.0 0.0 0.0 1.0
Φ1 0.9954 −0.0013 0.0118 −0.0071
Φ2
0.3025 0.2532 0.5818 0.2920
variational subspaces may also be obtained by rotation of the canonical RHF orbitals without an optimization procedure (see Appendix A). Analogous rules apply for open-shell systems.14 Note that the variational subspaces may be described as “partially” disjoint in the following manner. Let S´ be the set of degrees of freedom shared by all OLCAO. For M = n, ˜ S´ is the null set, and the variational subspaces are disjoint. The existence of a non-null S´ (see Tables I and II, below) makes the variational subspaces partially disjoint. The situation, Mi < M − n˜ + 1, is also valid, allowing the OLCAO expansions to be truncated relative to the full molecular basis set. Since the truncated model does not reproduce ERHF, the level of theory is unique to VSVB, with energy EVSVB. Otherwise, with NSC > 0, VSVB may be superior to RHF if, for example, singlet-coupled pairs of distinct OLCAO are allowed to share the same variational subspace. When optimized, such orbital pairs typically become polarized, reflecting a correlation. Thus, VSVB allows SCVB to be embedded in a way that addresses the shortcomings in  noted earlier. SCVB can be regarded as a special case of VSVB where NDOCC = 0, Mi = M, for all i, and the singly occupied OLCAO must be coupled to an overall spin using multiple determinants. While, in principle, VSHF could be used to obtain ERHF, it might not be the most efficient way (although VSHF offers computational advantages—see Sec. VII), rather VSHF provides validation and a reference point for VSVB, as will be shown next. III. EXAMPLE OF THE BERYLLIUM ATOM
The use of VSVB is illustrated with the 1 S ground state of the beryllium atom, where the minimal wave function involves TABLE II. Expansions of the two occupied orbitals to model the 1 S ground state of Be over a basis of five functions, for a different choice of variational subspace to that given in Table I. Variational subspaces AO
Φ1
s1 s2 s3 s4 s5
• • • •
Φ2
• • • •
Optimized Φ1 0.9954 0.0024 0.0149 0.0036
Φ2
0.3025 0.2532 0.5818 0.2920
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-4
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
just two doubly occupied orbitals, Φ1 and Φ2 (n˜ = 2). A computer program was developed to perform the calculations (see Appendix B for details). The basis set consists of the five s-type functions from the cc-pVQZ basis for Be,16 denoted “s i ,” i = 1-5, in which the first two are contracted over eight primitive gaussian functions and the remaining three are uncontracted. The VSHF calculation is summarized in Table I, where the second and third columns use dots to indicate that a given basis function is in the expansion of Φ1 and Φ2, respectively. In Table I, s1 is chosen as the unique degree of freedom for Φ1 while s5 is the unique degree of freedom for Φ2. A simple initial guess is given in columns four and five weighting the unique degrees of freedom with unity while other degrees of freedom begin with weights of zero. Following the optimization of the total energy to 10−8 AU, the final definitions of Φ1 and Φ2, which approximate the 1s and 2s shells of Be, are given in columns six and seven, respectively, to four figures. The energy, EVSHF = −14.572 968 12 AU, was verified to match the RHF energy by performing a standard RHF calculation using GAMESS.17 The overlap integral of Φ1 and Φ2, ⟨Φ1|Φ2⟩, has the value 0.0042. Table II summarizes a similar calculation to that shown in Table I in order to illustrate a different variational subspace, where s4 is chosen as the unique degree of freedom for Φ2, instead of s5. For brevity, Table II omits the wave function guess. Table II shows that the different choice of variational subspace causes an adjustment of the weights for Φ1, with ⟨Φ1|Φ2 = 0.0080⟩, while the total energy, EVSHF = −14.572 968 12 AU, matches that obtained earlier, and the impact on Φ2 is negligible. An advantage of VSVB is that it allows the HF-quality VSHF model to be improved by embedding SCVB. Table III gives examples of where the doubly occupied Φ1 and Φ2 are “split” into correlated orbital pairs. The split may be achieved by first replacing the original doubly occupied orbital in the wave function with two singly occupied copies, coupled to a singlet using a form analogous to (13), then applying opposing perturbations to the LCAO weights and re-optimizing. Since the resulting VSVB wave functions are qualitatively superior to HF with Mi = M − n˜ + 1, their energies are all lower than ERHF. To aid in the chemical interpretation of Table III, the orbitals are now labeled by the 1s and 2s atomic shells of Be they approximate; an underline indicates beta spin, as in (6), parentheses indicate coupling to singlet spin based on (13), and primes distinguish the members of singlet pairs. The relaxation energy relative to RHF is given in the right column.
The column labeled Ndet gives the number of determinants in the wave function. To provide a reference point, the first row of Table III represents the optimized wave function of Table I in the chemical notation. The wave function in the second row splits the valence shell 2s orbital into a spin-coupled pair, bringing radial correlations. The third row shows the analogous modification for the core 1s shell, where the relaxation is larger due to the core shell being more tightly bound, while the fourth row applies the splitting to both shells. It can be seen that the relaxations in the two shells are nearly additive. The final row allows comparison with the SCVB model in which full expansions for all orbitals are used and no degrees of freedom are omitted. The SCVB energy is slightly lower (by > n, ˜ VSVB requires significantly less memory than SCF. In summary, the principal computational advantages compared to SCF are as follows. 1. Wave function guess calculation not imperative. 2. Eigenproblems typically small and localized, rather than global. 3. Lower memory requirements. Such advantages are strongly aligned with the current trend toward dense architectures exemplified by accelerators and co-processors. Furthermore, it follows from SCVB theory that the single configuration VSVB wave function can be readily extended to model 1. bond formation, majority (or all) of the nuclear potential energy surface; 2. open- and closed-shell systems; 3. resonances (multiple spin couplings) and degenerate states; 4. ground and excited states.12–14 Thus, VSVB provides a versatile and scalable chemical model applicable to problems with hundreds of atoms in which overlapping, localized functions are used throughout. In this sense, VSVB obviates the N! problem.
VIII. FUTURE WORK
As mentioned in Sec. II, the level of theory equivalent to VSVB depends on many details of the specific application and, thus, may be hard to estimate. However, as with SCVB, the maximum level that may be attained is close to CASSCF.9 Further improvement will require an account of dynamic electron correlation effects. The natural strategy for a model based on locality such as VSVB is to begin with terms linear in the inter-electronic distances (“r 12”) in the wave function. Continuing the theme of locality, cost scaling may be improved with approximations to the asymptotic forms of the coulomb integrals over OLCAO, while storage for computing the density matrix elements in (3) may be reduced to asymptotic using a theorem arising from overlap decay.8 Meanwhile, the greater arithmetic load of VSVB compared to SCF can be mitigated through significant levels of concurrency, as mentioned in Sec. VI, together with vectorization. Major forthcoming topics include 1. use of explicit density normalization to reduce memory scaling to asymptotic, 2. the methodology for open-shell systems, 3. long-range coulomb integral approximations,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-8
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
4. parallel scalability and vectorization, 5. improve convergence by incorporation of second-order methods (e.g., Newton-Raphson), and 6. incorporation of dynamic electron correlation effects based on r 12.
ACKNOWLEDGMENTS
The author wishes to thank Brian Sutcliffe at Université Libre de Bruxelles, Ray Bair and Spencer Pruitt at Argonne National Laboratory, and Fred Manby at the University of Bristol (UK), for helpful comments. The use of MacMolPlot19 for images and GAMESS17 for RHF calculations is acknowledged. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DEAC02-06CH11357. This work is dedicated to the memory of Dr. Graham Doggett.
APPENDIX A: ROTATION METHOD
(A1)
where θ = tan−1
Cµi Cµ j
The present work was performed using a suite of software tools written by the author. The code for computing integrals is based on the Rys polynomial method.20 The author is grateful for the existence of software for solving common mathematical problems. Eigenproblems were solved using the recommended sequence of “eispack” routines.21 Determinants (for density matrix elements) were computed using a “LU decomposition” algorithm due to Dongarra.22 The total number of lines in the energy program is currently 7450, including comments. APPENDIX C: ETHANE CLUSTERS
In the 91-ethane cluster, the experimental lattice parameters for the cubic structure of solid ethane were used.23 The observed level of crystal disorder was approximated by choosing the three Cartesian angles about the mid-point of each CC bond relative to the reference frame to be random numbers between zero and 0.05 radians.24 1C.
The rotation method may proceed as follows. Step 1. For each canonical occupied RHF orbital, Φi , choose a unique degree of freedom, µ, subject to Cµi , 0. Step 2. From Φi , proceed to remove the degrees of freedom of Φ j , j , i, by rotation, as follows: ϕi′ = cos θϕi − sin θϕ j , ϕ ′j = sin θϕi + cos θϕ j ,
APPENDIX B: SOFTWARE IMPLEMENTATION
(A2)
and µ labels the unique degree of freedom of Φ j . In principle, the only condition upon j is the numerical stability of (8), provided the required degree of freedom is referenced. Step 3. Repeat step 2, for the remaining j , i. Step 4. Repeat steps 2 and 3 iteratively, as subsequent rotations may undo the effects of previous ones, until all weights in question are sufficiently small for use. Step 5. Repeat steps 2-4 for all remaining occupied orbitals, beginning with the canonical RHF orbitals each time, and saving the final Φi . Since ERHF is invariant to rotations of the doubly occupied orbitals, the combined final orbitals yield EVSHF = ERHF. The rotation method is computationally trivial compared to the VSHF method in Secs. I–VIII. Although both routes yield ERHF, neither is guaranteed to generate satisfactory orbitals and visualization is strongly recommended before any orbitals are used. In this respect, VSHF may yield more reliable results by virtue of needing variational stability to complete.
C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951). G. Fedorov and K. Kitaura, J. Phys. Chem. A 111, 6904 (2007). 3D. G. Fedorov and K. Kitaura, in Modern Methods for Theoretical Physical Chemistry and Biopolymers, edited by E. B. Starikov, J. P. Lewis, and S. Tanaka (Elsevier, Amsterdam, 2006), pp. 3-38. 4D. L. Cooper, J. Gerratt, and M. Raimondi, Chem. Rev. 91, 929 (1991), and references therein. 5Z. Rashid and J. H. v. Lenthe, J. Chem. Phys. 138, 54105 (2013). 6D. Danovich, P. C. Hiberty, W. Wu, H. S. Rzepa, and S. Shaik, Chem. Eur. J. 20, 6220 (2014); C. Angeli, R. Cimiraglia, and J.-P. Malrieu, J. Chem. Educ. 85, 150 (2008). 7D. L. Cooper and P. B. Karadakov, Int. Rev. Phys. Chem. 28, 169 (2009). 8G. D. Fletcher (unpublished). 9T. Thorsteinsson, D. L. Cooper, J. Gerratt, P. B. Karadov, and M. Raimondi, Theor. Chim. Acta 93, 343 (1996). 10G. E. Scuseria, C. A. Jimenez-Hoyos, T. M. Henderson, K. Samanta, and J. K. Ellis, J. Chem. Phys. 135, 124108 (2011). 11D. D. Kemp, J. M. Rintelman, M. S. Gordon, and J. H. Jensen, Theor. Chem. Acc. 125, 481 (2010). 12G. D. Fletcher, G. Doggett, and A. S. Howard, Phys. Rev. A 46, 5459 (1992). 13G. Doggett and G. D. Fletcher, J. Mol. Struct.: THEOCHEM 260, 313 (1992); G. Doggett, G. D. Fletcher, and F. R. Manby, J. Mol. Struct. 300, 191 (1992); F. R. Manby, G. Doggett, and G. D. Fletcher, J. Mol. Struct.: THEOCHEM 343, 63 (1995). 14G. D. Fletcher, Ph.D. thesis, University of York, UK, 1993. 15R. Pauncz, Spin Eigenfunctions (Plenum, New York, 1979). 16T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 17M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. J. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 18J. J. P. Stewart, J. Comput. Chem. 10, 209 (1989). 19B. M. Bode and M. S. Gordon, J. Mol. Graphics Modell. 16, 133 (1999). 20M. Dupuis, J. Rys, and H. F. King, J. Chem. Phys. 65, 111 (1976). 21B. S. Garbow, Comput. Phys. Commun. 7, 179 (1974). 22J. Dongarra, ACM Trans. Math. Software 10, 219 (1984). 23G. J. H. van Nes and A. Vos, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 34, 1947 (1978). 24See supplementary material at http://dx.doi.org/10.1063/1.4916743 for the atomic coordinates of the structure. 2D.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.
The variational subspace valence bond method Graham D. Fletcher Argonne National Laboratory, 9700 South Cass Ave., Lemont, Illinois 60439, USA
(Received 19 January 2015; accepted 23 March 2015; published online 7 April 2015) The variational subspace valence bond (VSVB) method based on overlapping orbitals is introduced. VSVB provides variational support against collapse for the optimization of overlapping linear combinations of atomic orbitals (OLCAOs) using modified orbital expansions, without recourse to orthogonalization. OLCAO have the advantage of being naturally localized, chemically intuitive (to individually model bonds and lone pairs, for example), and transferrable between different molecular systems. Such features are exploited to avoid key computational bottlenecks. Since the OLCAO can be doubly occupied, VSVB can access very large problems, and calculations on systems with several hundred atoms are presented. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4916743]
I. INTRODUCTION
In molecular electronic structure theory, Hartree-Fock (HF)1 provides the starting point for many ab initio methods. One purpose of HF is to provide a set of molecular orbitals for further analysis. A key property of molecular orbitals is their mutual orthogonality which has two major advantages: (1) it obviates the need to evaluate N! terms arising from a fermionic wave function, where N is the number of electrons and (2) it provides orbitals with variational support, preventing them from coalescing on the lowest energy solution to the Fock equations, yielding the shell structure of chemistry. However, orthogonality has several disadvantages that ultimately hinder applications of mainstream quantum chemistry to large problems. The tendency for molecular orbitals to be delocalized over the entire molecular system incurs high-order data structures for intermediate quantities that give rise to bottlenecks in storage (and/or communication for distributed, parallel processing), adversely impacting both efficient concurrency and the cost as a function of problem size. Molecular orbitals must be generated simultaneously, with further consequences for concurrency, and obtaining their initial guess requires a non-trivial energy calculation, often of low parallel scalability. Finally, molecular orbitals demand multiple spatial configurations in order to model excited states and regions of the nuclear potential energy surface away from the equilibrium geometry (consider, for example, the need for “anti-bonding” orbitals when breaking chemical bonds). Nevertheless, interest remains high in scaling quantumbased methods up to systems with hundreds to thousands of atoms, reflecting the impact of quantum effects well into the mesoscale, where ab initio methods obviate the need to anticipate exactly where and when the quantum effects will be important. In this aim, fragmentation approaches achieve both cost and concurrency scalability by subdividing a system into manageable parts. The fragment molecular orbital (FMO) method2 is, arguably, the most advanced of such approaches. However, the success of FMO is critically dependent upon the availability of a chemically intuitive fragmentation. The definition of fragment boundaries across chemical bonds can 0021-9606/2015/142(13)/134112/8/$30.00
be a source of instability3 in which orthogonality plays a central role, limiting this situation to small basis sets. The above issues cannot be entirely avoided so long as orthogonality is part of the model. Not surprisingly, then, the situation can be largely obviated by using overlapping orbitals. It is well-known from spin coupled valence bond (SCVB) theory that the wave function for a single product of overlapping orbitals can model many classes of chemical phenomena (see the work of Cooper and co-workers,4 and also van Lenthe and co-workers5). Overlapping orbitals also turn out to be “naturally” localized. For example, the minimal definition of a bonding orbital involves basis functions on two atomic centers. Thus, localization offers computational advantages because the orbital definitions rapidly become independent of the system size as the system grows. Localization, along with advances in computer hardware, ought to bring many more problems within the scope of methods based purely on overlapping orbitals, especially considering that determinants can be used to compute N! terms at a cost of order-N 3. However, the number of determinants remains an obstacle and is linked to the issue of variational support. The use of overlapping orbitals was first explored in the valence bond (VB) method in which determinant-based “structures” (being spin eigenfunctions) are built from atomic orbitals (AOs) and an eigenproblem in the structure basis is solved. VB tends to need many structures because the AO, individually, represent the chemistry rather poorly. Nevertheless, VB has retained a key role in chemical interpretation.6 SCVB solves the aforementioned problem using a single structure of overlapping linear combinations of AO (OLCAO), since the latter can directly model chemical objects such as bonds and lone pairs. The primary computational task of SCVB is to optimize the OLCAO, whose chemical forms then arise naturally by minimizing the energy. However, the number of determinants needed to spincouple N electrons tends to increase exponentially with N. The attempt to reduce the number of determinants by doubly occupying the OLCAO fails due to orbital coalescence (see, for example, Ref. 7, where it is termed “self-annihilation”). Furthermore, SCVB is reliant upon the antisymmetrizer for
142, 134112-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-2
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
variational support and, in general, will encounter coalescence for problems larger than eight or nine electrons, regardless of double-occupancy.8 The strategy of embedding SCVB within a complete active space self-consistent field (CASSCF) wave function ensures that the OLCAO are supported against collapse by a lower manifold of orthogonal orbitals.9 However, the computational limitations of SCVB become, essentially, those of CASSCF. Elsewhere, orbital overlap might be a component of a method, even if it is not the foundation.10 The effective fragment potential method of Jensen and co-workers is based on molecular orbitals, where the overlap between different fragments is approximated in order to model exchangerepulsion effects.11 In FMO, the overlap between orbitals of different fragments is purely incidental.2 On the other hand, the variational support of OLCAO, without recourse to orthogonalization, could access many computational advantages and extend their use far beyond what is currently possible. Such an approach is described in Secs. II–VIII. Notethat the present work significantly extends and generalizes earlier work12–14 and clarifies its connection with HF.
II. THEORY AND METHOD
For reasons that will become clear later, the current method will be called the variational subspace valence bond (VSVB) method. While VSVB shares many conceptual ingredients with SCVB,4,12 it can give very different results, as will be shown in Secs. III–VIII. In its most general form, the VSVB wave function is multi-structural, i ′ = Ci ΨVSVB . (1) ΨVSVB i
However, the focus of the present work is on the single configuration form and, hereafter, the superscript is dropped. The main departure from traditional SCVB is that OLCAO are allowed to be doubly occupied, with zero overall spin arising from them. Thus, the single-configuration VSVB wave function for N electrons is a product of “spin coupled” and “doubly occupied” wave functions, ˆ N ΞN ΨVSVB = AΦ , SC DOCC
(2)
where Aˆ is the antisymmetrizer, NSC is the number of spin coupled electrons, and NDOCC is the number of electrons occupying doubly occupied orbitals, with N = NSC + NDOCC .
(3)
The “SCVB” wave function involves a product of NSC orbitals coupled to an overall spin, S, NSC Φ NSC = *. Φi (i)+/ ΘSNSC. (4) , i=1 The spin coupling in (4) may be generally expanded over a basis of spin eigenfunctions ΘSNSC = d i ΘS,i (5) NSC i
(see Ref. 15, for example). As with SCVB, multiple spin couplings are an important component of VSVB although, for simplicity, Secs. III–VIII give examples with single spin couplings. The “doubly occupied” wave function is given by Ξ NDOCC =
NSC+N DOCC /2
¯ i ( j), Φi ( j − 1)Φ
(6)
i=NSC+1
where j = NSC + 2i, an overbar indicates β-spin, and its absence indicates α-spin. A useful feature of VSVB is that either (4) or (6) is allowed to be null, that is, NSC, NDOCC ≥ 0.
(7)
The orbitals in both (4) and (6) have LCAO form Φi =
Mi
Cµi η µ ,
(8)
µ ∈{η}
where η µ are the AO and Ci are their weights. Each LCAO is expanded over a subset of the molecular basis functions, {η}, in a manner described below, where Mi is the size of the subset. The OLCAO have the property
Φi Φ j , 0,i , j. (9) The molecular total energy may be expressed in terms of integrals and densities as follows: ΨVSVB Hˆ ΨVSVB EVSVB = ⟨ΨVSVB |ΨVSVB ⟩ = h i j Di j + gijkl Dijkl + Vnuc, (10) i, j
i, j,k,l
where h, g are the 1, 2-electron integrals over OLCAO weighted by densities, D, of matching order, and Vnuc is the nuclear repulsion energy. The D may be expressed in terms of the determinants of matrices containing the overlap integrals (9). The weights in (8) are determined variationally, subject only to the requirement of normalization. From a suitable guess at Ci , the ith OLCAO can be improved by setting the appropriate first derivatives of (10) to zero and solving H(i)Ci = S(i)Ci Ei ,
(11)
where the elements of the Hamiltonian matrix are (i) ˆ 1(1)Φ2(2) · · · η µ (i) · · · Φ N (N)ΘSN H µν = AΦ × Hˆ Φ1(1)Φ2(2) · · · η ν (i) · · · Φ N (N)ΘSN . (12) Hˆ is the Born-Oppenheimer molecular Hamiltonian. Analogous expressions apply for the unit operator. The order of (11) is Mi . The improved orbital corresponds to the lowest root of (11) for the ground state. Higher roots can be used to excite selected electrons because the solutions of (11) correspond to orthogonal N-electron states.13,14 Analogous expressions apply for the optimization of spin coupling weights (5). The simple first-order optimization method outlined in Eqs. (11) and (12) is intended primarily for illustrative purposes. In practice, the experience from SCVB4 shows that
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-3
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
second-order methods can significantly reduce the number of iterations needed to achieve convergence because the wave function parameters often interact strongly (for example, the OLCAO of a spin coupled pair). If S = NSC/2, (2) can be written as a single determinant. Alternatively, if two spatially distinct orbitals are to be singletcoupled, a spin function such as 1 Θ20 = √ (α(1) β(2) − β(1)α(2)) 2
(13)
may be used, needing two determinants. For n such pairs, 2n determinants would be needed. Given that the majority of chemical processes rearrange few electrons, n need not typically be large, provided OLCAO are used. However, the optimization of doubly occupied OLCAO remains to be addressed for which the variational properties of HF must be examined. The case of closed-shell, or Restricted HF (RHF), provides the simplest starting point. The RHF energy, ERHF, requires that the variational subspace of each occupied orbital contains a unique degree of freedom not found in the others. Thus, for a system with M basis functions and n˜ = NDOCC occupied orbitals, the total number of degrees of freedom for the ith orbital, Mi , is M − n˜ + 1. The latter condition is implicitly satisfied by orthogonality. In the context of OLCAO, recognizing that the weights of the normalized functions of the molecular basis set constitute equivalent units of variational freedom, the orbital expansions in (8) can be chosen to satisfy the same criterion without assuming orthogonality.14 Compared to RHF, the Mi may be preserved by explicitly removing n˜ − 1 degrees of freedom from each OLCAO expansion (8) while simultaneously restoring n˜ − 1 degrees of freedom by allowing overlap. Since the structure of the degrees of freedom is preserved from RHF, the model quality is exactly that of RHF. In so far as ERHF is obtained, the OLCAO are variationally supported. When Mi = M − n˜ + 1, the method will be referred to as variational subspace HartreeFock (VSHF), with EVSHF = ERHF. Note that terms such as “degrees of freedom,” “AO,” and “basis function” are, in the present work, synonymous. However, the distinctions can be useful if, for example, a degree of freedom itself corresponds to a LCAO, as may be the case when forming a symmetryadapted basis set. It is necessary to decide which degrees of freedom best represent the intended OLCAO types (core, lone pair, bonding, etc.), and declare them to be the unique degrees of freedom for each OLCAO. Where there are equivalent centers, a symmetry-adapted basis might be convenient in this respect, though this is not a requirement of VSHF. The latter choices define the variational subspaces and, therefore, the OLCAO expansions (8). In VSHF, n˜ = NDOCC, though, in general, n˜ counts the distinct variational subspaces. When solving (11), collapse is prevented by the non-coincident nature of the variational subspaces. During optimization, the energy penalty associated with OLCAO abandoning their unique degrees of freedom discourages coalescence. Thus, the use of variational subspaces is analogous to orthogonality. As ever, numerical success depends on the suitability of the molecular basis set and the wave function guess. However, the
TABLE I. Expansions of two doubly occupied orbitals to model the 1 S ground state of Be over a basis of five functions. The guess and optimum weights are shown. Variational subspaces AO
Φ1
s1 s2 s3 s4 s5
• • • •
Guess
Φ2
Φ1 1.0 0.0 0.0 0.0
• • • •
Optimized
Φ2
0.0 0.0 0.0 1.0
Φ1 0.9954 −0.0013 0.0118 −0.0071
Φ2
0.3025 0.2532 0.5818 0.2920
variational subspaces may also be obtained by rotation of the canonical RHF orbitals without an optimization procedure (see Appendix A). Analogous rules apply for open-shell systems.14 Note that the variational subspaces may be described as “partially” disjoint in the following manner. Let S´ be the set of degrees of freedom shared by all OLCAO. For M = n, ˜ S´ is the null set, and the variational subspaces are disjoint. The existence of a non-null S´ (see Tables I and II, below) makes the variational subspaces partially disjoint. The situation, Mi < M − n˜ + 1, is also valid, allowing the OLCAO expansions to be truncated relative to the full molecular basis set. Since the truncated model does not reproduce ERHF, the level of theory is unique to VSVB, with energy EVSVB. Otherwise, with NSC > 0, VSVB may be superior to RHF if, for example, singlet-coupled pairs of distinct OLCAO are allowed to share the same variational subspace. When optimized, such orbital pairs typically become polarized, reflecting a correlation. Thus, VSVB allows SCVB to be embedded in a way that addresses the shortcomings in  noted earlier. SCVB can be regarded as a special case of VSVB where NDOCC = 0, Mi = M, for all i, and the singly occupied OLCAO must be coupled to an overall spin using multiple determinants. While, in principle, VSHF could be used to obtain ERHF, it might not be the most efficient way (although VSHF offers computational advantages—see Sec. VII), rather VSHF provides validation and a reference point for VSVB, as will be shown next. III. EXAMPLE OF THE BERYLLIUM ATOM
The use of VSVB is illustrated with the 1 S ground state of the beryllium atom, where the minimal wave function involves TABLE II. Expansions of the two occupied orbitals to model the 1 S ground state of Be over a basis of five functions, for a different choice of variational subspace to that given in Table I. Variational subspaces AO
Φ1
s1 s2 s3 s4 s5
• • • •
Φ2
• • • •
Optimized Φ1 0.9954 0.0024 0.0149 0.0036
Φ2
0.3025 0.2532 0.5818 0.2920
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-4
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
just two doubly occupied orbitals, Φ1 and Φ2 (n˜ = 2). A computer program was developed to perform the calculations (see Appendix B for details). The basis set consists of the five s-type functions from the cc-pVQZ basis for Be,16 denoted “s i ,” i = 1-5, in which the first two are contracted over eight primitive gaussian functions and the remaining three are uncontracted. The VSHF calculation is summarized in Table I, where the second and third columns use dots to indicate that a given basis function is in the expansion of Φ1 and Φ2, respectively. In Table I, s1 is chosen as the unique degree of freedom for Φ1 while s5 is the unique degree of freedom for Φ2. A simple initial guess is given in columns four and five weighting the unique degrees of freedom with unity while other degrees of freedom begin with weights of zero. Following the optimization of the total energy to 10−8 AU, the final definitions of Φ1 and Φ2, which approximate the 1s and 2s shells of Be, are given in columns six and seven, respectively, to four figures. The energy, EVSHF = −14.572 968 12 AU, was verified to match the RHF energy by performing a standard RHF calculation using GAMESS.17 The overlap integral of Φ1 and Φ2, ⟨Φ1|Φ2⟩, has the value 0.0042. Table II summarizes a similar calculation to that shown in Table I in order to illustrate a different variational subspace, where s4 is chosen as the unique degree of freedom for Φ2, instead of s5. For brevity, Table II omits the wave function guess. Table II shows that the different choice of variational subspace causes an adjustment of the weights for Φ1, with ⟨Φ1|Φ2 = 0.0080⟩, while the total energy, EVSHF = −14.572 968 12 AU, matches that obtained earlier, and the impact on Φ2 is negligible. An advantage of VSVB is that it allows the HF-quality VSHF model to be improved by embedding SCVB. Table III gives examples of where the doubly occupied Φ1 and Φ2 are “split” into correlated orbital pairs. The split may be achieved by first replacing the original doubly occupied orbital in the wave function with two singly occupied copies, coupled to a singlet using a form analogous to (13), then applying opposing perturbations to the LCAO weights and re-optimizing. Since the resulting VSVB wave functions are qualitatively superior to HF with Mi = M − n˜ + 1, their energies are all lower than ERHF. To aid in the chemical interpretation of Table III, the orbitals are now labeled by the 1s and 2s atomic shells of Be they approximate; an underline indicates beta spin, as in (6), parentheses indicate coupling to singlet spin based on (13), and primes distinguish the members of singlet pairs. The relaxation energy relative to RHF is given in the right column.
The column labeled Ndet gives the number of determinants in the wave function. To provide a reference point, the first row of Table III represents the optimized wave function of Table I in the chemical notation. The wave function in the second row splits the valence shell 2s orbital into a spin-coupled pair, bringing radial correlations. The third row shows the analogous modification for the core 1s shell, where the relaxation is larger due to the core shell being more tightly bound, while the fourth row applies the splitting to both shells. It can be seen that the relaxations in the two shells are nearly additive. The final row allows comparison with the SCVB model in which full expansions for all orbitals are used and no degrees of freedom are omitted. The SCVB energy is slightly lower (by > n, ˜ VSVB requires significantly less memory than SCF. In summary, the principal computational advantages compared to SCF are as follows. 1. Wave function guess calculation not imperative. 2. Eigenproblems typically small and localized, rather than global. 3. Lower memory requirements. Such advantages are strongly aligned with the current trend toward dense architectures exemplified by accelerators and co-processors. Furthermore, it follows from SCVB theory that the single configuration VSVB wave function can be readily extended to model 1. bond formation, majority (or all) of the nuclear potential energy surface; 2. open- and closed-shell systems; 3. resonances (multiple spin couplings) and degenerate states; 4. ground and excited states.12–14 Thus, VSVB provides a versatile and scalable chemical model applicable to problems with hundreds of atoms in which overlapping, localized functions are used throughout. In this sense, VSVB obviates the N! problem.
VIII. FUTURE WORK
As mentioned in Sec. II, the level of theory equivalent to VSVB depends on many details of the specific application and, thus, may be hard to estimate. However, as with SCVB, the maximum level that may be attained is close to CASSCF.9 Further improvement will require an account of dynamic electron correlation effects. The natural strategy for a model based on locality such as VSVB is to begin with terms linear in the inter-electronic distances (“r 12”) in the wave function. Continuing the theme of locality, cost scaling may be improved with approximations to the asymptotic forms of the coulomb integrals over OLCAO, while storage for computing the density matrix elements in (3) may be reduced to asymptotic using a theorem arising from overlap decay.8 Meanwhile, the greater arithmetic load of VSVB compared to SCF can be mitigated through significant levels of concurrency, as mentioned in Sec. VI, together with vectorization. Major forthcoming topics include 1. use of explicit density normalization to reduce memory scaling to asymptotic, 2. the methodology for open-shell systems, 3. long-range coulomb integral approximations,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
134112-8
Graham D. Fletcher
J. Chem. Phys. 142, 134112 (2015)
4. parallel scalability and vectorization, 5. improve convergence by incorporation of second-order methods (e.g., Newton-Raphson), and 6. incorporation of dynamic electron correlation effects based on r 12.
ACKNOWLEDGMENTS
The author wishes to thank Brian Sutcliffe at Université Libre de Bruxelles, Ray Bair and Spencer Pruitt at Argonne National Laboratory, and Fred Manby at the University of Bristol (UK), for helpful comments. The use of MacMolPlot19 for images and GAMESS17 for RHF calculations is acknowledged. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DEAC02-06CH11357. This work is dedicated to the memory of Dr. Graham Doggett.
APPENDIX A: ROTATION METHOD
(A1)
where θ = tan−1
Cµi Cµ j
The present work was performed using a suite of software tools written by the author. The code for computing integrals is based on the Rys polynomial method.20 The author is grateful for the existence of software for solving common mathematical problems. Eigenproblems were solved using the recommended sequence of “eispack” routines.21 Determinants (for density matrix elements) were computed using a “LU decomposition” algorithm due to Dongarra.22 The total number of lines in the energy program is currently 7450, including comments. APPENDIX C: ETHANE CLUSTERS
In the 91-ethane cluster, the experimental lattice parameters for the cubic structure of solid ethane were used.23 The observed level of crystal disorder was approximated by choosing the three Cartesian angles about the mid-point of each CC bond relative to the reference frame to be random numbers between zero and 0.05 radians.24 1C.
The rotation method may proceed as follows. Step 1. For each canonical occupied RHF orbital, Φi , choose a unique degree of freedom, µ, subject to Cµi , 0. Step 2. From Φi , proceed to remove the degrees of freedom of Φ j , j , i, by rotation, as follows: ϕi′ = cos θϕi − sin θϕ j , ϕ ′j = sin θϕi + cos θϕ j ,
APPENDIX B: SOFTWARE IMPLEMENTATION
(A2)
and µ labels the unique degree of freedom of Φ j . In principle, the only condition upon j is the numerical stability of (8), provided the required degree of freedom is referenced. Step 3. Repeat step 2, for the remaining j , i. Step 4. Repeat steps 2 and 3 iteratively, as subsequent rotations may undo the effects of previous ones, until all weights in question are sufficiently small for use. Step 5. Repeat steps 2-4 for all remaining occupied orbitals, beginning with the canonical RHF orbitals each time, and saving the final Φi . Since ERHF is invariant to rotations of the doubly occupied orbitals, the combined final orbitals yield EVSHF = ERHF. The rotation method is computationally trivial compared to the VSHF method in Secs. I–VIII. Although both routes yield ERHF, neither is guaranteed to generate satisfactory orbitals and visualization is strongly recommended before any orbitals are used. In this respect, VSHF may yield more reliable results by virtue of needing variational stability to complete.
C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951). G. Fedorov and K. Kitaura, J. Phys. Chem. A 111, 6904 (2007). 3D. G. Fedorov and K. Kitaura, in Modern Methods for Theoretical Physical Chemistry and Biopolymers, edited by E. B. Starikov, J. P. Lewis, and S. Tanaka (Elsevier, Amsterdam, 2006), pp. 3-38. 4D. L. Cooper, J. Gerratt, and M. Raimondi, Chem. Rev. 91, 929 (1991), and references therein. 5Z. Rashid and J. H. v. Lenthe, J. Chem. Phys. 138, 54105 (2013). 6D. Danovich, P. C. Hiberty, W. Wu, H. S. Rzepa, and S. Shaik, Chem. Eur. J. 20, 6220 (2014); C. Angeli, R. Cimiraglia, and J.-P. Malrieu, J. Chem. Educ. 85, 150 (2008). 7D. L. Cooper and P. B. Karadakov, Int. Rev. Phys. Chem. 28, 169 (2009). 8G. D. Fletcher (unpublished). 9T. Thorsteinsson, D. L. Cooper, J. Gerratt, P. B. Karadov, and M. Raimondi, Theor. Chim. Acta 93, 343 (1996). 10G. E. Scuseria, C. A. Jimenez-Hoyos, T. M. Henderson, K. Samanta, and J. K. Ellis, J. Chem. Phys. 135, 124108 (2011). 11D. D. Kemp, J. M. Rintelman, M. S. Gordon, and J. H. Jensen, Theor. Chem. Acc. 125, 481 (2010). 12G. D. Fletcher, G. Doggett, and A. S. Howard, Phys. Rev. A 46, 5459 (1992). 13G. Doggett and G. D. Fletcher, J. Mol. Struct.: THEOCHEM 260, 313 (1992); G. Doggett, G. D. Fletcher, and F. R. Manby, J. Mol. Struct. 300, 191 (1992); F. R. Manby, G. Doggett, and G. D. Fletcher, J. Mol. Struct.: THEOCHEM 343, 63 (1995). 14G. D. Fletcher, Ph.D. thesis, University of York, UK, 1993. 15R. Pauncz, Spin Eigenfunctions (Plenum, New York, 1979). 16T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 17M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. J. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 18J. J. P. Stewart, J. Comput. Chem. 10, 209 (1989). 19B. M. Bode and M. S. Gordon, J. Mol. Graphics Modell. 16, 133 (1999). 20M. Dupuis, J. Rys, and H. F. King, J. Chem. Phys. 65, 111 (1976). 21B. S. Garbow, Comput. Phys. Commun. 7, 179 (1974). 22J. Dongarra, ACM Trans. Math. Software 10, 219 (1984). 23G. J. H. van Nes and A. Vos, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 34, 1947 (1978). 24See supplementary material at http://dx.doi.org/10.1063/1.4916743 for the atomic coordinates of the structure. 2D.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:11:11
Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.
