<HTML><HEAD><TITLE>Full-potential program lmf (v6.17) </TITLE><STYLE><!--P{font-size : medium;} --></STYLE></HEAD><BODY bgcolor="#cccccc">
<H1><B><FONT size="+2">Full-potential program lmf (v6.17)</FONT></B></H1>

<H2><A name="section1"><FONT size="+1">
1.  Overview of the full-potential method
</FONT></A></H2>

<P>

The full-potential program <B>lmf</B> was originally adapted from a program
<B>nfp</B> written by M. Methfessel and M van Schilfgaarde.  The
method is described in some detail in the following reference (click
<B><A href="nfp-doc.ps.gz">here</A></B> to view postscript file): <UL>
M. Methfessel, Mark van Schilfgaarde, and R. A. Casali, ``A
full-potential LMTO method based on smooth Hankel functions,'' in
  <I>Electronic Structure and Physical Properties of Solids: The Uses
  of the LMTO Method</I>, Lecture Notes in Physics
  <B>535</B>. H. Dreysse, ed. (Springer-Verlag, Berlin) 2000.<BR>
</UL>

The following features are unique to this method:
<BR>

<DL>
<DT><STRONG>Smoothed Hankel functions</STRONG>
<DD>
  The envelope functions are smoothed generalizations of Hankel functions
  that found in LMTO programs.  Unlike the normal Hankel functions, the
  smoothed versions are regular at the origin.  They are a significantly
  better choice of basis than the customary LMTO basis set.  However, the
  smoothing introduces complications because the augmentation of a smoothed
  Hankel function is less straightforward than of a normal Hankel.  The
  envelope functions are not screened into a tight-binding representation,
  as in the second-generation and later generation LMTO methods; thus wave
  functions are evaluated by Ewald summation.  Regrettably, this makes the
  method rather slow.  A screened version is in progress.

<DT><STRONG>Local Orbitals</STRONG>
<DD>
  This package extends the linear method through the use of local
  orbitals.  Augmented methods substitute radial solutions of the
  Schrodinger equation with combinations of partial waves of angular
  quantum number <i>l</i> inside the augmentation region.  Linear
  methods used a fixed radial function (more precisely, pair of
  functions), which has validity over only a certain energy window.
  With local orbitals, a third radial function is added to the basis,
  which greatly extends the energy window over which energy
  eigenvalues can be calculated.  It is necessary, for example, to
  obtain the reliable LDA bandgap in GaAs, as the Ga 3<i>d</i> and
  4<i>d</i> partial waves are both important.  To include local
  orbitals in the basis, see <B><A href="#localorbitals">here</A></B>
  for documentation.
  
<DT><STRONG>Augmented Plane Waves</STRONG>
<DD>

  In July 2008 Takao Kotani added APWs as additional envelope
  functions, which can increase the flexibility of the basis.  One can
  view this package as an extension of a conventional LAPW method,
  enabling through the use of a few MTOs a much faster convergence in
  energy cutoff of APWs.  Alternatively, it can be viewed as an
  extension of the original MTO method, With APWs, a good quality
  basis is easier to construct than it is with a strictly with the MTO
  basis.  Also APWs offer a systematic way of converging the accuracy
  of the combined MTO+APW basis systematic and reliable manner, to
  almost arbitarily high accuracy.  This is particularly important
  when reliable eigenvalues far above the Fermi level are needed, and
  to check the accuracy of a given MTO basis.   To include APWs
  in the basis, see <B><A href="#apw">here</A></B>
  for documentation.

<DT><STRONG>Augmentation and Representation of the charge density</STRONG>
<DD>

 The charge density representation is unique to this method.  It
 consists of three parts: a smooth density <i>n</i><sub>0</sub>
 carried on a uniform mesh, defined everywhere in space (<i>n</i><sub>0</sub> is not
 augmented, as occurs with conventional augmentation), the true
 density <i>n</i><sub>1</sub> expressed in terms of spherical harmonics
 <I>Y</I><sub>lm</sub> inside each augmentation sphere, and a
 one-center expansion of the smooth density <i>n</i><sub>2</sub>, also
 expanded in terms of harmonics <I>Y</I><sub>lm</sub> inside each
 augmentation sphere.  The total density is expressed as in the
 "threefold represesentation"
   <i>n</i> = <i>n</i><sub>0</sub> + <i>n</i><sub>1</sub> - <i>n</i><sub>2</sub>.

 This turns out to be an extremely useful way to carry out the
 augmentation procedure.  Quantities converge much more rapidly in the
 <i>l</i>-truncation in the augmentation sphere than occurs in
 conventional augmented wave constructions.  (The analysis is a little
 subtle; see the reference at the start of this document.)

</DL>


<P>The full-potential program builds on the basic package (documented
   in file <B><A href="lmto.html">lmto.html</A></B>) which contains an
   implementation of a tight-binding LMTO program in the Atomic
   Spheres Approximation (ASA), and shares most things in common with
   it, including a number of auxiliary programs useful to both ASA and
   this one.  For example, both methods are linear augmented-wave
   methods, and the wave functions inside the augmentation spheres are
   equivalent in the two cases.  The reader is strongly advised to
   review that documentation before proceeding with this one.  There
   are some differences in input between the two programs.  Tokens not
   used by the fp programs are marked as ASA-specific in file <B><A
   href="lmto.html">lmto.html</A></B>; others that are specific to
   this method are documented in <A href="#section3">Section 3</A>.
   See <A href="FPtutorial.html">FPtutorial.html</A> for a tutorial
   that helps you build an input file, and explains the output, and
   <A href="ASAtutorial.html">ASAtutorial.html</A> for a corresponding
   tutorial for the ASA programs.

<P>
 One important difference between is that the FP method has no neat
 parameterization of total density in terms of <A
 href="lmto.html#section3">the ASA energy moments Q0,Q1,Q2</A>, or the
 representation of the potential by a few potential parameters, as in
 the ASA (see ASA overview in <B><A
 href="lmto.html">lmto.html</A></B>).  However, the basis within the
 augmentation spheres is defined from the spherical average of the
 potential, just as in the ASA, and the linearization proceeds in the
 same way.  Program <B>lmf</B> uses the same ``<A
 href="lmto.html#section3">continuously variable principal quantum
 numbers</A>'' P described to establish a mapping between the
 linearization energy and logarithmic derivative at the MT boundary,
 and to float the linearization energy to band center-of-gravity.
 Thus, the description of P and IDMOD equally apply to the ASA and
 this program.

<P>
A second important difference is that the basis set is more
complicated, and in its current form, the user must choose parameters
defining the basis.  This complication is the most onerous part of the
present method, and there are plans for redesign; but at present is
necessary to treat the interstital region reliably.  The inputs are
described below; see also <A href="FPoptbas.html">FPoptbas.html</A>
for a tutorial about choosing the basis optimally.  Note also with
the incorporation of the <B><A href="#apw">APW basis</A></B>,
other options are available.

