# The Atomic Spheres Approximation

### Table of Contents

- Introduction
- The ASA Suite
- Structure of the ASA
- Augmentation sphere boundary conditions and continuous principal quantum numbers
- Generation of the sphere potential and energy moments Q
- Selection of Sphere Radii
- Downfolding in the ASA
- Other Resources

### Introduction

The Questaal package has three codes that implement DFT in the Atomic Spheres Approximation (ASA). Formulated by O. K. Andersen in the 1970’s to handle transition metals, the ASA overlaps the augmentation spheres so that the interstitial volume is zero (there is a geometry violation). Moreover, the potential is assumed to be spherically symmetric inside the spheres.

The ASA is very efficient, but its range of validity is limited. This is because the interstitial is omitted so spheres must fill space. Hence, there is a geometry violation that becomes severe if the spheres overlap too much. It works best for close-packed systems, and still remains one of the best and most highly efficient approaches to studying magnetic properties of transition metals and their alloys. The ASA package has a non-collinear framework and a fully relativistic Dirac branch.

Questaal’s implementation uses the “tight-binding” form of LMTO, sometimes called “second generation,” a linear transformation” of the original basis set that makes Hankel functions short ranged.

**Note** There is also a non self-consistent implementation of Anderen’s most recent basis, the ‘NMTO’. This code should largely be regarded as experimental, as there are practical pitfalls associated with it that haven’t been fully worked out.

### The ASA Suite

Questaal has three implementations of the ASA:

*lm*: a band method whose function is similar to the full-potential program*lmf*. It is interesting to compare the ASA band structure to the FP one, e.g. in PbTe. To get started with this code, see the Introductory tutorial.*lmgf*: is a crystal code similar to*lm*, but it uses a Green’s function formalism. An extra energy integration (in addition to the k integration) is required, which makes the program somewhat slower. However it has features*lm*does not: it can calculate magnetic exchange interactions and some other properties of linear response. This code can include spin-orbit coupling perturbatively, and has a fully relativistic Dirac formulation. It also implements the Coherent Potential Approximation, either for the study of alloys, or for disordered local moments, or a combination of the two.*lmgf*has an Introductory tutorial.*lmpg*: is an analog of*lmgf*for layered systems. Periodic boundary conditions are used in two dimensions, and a Principal Layer technique is used for the third dimension. This is advantageous because (1) boundary conditions in this dimension semi-infinite leads, corresponding to layered systems and (2) the computation time scales only linearly in the number of principal layers. It can be used in a self-consistent framework, and also to calculate transmission using Landauer-Buttiker theory. There is a non-equilibrium Keldysh formulation of the ASA hamiltonian of the theory described in this paper.

### Structure of the ASA

The ASA is like other augmented wave methods which divide into an “atomic’’ part which makes matrix elements and a “band’’ part which generates bands, densities-of-states, etc. The ASA makes two simplifications to the atomic part that make the method highlty efficient:

- The
**nonspherical part**of the density and potential are neglected. - The spheres are
**overlapped so that they fill space**.

The net interstitial volume is zero, and in the pure ASA it is neglected all together.

Both atomic and band parts become simpler than in full potential methods. Matrix elements of the potential become quite simple and reduce to a few parameters (the potential parameters). The band part need only generate the lowest three energy moments , , and of the density as described below, for each *l*. , , and , with a boundary condition for the partial wave, usually specified by continuous principal quantum number , are sufficient for the atomic part to construct a density and make potential parameters. In the self-consistency cycle the atomic part takes moments and generates potential parameters; the band part takes potential parameters and generates moments.

Self-consistency proceeds by alternating between the solid part and atomic part, generating moments, then potential parameters, then moments again until the process converges. The program can be started either with the band part, specifying potential parameters, or with the atomic part, specifying the moments.

### Augmentation sphere boundary conditions and continuous principal quantum numbers

This section refers to the parameterization of a partial wave inside an augmentation sphere. It applies to all the all the codes in the Questaal package that integrate partial waves, not just the ASA. As described below, the energy-dependence is parameterized through a “continously varying principal quantum number” *P _{l}*. There is a one-to-one correspondence between

*P*and for a fixed potential; but

_{l}*P*is more universal (it does not depend on constant potential shifts, for example) and is thus easier to work with.

