The empirical tight-binding code tbe

Summary

This page documents tbe, a code that evaluate the electronic structure in an empirical tight-binding framework, that is where the hamiltonian matrix elements are given as input.

It has several enhancements to a basis tight-binding scheme. In addition to standard features of tight-binding hamiltonians (Molecular dynamics and statics) it has several novel features. The is a facility for self consistency and for magnetic tight binding. This code also interfaces with Cheriotti’s i-pi, for prosecuting exotic molecular dynamics such as Feynman’s path-integral theory.

Background, and specification of matrix elements
- Example ME category
The TB and START categories
Self consistent, polarisable ion, tight binding
Magnetic tight binding
Molecular dynamics and statics
Band and DOS plotting
tbe-specific command line switches
References

Background, and specification of matrix elements

Initially the tbe program was written to make quick and easy tight-binding (TB) calculations in the spirit of the early work by Jim Chadi [1] and Walter Harrison [2]. Later, after the invention of the polarisable ion tight binding theory [3] the code became much more sophisticated and now includes charge and spin self consistency, molecular statics and dynamics. The theory has now been described in detail in Mike Finnis’s book [4] and in Tony Paxton’s lecture notes [5]. Magnetic tight binding is explained in Paxton and Finnis [6].

The empirical tight binding was built around the ASA package so as to exploit the sparse matrix handling, symmetrisation, BZ integration and so on. As a result, some ASA tokens are sought from the usual input file ctrl.ext, even if they are not used. An example is the P token in START.

Hamiltonian matrix elements are ready via two passes through the control file after the driver program, tbzint, has been called. There is not unlimited scope for specification of the hamiltonian; modern environment dependent matrix elements or complicated functions of distance are not catered for. In such cases the user will have to extend the code, possibly by enabling a disc read of tabulated functions, or an adaptation of the segment tbham.f.

Diagonal matrix elements are read from the START category. Off-diagonal matrix elements of the hamiltonian, overlap and the pair potential are read from ctrl.ext in the ME category, having the following syntax.

