## Density functional theory (DFT)

Ground-state properties are functions of the charge density

Density function theory begins with a theorem by Hohenberg and Kohn (1964, later generalized by Levy 1979), which states that all ground-state properties are functions of the charge density . Specifically, the total energy Et may be written as:

#### Eq. 1

where T [] is the kinetic energy of a system of noninteracting particles of density , U [] is the classical electrostatic energy due to Coulombic interactions, and Exc [] includes all many-body contributions to the total energy, in particular the exchange and correlation energies. Eq. 1 is written to emphasize the explicit dependence of these quantities on (in subsequent equations, this dependence is not always indicated).

As in other molecular orbital methods (Roothaan 1951, Slater 1972, Dewar 1983), the wavefunction is taken to be an antisymmetrized product (Slater determinant) of one-particle functions, that is, molecular orbitals (MOs):

#### Eq. 2

The molecular orbitals also must be orthonormal:

#### Eq. 3

The charge density summed over all molecular orbitals

In this case, the charge density is given by the simple sum:

#### Eq. 4

Spin-restricted and -unrestricted calculations

where the sum goes over all occupied MOs i. The density obtained from this expression is also known as the charge density. The MOs may be occupied by spin-up (alpha) electrons or by spin-down (beta) electrons. Using the same i for both alpha and beta electrons is known as a spin-restricted calculation; using different i for alpha and beta electrons results in a spin-unrestricted or spin-polarized calculation. In the unrestricted case, it is possible to form two different charge densities: one for alpha MOs and one for beta MOs. Their sum gives the total charge density and their difference gives the spin density--the amount of excess alpha over beta spin. This is analogous to restricted and unrestricted Hartree-Fock calculations (Pople and Nesbet 1954).

From the wavefunctions (Eq. 1) and the charge density (Eq. 4), the energy components can be written (in atomic units) as:

#### Eq. 6

In Eq. 6, Z refers to the charge on nucleus of an N-atom system. The first term, VN , represents the electron-nucleus attraction. The second, Ve/2, represents the electron-electron repulsion. The final term, VNN , represents the nucleus-nucleus repulsion.

The final term in Eq. 1, the exchange-correlation energy, requires some approximation for this method to be computationally tractable. A simple and surprisingly good approximation is the local density approximation, which is based on the known exchange-correlation energy of the uniform electron gas (Hedin and Lundqvist 1971, Ceperley and Alder 1980, Lundqvist and March 1983). Analytical representations have been made by several researchers (Hedin and Lundqvist 1971, Ceperley and Alder 1980, von Barth and Hedin 1972, Vosko et al. 1980, Perdew and Wang 1992). The local density approximation assumes that the charge density varies slowly on an atomic scale (i.e., each region of a molecule actually looks like a uniform electron gas). The total exchange-correlation energy can be obtained by integrating the uniform electron gas result:

#### Eq. 7

where Exc [] is the exchange-correlation energy per particle in a uniform electron gas and is the number of particles.

This implementation uses the form derived by Vosko et al. (1980), denoted the VWN functional, as the default. Other local spin-density functionals, developed by Von Barth and Hedin (BH), Janak, Morruzi, and Williams (JMW), and the latest Perdew and Wang work (PW), are also available when a nonstandard type of DMol calculation is chosen.

The next step in improving the local spin-density (LSD) model is to take into account the inhomogeneity of the electron gas which naturally occurs in any molecular system. This can be accomplished by a density gradient expansion, sometimes referred as the nonlocal spin-density approximation (NLSD). Over the past few years it has been well documented that the gradient corrected exchange-correlation energy Exc[, d ()] is necessary for studying the thermochemistry of molecular processes (see recent reviews by Ziegler 1991, Labanowski and Andzelm 1991, Politzer and Seminario 1995).

In the present implementation, several commonly used NLSD functionals are available. The default is the Perdew and Wang (PW) generalized gradient approximation for the correlation functional and the Becke (B) gradient corrected exchange functional. This level of approximation for NLSD Hamiltonian is referred to as the BP functional. In addition, the gradient corrected correlation functional of Lee, Yang, and Parr (LYP) is available.

The total energy can now be written as:

#### Eq. 8

To determine the actual energy, variations in Et must be optimized with respect to variations in , subject to the orthonormality constraints in Eq. 2 (Kohn and Sham 1965):

#### Eq. 9

Kohn-Sham equations

where ij are Lagrangian multipliers. This process leads to a set of coupled equations first proposed by Kohn and Sham (1965):

#### Eq. 10

The term µxc is the exchange-correlation potential, which results from differentiating Exc. For the local-spin density approximation, the potential µxc is:

#### Eq. 11

Use of the eigenvalues of Eq. 10 leads to a reformulation of the energy expression:

#### Eq. 12

Convenience of expanding MOs in terms of AOs, or basis functions

In practice, it is convenient to expand the MOs in terms of atomic orbitals (AOs):

#### Eq. 13

The atomic orbitals µ are called the atomic basis functions, and the C are the MO expansion coefficients. Because linear combinations of AOs are used, the C are also sometimes called the LCAO-MO coefficients. Several choices are possible for the basis set, including Gaussian (Andzelm et al. 1989), Slater (Versluis and Ziegler 1988), and numerical orbitals. In this implementation, numerical orbitals are used--these are discussed at some length under Numerical basis sets.

where:

and:

#### Eq. 16

SCF procedure

Because H depends upon C, Eq. 14 must be solved by an iterative technique. This can be done by the following procedure:

For an organic molecule, about ten iterations are typically required to obtain convergence at new - old < 10-6. For metallic systems, many more iterations are frequently required.

### Numerical basis sets

Atomic basis sets are generated numerically