</P>
<H2><A name="programs"></A></H2>
<P> The FP package adds two executable programs to the basic one:
<UL><LI>
 Program <B>lmf</B> is the program used for self-consistent full-potential
 calculations.  It requires a starting density, which it obtains either
 from a binary restart file `rst' (typically generated by a prior
 invocation of <B>lmf</B>) or a formatted (almost) equivalent of it `rsta', or
 from a superposition of free-atom densities (file `atm', which is created
 as follows).
 <BR> <BR>
  <LI>
 Program <B>lmfa</B> makes each species self-consistent for the free atom and writes
 the local density, plus a fit to the tail beyond the MT radius (fit as a
 linear combination of sm-Hankel functions) to file `atm'.  Also appended to
 file `atm' is a fit to the free-atom wave functions for r>rmt.  In version
 6.9n, this is where the envelope wave functions are specified, and as the
 functions are an excellent choice for a free atom, they are not tailored to
 the solid, and you are encouraged to tinker with the basis (described in
 the next Section) before proceeding to carry out full calculations.
</UL>

<HR>
<H2><A name="section2"><FONT size="+1">
2.  Capabilities of the package
</FONT></A></H2>
<P>

<B>lmf</B> in an implementation of the local-density approximation,
and generates density-of states, and total energy as a part of the
self-consistency cycle.

<P> Other extensions are:

<P> <FONT color="#bb3300"><I>Energy bands</I>.</FONT>
You can run <B>lmf</B> in band mode to generate energy
bands along lines or planes for, generating, e.g. Fermi surfaces.
<I>k</I>-point specifications and invocation is identical to that
described in the basic <B><A href="lmto.html#section9lm">ASA
suite</A></B>.  Particularly useful are the <A
href="ASAtutorial.html#bandweights">color weights.</A>

<P> <FONT color="#bb3300"><I>Partial DOS and Mulliken analysis</I>.</FONT>
<B>lmf</B> can generate partial dos within augmentation spheres, and
construct a Mulliken analysis.  DOS can be resolved by site, by site and
<I>l</I>, or by site and <I>lm</I>.  These options are invoked through
<B><A href="fp.html#section4">command-line switches</A></B>.  For
illustrations, invoke
<pre>
   fp/test/test.fp co 2
   fp/test/test.fp fe 2
</pre>

<P> <FONT color="#bb3300"><I>Charge density</I>.</FONT>

<B>lmf</B> can generate the charge density (smooth or not, with or
without cores), and the contribution to the density from a selected
window of states.  See <tt>--wden</tt> and <tt>--window</tt> in
<B><A href="fp.html#section4">command-line switches</A></B>

<P> <FONT color="#bb3300"><I>Core-level spectroscopy</I>.</FONT>
<B>lmf</B> can generate EELS spectra, which involve matrix elements
between core and valence electrons.  The EELS option is invoked with a <B><A href="fp.html#section4">command-line switches</A></B>
For illustrations, invoke
<pre>
   fp/test/test.fp fe 2
   fp/test/test.fp crn 2
</pre>

<P> <FONT color="#bb3300"><I>LDA+U</I>.</FONT>

 The LDA+U functional was built into <B>lmf</B> in v6.15 and later, by
 Walter Lambrecht.  LDA+U needs in addition to the LDA parameters
 <i>U</i> and <i>J</i> for selected orbitals, which are empirical.
 The LDA+U constructs an additional potential for a particular
 <I>l</I> subblock (<I>m</I>=<I>-l</I>..<I>l</I>) from the <I>U</I>
 supplied by the user, and the density-matrix, which is generated by
 <B>lmf</B>.  

 Two modifications must be added to the input, which are described 
 <A href="#ldaplusu">here</A>.

 In a strictly LDA calculation, complete information is carried by the
 density, contained in the restart file, <tt>rst.ext</tt>. In the
 LDA+U case, complete information is carried by  density and
 on-site density matrices, which are contained in file <tt>dmats.ext</tt>.

 An example that illustrates LDA+U method is ErAs, which you can run by
<pre>
   fp/test/test.fp eras
</pre>

ErAs is an interesting case because LDA puts all 4 minority <I>f</I>
electrons in a single extremely narrow band at the Fermi level.  In
LDA+U the minority <I>f</I> is spit into a 4- and 3-manifold; see PRB
67, 035104 (2003).

<P> <FONT color="#bb3300"><I>GW</I>.</FONT>
<B>lmf</B> is designed to work in coordination with a GW package by
T. Kotani (the GW package comes comes separately).  <B>lmf</B> acts
both as a driver for the GW package and also can be used in a
self-consistent GW cycle.  An extra driver <B>lmfgw</B> is compiled
as part of this extensions.  Use of this driver is described in the
<B><A href="gw.html">GW driver documentation</A></B>.  You need the
extension package GW.version.tar.gz.  Also you will need the GW package itself.
For illustrations of the driver invoke
<pre>
   gw/test/test.gw si
   gw/test/test.gw gas
</pre>

<P> <FONT color="#bb3300"><I>Spin-Orbit coupling</I>.</FONT>

Starting with v6.15, <B>lmf</B> can include spin-orbit coupling
perturbatively to the hamiltonian (courtesy of A. Chantis).  It is
possible to include the full L.S, or just the Lz.Sz part.  The
coupling works in conjuction with the self-energy generated by the GW method.
For illustrations of both L.S and Lz.Sz invoke
<pre>
   fp/test/test.fp felz
   fp/test/test.fp gasls
</pre>


<P> <FONT color="#bb3300"><I>Local orbitals</I>.</FONT> 

 In v6.12 and later, <B><A href="#localorbitals">local
 orbitals</A></B> may be added to the basis set. 
 These orbitals are important when energy bands over a very wide
 energy window are required, when high accuracy is needed for shallow
 (semi)core states, or for energies far above the Fermi energy.
 Examples of the former occur in oxides: bond lengths are small and
 cations with shallow <I>p</I> orbitals extend somewhat beyond the
 augmentation radius. 

<P>

 Local orbitals are presented in one of two flavors.  The first,
 conventional type of local orbital is constructed by solving the
 radial Schrodinger at a different linearization energy than the usual
 valence states, and then subtracting off a particular linear of the
 valence wave function phi and energy derivative phidot such that the
 local orbital's value and slope vanish at the augmentation
 radius. Thus

 <BR> <UL>

   <LI> A local orbital is strictly confined to the augmentation
         sphere, and has no envelope function at all.

   <LI> When taken in linear combination with the valence augmentation
         functions, it can solve the Schrodinger equation exactly for
         linearization the energy chosen (that is, for the spherical
         potential that defines the wave functions).

   <LI> It turns the <I>linear</I> method into a <I>quadratic</I> one.
        Of course, it would be possible, and more accurate, to
        construct a LMTO orbital of a different principal quantum
        number complete with envelope function; howewer, there is a
        corresponding loss in efficiency because the additional matrix
        elements of the envelope function must be evaluated.

   <LI> for states with energies in the vicinity of the local orbital
        energy, the envelope functions of the regular valence states
        combine with the local orbital to make the interstitial part
        of the wave function. 

</UL>

Example that illustrate local orbitals of this type are
<pre>
   fp/test/test.fp gas
   fp/test/test.fp cu
</pre>
The Ga 3<I>d</I/> semicore the high-lying As 5<I>s</I/> state are
included as local orbitals.  In the Cu case, the high-lying Cu
4<I>d</I/> is included, which is important in GW calculations.

<P>

The second, extended, kind of local orbital can only be used for
semicore states.  Instead of artificially subtracting off some linear
of the phi and phidot to make the orbital vanish at the augmentation
radius, a smooth Hankel tail is attached to the orbital.  The
smoothing of tail is constructed to match as well as possible the
kinetic energy of the semicore state.  This type of orbital has the
advantage that the valence envelope function need not `carry' the tail
of the semicore state.  Its drawback is that more things can "go
wrong," namely it may fail to do a good job of fitting the kinetic
energy.

an example that illustrates the second kind of local orbital is
<pre>
  fp/test/test.fp srtio3 
</pre>
The Sr 4<I>p</I/> and Ti 3<I>p</I/> semicore states are included as
local orbitals.  In the first part of the test, they are included as
local orbitals of the first type; then the last step is recalculated
using local orbitals of the second type.

<P> <FONT color="#bb3300"><I>Floating Orbitals</I>.</FONT> 

 In v6.15 and later, <B><A href="#floatingorbitals">floating
 orbitals</A></B> may be added to the basis set.  These orbitals can
 be important when very accurate calculations are needed in open
 systems, e.g. when reliable energy bands are needed for a wide energy
 window.  These orbitals differ from the usual smooth Hankels in that
 they are not centered at an atom.  They are augmented just as the
 other orbitals, but there is no "head" sphere.

<P> The following illustrate the inclusion of floating orbitals in the basis:

<pre>
   fp/test/test.fp te
   fp/test/test.fp gaslc
</pre>

<P> <FONT color="#bb3300"><I>Augmented Plane Waves</I>.</FONT> 

 In v6.17 and later, <B><A href="#apw">Augmented Plane Waves</A></B>
 can also be included in the basis.  They play the same role as
 floating orbitals do, but APWs are superior because they are simpler
 to use (there are no parameters and they do not need to be located at
 any particular site), and the control over convergence is more
 systematic.

<P> The following illustrate the inclusion of APWs in the basis:

<pre>
   fp/test/test.fp te
   fp/test/test.fp srtio3
   fp/test/test.fp felz 4
</pre>

<P> <FONT color="#bb3300"><I>Parallel Implementation</I>.</FONT> 

Two separate parallel versions of <B>lmf</B> have been made (courtesy
of A. T. Paxton).  One parallelizes over <I>k</I>-points, and is the
most efficient for scaling; the other parallelizes many points at a
lower level.  Installation is not automatic, however.
See <B><A href="MPI-instructions-for-lmf">instructions</A></B> for installation.

<HR>
<H2><A name="section3"><FONT size="+1">
3.  Additional input for the full-potential programs
</FONT></A></H2>
<P>
The full-potential reads from the same input file as the basic package,
documented in <B><A href="lmto.html">lmto.html</A></B>. File <B><A
href="input-file-style.txt ">input-file-style.txt</A></B> describes in
generic terms how the input file is grouped into categories and tokens;
<B><A href="lmto.html">lmto.html</A></B> documents tokens shared by both the basic package and this one.
This section documents tokens specific to <B>lmf</B>.

<H2><A name="spec"></A></H2>
<UL>
<LI>Category SPEC (fp extensions; see also
    <A href="lmto.html#SPEC">SPEC</A> in the ASA documentation)

<BR><BR>

First, the <A href="lmto.html#SPEC"> muffin-tin radii</A> should be
taken somewhat smaller than in the ASA.  There is no volume-filling
condition as found in the ASA.  Strictly speaking, there should be no
sphere overlaps, though in practice overlaps up to about 10% seem to
engender minimal error, because of the way the augmented potential is
constructed.  See section on "<A href="lmto.html#section8">Selection of Sphere
Radii</A>" for a discussion of choosing sphere radii and a way to
automatically find good choices for them.

<BR> <BR>
  The most important, and also the most difficult new feature of this
  method has to do with specification of
  <FONT color="#bb3300"><I>basis functions</I></FONT>.
  An envelope function takes the form of a sm-Hankel, which is
  specified by a smoothing radius <FONT size="+1"><tt>RSMH</tt></FONT> and an energy <FONT size="+1"><tt>EH</tt></FONT>.  Each of
  these may be <I>l</I>-dependent.  <B>lmfa</B> finds the <I>optimum
  choice</I> for each <I>l</I> and appends the results to the `atm'
  file, albeit for the <I>free atom</I> case.)  For example, in the
  copt test case something like the following lines are appended to
  file `atm'

 <BR> <UL>
   BASIS:
 <BR> <UL>
     Co RSMH= 2.425 2.717 1.017 EH= -0.360 -0.200 -0.222
 <BR>
     Pt RSMH= 2.461 3.042 1.085 EH= -0.441 -0.200 -0.200
  </UL>  </UL> <BR>

 In this case the basis on the Co site consists of an <I>s</I> orbital with
 (rsmh,eh)=(2.425,-0.36); a <I>p</I> orbital with (rsmh,eh)=(2.717,-0.2);
 and a <I>d</I> orbital with (rsmh,eh)=(1.0717,-0.222).  <B>lmf</B> requires a
 basis specification of that style in the input file, for each species in
 the SPEC (or CLASS) category.  Often you can just "cut and paste" the
 lines from the bottom of the atm file, and pluck them straightaway in the
 input file.  But of course what is optimum for the free atom may be very
 different from the solid case, and the user is strongly
 recommended to optimize the basis set before doing any systematic
 calculations, as described below.

 You may specify a single sm-H per <I>l</I>-channel or two sm-H.  Here is an
 a segment of an input file with a two sm-H declaration:
 <BR> <UL>
   SPEC   ATOM=A
 <BR> <UL>
           RSMH= 2,2,2       EH=-.1,-.1,-.1
 <BR>
          RSMH2=1.5,1.5,1.5 EH2=-1,-1,-1
  </UL>  </UL> <BR>


 You are strongly encouraged not to just take what <B>lmfa</B> supplies, but to spend a
 little time tinkering with these parameters to minimize the energy.  A quick
 way to do that is to use the HF=t option, which always used overlapped
 free-atomic densities and makes only one band pass.  Some rules of thumb:

<OL>
<LI> The total energy is more sensitive to the choice of smoothing radii than
     the energies; usually there is little benefit in optimizing the
     energies.   The main exception to this rule occurs for deep orbitals,
     such as the <I>s</I> orbital on a semiconductor anion.
     <BR>

<LI> For <I>sp</I> orbitals, a good rule of thumb is to use <FONT size="+1"><tt>RSMH</tt></FONT>=2/3
     muffin-tin radius.  Sometimes (particularly in cases with <I>d</I>
     electrons) a somewhat better choice is <FONT size="+1"><tt>RSMH</tt></FONT> slightly greater than RMT.
     <B>lmfa</B> usually produces such an <FONT size="+1"><tt>RSMH</tt></FONT> for <I>sp</I> orbitals in
     transition metals.  For example, rsmh for the <I>d</I> orbital in the
     copt example above are pretty close to optimum for the solid.
     <BR>

<LI> Occupied <I>d</I> orbitals are atomic-like and require a small <FONT size="+1"><tt>RSMH</tt></FONT>.  Usually
     the value produced by <B>lmfa</B> is pretty close to the optimum one.
     <BR>

<LI> Assembling the hamiltonian is more efficient when you use the same rsmh,eh
     for consecutive <I>l</I> channels, because it can better vectorize loops.
     However, there is a tradeoff in that the basis set may be poorer when
     it is not tailored.  If only a mRy or so is gained by splitting rsmh,eh,
     for consecutive <I>l's</I>, it is preferable to keep them the same.

</OL>

<BR>
Starting with v6.14, <strong>lmf</strong> offers some ability to find
optimum choices for parameters <FONT size="+1"><tt>RSMH</tt></FONT> and <FONT size="+1"><tt>EH</tt></FONT>. The optimization proceeds
by finding values for the <FONT size="+1"><tt>RSMH</tt></FONT> and <FONT size="+1"><tt>EH</tt></FONT> you specify which minimize the
total energy.  The optimization is a little primitive but seems to
work well.  There is a tutorial <A href="FPoptbas.html">FPoptbas.html</A> showing how to invoke
<strong>lmf</strong> to automatically find a reasonably optimized
minimal basis in SrTiO3.

<H2><A name="floatingorbitals"></A></H2>
<BR> <P> <FONT color="#bb3300"><I>Floating orbitals</I>.</FONT>

To specify a floating orbital, add a new species with atomic number 0
and augmentation radius 0.  No augmentation parameters (<FONT
size="+1"><tt>LMXA,LMXL,</tt></FONT> etc) are used, bu you must also
specify the envelope function.  In the following example
<pre>
 SPEC 
     .. 
     ATOM=X Z=0 R=0 LMX=1 RSMH=1.5,1.5,1.5,1.5 EH=-.3,-.3,-.3,-.3
</pre>

species <FONT size="+1"><tt>X</tt></FONT> is specified as a floating orbital.  Although
parameters <FONT size="+1"><tt>RSMH</tt></FONT> and <FONT size="+1"><tt>EH</tt></FONT> 
are supplied for <I>spdf</I> orbitals, the <FONT size="+1"><tt>LMX=1</tt></FONT> 
token specifies that just <I>s</I> and <I>p</I> functions are to used.

In the SITE category, you specify where to place the orbitals, e.g.
<pre>
 SITE
     .. 
     ATOM=X POS=1/2 1/2 1/2
     ATOM=X POS=3/4 3/4 3/4
</pre>


You can in principle put the floating orbitals anywhere at all.  In
practice, it makes little sense to put them too close to the regular
orbitals or too close to each other, numerical instabilities can
arise if you do so.

For an example that illustrates the floating orbitals, try elemental
Te, which has a very open structure:
<pre>
   fp/test/test.fp te
</pre>
Nine sites with floating sp orbitals are added to the three actual
atoms, increasing the basis from from 39 to 75 orbitals.  The test
shows that addition of the floating orbitals lowers the total energy
by 4 mRy and minimally affects the forces.  (Note that this test also
includes APWs, which in this case offers a more convenient and efficient
method to make the basis complete.)

<H2><A name="localorbitals"></A></H2>
<BR> <P> <FONT color="#bb3300"><I>Local orbitals</I>.</FONT>

 The local orbital is specified with the <FONT
 size="+1"><tt>PZ=</tt></FONT> token, which is its analog of the <A
 href="lmto.html#section3"> continuously variable principal quantum
 number <FONT size="+1"><tt>P</tt></FONT></A>.  Conventional local
 orbitals may be used either for semicore or high-lying states.  In
 this case, specify the integer part of <FONT
 size="+1"><tt>PZ</tt></FONT> must either as one less than the integer
 part of <FONT size="+1"><tt>P</tt></FONT> for a semicore state or one
 greater for a high-lying state.  As for the fractional part, it is
 recommended that you specify 0.9 or so for the semicore case; while
 for high-lying states 0.3 or so is reasonable (a larger number
 specifies a higher linearization energy).

  <P> <FONT color="#005500"><I>Extended local orbitals:</I>.</FONT>
  Extended local orbitals have a smooth Hankel envelope attached to
  the radial wave function which spills out into the interstitial.
  The tail is matched value and slope to the radial function, but
  <strong>lmf</strong> will try to vary both <FONT
  size="+1"><tt>RSMH</tt></FONT> and <FONT
  size="+1"><tt>EH</tt></FONT> to match the kinetic energy as well.
  To specify that a semicore local orbital is of the extended type,
  simply add 10 to token <FONT size="+1"><tt>PZ=</tt></FONT>.

  <P> <FONT color="#005500"><I>Hazard for conventional local
  orbitals:</I>.</FONT> Conventional local orbitals applied to
  semicore states must rely on the valence envelope functions for
  representation of these deep states in the interstitial.  Therefore be
  advised to choose the smoothing radius for at least one envelope
  function rather small (typically around 1 or so), with a
  corresponding energy rather negative (typically around -1 or so).

  <P> <FONT color="#005500"><I>Hazard for extended local orbitals:</I>.</FONT> 

  It may happen that <strong>lmf</strong> cannot find any reasonable
  combination of <FONT size="+1"><tt>RSMH</tt></FONT> and <FONT
  size="+1"><tt>EH</tt></FONT> that matches the envelope's value and
  slope to the radial wave functions.  Since this matching condition is required,
  <strong>lmf</strong> stops with a message like
<pre>
    Fit local orbitals to sm hankels, species Sr, rmt=3.1
     l  type    Pnu      Eval        K.E.       Rsm       Eh      Q(r>rmt)    Fit K.E.
    Exit -1 : mtchre : failed to match phi'/phi=-5.127 to envelope, l=1
</pre>
  This usually happens when the fractional part of <FONT
  size="+1"><tt>PZ</tt></FONT> is too large (it was set to 0.98 in
  this case).  The solution is to change the fractional part of <FONT
  size="+1"><tt>PZ</tt></FONT>.  

 <P> <strong>lmf</strong> will try to vary both <FONT size="+1"><tt>RSMH</tt></FONT> and <FONT size="+1"><tt>EH</tt></FONT>
  to match (value,slope,kinetic energy), and it may happen that <strong>lmf</strong> can match the first two but not 
  the kinetic energy of the envelope to the radial wave function.
  This doesn't cause an error but it can reduce the accuracy of the
  orbital.  It is seen in the following output:

<pre>
     l  type    Pnu      Eval        K.E.       Rsm       Eh      Q(r>rmt)    Fit K.E.
     1  low    4.900   -1.708558   -0.923699   1.00000  -0.30708   0.07250   -0.289954*
</pre>

The true and fit K.E. differ.  The solution is again to alter PZ (4.90
in this case).  It is not essential, however.  If <FONT
size="+1"><tt>PZ</tt></FONT> floats in the course of a self-consistent
cycle, <strong>lmf</strong> usually picks a good value for PZ.

<P>

 At present, the <FONT
 size="+1"><tt>P</tt></FONT> with the higher principal quantum number
 is automatically frozen to its input value (that is, <FONT
 size="+1"><tt>P</tt></FONT> is frozen when <FONT
 size="+1"><tt>PZ</tt></FONT> is a semicore orbital, and <FONT
 size="+1"><tt>PZ</tt></FONT> is frozen when it is a high-lying
 orbital).  Thus, only the <FONT size="+1"><tt>P</tt></FONT>
 corresponding to the deepest state is allowed to float (see <A
 href="lmto.html#section3"> IDMOD description</A>).

<H2><A name="apw"></A></H2>
<BR> <P> <FONT color="#bb3300"><I>APWs</I>.</FONT> 

It is very simple to add augmented plane waves to the basis.
Add these tokens to category <FONT size="+1"><tt>HAM</tt></FONT>:

<pre>
        PWMODE=1 PWEMIN=# PWEMAX=#
</pre>


<FONT size="+1"><tt>PWMODE</tt></FONT>=1 tells <b>lmf</b> to include
APWs in the basis (The default is <FONT
size="+1"><tt>PWMODE</tt></FONT>=0 which means no APWs will be
included).  <FONT size="+1"><tt>PWEMAX</tt></FONT> is an energy
cutoff (in Rydbergs): G-vectors whose energy fall below <FONT
size="+1"><tt>PWEMAX</tt></FONT> will be included.  Because MTOs can
also be included, and they already rather accurately describe low-energy states,
it makes sense to include a minimum-energy cutoff as well,
i.e. include APWs that satisfy <FONT size="+1"><tt>PWEMIN < E < PWEMAX
</tt></FONT>.

<P>

There are two further considerations the user should take note of.
First, when many APWs are used, the overlap matrix can become nearly
singular.  When this happens the secular matrix becomes unstable,
which can lead to somewhat wrong or even nonsensical results.  To
protect against this, add to <FONT size="+1"><tt>HAM</tt></FONT> the
following:
<pre>
 OVEPS=#
</pre>

If # is a positive number, <b>lmf</b> will diagonalize the overlap matrix,
and eliminate the subspace which has eigenvalues smaller than #.  This
procedure will largely eliminate these instabilities.  A reasonable
choice is <FONT size="+1"><tt>OVEPS=1E-7</tt></FONT>.

<P>

Second, high-energy APWs oscillate much more rapidly than MTOs do.
That means polynomials of higher order may be needed to augment the envelope
functions.  The order of polynomial cutoff is set by (species-specific) token
<FONT size="+1"><tt>KMXA</tt></FONT> in
the <FONT size="+1"><tt>SPEC</tt></FONT> category. When using  
<FONT size="+1"><tt>PWEMAX</tt></FONT> larger than 3 or 4,
take care that the results are converged wrt
<FONT size="+1"><tt>KMXA</tt></FONT>.  As a rough guide, the
dimensionless integer <FONT size="+1"><tt>KMXA</tt></FONT> should be
approximately as large as <FONT size="+1"><tt>PWEMAX</tt></FONT> is
(in Rydbergs).

<P>

You can see the interplay between MTOs and APWs in the Figure, which
shows the total energy <i>E</i> in SrTiO<sub>3</sub> as a function of the
number of APWs included in the basis. (The upper cutoff corresponds to
<FONT size="+1"><tt>PWEMAX</tt></FONT>=10.) 

<IMG
 WIDTH="361" HEIGHT="312" ALIGN="LEFT" BORDER="0"
 SRC="srtio3.png"
 ALT="\includegraphics[angle=0,width=.65\textwidth,clip]{srtio3.eps}">

<P>

Curve "1" includes just
<i>s</i> and <i>p</i> MT orbitals on the O sites, sufficient for a
crude representation of the valence bands (no orbitals are available
for conduction bands).  Also included are local orbitals to represent
Sr 4<i>p</i> and Ti 3<i>p</i> semicore states, as they are too extended to
be appoximated as core states.  A large number of APWs is need to get a good
total energy:  something like 150 APWs are needed to converge 
<i>E</i> to within about 50 mRy of the converged result.  (Had the O
<i>s</i> and <i>p</i> not been present, many more APWs would have been
required.)  

<P>

Curve "4" corresponds to an extreme tight-binding basis, consisting of
<i>sp</i> orbitals on Sr,
<i>spd</i> orbitals on Ti (the conduction band is mainly Ti <i>d</i>),
and <i>sp</i> orbitals on O.  The total energy of the MTO basis alone
(no APWs) is rather crude --- more than 200 mRy underbound.  However,
only 25 orbitals (plus 6 for the local orbitals) are included in this
basis.  The energy drops rapidly as low-energy APWs are included:
adding about 40 APWs is suffcient to converge <i>E</i> to about 50 mRy.

As more APWs are added,  as more APWs are added, the gain in energy
becomes more gradual; indeed convergence is very slow for large <i>E</i>.

<P>

Curve "5" differs from curve "4" only in that a Sr <i>d</i> orbital
was added.  With the addition of these 5 orbitals, the MTO-only basis
is already rather reasonable.  This would be the smallest acceptable
MTO-basis.  As in the Curve "4" basis, there is initially a rapid gain
in energy as the first few APWs are added, followed by a progessively
slower gain in energy as more APWs are added.  The blue curve close to
Curve "5" uses the same MTO basis, but <FONT
size="+1"><tt>PWEMIN</tt></FONT>=1.  Interestingly, the blue and black
curves almost superimpose on each other.  It means the gain in energy 
seems to depend on the number of APWs, but not details of their shape.

<P>

Curve "6" is a standard LMTO minimum basis: <i>spd</i> orbitals on all
atoms.  Comparing Curve "5" or Curve "6" to Curve "1" shows that the
MTO basis is vastly more efficient than the APW basis in converging
the total energy.  This is true until a minimum basis is reached.
Beyond this point, the gain APWs and more MTOs improve the total
energy with approximately the same efficiency, as the next tests show.


<P>

Curve "8" is a standard LMTO larger basis: <i>spdspd</i> orbitals on
Sr and Ti, and <i>spdsp</i> on O.  Comparing curves "6" and '8" shows
that the efficacy of any one orbital added to to the standard MTO
minimum basis is rather similar in the APW and MTO cases.  Thus,
increasing the MTO basis from 51 to 81 orbitals in the MTO basis
lowers the energy by 33 mRy; adding 33 APWs to the minimum basis
(curve "6") lowers the energy by 36 mRy.  APWs, however, have the
advantage that they are simpler to include.

<P>

Curve "11" enlarges the MTO basis still more, with Sr: <i>spdspd</i>,
Ti: <i>spdspd</i>, O: <i>spdspd</i>.  Also local orbitals are used to
represent the high-lying Ti 4<i>d</i> and O 3<i>s</i> and 3<i>p</i>
states.  For occupied states, these orbitals have little effect,
but they are important for unoccupied states higher than about 2 Ry.

Comparing the last curves, it appears that finally <i>E</i> is at last
converging to an absolute minimum at around -2.760 Ry.  The instabilities
in the overlap matrix mentioned above are largely controlled, but there
is a small uncertainty in the absolute energy on the order of 1 mRy.

It is clear even from inspection of Curve "1", that a pure APW basis
would require an extremely large number of PWs to reach this level of
convergence.

<A name="ldaplusu"></A>
<P> <FONT color="#bb3300"><I>LDA+U</I>.</FONT>

To get started, LDA+U needs in addition to the LDA input two modifications:
<UL>
          Parameters <I>U</I> and <I>J</I> 
<BR>
          Density-matrix, stored in file <FONT size="+1"><tt>dmats.ext</tt></FONT>.
         (In the course of self-consistency, it is made by <B>lmf</B>.)
</UL>

For the first, a line must be added to the input file for each species
that contains
<i>U</i> and <i>J</i>.  Consider the example taken from the ErAs test case:
<pre>
        IDU= 0 0 2 2 UH= 0 0 0.1 0.632 JH= 0 0 0 0.055
</pre>
The IDU token tells <B>lmf</B> that no <I>U</I> is to be added to the
<I>s</I> or <I>p</I> channels, but that a <I>U</I> is to be added to
the <I>d</I> and <I>f</I> channels.  IDU=2 specifies LDA+U functional
style 2; this is the "Fully Localized Limit" described in Liechtenstein, PRB 52,
R5467 (1995)).  IDU=1 specifies the "Around Mean Field" functional (Petukhov, PRB 67,
153106 (2003)).

