<HTML><HEAD><TITLE>GW driver package (v6.16) </TITLE><STYLE><!--P{font-size : medium;} --></STYLE></HEAD><BODY bgcolor="#cccccc">
<H1><B><FONT size="+2">GW driver package (v6.16)</FONT></B></H1>

This package coordinates with an all-electron <i>GW</i> package
written by T. Kotani, including an implementation of <b>quasiparticle
self-consistent <i>GW</i></b> (<i>QSGW</i>) theory developed by the
authors.  This package contains a driver <B>lmfgw</B> that generates
input the
<i>GW</i> code requires.  The <i>GW</i> code is a separate package,
and is not included here, though <i>GW</i> manual man033_ver1.pdf can
be found in this directory.

Additionally program <B><A href="fp.html">lmf</A></B> reads the
(energy-independent) self-energy as generated by the
<i>QSGW</i> implementation of the <i>GW</i> code.  This enables
<b>lmf</b> to do the same kinds of calculations as it performs with
the LDA potential, but substituting the <i>QSGW</i> self-energy for
the LDA exchange potential.  The
<i>QSGW</i> theory is described in some detail in
Phys. Rev. B<b>76</b>, 165106 (2007), including details of how
<b>lmf</b> is modified to use the <i>QSGW</i> self-energy, and how it
coordinates with the <i>GW</i> package.  Some <A
href="#section3">modifications to the input file</A></B> are required
for <b>lmf</b> to read the <i>QSGW</i> self-energy.

<H2><A name="section1"><FONT size="+1">
1.  lmfgw
</FONT></A></H2>

<P>

The GW driver package contains an executable <B>lmfgw</B>, which acts
to set up input needed by the GW package.  Apart from structural
information, it supplies eigenfunctions and eigenvalues, matrix
elements of the LDA potential, etc..  <B>lmfgw</B> works either with
standard 1-shot <i>GW</i> or with <i>QSGW</i>.

<P>

There is also a feedback to the FP program <B><A
href="fp.html">lmf</A></B>.  In the <i>QSGW</i> mode, a (static)
self-energy <font face="Symbol"><font size="-2">S</font></font> is generated.  When placed
into file <FONT size="+1"><tt>sigm.ext</tt></font> the full-potential
program <B>lmf</B> and the GW driver <B>lmfgw</B> read this potential
and add it as an additional (static) potential to the LDA potential.
(In practice,
<font face="Symbol">S - </font> <i>V</i><sub>xc</sub><sup>LDA</sup> is
stored in <FONT size="+1"><tt>sigm.ext</tt></font>.
See <A href="#section1.4">Reading in the self-energy <font face="Symbol"><font size="-2">S</font></font>
</A>.)

<P>

Apart from these linkages, the <i>GW</i> package is completely
separate from the FP-LMTO program.  The <i>GW</i> package does not
depend on the input files <B>lmf</B> and <B>lmfgw</B> use (e.g. <FONT
size="+1"><tt>ctrl.ext</tt></font> and <FONT
size="+1"><tt>rst.ext</tt></font>).

The <i>GW</i> requires various kinds of input, most all of which 
is supplied automatically by <b>lmfgw</b>.  It has one input file,
<FONT size="+1"><tt>GWinput</tt></font>, which the user must create.
<b>lmfgw</b> can generate a default <FONT
size="+1"><tt>GWinput</tt></font>, which is very useful, since  <FONT
size="+1"><tt>GWinput</tt></font> is rather complicated.  By
invoking <B>lmfgw</B> in a special mode (job -1), it will generate a
default file, which you can then edit.  <b>lmfgw</b> will read an an (optional) category
<FONT size="+1"><tt>GW</tt></font> in <FONT
size="+1"><tt>ctrl.ext</tt></font> to modify the defaults that affect
the generation of <FONT size="+1"><tt>GWinput</tt></font>.

<P>

<b>lmfgw</b> has some other 'job modes' that are read from standard input.
<UL>
<LI>
`Job 0' creates init mode; creates input files the <i>GW</i> package requires:
<FONT size="+1"><tt>SYMOPS, LATTC, CLASS, NLAindx, ldima</tt></font>
<LI>
`Job 1' generates eigenfunctions, eigenvalues, and wave function information
required by the <i>GW</i> package.
</LI> </UL>

These steps must be peformed in a particular order to work with the
<i>GW</i> package.  For details, read the 
<B><A href="man033_ver1.pdf#"><i>GW</i> manual</A></B>.
Alternatively, scripts
<B><A href="#section2">lmgw</B></A> and <B><A
href="#section2">lmgw1</B></A> will carry out 1-shot <i>GW</i> calculations
and automatically handle linkages between the two packages.
These scripts are described in a
<B><A href="#section2">Section 2</B></A>.  Another script <A
href="#section3"><b>lmgwsc</b></A> carries out <i>QSGW</i>
calculations.  It can be substituted for <FONT
size="+1"><tt>gwsc1shot</tt></font> and <FONT
size="+1"><tt>gwsc</tt></font> supplied in the <i>GW</i> package.

<p>

Finally this package contains a second program, <b>lmf2gw</b>, which translates output
of <b>lmf</b> into a format the <i>GW</i> package can read.  It should be compiled automatically
with the <FONT size="+1"><tt>--enable-gw</tt></font> switch when you invoke <b>configure</b>.

<P>

<H2><A name="section1.1"><FONT size="+0">
Category <b><FONT size="+2"><tt>GW</tt></font></b> and its tokens
</FONT></A></H2>

<P>

As noted, this category only affects <b>lmfgw</b> when it creates a default <FONT
size="+1"><tt>GWinput</tt></font> for the user (job -1).  Invoke <b>lmfgw</b> as usual, 
e.g. 
<pre>
  <b>lmfgw</b> ext
</pre>
The standard input will prompt you for a job.  Enter -1.
<b>lmfgw</b> will create a file <FONT size="+1"><tt>GWinput</tt></font>. 

<p>

Tokens in the following 
<b><FONT size="+2"><tt>GW</tt></font></b> category in your control file will affect the contents of
<FONT size="+1"><tt>GWinput</tt></font>.