_{l}Linear augmented wave methods almost invariably construct the basis set inside augmentation spheres from the spherical part of the potential. (In the ASA the potential is always spherical, but even in the full-potential case, the partial waves that make up the basis set are constructed from the spherical part of the potential.) For a fixed spherical potential, the Schrödinger equation separates into an angular part (whose solutions are spherical harmonics) and a radial part with quantum number *l*. Solution to the radial Schrödinger equation (aka “partial wave”) and its energy are uniquely determined by the boundary condition at the augmentation radius *s*. More precisely, is called a *partial wave* since it is only a partial solution to the full Schrödinger equation. Partial waves must be matched to the envelope function at the augmentation sphere radius; the condition that all partial waves match smoothly and differentiably at all surfaces is the quantization condition that determines allowed eigenvalues. Linear methods in fact require the partial wave and its energy derivative (or possibly at two different linearization energies.

The boundary condition can be given through the “logarithmic derivative function”

The energy fixes , or alternatively can be specified which fixes . is a cotangent-like function: it decreases monotonically from over a finite window of energy, after which it starts again at . There is thus a multiplicity of energies for a given , one branch for each principal quantum number.

For that reason the Questaal suite uses a “logarithmic derivative parameter” or a “continuous principal quantum number”

*P _{l}* increases smoothly and monotically with energy, acquiring an extra integer each time a new node appears. This construction is due to Michael Methfessel.

*Note:*there is a one-to-one correspondence to and the energy of the partial wave,

*ε*.

_{l}*Note:* should not be confused with O.K. Andersen’s “Potential functions,” which *also* characterize, but in a different manner, the correspondence between partial wave energy and slope at the augmentation boundary. It is unfortunate that these distinct but related functions have the same symbol.

##### Continuous principal quantum number for core levels and free electrons

*Core levels*- A core state is exponentially decaying as it approaches
*s*; therefore its logarithmic derivative*D*is approximately_{l}*s*/*ε*, which is large and negative._{l}

Using the fact that arctan(*x*→-∞)/*π*→-1/2, the fractional part of*P*is large and close to one._{l} *Free electrons*- In the absence of a potential the partial wave has the shape
*r*. Thus for free electrons, . This sets a reasonable lower bound to . For^{l}*l*=0, the (fractional part) and for*l*=1,2,3,4,5 is 0.25, 0.15, 0.10, 0.08, and 0.06, respectively.

In summary, the fractional part of *P _{l}* is close to one for deep states far below the Fermi level, for states far above the Fermi level it is small, at least for

*l*>0. As

*ε*increases from −∞ to ∞,

*P*changes in a continous way, acquiring an extra integer each time a new node appears.

##### Floating linearization energy

The Taylor series in about a linearization energy is an approximation. Since the occupied states determine the density and thus the potential in density-functional and Hartree-Fock theories, it is important that this energy be suitably chosen to make the method as accurate as possible.

The Questaal codes will normally float each (or the corresponding *P _{l}*) to the band “center of gravity” in the course of the self-consistency cycle. It estimates this by estimating that causes the first energy moment of the partial wave to vanish:

is coefficient specifying how much contributes to the density at energy .

*Note:* cases can arise where the automatic floating algorithms cause problems. See this page for a discussion.

In general the *n*^{th} energy moment for partial wave of quantum number *l* is

Of particular importance for the ASA are first three energy moments , , .

##### Numerics in solving the radial Schrodinger equation

The Schrodinger equation in central potential reads

It can be written this way because is separable in spherical coordinates. Thus *l* and *m* are a good quantum numbers: the angular part of the Laplacian can be integrated out with spherical harmonics as solutions. In the Questaal suite we use real analogs of the spherical harmonics, as described here.

The radial Schrodinger equation remains to be integrated

The Questaal codes do not solve the Schrodinger equation, but a scalar (two-component) approximation to the Dirac equation. The Scalar-Dirac equation is integrated numerically on a radial mesh of points *i*=1…*N*. Since valence states must be orthogonal to core states, varies rapidly near the nucleus and a shifted logarithmic mesh is used. Point *i* is given by

Near the nucleus mesh points are linearly spaced by *ab*. For *r _{i}* large compared to

*a*, the mesh points are spaced exponentially (equally spaced on a log scale, spacing

*a*). The first point falls at the origin and the last at the augmentation radius

**R**.

Three numbers are used to specify the mesh: augmentation radius, the number of points inside it, and the spacing parameter *a* (*b* can be determined from them). These can be specified in the input file as **SPEC_ATOM_R**, **SPEC_ATOM_NR**, **SPEC_ATOM_A** though usually you can rely on default values.

The integration is performed outwards from the origin to some matching point, and inwards from **R** to the same point. Then (or more commonly the logarithmic derivative ) is varied until the discontinuity in the slope vanishes, and the partial wave has the correct number of nodes. This inward/outward integration technique is suitable for stiff differential equations. The radial Schrodinger is in fact stiff for states well localized inside **R** (deep valence states anre core states) since must decay exponentially when it is localized.

The ASA codes have the facility to integrate the actual four-component Dirac equation; *lmf* does not yet, but *lmfa* can generate the core using the true Dirac equation.

### Generation of the sphere potential and energy moments Q

Because the method is a linear one, and because the density is constrained to be spherical, only three functions can carry charge inside a sphere per *l* channel (*φ _{l}*

^{2},

*φ*×(

_{l}*dφ*/

_{l}*dε*), (

*dφ*/

_{l}*dε*)

