The ASA Green's function program lmgf
Table of Contents
 Purpose
 Structure of Green’s function program
 Energy Contours, Potential Shifts and the Determination of the Fermi Level
 GF specific input
 lmgfspecific commandline arguments
 Test cases and examples
Purpose
This package implements the ASA local spindensity approximation using Green’s functions. The Green’s functions are contructed by approximating KKR multiplescattering theory with an analytic potential function. 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. With this package a new program, lmgf is added to the suite of executables. lmgf plays approximately the same role as the LMTOASA band program lm: you can use lmgf to make a selfconsistent density as you can do with lm. A potential is generated from energy moments $Q_0$, $Q_1$, and $Q_2$, 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. The most imprortant ones are:
 Calculate magnetic exchange interactions
 Calculate magnetic susceptibility (spinspin, spinorbit, orbitorbit parts)
 Calculate properties of disordered materials, either chemically disordered or spin disorder from finite temperature, within the Coherent Potential Approximation, or CPA.
 Calculate the ASA static susceptibility at $q=0$ to help converge calculations to selfconsistency.
Structure of Green’s function program
lmgf runs in much the same way as lm. 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, densitymatrix, etc., in the same way as bndasa does.
In contrast to band methods (implemented in lm) where the Hamiltonian $H$ 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 $G$ for a particular energy.
This fact highlights the strengths and weaknesses of a Green’s function approach. Energyintegrated properties such as the moments, must be obtained by integrating over energy. Calculating $G$ 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 energydependence is needed anyway. CPA theory yields an energydependent 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 selfconsistent cycle. Finally, an estimate for the Fermi level $E_F$ is determined using a from Pade approximation. If the original guess for $E_F$ 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.
Energy Contours, Potential Shifts and the Determination of the Fermi Level
For energyintegrated properties, 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_{F} is not known but must be fixed by charge neutrality. Thus $E_F$ must be guessed at and iteratively refined until the charge neutrality condition is satisfied. lmgf does not vary E_{F}; the user specifies it at the outset. Instead lmgf looks for a constant potential shift that satisfies charge neutrality; this must be searched for iteratively. Both the potential shift and $E_F$ are maintained in a file vshft.ext. Inspection of vshft.ext may look unecessarily complicated; it’s because you can use the file to add sitedependent shifts. 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.
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 $k$ and energypoints are summed over and the deviation from charge neutrality is determined, the code will attempt to find the potential shift that fixes charge neutrality.
It finds the Fermi level in one of two ways:

Using a Pade approximant, lmgf interpolates the diagonal elements of G. The interpolation is used to evaluate the GF on the starting elliptical contour shifted rigidly by a constant, and the shift is iterated until the chargeneutrality condition is satisfied. At this stage, there are two possibilities:
1. repeat the integration of $G$ over $k$ and the energy contour with the constant shift added to the potential.
2. Assume that the Padeapproximant to the diagonal $G$ is a sufficiently good estimate for the actual $G$.
If the potential shift is larger than a userspecifed tolerance (see padtol in GF_GFOPTS below), option 1 is taken and the Pade shift reevaluated. A new Pade estimate is made for the potential shift requiring charge neutrality, and it is tested once against the userspecified tolerance.
When the shift falls below the tolerance, option 2 is taken and lmgf proceeds to the next step. The user is advised to monitor these shifts and the deviation from charge neutrality.

