<HTML><HEAD><TITLE>ASA Green's function package (v6.14) </TITLE><STYLE><!--P{font-size : medium;} --></STYLE></HEAD><BODY bgcolor="#cccccc">
<H1><B><FONT size="+2">ASA Green's function package (v6.14)</FONT></B></H1>

This package implements the ASA local spin-density approximation using
Green's functions.  It adds a program <B>lmgf</B> to the ASA suite, which
plays approximately the same role as the LMTO-ASA band program <B>lm</B>.

<H2><A name="section1"><FONT size="+1">
Structure of Green's function program
</FONT></A></H2>

<B>lmgf</B> runs in much the same way as program <B>lm</B>.

The the band pass routine of <B>lm</B>, bndasa, generates eigenvalues
and eigenvectors, which can in turn generate quantities of interest.
<B>bndasa</B> is replaced by a Green's function routine, <B>gfasa</B>.
<B>gfasa</B> can generate output moments DOS, density-matrix, etc., in
the same way as <B>bndasa</B>.

<P>

<B>gfasa</B> loops over some energy contour; see <A href="gf.html#section3">below</A> how to specify an energy contour
and <A href="../gf/emesh.f">../gf/emesh.f</A> for the routine that generates it.
<B>gfasa</B> then accumulates various kinds of data, such
as the energy moments used in an ASA self-consistent cycle; see below
for what kinds of things you can create right now.  Two particularly
useful features in <B>lmgf</B> are the efficient computation of the
screened exchange potential (which can also be done by <B>lm</B>, but
far less efficiently), and the computation of magnetic exchange interactions.

</P><P>

For each energy point, the BZ integration is accomplished by routine
in <A href="../gf/gfibz.f">../gf/gfibz.f</A>, which loops over all
irreducible points, generating the "scattering path operator" <I>g</I>
and the corresponding <I>g</I> for all the points in the star of <I>k</I> to
generate a properly symmetrized <I>g</I>.  Within the ASA,
second-generation LMTO, <I>g</I> is converted to proper Green's
function <I>G</I>, corresponding to the orthogonal gamma
representation by an energy scaling.  The scaling is carried out in
routine <A href="../gf/gfg2g.f">../gf/gfg2g.f</A>.  (<I>g</I> and the
converstion to <I>G</I> will change with the construction of NMTO's).

Next the various integrated quantities sought are assembled (done
by <A href="../gf/gfidos.f">../gf/gfidos.f</A>).  Also the program may
may be extrapolated to other points using a Pade approximation. 

Finally an estimate for the Fermi level is determined (see next
section), and if the original estimate is sufficiently close, the
cycle is finished as in program <B>lm</B>.  If the estimate is too far
off, a new energy mesh is taken and the process repeated.

<H2><A name="section2"><FONT size="+1">
Energy Contours, Potential Shifts and the Determination of the Fermi Level
</FONT></A></H2>

For self-consistency, a very fine energy mesh would be required if the
energy integration were carried out close to the real axis.  Instead,
the integration contour is deformed into an elliptical path in the
complex plane.

A Green's function method requires integration to the Fermi level,
which is determined by charge neutrality.  In these codes, the Fermi
level is specified by the user and the crystal is shifted by a
constant potential to obtain neutrality.  The Fermi level and constant
potential are stored in array <I>vshft</I>, and permanently on disk in
file <B>vshft</B> described below; routine <A
href="../subs/iovshf.f">../subs/iovshf.f</A>) does the I/O.  (Actually the
freedom to shift potentials is more general and can accomodate
potential shifts at separate sites, useful in non-self-consistent or
limited self-consistent calculations.)  Inspection of this array
<I>vshft</I> in the <A href="../gf/gfasa.f">source code</A>) shows
some apparently unecessary complications, which arise because the
crystal GF codes share routines in common with layer GF codes, which
have more complicated treatment of the electrostatics.

<P></P>

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.

<ul>
 <B>Metal case</B> (set by BZ METAL=1): once the <I>k</I>- 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.  It does this in one of two ways:

<P></P><ul>

 - Pade approximant to interpolate the diagonal elements of
 the Green's function.  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 charge-neutrality condition is
 satisfied.  In this mode, when the potential shift exceeds a certain
 tolerance, <B>gfasa</B> will repeat the cycle over energy and
 <I>k</I>-points and energy with the Pade-estimate for the potential
 shift.  This tolerance can be specified by the user using option
 padtol in token GFOPTS= below.  The user is advised to monitor these
 shifts and the deviation from charge neutrality.

<P></P>

 - Integration along 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.

<P></P></ul>

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

<P></P>

<BR>
 <B>Nonmetal case</B> (set by BZ METAL=0): <B>lmgf</B> will not attempt to
 shift the potential, or ensure charge neutrality.  The user is
 cautioned to 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.

</ul>

<H2><A name="section3"><FONT size="+1">
GF-specific input
</FONT></A></H2>

<H2><A name="section3.0"><FONT size="+0">
Energy integration contour
</FONT></A></H2>

Green's functions require an energy integration contour.
the energy contour is specified in category BZ as

<pre>
  EMESH= nz mode emin emax [other args, depending mode]
</pre>
where
<pre>
         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)
</pre>

<br>
  Right now there are the following contours:

<P>
<ul>
  <B>mode=0</B>: a uniform mesh of points 
          between emin and emax, with a constant imaginary component.
<pre>
  EMESH= nz 0 emin emax Im-z [<A href="#section3.3">... + possible args for layer geometry</A>.]

         Im-z is the (constant) imaginary component.
</pre>
  This mode is generally not recommended for self-consistent cycles
  because the GF has a lot of structure close to the real axis (small
  Im-z), while shifting off the real axis introduces errors.  It is
  used, however, in other contexts, e.g. <A href="#section3.3">transport</A>.

</ul>

<H2><A name="section3.1"><FONT size="+0"></FONT></A></H2>

<ul>
  <B>mode=10</B>: a Gaussian quadrature on an ellipse. 

<pre>
  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 Ef.
</pre>
<BR>
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 <A href="gf.html#section2">section 2</A>.  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.
</ul>

<ul>
  <B>mode=310</B>: 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.

<pre>
  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.
</pre>

<BR><BR>
  <B>mode=2</B>: is the same contour as <A href="#section3.0">mode=0</A>.  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.

</ul>

<H2><A name="section3.2"><FONT size="+0"></FONT></A></H2>

<ul>
  <B>mode=110</B>: 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 <B>lmpg</B>.)

Thus the integration proceeds in two parts: first an integration on an
elliptical path is taken to the left Fermi level (as in <A href="#section3.1">mode=10</A>).
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 <A href="#section3.0">mode=0</A>.  For the nonequilibrium contour, three additional pieces of
information must be supplied:
<pre>
  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 Im-z on the nonequilibrium contour
</pre>
<br>
The mesh is specified as 
<pre>
  EMESH= nz 110 emin ef(L) ecc eps nzne vne delne [del00]
</pre>

The last argument plays the role of delne specifically for computing
the self-energy that determines the end boundary conditions.  There is
an incompatibility in the requirements for Im-z in the central and end
regions. the same incompatibility applies to transport and is
discussed [<A href="#section3.3">below]</A>.

<H2><A name="section3.3"><FONT size="+0">Modifications of energy contour for layer geometry</FONT></A></H2>

<ul>
When computing transmission coefficients via the Landauer-Buttiker
formalism, one chooses a contour as in <A
href="#section3.1">mode=0</A>.  However, a there is a problem in how
to choose <tty>Im-z</tty>.  A small Im-z is needed for a
reliable calculation of the transmission coefficient, but using a
small 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 accomodate these conflicting requirements, a surface-specific
<tty>Im-z</tty> should be used, called <tty>del00</tty>.
The <tty>mode=0</tty> mesh is specified as 
<pre>
  EMESH= nz 0 emin emax delta xx xx xx xx del00
</pre>
<tty>delta</tty> is <tty>Im-z</tty> for the central region;
<tty>del00</tty> is <tty>Im-z</tty> for the surfaces.<br>
Entries <tty>xx</tty> have no meaning but are put there for
compatibility with the contour used in nonequilibrium
calculations.  (A similar situation applies to the <A
href="#section3.1">nonequilibrium part of the contour</A>).