^{2}) and therefore the properties of a sphere, for a spherical potential and a linear method are completely determined by the first three energy moments

*Q*

_{0},

*Q*

_{1}, and

*Q*

_{2}of the density of states for each

*l*channel, which are called the atomic number and the boundary conditions at the surface of the sphere. In some sense these numbers are ‘‘fundamental’’ to a sphere; the atomic program will generate a self-consistent potential for a specified set of

*Q*

_{0},

*Q*

_{1},

*Q*

_{2}and boundary conditions, specified in the Questaal package through the continuous principal quantum number

*P*described in the previous section. This simplification depends on assumption of spherical densities, and is unique to the ASA. Information specifying the potential is carried compactly in the four numbers

*P, Q*in each

_{0…2}*l*channel.

This is a generalization of the free-atom case where the atomic density is determined by the zeroth moment *Q*_{0} in each *l* channel and the boundary condition that *φ _{l}* decay as

*r*→∞. Only

*Q*

_{0}is needed because the atomic level is sharp, having no energy dispersion. Also the boundary condition is fixed by the requirement that

*φ*is integrable.

*Potential Parameters*

Once a potential is specified (implicitly through *P, Q*_{0,1,2} ), “potential parameters” can be generated from Wronskians of the partial waves and at the muffin tin boundary. They are a compact representation of information needed to specify the hamiltonian. A description of how the parameters are generated and their significance is too involved to be described in this overview, but see Other Resources. The most important parameters are the “band center of gravity” *C _{l}* and the bandwidth Δ

*.*

_{l}*C*describes the band center, and is the analog of the on-site matrix element (or atomic level in the free atom)_{l}- Δ
characterises the width of the partial wave, i.e. approximately the maximum and minimum values a partial wave would take in the absence of hybridization with other atoms._{l} - γ
is a “phase shift” parameter that is closely associated with the “screening transformation” to a tight-binding form._{l}

Another useful quantity is the “small parameter” *p*, which tells you the energy window over which the partial wave is well described by the linear method These parameters build the potential function which forms the basis of multiple scattering theory in Green’s function context (implmented in *lmgf* and in *lmpg*) and the energy-independent ASA hamiltonian in the bands context (implemented in program *lm*).

See the book by Turek at al for a thorough and readable published reference of the potential parameters, and how Green’s functions are constructed, and their connection to energy-independent hamiltonians suitable for band band calculations.

For a connection between *C*, Δ, and scattering phase shifts, see downfolding below.

To generate bands and an output charge density (in the form of moments *Q*_{0,1,2}), only potential parameters are required. Nevertheless it is more common to start from the moments because rough values for them can easily be guessed. The ASA codes will assume default values (*Q*_{0} = occupation of the free atom, *Q*_{1} = *Q*_{2} = 0), which most of the time is good enough to reach self-consistency. These codes also have a lookup table for default values of *P* described above.

### Selection of Sphere Radii

One of the biggest inconveniences as a result of partitioning space augmented-wave methods, is that you must select the partitioning (augmenation sphere radius).