<I>U</I>=0.1 Ry is included on the
<I>d</I> orbital, and <I>U</I>=0.632 is included on the <I>f</I>
orbital.  Additionally <I>J</I>=0.055 is put on the <I>f</I> orbital.

<P>

The density-matrix is read from and written to a file
<FONT size="+1"><tt>dmats.ext</tt></FONT>.  Two density-matrices (1 for each spin) are written
to this file in a (2<I>l</I>+1) by (2<I>l</I>+1) block for each
<I>l</I> block for which a <I>U</I> is defined.  
<FONT size="+1"><tt>dmats.ext</tt></FONT> is
an ASCII file which you can read, and it's quite useful to interpret
what's going on.  The diagonal parts are the occupation numbers and
are the most important.  Note that the file may be stored in either
spherical harmonics or real harmonics, depending on how 
<FONT size="+1"><tt>SHARM=</tt></FONT> is set in the
<FONT size="+1"><tt>OPTIONS</tt></FONT> category.

<P>

Usually it is too tedious to supply 
<FONT size="+1"><tt>dmats.ext</tt></FONT> itself, especially since you
typically won't know what to choose.  Instead, you can supply an
"occupation numbers" file <FONT size="+1"><tt>occnum.ext</tt></FONT>,
which is a starting guess for the density-matrix (its diagonal part).
<FONT size="+1"><tt>occnum.ext</tt></FONT> has one line of
(2<I>l</I>+1) numbers for the occupation numbers of the first spin,
following by a line with the occupation numbers for the second spin.

