# The ASA Crystal Green's function program lmgf

### Summary

This package implements the local spin-density-functional theory, in the Atomic Spheres Approximation using Green’s functions. The implementation is built into code **lmgf**, which plays approximately the same role as the LMTO-ASA band program **lm**. However it has functionality **lm** does not. It can:

- Calculate magnetic exchange interactions
- Calculate magnetic susceptibility (spin-spin, spin-orbit, orbit-orbit parts)
- Calculate properties of disordered materials, either chemically disordered or spin disorder from finite temperature, within the Coherent Potential Approximation [1], or CPA.
- Calculate the ASA static susceptibility at to help converge calculations to self-consistency.

Details about how the ASA works, and Green’s functions in particular can be found in references section of ASA overview; see in particular Turek’s book.

See Other Resources for tutorials and related programs.

### Table of Contents

- Summary
- Introduction
- Structure of Green’s function program
- Energy Contours, Potential Shifts and Determination of the Fermi Level
- GF specific input
- lmgf-specific command-line arguments
- Test cases and examples
- The Coherent Potential Approximation
- Spectral function calculations with lmgf
- Footnotes
- Other Resources
- References

### Introduction

The Green’s functions are constructed by approximating KKR multiple-scattering theory with an analytic potential function, described below. The approximation to KKR is essentially similar to the linear approximation employed in band methods such as LMTO and LAPW. It can be shown that this approximation is nearly equivalent to the LMTO Hamiltonian without the “combined correction” term. Implementation of the Green’s function code is accomplished through **lmgf**. **lmgf** plays approximately the same role as the LMTO-ASA band program **lm**: you can use **lmgf** to make a self-consistent density as you can do with **lm**. A unique potential is generated from energy moments , , and , in the same way as the **lm** code. **lmgf** is a Green’s function method: Green’s functions have less information than wave functions, so in one sense the things you can do with **lmgf** are more limited: you cannot make the bands directly, for example. However, **lmgf** enables you to do things you cannot with **lm**, as described at the beginning of this document.

#### Potential functions

The scattering properties of a sphere whose potential is spherically symmetric can be encapsulated in terms of a phase shift *η _{l}* of a wave scattering off the sphere. Each angular momentum

*l*has its own phase shift, and it depends on energy. Alternatively

*η*can be defined in terms of the potential function

_{l}Here and are Wronskians of Hankel and Bessel functions and with partial waves at the muffin-tin sphere radius. Hankel and Bessel functions are solutions to the Schrodinger equation in the flat interstitial part of a muffin-tin potential, and Wronskians are used to match a linear combination of and to the value and the slope of . making the wave continuous and differentiable.

By *linearizing* the partial wave , the energy-dependent Hamiltonian becomes energy-independent. The Wronskians become fixed numbers, and reduce to four types thos if *H* and *J* with and its energy derivative, . These Wronkians are usually expressed in terms of the physically more interpretable potential parameters, the most important of which are the band center *C _{l}* and bandwith

*Δ*, and also

_{l}This simplification offers a huge advantage, which is why linear methods are ubiquitous in electronic structure theory. In LMTO, it is customary to use the ( is the energy of the Hankel and Bessel envelope functions) KKR phase shift in the following parameterization:

, and are potential parameters. they are explained in detail in the book by Turek et al; see also the description of downfolding) calculated from the partial waves inside augmentation spheres at the linearization energy.

The “scattering path operator” , essentially the inverse of the Green’s function apart from a scaling, s given by and the structure matrix . The latter are are structure constants, independent of potential, that relate the expansion Hankel function centered around a remote site in Bessel functions around another site:

Eigenstate appear where there are poles in *G* or when

Finding poles must be done by an searching in energy (a nonlinear eigenvalue problem). But it is easily shown with linearization the eigenvalue condition becomes a linear algebraic eigenvalue problem, which can be solved by band methods. However, when thee potential is energy-dependent, as it is, e.g., in the CPA, the hamiltonian cannot be so so simplified.

Green’s function methods must resolve by energy so advantage is gained by linearization is much less (the KKR method is the LMTO Green’s function method without the linear approximation). But **lmgf** makes the linear approximation anyway, for consistency with the rest of the Questaal suite.

*Note:* The potential function is easily confused with the “continuous principal quantum number” which bears the same symbol and has a similar purpose. This is unfortunate; you have to infer which is meant from the context.

### Structure of Green’s function program

**lmgf** runs in much the same way as **lm**, at least in its primary mode, **MODE=1**. The band pass routine of **lm**, **bndasa.f**, generates the eigenvalues and eigenvectors, which can in turn generate the quantities of interest. **bndasa** is replaced by a Green’s function routine, **gfasa**. **gfasa** can generate output moments, DOS, density-matrix, etc., in the same way as **bndasa** does.

In contrast to band methods (implemented in **lm**) where the Hamiltonian is energy independent and all the bands are generated by diagonalizing H, Green’s functions are calculated for a specific energy; information is extracted from for a particular energy.