The basis functions m are given numerically as values on an atomic-centered spherical-polar mesh, rather than as analytical functions (i.e., Gaussian orbitals). The angular portion of each function is the appropriate spherical harmonic Ylm(,). The radial portion F(r) is obtained by solving the atomic DFT equations numerically. A reasonable level of accuracy is usually obtained by using a range of 300 radial points from the nucleus to an outer distance of 10 Bohrs ( 5.3 Å). Radial functions are stored as a set of cubic spline coefficients for each of the 300 sections, so that F(r) is actually piecewise analytic. This is an important consideration for generating analytic energy gradients. In addition to the basis sets, the -2/ 2 terms required for evaluation of the kinetic energy are also stored as spline coefficients.

The use of the exact DFT spherical atomic orbitals has several advantages. For one, the molecule can be dissociated exactly to its constituent atoms (within the DFT context). Because of the quality of these orbitals, basis set superposition effects (Delley 1990) are minimized, and an excellent description of even weak bonds is possible.

Greater variational freedom is obtained by providing larger basis sets. Generation of an entire second set of functions results in doubling the basis set size; this is referred to as a double-numerical (DN) set. This notation is used by analogy with Gaussian double-zeta (DZ) sets, but the N is used to emphasize the numerical nature of these orbitals. Additional basis functions, including polarization, are obtained by several procedures:

For first-row atoms, functions from +2 ions provide a reasonable double basis set. A hydrogenic 3d orbital obtained for a nucleus of Z = 5 provides a good polarization function for each of these atoms. A hydrogenic 2p function for Z = 1.3 is used for hydrogen. The use of various nuclear charges to generate polarization functions is analogous to the variation of zeta in Gaussian basis sets. For metals, 4p polarization functions are generated by solving the atomic equations for a 4s 4p excited state. Basis set quality has been analyzed in detail by Delley (1990).

### Numerical integration

Evaluation of the integrals in Eqs. 15 and 16 must be accomplished by a 3D numerical integration procedure, due to the nature of the basis functions. The matrix elements need to be approximated by the finite sums:

#### Eq. 18

The sums run over several numerical integration points ri. The term Heff (ri) represents the value of the integrand of Eq. 15 at point ri and w(ri) represents a weight associated with each mesh point. Increasing the number of mesh points improves the numerical precision of the integral but also results in additional computational cost.

Careful selection of a set of integration points is important for the quality of the calculation (Delley 1990, Ellis and Painter 1968, Boerrigter et al. 1988). In general, the grid used to generate the atomic basis set is not suitable for molecular calculations. The grid used for atomic basis sets can take advantage of spherical symmetry, which greatly simplifies matters. For molecules, it is necessary to be able to treat the rapid oscillations of the molecular orbitals near the nuclei and to avoid integration of the nuclei themselves because of the nuclear cusps (Delley 1990).

The integration points are generated in a spherical pattern around each atomic center. Radial points are typically taken from the nucleus to an outer distance of 10 Bohrs ( 5.3 Å). The number of radial points within this distance is designed to scale with increasing atomic number. For example, Fe requires more points than C. The typical number of radial points NR for a nucleus of charge Z is:

#### Eq. 19

This number may, of course, be manually adjusted to accommodate the required precision or allowed cost of a calculation. The spacing between points is logarithmic--points are spaced more closely near the nucleus where oscillations in the wavefunction are more rapid.

Angular integration points N are generated at each of the NR radial points, creating a series of shells around each nucleus. Angular points are selected by schemes designed to yield points ri and weights w(ri), which could yield exact angular integration for a spherical harmonic with a given value of l. Such quadrature schemes (Stroud 1971, Lebedev 1975 and 1977, Konyaev 1979) are available for functions up to l = 41. Alternatively, a product-Gauss rule in cos and may be used for arbitrary values of l (Stroud 1971). The product-Gauss methods use (l + 1) points on each shell, while the quadrature methods use more points:

l N
5   14
7   26
11   50
17   110
23   194
29   302
35   434
41   590

In practice, the quadrature method is used to generate the integration grid, and the product method is used as a check of numerical precision. On each radial shell, the sum of the atomic densities and the atomic electrostatic potentials is computed using both angular schemes. If the difference between the results of the two schemes is too large, then the next highest l value is used to generate a set of new points. This assures a consistent level of numerical precision throughout the integration grid. Typically, angular sampling increases when the difference between the two schemes is greater than 10-4. This yields 1000 points per atom and offers a good compromise between computational cost and numerical precision.

Partition functions are used to increase the convergence of the numerical integration and to avoid integrating over nuclear cusps (Delley 1990, Hirshfeld 1977, Becke 1988). A partition function a is defined as:

#### Eq. 20

where is an atom index and g(r - R) is a function that typically is large for small r - R and small for large r - R (i.e., larger near the nucleus). Integrals are rewritten using partition functions as:

#### Eq. 21

which is further reduced to a sum over 3D integration points:

#### Eq. 22

Preferred partition function

In practice, the partition functions are combined with the integration weights w(ri) to simplify the computation. The preferred choice for a partition function is:

#### Eq. 23

where r = |ri - R|, r0 = 0.5, and a is the atomic charge density for atom . Other options for partition functions include:

### Evaluation of the effective potential

Evaluation of the exchange-correlation energy xc and potential µxc is accomplished with Eqs. 52, 55, and 57. This requires numerical evaluation of the charge density (r) at many points in space (i.e., xc and µxc are tabulated numerically). This restriction actually applies to most density function methods, even if analytical basis functions are used (Andzelm et al. 1989, Versluis and Ziegler 1988). The use of numerical basis functions facilitates this process, since all required quantities are already available on a grid of adequate numerical precision. An alternative approach (Baerends et al. 1973) is to fit the charge density to an analytic multipolar expansion via a least-squares fitting procedure. This simplifies the evaluation of xc and µxc, but still requires the use of a numerical grid for the least-squares fitting.

The Coulombic potential is evaluated by solving the Poisson equation for the charge density:

#### Eq. 28

rather than by explicitly evaluating the Coulombic term as:

#### Eq. 29

In the this approach, the Poisson equation is solved in a completely numerical (non-basis set) approach (Delley 1990). This provides greater numerical precision, since the evaluation of Ve is essentially exact once the form of (r) has been specified. Such a method requires specification of an analytic form of (r), as discussed above. However, rather than use a least-squares fitting procedure, a projection scheme is used. The charge density is first partitioned into atomic densities and then decomposed into multipolar components. Appropriate partition functions can ensure that such an expansion is rapidly convergent.