In this test:
<pre>
   fp/test/test.fp eras
</pre>

the script assumes a particular starting spin configuration through
the occupation number file <FONT size="+1"><tt>occnum.eras</tt></FONT>
it uses.  Er has 11 <i>f</i> electrons, 7 of which go into the
majority channel and 4 into the minority channel.  There is some
choice in which <i>m</i> states to fill and which to keep empty.  A key point
is that the self-consistent solution you end up with will depend on this
choice.   The
ErAs test uses the following input file for <FONT
size="+1"><tt>occnum.eras</tt></FONT> :
<pre>
  0 0 0 0 0
  0 0 0 0 0
  1 1 1 1 1 1 1
  0 1 1 1 1 0 0
</pre>

The first and second lines are occupation numbers for the majority and
minority <i>d</i> channel; the third corresponds to the majority
<i>f</i> channel where all states are taken to be filled.
The last line corresponds to the minority <i>f</i> channel.  In this case,
<i>m</i>=-2,-1,0,1 are filled and <i>m</i>=-3,2,3 are empty.
As the script notes, different choices of starting occupation numbers
lead to different self-consistent solutions.  The one with the lowest
energy is that which satisfies Hund's rule (<i>m</i>=0,1,2,3 filled
and <i>m</i>=-3,-2,-1 empty).

