{'fraction_non_alphanumeric': 0.04947033527951614, 'fraction_numerical': 0.012690976470999546, 'mean_word_length': 4.517746913580247, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonvanishing quadrature derivatives in the analytical gradients of density functional energies in crystals and helices

It is shown that the quadrature derivatives in some analytical gradients of energies evaluated with a multi-centre radial-angular grid do not vanish even in the limit of an infinitely dense grid, causing severe errors when neglected. The gradients in question are those with respect to a lattice constant of a crystal or to the helical angle of a chain with screw axis symmetry. This is in contrast with the quadrature derivatives in atomic gradients, which can be made arbitrarily small by grid extension. The disparate behaviour is traced to whether the grid points depend on the coordinate with respect to which the derivative of energy is taken. Whereas the nonvanishing quadrature derivative in the lattice-constant gradient is identified as the surface integral arising from an expanding integration domain, the analytical origin of the nonvanishing quadrature derivative in the helical-angle gradient remains unknown.


I. INTRODUCTION
The analytical energy gradient method of density functional theory (DFT) [1,2] is one of the most heavily used computational machineries in chemistry. It serves as essential infrastructure of chemical research, upon which its accuracy, efficiency, and stability have great impact.
The method is based on Pulay's formalism [3] of the analytical energy gradients for Hartree-Fock (HF) theory. In DFT, an additional term emerges from the exchange-correlation (XC) energy, which is often evaluated by quadrature using Becke's interlocking multi-atom-centre grid [4] (see also Satoko [5]). Its formulation was first reported by Delley [6], who not only established a method for computing these XC energy gradients but also pointed out the existence of a correction due to the variations in grid coordinates and weights. Delley and coworkers [6,7] deemed this correction negligible in atomic gradients. Johnson, Gill, and Pople [8] reported an implementation fully accounting for these quadrature derivatives. Johnson and Frisch [9], Malagoli and Baker [10], and more recently Shang and Yang [11] all concluded that the quadrature derivatives are nonnegligible in vibrational frequencies, although they are generally small and can be made arbitrarily smaller by increasing grid size.
The objective of this article is to point out surprising cases where the quadrature derivatives are large and nonvanishing even in the limit of an infinitely dense grid, causing severe errors when they are neglected. This problem manifests itself in the atomic-basis-set DFT for crystals and helices [12,13].

II. QUADRATURE DERIVATIVES
Let us take an infinite helical polymer (screw axis along x with lattice constant or "rise" a and helical angle or "twist" θ [14]) as our example [12,13]. The atomic positions in the nth helical repeat unit cell (hereafter simply "cell") are rotated anticlockwise by nθ. Therefore, the coordinates of the Ith atom * sohirata@illinois.edu in the nth cell are related to those in the zeroth cell by X I(n) = X I(0) + na, (1) Y I(n) = Y I(0) cos nθ − Z I(0) sin nθ, (2) Z I(n) = Y I(0) sin nθ + Z I(0) cos nθ. ( The formulation for this system encompasses molecules (where nonvanishing quadrature derivatives do not occur), linear (nonhelical) polymers, and two-and three-dimensional crystals. Spin-restricted hybrid DFT formalism [12,13,15,16] for the energy per cell is given in Appendix A. For a more complete formalism of the analytical energy gradient with respect to θ, see Appendix B. Here, we focus our attention on the quadrature derivatives in the XC energy gradients. Atomic units are used throughout this article unless otherwise stated. The XC energy per cell is defined by E XC = a 0 dx dydz f XC (ρ α , ρ β , γ αα , γ αβ , γ ββ ), (4) where ρ's are the electron spin densities and γ's are the spin density gradient invariants [8]. This integral is usually evaluated by quadrature with Becke's interlocking multi-atomcentre grid [4] (see also Satoko [5]) as E XC ≈ I g(I) w I (r g(I) ) f XC (r g(I) ), where I runs over all atoms in the zeroth cell, g(I) over all atom-centre grid points of the Ith atom, r g(I) is the coordinates of the g(I)th grid point, and w I (r g(I) ) is the quadrature weight. It should be stressed that the integration domain is always the zeroth cell (although its boundaries are "fuzzy" [4]) and the grid points and weights are defined only within that cell.
Since the XC energy depends on density ρ, atom-centre basis functions, and atom-centre grid, its analytical gradient with respect to a general coordinate ξ also consists of three terms: The first term with superscript '(b)' involves the direct derivatives of the basis functions (not through variation in grid coordinates). The second term with '(d)' comes from the derivatives of the density matrix elements, which are absorbed in the so-called Pulay force [3]. See Appendix B for the first two terms. The third term with '(q)' is the quadrature derivatives, which consist of the derivatives of grid weights and the derivatives with respect to grid coordinates.

A. In-phase atomic gradients
If ξ is an in-phase atomic coordinate, the sum of the first two terms of Eq. (6) (evaluated by the well-known method [6,8,12]; see Appendix B) is nearly equal to the total XC energy gradient [6] for a typical grid size. Small errors due to the neglect of the quadrature derivatives (the third term) are on the order of 10 −5 a.u. for a grid used in the preceding energy calculation, which are nevertheless noticeable as a small violation of the translational invariance [8]. The quadrature derivatives erase these small errors and restore the translational invariance at near machine precision (see Sec. III).
Taking the x coordinate of atom I in the zeroth cell for ξ, the quadrature derivative is written as ∂w J (r g(J) ) ∂X I(0) f XC (r g(J) ) + g(I) w I (r g(I) ) ∂ f XC (r g(I) ) ∂x .
The first and second terms in the right-hand sides are respectively the derivative of grid weights and the derivative with respect to grid coordinates [6]. They have a large and nearly equal magnitude with opposite signs, cancelling each other nearly exactly. This is why the quadrature derivatives are negligibly small in atomic gradients. The programmable expression of the first term can be found in Ref. [8]. It is however more easily implemented by directly applying the chain rule of differentiation to the code that calculates the quadrature weights. The second term is also easily programmed as a minor modification to the code for x atomic gradients (see Appendix B).

B. Lattice-constant gradient
Next, let us consider the gradient with respect to lattice constant a. Generally, such a gradient should be obtained as a by-product of atomic gradients along x axis by virtue of the following relationship [15]: where n runs over the same range of cell indexes adopted for the lattice sums [17] in the preceding energy calculation, I over all atoms in a cell, and Eq. (1) was used. This gives the correct lattice-constant gradient of the HF energy, as pointed out by Teramae et al. [15]. We found [12] that the above formula does not work for the atomic-basis-set DFT. This is because the integration domain of the XC energy depends on a [see Eq. (4)] and the formula ignores such dependence [12]. Instead of Eq. (6), we had to begin with the following breakdown of the lattice-constant gradient: where the last term with superscript '(s)' stands for a surface integral arising from the explicit dependence of the integration domain on a. Specifically, it is the integral of the XC function over the cell boundary at x = a, which is substantial and should not be neglected. In principle, this can be evaluated by two-dimensional quadrature, but we used the following mathematical trick to convert it into a volume integral and evaluated it with the same three-dimensional quadrature of Becke [12]: where the penultimate equality used Gauss' theorem. In this article, we call this term the surface-integral correction, although it is evaluated as the volume integral. Comparing Eqs. (6) and (12), we can identify the surface integral with the quadrature derivative (albeit not identically but closely; see Sec. III), which is nonvanishing regardless of grid size. The quadrature derivative with respect to a is written as ∂E XC ∂a (q) = I g(I) ∂w I (r g(I) ) ∂a f XC (r g(I) ) + I g(I) w I (r g(I) ) ∂ f XC (r g(I) ) ∂a ∂w I (r g(I) ) ∂a f XC (r g(I) ) + I g(I) w I (r g(I) ) ∂r g(I) ∂a ∂ f XC (r g(I) ) ∂r g(I) = I g(I) ∂w I (r g(I) ) ∂a f XC (r g(I) ). (20) Since no grid points in the zeroth cell move with a (recall E XC is an integral within the volume of the zeroth cell), the second term (the derivatives with respect to grid coordinates) vanishes identically. The first term survives because Becke's fuzzy cell weights for grid points in the zeroth cell are influenced by grid points of nearby atoms, which may move with a. Consequently, near exact cancellation of the two terms does not occur in the quadrature derivatives with respect to lattice constant. The quadrature derivative in the lattice-constant gradient remains large and should not be neglected.

C. Helical-angle gradient
Invoking the same strategy embodied by Eq. (11), we may be led to an identity, This appears to suggest that the helical-angle gradient should be obtainable as a by-product of atomic gradients along y and z axes. This turns out to be false even at the HF level. What is erroneously neglected here [but not in Eq. (11)] is the fact that the basis functions not only shift their centres but also rotate with the helix (and only with these translating-rotating basis functions does the Bloch theorem hold). See Appendix B for the helical-angle gradient formalism that takes into account all θ dependencies including the basis function rotation. ∂w I (r g(I) ) ∂θ f XC (r g(I) ) + I g(I) w I (r g(I) ) ∂ f XC (r g(I) ) ∂θ ∂w I (r g(I) ) ∂θ f XC (r g(I) ) + I g(I) w I (r g(I) ) ∂r g(I) ∂θ ∂ f XC (r g(I) ) ∂r g(I) = I g(I) which also lacks the contribution arising from the derivatives with respect to grid coordinates because the grid is defined only in the zeroth cell and its coordinates do not depend on θ (but its weights can vary with θ). It is therefore nonnegligible no matter how dense the grid is. The nonvanishing quadrature derivative implies that the helical-angle gradient also has an additional term analogous to the surface integral in the lattice-constant gradient. This is not self-evident because the integration domain (i.e., the cell volume) does not expand or contract with θ, but it is only twisted relative to the adjacent domains. We have been unable to identify its mathematical form or analytical origin. Viewed chemically, the variation in θ has the similar effect of lengthening or shortening the chemical bond across the boundary between the zeroth and first cells as does the variation in a. Therefore, if we take as the integration domain a tube enclosing all chemical bonds that is terminated by a plane perpendicular to the bond across the cell boundary, the helical-angle gradient should have a surface integral on the terminal plane, which is analogous to the surface integral in the lattice-constant gradient. Viewed mathematically, however, the concept of chemical bond seems incongruous and the origin of the additional term remains a mystery, although its numerical value can be easily determined as the quadrature derivative. In this sense, taking into account the quadrature derivatives of atom-centre grids seems to be the most general and reliable method of formulating analytical derivatives even when a different numerical integration method is being employed.

III. NUMERICAL ANALYSIS
Some of the B3LYP/6-31G and B3LYP/cc-pVDZ energy gradients [18] of helical polyethylene are given in Tables I and  II, respectively. The structure of the helical repeat unit cell, which is a CH 2 unit, along with the lattice constant or "rise" a and helical angle or "twist" θ, are given in Table III. Note that our B3LYP functional used the VWN5 local correlation functional. The Namur cutoff criterion was adopted with the short-range cutoff of S = 10 and long-range cutoff of L = 30 without a multipole-expansion correction [17]. The number of wave vector (k) sampling points was 100 in the reciprocal unit TABLE I. Energy gradients of helical polyethylene computed by B3LYP/6-31G (using VWN5). Short-range (S ) and long-range (L) lattice sum cutoffs of 10 and 30 were used. The number of k points (K) in the reciprocal unit cell was 100. derivatives" means the analytical gradient calculations without the quadrature derivatives or surface-integral correction. "+ Surface integral" is the "No quad. derivatives" calculation plus the surface-integral correction for the lattice-constant gradient. "Quad. derivatives" stands for the analytical gradient calculations including the quadrature derivatives. "Finite differences" refers to the numerical gradient calculations using a symmetric two-point finite-difference formula with a displacement of ±0.001 bohr or ±0.01 • . b The number of Gauss-Chebyshev radial grid points per atom × the number of Lebedev angular grid points per atom. '8656' refers to a mixed grid of 25 × 302 + 13 × 50 + 12 × 38. c The absolute value of the sum of the three in-phase atomic gradients along x axis, which reports the deviation from the perfect translational invariance. TABLE II. Energy gradients of helical polyethylene computed by B3LYP/cc-pVDZ (using VWN5 and spherical d orbitals). Short-range (S ) and long-range (L) lattice sum cutoffs of 10 and 30 were used. The number of k points (K) in the reciprocal unit cell was 100.  III. The Cartesian coordinates (in bohr) of atoms, the lattice constant or "rise" a (in bohr), and the helical angle or "twist" θ (in degrees) of the CH 2 repeat unit in helical polyethylene.
cell. The integral screening threshold for the two-electronintegral evaluation [19] was 1 × 10 −11 a.u. All calculations were performed with the polymer programme [20].
The grid is defined only in the zeroth cell and consists of interlocking atom-centre grids [4], each of which is a product of the Gauss-Chebyshev radial grid [4] and Lebedev angular grid [21][22][23][24][25][26]. Their orientations are fixed in some user-defined coordinates (cf. Eq. (5) of Ref. [8]) rather than in the molecular axes (cf. Eq. (7) of Ref. [9]), thus making the energy and gradients vary slightly with the user's arbitrary choice of the coordinates since the Lebedev grids are not rotationally invariant. This orientation-dependence is minimized by inclusion of the quadrature derivatives [9] and eventually eliminated in the limit of an infinitely dense grid. Conclusions of this article are unchanged by the lack of rotational invariance of the grids.
In-phase atomic gradients change noticeably upon inclusion of the quadrature derivatives for the small (25 × 86) grid consisting of 25 radial points × 86 angular points per atom. With increasing grid size, the differences between gradients with and without the quadrature derivatives taper below 2 × 10 −7 a.u. Concomitantly, the deviations from the translational invariance (the absolute value of the sum of the x atomic gradients) without the quadrature derivatives are on the order of 10 −4 a.u. or more for the 25 × 86 grid, but they decrease to 10 −8 a.u. or less for the 300 × 974 grid. With the quadrature derivatives, the deviations become on the order of the integral screening threshold regardless of grid size.
Curiously, the atomic gradients without the quadrature derivatives seem consistently closer to the values in the limit of an infinitely dense grid. However, the atomic gradients with the quadrature derivatives are within 1 ×10 −7 a.u. of the finitedifference benchmarks regardless of grid size, and are therefore more precise. With a grid that is sufficiently large for energy calculations (the 8565 or 100 × 302 grid), the differences are negligible for atomic gradients, if not for frequency calculation [9][10][11].
The lattice-constant gradient without the quadrature derivatives or surface-integral correction differs from the correct value by an order of magnitude and has the wrong sign even for the gigantic 300 × 974 grid. The addition of the surface integral brings it within 2 × 10 −4 a.u. of the correct value for the 25 × 86 grid and within less than 10 −7 a.u. for the 300 × 974 grid. On the other hand, the deviations of the lattice-constant gradient with the quadrature derivatives from the finite-difference benchmarks are on the order of 10 −7 a.u. for the 8656 or 25 × 86 grid (these slight deviations are due to the finite-difference approximation). The more precise agreement between the latter two methods indicates that the surface integral method is slightly approximate and the method including the quadrature derivatives should always be trusted, although the error in the surface integral is negligible. Not including either the surface-integral correction or quadrature derivatives will lead to unacceptably large errors in the latticeconstant gradient, rendering the method useless.
The helical-angle gradient without the quadrature derivatives is wrong by the same order of magnitude as the gradient itself regardless of grid size for this particular geometry. The value with the quadrature derivatives agree with the finitedifference benchmark within a few 10 −7 a.u. for the 8656 or 25 × 86 grid, and these small errors are ascribed to the finitedifference approximation. It is therefore crucial to include the quadrature derivatives in this gradient also. In atomic units, energy per helical repeat unit cell of an infinitely extended helix is written as which encompasses the spin-restricted HF theory (c = 1 and E XC = 0) and a pure local or gradient-corrected DFT (c = 0 and E XC 0) as well as a hybrid DFT such as B3LYP [18] (c 0 and E XC 0). Generalizing it to the spin-unrestricted case should be straightforward.
The first term is the nuclear-repulsion energy, which is written as where I and J run over all nuclei with nuclear charges of Z I and Z J , and n runs over cells (−L ≤ n ≤ L) with L standing for the long-range lattice sum cutoff [17]. Summands with zero denominator are excluded from the summation. R I(n) is the position of the Ith nucleus in the nth cell, whose elements are defined by Eqs. (1)-(3). The nth cell is twisted anticlockwise by nθ relative to the zeroth cell. The second term is the kinetic energy per cell expressed as whereP µ(0)ν(n) is the density matrix element (defined at the end of this Appendix) for the µth contracted Gaussian-type orbital (cgs) in the zeroth cell and the νth cgs in the nth cell that translates and rotates with the helix (defined in the next paragraph). S is the short-range lattice sum cutoff [17].T µ(0)ν(n) is the corresponding kinetic-energy integral, defined bỹ whereχ ν(n) (r) is the translating-rotating νth cgs in the nth cell. Generally, tilde indicates screw-axis-symmetry-adapted quantities throughout this article. Specifically, a screw-axis-symmetry-adapted cgs is given as a linear combination of the also screw-axis-symmetry-adapted primitive Gaussian-type orbital (pgs). It is therefore defined byχ where c κα is the product of contraction coefficient and normalization factor of the αth pgs, ζ α is the exponent of the αth pgs, and R κ(n) is the coordinates of the nucleus in the nth cell on which the κth cgs and all of its constituent pgs are centred. The coordinates (x (n) ,ỹ (n) ,z (n) ) are the nth-cell-fixed Cartesian coordinates that rotate with the helix and are centred on R κ(n) = (X κ(n) , Y κ(n) , Z κ(n) ). They are related to the space-fixed Cartesian coordinates r = (x, y, z) bỹ Only with these screw-axis-symmetry-adapted (i.e., translating-rotating) basis functions does the Bloch theorem [27] hold, guaranteeing the translational invariance of matrices with tilde:T µ(m)ν(n) =T µ(0)ν(n−m) , (A10) P µ(m)ν(n) =P µ(0)ν(n−m) . (A11) However, various molecular integrals are more conveniently evaluated [19] for cgs in the space-fixed Cartesian coordinates without tilde: (A12) The kinetic-energy-integral matrix in these non-symmetryadapted cgs is defined by and computed as using the standard algorithm such as the Obara-Saika method [19]. Generally, quantities without tilde designate those in the basis of the space-fixed Cartesian coordinates in this article. The screw-axis-symmetry-adapted kinetic-energy integrals are obtained from the non-symmetry-adapted ones by the transformation,T where R νν ′ (n) is the rotation matrix for the nth cell. The latter is a unit matrix in which the 3-by-3 block in a normalized p subshell (ordered p x , p y , p z ) is replaced by with c = cos nθ and s = sin nθ, the 6-by-6 block in a normalized Cartesian d subshell (ordered d xx , d yy , d zz , d xy , d yz , d zx ) is substituted by and the 10-by-10 block in a normalized Cartesian f subshell The third term of Eq. (A1) is the nuclear-attraction energy given by where the nuclear-attraction integrals over the translatingrotating basis functions are obtained by the transformation, The nuclear-attraction integrals are computed in the spacefixed Cartesian coordinates as (A20) The subsequent two terms in Eq. (A1) are the Coulomb and exchange energies given by where the so-called J and K integrals in the translatingrotating basis functions are obtained by the same transformation,J The J and K integrals in the space-fixed Cartesian coordinates are, in turn, built from two-electron integrals: with in the Mulliken notation. The density matrix in the space-fixed Cartesian coordinates entering the above equations is obtained by the back-transformation, which is no longer expected to possess the translational invari-ance, i.e., The penultimate term of Eq. (A1), E XC , is the exchangecorrelation (XC) energy evaluated by quadrature [see Eq. (5)]. Integrand f XC has the same screw axis symmetry as the helix and so the XC energy per cell can be evaluated as its integral over the zeroth cell volume. In practice, an infinitely extended f XC is divided into cellwise contributions in a gradually overlapping manner by Becke's fuzzy cell scheme.
The electron density and density gradient invariant are defined only within the zeroth cell (where the space-fixed and translating-rotating Cartesian coordinates coincide) and evaluated explicitly with the translating-rotating basis functions.
where Eqs. (A8) and (A9) were used. The XC integrals, necessary in the self-consistent field pro-cedure (see below), are also evaluated by quadrature.

I g(I)
w i (r g(I) )∇χ * µ(m) (r g(I) ) · ∇ρ(r g(I) ) × u XC (r g(I) )χ ν(m+n) (r g(I) ), where the XC potentials are defined only within the zeroth cell as The derivatives of the translating-rotating basis functions with respect to the space-fixed Cartesian coordinates are evaluated as The derivatives in the right-hand sides are computed on each grid point explicitly. In Eqs. (A30)-(A36), two translating-rotating basis functions are generally centred in different cells, while v XC is infinitely extended and periodic. Therefore, in principle, one needs to integrate the respective functions in all cells, using grid that is translated and rotated along the helix. An equivalent integration is achieved by defining grid only within the zeroth cell and translate and rotate the two basis functions along the helix. We have chosen the latter method and the above equations reflect this choice.
The last term, E MPE , stands for the multipole-expansion correction [17]. In the leading order, it is given by where (µ x , µ y , µ z ) is related to the zeroth-cell dipole moment, and ∞ is mimicked by a large integer (10,000 in our implementation). Note that the charge-quadruple interactions, which also decay as r −3 and were recommended for inclusion by Delhalle et al. [17], are excluded in this study, so that the energy is variational and consistent with the analytical gradients (see Appendix B). The corresponding correction to the Kohn-Sham Hamiltonian (or Fock) matrix elements arẽ Integrals x µ(0)ν(n) , y µ(0)ν(n) , and z µ(0)ν(n) are the moment integrals of the space-fixed Cartesian coordinates in the nonrotating basis functions. For instance, a y moment integral is defined by and evaluated as which, in turn, gives the y "dipole moment" for the zeroth cell as where the density matrix in the space-fixed Cartesian coordinates is given by Eq. (A28). The Kohn-Sham Hamiltonian (or Fock) and overlap matrices in the translating-rotating basis functions are defined bỹ andS with S µ(0)ν(n) = dr χ * µ(0) (r)χ ν(n) (r). (A49) The dynamical Fock and overlap matrices are then constructed asF where k is the wave vector (−π/a ≤ k < π/a). At this point, the dynamical matrices in the spherical d and f translatingrotating basis functions (if applicable) can be obtained by the transformations,F The transformation matrix is a unit matrix in which the 6-by-6 block in a normalized Cartesian d subshell (in the same order as before) is replaced by or the 10-by-10 block in a normalized Cartesian f subshell (in the same order as before) is swapped by The generalized eigenvalue equation of the form, is solved with the Löwdin orthogonalization method [28] (at each of the K evenly spaced wave vector k in the reciprocal unit cell) for the energy bands ǫ m (k) and Bloch orbital coeffi-cientsC νm (k) in the spherical d and f translating-rotating basis functions. The latter are transformed to those in the Cartesian d and f translating-rotating basis functions bỹ The density matrix in the Cartesian d and f translatingrotating basis functions is then given bỹ with f m (k) being the Fermi-Dirac distribution function, where β = (k B T ) −1 and the Fermi energy ǫ F is determined by the electro-neutrality condition: Appendix B: Energy gradients of a helical polymer with respect to the helical angle We shall consider the gradient of the total energy per cell with respect to the helical angle θ as the most complex example: Gradients with respect to the in-phase atomic coordinates X I(0) , Y I(0) , Z I(0) , or to the lattice constant a can be inferred from it. The first term in the right-hand side of the above equation is the Hellmann-Feynman force, the second term has to do with the derivatives of basis functions (but not through grid coordinates), the third term the density matrix derivatives, and the fourth term the quadrature derivatives. The Hellmann-Feynman force consists of two terms: each of which is easily formulated and will not be repeated here. The density matrix derivatives are consolidated into the so-called Pulay force by virtue of the variational nature of the HF or DFT energy. It is given by ∂E ∂θ with the energy-weighted density matrix, The quadrature derivatives exist only in the XC energy gradient and are extensively discussed in the main text.

∂E ∂θ
Therefore, in the remainder of this Appendix, we shall focus on the basis function derivatives: . (B6) The basis function derivative of the kinetic energy with respect to θ is given by In the last expression, the first term accounts for the basis function rotation, while the second for the basis function centre shift. Equation (22) is wrong for it neglects the first term. The derivative of the kinetic-energy integral in the space-fixed Cartesian coordinates with respect to θ is computed as where we used the fact that the basis functions in the spacefixed Cartesian coordinates depend on θ only through their centres and that ∂/∂Y = −∂/∂y, etc. The molecular integral derivatives appearing in the last expression are evaluated by the standard method [19]. The corresponding derivatives of S µ(0)ν(n) andÑ µ(0)ν(n) are obtained in the same way.
The basis function derivative of the Coulomb energy is given by The basis function derivative explicitly neglects the variation of density matrix elements. Hence, what is meant by the basis function derivative of the density matrix in the space-fixed Cartesian coordinates is as follows: which includes only the effects of basis function rotations on the density matrix, but not the change in the density matrix (P) itself. The derivative of a two-electron integral with respect to θ is given by ∂ ∂θ (µ(0)ν(n)|κ(l)λ(l + m)) = (−nY ν(0) sin nθ − nZ ν(0) cos nθ) µ(0) − ∂ν(n) ∂y |κ(l)λ(l + m) The corresponding derivative of the exchange energy can be obtained in the same way.
The basis function derivative of the XC energy reads I g(I) w I (r g(I) ) f XC (r g(I) ) (B17) = nuc.