The density obtained in this way is called the model density. The term alm r - R gives the model density for the multipolar component lm on atom for a shell at r - R distance from the nucleus:

#### Eq. 30

Note that the partition function p used for decomposition of the density is in general not the same as that used to improve the numerical integration in Eq. 22. The total model density is obtained from the summation over all alm:

#### Eq. 31

Effect of angular truncation on precision of model charge density

The total model charge density is, in general, not equal to the orbital density because of angular truncations. However, the flexibility of this model charge density is superior to that obtained with fitting procedures. The degree of angular truncation can be specified as an input parameter. Typically, a value of l one greater than that in the atomic basis provides sufficient precision; for example, the use of l = 3 truncation when p functions are present in the basis or l = 2 truncation if only p functions are used.

The Coulombic potential for each component is calculated using the Green function of the Laplacian (Delley 1990):

#### Eq. 32

The total potential

and the total potential is given by:

### Computational self-consistent field procedure

Interpolating the numerical atomic bases onto the molecular grid

Before the self-consistent field (SCF) procedure can be used, a step is required that is analogous to the evaluation of integrals over atomic orbitals, such as in Hartree-Fock methods. This is the interpolation of the numerical atomic basis set onto the molecular grid. Neglecting symmetry and any frozen-core approximations, this step requires a computational effort on the order of N X P for N atomic orbitals and P integration points. The basis set is controlled by the user, the number of points typically being on the order of 1000 points per atom. The overlap matrix (Eq. 16) and the constant portion of the effective Hamiltonian (Eq. 15) (kinetic and nuclear attraction terms) are constructed at this time, and each requires N X (N + 1) X P operations. The interpolated values can be stored externally and read as they are required. Alternatively, these data can be generated as needed, obviating the need for storage--this is termed a direct SCF procedure, by analogy with the direct Hartree-Fock method (Almlöf et al. 1982).

In practice, it is more convenient to skip Steps 1 and 2 (choosing an initial C and constructing an initial set of i, (see Choose an initial set of Cim.)) and to begin instead with an initial constructed from the superposition of atomic densities (quantities that are readily constructed from the numerical atomic basis set). In the SCF procedure, reconstruction of the new density requires a computational effort on the order of N X No X P, where No is the number of occupied orbitals.

Once (r) is known, evaluation of the exchange-correlation potential µxc requires only P operations. Construction of the Coulombic potential requires only P X M effort, where M is the number of multipolar functions. M is typically on the order of 9 functions atom-1 (l = 2) or 16 functions atom-1 (l = 3). Neither of these steps is especially time consuming.

Construction of the Hamiltonian matrix elements (Eq. 14) is among the most time-consuming steps, requiring N X (N + 1) X P operations each iteration. Solution of the secular equation is also time consuming, requiring N3 operations. This can be reduced by solving only for the eigenvectors corresponding to occupied orbitals to N2 X No.

For large systems, the computational cost does not necessarily grow as rapidly as implied by the above comments. Since the atomic basis functions have a finite extent (10 Bohrs), only a limited number of points contribute to each matrix element and P eventually converges to a constant. In addition, construction of the density and the secular matrix can be accomplished using sparse matrix multiplier routines, further reducing the computational cost.

Construction of a new density follows solution of the secular equation. Damping is usually required to ensure smooth convergence. In the current method, simple damping is possible:

#### Eq. 34

where d is the damping factor, old is the density that was used to construct the secular matrix, new is the density constructed from the new MO coefficients without damping, and is the density that is actually used in the next iteration. An interpolation/extrapolation scheme is available. This technique constructs an effective vector from old to new for the current iteration and for the previous iteration. The effective vectors are generally skew vectors. The point of closest approach between old and new is used to extrapolate the actual density for the next iteration. The Pulay DIIS method (Pulay 1982) has been implemented.

The direct inversion of the iterative subspace (or DIIS) method developed by Pulay (Pulay, 1982) has been implemented in DMol as a mechanism to speed up SCF convergence. The method rests on the suitable definition of an error vector that is zero when convergence is achieved and on performing a linear combination of a set of error vectors sampled along iterations that produces a new error vector with minimal norm. (This algorithm was present in previous versions of DMol, yet restricted to interpolation of at most two error vectors.) The DIIS method is much more powerful if one allows for a slightly larger dimension of iterative subspace. The default dimension currently adopted in DMol is 4. An upper limit of 10 was imposed in the current DIIS scheme, which is motivated by the fact that linear dependencies between error vectors and therefore the amount of redundant information increases with the size of the vector space.

The error vector for the DIIS procedure is defined as the difference between the input radial densities and the projected output radial densities (see Eq. 30 for the definition of the model density):

#### Eq. 35

The final model density is a linear combination of the densities at each SCF iteration i:

#### Eq. 36

with the constraint that . The Ci coefficients are calculated by minimizing the norm:

where .

#### Eq. 38

Effect of auxiliary density approximation on accuracy of calculated total energy

To minimize the impact of this difference a total-energy formula is used that is second-order in (Delley et al. 1983, Delley 1990, Delley 1991). The default starting density in DMol is the sum of the spherical atom densities. The total energy calculated in the first SCF cycle is thus the so-called Harris approximation (Harris 1985). The atomic dissociation energy in the first cycle is usually overestimated, since the electrostatic error term for the total energy:

#### Eq. 39

is negative definite. The less important second-order term for local density functionals is positive definite, which leads to a slight overestimation of the total energy during the SCF iterations.

#### Eq. 40

where is the density for spin alpha and beta, respectively; ni are the occupations of orbitals with the orbital and spin labels i, ; i is the corresponding eigenvalue, which has been calculated using the static V~e and exchange-correlation potential µ~xc, arising from the spin densities ~. Exc is the exchange-correlation functional (local or nonlocal).

### Energy gradients