The inconvenience is particularly pronounced in the Atomic Spheres Approximation, because the ASA takes only the spherical part of the potential inside each augmentation sphere; and also the ASA neglects the interstitial. This omission, because the system must conserve charge, implies that the spheres must fill space. The best we can do is make sum-of-sphere volumes equal the unit cell volume, but this requires that spheres overlap which means some points are double counted and others not at all. Such geometry violations can become severe unless sites are reasonably close packed. If they are not, you will need to alleviate the problem by including additional fictitious atoms “empty spheres” (fictitious atoms with atomic number zero), as was done in this tutorial. Because of the space-filling constraint, the criteria in selecting sphere radii is somewhat different from the full-potential case. If these extra sites are needed, Questaal has some algorithms to automatically help you locate these empty sites, as described below.

For either the ASA or FP, the Questaal package has several tools to help you select radii automatically.

- The input file maker,
*blm*, automatically selects them. Many tutorials, such as the basic*lmf*and*lm*tutorials, start with*blm*. - If you want to modify a
*ctrl*file you already have, the geometry checker*lmchk*will find radii (`lmchk --getwsr`

) - Questaal programs can rescale preselected sphere radii up to a specified volume within constraints you supply.

Finding sphere radii automatically is relatively straightforward in the FP case; the ASA can be tricky because of the space-filling requirement. See this page for some general considerations and click on the line below for observations mostly specific to the ASA.

- Geometry violation of overlapping spheres
- Overlapping spheres count some parts of space twice and others not at all. The ASA band code
*lm*, has a “combined correction” term that partially reverses this error, but not completely. The Green’s function codes*lmgf*and*lmpg*do not have this term. - ASA Requirement for space-filling spheres
- The ASA functional requires that the sum-of-sphere volumes equals the cell volume. More precisely, the density is carried by the spheres (superposition of spherically symmetrical sphere densities). This criterion mitigates directly against the preceding one. The more closely packed a system is the better suited the ASA. For open systems, you must add “empty spheres:” fictitious atoms of zero atomic number that enable space to be filled without too large a geometry violation.
- Overlapping spheres wth
*lmf* - The full-potential code
*lmf*has a unique augmentation, constructed so that the sphere contributions vanish quadratically for radii approaching the MT radius. Overlap errors tend to be small until overlaps reach about 10% of the internuclear distance. It has been found empirically, however, that self-consistency proceeds more slowly when spheres overlap. - Large sphere radii assign more volume to augmented functions
- Augmented wave functions are very accurate, and the more space covered by them the more reliable the basis set.
*l*-convergence is most rapid for small sphere radii- The larger the sphere radius, the slower the convergence with
*l*, because the angular momentum increases rapidly with*l*. - Larger spheres better contain shallow semicore states
- Ideally the core is completely localized within augmentation spheres. Particularly in the full-potential case where spheres overlap less than in the ASA, shallow semicore states can be an issue. In the FP case, you can always add a local orbital to address this problem.
- MT potentials are exactly solvable
- The KKR method is essentially exact for a MT potential, i.e. one that is spherical inside augmentation spheres and approximately flat in the interstitial. The LMTO basis starts from the KKR basis; thus a partitioning of space which best resembles a MT potential is the best choice. The automatic sphere radii algorithms try to select radii that make the intersitial potential flat.

#### Algorithm to automatically determine sphere radii

The ideal choice of sphere radii best approximates a potential that is spherical within the augmentation spheres and flat outside. *blm* and *lmchk* use an algorithm that makes a reasonable initial choice: they compute the (electrostatic) potential obtained from overlapping free-atom densities along all connecting vectors between a given site and its relatively near neighbors. The augmentation radius is taken as the first potential maximum along any ray. This choice is a pretty reasonable estimate for the potential being approximately spherical inside. Also, for a completely symmetric bond, the potential maximum will fall exactly midway between the bond, so for that case the two sphere radii will exactly touch and have equal potentials.

*blm*uses this algorithm automatically*lmchk*will run the algorithm and print the results to stdout if you invoke it with`lmchk --getwsr ...`

You must use your editor to copy each radius to the appropriate**SPEC_ATOM_RMAX**tag in the ctrl file.

There is at present no automatic way to incorporate radii generated by*lmchk*’s output into the input system.

#### Autoscaling of sphere radii

Questaal programs can scale sphere radii as large as possible within constraints you supply. This option iteratively adjusts sphere radii as large as it can within certain constraints, or until the aggregate sphere volumes equal the target you set. To autoscale sphere radii, set **SPEC_SCLWSR**=*f*, where *f* is the target aggregate sphere volumes as a fraction of the cell volume. For the ASA, *f* should be 1. You can use it for the FP codes too, but it usually isn’t necessary.