<br><br>
The mesh for <A href="#section3.1">self-consistent nonequilibrium calculations</A> is
<pre>
  EMESH= nz 110 emin ef(L) ecc eps nzne vne delne del00
</pre>

</ul>
</ul>

<H2><A name="section3.4"><FONT size="+0">
Green's function category
</FONT></A></H2>

<B>lmgf</B> requires a GF-specific category.

<pre>
GF  MODE=1 GFOPTS=options
</pre>

<P>
Token <B><tt>MODE=</tt></B> controls what <B>lmgf</B> calculates.  Options are 
<tt>MODE=1</tt>, <tt>MODE=10</tt>, <tt>MODE=11</tt>, described below.


<ul>
  <B><tt>MODE=1</tt></B> goes through the usual cycle calling gfasa.  It performs a
  function analogous to bndasa in the band program, generating output
  density, moments, and optionally other quantities such as
  density-of-states.
<BR>
  Taken with the special integration contour mode=2 (see <tt>EMESH</tt> above), N(E)
  and its integral are computed and tabulated over the window
  specified.  With the following sample input segment:

<BR>
<pre>
  % const ef=-0.025725
  BZ      EMESH=5 2 {ef} {ef+.002*4} .001 0

  the integration would be tabulated for five points
  ef, ef+.002, ef+.004, ef+.006, ef+.008 like so (spin-polarized 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
</pre>


  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.
<P>
  <B><tt>MODE=10</tt></B> invokes a special branch that computes magnetic exchange
  interactions using a linear response technique.  In particular, J_ij is
  computed for pairs of sites (i,j), where the J's are the parameters in
  the Heisenberg hamiltonian
<pre>
    E(s_i, s_j) = sum_ij J_ij s_i . s_j
</pre>

  Thus, the J's are coefficients to energy changes for small rotations
  of the spins. It can be computed from a change in the band energy,
  and for small rotations can be done analytically.
<P>
  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.
<P>
  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 is 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 switch
<pre>
   --sites:pair:site-list
</pre>
  See <A href="#section4">command-line arguments</A></B> for more details.
<P>
  <B><tt>MODE=11</tt></B> is an exchange branch that is run after <tt>MODE=10</tt>.  It prints out
  the J_ij and does several other analyses.  Switch
<pre>
   --sites:pair:site-list
</pre>
  also applies to this mode; see see <A href="#section4">command-line arguments</A>.

</ul>

Token <B><tt>GFOPTS= option-list</tt></B> causes gfasa to do a variety of different things.

<pre>
  Options are entered as a series of strings string1;string2;...  
  Option strings are:

    emom   generate the output ASA moments, needed for self-consistency
    idos   make integrated properties, such as the sum of one-electron energies
    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 a file 'dmat'
    pdos   Make the partial density of states (this has never been checked)
    p3     Use third order potential functions
    padtol Set the tolerance for maximum potential shift permissible by
           Pade interpolation, as described above
    frzvc  Suppress saving the constant potential shift used 
           to determine charge neutrality
</pre>

<H2><A name="section4"><FONT size="+1">
lmgf-specific command-line arguments
</FONT></A></H2>

<pre>

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

</pre>
  The following are specific to the exchange calculation modes 10 and 11:
<pre>
  --sites[:pair]:site-list  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.
</pre>
  See <A href="lmto.html#IntegerListSyntax">Syntax of Integer
  Lists</A> for the syntax of `site-list'.
<BR>
  Example:  running lmgf using <tt>MODE=10</tt> and this command line argument
<pre>
          --sites:pair:1,3,5,7
</pre>
  generates J connecting sites 1, 3, 5 and 7 to all neighbors.
<BR><BR>
  Running lmgf using <tt>MODE=11</tt> with the same <tt>--sites</tt> switch
  argument will print out the exchanges just between pairs of these sites.
<BR><BR>
  Running lmgf using <tt>MODE=11</tt> without any <tt>--sites</tt>
  argument will print out the exchanges between these sites and all neighbors.
<BR><BR>

<pre>
  --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

</pre>

<H2><A name="section5"><FONT size="+1">
Installation, test cases, and examples
</FONT></A></H2>

See <B><A href="../gf/README">../gf/README</A></B>