Predicting chemical structure

The ability to evaluate the derivative of the total energy with respect to geometric changes is critical for the study of chemical systems. Without the first derivatives, a laborious point-by-point procedure is required, which is taxing to both computer and human resources. The availability of analytic energy derivatives for Hartree-Fock (Pulay 1969), CI (Brooks et al. 1980), and MBPT (Pople et al. 1979) theories (to name just a few) has made these remarkably successful methods for predicting chemical structures.

The derivative of the total energy in Eq. 12 with respect to a nuclear perturbation in direction a (x, y, or z) of atom may be written as:

#### Eq. 41

where the derivative density is defined as:

with:

#### Eq. 43

Derivative of the basis function

The derivative of the basis function m with respect to the perturbation a can be computed analytically because of the representation of the numerical basis sets. The angular portion of m is a spherical harmonic function, which is analytic and easily differentiated. The radial portion is represented by several cubic splines, each of which is also differentiable.

The derivatives of the eigenvalues can be obtained from Eq. 10. Multiplying by i and integrating gives an expression for i:

#### Eq. 44

Differentiating and rearranging yields:

#### Eq. 45

The terms in Eqs. 45 and 49 involving Z represent the Hellmann-Feynman force (Hellmann 1939, Feynman 1939), which gives the derivative in the absence of any orbital relaxation.

Substituting Eq. 45 into Eq. 41 yields:

#### Eq. 46

Where Eat is the Hellmann-Feynman term. Now (Andzelm et al. 1989):

#### Eq. 47

and recalling the definition in Eq. 11, we can also write:

#### Eq. 48

The final equation for the derivative of the energy

Therefore, the terms in Eq. 46 involving xc cancel. In addition, if Eq. 29 is used to construct the charge density, then the last two terms in Eq. 46 also cancel, leaving:

#### Eq. 49

which is formally the same as the equation derived by other researchers (Andzelm et al. 1989, Versluis and Ziegler 1988). In practice, however, it is necessary to compute both Va/2 and aV/2, because the model density from Eq. 31 is not exactly equal to the numerical charge density computed from Eq. 4.

The time required to evaluate all 3N first derivatives for an N-atom system is typically the same as the time required for 3-4 SCF iterations. If convergence is achieved in, say, 10-12 iterations, then 25-30% additional time is required to evaluate the derivatives. Others have obtained similar results (Versluis and Ziegler 1988).

Because of numerical precision, two potential problems have been observed in evaluating gradients. First, the energy minimum does not correspond exactly to the point with zero derivative (Versluis and Ziegler 1988). The gradients are typically about 10-4 at the energy minimum. A second important point is that the sum of the gradients is not always zero, as it must be for translational invariance. The sum can be as high as 10-3 if the calculation is very poor. Increasing the quality of the integration mesh and the number of multipolar functions in the model density can reduce this to about 3.0 X 10-5. This magnitude of error seems to be permissible for geometry optimizations: the error introduced in the geometry is typically on the order of 0.001 Å. Only for very flat potential energy surfaces should this be a problem.

Currently, the geometry is optimized using both Cartesian and internal coordinates (see Geometry optimization -- The OPTIMIZE suite of algorithms). When the geometry is optimized under conditions of symmetry, only forces that maintain molecular symmetry are evaluated, resulting in considerable time savings. Even in the absence of symmetry, certain forces can be omitted from the calculation, resulting in faster calculations. For example, for a substrate adsorbed on a metal cluster, it is possible to compute gradients for the substrate only and perform no optimization on the metal.

### Electron gas model

The equations actually used in DMol

The parameterized equations for the electron gas exchange-correlation energy as used in DMol are presented explicitly in this section. An excellent review of the electron gas model can be found in Appendix E of Parr and Yang (1989). As mentioned earlier, the implementation in DMol is based on the work of von Barth and Hedin (1972).

#### Eq. 50

The exchange energy

Then the exchange energy is given by:

where:

#### Eq. 52

Spin-restricted computations

The correlation contribution depends upon whether this is a spin-restricted or unrestricted computation. For a restricted computation, the correlation contribution is given by:

where:

co = 0.0225,

r0 = 21,

and:

#### Eq. 54

Spin-unrestricted computations

For a spin-unrestricted computation, the expression is:

where:

and:

## Point-group symmetry

DMol supports most of the chemically important symmetry point groups. If, however, a system is selected with a point group that is not supported, DMol automatically switches to the highest-order subgroup. The Cn, Cnh, and S2n groups are not supported yet, because the complex character tables are not available in DMol. The rotational groups Cv and Dh are also not supported yet, but such systems can be computed in C6v and D6h, respectively, without major loss in efficiency. The supported groups are:

```	C1	Cs	Ci
C2	C2v     C2h      D2      D2h     D2d
C3v		D3	D3h	D3d
C4v		D4	D4h	D4d
C5v		D5	D5h	D6d
C6v		D6	D6h	D5d
Td	O	Oh	I	Ih
```

## Geometry optimization -- The OPTIMIZE suite of algorithms

### Introduction

DMol includes a powerful suite of algorithms for geometry optimization, referred to collectively as OPTIMIZE. Within the Insight environment, you can control OPTIMIZE via parameters in the Optimize/Opt_Parameters command.

After this introductory section, the theory behind the principal algorithms in OPTIMIZE is presented, starting under Theory and implementation. (This section can be skipped by those whose sole interest is how to use the program.) Methodology--Insight, includes a discussion of input parameters, options, and how to use OPTIMIZE within the Insight interface. Lessons demonstrating how to use the Optimize/Opt_Parameters command are described in Tutorial--The Insight Environment.

OPTIMIZE is a general geometry-optimization package for locating both minima and transition states on a potential energy surface. It can optimize in Cartesian coordinates or in a nonredundant set of internal coordinates that are generated automatically from input Cartesian coordinates. It also handles fixed constraints on distances, bond angles, and dihedral angles in Cartesian or (where appropriate) internal coordinates.