<P>

The occupation numbers are by default correpond to spherical
harmonic representations of <i>Y</i><sub><i>lm</i></sub>.
If you want to define the occupation numbers in real harmonics, put
<pre>
% real
</pre>
on the first line of <FONT size="+1"><tt>occnum.ext</tt></FONT>.

<P>

Note that in the LDA case, complete information is contained in the
density, stored in file <FONT size="+1"><tt>rst.ext</tt></FONT>.  In
the LDA+U case, complete information is contained in the combination of
<FONT size="+1"><tt>rst.ext</tt></FONT> and <FONT
size="+1"><tt>dmats.ext</tt></FONT>.

<P>

The density-matrix is at present mixed independently of the charge,
with linear mixing.  (This will likely change in future) To do the
mixing there is a special-purpose parameter in the <A href="#MIX">MIX</A> category

<pre>
  UMIX=# 
</pre>

# is a parameter between 0 and 1.  It plays the same role for the
density-matrix as the <A href="lmto.html#MIX">mixing beta</A> 
plays for the mixing of the regular density.

<P> There is additionally a tolerance parameter

<pre>
  TOLU=#    
</pre>

that tells <B>lmf</B> to stops mixing dmats when its rms change falls
below <FONT size="+1"><tt>TOLU</tt></FONT>.  Usually it's not
necessary, and setting <FONT size="+1"><tt>TOLU</tt></FONT>=0 (or
leaving it out) means it plays no role.

<H2><A name="corehole"></A></H2>
<P> <FONT color="#bb3300"><I>Partially occupied core holes</I>.</FONT> 

Calculations involving partial core hole occupancy are useful in the
context of Slater transition-state theory, which undoes most of the
error in the LDA description of the core hole eigenvalue (see
example J. Phys. Cond. Mat. 12, 729 (2000)).

You specify which orbital in which species is to be treated as a
partially occupied by core by adding a token <tt> C-HOLE= </tt> to the
<tt> SPEC </tt> category.  You also have to specify what the partial
occupation is, which you do with token <tt> C-HOLE= </tt>.  An example
is the N 1<I>s</I> core.  A core hole of 1 electron can be put in by
<pre>
     C-HOLE=1s C-HQ=-1
</pre>
Note the sign of the charge.  The number refers to the excess electron charge.
to put in a hole, use a negative charge.  As an illustration in CrN,
you can run the test case
<pre>
     fp/test/test.fp crn
</pre>

<P> Core holes are also useful as an approximate workaround in the LDA
context to deal with (almost) nonbonding <I>f</I> electrons.  In most
4<I>f</I> systems, the <I>f</I> states get shifted away from the Fermi
level, even though the LDA typically is unable to do this (except for
Gd), because it lacks a nonlocal exchange as in LDA+U.  An approximate
workaround is to treat the 4<I>f</I> electrons as core.

For Gd in particular, the 7 majority states should be filled, while
the 7 minority states empty.  Thus a core hole of -7 is required, but
it is necessary to further specify the spin polarization of the core.
This is accomplished with a second argument.  For the Gd case (4<I>f</I> core, -7
excess electrons, with the core magnetic moment +7),  use
<pre>
    C-HOLE=4f C-HQ=-7,7 
</pre>
As an illustration, try
<pre>
    fp/test/test.fp gdn
</pre>

<H2><A name="rsteditor"></A></H2>
<P> <FONT color="#bb3300"><I>Restart file editor</I>.</FONT> 

Not documented.