This fact highlights the strengths and weaknesses of a Green’s function approach. Energy-integrated properties such as the moments, must be obtained by integrating over energy. Calculating explicitly at a family of energies is more cumbersome than diaonalizing a Hamiltonian. On the other hand, Green’s function methods are naturally suited to contexts where the energy-dependence is needed anyway. CPA theory yields an energy-dependent potential; Green’s functions are a natural way to implement it. Similarly, noninteracting susceptibilities can be expressed as G×G (`×’ implies either convolution or product, depending on the space you are working in).

**lmgf** always loops over some energy contour; what contour you use depends on the context as described below. **gfasa** accumulates various kinds of data for each mesh point, such as the point’s contribution from energy moments used in an ASA self-consistent cycle. Finally, an estimate for the Fermi level is determined from a Pade approximation. If the original guess for is sufficiently close, the cycle is finished as in **lm**. If the estimate is too far off, a new energy mesh is taken and the process is repeated.

In addition to its primary mode (**MODE=1**), there are other modes, notably **MODE=10** and **MODE=11** for computing magnetic exchange interactions within ASA-linear response.

### Energy Contours, Potential Shifts and Determination of the Fermi Level

For energy-integrated properties (**MODE=1**), a very fine energy mesh would be required if the energy integration was carried out close to the real axis. It is much more efficient to deform the integration contour into an elliptical path in the complex plane, approaching the real axis only at the lower and upper integration limits.

To integrate quantities over occupied states, integration to the Fermi level *E _{F}* is required.

*E*is not known but must be fixed by charge neutrality. Thus must be guessed at and iteratively refined until the charge neutrality condition is satisfied.

_{F}**lmgf**does not vary

*E*; the user specifies it at the outset. Instead

_{F}**lmgf**looks for a global constant potential shift

**vconst**added to the one-particle hamiltonian that allows it to satisfy charge neutrality.

**vconst**must be found by an iterative search, as described in the subsection below. Both

**vconst**and

*E*are maintained in a file

_{F}*vshft.ext*. Inspect this file and you may find it unecessarily complicated; it’s because you can also use it to add site-dependent shifts.

**vconst**and

*E*are stored in the first line.

_{F}*vshft.ext*is also used by the layer Green’s function code

**lmpg**, which requires extra information about shifts on the left and right leads; this additional information is also stored in

*vshft*.

Metals and nonmetals are distinguished in that in the latter case, there is no DOS in the gap and therefore the Fermi level (or potential shift) cannot be specified precisely.

**Metal case** (set by **BZ_METAL=1**): once the - and energy-points are summed over and the deviation from charge neutrality is determined, the code will attempt to find the potential shift that fixes charge neutrality.

This tutorial offers a working example.

**Nonmetal case** (set by **BZ METAL=0**): **lmgf** will not attempt to shift the potential, or ensure charge neutrality. The user is cautioned to pay rather closer attention to deviations from charge neutrality. It can happen because of numerical integration errors, or because your assumed Fermi level does not fall within the gap. You can use **METAL=1** even if the material is a nonmetal, but be advised that this scheme is not foolproof: the vanishing DOS at *E _{F}* puts a strain on the iterative search algorithm.

##### Some details concerning how mode 1 works

For each energy point, the BZ integration is accomplished by routine in **gf/gfibz.f**, which loops over all irreducible points, generating the “scattering path operator” and the corresponding for all the points in the star of to generate a properly symmetrized . Within the ASA, second-generation LMTO, is converted to proper Green’s function *G*, by an energy scaling. The scaling is carried out in routine **gf/gfg2g.f**. Next the various integrated quantities sought are assembled (done by **gf/gfidos.f**). **vconst** needed satisfy charge neutrality is found by a nested iterative procedure explained in the following paragraph. The shift is stored in *vshft.ext*, as noted above; I/O to this file is handled by routine **subs/iovshf.f**.

The search for charge neutrality proceeds iteratively inside two nested loops. In the outer loop, sphere moments are integrated from for a given trial potential **vconst**, which yields deviation **Δq** from charge neutrality in the entire system. **vconst** is a global constant potential shift of the one-particle hamiltonian, and it must be adjusted until **Δq** is zero. (Alternatively the Fermi level may be adjusted keeping **vconst** fixed; see **shftef** in **GFOPTS** below). From various data (initially, density-of-states at *E _{F}*; later, data pairs (

**vconst**,

**Δq**) from prior iterations) are used to estimate a new

**vconst**. The procedure is iterated until the change in

**vconst**falls below a tolerance

**padtol**.

It is possible to find a new **vconst** using the outer loop alone. The outer loop is simply repeated for new values of **vconst** until tolerance is reached. This is expensive, however, and there is additionally an inner “Pade” loop to accelerate convergence. The diagonal elements of *G* are fit to a Pade polynomial, and the Pade approximation to *G* is used in the iterative procedure described above. Once the neutrality point is found within the Pade approximation, the outer loop is repeated until the change in **vconst** falls below **padtol**. After convergence the search terminates and the updated **vconst** is written to file *vshft.ext*.

There is a second (optional) tolerance **qtol**, that modifies how the search terminates. If the **qtol** tag is missing or its value is zero, this tag has no effect. Otherwise when **Δq < qtol** the internal search for **vconst** ends. Specifically:

- A test is made before the inner loop starts. If
**Δq < qtol**the search terminates without any Pade correction. - Otherwise the search proeeds with inner loop (i.e. using a Pade estimate for
*G*). If**Δq < qtol**for the Pade estimate for*G*, the inner loop terminates and the search reverts to the outer loop.

It finds the Fermi level in one of three ways:

User has specified a nonzero

**qtol**, and also**Δq < qtol**for*G*calculated on the elliptical contour.*G*is accepted as is, and the search terminates.Using a Pade approximant,

**lmgf**interpolates the diagonal elements of G. The interpolation is used to estimate*G*on the starting elliptical contour shifted rigidly by**vconst**, and the shift is iterated until the charge-neutrality condition is satisfied. At this stage, there are two possibilities:1. repeat the integration of over and the energy contour with a new trial

**vconst**. This happens when the change in**vconst**exceeds**padtol**.2. Assume that the Pade-approximant to the diagonal

*G*is a sufficiently good estimate for the actual*G*. If the change in**vconst**is less than**padtol**, the search terminates and the last Pade approximation is taken for the diagonal part of*G*.The contour is continued on a set of uniformly spaced points on the real axis starting from the Fermi level. A trapezoidal rule is used (or Simpson’s rule using a Pade approximate for the midpoint), Integrals for the moments are accumulated until charge neutrality is found. There is no iterative scheme as with the Pade approximation. This option tends to be a little less accurate than the Pade, but somewhat more stable as it is less susceptible to interpolation errors.

One last comment about **vconst**: by default the program will save the potential shift to use in the next iteration. You can suppress this save (see **frzvc**), which again can be less accurate, but more stable. In particular, if you are working with an insulator where stability can be an issue (determination of the Fermi level is somewhat ill conditioned), a stable procedure is to use this option together with second energy integration scheme described above (the integration contour on the real axis).

### GF specific input

#### Energy integration

Green’s functions are always performed on some energy contour, which is discretized into a mesh of points in the complex energy plane. (A description of the various kinds of contours this code uses is documented in the comments to *gf/emesh.f*.) is “spikey” for energies on the real axis (it has poles where there are eigenstates).

##### Energy contours and the EMESH token

*Note:* the **EMESH** tag applies to both **lmgf** and **lmpg**.

To compute energy-integrated properties such as magnetic moments or the static susceptibility, the calculation is most efficiently done by deforming the contour in an ellipse in the complex plane.

At other times you want properties on the real axis, e.g. density-of-states or spectral functions. You specify the contour in category **BZ** as:

```
EMESH= nz mode emin emax [other args which depend on mode]
```

**nz**number of energy points**emin,emax**the energy window (emax is usually the Fermi level when integrating over the BZ)**mode**specifies the kind of contour:**mode=10**: a Gaussian quadrature on an ellipse.- This is the standard mode for integrating over the occupied states.

`EMESH= nz 10 emin emax ecc eps`

**ecc**is the eccentricity of the ellipse, ranging from 0 (circle) to 1 (line)**eps**is a ‘bunching’ parameter that, as made larger, tends to bunch points near**emax**.

As a rule,**eps=0**is good, or maybe**eps=0.5**to emphasize points near the Fermi level.

After the integration is completed, there will be some deviation from charge neutrality, because

**emax**will not exactly correspond to the Fermi level. This deviation is ignored if**METAL=0**; otherwise**lmgf**or**lmpg**has to go through an iterative procedue to find the charge neutrality point. There is a tutorial demonstrating the process.**mode=110**: is a contour specific to nonequilibrium principal layer Green’s function calculations.- Use this mode for iterating to self-consistency in the non-equilibrium case.

The nonequilibrium Green’s function requires additional information for the energy window between the left and right leads. (The nonequilibrium Green’s function is implemented for the layer geometry in

**lmpg**.) Thus the integration proceeds in two parts: first an integration on an elliptical path is taken to the left Fermi level (as in**mode=10**). Then an integration over is performed on the nonequilibrium contour, i.e. the energy window from the left to the right Fermi level. This integration is performed on a uniform mesh close to the real axis, as in**mode=0**. For the nonequilibrium contour, three additional pieces of information must be supplied:**nzne**number of (uniformly spaced energy points on the nonequilibrium contour**vne**difference in fermi energies of right and left leads,*E*(R)minus;E_{F}_{F}</i>(L)**delne**Im-z on the nonequilibrium contour

Parameters are specified as

`EMESH= nz 110 emin ef(L) ecc eps nzne vne delne [delend]`

The last argument

**delend**plays the role of**delne**specifically for computing the self-energy of the left and right leads. There is an incompatibility in the requirements for in the central and end regions. See**Footnote 1**.**mode=0**: a uniformly space mesh of points between**emin**and**emax**, with a constant imaginary component.- Use this mode when you want quantities resolved by energy on the real axis, e.g. spectral functions or transmission in the
**lmpg**case. **mode=1**: is the same contour as**mode=0.**- The only difference is that the sign of Im z is flipped.
**mode=2**: is the same contour as**mode=0.**- The only difference is that weights for each energy are set to unity instead of the spacing between energy points. Whether to choose mode 0 or mode 2 depends on if the energy-resolved property of interest, .e.g. density-of-states, is to be given as is, or to be weighted for an integration over energy.

`EMESH= nz 0|1|2 emin emax Im-z [... + possible args for layer geometry.]`

Im-z is the (constant) imaginary component.

For the Principal layer code

**lmpg**, use one of:`EMESH= nz 0 emin emax delta xx xx xx xx delend [equilibrium] EMESH= nz 110 emin ef(L) ecc eps nzne vne delne delend [non-equilibrium]`

**emin,emax**the energy window (emax is usually the Fermi level when integrating over the BZ)**delta**is Im*z*for the central region**xx**has no meaning but are present for compatibility with the contour used in nonequilibrium calculations.**nzne**,**vne**,**delne**are for nonequilibrium calculations (see mode**110**).**delend**is used to make the self-energy of the leads; see**Footnote 1**.

**mode=310**: Alternative to Pade approximant in finding the Fermi level.- Acts like
**mode 10**for the elliptical contour to**emax**, then switches to a uniform mesh to fine the Fermi level.

This mode integrates with a Gaussian quadrature on an ellipse to a trial

**emax**, as in**mode 10**. However, the search for the Fermi level is not done by Pade approximant, as in**mode 10**. Instead, a second integration proceeds along a uniform mesh from emax to some (Fermi) energy which satisfies charge neutrality. This procedure is not iterative.`EMESH= nz 310 emin emax e1 e2 delz`

+e1 and e2 are just as in mode 10 +delz is the spacing between energy points for the second integration on the uniform mesh.

#### Green’s function category

**lmgf** requires a GF-specific category.

```
GF MODE=1 GFOPTS=options
```

##### The GF_MODE token

**MODE= n** controls what

**lmgf**calculates. Options are

**MODE=1**,

**MODE=10**,

**MODE=11**,

**MODE=26**, described below.

**MODE=1** goes through the self-consistency cycle, calling **gfasa**. It performs a function analogous to **bndasa** in the band program, generating output density, moments, and other quantities such as density-of-states.

Taken with the special integration contour **mode=2** (see **EMESH** above), the density-of-states and its integral are computed and tabulated over the window specified.

**MODE=10** invokes a special branch that computes magnetic exchange interactions using a linear response technique. (The source code has its entry point in **gf/exasa.f**.)

In particular, is computed for pairs of sites , where the *J*’s are the parameters in the Heisenberg Hamiltonian

Thus, the *J*’s are coefficients to energy changes for small rotations of the spins. They can be computed from a change in the band energy; changes from small rotations are done analytically.

Taken with the usual elliptical integration contour, the *J*’s are computed by energy integration to the Fermi level. Taken with the special integration contour **mode=2** (see EMESH above), **dJ/dE** is computed instead. There is a shell script

```
gf/test/getJq0z
```

(invoke with no arguments to see usage) that will collect some of the ouput for you into tables. The data are collected into file **dj0dz**. For an example illustrating this mode, invoke

```
gf/test/test.gf co 5
```

This test computes the exchange coupling both for the usual elliptical contour and resolves the energy-dependence of *J* in a small window near the Fermi level.

Often only some atoms are magnetic, and all that is desired are the exchange parameters *J* connecting a partial list of sites to its neighbors. This can be useful, even essential, for large systems because it can be very expensive both in time and memory to compute exchange interactions for all pairs. To compute exchanges only for a list of sites, use command-line argument

```
--sites:pair:site-list
```

For more details, see command-line arguments invoked with **lmgf**.

**Caution**. **lmgf** reads and writes a potential shift file *vshft.ext* which shifts site potential by a constant to cause the Fermi level to match what is specified by the input. This shift also gets added into the atom file; potential **VES** in line **PPAR** is adjusted. When calculating exchange interactions, *vshft.ext* is not read. However, the shift is preserved because they are held in the potential parameters section of the atom file. But if you run the atom part **lm** or **lmgf** and remake the PP’s from the moments (**START BEGMOM=1**), this causes estat potential to be remade, but the sphere program does not add the contents **vshft** (it is done at the start of the Green’s function calculation). The exchange parameters should be evaluated with the potential parameters generated by **lmgf**. If they are alternatively evaluated from the atom files generated by **lm**, the Fermi level needs to be aligned to the Fermi level of **lm** (or close to it; there are slight differences between Fermi levels generated by **lm** and by **lmgf**).

**MODE=11** is an exchange branch that is run after **MODE=10**. It prints out the and does several other analyses.

Switch `--sites:pair:site-list`

applies to mode 11 as well as mode 10; see command-line arguments.

##### The GF_GFOPTS token

The **GF** category has a token **GFOPTS=tag;tag;…**, which causes **lmgf** to perform a variety of special purpose functions.

Options are entered as a sequence of tags delimited by a semicolons: **tag1;tag2;…** .

Tag | Purpose |

emom | generate the output ASA moments, needed for self-consistency |

noemom | suppress generation of the output ASA moments |

idos | make integrated properties, such as the sum of one-electron energies |

noidos | reverse of idos |

pdos | Make the partial density of states (this has not been checked recently). |

dmat | make the density-matrix G_{RL,R’L’} |

sdmat | make the site-diagonal density-matrix G_{RL,RL’}. The density matrix is written to dmat.ext. |

p1 | Use first order potential functions (rarely used) |

p3 | Use third order potential functions |

pz | Exact potential functions like KKR (some unreseolved problems; not recommended). |

shftef | Find charge neutrality point by shifting the Fermi level, rather than adding a constant potential shift |

frzvc | Suppress saving the constant potential shift used to determine charge neutrality Δq. |

padtol=# | Set the tolerance for maximum potential shift permissible by Pade interpolation, as described above |

qtol=# | Set the tolerance for deviation from charge neutrality. If tag is missing or # is zero, this switch has no effect. |

Otherwise qtol is used as another check for deviation from charge neutrality Δq. When Δq < #, code exits search for neutral point. | |

If search falls within a Pade interpolation, search reverts to the outer iterative search calculating the full G. | |

If search precedes Pade interpolation, the search terminates and the existing combination of Fermi level/potential shift is accepted. |

The following are specific to the the layer code

nclead | Allow leads to be noncollinear |

declead | Calculate g in the leads with decimation (this is usually the default) |

sreslead | Resolves transmission by spin components in the leads. Only applies if spins in the lead are coupled. |

refinegs=# | Use embedding to iteratively refine the surface g in the lead, after it has been made by decimation. 0 => no refinement |

The following are specific to the CPA:

omgtol | Tolerance in the Omega potential, CPA self-consistency |

omgmix | How much of prior iterations to mix, CPA self-consistency |

nitmax | Maximum number of iterations for CPA self-consistency |

lotf | Learn on-the-fly |

specfun | Make spectral function |

specfr | Make spectral function resolved by site |

specfrl | Make spectral function resolved by site and l |

dmsv | Record density matrix to file |

dz | Shift Omega potential by dz |

sfrot | Rotate the GF in spin space by this angle around the y axis |

refinegs=# | Use an iterative procedure to refine surface g after it has been made |

### lmgf-specific command-line arguments

```
-ef=# overrides upper limit of energy integration (Fermi level) and assigns to #
```

The following are specific to the exchange calculation **modes 10** and **11**:

```
--sites[:pair]:_site-list_ Make the exchange parameters J_ij only for sites in _site-list_, which is a standard [Questaal integer list](/docs/numerics/integerlists).
```

*Example*: running **lmgf** using **MODE=10** with this command line argument

```
--sites:pair:1,3,5,7
```

generates *J* connecting sites 1, 3, 5 and 7 to all neighbors.

This switch also works in **mode 11**, but the meeaning is different: **lmgf** will print out the exchange interacations only from between in

Running **lmgf** using **MODE=11** with the same **--sites** argument will print out the exchanges just between pairs of these sites.

Running **lmgf** using **MODE=11** without any **--sites** argument will print out the exchanges between these sites and all neighbors.

```
--wrsj[:fn=name][:scl=#][:tol=#] (mode 11 only)
Writes the Heisenberg exchange parameters in a standard format, suitable for use
in spin dynamics simulations.
fn=name writes to file 'name' (default name is rsj)
scl=# scales the parameters by #
tol=# writes only parameters with energy > tol
--rcut=#
Truncates the range of the R.S exchange parameters ...
useful to assist in the determination of the effect distant neighbors.
--2xmsh
When integrating over the BZ to estimate Tc from Tablikov formula, this option doubles the k-mesh.
Can be helpful in testing k-convergence of the singular q->0 limit entering into the formula.
--amoms=mom1,mom2,...
--amom=mom1,mom2,...
This switch overrides ASA moments (which are automatically generated).
The first switch reads a vector of nbas moments, one for each site.
The first switch reads a vector of nclass moments, one for each class.
```

Sphere magnetic moments are tabulated in the printout at the end of **mode 10**, and the start of **mode 11**. If you are importing exchange parameters (file *jr.ext* , e.g. from the full-potential code, you will want to supply the moments calculated from that program.)

### Test cases and examples

This script:

```
gf/test/test.gf --all
```

carries out a number of tests, which also demonstrate various branches of the code. To see the materials and corresponding tests try

```
gf/test/test.gf --list
```

### The Coherent Potential Approximation

The CPA implementation for substitutional alloys and for spin disorder follows the formulation explained in these References [1,2,3]. Particularly, see the description of the numerical implementation in Turek’s book.

CPA self-consistency is based on iterating the coherent interactor Ω, which is a spin-dependent single-site matrix defined for each CPA site at each complex energy point. The linear mixing of Ω can be interleaved with charge mixing steps. However, experience shows that much faster convergence can be achieved by iterating Ω at each z-point until its misfit reaches a sufficiently low tolerance (say, 1d-3), between charge mixing steps. In addition, it is better to skip charge mixing if sufficiently accurate charge-neutrality has not been achieved (the reason being that Ω is not Pade-adjusted). There are a few parameters controlling Ω convergence, which are summarized below along with the recommended settings that work quite well in most cases. The Ω matrices are recorded in files **omegaN.ext**, where *N* is the number of the CPA site. A human-readable version (with fewer decimal digits) is recorded in **om-hrN.ext**.

#### CPA-specific input

To turn on chemical and/or magnetic CPA, additons are required to the **SPEC** and **GF** categories in the ctrl file.

##### SPEC category

**Chemical Disorder**. Additional species must be defined for chemical CPA, and their concentrations.

```
SPEC ATOM CPA= and C= together turn on chemical CPA for a particular species.
```

They specify which species are to be alloyed with this species, and the concentrations of the other species. For example,

```
SPEC ATOM=Fe ... CPA=1 4 5 C=0.5 0.3 0.2
```

specifies that species Fe (whenever it appears in the basis (defined in **SITE** category) in fact refers to a disordered site composed of three kinds of elements. Numbers following **CPA=** refer to indices in the **SPEC** category: thus “**CPA=1 4 5**” indicate that the three elements to be identified with sites referring to this species are the 1st, 4th, and 5th species declared in the **SPEC** category. **C=** indicates the concentrations of each species; the concentrations must sum to 1. In the example given, sites with species label Fe are composite elements with with 50% of species 1, 30% of species 4 and 20% of species 5 (up to 10 species may be given).

A CPA species may refer to itself. For example, if the Fe species above is the first species to be read from the ctrl file, then **CPA=1** refers to itself. All other parameters like *Z*, *R*, will be taken from this species.

**Spin Disorder**. No additional species are required, but the number of orientations must be specified.

```
SPEC ATOM NTHET= turns on spin disorder for a particular atom type.
```

A species with non-zero **NTHET** can be listed as a CPA component, and it will be included as NTHET components with different directions of the local moment.

**NTHET=2** specifies that there will be two CPA-DLM components with polar angles 0 and π. **NTHET=N** with specifies a vector-DLM model, for which N polar angles for the local moment direction are selected using the Gaussian quadrature for the sphere. (Axial symmetry is always assumed and the integral over the azimuthal angle is taken analytically.)

**Combined Chemical and Spin Disorder**. Either spin or chemical disorder may be specified; they may also be included simultaneously. If only **CPA=** is chosen, that species will be treated with chemical, not spin, disorder. If only **NTHET=** is chosen, that species be treated with spin disorder only. Specifying both means that the CPA will include both chemical and spin disorder. For example, in the above example for CPA, if **SPEC ATOM=Fe** includes a tag **NTHET=2** (while species 4 and 5 have **NTHET=0**), species Fe describes a CPA site with 4 components: 25% Fe, 25% Fe with a reversed local moment, 30% species 4 and 20% species 5.

##### GF category

The following token turns on the CPA and/or DLM:

```
GF DLM= controls what is being calculated.
```

At present, these values are supported:

```
DLM=12: normal CPA calculation; both charges and Ω's are iterated
DLM=32: no charge self-consistency; only CPA it iterated until Ω reaches
prescribed tolerance for each z-point.
DLM=112: special-purpose experimental branch (not documented)
```

The following are optional inputs:

```
GF BXY=1 turns on the self-consistent determination of the
constraining fields for vector DLM calculations.
GF TEMP= supplies the spin temperature (not implemented yet)
```

Self-consistency in Ω is controlled by the following tags supplied in GF GFOPTS:

```
lotf if present, Ω is iterated at each z-point until converged to omgtol (recommended)
nitmax= maximum number of Ω iterations (30 is usually sufficient)
omgmix= linear mixing parameter for Ω (0.4 works well in most cases)
omgtol= tolerance for Ω
padtol= same meaning as usual, but note that Ω is not mixed unless padtol is reached
(1d-3 is recommended for all CPA calculations)
dz= special branch, in which z-points are shifted by dz along the real axis (experimental)
```

Recommended options:

```
GF GFOPTS=[...];omgmix=0.4;padtol=1d-3;omgtol=1d-3;lotf;nitmax=30
```

#### Compatibility with other features

Downfolding is supported. Note, however, that downfolding applies to the crystal Green’s function and not to individual CPA components. The downfolding options are taken from the first species appearing in the CPA list. Gamma representation is supported with a caveat. CPA does not allow random structure constants, which means that the screening parameters must be the same for all components on the same CPA site. In the present implementation, the screening parameters are taken from the first class listed for the given CPA site (for a DLM site this is angle #1).

*LDA+U* **is not** supported, and density matrices are not calculated for the components on the CPA site. However, the modes **IDU=4** and **IDU=5** are supported. The U and J parameters for these modes are taken from the first species appearing in the CPA list.

Broyden mixing for charged works fine if omgtol is set to a sufficiently low value. If Broyden mixing seems to act strangely, try to reduce omgtol. Charge self-consistency in CPA may sometimes be difficult for impurities with low concentrations. (Note that an isolated impurity can be described by adding it as a CPA component with zero concentration.)

#### Atomic files

It is important to understand the atomic file handling with CPA. For a CPA site (say, species Fe) the code creates an atomic file per each CPA component. In the above example with **SPEC ATOM=Fe … NTHET=2 CPA=1 4 5** there will be four atomic files: **fe#1.ext** for Fe with angle 0, **fe#2.ext** for Fe with angle π, **fe#3.ext** for species type 4, and **fe#4.ext** for species type 5. Note that fe#3 and **fe#4** will not actually correspond to Fe atoms, but to those described by species 4 and 5. Because convergence can be delicate, it is always recommended to copy appropriately prepared atomic files before attempting a CPA calculation. In the above example, converge a Fe atom and copy the atom file to **fe#1.ex**t and **fe#2.ext**; then converge species type 4 and copy it to **fe#3.ext**, and so on. For DLM with **NTHET=N**, make N copies of the atomic file: say, **fe#1.ext**, **fe#2.ext**, …, **fe#N.ext**.

#### Outputs

At the beginning of the run, some debugging information is printed, listing the indexing for the CPA sites. **DLMWGTS** lists the polar angles (0 for non-DLM classes) and weights for all CPA classes (this is also for debugging purposes). **GETZV** prints the total valence charge, which in CPA is generally not integer. Output for each CPA component includes the usual information (charge, local moment, etc.). Exchange constants J0 are automatically calculated for all CPA components using the linear response formula from Liechtenstein et al. (it can not be disabled, but the computational cost in any case negligible). Off-diagonal local moments and constraining fields are always printed out, even if DLM is not used. These include the diagonal local moment as well. All these moments are output, unmixed values. In the self-consistent state the z-component should equal to the input moment.

At the end of the iteration Ω is mixed, and its misfit for each CPA site is printed out (see “*Mixed Omega for site …*”) The total energy is correctly calculated and printed out as ehk, as usual.

#### Partial densities of states

Partial DOS can be calculated as usual using contour type 2 and adding pdos to the **GFOPTS** tag. Note that in this case Ω needs to be converged anew at each point of the new contour. This destroys the old converged Ω file, so it is recommended to create a separate directory for a DOS calculation. The file **dos.ext** contains the usual information, but the data for CPA sites are averaged over components. The partial DOS for all components are separately recorded in files **dosN.ext**, where N is the number of the CPA site. The format of this file is the same as that of dos.ext, as if it described a system with M sites (where M is the number of CPA components). For example, for a binary CPA on site 2 with spd basis, file **dos2.ext** contains channels 1:6 for the first CPA component and channels 7:12 for the second CPA component. This file can be processed using **pldos**, as a conventional dos file.

#### Spectral Functions

**lmgf** can generate spectral functions. It is very useful way to see the broadening of states from disorder, and you can plot energy bands with it. This document explains how to make them and draw energy bands.

#### CPA Test case

To familiarize yourself with a CPA case you go through the CPA tutorial.

You can also run the following test case

```
~/lm/gf/test/test.gf fe2b
```

### Spectral function calculations with lmgf

Spectral function have been implemented in **lmgf** **v7.10** by *Bhalchandra Pujari* and *Kirill Belashchenko* (*belashchenko@unl.edu*). The details of the theory in the CPA case can be found in Ref [1].

#### How to calculate spectral functions?

The spectral function can be calculated both with and without CPA. The calculation is performed in 3 steps:

- Charge self consistency. (The spectral function can be calculated for any potential, but it is usual to work with the self-consistent one).
- CPA self-consistency in the coherent interactor (CPA only). Since is energy dependent, it has to be calculated for the energy points where the spectral function is needed. For drawing spectral functions this is usually a uniform mesh of points close to the real axis.
- Calculation of the spectral function on some contour, usually a uniform mesh close to the real axis.

**Charge self consistency** is performed in the usual manner, for example with the following options:

```
BZ EMESH=31 10 -.9 0 .5 .0
GF MODE=1 DLM=12 GFOPTS=p3;omgmix=1.0;padtol=1d-3;omgtol=1d-5;lotf;nitmax=50
```

Note the **EMESH mode** (contour type) is elliptical (type 10). If CPA is used, the coherent interactors for all CPA sites are also iterated to self-consistency during this calculation, but this is done for the complex energy points on the elliptical contour. The following additional step is needed in this case to obtain self-consistent at those points where the spectral function will be calculated.

**Omega self-consistency** is turned on by setting **DLM=32**. In this mode only Ω for each CPA site is converged, while the atomic charges are left unchanged. It is important to converge Ω to high precision. Typically **omgtol=1d-6** is a good criterion. The contour type should be set to 2. Example input for this step is

```
BZ EMESH=150 2 -.25 .25 .0005 0
GF MODE=1 DLM=32 GFOPTS=p3;omgmix=1.0;padtol=1d-3;omgtol=1d-6;lotf;nitmax=50
```

The highlighted parameters are of particular importance. lotf is required to iterate for convergence (and it is recommended to keep it enabled in all CPA calculations, including charge self-consistency). It is also necessary to monitor the output file (set **--pr41**) and make sure that the required precision has been achieved for all energy points. If convergence appears to be problematic, try to start with a larger imaginary part for the complex energy or reduce the mixing parameter omgmix.

Calculation of the spectral function should be done with **EMESH** set to the same mesh as used for self-consistency, e.g.

```
BZ EMESH=150 2 -.25 .25 .0005 0
GF MODE=1 DLM=12 GFOPTS=p3;omgmix=1.0;padtol=1d-3;specfun
```

**Important: If there are sites treated in CPA, the contour specified by EMESH should be kept exactly as in the previous step when CPA self-consistency was performed.**

In order to start the calculation, invoke **lmgf** with the **--band** flag referring to the symmetry-line file (same format as used for band structure calculation with **lm**):

```
lmgf «sys» --band:fn=syml
```

where *«sys»* is the extension of the **ctrl** file. Once completed, the program will generate a **spf.«sys»** file containing the complete spectral function along the lines given in the **syml.«sys»** file. Other options included with **--band** are currently not used.

*Note:* Spectral function calculations can be run with MPI parallelization.

#### Plotting the spectral function

Upon successful completion of calculation **spf.«sys»** file will be generated. The file has the following format:

```
ZP KP SF_UP SF_DN
.
.
-0.250 0.137 5.708839 4.742153
-0.250 0.160 10.658114 4.844647
.
.
```

where *ZP* are the points on the energy contour, *KP* are the points of the K-mesh and SF_UP, SF_DN are the spectral functions for Up and Down channel respectively. The very first line of the file indicate the location of the high symmetry points of the Brillouin zone. User can utilize this information to visualize the spectral function using any desired graphics package. A small bash utility, **SpectralFunction.sh**, is given for the sake of convenience. This bash script uses Gnuplot to view and save the spectral function. Users with standard Linux/unix distro should be able to use it without special prerequisites. Typical output of the script is shown in the figure. Please see SpectralFunction.sh -h for usage. (The file **spf.«sys»** should be renamed **specfun.«sys»** to be read by this script. See below why.)

#### Site-resolved spectral functions

Site-resolved spectral functions can be obtained by including the option **specfr** in **GFOPTS**. (It supersedes the option specfun.) A series of output files will then be generated, which contain the spectral function on each basis site . The file names are **spfN.«sys»**, where N is the index of the site. The format of this file is the same as for **spf.«sys»**, and it can be also plotted using the SpectralFunction.sh script (first rename to **specfun.«sys»**).

### Footnotes

^{1} When computing Green’s functions near the real axis via the Landauer-Buttiker formalism for transmission through the active region, or the nonequilibrium part of the contour in nonequilibrium calculatations, or in special modes that search for the Fermi energy by integrating points on the real axis, there is a problem in how to choose . A small is needed for a reliable calculation of the transmission coefficient, but choosing a small to determine the surface Green’s function may not succeed because *G* can become long range and the iterative cycle used to generate it may not be stable. To accommodate these conflicting requirements, a surface-specific should be used. It is entered as an element **delend** in **EMESH**.

The mode=0 mesh is specified as

When computing transmission coefficients via the Landauer-Buttiker formalism, one chooses a contour as in mode=0. However, there is a problem in how to choose Imz{\mathrm{Im}}\, zImz. A small Imz{\mathrm{Im}}\, zImz is needed for a reliable calculation of the transmission coefficient, but using a small Imz{\mathrm{Im}}\, zImz to determine the surface Green’s function may not succeed because the GF can become long range and the iterative cycle used to generate it may not be stable. To accommodate these conflicting requirements, a surface-specific Imz{\mathrm{Im}}\, zImz should be used, called **delend**. The mode=0 mesh is specified as

### Other Resources

An overview of the Atomic Spheres Approximation can be found here.

If you haven’t already done so, you are advised to go through the tutorial for the band code **lm**) before doing this one. **lm** and **lmgf** share many feature in common, but **lm** is somewhat easier to use.

See this page for an Introductory tutorial for **lmgf**, and this page for a tutorial explaining how **lmgf** implements the CPA in practice.

There is a related program **lmpg**, that has a similar function but is designed for layer geometries, enabling Landauer-Buttiker transport.

### References

- I. Turek et al.,
*Electronic strucure of disordered alloys, surfaces and interfaces*(Kluwer, Boston, 1996). - J. Kudrnovsky and V. Drchal, Phys. Rev. B 41, 7515 (1990).
- J. Kudrnovsky, V. Drchal, and J. Masek, Phys. Rev. B 35, 2487 (1987).