The process is iterative, with repeated calculations of energies and gradients and calculations or estimations of Hessians in every optimization cycle until convergence is attained (see scheme in Figure 1). The whole art of geometry optimization lies in calculating the step h so as to converge in as few such cycles as possible.

#### Figure 1     . Schematic of geometry optimization by OPTIMIZE set of algorithms

OPTIMIZE is designed to operate with minimal user input. All you need to provide is the initial geometry in Cartesian coordinates (obtained from the Insight or Discover programs or from an appropriate database), the type of stationary point sought (minimum or transition state), and details of any imposed constraints. All decisions as to the optimization strategy--how to handle the constraints, whether to use internal coordinates, which optimization algorithm to use--are made by OPTIMIZE.

The heart of the program (for both minimization and transition-state searches) is the EF (eigenvector-following) algorithm (Baker 1986). The Hessian mode-following option incorporated into this algorithm is capable of locating transition states by walking uphill from the associated minima. By following the lowest Hessian mode, the EF algorithm can locate transition states, starting from any reasonable input geometry and Hessian.

An additional option available for minimization is GDIIS, which is based on the well known DIIS technique for accelerating SCF convergence (Pulay 1982).

The strategy adopted for constrained optimization depends on the starting geometry and the nature of the constraints. Constraints can be handled easily in internal coordinates, provided that (1) the constrained parameter (distance, angle, or dihedral) is a natural part of the coordinate set and (2) the constraint is rigorously satisfied in the starting structure. If both (1) and (2) hold for all desired constraints, then OPTIMIZE carries out the optimization in internal coordinates. Otherwise, the constrained optimization is performed in Cartesian coordinates.

Traditional wisdom has it that optimization in Cartesian coordinates is inefficient relative to internal coordinates; however, recent work (Baker and Hehre 1991) has clearly demonstrated that, if a reasonable estimate of the Hessian matrix is available (e.g., from a molecular mechanics forcefield) at the starting geometry, then optimization directly in Cartesian coordinates is as efficient as an internal coordinate optimization. In particular, constrained optimization can be handled in Cartesian coordinates as efficiently as with a Z-matrix, with the additional advantages that any distance, angle, or dihedral constraint between any atoms in the molecule can be dealt with (i.e., there is no formal connectivity requirement), and the desired constraint does not have to be satisfied in the starting structure.

OPTIMIZE incorporates a very accurate and efficient Lagrange multiplier algorithm for handling constraints in Cartesian coordinates, with a more robust (but less efficient) penalty function algorithm as a backup. Both algorithms are suitably modified versions of the basic EF algorithm (Baker 1992). The Lagrange multiplier code can locate constrained transition states, as well as minima.

### Theory and implementation

#### The EF algorithm and mode following

Shifting the Newton- Raphson step to favor optimization along an eigenmode

Mode following is a powerful technique for geometry optimization. It involves modifying the standard Newton-Raphson step:

#### Eq. 58

by introducing a shift parameter so that (Cerjan and Miller 1981):

#### Eq. 59

In terms of a diagonal Hessian representation, this can be written:

#### Eq. 60

where the ui and bi are the eigenvectors and eigenvalues of the Hessian matrix H, and = uti g is the component of g along the local eigenmode ui. Scaling the Newton-Raphson step in this way has the effect of directing the step to lie primarily (but not exclusively) along one of the local eigenmodes, depending on the value chosen for .

Various recipes for choosing a suitable shift parameter exist: the EF algorithm utilizes a rational function approximation to the energy, yielding an eigenvalue equation of the form (Banerjee et al. 1985):

#### Eq. 61

from which a suitable can be obtained. This RFO matrix equation has the following important properties:

EF enables both minimization and transition-state optimization

Property 3--the separability of the positive and negative Hessian eigenvalues--allows two shift parameters p and n to be used, one for modes along which the energy is to be maximized and the other for which it is minimized. Specifically, for a transition state (a saddlepoint of order 1) in terms of the Hessian eigenmodes, we have the two matrix equations:

#### Eq. 63

where it is assumed that maximization is along the lowest Hessian mode bi. Note that p is the highest eigenvalue of Eq. 62--it is always positive and approaches zero at convergence--while n is the lowest eigenvalue of Eq. 63--it is always negative and again approaches zero at convergence.

#### Eq. 64

This step may be further scaled down if it is considered too long. For minimization, only one shift parameter n would be used and this would act on all modes. It is often possible to locate different transition states from the same starting structure by maximizing along a mode other than the lowest (hence "mode following").

#### Constrained optimization

Lagrange multipliers as constraints

The essential problem in constrained optimization is to minimize a function of, say, n variables F(x) subject to a series of m constraints of the form Ci (x) = 0, (i = 1 ... m). This can be handled by introducing the Lagrangian function (Fletcher 1981):

#### Eq. 65

which replaces the function F(x) in the unconstrained case. Here, the i are the so-called Lagrange multipliers, one for each constraint Ci(x). Taking the derivative of Eq. 65 with respect to x and gives:

and:

#### Eq. 67

At a stationary point of the Lagrangian function, we have L = 0, that is, all L/xj = 0 and all L/i = 0. This latter condition means that all Ci (x) = 0, and so all constraints are satisfied. Hence, finding a set of values (x, ) for which L = 0 gives a possible solution to the constrained optimization problem in precisely the same thing as finding an x for which g = F = 0 gives a solution to the corresponding unconstrained problem.

We can implement mode following in constrained optimization by simply adopting Eq. 61, but with H replaced by 2L and g replaced by L. However, it is important to realize that each constraint introduces an additional mode to the Lagrangian Hessian (2L), which has negative curvature (a negative Hessian eigenvalue). Thus, when considering minimization with m constraints, you should look for a stationary point of the Lagrangian function whose Hessian has m negative eigenvalues, that is, for a saddle point of order m.

Insofar as mode following is concerned then, assuming a diagonal Lagrangian Hessian representation, Eqs. 62 and 63 for an unconstrained system should be replaced by the following for a constrained system:

#### Eq. 69

