| |
|
|
|
This chapter contains a description of some of the theory behind the methods used in Jaguar. Section 7.1 describes the pseudospectral method itself. Section 7.2, Section 7.3, and Section 7.4 describe GVB, GVB-RCI, and LMP2 calculations and how the pseudospectral method improves computational scaling and efficiency for these methods. Section 7.5 contains a brief description of density functional theory. Chapter 3 includes information about performing Jaguar calculations using the techniques described here.
Like conventional ab initio electronic structure codes, Jaguar solves the Schrödinger equation iteratively, using self-consistent field methods to calculate the lowest-energy wavefunction within the space spanned by the selected basis set. For calculations on large molecules, both conventional and pseudospectral techniques must recalculate key integral terms for each SCF iteration, since storage costs for these terms are prohibitive.
However, most of the fundamental integrals calculated in the pseudospectral method [19] are computed in physical space, on a grid, rather than in the spectral space defined by the basis functions. The pseudospectral method takes the density matrix from the wavefunction guess used at the beginning of each SCF iteration and the values of the integrals upon the grid points and manipulates them to produce the necessary operators on the grid, then assembles the Fock matrix by transforming these components back into spectral space, where the Fock matrix is employed in the usual way to generate the wavefunction for the next iteration.
For medium and large molecules, the additional overhead the pseudospectral method requires to compute the information needed for the transformation between physical and spectral space is vastly outweighed by the advantages of evaluating the integrals in physical space. The matrix needed for the transform from physical to basis set space [7] can be assembled before the SCF iterations by calculating the least-squares operator Q, which is given by the equation
where S is the usual analytic overlap matrix, R is a matrix of fitting functions evaluated at the grid points, and w is a diagonal matrix of grid weights. The fitting functions used to construct the matrix R include both basis functions and dealiasing functions, which are chosen in order to span the function space represented by the grid more completely than the basis functions alone. The operator Q can be calculated for the relevant basis functions using several different sets of grid points, where each set of points defines a grid type, ranging from coarse to ultrafine depending on the number of grid points.
In practice, not all possible
elements are calculated for each basis function i and each grid point g, because most basis functions drop off sharply enough that they have no significant value on some or most grid points. These functions are classified as short-range functions and are grouped together by atom, while the remaining functions are classified as long-range functions, which are all considered to be in one single group [13].
Since Q does not depend on the wavefunction itself, it can be fully computed before the SCF procedure. However, since the Q for each grid type contains Nbasis x Ngrid elements, where Nbasis is the number of basis functions and Ngrid the number of grid points (which is generally larger than Nbasis), we sometimes reduce memory demands by only computing and storing the Nbasis x Nfit matrix
in the program rwr, for cases where the Q for that grid type is only needed for one SCF iteration. We then assemble the full Q during the SCF iteration for which it is needed.
After the Jaguar program rwr has generated the Q or
matrix, the program scf takes the initial orbitals and iteratively modifies them with the pseudospectral method until convergence. This process involves calculating the values of the necessary integrals on the grid points, and actually assembling the Fock matrix from the computed information. The three-center, one-electron pseudospectral integrals on the grid points are defined by
where
and
are basis functions and the index g represents a grid point. These integrals are calculated for all combinations of basis functions and grid points not eliminated by cutoffs, and the Fock matrix is assembled from its Coulomb and exchange matrix components
and
, which are calculated in physical space and transformed back into spectral space by the following equations:
where D is the usual spectral space density matrix,
is the value of the function j at grid point g, and
is given by Equation (2). The grid points used for each SCF iteration are determined by the grid type (coarse, medium, fine, or ultrafine) chosen for that iteration. The number of arithmetic operations involved in the assembly of the matrices J and K in Equation (3a) and Equation (3b) scales formally as N3, as opposed to the N4 scaling for the matrix assembly in the conventional spectral space algorithm.
Jaguar actually uses the pseudospectral method described above for the majority of the computationally intensive two-electron integral terms, but calculates the one-electron and some of the largest and most efficiently computed two-electron terms analytically [13]. More specifically, for the Coulomb matrix elements, we calculate the analytic terms
for cases in which i, j, k, and l meet certain cutoff criteria and the term (ij | kl) is of the form (aa | aa), (aa | ab), (aa | bb), (ab | ab), or (aa | bc), where a, b, and c indicate the atom upon which the function is centered. Similar correction terms are computed for the exchange operator, as detailed in ref. [13]. The corresponding pseudospectral terms, as dictated by Equation (3a) and Equation (3b) for the appropriate i, j, k, and l choices, must be subtracted from the pseudospectral J and K elements as well. This combined pseudospectral/analytic approach allows Jaguar to take advantage of the strengths of both methods, since it can largely maintain the pseudospectral method speedups for a particular grid, and can also use a coarser grid than a purely numerical calculation would allow.
The pseudospectral method has also been extended to electron correlation methods, with a particular focus on Generalized Valence Bond (GVB) [20] calculations. Highly refined GVB initial guess [14] and convergence [11] algorithms have been automated within Jaguar, allowing the scaling advantages resulting from the pseudospectral method to be maintained for GVB calculations. The method yields very accurate excitation energies, rotational barriers, and bond energies for many molecules, and GVB calculations with Jaguar are typically 10 to 100 times more efficient than the best conventional GVB programs, even for molecules as small as ten atoms [6].
In the GVB approach, each bond or other electron pair is described by two non-orthogonal orbitals, whose contributions to the bond description are obtained variationally. The bond description can thus change smoothly from a description with two atomic-like orbitals at large bond distances to a description with bond-like orbitals at short distances. This improvement over Hartree-Fock, which treats bonds as having equal amounts of covalent and ionic character, allows GVB to describe charge transfer reactions and bond breaking and formation accurately, and also gives better results for other molecular properties than an HF treatment alone can provide.
The goal of a GVB calculation, then, is to obtain pairs of GVB orbitals ypa and ypb, where p ranges from 1 to the number of GVB pairs Ngvb, that lead to a minimum energy for the molecular wavefunction
For a given p, the orbitals ypa and ypb form a pair that describes a particular bond or other pair of electrons. Under the perfect pairing restriction, the GVB orbitals within a pair are not orthogonal, although they are each orthogonal to all GVB orbitals in other pairs. For computational purposes, it is useful to form orthogonal GVB natural orbitals ypg and ypu from the GVB orbitals ypa and ypb and their overlap Sp, as follows:
The ypg orbitals generally have a bonding character, while the ypu orbitals are anti-bonding. The contribution to the GVB wavefunction from each pair is given by
where the GVB configuration interaction (CI) coefficients Cpg and Cpu satisfy the following equations:
Solving for the optimal GVB orbitals is therefore a matter of determining both the GVB natural orbitals and the GVB CI coefficients that minimize the energy of the GVB wavefunction. This energy is given by the equation
where m and n range over all GVB natural orbitals (bonding and anti-bonding), and where these orbitals are expanded in terms of the basis functions, as shown here:
The terms hmm, Jmn, and Kmn are defined by:
and the quantities amn and bmn obey the following rules:
Examining the variation of the energy E with respect to the basis set coefficients c gives the equations for the Fock operator corresponding to each GVB natural orbital:
Each orbital's Fock operator thus depends on the other orbitals' Coulomb and exchange operators.
At the beginning of each SCF iteration, the scf program is provided with a set of proposed natural orbitals and a set of CI coefficients that dictate the contribution of each natural orbital to the GVB orbitals. For that set of GVB natural orbitals, the program first solves for revised CI coefficients by evaluating the Coulomb and exchange matrix elements for those orbitals and diagonalizing the two-by-two operators Yp in the basis of the two natural orbitals in pair p, as described by these equations:
In practice, since the CI coefficients mutually interdependent, they are determined using a self-consistent iterative procedure.
Next, holding the CI coefficients fixed, the program evaluates the energy and the Fock matrix and adjusts the basis set coefficients describing the GVB natural orbitals accordingly, in basically the same manner used for the usual HF treatment. The revised orbitals and CI coefficients are then used in the next SCF iteration, and the process continues until both the GVB natural orbitals and the CI coefficients have converged.
The GVB treatment can also be applied to open shell cases, or restricted to certain electron pairs. These variations are described in reference [20], which also provides much more detail about the GVB methods and equations. The ability to restrict the use of GVB to particular electron pairs is an important strength of the method. This feature allows computationally inexpensive correlation of critical regions in very large molecules.
The GVB-RCI (restricted configuration interaction) wavefunction is the simplest multideterminantal reference wavefunction which properly dissociates to open shell fragments regardless of the spin multiplicity of the fragments. A critical advantage of GVB-RCI is that the GVB and RCI computations can be confined to a localized region of the molecule. The GVB-RCI method is therefore particularly useful for evaluating bond energies and bond formation and breaking, as well as for studies of open shell radicals and other systems for which it is important to avoid spin contamination problems.
The version of GVB-RCI within Jaguar uses pseudospectral numerical methods and a novel internal contraction scheme in which a GVB-PP wavefunction is used as a correlated mean field reference state [12]. This implementation of GVB-RCI can be used to generate highly accurate GVB-RCI wavefunctions, with energies within about 0.1 kcal/mole of results from all-analytical integral calculations [12]. The internal contraction scheme used restricts the number of CI coefficients in the RCI calculation to ~n3, where n is the number of GVB pairs, yet is in excellent agreement with a fully uncontracted CI which by contrast would contain 2nn3 CI coefficients (the number of uncontracted determinants).
The GVB-RCI program within Jaguar generates a correlated wavefunction from intra-pair excitations of the GVB reference wavefunction described in Section 7.2, using a highly effective contraction procedure to reduce the length of the CI expansions. The program employs the pseudospectral method to speed up integral evaluation, and systematically includes the most important configurations to make the calculation more practical, with minimal loss of accuracy relative to the fully uncontracted expansion.
The spatial states for an RCI pair are constructed from the same natural orbitals as those used for the GVB reference wavefunction, ypg and ypu , but in addition to the GVB spatial state from Equation (6), rewritten here:
the RCI spatial states include the orthogonal complements xp1 and xp2:
Just as the GVB method allows the user to correlate particular electron pairs for maximal efficiency, the RCI treatment can be applied to any user-specified subset of the GVB pairs. We then use a GVB mean field procedure to evaluate a Coulomb-exchange mean field operator describing the effect of the non-excited GVB pairs on the RCI pairs. This treatment effectively reduces the two-electron part of the Hamiltonian to the space of the RCI coordinates. Even for cases with many RCI pairs, we restrict the configurations to those with only a small number of excitations and use the mean field treatment for each configuration's calculation.
The RCI spatial states xp1 add an extra complication to the necessary evaluation of Coulomb and exchange matrix elements using the natural orbitals ypg and ypu . For the GVB case, it is sufficient to compute the following matrix elements (corresponding to Equation (10b) and Equation (10c)):
where the m and n can each be any natural orbital. For the RCI pairs, on the other hand, we need all matrix elements of the form:
where the a and b natural orbitals are from the same RCI pair p (and may be the same natural orbital), while the g and d natural orbitals are from the same RCI pair with index q. The complicated part of the calculation of the Coulomb and exchange operators, then, is evaluating matrix elements in atomic orbital (AO) space and using the AO-space matrix elements to produce the matrix elements in the natural orbital space, a process that normally requires a four-index transform.
By using the pseudospectral method, however, Jaguar reduces the scaling of the evaluation of each Coulomb or exchange matrix operator in basis function space from N4 to N3, and solves for the necessary matrix elements with a two-index transform rather than an expensive four-index one. For simplicity, we will describe this process for the Coulomb matrix elements only; the equations for K are similar. We first evaluate the usual three-center, one-electron integrals Aklg in basis function space (see Equation (2)). We then evaluate Jgd in physical space for all gd corresponding to orbital products of each RCI pair, ypgypg , ypgypu , and ypuypu , using the equation
Next, we transform back into spectral space to get
, where i and j are basis functions, using the pseudospectral method in the usual manner described in Section 7.1:
where Q is the pseudospectral least-squares operator and Rjg is the value of the basis function j at grid point g. We perform a final two-index transform,
to obtain the matrix elements in the natural orbital basis.
When Jaguar has obtained all Coulomb and exchange operators, it performs an iterative diagonalization of the Hamiltonian to obtain the RCI coefficients. We employ the Davidson method for this step.
Second order Møller-Plesset perturbation theory (MP2) is perhaps the most widely used ab initio electron correlation methodology, recovering a large fraction of the correlation energy at a relatively low computational cost. The method greatly improves Hartree-Fock treatments of properties such as transition states, dispersion interactions, hydrogen bonding, and conformational energies. However, the scaling of conventional MP2 algorithms with system size is formally nN4, where N is the number of basis functions and n the number of occupied orbitals, due to the necessity of carrying out a four index transform from atomic basis functions to molecular orbitals. In principle, it is possible to reduce this scaling by using integral cutoffs, as for Hartree-Fock calculations. However, the reduction is noticeably less effective in MP2, particularly for the large, correlation-consistent basis sets that are required for accurate correlation effects on observable quantities. Thus, MP2 techniques have traditionally been used primarily for small molecules.
Several years ago, Pulay and coworkers [35, 36] formulated a version of MP2 in which the occupied orbitals are first localized (e.g., via Boys localization [38]) and the virtual space correlating such orbitals are then truncated to a local space, built from the atomic basis functions on the local atomic centers orthogonalized to the occupied space. Another critical advantage of LMP2 (as for other localized correlation methods such as GVB and GVB-RCI) is that one can very precisely control which region of the molecule is correlated, reducing CPU costs enormously. The method has been shown to yield an accuracy for relative energies that is, if anything, superior to conventional MP2, due to elimination of basis set superposition error [37]. However, localized MP2 implementations in conventional electronic structure codes have not yet led to substantial reductions in CPU time, since the first few steps of the necessary four-index transform are unaffected by localization of the occupied orbitals, and the localized orbitals have tails that extend throughout the molecule.
We have carried out extensive tests demonstrating the accuracy and computational efficiency of the pseudospectral implementation of LMP2, as detailed in ref. [16]. In the pseudospectral approach, we assemble two-electron integrals over molecular orbitals directly and are thus able to fully profit from the huge reduction in the size of the virtual space in Pulay's theory. Formally, the PS implementation of LMP2 scales as nN3; however, various types of cutoffs and multigrid procedures can reduce this to ~N2. In fact, for calculations involving both the 631G** and Dunning ccpVTZ basis sets, we find a scaling ~N2.7 with system size.
The physical idea behind the LMP2 method is that if the molecular orbitals are transformed so that they are localized on bonds or electron pairs, correlation among the occupied pairs can be described by the local orbital pairs and their respective local pair virtual spaces defined from the atomic orbitals on the relevant atom or pair of atoms. The localized orbitals can be generated by any unitary transformation of the canonical orbitals. For LMP2, we use Boys-localized [38] orbitals, for which the term
is maximized. The local virtual space for each atom is defined by orthogonalizing its atomic basis functions against the localized molecular orbitals. The correlating orbitals included in the local virtual space are thus mostly near the atom itself, but because of the orthogonalization procedure, they are not particularly well localized.
The Jaguar LMP2 program uses Pulay's method [35, 36, 37] to expand the first order wavefunction correction
as a linear combination of determinants formed by exciting electrons from localized orbitals i and j to local virtual space correlation orbitals p and q:
For local MP2, we must iteratively solve the following equation, which has been derived in detail by Pulay and Sæbo, for the coefficients
:
Here F is the Fock matrix and S is the overlap matrix. The exchange matrix Kij is restricted to the dimensions of the virtual space corresponding to the occupied localized molecular orbitals i and j. The simplest updating scheme for the Cij coefficients is to obtain updated coefficients
iteratively from the equation:
where ei and ej are the matrix elements Fii and Fjj in the localized molecular orbital basis and ep and eq are the eigenvalues of the Fock matrix in the local virtual basis.
From the Cij coefficients and the exchange matrices Kij, Jaguar computes the second order energy correction
from the equations:
where the bracket in Equation (21a) denotes a trace and dij is 1 if i = j and 0 otherwise. Computing the exchange matrix elements for Equation (21a) is approximately 80% of the work for an energy correction computation, while generating the Cij coefficients comprises about 20% of the work.
Jaguar performs localized MP2 calculations using pseudospectral methods, evaluating integrals over grid points in physical space in a manner similar to that described for HF and GVB calculations in Section 7.1 and Section 7.2. The two-electron exchange integrals needed for Equation (21a) are evaluated over grid points g as follows:
where Qig is the least squares fitting operator for molecular orbital i on grid point g, Rpg is the physical space representation of virtual orbital p, and Ajqg is the three-center, one-electron integral over the occupied molecular orbital j and the local virtual orbital q. The last term is related to the three-center, one-electron integrals in atomic orbital space, Aklg, described in Equation (2), by
Therefore, obtaining the exchange integral values requires solving for the integrals Aklg in atomic orbital space; transforming these as:
performing the second transformation step to yield the integrals in molecular orbital space:
and forming the necessary
elements using Equation (22). Jaguar's local MP2 module also includes analytical corrections similar to those described earlier for Hartree-Fock and GVB calculations, and a length scales algorithm, both of which are explained in ref. [13].
Density functional theory (DFT) is based on the Hohenberg-Kohn theorem [92], which states that the exact energy of a system can be expressed as a functional depending only on the electron density. In the Kohn-Sham implementation of DFT [93], this density is expressed in terms of Kohn-Sham orbitals {yi}:
similarly to the density expression used for Hartree-Fock SCF calculations. (For simplicity, we consider only closed shell systems in this overview of the method.)
The Kohn-Sham orbitals are expressed as a linear combination of basis functions
, and the coefficients for this expansion are solved iteratively using a self-consistent field method, as for Hartree-Fock. However, DFT includes exchange and/or correlation density functionals within the Fock matrix used for the SCF procedure. For DFT calculations, the Hartree-Fock exchange term Kij in the Fock matrix is replaced by the exchange-correlation potential matrix elements
:
where
is an exchange-correlation functional and g is
.
The exchange-correlation functional
is usually separated into exchange and correlation functional components that are local or non-local in the density:
Under the local density approximation (LDA), the non-local functionals
and
are ignored; when either or both of these terms are included, the generalized gradient approximation (GGA), also known as the non-local density approximation (NLDA), applies. The local and non-local exchange and correlation functionals available within Jaguar are described in section 3.1 and its references.
The electronic ground state energy E0 is given by
(in Hartree atomic units), where Vnuc is the nuclear potential and J is the Coulomb potential. Therefore, for a given exchange-correlation functional, it is possible to solve iteratively for Kohn-Sham orbitals
and the resulting density r to yield a final DFT energy.
A more detailed description of density functional theory can be found in references [94] and [95].
|
Schrödinger, Inc. http://www.schrodinger.com Voice: (503) 299-1150 Fax: (503) 299-4532 help@schrodinger.com Last updated: Thu Oct 11 19:10:38 2001 |
| |
|
|
|
.