<H3><A name="spectable"></A>
lmf-specific tokens in the <tt> SPECIES </tt> category
</H3>

 The following table documents the remaining <B>lmf</B>-specific input for SPEC
 category.  Most are optional.  Symbols in the ``Arguments'' column refer
 to the data expected following the token.  `l', `i', `r', and `c' refer
 to logical, integer, floating-point and character data, respectively.

<TABLE border="1">
  <TBODY>
   <TR>
     <TD bgcolor="#000000"><FONT color="#ffffff">Token </FONT></TD>
     <TD bgcolor="#000000"><FONT color="#ffffff">Arguments
     <TD bgcolor="#000000"><FONT color="#ffffff">Explanation</FONT></TD>
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">LMXA=</TD>
   <TD bgcolor="#feffc1"><PRE> i</PRE></TD>
   <TD bgcolor="#c5c1ff">
           angular momentum <I>l</I>-cutoff for projection of wave functions tails
           centered at other sites in this sphere.
           Must be at least the basis <I>l</I>-cutoff (specified by LMX=).
   <BR>
           Default:NL-1
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">KMXA=</TD>
   <TD bgcolor="#feffc1"><PRE> i</PRE></TD>
   <TD bgcolor="#c5c1ff">
           polynomial cutoff for projection of wave functions in sphere.
           Smoothed Hankels are expanded in polynomials around other sites
           instead of Bessel functions as in the case of normal Hankels.
   <BR>
           Default:3
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">RSMA=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           smoothing radius for projection of smoothed Hankel tails into
           augmentation spheres.  sm-Hankels are expanded in polynomials by
           integrating with Gaussians at that site.  Thus RSMA=0 =>
           sm-Hankels Taylor series expansion about the origin.  For large
           KMXA the choice is irrelevant, but RSMA is best chosen that
           maximizes the convergence of sm-Hankels with KMXA.
   <BR>
           Default: R*0.4, with R=augmentation (muffin-tin) radius.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">LMXL=</TD>
   <TD bgcolor="#feffc1"><PRE> i</PRE></TD>
   <TD bgcolor="#c5c1ff">
           angular momentum <I>l</I>-cutoff for explicit representation of local
           charge on a radial mesh.
   <BR>
           Default:NL-1
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">RSMG=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           smoothing radius for Gaussians added to sphere densities to correct
           multipole moments needed for electrostatics.  Value should be as
           large as possible but small enough that the Gaussian doesn't spill
           out significantly beyond rmt.
   <BR>
           Default: R/4, with R=augmentation (muffin-tin) radius.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">LFOCA=</TD>
   <TD bgcolor="#feffc1"><PRE> i</PRE></TD>
   <TD bgcolor="#c5c1ff">
           prescribes how core is treated.
<DL><DD>
           0 confines core to within RMT.  Usually the least accurate.
<DD>
           1 treats the core as frozen but lets it spill into the interstitial
<DD>
           2 same as 1, but interstitial contribution to vxc treated perturbatively.
</DL>
           Default: 1
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">RFOCA=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           smoothing radius fitting tails of core density.  A large radius
           produces smoother interstitial charge, but less accurate fit.
   <BR>
           Default: R*0.4, with R=augmentation (muffin-tin) radius.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">RSMFA=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           smoothing radius for tails of free-atom charge density.
           Irrelevant except first iteration only (non-self-consistent harris).
           A large radius produces smoother interstitial charge, but somewhat
           less accurate fit.
   <BR>
           Default: R/2, with R=augmentation (muffin-tin) radius.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">C-HOLE=</TD>
   <TD bgcolor="#feffc1"><PRE> c</PRE></TD>
   <TD bgcolor="#c5c1ff">
           channel for core hole.  You can force partial core occupation.
           Syntax consists of two characters, the principal quantum number and
           the second one of `s',`d',`d',`f' for the <I>l</I> quantum number, e.g. `2s'
           See <A href="fp.html#corehole">Partially occupied core holes</A> for description and examples. <BR>
           NB: this also works with the ASA code <b>lmf</b>.
   <BR>
           Default: nothing
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">C-HQ=</TD>
   <TD bgcolor="#feffc1"><PRE> r[,r]</PRE></TD>
   <TD bgcolor="#c5c1ff">
           First number specifies the number of electrons to remove from the <I>l</I>
           channel specified by C-HOLE=.  <BR>
           Second (optional) number specifies the hole magnetic moment.  <BR>
           See <A href="fp.html#corehole">Partially occupied core holes</A> for description and examples.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">P=</TD>
   <TD bgcolor="#feffc1"><PRE> r,r,...</PRE></TD>
   <TD bgcolor="#c5c1ff">
           starting values for potential functions, one for each of <I>l</I>=0..LMXA
   <BR>
           Default: taken from an internal table.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">PZ=</TD>
   <TD bgcolor="#feffc1"><PRE> r,r,...</PRE></TD>
   <TD bgcolor="#c5c1ff">
           starting values for local orbital's potential functions,
	   one for each of <I>l</I>=0..LMX.  Setting PZ=0 for any
	   <I>l</I> means that no local orbital is specified for this
	   <I>l</I>.  Each integer part of PZ must be either one less
	   than P (semicore state) or one greater (high-lying state).
   <BR>
           Default: zero
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">Q=</TD>
   <TD bgcolor="#feffc1"><PRE> r,r,...</PRE></TD>
   <TD bgcolor="#c5c1ff">
           charges for each <I>l</I>-channel making up free-atom density
   <BR>
           Default: taken from an internal table.
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">MMOM=</TD>
   <TD bgcolor="#feffc1"><PRE> r,r,...</PRE></TD>
   <TD bgcolor="#c5c1ff">
           magnetic moments for each <I>l</I>-channel making up free-atom density
           Relevant only for the spin-polarized case.
   <BR>
           Default: 0,0,...
   </TR>
  </TBODY>
</TABLE>
<BR><BR>

<H2><A name="HAM"></A></H2>
<LI>Category `HAM' has several <B>lmf</B>-specific additions to the
    <A href="lmto.html#HAM"</A>HAM</A> category documented for the ASA.

<TABLE border="1">
  <TBODY>
    <TR>
      <TD bgcolor="#000000"><FONT color="#ffffff">Token </FONT></TD>
      <TD bgcolor="#000000"><FONT color="#ffffff">Arguments
      <TD bgcolor="#000000"><FONT color="#ffffff">Explanation</FONT></TD>
    </TR>

    <TR>
      <TD bgcolor="#c5c1ff">GMAX=</TD>
      <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
      <TD bgcolor="#c5c1ff">an energy cutoff used to create mesh for the interstitial
           density.  A uniform mesh is created that best `fits into'
           cutoff gmax.<BR>
      <FONT color="#ff3366">Required input; alternatively information may be specified by the following token.</FONT></TD>
    </TR>

    <TR>
      <TD bgcolor="#c5c1ff">FTMESH=</TD>
      <TD bgcolor="#feffc1"><pre> i1 [i2 i3]</pre></TD>
      <TD bgcolor="#c5c1ff">
           the number of divisions specifying the uniform mesh density
           along the three reciprocal lattice vectors.  The second
           and third arguments default to the value of the first one,
           if they are not specified.
      <FONT color="#ff3366">Required input; alternatively information may be specified by the preceding token.</FONT>.</TD>
    </TR>

    <TR>
      <TD bgcolor="#c5c1ff">TOL=</TD>
      <TD bgcolor="#feffc1"><pre> r</pre></TD>
      <TD bgcolor="#c5c1ff">specifies the precision to which the wave
      functions are expanded in a Fourier expansion.<BR>
      Default:10<sup>-6</sup></TD>
    </TR>

    <TR>
      <TD bgcolor="#c5c1ff"> FRZWF=</TD>
      <TD bgcolor="#feffc1"><pre> l</pre></TD>
      <TD bgcolor="#c5c1ff">specifies whether the augmented part of the wave functions be
           frozen.  Normally they are updated as the potential changes,
           but with FRZWF=t the potential used to make augmentation wave
           functions is frozen at what is read from the restart file (or
           free-atom potential if starting from superposing free atoms).<BR>
      Default:F</TD>
    </TR>

    <TR>
      <TD bgcolor="#c5c1ff">FORCES=</TD>
      <TD bgcolor="#feffc1"><pre> i</pre></TD>
      <TD bgcolor="#c5c1ff"> controls how forces are to be calculated, and how the
           second-order corrections are to be evaluated.  Through the
           variational principle, the total energy is correct to
           second order in deviations from self-consistency, but
           forces are correct only to first order.  To obtain forces
           to second order, it is necessary to know how the density
           would change with a (virtual) displacement of the
           core+nucleus, which requires a linear response treatment.
           <B>lmf</B> estimates this changes in one of two ways:
           <OL>
             <LI> the free-atom density is subtracted from the total
             density for nuclei centered at the original position and
             added back again at the (virtually) displaced position.
             <BR>For this ansatz, use FORCES=1.

             <LI> the core+nucleus is shifted and screened assuming a
             Lindhard dielectric response.
	     <BR>For this ansatz, use FORCES=12.  You also must
	     specify ELIND, below.
           </OL>
      Default:0</TD>
    </TR>
    <TR>
      <TD bgcolor="#c5c1ff"> ELIND= </TD>
      <TD bgcolor="#feffc1"><pre> r</pre></TD>
      <TD bgcolor="#c5c1ff">
           a parameter in the the Lindhard response function, namely the
           Fermi level for a free-electron gas relative to the bottom of
           the band.  You can specify this energy directly, by using a
           positive number for the parameter.  If you use instead a
           negative number, the program will choose a default value from
           the total number of valence electrons and assuming a
           free-electron gas, scale that default by the absolute value of
           the number you specify.  If you have a simple sp bonded system,
           the default value is an excellent choice.  If you have <I>d</I>
           or <I>f</I> electrons, it tends to overestimate the response.
           Use a something smaller, e.g. ELIND=-0.7.  ELIND is used in
           three contexts:
          <DL>
          <DD>
             (1) in the force correction term; see FORCES= above
          <DD>
             (2) to estimate a self-consistent density from the input and
                 output densities after a band pass
          <DD>
             (3) to estimate a reasonable smooth density from a starting
                 density after atoms are moved in a relaxation step.
          </DL>
      Default:-1<BR>
      </TD>
    </TR>
    <TR>
      <TD bgcolor="#c5c1ff"> SIGP: </TD>
      <TD bgcolor="#feffc1"><pre> r</pre></TD>
      <TD bgcolor="#c5c1ff">
           A collection of parameters use to interpolate the self-energy
           Sigma. Used in conjuction with the GW package. See gw.txt
           for description.
      Default:not used<BR>
      </TD>
    </TR>

  </TBODY>