The charge is integrated in a contour near the real axis subsequent to the elliptical contour. In this mode, the determination of the potential shift is accomplished by continuing the integration contour on the real axis starting from the originally estimated Fermi level. A trapezoidal rule is used (or Simpson’s rule using a Pade approximate for the midpoint), and new energy points are computed and integrals 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 the METAL case: by default the program will save the potential shift to use in the next iteration. You can suppress this save (see frzvc below), 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).
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.
Some details concerning how lmgf works internally
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” $g$ and the corresponding $g$ for all the points in the star of $k$ to generate a properly symmetrized $g$. Within the ASA, secondgeneration LMTO, $g$ is converted to proper Green’s function $G$, corresponding to the orthogonal gamma representation 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). The potential shift to satisfy charge neutrality is found, and stored in vshft.ext. The I/O is handled by routine subs/iovshf.f.
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.) $G$ is “spikey” for energies on the real axis (it has poles where there are eigenstates).
The BZMESH token
To compute energyintegrated 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. densityofstates or spectral functions. You specify the contour in category BZ as:
EMESH= nz mode emin emax [other args, depending mode]
where
nz number of energy points
mode specifies the kind of contour; see below
emin,emax are the energy window (emax is usually the Fermi level)
Right now there are the following contours:
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, e2=0 is good, or maybe e2=.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, the mesh is rigidly shifted by a constant amount, and the diagonal GF interpolated using a Pade approximant to the shifted mesh. The shifting+interpolation is iterated until charge neutrality is found, as described in section 2. If the rigid shift exceeds a specified tolerance, the Pade interpolation may be suspect. Thus, the entire cycle is repeated from scratch, on the shifted mesh where the shift is estimated by Pade.
mode=0: a uniform mesh of points between emin and emax, with a constant imaginary component.
EMESH= nz 0 emin emax Imz [... + possible args for layer geometry.]
Imz is the (constant) imaginary component.
This mode is generally not recommended for selfconsistent cycles because the GF has a lot of structure close to the real axis (small ${\mathrm{Im}}\, z$), while shifting off the real axis introduces errors. It is used, however, in other contexts, e.g. transport.
mode=110: is a contour input specific to nonequilibrium Green’s function.
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, ef(R)ef(L)
delne Imz on the nonequilibrium contour
The mesh is specified as
EMESH= nz 110 emin ef(L) ecc eps nzne vne delne [del00]
The last argument plays the role of delne specifically for computing the selfenergy of the left and right leads. There is an incompatibility in the requirements for $\mathrm{Im}\,z$ in the central and end regions. The same incompatibility applies to transport and is discussed in the following section.
mode=310: Alternative Pade approximant in finding Fermi level.
a Gaussian quadrature on an ellipse to a trial emax, as in mode 2. 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.
mode=2: is the same contour as mode=0. However, it is designed for cases when you want to resolve the energy dependence of some quantity, such as the DOS or magnetic exchange coupling. These are discussed in the GF category below.
Modifications of energy contour for layer geometry
When computing transmission coefficients via the LandauerButtiker formalism, one chooses a contour as in mode=0. However, there is a problem in how to choose ${\mathrm{Im}}\, z$. A small ${\mathrm{Im}}\, z$ is needed for a reliable calculation of the transmission coefficient, but using a small ${\mathrm{Im}}\, z$ 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 surfacespecific ${\mathrm{Im}}\, z$ should be used, called del00. The mode=0 mesh is specified as
EMESH= nz 0 emin emax delta xx xx xx xx del00
delta is Im z for the central region; del00 is ${\mathrm{Im}}\, z$ for the surfaces.
Entries xx have no meaning but are put there for compatibility with the contour used in nonequilibrium calculations. (A similar situation applies to the nonequilibrium part of the contour).
The mesh for selfconsistent nonequilibrium calculations is
EMESH= nz 110 emin ef(L) ecc eps nzne vne delne del00
Green’s function category
lmgf requires a GFspecific 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 selfconsistency cycle, calling gfasa. It performs a function analogous to bndasa in the band program, generating output density, moments, and other quantities such as densityofstates.
Taken with the special integration contour mode=2 (see EMESH above), the densityofstates $D(E)$ and its integral are computed and tabulated over the window specified.
With the following sample input segment:
% const ef=0.025725
BZ EMESH=5 2 {ef} {ef+.002*4} .001 0
The integration will be tabulated for five points ef, ef+.002, ef+.004, ef+.006, ef+.008 like so (spinpolarized case)
Re z Im z spin dos idos
0.025725 0.001000 1 13.55272 0.00000
0.025725 0.001000 2 10.38435 0.00000
0.025725 0.001000 t 23.93706 0.00000
0.023725 0.001000 1 9.17407 0.02273
0.023725 0.001000 2 4.13694 0.01452
0.023725 0.001000 t 13.31101 0.03725
0.021725 0.001000 1 15.33776 0.04724
0.021725 0.001000 2 7.42200 0.02608
0.021725 0.001000 t 22.75976 0.07332
0.019725 0.001000 1 19.58433 0.08216
0.019725 0.001000 2 7.52708 0.04103
0.019725 0.001000 t 27.11141 0.12319
0.017725 0.001000 1 20.83078 0.12258
0.017725 0.001000 2 9.31350 0.05787
0.017725 0.001000 t 30.14428 0.18045
If the partial DOS is generated, the usual tokens in the BZ category specifying the window (DOS=) and number of points (NPTS=) are overridden by the parameters in EMESH.
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, $J_{ij}$ is computed for pairs of sites $(i,j)$, where the J’s are the parameters in the Heisenberg hamiltonian
$E(s_i, s_j) = \sum_{ij} J_{ij} s_i s_j$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 energydependence 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 commandline argument
sites:pair:sitelist
For more details, see commandline 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 $J_{ij}$ and does several other analyses.
Switch sites:pair:sitelist
applies to mode 11 as well as mode 10; see commandline 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 selfconsistency 
noemom  suppress generation of the output ASA moments 
idos  make integrated properties, such as the sum of oneelectron energies 
noidos  reverse of idos 
dmat  make the densitymatrix G_{RL,R’L’} 
sdmat  make the sitediagonal densitymatrix G_{RL,RL’}. The density matrix is written to dmat.ext. 
pdos  Make the partial density of states (this has not been checked recently). 
p3  Use third order potential functions 
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 
padtol  Set the tolerance for maximum potential shift permissible by Pade interpolation, as described above 
The following are specific to the CPA:
omgtol  Tolerance in the Omega potential, CPA selfconsistency 
omgmix  How much of prior iterations to mix, CPA selfconsistency 
nitmax  Maximum number of iterations for CPA selfconsistency 
lotf  Learn onthefly 
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 
lmgfspecific commandline 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]:sitelist Make the exchange parameters J_ij only for sites i in the site list.
In mode 11, option :pair means that only parameters J_ij where both i and j are printed.
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. See Syntax of Integer Lists for the syntax of sitelist.
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 kmesh.
Can be helpful in testing kconvergence 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 fullpotential 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
Edit This Page