The overlap between spheres at sites and is defined as where is the augmentation radius for sphere and the distance between sites and . The constraints on come in the following flavors (all of them are imposed):

- Constraints on sphere overlaps
- There are constraints on sphere overlaps set through tags
**SPEC_OMAX1**and**SPEC_OMAX2**.

*o*/_{ij}*d*is constrained to be less than_{ij}**OMAX1**

*o*/ min(_{ij}*r*,_{i}*r*) is constrained to be less than_{j}**OMAX2** - Maximum sphere radius
- Cap the maximum sphere radius by setting
**SPEC_WSRMAX**

Lock sphere radii of specific species, by setting**SPEC_ATOM_CSTRMX**in any species you want to freeze.

#### Finding empty spheres

You can use “empty spheres” as a device to fill space. The ASA constraint that the aggregate sphere volume match the cell volume can only be realized with reasonable sphere overlaps (<18%) for fairly close-packed systems. The geometry violation for open systems, e.g. Si, is too severe. To use the ASA for such systems you must pack the volume with “atoms” with atomic number zero (“empty spheres”). This can be tedious, but *blm* and *lmchk* have an automatic “empty sphere” finder that can greatly facilitate this step. See this tutorial.

##### Assigning lower priority to resizing empty spheres

Particularly in the ASA, empty spheres are often needed to get reasonable sphere packing. However, it is reasonable that their radii should be scaled after the real spheres are rescaled. You can tell the resizer to do this through the 10’s digit of tag **SPEC_SCLWSR**. The 10’s digit behaves like a flag to cause the resizer to treat empty spheres on a different footing from real atoms.

- Add 10 to
**SPEC_SCLWSR**to initially scale real atoms (those with Z>0) first. The scaling is done using radii of size zero for all empty spheres. After this initial scaling, the resizer will proceed rescaling all the spheres. - Add 20 to
**SPEC_SCLWSR**is similar to adding 10. However, the final rescaling applies only to the empty spheres; the real atoms’ spheres change only in the first scaling, without reference to the empty spheres.

### Downfolding in the ASA

The *lm* and *lmgf* codes have procedures for constructing minimal basis by downfolding orbitals whose band center of gravity is far above the Fermi level.

Downfolding is a procedure for constructing minimal basis sets and for avoiding ghost bands. The best description is in Ole Andersen’s Varenna Notes (section 4.12), and for the stout-hearted, there is a full account in Lambrecht and Andersen, Phys. Rev. B, 34, 2439 (1986). It is implemented in the ASA including the “combined correction” term. We include in this documentation a plain TeX source file of notes explaining in some detail how downfolding is implemented into the *lm* code.

One way to look at the scheme is the following. When an electron encounters an atomic sphere, the scattering it experiences can be described in terms of its phase shift, *η _{l}*. The tangent of the phase shift is a property of the scattering potential and the angular momentum of the electron,

*l*, and is a function of energy. Some electrons are weakly scattered, while others (for example

*d*electrons in transition metals) may scatter strongly, especially when their energy is close to the resonant energy

*E*. In a linear method, the phase shift is parameterized:

_{l}so that one can construct an energy-independent hamiltonian. In LMTO, it is customary to use the *κ*=0 (*κ*^{2} is the envelope function energy) KKR phase shift in the following parameterization:

which is correct to second order in . Potential parameters *C* and Δ are readily identified by comparing the two equations: Δ is the width *W* of the resonance, and *C* is the band center. γ is the second order distortion parameter, which can be seen to add a constant background to the phase shift. In practice, one can also include third order terms using the small parameter *p*; see Varenna notes.

*Note:* In (1), *P _{l}* is indeed Andersen’s “potential function,”

*not*the “continuously varying principal quantum number”! 1/

*P*is the LMTO analogue of the tangent of the phase shift in multiple scattering and KKR theories.