<pre>
  GCUTB=   plane-wave cutoff for basis envelope functions

  GCUTX=   Plane-wave cutoff for interstitial response function

  QOFFP=   k-point offset parameter for special BZ integration for xi
           See man-gw.tex for description

  NKABC=   No. <I>k</I>-points along each of 3 lattice vectors

  NBAND=   No. bands to include in calculating response function

  NIME=    Number of imaginary energy points for energy contour

  DELRE=   Energy mesh spacing on real axis

  DELTA=   Imaginary energy added to denominator in energy integration of x0.
           Delta should be negative

  DELTAW=  Width for finite difference in numerical differentiaion of
           self-energy, needed to compute the Z factor.

  PBTOL=   retain product basis with overlap eigenvalues>PBTOL
</pre>

<H2><A name="section2"><FONT size="+0">
2. One-shot <i>GW</i> calculations
</FONT></A></H2>

The standard distribution has a script <FONT
size="+1"><tt>gw_lmfh</tt></font> to do one-shot calculations.  It is
described in the <i>GW</i> manual, e.g. man033_ver1.pdf.

This script is very simple, and its functionality is limited.  An
alternative, there are two scripts <b>lmgw</b> and <b>lmgw1</b> that
facilitate calculations and the subsequent extraction of data.

<p>

<FONT size="+1"><tt>lmgw</tt></font>
is a multipurpose script that exercises the <i>GW</i> codes in
many kinds of modes; in particular it can be used in one-shot calculations.
Invoke lmgw without arguments to see what options are available.
The following is a typical invocation:
<pre>
   lmgw --insul=19 --ht srtio3
</pre>

It calculates the QP energies as a 1-shot perturbation given an LDA
pair of input files <FONT size="+1"><tt>ctrl.ext</tt></font> and
<FONT size="+1"><tt>rst.ext</tt></font>,
and file <FONT size="+1"><tt>GWinput</tt></font>.  There is an optional file 
<FONT size="+1"><tt>`switches-for-lm,'  </tt></font> that 
<b>lmgw</b> will use, if it is supplied.  File
<FONT size="+1"><tt>switches-for-lm</tt></font> should contain only one line.
Its contents are passed to <b>lmf</b> and <b>lmfgw</b> as command-line arguments.

<P>

<b>lmgw</b> will carry out all the steps required for a one-shot
calculation.  The main results are kept in files <FONT
size="+1"><tt>TOTE2.UP</tt></font> (and <FONT
size="+1"><tt>TOTE2.DN</tt></font> in spin-polarized cases) and in a
more readable format in files <FONT size="+1"><tt>QPU</tt></font> (and
<FONT size="+1"><tt>QPD</tt></font> for spin-polarized calculations).

<P>

<b>lmgw1</b> is a higher level script that calls <b>lmgw</b>.  It
greatly facilitates the exercising of <b>lmgw</b> by copying files in
and out of directories to enable batch runs.  <b>lmgw1</b> calls extra
scripts:
<UL>
<LI>
<b>infgw</b> extracts timing and other data from output files, and
extracts QP levels into an easy-to-extract format.  Using this
script, <b>lmgw1</b> extract QP data and stores it into file dat.
<LI>
<b>savegwfiles</b> copies selected files in and out of subdirectories.
<LI>
<b>gw-extract-prodbas-and-time-from-output</b> need not be installed.
If it is, it prints out out the size of the LMTO and product bases,
and execution times.
<LI>
<b>xqp</b> is not normally invoked by <b>lmgw1</b>, so you need not
have it installed.  <b>xqp</b> is very useful post-processing script,
however: it formats the results in an easy-to-read manner.  It has
several modes, e.g. it can compare the differences in QP levels from
data in two directories.  For <b>xqp</b> to work, you must have
installed the matrix calculator program
<b>mc</b>. 
</LI> </UL>

<FONT size="+1"><i>Example</i></FONT>.  Assess the effect of floating
orbitals on the QP states in SRTiO<sub>3</sub>.

<P>
Prepare two directories, <FONT size="+1"><tt>NOFLOAT</tt></font> and
<FONT size="+1"><tt>FLOAT</tt></font>.  Each directory should contain
its respective <FONT size="+1"><tt>GWinput'</tt></font>, <FONT
size="+1"><tt>switches-for-lm'</tt></font>, (or distinct files <FONT
size="+1"><tt>ctrl.srtio3</tt></font>) and distinct <FONT
size="+1"><tt>rst.srtio3</tt></font>).  The two directories serve as
distinct input streams, in this case to two different LMTO basis sets:
one containing floating orbitals and the other not.  (Obviously you
can make other kinds of variations such as vary entries in <FONT
size="+1"><tt>GWinput'</tt></font>).

<P> Do the 1-shot calculations:

<pre>
   lmgw1 --autoht --infgwx --insul=19 -job= FLOAT srtio3
   lmgw1 --autoht --infgwx --insul=19 -job= NOFLOAT srtio3
</pre>

<FONT size="+1"><tt>--insul=19 </tt></font> tells <b>lmgw</b> to
define the zero at the 19th band.  (The <i>GW</i> code always uses a
particular level at the <font face="Symbol">G</font> point,
i.e. <b>q</b>=0, in this case the level 19, for the reference).

<p>
Invoke <b>lmgw1</b> with no arguments to see what the rest of the switches mean.

<p>
<b>lmgw1</b> puts new files into directories
<FONT size="+1"><tt>NOFLOAT</tt></font> and
<FONT size="+1"><tt>FLOAT</tt></font>, e.g files:
<pre>
   ctrl.preprocessed.srtio3 GWinput switches-for-lm TOTE.UP dat job QPU TOTE2.UP version rst.srtio3
</pre>

Use <b>xqp</b> to retrieve data in a pretty format.  This command
generates QP data for levels 1,2,6,9,12,15,18,20 in 3 groups of 3:
<pre>
   xqp -lda -lst=1,2,6,9,12,15,18,20 -eshft=-.36,-.36,-.35 -qp=0.00,0.00,0.00 -qp=0.50,0.50,0.50 -qp=0.50,0.00,0.00 FLOAT
</pre>
Output:
<pre>
    qp=0.00,0.00,0.00      qp=0.50,0.50,0.50   qp=0.50,0.00,0.00 
    eQP  eQP(noZ) eLDA    eQP  eQP(noZ) eLDA    eQP eQP(noZ)  eLDA
  -32.70 -32.87 -32.49  -33.43 -34.46 -32.43  -33.14 -33.64 -32.57
  -33.43 -34.49 -32.36  -33.43 -34.46 -32.43  -33.42 -34.48 -32.36
  -16.99 -17.66 -16.13  -17.04 -17.66 -15.74  -16.97 -17.63 -16.09
  -15.65 -15.91 -14.53  -15.99 -16.24 -14.82  -15.59 -15.85 -14.48
   -3.00  -3.09  -2.85   -3.97  -3.87  -4.21   -3.01  -3.08  -2.89
   -1.25  -1.27  -1.21   -3.65  -3.59  -3.79   -2.30  -2.25  -2.41
   -0.36  -0.36  -0.35    0.00   0.00   0.00   -0.77  -0.79  -0.73
    3.81   4.53   1.77    6.41   7.39   4.48    4.02   4.77   1.95
</pre>

The <FONT size="+1"><tt>-eshft</tt></font> switch was necessary here
because the QPU generator shifts the energy zero around by assuming
the VBM is at Gamma; for SrTiO<sub>3</sub> the VBM falls at
(1/2,1/2,1/2).

To see the difference in the two QP calculations, invoke <b>xqp</b> with two
directories as:

<pre>
   xqp -lda -lst=1,2,6,9,12,15,18,20 -qp=0.00,0.00,0.00 -qp=0.50,0.50,0.50 -qp=0.50,0.00,0.00 FLOAT NOFLOAT
</pre>

Output:
<pre>
    qp=0.00,0.00,0.00      qp=0.50,0.50,0.50   qp=0.50,0.00,0.00 
    eQP  eQP(noZ) eLDA    eQP  eQP(noZ) eLDA    eQP eQP(noZ)  eLDA
    0.09   0.13   0.01    0.06   0.09   0.01    0.08   0.11   0.01
    0.06   0.08   0.01    0.06   0.09   0.01    0.06   0.08   0.00
    0.06   0.06   0.00    0.06   0.07   0.00    0.05   0.07   0.00
    0.05   0.07   0.01    0.05   0.07   0.01    0.05   0.07   0.01
    0.01   0.01   0.00    0.01   0.00   0.01    0.01   0.00   0.00
   -0.02  -0.02  -0.01    0.01   0.00   0.01    0.00  -0.01   0.01
    0.00   0.00   0.00   -0.01  -0.02   0.00    0.00  -0.01   0.00
   -0.03  -0.03  -0.01   -0.02  -0.03   0.00   -0.02  -0.04   0.00
</pre>

In the SrTiO<sub>3</sub> case, floating orbitals don't matter much.

<p>

Run <b>xqp</b> with no arguments to see other options.  


<H2><A name="section3"><FONT size="+0">
3. Quasiparticle self-consistent <i>GW</i> calculations
</FONT></A></H2>

To performed <i>QSGW</i> calculations, simple shell scripts
<b>gwsc1shot</b> and <b>gwsc</b> supplied in the <i>GW</i> package,
and how to use them is documented in the <B><A
href="man033_ver1.pdf#"><i>GW</i> manual</A></B>.

Alternatively, a high-level script <b>lmgwsc</b> can be used.  It
calls script <A href="#section2"><B>lmgw1</B></A> for a single
iteration, then performs some checks, and calls <b>lmgw1</b> again,
cycling through iterations.  <b>lmgwsc</b> is useful because it has
many more options than <b>gwsc</b>. For example it monitors the RMS
change in <font face="Symbol"><font size="-2">S</font></font> and will stop iterating when a
certain tolerance is met. <b>lmgwsc</b> `iteration 0' corresponds to
<b>gwsc1shot</b>, and 
<b>lmgwsc</b> corresponds to <b>gwsc</b> in subsequent iterations.
Invoke <b>lmgwsc</b> without arguments to see the various options.