</TABLE>
<BR><BR>

<H2><A name="MIX"></A></H2>
<LI>Category `MIX' has the following <B>lmf</B>-specific tokens

<TABLE border="1">
  <TBODY>
   <TR>
     <TD bgcolor="#000000"><FONT color="#ffffff">Token </FONT></TD>
     <TD bgcolor="#000000"><FONT color="#ffffff">Arguments
     <TD bgcolor="#000000"><FONT color="#ffffff">Explanation</FONT></TD>
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">CONV=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           maximum allowed change in energy between iterations to reach
           self-consistency.
   <BR>
           Default: 10<sup>-4</sup>
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">CONVC=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           maximum allowed change in output-input density to reach
           self-consistency.
   <BR>
           Default: 10<SUP>-4</SUP>
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">UMIX=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           (specific to LDA+U) linear <A href="lmto.html#MIX">mixing beta</A> for density matrix
   <BR>
           Default: 1
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">TOLU=</TD>
   <TD bgcolor="#feffc1"><PRE> r</PRE></TD>
   <TD bgcolor="#c5c1ff">
           (specific to LDA+U) tolerance in density-matrix.  It is not
           updated in the course of self-consistency when the RMS change in
           output-input density matrix falls below TOLU
   <BR>
           Default: 0
   </TR>


  </TBODY>
</TABLE>
<BR><BR>

<LI>Category `DYN' is used for molecular statics and dynamics.
    (dynamics is not implemented).

<TABLE border="1">
  <TBODY>
   <TR>
     <TD bgcolor="#000000"><FONT color="#ffffff">Token </FONT></TD>
     <TD bgcolor="#000000"><FONT color="#ffffff">Arguments
     <TD bgcolor="#000000"><FONT color="#ffffff">Explanation</FONT></TD>
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">MSTAT:</TD>
   <TD bgcolor="#feffc1"><PRE> l1,l2,r3,r4,r5</PRE></TD>
   <TD bgcolor="#c5c1ff">
           a collection of parameters specifying how site positions are
           relaxed given the internuclear forces
           <pre> l1=0: skip relaxation of atom positions<BR>    4: relax using conjugate-gradients<BR>    5: relax using Fletcher-Powell<BR>    6: relax using Broyden</pre>
           <pre> l2 T: read hessian from disk<BR>    F: assume initial hessian=1</pre>
           <pre> r3 convergence criterion in displacements<BR>    r3=0 means not to use this criterion</pre>
           <pre> r4 convergence criterion in forces<BR>    r4=0 means not to use this criterion<BR>    NOTE: you should have either r3>0 or r4>0</pre>
           <pre> r5 initial step length, in units of the lattice parameter a</pre>
           Default: no relaxation
   </TR>

   <TR>
   <TD BGCOLOR="#C5C1FF">NIT=</TD>
   <TD bgcolor="#feffc1"><PRE> i</PRE></TD>
   <TD bgcolor="#c5c1ff">
           maximum number of relaxation steps
   <BR>
           Default: 1
   </TR>

  </TBODY>
</TABLE>

<HR>
In molecular statics calculations, <B>lmf</B> will iterate at a fixed
set of positions until it meets the self-consistency criteria (CONV=
and CONVC= above), or until the maximum number of iterations is
reached in the self-consistency cycle (START NIT=).  then it will move
the atoms using the forces it has using a conjugate-gradients,
Fletcher-Powell, or Broyden technique.  This cycle is repeated until
the forces and/or displacements are converged to the specified
tolerance, or until atoms have been moved NIT times.

<BR><BR>

See <B><A href="force-minimization.html"><tt>force-minimization.html</tt></A></B>
for more description on relaxation of atoms, and strategies to
accelerate this procedure.



<P>
<H2><A name="section4"><FONT size="+1">
4.  Command-line switches
</FONT></A></H2>
See <A href="lmto.html#section9">lmto.html</A>
for a general discussion of command-line switches, and those common
to most programs.
<BR><BR>

<STRONG>Switches special to lmf</STRONG>

