The Atomic Spheres Approximation

Table of Contents

• Introduction
• The ASA Suite
• Structure of the ASA
• Augmentation sphere boundary conditions and continuous principal quantum numbers
  • Continuous principal quantum number for core levels and free electrons
  • Floating linearization energy
  • Numerics in solving the radial Schrodinger equation
• Generation of the sphere potential and energy moments Q
  • Potential Parameters
• Selection of Sphere Radii
  • Algorithm to automatically determine sphere radii
  • Autoscaling of sphere radii
  • Finding empty spheres
    • Assigning lower priority to resizing empty spheres
• 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:

1. The nonspherical part of the density and potential are neglected.
2. 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 $Q_0$, $Q_1$, and $Q_2$ of the density as described below, for each l. $Q_{0l}$, $Q_{1l}$, and $Q_{2l}$, with a boundary condition for the partial wave, usually specified by continuous principal quantum number $P_l$, 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 $\phi_l(\varepsilon,r)$ 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_l and $\varepsilon$ for a fixed potential; but P_l is more universal (it does not depend on constant potential shifts, for example) and is thus easier to work with.

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”) $\phi_l$ and its energy are uniquely determined by the boundary condition at the augmentation radius s. More precisely, $\phi_l$ 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 $\phi$ and its energy derivative $\dot\phi$ (or possibly $\phi$ at two different linearization energies.

The boundary condition can be given through the “logarithmic derivative function” $$D_l(\varepsilon) \equiv D\{\phi_l(\varepsilon)\} = \left({\frac{d\ln\phi_l(\varepsilon,r)}{d\ln r}} \right)_{s} = \left( {\frac{r}{\phi_l(\varepsilon,r)}} {\frac{d\phi_l(\varepsilon,r)}{dr}} \right)_{s} .$$

The energy $\varepsilon$ fixes $D$, or alternatively $D$ can be specified which fixes $\varepsilon$. $D_l$ is a cotangent-like function: it decreases monotonically from $(+\infty,-\infty)$ over a finite window of energy, after which it starts again at $+\infty$. There is thus a multiplicity of energies for a given $D_l$, 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 = 0.5 - \arctan(D_l)/\pi + \mathrm{(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 $P_l$ and the energy of the partial wave, ε_l.

Note: $P_l$ 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_l is approximately s/ε_l, which is large and negative.
Using the fact that arctan(x→-∞)/π→-1/2, the fractional part of P_l is large and close to one.

Free electrons

In the absence of a potential the partial wave has the shape r^l. Thus for free electrons, ${\mathrm{frac}}[P_l^{\mathrm{free}}] = 0.5 - \arctan(l)/\pi$. This sets a reasonable lower bound to $P_l$. For l=0, the $P_l^{\mathrm{free}}{=}1/2$ (fractional part) and for l=1,2,3,4,5 $P_l^{\mathrm{free}}$ 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 $\phi_l(\varepsilon,r)$ about a linearization energy $\varepsilon_\nu$ 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 $\varepsilon_\nu(l)$ (or the corresponding P_l) to the band “center of gravity” in the course of the self-consistency cycle. It estimates this by estimating $\varepsilon_\nu(l)$ that causes the first energy moment of the partial wave to vanish: $$Q_1 = \int_\mathrm{occ} w_l(\varepsilon) {\phi _l (\varepsilon ,r)(\varepsilon - \varepsilon _\nu )\phi _l (\varepsilon ,r)} d^3 {\mathbf{r}} = 0$$

$w_l$ is coefficient specifying how much $\phi_l(\varepsilon)$ contributes to the density at energy $\varepsilon$.

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 $$Q_{nl} = \int_\mathrm{occ} w_l(\varepsilon) {\phi _l (\varepsilon ,r)(\varepsilon - \varepsilon _\nu )^n\phi _l (\varepsilon ,r)} d^3 {\mathbf{r}} = 0$$

Of particular importance for the ASA are first three energy moments $Q_{0l}$, $Q_{1l}$, $Q_{2l}$.

Numerics in solving the radial Schrodinger 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², φ_l×(dφ_l/dε), (dφ_l /dε)²) 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₀, Q₁, and Q₂ 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₀, Q₁, Q₂ 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_0…2 in each l channel.

This is a generalization of the free-atom case where the atomic density is determined by the zeroth moment Q₀ in each l channel and the boundary condition that φ_l decay as r→∞. Only Q₀ 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 $\phi$ and $\dot\phi$ 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_l 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.

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₀ = occupation of the free atom, Q₁ = Q₂ = 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.

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 $i$ and $j$ is defined as $o_{ij} = (s_i + s_j - d_{ij})$ where $s_i$ is the augmentation radius for sphere $i$ and $d_{ij}$ the distance between sites $i$ and $j$. The constraints on $o_{ij}$ 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_ij is constrained to be less than OMAX1
  o_ij / min(r_i,r_j) is constrained to be less than 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.

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:

1. 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)

2. 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)

3. 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).

4. 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.