ME    is a category that tells the program what hamiltonian matrix
      to use. It begins with an integer, which defines a mode
      to specify the form of matrix elements.
      `mode' is one of:
         0  to indicate that the hopping integrals are fixed
         1  to use the Harrison universal matrix elements for
            semiconductors
         2  for exponential decay
         3  for inverse power decay
         4  general power * exponential decay
         5  for Goodwin Skinner Pettifor (GSP) scaling

Each mode has a group of numbers of input, whose meaning (and number) depend on the mode. But for each mode, the hamiltonian is specified by a pattern like this:

ME  mode
    name1  name2 | (set-of-numbers-specifying-matrix-elements)

name1 and name2 indicate the two species of atoms as referenced in the SPEC category of the ctrl file; the bar (|) indicates that a set of numbers will follow; how many numbers and what they mean depends on mode.

Mode 0: matrix elements are fixed to constant values, independent of, e.g. bond length d. The hamiltonian matrix elements are given as a string of numbers. Following the bar (|) 10 numbers are read in:

     ssσ, spσ, ppσ, ppπ, sdσ, pdσ, pdπ, ddσ, ddπ, ddδ

These are integrals V_pair, read in the order shown, connecting pairs of atoms with species name1 and name2 respectively.

Mode 1: matrix elements follow Harrison’s universal form as described in his book [2]. That is: $V_\mathrm{pair} = \eta_\mathrm{pair} \hbar^2 / md^2$ where

      η_ssσ   = -1.32 eV
      η_spσ   = +1.42 eV
      η_ppσ   = +2.22 eV
      η_ppπ   = -0.63 eV

The user does not specify any rules for matrix elements (they are taken to be universal). Note that this mode is suitable to sp semiconductors only.

Modes 2 and 3: additional information is supplied about the dependence of the matrix element on bond length d. Mode 2 uses an exponential decay; Mode 3 a power law decay.

For each matrix element two coefficients (a and b) are required to define a matrix element: $V_\mathrm{pair} = a e^{-b d} \quad \mathrm{(exponential\ decay)} \qquad V_\mathrm{pair} = a d^{-b} \quad \mathrm{(power\ law\ decay)}$

Coefficients a are specified as a vector of 10 (or fewer) numbers, in the same format as Mode 0 above.
The decay parameters b are specified by a vector following token DECAY=.

Note: For ME modes 2 and 3, DECAY= expects to read the quantities b as positive numbers. tbe will make “sanity checks” on the parameters to identify any sign errors.

For example, in the canonical d band model of Spanjaard and Desjonqueres [7], we would write:

ME      2
        1 1  | 0 0 0 0 0 0 0 -fd*6 fd*4 -fd
               DECAY=0 0 0 0 0 0 0 q q q

You can specify a global default decay coefficient. If an explicit specification is omitted for a particular pair, the default value is used. Specify the default parameter with token DECAY0=, after the ME mode specification. The Spanjaard and Desjonqueres [7] rule could thus have been specifed as:

ME      2  DECAY0=q
        1 1  | 0 0 0 0 0 0 0 -fd*6 fd*4 -fd

ME      2  DECAY0=q
        1 1  | 0 0 0 0 0 0 0 -fd*6 fd*4 -fd

Mode 4: the matrix element takes the form $V_\mathrm{pair} = \sum\nolimits_{k=1{\ldots} 3} a_k d^{b_k} e^{-c_k d}$

where d is bond length. Thus, each matrix element requires 9 numbers as input. The nine numbers after the vertical bar (|) are the three groups of a b c parameters.

Mode 5: The Goodwin-Skinner-Pettifor matrix elements take the form $V_\mathrm{pair} = V (r_0/d)^n \exp[-n(d/r_c)^{n_c}-n(r_0/r_c)^{n_c}]$

In this mode, each matrix element requires 5 numbers, read in the order V n n_c r₀ r_c . For example, a hamiltonian with s and p orbitals would read:

ME  5
      1 1 | sss n nc r0 rc
            sps n nc r0 rc
            pps n nc r0 rc
            ppp n nc r0 rc

where sss=ssσ, sps=spσ, pps=ppσ, ppp=ppπ.

Empirical pair potential: Many tight-binding hamiltonians adopt a classical pair potential (energy) to stabilize the lattice at its equilibrium spacing. A pair potential may be specified by a shriek (!) ; the format follows the power-law-exponential decay format, i.e. $\sum\nolimits_{i=1{\ldots} 3} a_i d^{b_i} e^{-c_i d}$

The nine numbers after the shriek are the three groups of a b c parameters. In the canonical model of Spanjaard and Desjonqueres [7], the hamiltonian and pair potential are thus specified as:

ME      2 
        1 1  | 0 0 0 0 0 0 0 -fd*6 fd*4 -fd  
               DECAY=0 0 0 0 0 0 0 q q q
        1 1  ! b 0 p   0 0 0    0 0 0       

If the power exponents are positive, then the pair potential does not take the above form but Chadi’s form [1] instead, namely a₁ε + a₂ε₂: the third number in each set is the equilibrium bond length.

A line such as this

    ! A 1 -1 n nc r0 rc 0 0

implements a pair potential of the GSP form. The last two numbers should be zero for the standard GWP potential. If they are nonzero, an exponential or power law is added to the GSP potential depending whether the last number is positive (last number is decay) or negative (last number is the power). The second-last coefficient is the prefactor.

Additional syntax enables the input of further parameters. Here is a complete list.

ME    [...]
      name1 name2 | rule for Hamiltonian
      name1 name2 ! rule for the classical repulsive pair potential
      name1 name2 @ rule for overlap
      name1 name2 & rule for crystal field
      name1 name2 % rule for crystal field-overlap

Of the last three all take exactly the same number of parameters as the hamiltonian. (In general an overlap matrix element will take the opposite sign to its corresponding hamiltonian matrix element.) The last two implement Jim Chadi’s empirical crystal field scheme [8]; to a large extent this is now superseded by self consistent tight binding with polarisable ions.

Note that if the basis is non orthogonal, then the overlap matrix elements must also have their scaling specified by DECAY= and CUT= tokens.

The ME category also accepts these additional tokens, which may be specified differently for each species pair.

MEMODE=    In this way the mode can be different for different species pairs
POLY=      The degree of the cut-off polynomial (4 or 5)
CUTMOD=    cutoff mode							     
	   0 no cutoff							     
	   1 augmentative cutoff: all functions augmented with the cutoff polynomial 
	   2 multiplicative cutoff: all functions multiplied by the cutoff polynomial
	   (see poly45.f for more documentation)                                  
CUTPP=     cut off radii, r₁ and r₁, the
           cut-off polynomial is defined on [r₁, r₁]
           so as to match value and slope at r₁ and to have zero
           value and slope at r₂. r₁ and r₁ are expected
           in either bohr (a.u.) or units of ALAT depending on the value
           of the token SCALE= in category **TB**.
PPMODE=    0 no pair potential
           10 a sum of power law times exponential (up to three terms)
           20 Quadratic pair potential (see Ref [[1]](#references))
           21 sum of a polynomial defined in [0 < r < r₂] and two power law functions
           30 GSP
           32 GSP + power law decay 
           33 GSP + (a1*exp(-2r) + a2/r^6)*left_cutoff
           (see code segments _makvp.f_{:.path} and _vppder.f_{:.path} for full documentation)

Example ME category

Below is an example for a case of two species “Fe” and “C” in which the values have been set in the CONST category. Note the following

If MEMODE= is set, there still needs to be a global mode set as the first number after the “ME” string in the ME category.
It is not enough to give the Fe C pair parameters, the C Fe parameters must also be given. Indeed they do not need to be the same. In this way the typical case in semiconductors can be covered: say in GaAs, the s(Ga)-p(As) need not be equal to p(Ga)-s(As); however only spσ and not psσ are specfied in each input line. tbe will take care of symmetrising the hamiltonian in this case.

ME      2
        Fe Fe MEMODE=2 PPMODE=10 POLY=5 CUTMOD=cutmod CUTPP=r1pp rcpp
            | fss fsp fpp -fpp/2 fsd fpd -fpd/sqrt(3)
                                                 fdds fddp fddd
        DECAY=qss qsp qpp qpp    qsd qpd  qpd    qdds qddp qddd
        CUT=r1ss rcss 0 0 0 0 0 0 r1sd rcsd 0 0 0 0 r1dd rcdd r1dd rcdd r1dd rcdd
            @ oss osp opp -opp/2 osd opd -opd/sqrt(3)
                                                 odds oddp oddd
        DECAY=qoss qsp qpp qpp    qosd qpd  qpd    qdds qddp qddd
        CUT=r1ss rcss 0 0 0 0 0 0 r1sd rcsd 0 0 0 0 r1dd rcdd r1dd rcdd r1dd rcdd
            ! b0 m0 p0   b1 m1 p1   0 0 0

        C C MEMODE=CCmode PPMODE=10 POLY=5 CUTMOD=cutmod CUTPP=r1CC rcCC
            | fCCsss fCCsps fCCpps fCCppp 0 0 0 0 0 0
        DECAY=qssCC   qspCC    qppCC    qppCC     0 0 0 0 0 0
        CUT=  r1CC rcCC r1CC rcCC r1CC rcCC r1CC rcCC
              0 0  0 0  0 0  0 0  0 0  0 0
            @ oCCsss oCCsps oCCpps oCCppp 0 0 0 0 0 0
        DECAY=qssCC   qspCC    qppCC    qppCC     0 0 0 0 0 0
        CUT=  r1CC rcCC r1CC rcCC r1CC rcCC r1CC rcCC
              0 0  0 0  0 0  0 0  0 0  0 0
            ! bCC mCC pCC  0 0 0    0 0 0

        Fe C MEMODE=2 PPMODE=10 POLY=5 CUTMOD=cutmod CUTPP=r1CFpp rcCFpp
            | fCFss  fCFsp  0  0 fCFsd fCFpds fCFpdp 0 0 0
        DECAY=qCFss  qCFsp  0  0 qCFsd qCFpds qCFpdp 0 0 0
        CUT=  r1CFs rcCFs r1CFsp rcCFsp 0 0 0 0
              r1CFd rcCFd r1CFd  rcCFd r1CFd rcCFd 0 0 0 0 0 0
            @ oCFss  oCFsp  0  0 oCFsd oCFpds oCFpdp 0 0 0
        DECAY=qoCFss qoCFsp  0  0 qoCFsd qoCFpds qoCFpdp 0 0 0
        CUT=  r1CFs rcCFs r1CFsp rcCFsp 0 0 0 0
              r1CFd rcCFd r1CFd  rcCFd r1CFd rcCFd 0 0 0 0 0 0
            ! b0CF n0CF q0CF   b1CF n1CF q1CF  0 0 0
        C Fe MEMODE=2 PPMODE=10 POLY=5 CUTMOD=cutmod CUTPP=r1CFpp rcCFpp
            | fCFss  fCFsp  0  0 fCFsd fCFpds fCFpdp 0 0 0
        DECAY=qCFss  qCFsp  0  0 qCFsd qCFpds qCFpdp 0 0 0
        CUT=  r1CFs rcCFs r1CFsp rcCFsp 0 0 0 0
              r1CFd rcCFd r1CFd  rcCFd r1CFd rcCFd 0 0 0 0 0 0
            @ oCFss  oCFsp  0  0 oCFsd oCFpds oCFpdp 0 0 0
        DECAY=qoCFss qoCFsp  0  0 qoCFsd qoCFpds qoCFpdp 0 0 0
        CUT=  r1CFs rcCFs r1CFsp rcCFsp 0 0 0 0
              r1CFd rcCFd r1CFd  rcCFd r1CFd rcCFd 0 0 0 0 0 0
            ! b0CF n0CF q0CF   b1CF n1CF q1CF  0 0 0

See Table of Contents

The TB and START categories

Most species-independent, tbe-specific input is placed in category TB. A brief description of the tokens and their contents can be obtained by runnning

tbe --input

The START category is used to read in diagonal hamiltonian matrix elements and the free-atom charges where needed. See the general description for P and for Q. A sample input specific to ASA can be found here.

tbe reads (P,Q) as a set in the ASA format; but the meanings are different. The P parameters must be set for consistency with the ASA, but they are not used. The Q parameters which in the ASA correspond to the 0^th, 1^st, and 2^nd energy moments of the density, correspond in tbe to the number of electrons, the on-site hamiltonian elements and the Hubbard U for each class.

tbe uses a hierarchy of choices to determine the total number of electrons: if ZVAL is supplied (see the BZ category) it will use this number; if not then it will be taken from START as the total 0^th moment; if these are all zero, tbe will use the (compulsory) atomic numbers, Z= in category SPEC.

See Table of Contents

Self consistent, polarisable ion, tight binding

The theory is described in detail in in Mike Finnis’s book [4] and in Tony Paxton’s lecture notes [5]. This branch of tbe is designed to implement the scheme exactly as written down there.

It is invoked by setting UL=T in the TB category. Polarisability of each species is specified through parameters called Δ_l’l’‘l They are entered in ctrl.ext at the token QPOL= in category SPEC which takes up to ten numbers of which the first seven are

        Δ₀₁₁ = Δ₁₀₁
        Δ₁₁₂
        Δ₀₂₂ = Δ₂₀₂
        Δ₁₂₁ = Δ₂₁₁
        Δ₂₂₂
        Δ₁₂₃ = Δ₂₁₃
        Δ₂₂₄

The Hubbard U are read from category START.

Note on units: in non self consistent mode tbe can in principle take input in any consistent set of units (say, eV and angstrom); however if UL=T is set, then input must be in atomic Rydberg units.

Token IODEL= amounts to a “restart” switch: tbe writes the increments of the hamiltonian or the multipole moments of the charge to disc and these can be read back in at the start of a new calculation.

Mixing is implemented in three modes.

The default is to mix the increments to the hamiltonian
The command line option –mxq enables mixing of the charge
The command line option –mxq2 enables mixing of the charge and magnetic moment separately

Mixing parameters are set exactly as in ASA using the category ITER, the only option being Anderson mixing.

Convergence is specified with respect to the rms difference in charge through the token CNVG= in category START, with the maximum number of iterations following the token NIT=, also in category START.

See Table of Contents

Magnetic tight binding

The most straightforward extension of the tight binding theory to magnetic systems, consistent with the local spin density approximation, is to add a term to the second order hamiltonian, $-(1/2) \sum\nolimits_\mathbf{R} J_\mathbf{R} [\delta (q_\mathbf{R}^\uparrow)^2 + (\delta q_\mathbf{R}^\downarrow)^2]$

which results in an additional term in the effective potential (see Ref [10]) $-(1/4)\sum\nolimits_\mathbf{R} J_\mathbf{R} m_\mathbf{R}$

in which $m_\mathbf{R} = \delta q_\mathbf{R}^\uparrow - \delta q_\mathbf{R}^\downarrow$ is the magnetic moment at site R. This is a spin-dependent, but orbital independent potential as long as the Stoner parameter J is averaged over the basis orbitals at site R. The orbital and spin dependent equivalent of LDA+U is TB+U and is still under development.

The Stoner parameters are taken from ctrl.ext in category SPEC at the token JH= which takes nl values which are the J parameters for each l-channel. This is consistent with the syntax used in the LDA+U as implemented in lmf. As with the Hubbard U these are averaged over the site internally.

Mixing can be tricky in magnetic metals and further work may be needed. By setting the command line option –mxq2 it is possible to use the functionality of the MIX structure in the category ITER to make a repeating sequence of mixing charge and magnetic moment separately or together using different Anderson β parameters: token b= sets the mixing for the charge and token bv= sets the mixing for the moment. If either of these are zero for a short sequence then only the other quantity is mixed. The other two schemes, <A HREF=#mixing>above</A>, are implemented except in the case of a non orthogonal basis, in which case only charge mixing is possible. The magnetic symmetry needs to be broken at the start of a calculation (unless the restart option IODEL=T is used and a legal restart file from a magnetic state can be found on disc). For potential mixing (the default) the up and down on-site energies must be split and this is done in ctrl.ext using the token DELTA= in category SITE. For example in the test file ctrl.fecr there are the lines

SITE    ATOM=Cr POS= 0 0 0         DELTA=0 0 -deltaCr 0 0 deltaCr
        ATOM=Fe POS= 1 0 0+dz+ddz  DELTA=0 0 -deltaFe 0 0 deltaFe

which split the the d-electron on-site energies by plus and minus delta; NB this is a ferromagnetic starting state: reversing the signs on one of the atoms would provide an antiferromagnetic starting state. Self consistency will flip these back if there is no local antiferromagnetic solution. In the case of charge mixing, invoked by the command line option –mxq[2] the input charge must have a magnetic moment, achieved by adjusting the input 0th moments in category START as described <A HREF=#moments>above</A>. Here are the lines in the START category from the test file ctrl.fecr. (Again, an antiferromagnetic starting state could be made by using momCr={nsp==1?0:-3}, say.

%const nsp=2 spd=1 NdFe=6 NdCr=4
CONST   q0s={spd?1:0} q0p=0 q0dFe={spd?7:NdFe} q0dCr={spd?5:NdCr}
        esFe=0.2 epFe=0.45 edFe=-0.01 momFe={nsp==1?0:3}
        esCr=0.2 epCr=0.45 edCr=0.01  momCr={nsp==1?0:3}
        U=1 Us=U Up=U UdFe=U UdCr=U
START   CNTROL=T NIT={nitq} CNVG=qtol
        ATOM=Fe   P= 4 4 3 4 4 3
                  Q= q0s/{nsp}            esFe   Us
                     q0p/{nsp}            epFe   Up
                     (q0dFe+momFe)/{nsp}  edFe  UdFe
                     q0s/{nsp}            esFe   Us
                     q0p/{nsp}            epFe   Up
                     (q0dFe-momFe)/{nsp}  edFe  UdFe
        ATOM=Cr   P= 4 4 3 4 4 3
                  Q= q0s/{nsp}            esCr   Us
                     q0p/{nsp}            epCr   Up
                     (q0dCr+momCr)/{nsp}  edCr  UdCr
                     q0s/{nsp}            esCr   Us
                     q0p/{nsp}            epCr   Up
                     (q0dCr-momCr)/{nsp}  edCr  UdCr

Note the use of the value of the token NSPIN={nsp} from category OPTIONS to split the spins in the spin polarised case, yet to provide the correct input moments in the non spin polarised case also (nsp=1).

It’s important to understand the role of the token NSPIN={nsp} in category OPTIONS as interpreted by tbe. If token UL=F in category TB then tbe will do non self consistent tight binding, in which case NSPIN=2 turns on the empirical spin orbit coupled branch of the code (not documented). Alternatively if token UL=T in category TB, tbe will do self consistent tight binding which does not include empirical spin orbit coupling; in that case the token NSPIN=2 in category OPTIONS tells tbe to make a magnetic tight binding calculation.

See Table of Contents

Molecular dynamics and statics

Conventions are in line with the input syntax in lmf and the mol code. Molecular dynamics NVE, NVT and NpT ensembles are implemented in the code segment verlet.f90 using exactly the reversible integrators employing Liouville operators published by Martyna et al. [9]. Here is the relevant section from invoking tbe –input

 --- Parameters for dynamics and statics ---
 DYN_NIT           opt    i4       1,  1          default = 1
   maximum number of relaxation steps (statics) or time steps (dynamics)
 DYN_MSTAT         opt    ---
   Parameters for molecular statics
 DYN_MSTAT_MODE    reqd   i4       1,  1          default = 0
   0: no relaxation  4: conjugate gradients  5: Fletcher-Powell  6: Broyden
 DYN_MSTAT_HESS    opt    lg       1,  1          default = T
   Read hessian matrix
 DYN_MSTAT_XTOL    opt    r8       1,  1          default = 1e-3
   Convergence criterion in displacements
   XTOL>0: use length;  <0: use max val;  =0: do not use
 DYN_MSTAT_GTOL    opt    r8       1,  1          default = 0
   Convergence criterion in gradients
   GTOL>0: use length;  <0: use max val;  =0: do not use
 DYN_MSTAT_STEP    opt    r8       1,  1          default = 0.015
   Initial (and maximum) step length
 DYN_MSTAT_NKILL   opt    i4       1,  1          default = 0
   Remove hessian after NKILL iter
 DYN_MSTAT_PDEF    opt    r8v      6,  1          default = 0 0 0 ...
   Lattice deformation modes
 DYN_MD            opt    ---
   Parameters for molecular dynamics
 DYN_MD_MODE       reqd   i4       1,  1          default = 0
   0: no MD 1: NVE 2: NVT 3: NPT
 DYN_MD_TSTEP      opt    r8       1,  1          default = 20.671
   Time step (a.u.)
 DYN_MD_TEMP       opt    r8       1,  1          default = 0.00189999
   Temperature (a.u.)
 DYN_MD_TAUP       opt    r8       1,  1          default = 206.71
   Thermostat relaxation time (a.u.)
 DYN_MD_TIME       reqd   r8       1,  1          default = 20671000
   Total MD time (a.u.)
 DYN_MD_TAUB       opt    r8       1,  1          default = 2067.1
   Barostat relaxation time (a.u.)
 DYN_MD_P          opt    r8       1,  1          default = 0
   External pressure
 NB: 1 deg.K = 6.3333e-6 a.u.; 1 fs = 20.67098 a.u.

Category DYN refers to both statics and dynamics which are chosen by reference to the token MODE set in either the MD or MSTAT structures. MODE=1-3 is reserved for the three MD ensembles and MODE=4-6 sets the energy minimisation algorithm: Conjugate Gradients, Fletcher-Powell and Broyden respectively. <A HREF=#units>Units</A> must be in atomic Rydbergs: useful conversion factors are included in the segment of output (above). Here’s an example few lines from a typical ctrl file:

DYN
% const dyn=0 temp=300 taup=10 taub=100 time=100000 tstep=0.5
% const hess=F relax=0 nit=50 xtol=1d-3 gtol=5d-4 step=0.01 nkill=100 nitf=50
% const fs=0.048377 K=1/0.6333328d-5
% if dyn==1|dyn==2|dyn==3
        MD[MODE={dyn} TSTEP={tstep/fs} TEMP={temp/K} TAUP={taup/fs}
           TIME={time/fs} TAUB={taub/fs} P=0]
% elseif relax>0
        MSTAT[MODE={relax} HESS={hess} XTOL={xtol} GTOL={gtol}
              STEP={step} NKILL={nkill}] NIT={nitf}
% endif

For MD the mass of the atom is read from category SPEC. As always units are Rydberg a.u. so for example the category may reference a carbon atom as

% const amass=1.09716d-3
        ATOM=C Z=6 R=R/W AMASS=12.0107/{amass}

Band and DOS plotting

The command line option –band[options] turns on the band plotting. The options are as described in ASA documentation. In the case of non self consistent tight binding, tbe –band will only plot bands; if UL=T is set then tbe will make the hamiltonian self consistent before plotting bands. It will be quicker in that case to invoke tbe with IODEL=T. In either case tbe looks for the Fermi energy on disc. This may be overidden on the command line, as in the ASA, with ~ef=#.

Fermi surface plotting and Brillouin zone mapping are not currently implemented in tbe.

For plotting the total DOS, the switch SAVDOS= in category BZ dumps out a file that can be read by the program pldos.f which is in the distribution. There is also an option to make Mulliken decomposed densities which interfaces with the lmdos utility.

The procedure is to set the MULL= token to the required value in category BZ. MULL= takes an integer whose one’s digit translates to a parameter mull1 that describes how the charge on site is deconstructed (by atom, l or lm); whose ten’s digit enables bond order output (mull2); and whose 100’s digit sets a parameter imfil that points to an ASCII input file holding details of what decomposition is required. The specification of these three parameters is as follows (see also the segment mkwttb.f)

   For (mull=0,1; imfil=0) all the partial DOS channels are generated
     automatically.  All of the angular momentum channels for a given
     atom appear consecutively in the DOS file (see below for specific
     ordering).  The atom ordering follows the class indices for
     (mull=0) and the basis indices for (mull=1).  In the case of
     (nsp=2) all of the spin up channels for a given atom or class
     appear consecutively followed by all of the spin down channels
     for the same atom or class.

  For (mull=0,1; imfil=1) the partial DOS channels are specified in
    the file MULL using the class indices (mull=0) or the basis
    indices (mull=1) and the angular momentum indices below.  The
    format of the lines in this free format file is:

      <ang mom ind 1> <atom ind 1>
      ...

    In the case of (nsp=2) the choice of spin up or down is specified
    in the file MULL by choosing the appropriate lm index (see below).

  For (mull=10,11; imfil=0) the bond charge is automatically
    generated for all atom pairs.  The ordering of the channels is
    done according to the class indices for (mull=10) and the basis
    indices for (mull=11).  The ordering is as follows:

         (mull=11)            (mull=10)

      Atom 1   Atom 2      Atom 1   Atom 2
      ------   ------      ------   ------
         1        1           1        1    On-site contribution
         1        2           1        1    All others (if existent)
         1        3           1        2
       ...      ...         ...      ...
         2        2           2        2    On-site contribution
         2        3           2        2    All others (if existent)
         2        4           2        3
       ...      ...         ...      ...

  For (mull=10,11, imfil=1) the bond charge DOS channels are specified
    in the file MULL using the class indices (mull=10) or the basis
    indices (mull=11).  The format of the lines in this free format
    file is:

      <atom ind 1> <atom ind 2> <ionsit>
      ...

     The switch (ionsit) must always be present and must be equal to
     zero unless [mod(mull,10)=0] and the two atom class indices are
     the same.  In this case (ionsit=1) refers to the on-site
     contribution and (ionsit=0) refers to all other contributions
     (if they exist).

   For (mull=20,21) and (mull=30,31) the atoms are indexed according
     to the class index (mull=20,30) or the basis index (mull=21,31).
     Ths DOS channels are specified in the file MULL using these atom
     indices and the angular momentum indices indicated below.  The
     format of the lines in this file is (free format):

       <lm ind 1> <atom ind 1> <lm ind 2> <atom ind 2> <ionsit>
       ...

     The switch (ionsit) is the same as for (mull=10,11; imfil=1)

     In the case of (nsp=2) the spin combination is specified in the
     file MULL by choosing the appropriate lm indices (see below).

   The angular momentum channels for each atom or class for
     (mull=0,20,21) are indexed as follows:
       lm = 1 : s
       lm = 2 : p (summed over m = -1, 0, 1)
       lm = 3 : d (summed over m = -2, -1, 0, 1, 2)
       ...

     In the case of (nsp=2) all of the spin up lm indices appear
     consecutively followed by all of the spin down lm indices
     (i.e. if nl=2 then lm=4 corresponds to spin down p).

   The angular momentum channels for each atom or class for
     (mull=1,30,31) are indexed as follows:
       lm = 1     : s
       lm = 2,3,4 : p_x, p_y, p_z
       lm = 5-9   : d_xy, d_yz, d_zx, d_(x²-y²), d_(3z²-r²)
       ...

     In the case of (nsp=2) all of the spin up lm indices appear
     consecutively followed by all of the spin down lm indices
     (i.e. if nl=2 then lm=6 corresponds to spin down p_x).

   The DOS channels are symmetrized if needed [ng > 0 and
     mod(mull,10) = 1] for (mull=1,11,21).  No symmetrization is done
     for (mull=30,31) and for (mull=1, L > 2) and therefore the full
     BZ should be used (no special points) or symmetry-equivalent
     channels should be averaged.

Once tbe has run, construct the DOS using the lmdos utility:

lmdos --dos:tbdos:[other switches ..] ext

In the spin polarised tight binding, use

lmdos --dos:tbu:[other switches ..] ext

When invoking lmdos, take care to use the same command line switches as used to invoke tbe, if for example the number of k-points or spins is different from that specified in ctrl.ext. The pldos utility may used to produce a file dosp.dat

pldos -fplot dos.ext

Draw a figure with the fplot utility, or, since dosp.dat contains columns of data, use your own favorite plotting package.

fplot -disp -pr10 -f plot.dos

tbe-specific command line switches

See this page

References

D. J. Chadi, Phys. Rev. Lett., 41, 1062 (1978)
W. A. Harrison, Electronic structure and the properties of solids, the physics of the chemical bond, (W. H. Freeman, San Francisco) 1980 (now published by Dover)
M. W. Finnis, A. T. Paxton, M. Methfessel and M. van Schilfgaarde, Phys. Rev. Lett. 81, 5149 (1998)
M. W. Finnis, Interatomic forces in condensed matter, (OUP, Oxford) 2003
A. T. Paxton, in Multiscale Simulation Methods in Molecular Sciences, NIC series, Vol. 42, edited by J. Grotendorst, N. Attig, S. Bluegel, and D. Marx (Institute for Advanced Simulation, Forschungszentrum Julich) 2009 pp. 145-174. Available on-line at [http://webarchiv.fz-juelich.de/nic-series/volume42/volume42.html[(http://webarchiv.fz-juelich.de/nic-series/volume42/volume42.html).
A. T. Paxton and M. W. Finnis, Phys. Rev. B, 77, 024428 (2008)
D. Spanjaard and M. C. Desjonqueres, Phys. Rev. B, 30, 4822 (1984)
D. J. Chadi, in Atomistic simulation of materials: beyond pair potentials,, ed. V. Vitek and D. J. Srolovitz, (Plenum Press, New York), p. 309, 1989 (ISBN-l3: 978-1-4684-5705-6, e-ISBN-13: 978-1-4684-5703-2 DOl: 10.1007/ 978-1-4684-5703-2)
G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Molecular Physics, 87, 1117 (1996)
J. Magnetism and Magn. Materials</i>., 15-18, 847 (1980))