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; this is 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.

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” Pl. There is a one-to-one correspondence between Pl and for a fixed potential; but Pl 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”) 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”

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

Free electrons

In the absence of a potential the partial wave has the shape rl. Thus for free electrons, . This sets a reasonable lower bound to . For 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 Pl 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 Pl) 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 nth energy moment for partial wave of quantum number l is

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 ri 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 (φl2, φ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 momentsQ0, Q1, and Q2 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 Q0, Q1, Q2 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, Q0…2 in each l channel.

This is a generalization of the free-atom case where the atomic density is determined by the zeroth moment Q0 in each l channel and the boundary condition that φl decay as r→∞. Only Q0 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, Q0,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” Cl and the bandwidth Δl.

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

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 Q0,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 (Q0 = occupation of the free atom, Q1 = Q2 = 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 for augmented-wave programs is the choice of sphere radius. Results are much more sensitive to choice of spheres in the ASA than in the full-potential case, in part because the energy functional (and potential) change with MT radii, whereas in the FP case, they do only weakly so. The ASA also has the additional constraint that the sum-of-sphere volumes equals the unit cell volume, so the criteria in selecting them is somewhat different.

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.

Geometry violation of overlapping spheres

Overlapping spheres count some parts of space twice and others not at all. The full-potential code 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. Note: the current GW implemntation doesn’t have this property: there, spheres should not overlap.

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.

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.oij / dij is constrained to be less than OMAX1oij / min(ri,rj) 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.

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 El. In a linear method, the phase shift is parameterized:

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.

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 El. (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 El since then the eigenvalues are exact at El, 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 El. This is called folding down about 1/P(El).

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(El) 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 sp3s* basis). Now try folding down the empty sphere s as well: any worse than Harrison’s minimal basis? (Try the deformation potentials!)

Other Resources

For an overview of linear methods, and their connection to pseudopotentials, see these lecture notes given at a CECAM workshop at Oxford in 2013. The 2nd generation potential parameters this package uses are particularly helpful because they refer to conceptually accessible quantities, such as the band-center parameter C, and the bandwidth parameter Δ, as the lecture notes briefly describe.