where now the bi are the eigenvalues of 2L with corresponding eigenvectors ui and = uti L. Constrained transition-state searches can be carried out by selecting one extra mode to be maximized in addition to the m constraint modes, that is, by searching for a saddlepoint of the Lagrangian function of order m + 1.

#### GDIIS

In the GDIIS method, geometries (xi) generated in previous optimization steps are linearly combined to find the "best" geometry on the current cycle (Császár and Pulay 1984):

#### Eq. 70

Finding appropriate coefficients for use in GDIIS method

The problem here, of course, is to find appropriate values for the coefficients Ci.

and:

#### Eq. 72

are satisfied, then the relation:

#### Eq. 73

also holds.

The true error vectors ei are, of course, unknown. However, they can be approximated by:

#### Eq. 74

where gi is the gradient vector corresponding to the geometry xi. Minimization of the norm of the residuum vector (Eq. 71), together with the constraint equation (Eq. 72), leads to a system of m + 1 linear equations:

#### Eq. 75

where Bij = ei ej is the scalar product of the error vectors ei and ej, and is a Lagrange multiplier.

The coefficients Ci determined from Eq. 75 are used to calculate an intermediate interpolated geometry:

#### Eq. 76

and its corresponding interpolated gradient:

#### Eq. 77

Relaxing the intermediate geometry

A new, independent geometry is generated by relaxing the interpolated geometry according to:

#### Eq. 78

Modifications of the original GDIIS algorithm

In the original GDIIS algorithm, the Hessian matrix is static, that is, the original starting Hessian remains unchanged during the entire optimization. However, updating the Hessian at each cycle generally results in more rapid convergence, and this is the default in OPTIMIZE. Other modifications to the original method include limiting the number of previous geometries used in Eq. 70 by neglecting earlier geometries and eliminating any geometries more than a certain distance (default = 0.3 au) from the current geometry.

## COSMO--solvation effects

The DMol module in release 4.0.0 of the Insight program includes some COSMO controls, which allow for the treatment of solvation effects.

#### Eq. 79

For further details of the COSMO theory see Klamt and Schüürmann (1993). COSMO provides the electrostatic contribution to the free energy of solvation. In addition, there are nonelectrostatic contributions to the total free energy of solvation, which describe the dispersion interactions and cavity formation effects.

### DMol/COSMO

The COSMO electrostatic energy is analogous in form to the DMol electrostatic energy, with the Coulombic operator replaced by the dielectric operator D = BA-1 B. From Eq. 40 we can write:

#### Eq. 80

where represents the auxiliary density, which is introduced to solve the Poisson equation for the electrostatic potential of the solute. This total energy is minimized, resulting in the Kohn-Sham equations for the molecular orbitals. The Kohn-Sham Hamiltonian now includes an electrostatic COSMO potential:

#### Eq. 81

This potential is present in every SCF cycle. This direct incorporation of the solvent effects within the SCF procedure is a major computational advantage of the COSMO scheme. Since the DMol/COSMO orbitals are obtained using the variational scheme, we can derive accurate analytic gradients with respect to the coordinates of the solute atoms. The complete theory is presented in Klamt and Schüürmann (1993) and Andzelm et al. (1995). The gradients include the forces between the solute charges Q and the screening charges q.

To summarize, a DMol/COSMO calculation begins with the construction of the cavity surface. The screening charges are then evaluated using Eq. 79 and the initial solute charges Q. This allows for calculation of the electrostatic COSMO potential. The process is repeated until DMol SCF convergence is achieved. The final total energy includes the DMol/COSMO electrostatic energy (Eq. 80). If geometry optimization is requested, DMol/COSMO gradients are evaluated, and the new geometry is calculated. The next optimization cycle begins with reconstruction of the cavity surface, and the process continues until the DMol optimization convergence criteria are met.

### Determination of the cavity surface (or solvent-accessible surface)

The surface is obtained as a superimposition of spheres centered at the atoms, discarding all parts lying on the interior part of the surface (Klamt and Schüürmann 1993). The spheres are represented by a discrete set of points, the so-called basic points. Eliminating the parts of the spheres that lie within the interior part of the molecule thus amounts to eliminating the basic grid points that lie in the interior of the molecule.

### Determination of nonelectrostatic contributions to the free energy of solvation

The free energy of solvation is calculated as:

#### Eq. 82

where Eo is the total DMol energy of the molecule in vacuo, E is the total DMol/COSMO energy of the molecule in solvent, and Gnonelectrostatic is the nonelectrostatic contribution due to dispersion and cavity formation effects. The nonelectrostatic contributions to the free energy of solvation are estimated from a linear interpolation of the free energies of hydration for linear-chain alkanes as a function of surface area:

#### Eq. 83    Gnonelectrostatic = A_Constant + B_Constant * surface_area

In the present implementation, the VWN potential, DNP basis set, and fine integration grid of DMol were used to calculate energies of methane and octane (C2h). The experimental values of the free energy of solvation are 1.9 and 2.9 kcal mol-1 for methane and octane, respectively (Ben-Naim and Marcus 1984). The default (see next section) set of COSMO parameters was chosen. The calculated surface areas were 38.4 and 104.1 Å2 for methane and octane, respectively. The following values of A_Constant and B_Constant were used:

The best A_Constant and B_Constant values to use depend on the choice of DMol and COSMO input parameters. Our tests (Andzelm et al. 1995) indicate that selection of other parameters, such as the nonlocal BP potential and the van der Waals radii of atoms, can influence Gnonelectrostatic by as much as 0.5 kcal mol-1.

## Electric field gradients

The electric field gradient at a nucleus provides information about the electronic environment of that atom. The property can be probed experimentally by measuring NMR line widths for molecular species that contain NMR-active isotopes with a nuclear quadrupole moment where the quadrupolar relaxation mechanism determines the observed line width. The fact that electric field gradients are a function of the chemical environment is also exploited in NQR (zero-field NMR, or commonly, nuclear quadrupole resonance) scanning to detect explosives (C&E News, 1996).