<H2><A name="section3.1"><FONT size="+0">
Additions to lmf to read <i>QSGW</i> self-energy <font face="Symbol"><font size="-2">S</font></font>.
</FONT></A></H2>


To read the the self-energy from file <FONT
size="+1"><tt>sigm.ext</tt></font>, you must add a token to category
<b><FONT size="+2"><tt>HAM</tt></font></b>.

<pre>
  RDSIG=   token specifying whether to read the self-energy generated
           by the GW package and add to to the LDA hamiltonian.
           The argument of RDSIG= consists of three digits:
           1s digit 1  : read sigma but to not symmetrize
                    2  : read sigma and symmetrize it

           10s digit interpolation mode; see Notes on interpolation of sigma
                    0  : Simple bloch sum of sigma
                    1  : (strongly recommended) approx high-energy states
                         with diagonal part of sigma.
                    2  : perturbation approach (not documented)
                    3  : Linear interpolation of sigma from neighboring points
                         (not documented)

          100s digit is specific to 10s digit = 2.  It specifies how many
               neighboring points should be included in the interpolation.
               Choosing 0 causes the program to choose the default value of 4.

        10000s digit has to do with symmetry operations and how
               lattice vectors are assembled.

               1 flags lmf that file <B>sigm</B>.ext contains no
                 symmetry operations.  Equivalent to --rsig:fbz.

               2 tells <B>lmf</B> not to augment the standard nk1*nk2*nk3
                 lattice translation vectors that would follow from
                 standard inverse fast fourier transform.  By default
                 some extra vectors are added to make the system
                 compatible with symmetry operations.  If you use this
                 option, you cannot symmetrize sigma.

               4 warning on file <I>k</I>-point mismatch.  When reading file
                 <B>sigm</B>.ext, the reader checks that the file <I>k</I>-point
                 is equivalent to the current <I>k</I>-point in the list of
                 irreducible <I>k</I>-point it has constructed.  If there is a
                 mismatch, the default is to abort.  Use this switch
                 when, e.g. you deform a lattice slightly and want to
                 use estimate sigma from the sigma for the original
                 lattice.
                 
                 Linear combinations of these switches are allowed.

  RSRNGE=  token specifying the maximum range of the r.s. sigma (units of alat)
           computed by fourier transform of sigma(k) generated by the GW package.

