{'fraction_non_alphanumeric': 0.029702429702429704, 'fraction_numerical': 0.006060606060606061, 'mean_word_length': 4.7084307308711235, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Full-Potential LMTO: Total Energy and Force Calculations

The essential features of a full potential electronic structure method using Linear Muffin-Tin Orbitals (LMTOs) are presented. The electron density and potential in the this method are represented with no inherent geometrical approximation. This method allows the calculation of total energies and forces with arbitrary accuracy while sacrificing much of the efficiency and physical content of approximate methods such as the LMTO-ASA method.


Introduction
This paper describes a particular implementation of a full-potential electronic structure method using Linear Muffin-Tin Orbitals (LMTO's) [1,2] as basis functions. There have been several "FP-LMTO" implementations [3,4,5,6,7]. The one described here has not been published in detail, although calculations performed with this method have been reported for quite some time. [3] There are many aspects to an electronic structure method. This paper is focussed on those aspects which enable a full potential treatment. Relatively small details pertaining to full-potential methods will be discussed while larger details having to do with, for example, relativity will not be.
The emphasis of a variational full-potential method is somewhat different from that of a method such as the LMTO-ASA method. The emphasis of the former is on the completeness of the basis while in the latter it is in the physical content (and interpretability) of the basis. These concepts are, of course, intimately related, but the emphasis is different.
The exposition here is for an infinite system periodic in three dimensions. This method has been implemented for two-dimensional systems, [8] but that will not be discussed here.

Notation
Papers on electronic structure methods unavoidably carry a high overhead in functional symbols and indices. It is simplest to define here, without motivation, the special symbols and functions that will be used in this paper, for future reference. These special functions (although not necessarily the symbols used here) have been used extensively in LMTO documentation and are largely due to Andersen. [1] Spherical harmonics: where Y is a spherical harmonic. [9] Bessel functions: K ℓ (κ, r) ≡ −κ ℓ+1 n ℓ (κr) − ij ℓ (κr) κ 2 < 0 n ℓ (κr) κ 2 > 0 (4) where L denotes ℓm and n ℓ and j ℓ are spherical Neumann and Bessel functions, respectively.
Geometry: For computational purposes, the crystal is divided into nonoverlapping spheres surrounding atomic sites (muffin-tin spheres) where the charge density and potential vary rapidly and the interstitial region between the spheres, where the charge density and potential vary slowly. This is the muffin-tin geometry used as an idealized potential and charge density in early electronic structure methods (KKR and APW). Here, the division is a computational one, and does not restrict the final shape of the charge density or potential. In the muffin-tin spheres, the basis functions, electron density, and potential are expanded in spherical waves; in the interstitial region, the basis functions, electron density, and potential are expanded in Fourier series. There are many relevant considerations in choosing muffin-tin radii. Assuming all expansions are taken to convergence, the density and potential depend on the muffin-tin radii only through the dependence of basis functions on the radii. As discussed below, basis functions have a different functional form inside the muffin-tin spheres, and the choice of muffin-tin radius affects this crossover. Hence, assuming the Hamiltonian is the same inside and outside the spheres (the treatment of relativity may affect this as discussed below), the muffin-tin radii are variational parameters and the optimum choice minimizes the total energy. If the basis is large enough however (suitably complete within and without the spheres), the energy is insensitive to the choice of radii. A reasonable choice results from choosing radii that are within both the minimum in charge density and the maximum in potential along a line between nearest neighbors. Relativistic effects are usually taken into account only in the muffin-tin spheres, in which case the Hamiltonian depends on the radii; hence when relativistic effects are important, the radii are not variational parameters.
In what follows, lattice positions are vectors R = Rn, integer multiples of a basis R. Atomic positions in the unit cell are denoted by τ . A set of atomic positions invariant under the point group of the lattice are said to be of the same symmetry type, t. Similarly, in the reciprocal lattice, vectors are g = Gn for the reciprocal basis G = 2πR −T . Brillouin zone (or reciprocal unit cell) vectors are denoted by k.
Symmetric functions: Within the muffin-tin region, functions invariant are expressed in harmonic series. If f (r) is such a function, at site τ In Equation (8a), D τ is a transformation to a coordinate system local to site τ ; the local coordinates of sites of the same type are related by an element of the crystal point group that takes one site into another. Expressed in this way, the functional form of D ht (Equation (8b)) depends only on symmetry type.
In the interstitial region, symmetric functions are expressed in Fourier series: The sum in Equation (9a) is over symmetry stars S of the reciprocal lattice.

Interstitial
In the interstitial region (symbolically I) between the muffin-tin spheres, bases are Bloch sums of spherical Hankel or Neumann functions: The rotation D τ in (10) takes the argument into a coordinate system local to each site τ . The purpose of this will be made evident later. The function on the right hand side of Equation (10) is sometimes called the envelope function. Notice the parameters, specifying a basis function, inherent in this definition. They are the site τ in the unit cell on which the spherical wave is based, the angular momentum parameters ℓ and m of the spherical wave with respect to its parent cell, and the kinetic energy κ 2 of the basis in the interstitial region. The angular momentum parameters specifying the basis set are chosen to represent the atomic states from which crystal eigenstates are derived. In the LMTO-ASA, it is usual to include ℓ bases one higher than the highest relevant band. In the method described here, this is rarely necessary, possibly because of the multiplicity of bases with the same angular momentum parameters. It is usual to use "multiple κ" basis sets, having all parameters except the tail parameter the same.
There appears to be no simple algorithm for choosing a good set of interstitial kinetic energy parameters. Schemes such as bracketing the relevant energy spectrum have been proposed. [7] The optimum set would minimize the total energy. This can be done but is time consuming even for relatively simple systems. It seems, however that parameter sets obtained in this way for simple systems in representative configurations can give good results when used for related systems over a broad pressure range. Thus good sets are arrived at through some experimentation. The choice can be important as it's possible to pick a set of parameters that will give very bad results, and the parameter set used in any new calculation should be always checked for stability.

Muffin Tins
In the muffin-tin spheres, bases are linear combinations of spherical waves matching continuously and differentiably to the envelope function at the muffin-tin sphere. The envelope function K may be expanded in a series of spherical Bessel functions about any site except it's center. A basis function on a muffin-tin sphere in the unit cell at R = 0 is therefore Full-Potential LMTO where r τ ≡ r − τ and B is equivalent to the KKR structure constant. [10] The unitary transformation applied to B rotates components into site-local coordinates from the left and right.
Equation (11) is compactly expressed by defining a two-component row vector K so that and a two component column vector S so that Then the value of a basis function on a muffin-tin boundary is expressed simply as The radial part a basis function inside a muffin-tin sphere is a linear combination of atomic like functions φ and their energy derivativesφ [1,2] matching continuously and differentiably to the radial function K in Equation (14). Collecting φ andφ in a row vector U(e, r) ≡ φ(e, r),φ(e, r) , a simple case of this matching condition may be expressed as U(e, s)Ω(e, κ) = K(κ, s) and U ′ (e, s)Ω(e, κ) = K ′ (κ, s), where Ω is a matrix of order 2.
The use of these radial functions in the method described here is different than that used by most other methods, however. For the broadest utility, a basis set must be flexible enough to describe energy levels derived from atomic states having different principle quantum numbers but the same angular momentum quantum number. For example, describing the properties of elemental actinides at any pressure requires a basis with both 6p and 7p character. Similarly, an adequate calculation of the structural properties of transition metal oxides requires both semi-core and valence s and p states on the transition metal ions. The description of the evolution of core states from localized to itinerant under pressure also requires multiple principle quantum numbers per ℓ value. It is usual in LMTO-based methods to perform calculations for the eigenstates and eigenvalues of "semi-core" and valence states separately, using a different basis set, with a single set of energy parameters {e ℓ }, for each "energy panel". This approach fails when energy panels overlap, and has the disadvantage that the set of eigenvectors is not an orthogonal set. The problem of "ghost bands" also arises. [2] In the method described here, bases corresponding to multiple principle quantum numbers are contained within a single, fully hybridizing basis set. This is accomplished simply by using functions φ andφ calculated with energies {e nℓ } corresponding to different principal quantum numbers n to describe the radial dependence of a basis in the muffin-tin spheres. The Hamiltonian matrix for an actinide, for example, will have elements ψ 6p H ψ 7p and the overlap matrix elements ψ 6p ψ 7p , We may formally express the radial part of basis i in a muffin-tin sphere by the function f (r) = n a i (nℓ)U(e nℓ , s)Ω(e nℓ , κ i ) but in practice it is sufficient to restrict the coefficients by a i (nℓ) = δ(n, n i ) so that the basis set (although not eigenvectors) will have pure principal quantum number "parentage". This method of expanding the energy range of a basis set has been used (and reported) extensively. Representative calculations in which this method was essential are described in Reference [11].
Thus another parameter specifying a basis function is the set of energy parameters {e tℓ } that will be used to calculate the radial basis functions φ tℓ andφ tℓ used to express the basis function in muffin-tin spheres of each symmetry type. A basis function in a muffin-tin sphere is therefore where e i means "use the energy parameter e nℓ corresponding to the principal quantum number specified for basis i" and The necessary cutoff in angular momentum has now been made explicit. The 2 × 2 matrix Ω matches U to K continuously and differentiably at the muffintin radius. Specifically, Ω is specified by In principle, and as programmed, each (τ ℓκ) basis can use its own unique energy set. It is more usual to use a common energy set for a set of basis states giving rise to bands of similar energy within the scope of a particular calculation. The configuration of the basis shown in Table 1 for example uses a set of energies for "semi-core" 6s and 6p bases, and another set of energies to represent "valence" bases. The calculation of energies in an energy parameter set is discussed below.
A parameter introduced in (15) is the angular momentum cut-off ℓ m . In most cases, a converged total energy is achieved with values ℓ m ∼ 6 − 8. Note that since a basis set generally contains functions based on spherical waves with ℓ ≤ 3, the KKR structure constant in (13) is rectangular. Table 1. Parameters for typical basis set for an elemental actinide: parent angular momentum parameter (ℓ), energy set for radial expansions (e-set), and the index of the kinetic energy in the interstitial region(κ-index). A typical set of κ 2 values, corresponding to the kinetic energy indices, is given at the bottom of the table.

Muffin-Tin Matrix Elements
The potential in a muffin-tin at τ has an expansion in linear combinations of spherical harmonics invariant under that part of the point group leaving τ invariant: The utility of referring bases and potentials in muffin-tin spheres to sitelocal coordinates is apparent in (18a). If the site local coordinates of sites are constructed so that D τ ′ = D τ Q −1 for some Q such that Qτ = τ ′ , then the harmonic functions D ht depend only on the symmetry type, rather than on each site. The normalization for the spherical harmonic in (18a Combining (15) and (18a), the potential matrix is The matrix element of the D ht is a sum over Gaunt coefficients: In electronic structure methods using muffin-tin orbitals, the muffin-tin energy parameters {e ℓ are usually taken from "ℓ-projected average energies". With multiple energy sets, this is a reasonable choice provided that the basis set, which uses separate sets, gives rise to bands well separated in energy. The ℓ-projected charge, integrated over a muffin-tin sphere, is a sum over cross terms between energy sets and must be made diagonal in some approximation for the resulting energyand ℓ-projected energies and charges to be representative.
Another criterion, particularly useful for states using different sets not well separated in energy or for states not having significant occupation is to maximize the completeness of the basis. To accomplish this, the energy parameter for the low energy state e ℓ (1) can be set to a set of projected energy averages, and the energy parameters for the same ℓ in higher energy sets may be chosen so that the radial function has one more node and the same logarithmic derivative at the muffin-tin radius, hence Although this usually generates energy parameters out of the range of occupied states (since the logarithmic derivative of semi-core states is usually large in magnitude and negative), this choice seems to give a total energy close to the minimum with respect to this parameter. This is an example of the difference mentioned in the introduction in emphasis between an accurate "basis-set" method and a method motivated by a physical model. The convergence of the harmonic expansion of the potential in a muffintin sphere (18a) depends, of course, on the basis, atomic constituents, and geometry. Using harmonics through ℓ hmax = 6 is usually sufficient, and it has never been necessary to go beyond ℓ hmax = 8.

Interstitial Matrix Elements
Overlap and Kinetic Energy: The interstitial overlap matrix can be easily obtained from an integral over the interstitial surface (the only non-zero contributions, in a crystal periodic in three dimensions, come from the surfaces of the muffin-tin spheres) and the kinetic energy is proportional to the overlap: where Basis functions on muffin-tin spheres are given in (14), hence In the limit κ 2 j → κ 2 i , the evaluation of (21) requires the derivative with respect to κ 2 of the structure constant.
Potential Matrix Elements: The greatest difference between LMTObased full-potential methods is in the way the matrix elements of the potential are calculated over the interstitial region. The method being described here uses a Fourier representation of basis functions and the interstitial potential to calculate these matrix elements. Other approaches for computing these elements are described in the literature. [4,5] A Fourier transform of the basis functions described in Section 2 would be too poorly convergent for practical use. However, the evaluation of the interstitial potential matrix requires only a correct treatment of basis functions and potential in the interstitial region. This degree of freedom can be used to design "pseudo basis-set", equal to the true basis in the interstitial region although not in the muffin-tin spheres, and have a Fourier transform which converges rapidly enough for practical use. We define this pseudo basis set byψ Since rapid Fourier convergence is the criterion for constructing the pseudobasis, it is useful to consider the Fourier integral of a Bloch function with wave-number k: where V c is the unit cell volume. Equation (24) is obtained by casting ∇ 2 +κ 2 on the plane wave then doing two partial integrations; surface terms vanish due to periodicity. From (24) it is evident that the Fourier integral of a pseudo-basis satisfying the first criterion (equal to the true basis in the interstitial region) may be obtained from integral over muffin-tin spheres. If in addition, the pseudo-basis is different from a Hankel function only in it's parent sphere, the Fourier integral is a finite integral over a single muffin-tin sphere. The problem then is to find a functionψ such that (∇ 2 +κ 2 )ψ has a rapidly convergent Fourier integral, vanishes outside a radius less than or equal to the parent muffin-tin radius for the basis, and has a value and slope equal to K at this radius. A good choice for such a function is obtained by solving for a radius s < s ti , and with with c ℓ chosen to match on to K at s. This is easily done analytically. The resulting Fourier transform is where N = ℓ i +n i +1. The subscript i has been purposely left off N and s (see below). These coefficients converge like 1/g n+4 , provided J N (|k+g|, s) achieves it's large argument behavior, and n can be chosen to optimize convergence. Weinert [12] used an analogous construction as tool to solve Poisson's equation. He proposed a criterion for the convergence of the Fourier serie (26) which amounts to choosing the exponent n in Equation (26) so that |k+g max |s would be greater than the position of the first node of J ℓ+n+1 . We find this criterion to be useful provided anisotropy in reciprocal space is accounted for. This is accomplished by using the minimum reciprocal lattice vector on the surface of maximal reciprocal lattice vectors, rather than simply using g max .
Notice that this criterion is a criterion for N = ℓ+n+1. The basis Fourier components are simplified, and the amount of information stored reduced, by simply using a single argument for all bases; i.e. all bases use the same value of N . It is also possible to use a single radius s, less than or equal to the smallest muffin-tin radius, since the only requirement is on the pseudo bases in the interstitial region. In practice, a few radii are desirable if large and small atoms are present in the same calculation, since small radii give less convergent Fourier coefficients. In any event, no more than a few radii are necessary to handle systems with many atoms. Notice also that local coordinates have been left out of (26). The resulting potential matrix may be easily rotated to local coordinates at the end of the calculation.
As expressed in (26), the Fourier components are products of phases e −i(k+g)·τ , which scale like the number of atoms squared (the size of the reciprocal lattice grid grows linearly with the number of atoms), and a function of lattice vectors and a few parameters, which scales linearly with the number of atoms. The phase factors are simple to calculate by accumulation and need not be stored.
The potential in the interstitial region is similarly obtained from a "pseudopotential"Ṽ that equals the true potential in the interstitial region and has rapidly converging Fourier coefficients: The sum in Equation (27b) is over stars S of the reciprocal lattice. Integrals over the interstitial region are performed by convoluting the potential with an interstitial region step function and integrating over the unit cell: The potential matrix element is calculated by convoluting the convoluted potential with a basis, and performing a direct product between convoluted and unconvoluted bases. If basis functions are calculated n 3 reciprocal lattice vectors, the interstitial potential will be calculated on (2n) 3 vectors. The convolution is exact if it is carried out on a lattice containing (4n) 3 vectors. The size of the set of reciprocal lattice vectors necessary to converge the total energy using this treatment of the interstitial region varies from between ∼ 150 -300 basis plane waves per atom, depending on the smoothness of the potential and the convergence required.
Another way of integrating over the interstitial region, more usual in sitecentered methods, is to integrate Fourier series over the unit cell and subtract the muffin-tin contributions with pseudo-bases and pseudo-potential expressed as an expansion in spherical waves. The convolution has an advantage in acting with a single representation, and, given a finite representation for bases and potential, the convolution may be done exactly.
Empty spheres are never used with this scheme. Bases, and the charge density and potential are calculated as accurately as necessary using the scheme described above and a basis set expanded with tail parameters and energy sets has proven to be flexible enough to accurately describe the contribution of the electronic states in the interstitial region.

Charge Density
When a solution to the wave equation at every physical energy is available, the charge density may be obtained from a set of energy-dependent coefficients.
The spherically symmetric charge density in a muffin-tin sphere, coupled with an ℓ − projected density of states, is an example. In a variational calculation, as is being described here, all that is available is a (variational) solution to the wave equation at a set of discreet energies, and the charge density must be obtained simply from the square of the eigenvectors, or equivalently from expectation values of occupation numbers.
Having calculated a set of eigenvalues and eigenvectors A of the generalized eigenvalue problem, the charge density in the interstitial region is n(r) where N S is the number of vectors in the reciprocal lattice star S. The square of the wave function is obtained by convoluting the Fourier components of ψ with A, Fourier transforming, and taking the modulus.
In the muffin-tin spheres the charge density is The process of calculation is evident in the sequence of equations.

Core States
Core states, even spherically symmetric complete shells, contribute non-muffintin components to the interstitial region and to muffin-tin spheres surrounding other sites. Whether it is essential to include this contribution depends on the size of the contribution, and any sizable contribution implies that there are states being treated as localized which aren't localized within the scope of the calculation. Nevertheless, confining states to the core is often useful, and including the core contribution to the full potential is not difficult. One possibility, the one used in this method, is to fit the part of the core electron density to a linear combination of Hankel functions, and expand this density in the interstitial region as a Fourier serie and in the muffin-tin spheres in a harmonic series, in the same way the basis functions are treated.

Coulomb Potential
The Coulomb potential is obtained by first calculating the Coulomb potential in the interstitial region, then, using the value of the interstitial potential on the muffin-tin sphere, calculating the potential in the spheres by a numerical Coulomb integral of the muffin-tin electron density for each harmonic. The interstitial Coulomb potential is calculated in a way similar to that suggested by Weinert [12]. Express the electron density as whereñ is the squared modulus of the pseudo-eigenvectors, which is equal to the true electron density in the interstitial region. The first term on the righthand side of (30) has, by construction, a convergent Fourier series. The second term is confined to muffin-tin spheres. To calculate the Coulomb potential in the interstitial region, this term may be replaced by any density also confined to the muffin-tin spheres and having the same multipole moments. If a charge density satisfies these requirements and also has a convergent Fourier series, the Coulomb potential in the interstitial region may be easily calculated from the combined Fourier series. Such a charge density can be constructed in a similar way to that detailed for the pseudo-bases. Construct a pseudo chargedensity satisfying This charge density has Fourier components where the multipole moments Q are defined by The Fourier componentsñ (p) (r) converge like 1/g n+2 provided j ℓ+n+1 attains it's asymptotic form. The exponent n is chosen using the same considerations as for the pseudo-basis set.
The Coulomb potential in the interstitial region is then given by From the Coulomb potential in the interstitial region follows the Coulomb potential on the surface of the muffin-tin spheres. The coulomb Potential inside the muffin-tin spheres is is the harmonic component of the potential on a sphere boundary.

Density Gradients
Gradients of the electron density are needed for the evaluation of gradient corrected density functionals. These functionals depend on invariants (with respect to the point group) constructed from density gradients (e.g. |∇n| 2 ). This reduces computation significantly in the muffin-tin spheres, for if f and g are invariant functions (i.e. f (r) = h f h (r)D h (r)), and d = ∇f · ∇g, then d(r) = h d h (r)D h (r) with where the set of parameters I is easily calculable from 3j and 6j coefficients and integrals over the harmonic functions D h , and and similarly for g. Gradients of the interstitial charge density, represented as a Fourier series, are poorly represented by differentiating the series term by term. A stable representation of the density gradient that converges well is obtained by defining the derivative as the difference between adjacent grid points, divided by twice the grid spacing as suggested by Lanczos. [13] This is equivalent to differentiating, term by term, the Lanczos-damped series for the charge density.

Symmetry
The set of internal forces acting on the atomic sites of a crystal is a symmetric, discrete function of atomic coordinates and has a spherical expansion on the crystal sites with the same coefficients as continuous symmetric functions (8a) and (8b). Since forces are vectors, their representation has ℓ = 1, and if a site has no invariant harmonics with ℓ = 1, there is no force on that site. So the force on an atomic site may be expressed as where the coefficients α are as in (8b), theê m are spherical unit vectors, [14] and U τ is the transformation to local coordinates for spherical vectors. A force calculation is, as much as possible, a calculation of the set {f ht }; The size of this set is often much smaller than three times the number of atoms. The displacements of atoms allowed by symmetry also have the form of (39): Minimizing the energy with respect to the atomic positions is a process of finding the set {δτ ht } that gives f ht = 0.

Force Calculations
The calculation of forces in an all-electron method has been nicely described by Yu et al. [15] for the LAPW method. In addition to the terms discussed in that paper, a force calculation using a site-centered basis has the additional, and significant, complication that the bases depend on atomic position not only through augmentation but also through parentage.
The contributions to the total force on a site in an all-electron calculation follow directly from a derivative of the LDA total energy with respect to atomic positions. The terms listed by Yu et al. are 1) a "Helmann-Feynman" term, ∂E/∂τ , which accounts for the explicit dependence of the energy functional on atomic positions, 2) an "Incomplete Basis Set" (IBS) term, which arises when derivatives of basis functions aren't contained in the space covered by the basis set, 3) a core-correction term, arising because core states are calculated using only the spherical average of the potential, and 4) a muffin-tin term, a surface term arising from the change in integration boundaries when atoms are moved and the discontinuity of the second derivative of basis functions across muffin-tin boundaries. There are two other terms to consider. The first arises when a calculation isn't fully self-consistent, and has the form − Vc (V out − V in )dn(r)/dτ , where V out and V in are output and  [16]. The LDA equilibrium volume is .958 Vexp; the GGA volume is 1.037 Vexp. The energy was also minimized with respect to the rhombohedral angle at each volume.
input potentials. The second term arises from the way in which Brillouin zone integrals are done. Whether by quadrature or linear interpolation, the result is a set of weights (occupations) multiplying quantities evaluated at discrete Brillouin zone points. The terms listed above do not take into account the change of weights with atomic positions. The evaluation of the IBS term in a method using site-centered bases is significantly more involved than in the LAPW method. This term has the form where the A are eigenvectors. Both LAPW and LMTO methods have a dependence on atomic positions through augmentation (the expansion of the basis set in atomic-like spherical waves) in the muffin-tin spheres, and both methods have an implicit dependence of basis functions on atomic positions through self-consistency, a term largely ignored and usually negligible. A sitecentered basis, however, depends on atomic positions also through it's parent site (the site it's centered on). The contribution from augmentation is fairly easily accounted for at the density stage of a calculation, after integrals over the Brillouin zone have been done. The parent contribution, however, requires evaluation at the part of the calculation where eigenvalues and vectors are obtained, which makes its calculation time consuming.
There are four types of contributions to dψ/dτ : wherep is the momentum operator −i∇. The first two terms, Equations (43) and (44), are parent terms, changes in a basis due to a change in the site the basis is centered on. The first term, Equation (43), is the derivative of the wave function in the interstitial region (Equation (10) with respect to its parent site. Since the gradient of a solution to the Helmholtz equation is a solution to the Helmholtz equation, matrix elements ψ i pψ j I and ψ i − ∇ 2 pψ j I are calculated as integrals over the surface of the muffin-tin spheres. As in Equation (22), when interstitial region tail parameters are the same, the evaluation requires κ 2 derivatives of structure functions. Working out this contribution proceeds as in Equation (22), although arriving at a finite form requires identities such as As I Fig. 2. Relaxation of a silicon 65 atom supercell containing a vacancy, a Si interstitial, and an As interstitial. Of the 106 internal coordinates in this cell, 104 were allowed to relax (2 coordinates were fixed to fix the center of mass of the crystal). The calculation used a simple Broyden's method to zero atomic forces.
Potential matrix elements ψ i V ψ j are calculated using Fourier series as in Sect. 3.2 with gradients taken as discussed after equation (38). The second term, equation (42), is the analog of the first term in the muffin-tin spheres; i.e., this term is the derivative of a basis with respect to its parent site evaluated in the muffin-tin spheres. This term requires the gradient with respect to atomic positions of the structure function B. This gradient is easily obtained from the structure function itself: If convergence with respect to ℓ on the left hand side of the structure function is sufficient for the energy, terms in ℓ max + 1 in Equation (48) may be neglected in evaluating forces. As stated above, the evaluation of these terms is somewhat time consuming. Examples of the use of forces for structural relaxation are given in Figures  1 and 2. Figure 1 shows deviations from ideal lattice positions calculated for rhombohedral BaTiO 3 as a function of volume compared to experiment [16]. The rhombohedral angle was also relaxed at each volume in this calculation. The Ti coordinate is a displacement along [111]. The oxygen displacements ∆x are along face diagonals while ∆z is toward the cell center. These calculations included Ti 3s and 3p and Ba 5s and 5p along with the usual valence bases in a single, fully hybridizing basis. At convergence, forces on internal coordinates were less than 1 mRy/Bohr. Figure 2 is a calculation of structural relaxation of As-vacancy-interstitial complex in Si. To a sixty-four atom Si supercell was added an As impurity at a tetrahedral interstitial position and a Si interstitial at an exchange position both surrounding a vacancy. The crystal, far from equilibrium, was then allowed to relax. Two internal coordinates (of a total of 106) were fixed to fix the center of mass of the crystal. The energy was minimized with respect to the other 104 internal coordinates by zeroing the forces (to with 1 mRy/Bohr). The forces were zeroed using a simple Broyden's method.

Conclusion
In this article we have described our highly accurate full-potential LMTO method for solving the Kohn-Sham equations. In particular, we have shown that by dividing the crystal space into non-overlapping "muffin-tin" spheres and an interstitial region, we can compute the charge density or the potential without any shape approximation, thus eliminating any need for empty spheres which are necessary in other LMTO implementations when the crystal is not closely packed. Another feature of our implementation is that we can describe multiple principle quantum numbers within a single, fully hybridized basis set. This is accomplished simply by using functions φ andφ calculated with energies {e nℓ } corresponding to different principal quantum numbers n to describe the radial dependence of a basis in the muffin-tin region. In the interstitial region our method uses "multiple κ" basis sets, for a better description of the interstitial charge density. Highly accurate charge density can be obtained by systematically increasing the number of variational parameters κ for each angular momentum of the basis set.
The potential in a muffin-tin sphere at τ has an expansion in linear combinations of spherical harmonics invariant under that part of the point group that leaves atomic positions invariant. The evaluation of the interstitial potential matrix only requires a correct treatment of basis functions (and potential) in the interstitial region. We have used this degree of freedom to design "pseudo-basis functions", equal to the true basis functions in the interstitial region and are smooth functions in the muffin-tin region, with the requirement that their Fourier transforms converge rapidly enough for practical use.
The set of internal forces acting on the atomic sites of a crystal is a symmetric, discrete function of atom coordinates and has a spherical expansion on the crystal sites with the same coefficients as continuous symmetric functions. The total force on a site is given by the derivative of the LDA total energy with respect to the atomic position. Our implementation of the forces is in many ways similar to that of Yu et al. for the LAPW method [15]. Because our basis set is a site-centered one, we are required to compute additional terms, which can be time consuming. These contributions to the forces are non existant in plane-wave based methods, such as the pseudo-potential method. In addition to the "Helmann-Feynman" term, which accounts for the explicit dependence of the energy functional A on atomic positions, the other contributions are: (1) an "Incomplete Basis Set" term, (2) a core-correction term, (3) a surface term arising from the change in integration boundaries when atoms are moved, (4) a term which arises when the calculation isn't fully self-consistent, and (5) a term arising from the way in which the Brillouin zone integrals are performed. We have showed that the forces are accurate enough to relax atomic structures. As examples, forces have been used to optimize the internal coordinates of rhombohedral BaTiO 3 as a function of volume and the geometry of a 65 atom As, vacancy, and interstitial defected Si supercell. Where experimental results are available, good agreement is obtained.