<PRE>
 --rs=#1,#2,#3,#4,#5  tells <B>lmf</B> what parts to read from the rst file.
                  #1= 0: do not read the restart file, but overlap
                         free-atom densities
                      1: read restart data from binary rst.ext
                      2: read restart data from ascii rsta.ext
                      3: same as 0, but also tells <B>lmf</B> to overlap
                         free-atom densities after a molecular statics
                         or molecular dynamics step.
                     11 or 12: same as 1, or 2 but additionally adjust
                         the mesh density for shifts in site positions
                         relative to those used in the generation of the
                         restart file.  Note: see --rs switch #3 below
                         for which site positions the program uses.  The
                         same principle for adjusting the density is used
			 as in computing corrections to the Helman-Feynman
                         forces; see token FORCES= in <A href="#HAM">HAM</A>.

                      Additionally, if you read from rst or rsta
                      you can add 10 to #1

                  #2= 0: exactly as #1, but except switches apply to
                         writing. Value zero suppresses writing.
                      1: write binary restart file rst (default)
                      2: write ascii restart file rsta
                      3: write binary file to rst.#, where # = iteration number

                  #3= 0: read site positions from restart file,
                         overwriting positions from input file
                         (this is the default)
                      1: ignore positions in restart file

                  #4= 0: read guess for Fermi level and window from
                         restart file, overwriting data from input
                         file.  This data is needed when the BZ
                         integration is performed by sampling.
			 (this is the default)
                      1: ignore data in restart file

                  #5= 0: read linearization pnu from restart file,
                         overwriting data from input file
			 (this is the default)
                      1: ignore pnu in restart file

                  The <B>lmf</B> output around lines
                     iors  : read restart file (binary, mesh density)
                  tells you what bits are read and what is ignored.
                  If not specified, <B>lmf</B> defaults to --rs=1,1,0,0,0

 --rdbasp[:fn]    tells the program to read basis parameters (input using
                  tokens RSMH=,EH=,RSMH2=,EH2=,PZ= for each atom in the
                  ctrl file) from file `basp' (or optionally file `fn').
                  Parameters input in this mode supersedes parameters
                  read from the ctrl file.  You can specify none, or any
                  set of the sets (RSMH=,EH= ; RSMH2=,EH2=; PZ=)
                  for each species.  Parameters not specified here
                  default to what was specified in the ctrl file.

                  This switch (and --optbas described below) is a useful
                  way to get started when you don't know what to choose
                  for parameters in the basis set.

 --optbas[:sort][:spec=name[,rs][,e][,l=###]...]
                  operates the program in a special mode to optimize
                  the total energy wrt the basis set. <STRONG>lmf</STRONG> makes
                  several band passes (not generating the output density
                  or adding to the save file), varying selected
                  parameters belonging to tokens RSMH= and EH= to
                  miniminize the total energy wrt these parameters.
                  Either the smoothing radius [,rs] or the energy [,e]
                  must be selected for optimization (you can select
                  both).  Select which <I>l</I> quantum numbers whose parameters
                  you want to optimize using `,l=##..', e.g. l=02 .  The
                  optimization routine is rather primitive, but it seems
                  to work reasonably well. See the <A href="FPoptbas.html">basis optimization tutorial</A>
                  for a more complete description and an example.

 --etot           is a special mode designed to be used in conjuction
                  with the GW suite.  It generates parameters for the
                  LDA total energy without disturbing the rst or mixing
                  files, or logging the energy in the save file.

 --rpos=filename  tells the program to read site positions from file
                  `filename' after the CTRL file has been read

 --wpos=filename  tells the program to write site positions to file
                  `filename' after a relaxation step.

 --band[:options] tells <B>lmf</B> to generate energy bands instead of make a
                  self-consistent calculation.  See
                  the <A href="FPtutorial.html#bands">FP tutorial</A> for an example and
                  <A href="../subs/suqlst.f">subs/suqlst.f</A> for a description of options.

 --pdos[:options] tells <B>lmf</B> to generate weights to enable calculation
                  of partial DOS inside augmentation spheres. See <B>lm</B>'s
                  <A href="lmto.html#section9lm">command-line switches</A> for description of the options.
                  NB: unlike the ASA, <B>lmf</B> does not by default save
                  information in the moments file to decompose the
                  dos into channels.  You must run <B>lmf</B> with the
		  --pdos switch set before invoking <B>lmdos</B>.

 --mull[:options] tells <B>lmf</B> to generate weights for Mulliken analysis.
                  --mull and --pdos may not be used in conjunction.  The
                  options to --mull are the same as those to --pdos;
                  see <A href="lmto.html#section9lm"><B>lm</B> command-line switches</A> for a description.

 --cls[options]   tells <B>lmf</B> to generate weights to compute matrix
                  elements and weights for core-level-spectroscopy.  See
                  <A href="../subs/suclst.f">subs/suclst.f</A> for a description of options.

 --wden[:options] writes one plane of the charge density to disk, on a
                  uniform of mesh of points.  Information for the
                  plane is specified by three groups of numbers: the
                  origin (i.e. a point through which the plane must
                  pass), a first direction vector with its number of
                  points, and a second direction vector with its
                  number of points.  Default values will be taken for
                  any of the three sets you do not specify.  The
                  density generated is the smooth density, augmented
                  by the local densities in a polynomial approximation
                  (see option core= below)

                  The options are specifications (see below) and
		  different options are separated by delimiters
		  (chosen to be `:' in this text; the delimiter
		  actually taken is the first character after `wden')

                  At present, there is no capability to interpolate
		  the smoothed density to an arbitrary plane, so you
		  are restricted to choosing a plane that has points
		  on the mesh.  Accordingly, all three groups of
		  numbers are given sets of integers, as will be
		  explained below.  Supposing your lattice vectors are
		  p1, p2 and p3, which the smooth mesh having (n1,n2,n3)
		  divisions.  Then the point (#1,#2,#3) corresponds to
		  the Cartesian coordinate
		     #1/n1 p1 + #2/n2 p2 + #3/n3 p3

                  Specify the origin (a point through which the plane
                  must pass) by
                     :o=#1,#2,#3
                  Default value: o=0,0,0.

                  Specify the direction vectors by
		      :l1=#1,#2,#3[,#4]
		      :l2=#1,#2,#3[,#4]

		  l1 and l2 specify the first and second direction
                  vectors, respectively.  The first three numbers
                  specify the orientation and the fourth specifies the
                  `length'.  #1,#2,#3 select the increments in mesh
                  points along each of the three lattice vectors that
                  define the direction vector.  Thus a direction
                  vector in Cartesian coordinates is
		     #1/n1 p1 + #2/n2 p2 + #3/n3 p3
		  The last number (#4) specifies how many points to take
		  in that direction and therefore corresponds to a length.
		  Default values:
                    l1=1,0,0,n1+1
                    l2=0,1,0,n2+1

		  Other options:

                  core=#   specifies how local densities is to be included.
                           Any local density added is expanded as
                           polynomial * gaussian, and added to the
                           smoothed mesh density.
		           #=0 includes core densities + nuclear contributions
		           #=1 includes core densities, no nuclear contributions
		           #=2 exclude core densities
		           #=-1 no local densities to be included (only interstitial)
		           #=-2 local density, with no smoothed part
		           #=-3 interstitial and local smoothed densities

                           Default: core=2

                  fn=name  specifies the file name for file I/O
                           The default name is `smrho'.

		  Example: use `:' as the delimiter, and suppose
		           n1=n2=48 and n3=120.  The specification
                             :fn=myrho:o=0,0,60:l1=1,1,0,49:l2=0,0,1,121
			   writes 'myrho.ext' a mesh (49,121) points.
			   The origin (first point) lies at (p3/2).
			   The first vector points along (p1+p2), and
			   has that length; the second vector points
			   along p3, and has that length.

 --window=#1,#2   accumulates the density in an energy window
                  specified by the limits #1,#2.  (This option is
                  intended to be used in conjunction with --wden).  If
                  invoked, <B>lmf</B> exits after a single band pass.

 --mixsig=#1[,#2] (used in conjunction with GW self-energy sigma)
                  <B>lmf</B> takes a linear combination of the
                  self-energy read from one or two files.  See the
                  <B><A href="gw.html#section1.5">GW documentation</A></B>.

 --wsig[options] (used in conjunction with GW self-energy sigma)
                  <B>lmf</B> writes self-energy to file <I>sigm2.ext</I> and exits.
                  See the <B><A href="gw.html#section1.6">GW documentation</A></B>.

 --oldvc          chooses nfp-style energy zero, which sets the cell
                  average of the potential to zero.  Normally the
                  average estat potential at the RMT boundary is
                  chosen to be the zero.  That puts the Fermi level
                  near zero, like in the ASA.
</PRE>
</BODY>
</HTML>

<BR><BR>
<H2><A name="section5"><FONT size="+1">
5. Frequently asked questions
</FONT></A></H2>


<BR><I>Q: I want to generate the total </I> DOS (SAVDOS=T, DOS=range) <I>over a specified range.
But <B>lmf</B> automatically resets the range of DOS to cover all the occupied states.
When I do NiO (with deep core), the range is very large. Is there a way to avoid this problem?</I>

<BR>
<BR><I>A</I>: This is a small bug, introduced to avoid problems in the sampling integration algorithm.
Use:
<pre>
  --no-fixef0
</pre>
and the range will not change.

<BR><BR>
<DT><B>Energy Functionals</B>

<BR>
<BR><I>Q: Should the Harris-Foulkes and Hohenberg-Kohn Sham
functionals agree at self-consistency?
</I>
<BR><BR><I>A:</I> Yes.  If they do not, either you have not reached self-consistency, 
or a tolerance is set too loosely in the evaluation of the matrix elements.
A likely culprit is that the mesh of points for the interstitial
integrals is not fine enough.

Here is part of an output where mesh of 14x14x14 divisions were used for
the interstitial density.
<PRE>

 GVLIST: gmax = 4.702 a.u. created 973 vectors of 2744 (35%)
         mesh has 14 x 14 x 14 divisions; length 0.588, 0.588, 0.668
 SGVSYM: 170 symmetry stars found for 973 reciprocal lattice vectors

...

 sugcut:  make orbital-dependent reciprocal vector cutoffs for tol= 1.00E-05
 spec      l    rsm    eh     gmax    last term    cutoff
  Te       0*   1.30  -0.10   5.220    9.30E-05     973*
  Te       1*   1.30  -0.10   5.597    4.36E-04     973*
  Te       2    1.30  -0.10   5.976    2.04E-03     973*
  Te       0*   1.40  -1.00   4.847    2.11E-05     973*
  Te       1    1.40  -1.00   5.182    9.90E-05     973*
  E        0*   1.50  -0.30   4.524    1.06E-05     847 
  E        1    1.50  -0.30   4.823    2.01E-05     973*
</PRE>

The table reflects how well each of the basis orbitals is converged in
a PW expansion, which is used to make matrix elements of the
interstitial potential.  Only the E <I>s</I> orbital is converged
below the specified tolerance (10<sup>-5</sup> ).  Convergence in the
Te <I>d</I> orbital is particularly poor; however, the Te <I>d</I> is
only minimally part of the bonding, and its poor convergence doesn't
affect the energy much.  In this case (3 Te atoms), the HF and HK
total energies differ by ~10<sup>-4</sup> Ry.  Below is a table
showing the convergence in HF energy as a function of the number of
divisions, and also the difference between HK and HF energies.
<PRE>
  n       ehf         ehf-ehk
  14 -14561.140624    0.000089
  15 -14561.133966    0.000017
  16 -14561.132309    0.000005
  20 -14561.131170    0.000000

</PRE>

<BR>
<BR><I>Q: What does this mean?</I>
<PRE>
   Exit -1 ICLBSJ: sought atom no. 1 in class 7 but only 0 atoms exist
</PRE>
<I>A:</I> Usually it means that there is an improper
synchronization between site and species data. (e.g. species are
defined that have no corresponding sites)

<BR><BR>
<DT><B>Questions about symmetry operations</B>
<BR>
<BR><I>Q: Strangely enough, <B>lmf</B> already calculates bands in bndfp and then when
making output density it runs into this problem with symprj:
</I>
<PRE>    Exit -1 symprj: no site mapped into first under op # </PRE>
<I>
Why does it not detect this problem earlier?
</I>
<BR>
<BR><I>A:</I> <B>lmf</B> has to symmetrize the output density by rotating all densities of
the same class into one site, averaging them, and restoring the
average densities by inverse rotations into each site of the class.
Until that point <B>lmf</B> doesn't need any information about classes.  This
happens when symmetry operations are specified that are incomptible
with the lattice, and should not occur if the program specifies the
symmetry operations automatically.
<BR>
<BR><I>Q:</I> what is tol in the following output line
<PRE>    SYMCRY: crystal invariant under 8 symmetry operations for tol=1e-5</PRE>
and can it be set somewhere to control how symmetry routines work?
<BR>
<BR><I>A:</I> tol is the fuzz the symmetry finder allows determine
what atom positions are related by symmetry.  You can tell <B>lmf</B> (or
other programs) to adjust the positions that are nearly coincident
under symmetry operations, but are not, by the command-line aurgment
<PRE>    --fixpos:# </PRE>
where # is the tolerance you specify for sites being `close enough' to
coincident.