</pre>

<H2><A name="section3.2"><FONT size="+0">
Reading in the self-energy <font face="Symbol"><font size="-2">S</font></font>
</FONT></A></H2>

The self-energy <font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) (actually <font
face="Symbol">S - </font> <i>V</i><sub>xc</sub><sup>LDA</sup>) is
generated by the <i>GW</i> package on a mesh of points in the
irreducible Brillouin zone (BZ).  When this information is stored in
file <FONT size="+1"><tt>sigm.ext</tt></font>, in can be read by
<b>lmf</b> and
<b>lmfgw</b>, and added to the LDA one-particle potential.  When
reading this file, <font face="Symbol"><font size="-2">S</font></font> is rotated to the
entire BZ and then inverse Fourier (Bloch) transformed into a
real-space form <font face="Symbol"><font size="-2">S</font></font>(<B>T</B>).  Thus, <font
face="Symbol">S</font>(<B>T</B>) can be computed at any <b>k</b>-point, so 
<b>lmf</b> can work in the just the same way in using a <i>QSGW</i>
potential as it does an LDA potential.  This is a crucial advantage of
the present implementation.   The output from <b>lmf</b> indicating that
the Bloch transformed has been effected looks like the following:
<pre>
   q-points in full BZ where sigma calculable ...
   BZMESH:  18 irreducible QP from 18 ( 3 3 2 )  shift= F F F
   Irr. qp for which sigma is calculated ...
   BZMESH:  6 irreducible QP from 18 ( 3 3 2 )  shift= F F F
   hft2rs created hrs:  ndhrs=29  max Re(hrs) = 1.21  max Im(hrs) = 6e-4
</pre>


The correctness of this inverse FT is checked by computing the Forward
FT of <font face="Symbol"><font size="-2">S</font></font>(<B>T</B>) at the known
<B>k</B>-points; and the difference in the two is printed in the lines
<pre>
   comparing s(iq) to interpolated sq for iq=   1  errmx= 5.6E-17 at i,j=  30  30
   comparing s(iq) to interpolated sq for iq=   2  errmx= 6.2E-17 at i,j=  34  34
   ...
</pre>

The error should be approximately that of machine precision.  Next, if
you so specify (see token <FONT size="+1"><tt>RDSIG</tt></font>),
<B>lmf</B> will symmetrize <font face="Symbol"><font size="-2">S</font></font>(<B>T</B>)
according to the given group operations; the difference in the file
<font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) and FT of the symmetrized <font
face="Symbol">S</font>(<B>T</B>) is checked once again.  Because of
small numerical errors in the generation of <font
face="Symbol">S</font>, the symmetrized <font face="Symbol"><font size="-2">S</font></font>
will differ a little from the unsymmetrized one.  <FONT
size="+1"><tt>"maximum error"</tt></font> in the following output
indicates the effect of symmetrization:
<pre>
   rsmsym: symmetrizing complex s(1..1704) using 8 group operations
   symstr: max asymmetry = 7.93e-7
   check FT s(R) against s(q) ... maximum error = 6.2e-5
   check FT s(R) against s(q) ... maximum error = 6.2e-5 spin 2
</pre>


<P>

<H2><A name="section3.3"><FONT size="+0">
Notes on interpolation of sigma
</FONT></A></H2>

It turns out that it is difficult to interpolate sigma to other
<B>k</B>-points.  The difficulty apparently rises from some numerical
instabilities the diagonal part of sigma (in the basis of LDA
eigenvectors) for higher energies, typically above 2 Ry.  A couple of
procedures have been developed to circumvent this problem.