The electric field gradient is defined as the gradient of the electric field and is therefore the second derivative (or Hessian) of the electrostatic potential at a given position:

#### Eq. 84

where R0 = position at which the derivative is evaluated, Q = electric field gradient, V = electrostatic potential, R = general position vector, and a, b = x, y, or z.

#### Eq. 85

The trace of Q defined as above yields (Poisson equation):

#### Eq. 86

with = the charge density.

#### Eq. 87

In DMol, the tensor Q is computed through expanding the electrostatic potential into a Taylor series with respect to each of the nuclei (if symmetry is imposed, only the symmetry-unique nuclei are visited). The matrix Q then appears as the coefficient of the second-order term in the Taylor expansion:

#### Eq. 88

Uses

The energy of an electrostatic quadrupole moment embedded in an electrostatic potential is determined by the electric field gradient at the position of the quadrupole. This interaction contributes to the NMR line width (as further outlined below).

The electric field gradient at a quadrupolar nucleus (spin I > 1/2) within a molecule determines the NMR longitudinal relaxation rate 1/T1 or the line width W1/2 (this assumes that the quadrupolar relaxation mechanism dominates):

#### Note

In Eq. 89, I = nuclear spin quantum number (e.g., I = 1 for 14N, I = 5/2 for 17O, I = 3/2 for 33S, I = 7/2 for 45Sc), = nuclear quadrupolar coupling constant and is = 2 Q qzz/h, where qzz = largest principal component of the EFG tensor q, (asymmetry parameter) = |qxx - qyy|/qzz, and the x, y, z axes are chosen such that || < 1. (Since enters into the relaxation time as 2, the sign of is not important and may be dropped (Bagno & Scorrano 1996).)

#### Eq. 90

where M = molecular weight, = solute density, and Na = Avogadro's number.

The EFG in atomic units can be converted to SI units with the conversion factor:

14N Q = 2.01 fm2
17
O Q = -2.558 fm2
33S Q = -6.78 fm2
45Sc Q = -22 fm2 (Q < 0 : oblate, Q > 0 : prolate)

Quadrupole interaction energy tensor:

## Fitting atomic point charges to the electrostatic potential (ESP)

Atomic multipole properties are often used to obtain the electrostatic parameters of classical forcefields. One of the most common approaches is to determine the atomic multipole properties by fitting them so as to reproduce the molecular electrostatic potential (ESP) (Singh & Kollman 1984). Numerous applications of ESP-fitted charges in simulations of biochemical systems prove the usefulness of this technique (Bakalarski et al. 1996, Bayly et al. 1993, Merz 1992, Grochowski & Lesyng in press). The ESP-derived charges can reproduce the intermolecular interaction properties of molecules well with a simple two-body additive potential.

#### Eq. 91

Practical implementation of the fitting technique in DMol requires first, calculation of the ESP on a three-dimensional rectangular grid that covers the molecule. The extension of the grid ri is determined by the atomic radii as will be explained in detail below.

#### Eq. 92

where wi are the point weights and , are the atom-centered charges and dipoles, respectively.

The grid points in Eq. 92 are selected based on the following criteria:

#### Eq. 93

where and are the internal and external radii of the atomic shells and depend on the atom type.

To make the results less sensitive to the selection of the grid, the concept of a layer border was introduced. The weights in Eq. 92 change smoothly across the border layer, as evident from the following formula:

#### Eq. 94

where and are the partial weights calculated with respect to all ESP centers in the system:

#### Eq. 96

where R is the "diffusion" width of the layer border.

## Optical absorption spectra

The measurement of optical absorption spectra yields information about the difference between energy levels and about the intensity of the transition between energy levels, which is related to the strength of the coupling with the external electric field inducing the transition.

If the initial state of lower energy is denoted as with energy Ei, and the final, excited state of higher energy is denoted as with energy Ef, an absorption band should be observed at an energy:

DMol calculates both the energy differences Ef - Ei and transition dipole moments IE-FieldF based upon the following assumptions:

The excitation of a single electron does not alter the spin of the electron, because the resulting matrix element IE-FieldF would otherwise vanish. Such a transition is called spin forbidden.

The algorithm for the computation of optical absorption spectra within DMol is:

```Loop over all occupied orbitals |i>
Loop over all unoccupied orbitals |a>
Compute Energy difference dE = Ea - Ei
Compute Matrix elements <i|x,y,z|a>
End Loop
End Loop```
where the occupied Kohn-Sham orbitals are i with orbital energies Ei, and the unoccupied Kohn-Sham orbitals are a with orbital energies Ea.

The computation of transition moments can be made more efficient if the molecular system under consideration possesses symmetry. If it does, selection rules can be applied to avoid computation of matrix elements that are known to be zero on the basis of symmetry arguments alone. The transitions that are not allowed due to symmetry are called symmetry-forbidden. It is possible, of course, that a matrix element is symmetry-allowed but turns out to be zero nevertheless. This happens, for example, when the calculation utilizes only a subgroup of the actual point group of the molecular system.

## Density of states analysis

Once a converged electron density has been calculated, there are several ways to analyze the results, in particular the wavefunction. A convenient way of displaying the molecular orbital spectrum is by constructing and plotting the density of states. For molecular systems, this is commonly done by graphing the molecular orbitals as a function of the MO eigenvalues. The degeneracy of the orbitals is then indicated by the height of the functions. The expression used is a simple sum of delta functions:

#### Eq. 97

For analysis of metallic systems, some sort of artificial broadening is usually applied to the DOS plot. This results in a better match with experimental data obtained from methods like UPS and XPS. Two common ways of doing this are Gaussian and Lorentzian broadening. For Gaussian broadening every energy level is convoluted with a function like:

#### Eq. 98

For Lorentzian broadening we use the function:

#### Eq. 99

The sigma parameter indicates the width of the peaks. For sigma values approaching zero, we obtain the delta representation from Eq. 97. The value of sigma is typically varied so as to best represent the available experimental data.

#### Eq. 100

