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

### Table of Contents

- Background, and specification of matrix elements
- 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:

MEmodename1 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: where

η= -1.32 eV_{ssσ}η= +1.42 eV_{spσ}η= +2.22 eV_{ppσ}η= -0.63 eV_{ppπ}

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:

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

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

*Mode 4*: the matrix element takes the form

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

In this mode, each matrix element requires 5 numbers, read in the order ** V n n _{c} r_{0} 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.

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 _{1}ε* +

*a*: the third number in each set is the equilibrium bond length.

_{2}ε_{2}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 (seepoly45.ffor more documentation) CUTPP= cut off radii,rand_{1}r, the cut-off polynomial is defined on [_{1}r,_{1}r] so as to match value and slope at_{1}rand to have zero value and slope at_{1}r._{2}rand_{1}rare expected in either bohr (a.u.) or units of_{1}ALATdepending 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**{: style="color: green"} and **vppder.f**{: style="color: green"} for full documentation)_{2}

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

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

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

Δ_{011}= Δ_{101}Δ_{112}Δ_{022}= Δ_{202}Δ_{121}= Δ_{211}Δ_{222}Δ_{123}= Δ_{213}Δ_{224}

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

### 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,

which results in an additional term in the effective potential (see Ref [10])

in which is the magnetic moment at site **R**. This is a spin-dependent, but orbital *in*dependent 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{: style=”color: green”} 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{: style=”color: green”}. (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.

### 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^{2}-y^{2}), d_(3z^{2}-r^{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=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))