<pre>
mode 1: Above a certain cutoff specified by the user, the off-diagonal
        part of sigma in the representation of LDA eigenstates is
        eliminated, and also the the diagonal part of sigma (sigii) is
        constrained for these states.  That is, starting
        representation is computed from the LDA hamiltonian, with the
        states ordered in ascending order of the LDA eigenvalues.  In
        this representation, the self-energy for states above a cutoff
        (you choose the cutoff either as those states above a
        specified LDA energy, or above a specied number of states) is
        taken to be diagonal, and the diagonal part constrained.

        There are a couple of options for the constraints (see mode
        below).  The default option is to contrain sigii to lie
        between a minimimum and maximum value.  You can also constrain
        to be either larger than, or equal to, a specified linear
        function of the LDA energy (asig+bsig*Elda) (NB: it turns out
        that sigma is reasonably linear in Elda).  In any case, these
        high-lying states are far enough away from the Fermi level
        that their effect should be small, and the result should
        depend very little on the choice of constraint. Another point
        to keep in mind is that approximating sigma for these states,
        you ensure that the LDA and quasiparticle eigenvectors for
        those states are the same.

        Finally, you can also elect to eliminate the off-diagonal
        parts of sigma for the lowest-lying states.  In this case,
        there is no option to constrain sigii.

        When using this mode, you specify the cutoffs and constraints
        using the SIGP: token.  The syntax of this token (category HAM) is
          SIGP:modsgp,nmin,emin,nmax,emax,asig,bsig,efit

        modsgp selects the kind of constraint:
              0 constrain sigii to be > asig+bsig*Elda
              1 constrain sigii to be = asig+bsig*Elda
              2 constrain sigii to be  asig < sigii < bsig
              3 constrain sigii as in modsgp 1.  The difference between
                modes 1 and 3 are merely informational; see below.
                default: 3

        nmin  is an integer specifying how many of the lowest-lying
              states are approximate sigma by discarding the
              off-diagonal parts in the basis of LDA functions.
              If nmin is zero, no low-lying states are approximated.
                default: 0

        emin  is an alternative way to specify approximation of low-
              lying sigma.  It is only used if nmin<0, which which case
              emin is an energy cutoff: states below emin are approximated.
                default: not used

        nmax  is an integer specifying which of the highest (LDA)
              states are to be approximated.  States higher than nmax
              have the off-diagonal part of sigma stripped; moreover
              (unlike the low-lying states) the diagonal part of sigma
              is constrained (see modsgp above).
              If nmax is <= 0, it is not used; see emax below.
                default: 0

        emax  is an alternative way to specify approximation of high-
              lying sigma.  It is only used if nmax<=0, which which case
              emax is an energy cutoff: states above emax are approximated.
                default: 2

        asig  is the first coefficient in the constraint.
              In the linear constraints (modes 0,1) it is the constant
              coefficient; for modsgp=2, it is the lower bound.
                default: 0.02 (good estimate for Si)

        bsig  is the second coefficient in the constraint.
              In the linear constraints (modes 0,1) it is the linear
              coefficient; for modsgp=2, it is the upper bound.
                default: 0.06 (good estimate for Si)

        efit  lower bound for least squares fit; see next section.

       *Estimation of asig and bsig:

        lmf will supply you information to find reasonable estimates for
        asig and bsig.  For modsgp<3, it will make a least-squares fit
        to sigii for states higher than efit. For modsgp=3, it will make
        a least-squares fit for states between efit and emax.  You must
        use the latter if you are going to only going to calculate sigma
        for states above some emax (set as in the last lines in the
        GWIN0 file).  For the former, you must invoke lmf ON THE MESH of
        k-points for which sigma is known (there appear to be fewer
        problems with interpolation on that mesh).  lmf accumulates the
        minimum, maximum, and least-squares fit for the as-given sigma
        for all the states above the cutoff.  (Look in the output for a
        line beginning with "hambls:".  Also, if you set the verbosity
        above 45, lmf will print out the calculated sigii for each of
        these states, together with the constrained value.)  For either
        style, lmf will write to file sigii.ext the data used to make
        the fit, and summarize the fit and the end of the file.

mode 2: is an experimental perturbation approach that is not documented.
        mode 1 works well enough that it is probably not needed.

mode 3: linearly interpolates sigma in the basis of the LDA wave
        functions, by finding k-points in the vicinity of the desired
        k-point and performing a linear interpolation.  The number of
        neighbors used in the interpolation can be 1,2,3, or 4.  The
        difficulty with this approach is that the interpolation is not
        smooth.  The experience to date suggests that mode 1 is the
        most reliable.  The default number of k-points used in the
        interpolation is 4.  You can specify a different number
        with the 100s digit of token <A href="#section3.1"><B><FONT
        size="+1"><tt>RDSIG</tt></font></B></A>).
</pre>


<H2><A name="section3.4"><FONT size="+0">
Notes on generating m-resolved DOS.
</FONT></A></H2>

In its current form, lmf does not allow you to generate m-resolved dos
with symmetry operations turned on, because they are not properly
symmetrized.  To make m-resolved dos, you must do the following:

<pre>
1.  First create a self-energy file (sigm.ext) without symmetry.  You
    can do this by invoking <A href="#section3.6"><B>lmf --wsig:fbz ...</B></A>
    this will cause lmf to write a new sigma file sigm2.ext with the
    same structure as sigm.ext (see below), but generated for points
    in entire BZ, which of course corresponds to the irreducible set
    with no symmetry operations.

2.  Turn off symmetry operations in your ctrl file, e.g. add a line
    SYMGRP E

3.  run lmf with
    --<B>pdos:mode=2:</B>...
    (see lmto.doc for a description of the options available to --<B>pdos</B>)

4.  run lmdos with the same --<B>pdos</B> switch as you ran lmf.
    File dos.ext will be created, containing DOS resolved by m .

</pre>

<H2><A name="section3.5"><FONT size="+0">
Taking a linear combination of two sigma files.
</FONT></A></H2>

You can read two different self-energies sigma from distinct files and
combine them in a linear combination.  Usually sigma is read from file
<B>sigm</B>.ext.  You can read sigma as some linear combination of <B>sigm</B>.ext
and <B>sigm1</B>.ext as
<pre>
    sigma = alpha*(sigm.ext) + beta*(sigm1.ext)
</pre>
by invoking lmf with
<pre>
   --mixsig=alpha,beta
</pre>
If second parameter is missing, the effect is to read in <B>sigm</B>.ext scaled
by factor alpha.

<H2><A name="section3.6"><FONT size="+0">
Output sigma file and modifications to sigma file.
</FONT></A></H2>

You can generate an output sigma file with the switch --<B>wsig</B>.
A new file, <FONT size="+1"><tt>sigm2.ext</tt></font> written, and 
<B>lmf</B> then exits.

You can use this switch to generate, for example, a properly
symmetrized self-energy (see <A href="#section3.1"><B><FONT
size="+1"><tt>RDSIG</tt></font></B></A>).

<P>