This does give a reasonable indication of the contribution of that AO to the MO, but a major disadvantage is that the values are not normalized. Adding up the partial DOS for all orbitals in the system does not add up to the total number of electrons of the system, because of the non-orthogonality of the basis functions on different atoms.

#### Eq. 101

where Dij is the AO density matrix. This definition allows the Mulliken density of states to become negative. The Loewdin analysis does not have this drawback, since, because of its definition, it always gives positive values:

#### Eq. 102

where Sif = i | i and S1/2 is the square root of the overlap matrix S. However, we should stress that these different population analysis methods are all somewhat arbitrary, due to the partitioning of S, and care should be taken in attributing too much physical significance to them.

## Mulliken and Mayer bond orders

The concept of bond order and valence indices is well established in chemistry. It allows for interpretation and deeper understanding of the results of DMol calculations using ideas familiar to chemists.

#### Eq. 103

and matrix Pµ and a set of atomic orbitals completely specify the charge density (Eq. 42).

#### Eq. 104

Summing (PS)µ contributions over all µ A, B , where A,B are centers, we can obtain PAB, which can be interpreted as the number of electrons associated with the bond A-B. This is the so-called Mulliken population analysis (Mulliken, 1955). The net charge associated with the atom is then given by:

#### Eq. 105

where ZA is the charge on the atomic nucleus A.

#### Eq. 106

where P, P are the density matrices for spin and .

#### Eq. 108

where PS = P - P.

## Statistical thermodynamics

### Theory

Knowledge of the energy levels of a model allows the computation of its thermodynamic properties using statistical thermodynamics (Herzberg 1945, McQuarrie 1976). The central quantity for statistical thermodynamic computations is the partition function Q of the system:

#### Eq. 109

where the sum runs over all states n of total energy En and degeneracy gn, k is the Boltzmann constant (k = R/Na), and T is the absolute temperature.

This approximation allows one to derive the macroscopic, thermodynamic properties of X from knowledge of the states of isolated molecules, which can be derived from quantum mechanical calculations, if accurate molecular mechanics forcefields are not available.

Dropping this approximation necessitates performing calculations on molecular assemblies in particular ensembles with the help of molecular dynamics or Monte Carlo simulations. These approaches typically require molecular mechanical forcefields, since quantum mechanical molecular dynamics and Monte Carlo simulations are still too demanding to be routine.

Electronic contribution

The separation of electronic contributions is justified where the Born-Oppenheimer approximation is justified for calculating potential energy surfaces.

#### Eq. 110

where Qtrans = translational contribution, Qrot = rotational contribution, and Qvib = vibrational contribution.

Quantum effects are important only for very light molecules at low temperatures (e.g., H2; ortho- and para-hydrogen). In these cases, the coupling between nuclear spin states and rotational states has to be considered.

With all these approximations, we arrive at the following final expressions for Qtrans, Qrot, and Qvib, taking the special cases of atoms and linear molecules into account.

The translational partition function is:

#### Eq. 111

where V = volume (for an ideal gas, related to pressure according to V = RT/p), m = molecular mass, k = Boltzmann factor, T = absolute temperature, and h = Planck's constant.

The vibrational partition function is, within the harmonic oscillator approximation:

#### Eq. 112

which is a product over all vibrations i with degeneracy di (di = 1 for non-degenerate vibrations) and frequency i (in wavenumber units); c is the velocity of light, which couples wavelength and frequency as c = .

#### Eq. 113

where k is the Boltzmann constant.

The computation of the rotational partition function under the assumption of a rigid nuclear framework (rigid rotor, ignoring couplings between rotational and vibrational motion) requires determination of the moments of inertia. The 3 X 3 moments of inertia tensor M is defined as:

#### Eq. 114

where mk = mass of atom k; rk = position vector of atom k relative to the center of mass; and a,b = Cartesian components x,y,z.

#### Eq. 115

According to equalities among the three moments of inertia IA, IB, and IC with respect to the principal axes, the molecule is classified as either:

Linear molecules are symmetric tops according to this classification, yet they have only two rotational degrees of freedom. Rotation with respect to the molecular axis does not affect the nuclear positions. For linear molecules, the partition function is:

#### Eq. 116

where B = h/(8 2 c IB), taking IA = 0, IB = IC.

#### Eq. 117

where A,B,C = h/(8 2 c I(A,B,C)).

#### Eq. 118

Approximate accounting for nuclear spin: the symmetry number

The coupling between nuclear spin and molecular rotation follows the same principle as the coupling between the electron spin and the electronic wavefunction. It is a quantum effect, but can be explained by a classical argument, resulting in the assignment of the so-called symmetry number.

#### Eq. 119

where SYM, the additional factor, is the symmetry number.

Given that electronically excited states are neglected, the ground state does contribute to the partition function if it is degenerate and if a special convention for the zero of the total energy is adopted.

### Thermodynamic quantities

Once the partition function Q is obtained, thermodynamic quantities can be derived. We have to emphasize again that these apply only to molecular gases (where intermolecular interactions are negligible) at temperatures that are neither too low (so that quantum effects are negligible) nor too high (so that the rigid rotor/harmonic oscillator approximation is applicable).

The thermodynamic relations that relate the partition function Q to these quantities are as shown in McQuarrie (1976) for linear polyatomic molecules and nonlinear polyatomic molecules.

The Helmholtz energy A and the heat capacity at constant volume Cv are given through the relationship:

Together with the equation of state, pV = NkT, this yields:

For linear polyatomic molecules

#### Eq. 124

Nonlinear polyatomic molecules

#### Eq. 129

where e = Euler constant (ln e = 1), N = Avogadro's number, k = Boltzmann constant, and = characteristic rotational temperature (P = A,B,C).

### Implementation

Thermodynamic quantities are calculated either as part of the quantum background job without additional user input (Turbomole, DMol) or as part of the analysis capabilities of the Insight Turbomole or DMol module, using the Analyze/Stat_Mech command.

Last updated September 25, 1997 at 03:12PM PDT.
Copyright © 1997, Molecular Simulations, Inc. All rights reserved.