For electrons that scatter only weakly, one can further approximate the hyperbola (1) by a linear function. This is exactly what happens if one throws away orbitals from the basis—one approximates 1/*P* for these “high” partial waves by a tangent drawn through the hyperbola (1) at the energy , which is the depth of a square well pseudopotential with the same scattering properties as the atomic sphere for energy *E _{l}*. (If the structure constants have been screened in these channels then the tangent goes through the potential parameter .)

The best possible way to treat weakly scattered electrons is to make the tangent at *E _{l}* since then the eigenvalues are exact at

*E*, and the wavefunctions are correct to zeroth order. The way to do this, is to change the representation of the hamiltonian before discarding the orbitals. The effect of using a representation beta is to shift the inverse potential function (1) rigidly by the amount beta. This is done by choosing a beta such that so that when the orbitals are subsequently discarded from the basis, this amounts to approximating their scattering by a linear tangent to 1/P at

_{l}*E*. This is called folding down about 1/P(

_{l}*E*).

_{l}If one merely wants to avoid ghost bands, then turn on **OPTIONS_ASA_ADNF** and keep an eye on which orbitals are being folded down. The automatic switch will choose to fold down about 1/P(*E _{l}*) or about the screening parameter depending on how weakly they are scattering. It is often useful to set the switches manually as the criteria in the automatic mode are set to cause no loss of accuracy. Very often one can fold down orbitals and save a lot of time with only a small error in the eigenvalues. Examples are

*p*orbitals in transition metals and

*d*orbitals in Al. In Fe, the

*f*orbitals must be folded down to avoid a ghost band. Another application is in constructing minimal basis sets. As an exercise try folding down orbitals in Si right down to

*s*and

*p*on the atoms and

*s*in the empty spheres (these are the analogues of Vogl’s

*sp*

^{3}s* basis). Now try folding down the empty sphere

*s*as well: any worse than Harrison’s minimal basis? (Try the deformation potentials!)

### Other Resources

Questaal’s standard methods paper has a detailed description of the ASA, summarizing decades of work by many authors, especially O. K. Andersen who invented linear methods. See:

Dimitar Pashov, Swagata Acharya, Walter R. L. Lambrecht, Jerome Jackson, Kirill D. Belashchenko, Athanasios Chantis, Francois Jamet, Mark van Schilfgaarde, *Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique*, Comp. Phys. Comm. **249**, 107065 (2020).

Other classic papers:

This classic paper established the framework for linear methods in band theory: O. K. Andersen, “Linear methods in band theory,” Phys. Rev. B12, 3060 (1975)

This paper lays out the framework for screening the LMTO basis into a tight-binding form: O. K. Andersen and O. Jepsen, “Explicit, First-Principles Tight-Binding Theory,” Phys. Rev. Lett. 53, 2571 (1984)

This book explains the ASA-Green’s function formalism, with some description of potential parameters: I. Turek et al., Electronic strucure of disordered alloys, surfaces and interfaces (Kluwer, Boston, 1996).

These papers lay the foundation for the Green’s function theory in the ASA: O. Gunnarsson, O. Jepsen, and O. K. Andersen, Phys. Rev B27, 7144 O. K. Andersen, Z. Pawlowska, and O. Jepsen, Phys. Rev. B 34, 5253 (1986)

These notes explain the “second generation” LMTO and the screening transformation

O.K. Andersen, A.V. Postnikov, and S. Yu. Savrasov, in “Applications of Multiple Scattering Theory to Materials Science,” eds. W.H. Butler, P.H. Dederichs, A. Gonis, and R.L. Weaver, MRS Symposia Proceedings No. 253 (Materials Research Society, Pittsburgh, 1992) pp 37-70.

O.K. Andersen, O. Jepsen and M. Sob, in Lecture Notes in Physics: Electronic Band Structure and Its Applications, eds. M. Yussouff (Springer-Verlag, Berlin, 1987).

These papers go beyond the 2nd generation LMTO, still within the ASA:

O.K. Andersen, O. Jepsen, and G. Krier in Lectures on Methods of Electronic Structure Calculations, edited by V. Kumar, O. K. Andersen, and A. Mookerjee (World Scientific Publishing Co., Singapore, 1994), pp. 63-124.

O.K. Andersen, C. Arcangeli, R.W. Tank, T. Saha-Dasgupta, G. Krier, O. Jepsen, and I. Dasgupta in Tight-Binding Approach to Computational Materials Science, Eds. L. Colombo, A. Gonis, P. Turchi, Mat. Res. Soc. Symp. Proc. Vol. 491 (Materials Research Society, Pittsburgh, 1998) pp 3-34.