The full syntax is, with available options:
<pre>
  --<B>rsig[:fbz][shftq[=#,#,#]]</B>
</pre>

<pre>
 :fbz   flags lmf that file <B>sigm</B>.ext contains no symmetry operations.  
        Equivalent to setting 10000s digit=1 in token RDSIG= .

 :shftq=#,#,# tells the program that <font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) stored in file 'sigm'
        was generated on a kmesh offset by #,#,#

 :shftq without arguments performs the same function.   The offset is taken from a default.

</pre>


Additionally you can make several modifications to the sigma file
before writing, though switches passed with the --<B>wsig</B> argument.

The full syntax is, with available options:
<pre>
  --<B>wsig[:fbz][:newkp][:onesp][shftq[=#,#,#]][:rot=rotation-string][:edit][:phase][:trans=#]</B>
</pre>
The delimiter `:' separating options can actually be any ascii character,
e.g. `^'.  (Sometimes you can't use `:' as the delimiter because it
may used in one of the arguments, and you will confuse the parser.)
The options do the following

<pre>
 :fbz   causes <B>lmf</B> to ignore symmetry operations and generate a sigma
        file for k-points in the entire BZ.  It is useful when generating a
        sigma file that allows fewer symmetry operations, e.g. when making the
        <A href="#section3.3"><I>m</I>-resolved density-of-states</A>, or creating a <A href="#section3.6">trial sigm for a sheared crystal</A>.
        Note: to <I>read</I> a sigm file generated by --wsig:fbz, you must set 
        the 10000's digit in token <A href="#section3.6">RDSIG</A>.

 :newkp causes <B>lmf</B> to use the mesh used by the lmf program.  By default
        <B>lmf</B> generates sigm on the mesh as defined in the sigm file.
        This switch enables the user to interpolate the sigma to a new <B>k</B> mesh.

 :onesp (for spin-polarized case) causes <B>lmf</B> to average the up- and down-
        self-energies, and write a nonmagnetic self-energy to disk.

 :shftq=#,#,# is used in conjunction with trans=#, where #=1,2,3; see below.
        It flags that an offset of #,#,# should be added to the standard <i>k</i> mesh
        when generating transformed <font face="Symbol"><font size="-2">S</font></font>.
        to avoid inadvertant swap of (nearly) degenerate LDA eigenvalues as the LDA basis
        set is changed.  (The scheme requires that eigenvalues in the two basis sets
        synchronize in the energy range of interest).

 :rot=rotation-string causes <B>lmf</B> to rotate sigm to a coordinate
        system, whose rotation is specified by 'rotation-string'.
        Here is an example :
             lmf si --wsig^rot=z:.1,y:2,z:.3
             cp sigm2.si sigm.si
             ... Edit the STRUC category to include ROT=z:.1,y:2,z:.3 and then do:
             lmf si --rs=101

        Points to note:

        1. Switch --rs=101 tells lmf to rotate the local charge density read from
        the restart file.

        2. Note the use of caret (^) as delimiter, since (:) is needed in the string

        3. <B><A href="rotations.html">rotations.html</A></B> describes how to specify a rotation.

        4. Be advised that the interpolation may not proceed in exactly
        the same way in the rotated and unrotated systems, owing to different
        sets of lattice vectors taken for large r.  To minimize these differences,
        work with a symmetrized self-energy as much as possible, and
        read the density (and lattice positions) from the same restart file.

        The example below follows the recommended procedure. The STRUC
	category of input file <B>ctrl</B>.si has this segment:
          %ifdef rot
        	  ROT=z:.1,y:2,z:.3
          %endif
        which applies a rotation when variable <I>rot</I> is nonzero.
        Also the HAM category has token <B>RDSIG</B>=<I>sig</I>. 

          ... first step symmetrizes sigma
          lmf si -vsig=12 -vrot=0 --wsig
          cp sigm2.si sigm.si

          ... this pass in the unrotated coordinated system
          lmf si -vsig=12 -vrot=0 

          ... should match this pass in the rotated coordinated system
          lmf si -vsig=12 --rs=101,0 -vrot=0 --wsig^rot=z:.1,y:2,z:.3
          cp sigm2.si sigm.si
          lmf si -vsig=12 --rs=101,0 -vrot=1 

 :trans=# is designed to enable transformations of <font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) 
        to new basis sets.

        trans=1 causes <B>lmf</B> to read <font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) (file sigm) in its usual form (orbital basis)
        and transform it to the LDA basis. <font face="Symbol"><font size="-2">S</font></font><sup>LDA</sup>(<B>k</B>) is written to file sigm2.

        trans=2 is similar, but the high-energy part of <font face="Symbol"><font size="-2">S</font></font><sup>LDA</sup>(<B>k</B>) is replaced before
        it is written, as described <A href="#section3.3">here</A>.

        trans=5  also substitutes the high-energy part of <font face="Symbol"><font size="-2">S</font></font><sup>LDA</sup>(<B>k</B>), but stores
        <font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) in the orbital basis.

        trans=3 causes <B>lmf</B> to generate and LDA eigenfunctions <i>Z</i><sup>LDA</sup>(<B>k</B>) and store them
        in file 'evec.'  <font face="Symbol"><font size="-2">S</font></font>(<B>k</B>) is not used in this mode; file sigm is not touched. 
        Note: you must run trans=1 first because it generates a list of k-points.
        
        trans=-1 reads <font face="Symbol"><font size="-2">S</font></font><sup>LDA</sup>(<B>k</B>) from file 'sigm2' and <i>Z</i><sup>LDA</sup>(<B>k</B>) from file 'evec'
        and generates <font face="Symbol"><font size="-2">S</font></font><sup>orb</sup>(<B>k</B>), i.e. <font face="Symbol"><fontsize="-2">S</font></font> in the usual orbital basis.
	The transformed file is written into file 'sigm3'.

        The following sequence should: make sigma(LDA) and evec, and then transform sigma(LDA)->sigm(orb)
        It should recover the original sigm file (apart from numerical imprecision).

        doing nothing (except for roundoff errors)
          <b>lmf</b> ext --wsig:shftq:fbz:trans=1
          <b>lmf</b> ext --wsig:shftq:fbz:trans=3
          <b>lmf</b> ext --wsig:trans=-1
          cp sigm3.ext sigm.ext
        ... run with:
          <b>lmf</b> ext --rsig:fbz:shftq 

        However, you can change the orbital basis between lmf invocations, thus
        transforming sigm to a new orbital basis.
        The following is the recommended way to remake sigm in a new basis
          lmf ext -vsig=12 --rs=1,0 -vnit=1 -vemaxs=2 --pr45 --wsig:trans=5
          cp sigm2.ext sigm.ext
          cp sigm.ext sigm.ext~
          <b>lmf</b> ext --wsig:shftq:fbz:trans=1
        ... change basis set in input file         
          <b>lmf</b> ext --wsig:shftq:fbz:trans=3
          <b>lmf</b> ext --wsig:trans=-1
          cp sigm3.ext sigm.ext
        ... run with:
          <b>lmf</b> ext --rsig:fbz:shftq         

         For an example, after running test 
           gw/test/test.gw si
         Do:
           cp gw/test/si/sigm.si .
           rdcmd '-f:#rdcmd:%2f' -cat:TSIGT --noerr ctrl.si > out.trans.si

 :edit  invokes an interactive editor that allows you to modify
        elements in the sigma file

 :phase causes <B>lmf</B> to scale rows and columns associated with particular
        sites in the sigm matrix by a <B>q</B>-dependent phase.  This is
        because in the standard FP program <B>lmf</B>, the structure matrix
        is computed by Ewald sums which contain extra site-dependent <B>R_j</B>
        and <B>q</B>-dependent phases 
             exp(i <B>q.R_j</B>)
        in when summing over sites <B>R_j</B> in the lattice.  The
        coordinates R_j are internally shifted automatically by some
        combination of lattice vectors, to shorten their length and
        avoid convergence problems in the Ewald sum.

        For that reason, when using the sigm file generated by
        rot=... , the user is strongly advised to construct the
        input file and its basis vectors in the unrotated coordinate
        system, and rotate the lattice and basis vectors through the
        <B>ROT=</B> token in the <A href="lmto.html#STRUC"><B>STRUC</B></A> category.
        (The syntax for this string is identical to :rot= above.)
        Additionally, the site positions should be read through the
        rst file, as they will not be further shortened by <B>lmf</B>.

        Use option :phase if for some reason you do change the
        shortening of lattice vectors.  For this option to work, the
        translation vectors must be supplied in a file
        <B>shorps</B>.ext.   <B>lmf</B> reads this file and scales
        sigm accordingly.

        You can create file <B>shorps</B>.ext by hand, or do it automatically.
        Either way, you must coordinate the positions <B>lmf</B>
        actually uses with the sigm file.

        To generate file <B>shorps</B>.ext automatically, you can invoke <B>lmf</B>
        with the switch --<B>shorps</B>, which, after reading site positions from
        file <B>rst</B>.ext, further shortens the site positions and saves
        the lattice translations used for shortening in <B>shorps</B>.ext.
        Take care to rewrite file <B>rst</B>.ext, e.g. by invoking
        <B>lmf --shorps ...</B> with no iterations.

        Once <B>shorps</B>.ext is created, you can generate a new sigma:
          <B>lmf --wsig:phase ...</B>
          <b>cp sigm2.ext sigm.ext</b>

        With both files <B>rst</B>.ext and <B>sigm</B>.ext altered,
        <B>lmf</B> should generate identically the same eigenvalues
        and output density.

</pre>

<H2><A name="section3.6"><FONT size="+0">
Example: adapting a sigma file to a symmetry-lowering shear distortion
</FONT></A></H2>

Two changes have to be made to the sigma file <FONT size="+1"><tt>sigm.ext</tt></font>.
<P>
1. Because the starting <FONT size="+1"><tt>sigm.ext</tt></font> has higher symmetry than the case
of interest, the new <FONT size="+1"><tt>sigm.ext</tt></font> may require more irreducible <I>k</I>-points.
<P>
2. The
sigma file reader will by default require the <I>k</I> points in the
file correspond to the <I>k</I> points generated by <B>lmf</B>.  This will no
longer be the case if you want to adapt <FONT size="+1"><tt>sigm.ext</tt></font> to a distorted mesh.

<P>

Supposing for definiteness that your input file has token <FONT
size="+1"><tt>RDSIG</tt></font> set by variable sig and that you want
to make a shear distortion along the [001], which you control by
variable tet:
<pre>
  HAM  RDSIG=sig
  STRUC SHEAR=0 0 1 tet
</pre>

<P>

To solve problem 1 above, generate a new <FONT size="+1"><tt>sigm.ext</tt></font> on a mesh without any symmetry operations:
<pre>
  lmf ... original arguments ... --wsig:fbz
  cp sigm2.ext sigm.ext
</pre>

<P>

Now you are in a position create a useable <FONT size="+1"><tt>sigm.ext</tt></font>.
Change your input file or alter command-line switches for the symmetry-changing distortion.  
You also must tell the sigma reader to (1) that the input <FONT size="+1"><tt>sigm.ext</tt></font> has no symmetry operations and 
to overlook the mismatch in <I>k</I>-points.  To this end, use token RDSIG=50011.
In the above example, you would not change the input file but invoke:
<pre>
  lmf ... original arguments ... -vsig=50011 -vtet=1.01 --wsig:newkp
  cp sigm2.ext sigm.ext
</pre>

<P>

This will create a <FONT size="+1"><tt>sigm.ext</tt></font> with the <I>k</I>-points appropriate to the distorted lattice.

<H2><A name="section3.7"><FONT size="+0">
Example: illustration of the sigma editor.
</FONT></A></H2>

In this example, we consider GdN in the zincblende structure, and zero
out all of the sigma except the site-diagonal <I>f-f</I> block.  This enables
us to compare what contributions are LDA+U-like, and what lie outside
the domain of the usual assumptions in LDA+U.  

<P>

The Gd basis in this example consists of two envelope functions and
and <I>f</I> local orbital.  Thus without some further reductions, it is not
particularly meaningful to say just what U is. 


<P>

Our plan is instead to include the <I>f-f</I> channels in all three kinds of
<I>f</I> blocks; this is presumed to be what a suitably chosen LDA+U would
generate.  The sigma for all other blocks will be zeroed out in this example.

<P>

The sigma editor modifies directly the real space sigma, so it is a
straightforward exercise to mask out all parts of sigma not in one of
the three subblocks, as the example shows.

To invoke the editor, run lmf as:
<pre>
  lmf ... original arguments ...  --wsig:edit
</pre>
You should see this appear on the screen:
<pre>
  Welcome to the sigma file editor.
  Option (enter ? for menu) : 
</pre>

The editor works on the R.S. sigma_RR', in one of two modes.  Either
you edit a range of pairs (RR'), in entire blocks, or you can work on
a particular RR' pair, where you can do things in more detail.  We
have to do both here, because we want to eliminate sigma_RR' unless R
and R' are the (on-site) Gd pair --- zero translation vector.

<P>

The editor starts out in the 'range' mode, with a default range
corresponding to all pairs.  The pairs are enumerated in a long list,
you must refer to the pair you want by a number.  To see what number
is appropriate, you can get an summary of the neighbor table by
invoking <pre> sum </pre> In this example, there are two atoms, Gd and
N.  The summary of the table looks like : 
<pre>
  Summary of pair table for R.S. sigma
  site  start   end  size
    1      1    553  553
    2    554   1106  553
</pre>
Atom 1 is Gd; there are many neighbors (this test was done for 6x6x6 k points).
The first one corresponds to the on-site Gd-Gd pair.  To see this explicitly, type 
<pre>
  show
</pre>
The editor will print out many lines, one for each pair.  A partial list of the table is:
<pre>
  ib  jb         --- r(jb)-r(ib) ---          d       -x-plat-   pair
   1   1   0.000000   0.000000   0.000000   0.0000     0  0  0     1
   1   2  -0.500000   0.000000   0.000000   0.5000     1  0 -1     2
...
   1   1   1.500000  -0.500000   1.000000   1.8708    -1  3  0   250
   1   1  -0.500000   1.500000  -1.000000   1.8708     1 -3  2   251

  ib  jb         --- r(jb)-r(ib) ---          d       -x-plat-   pair
   2   2   0.000000   0.000000   0.000000   0.0000     0  0  0   554
   2   1  -0.500000   0.000000   0.000000   0.5000     0 -1  0   555
   2   1   0.000000   0.500000   0.000000   0.5000     0 -1  1   556
...
</pre>
Pair 1 corresponds to ib=1 and jb=1 (Gd-Gd pair) with 0 distance
between them.  (Pair 554 corresponds to onsite N-N pair.  It is
interesting to see what the effective U is on the N p orbital also.)
So, it's pair 1 we want to keep, and discard all other pairs.

To do this, invoke 
<pre>
 mask 1,1
</pre>
This zeros out hrs for any pair not in the range (1..1), i.e. all pairs
except the first one.

We could quit at this point.  If we did so, we would include
<I>all</I> of the U parameters coupling Gd orbitals to themsevlves:
not only the Gd <I>f-f</I> blocks, but the <I>f-d</I> blocks, the
<I>p-d</I> blocks, and so on.

To restrict the sigma to just the <I>f-f</I> parts, we must specialize
to the first pair.  Tell the editor that you want to specialize to
that pair by invoking
<pre>
  pair 1
</pre>
You can get a summary of the structure of the hamiltonian for that pair by invoking 
<pre>
  sum
</pre>
(Note the 'sum' command takes a different meaning when editing a particular pair).
In this example, the hamiltonian for Gd has this structure:
<pre>
  Summary of hamiltonian subblocks for site   1
      kap  lstart  h(k)    h(rs)     size
                  start    range
       1      0      1     1,25       25
       2      0     26     26,41      16
       3      1     42     42,44       3
       3      3     45     45,51       7
</pre>
The three <I>f</I> blocks start are along the row and column dimensions (10,16), (35,41), (45,51).
To mask out the block except the 9 <I>f-f</I> subblocks, do
<pre>
  mask 10,16,10,16 10,16,35,41 10,16,45,51 35,41,10,16 35,41,35,41 35,41,45,51  45,51,10,16 45,51,35,41 45,51,45,51 
</pre>

You can see what the on-site sigma looks like now, by invoking 
<pre>
  show
</pre>
That command prints out a lot of data, which isn't convenient to look at.
You can tell the editor to display just the subblocks you are interested in:
<pre>
  show 10,16,10,16 10,16,35,41 10,16,45,51 35,41,10,16 35,41,35,41 35,41,45,51  45,51,10,16 45,51,35,41 45,51,45,51
</pre>

Finally, we are ready to quit.  Invoke
<pre>
  q
</pre>
and the editor will exit, saving the file into sigm2.ext.
In the usual way, you can do
<pre>
  cp sigm2.ext sigm.ext
</pre>
and you are ready to invoke lmf with a modified sigma.


<H2><A name="section4.1"><FONT size="+0">
4. Frequently asked questions
</FONT></A></H2>

<BR><I>Q</I>: When I do Ni with 10x10x10 divisions, I get

<pre>
 hft2rs created hrs:  ndhrs=39  max Re(hrs) = 0.771  max Im(hrs) = 2.33e-9

   i   j      diff              bloch sum                file value
  21   1    0.000174     -0.036707    0.000000     -0.036533    0.000000

 Exit -1 rdsigm: Bloch sum deviates more than allowed tolerance (tol=5e-6)
</pre>

What does this mean?  or how to avoid the problem?

<BR>
<BR><I>A:</I>

It almost always means that the range set up for the r.s. Bloch sigma
is not long enough to encompass all possible connecting vectors.  The
easiest way to fix it is to increase the range.  You need to add to
the HAM category something like

<pre>
HAM   RSRNGE=8
</pre>

In the output, look for where the sigma is read in.  You will see
something like this output

<pre>
 RDSIGM: read file sigm and create COMPLEX sigma(R) by FT ...
         Sigm will be approximated by:  diagonal Sigma for high and low states 
         Approximate sigma for energies E(lda)>2.5
         For high states Sigii = 0.000 + 0.080 * E(lda) 
         Linear fit to diagonal sigma for E(lda)>0 and E(lda)<2.5 
         sigm file has 20 irreducible QP: nk = ( 4 4 4 )  shift= F F F
 hft2rs: make neighbor table for r.s. hamiltonian using range = 5 * alat
 pairc:   8278 pairs total 4139 is max cluster size
 hft2rs: found 256 connecting vectors out of 256 possible for FFT
</pre>

<tt>range = 5 * alat</tt> means that RSRNGE=5 (default value).  If RSRNGE is
large enough, this line
<pre>
 hft2rs: found 256 connecting vectors out of 256 possible for FFT
</pre>
indicates that all possible connecting vectors were found (as they are
in this example); then there should be no problem.
