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:
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 Vpair, 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
η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:
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 nc r0 rc . 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 a1ε + a2ε2: 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, r1 and r1, the
cut-off polynomial is defined on [r1, r1]
so as to match value and slope at r1 and to have zero
value and slope at r2. r1 and r1 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 < r2] 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
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 0th, 1st, and 2nd 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 0th 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 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.
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_(x2-y2), d_(3z2-r2)
...
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))