Written by M. van Schilfgaarde, J. E. Klepeis and A. T. Paxton Bugs and comments to Tony.Paxton@QUB.ac.uk Documentation for Version 9.0 Oct 2008
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
#1 #2 | (set-of-numbers-specifying-matrix-elements)
#1 and #2 indicate the two classes of atoms; the bar (|) indicates
that a set of numbers will follow; how many numbers and what they mean depends on mode.
sssigma, spsigma, ppsigma, pppi,
sdsigma, pdsigma, pdpi, ddsigma,
ddpi and dddelta
These are integrals Vij, read in the order shown.
Vij = ηij hbar2 / md2
ηss&sigma = -1.32 eV
ηsp&sigma = +1.42 eV
ηpp&sigma = +2.22 eV
ηpp&pi = -0.63 eV
In this mode, the use does not specify any rules for matrix elements. Note that
this mode is suitable to sp semiconductors only.
V = a×e-b d (exponential decay) or V = a×d-b (power law decay)Coefficients a are specified as a vector of 10 (or fewer) numbers, in the same format as Mode 0 above.
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 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
V n nc r0 rc
For example, a hamiltonian with s and p orbitals would read:
ME 5
1 1 | sssigma n nc r0 rc
spsigma n nc r0 rc
ppsigma n nc r0 rc
pppi n nc r0 rc
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
! 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.
ME [...]
i-list j-list | rule for Hamiltonian
i-list j-list ! rule for the classical repulsive pair potential
i-list j-list @ rule for overlap
i-list j-list & rule for crystal field
i-list j-list % rule for crystal field-overlap
(i-list and j-list are explained below.) 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; to a large extent this is now superseded by
self consistent tight binding with polarisable ions.
Empirical spin orbit coupling is enabled by setting NSPIN=2; exactly two spin orbit coupling parameters being read in at the token VSO= in category SPEC. An example test file ctrl.fecr invokes all of the above features.
The token CLSTYL= allows more flexibility to what pairs are associated with a specified rule for matrix elements. This is useful if there is a large number of species. To invoke, ME has optional CLSTYL, which must follow memode:
ME memode [CLSTYL=0 or 1 or 2] [DECAY0=#]
i-list j-list | rule1
i-list j-list | rule2
...
CLSTYL=1 uses the original convention for i-list and
j-list, namely a list of class numbers using the usual syntax
in mkilst.f, eg i-list looks like 1,2,7:10
CLSTYL=2 takes any class satisfying an expression involving variables
ic (the class number) and z (the class atomic number), e.g.
i-list looks like, e.g.
ic<10&z==6
It this case any of the first ten classes with Z=6 will be included
in i-list.
CLSTYL=3, for unix systems, uses file names with wild cards. For each class with name 'nam' a file is created, file 'nam.ext' You enter i-list as a class name with wild cards, eg i-list looks like a* Any class whose corresponding file the unix shell finds by expanding the file 'i-list.ext' (in this case all classes beginning with 'a') are included in the list.
Each species of atom may be given basis function orbitals up to l=2, that is, s, p and d. To implement f orbitals would require an impractically large rewrite of the code. The token NL= specifies the maximum number of l's at any site, and the tokens IDXDN= are used to specify which orbitals are required for each species: as in the ASA IDXDN=1 includes an orbital in the basis and IDXDN=3 excludes it. Thus if NL=3, a species with IDXDN=1 3 1 will have an sd basis, while a species with IDXDN=3 1 3 has p orbitals only.
tbe --inputwith some additional remarks.
TB_OVLP opt lg 1, 1 default = F Non orthogonal tight-binding TB_CRYSF opt lg 1, 1 default = F Crystal field terms in hamiltonian TB_OVCF opt lg 1, 1 default = F Crystal field terms in overlap TB_ADDES opt lg 1, 1 default = F Add ebarLL' * sLL' to hLL' TB_RMAXH opt r8 1, 1 default = 0 Hamiltonian cut-off length in units of a TB_FORCES opt lg 1, 1 default = T Calculate forces TB_FIJ opt lg 1, 1 default = F To get forces by atom pair, e.g., for stresses TB_3PV opt lg 1, 1 default = F Calculate pressure TB_EVDISC opt lg 1, 1 default = T Can be F for insulators or to save space for metals TB_PAIR opt lg 1, 1 default = F Pair potential only TB_TRH opt lg 1, 1 default = F Calculate local projection of band energies TB_RHO opt lg 1, 1 default = F Calculate local projection of charges TB_U1 opt lg 1, 1 default = F For electrostatics with L>=0 (first model) TB_UL opt lg 1, 1 default = F For electrostatics with L>=0 TB_TBU opt lg 1, 1 default = F Tight-binding+U TB_GAMMA opt lg 1, 1 default = F Do gamma-point only TB_MOL opt lg 1, 1 default = F molecule: no PBCs
The START category is used to read in diagonal hamiltonian matrix elements and the free-atom charges where needed. See the general description of parameters P and Q in the ASA documentation, and here for the format of input tokens. 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 moments, 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 set in the BZ category it will use this number; if not then it will taken from START as the total 0th moment; if these are all zero, tbe will use the (compulsory) atomic numbers, Z= in category SPEC.
It is invoked by setting UL=T in category TB. Polarisability of each species is specified through parameters called &Deltal'l"l They are entered in ctrl at the token QPOL= in category SPEC which takes up to ten numbers of which the first seven are
&Delta011 = &Delta101
&Delta112
&Delta022 = &Delta202
&Delta121 = &Delta211
&Delta222
&Delta123 = &Delta213
&Delta224
internally, these are multiplied by The Hubbard U are read from category START (see above). Finnis' formulation does not allow for an l-dependent U, which is consistent with density functional theory. If U in START are not the same in each l-channel then they are averaged internally over the orbitals for which IDXDN=1 on that species. This amounts to setting the token NOUAVG=F in category TB. The case NOUAVG=T constitutes an orbital-dependent effective potential as found in theories such as SIC, LDA+U. This implementation is under development (Feb 2007) under the general heading TB+U and isn't ready yet.
Important 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 two modes.
Mixing parameters are set exactly as in ASA using the category MIX, the only option being Anderson mixing. A line like this will do
%const nmix=5 kmix=10 beta=1
MIX MODE=A{nmix},k={kmix},b={beta}
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.
-(1/2)&sumR JR ( &deltaqR(up)2 + &deltaqR(down)2 )
which results in an additional term in the effective potential (see Pettifor,
J. Magnetism and Magn. Materials., 15-18, 847 (1980))
-(1/4)&sumR JR mR
in which mR=&deltaqR(up) -
&deltaqR(down) 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 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 implementation in lmf (q.v.); 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. Currently the two schemes, above, 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 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 the input
charge must have a magnetic moment, achieved by adjusting the input 0th
moments in category START as described above. 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. Developers should note that internally this case causes variables nsp and nspc to become respectively 1 and 2. nspc=2 instructs the code to make a single diagonalisation at each k-point of a hamiltonian twice the dimension of the non spin polarised system. The cross terms contain the spin orbit coupling parameters and will also be needed in any future development that implements non collinear magnetism. 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. Now, as in standard scalar relativistic LSDA, the hamiltonian is set up and diagonalised twice at each k-point and the spins are uncoupled. (NB, fully relativistic LSDA with spin orbit coupling is now implemented in the lmf code (q.v.))
category DYN (optional): parameters for dynamics
token MSTAT: of cast struc, elements: relax hess xtol gtol step nkill (optional)
Parameters for molecular statics:
relax=0 (no relaxation) 4 (conjugate gradients) 5 (Fletcher-Powell)
6 (Broyden)
hess=T read hessian matrix
xtol=# convergence criterion in displacements
>0: use length; <0: use max val; =0: do not use
gtol=# convergence criterion in gradients
>0: use length; <0: use max val; =0: do not use
step=# step length
nkill=# Remove hessian after # iter
token MD: of cast struc, elements: dyn tstep temp rltime time (optional)
Parameters for molecular statics:
dyn=0 (no relaxation) 1 (molecular dynamics)
tstep=# time step (a.u.)
temp=# temperature (a.u.)
rltime=# Nose-Hoover relaxation time (a.u.)
time=# total MD time (a.u.)
NB: 1 deg.K = 6.3333e-6 a.u.; 1 fs = 20.67098 a.u.
NB: 1 deg.K = 6.333e-6 a.u.; 1 sec = 20.671e15 a.u.
token NIT= of cast integer (optional) :
maximum number of relaxation steps (MSTAT or MD)
Category DYN refers to both statics and dynamics which are chosen by reference
to tokens relax= and dyn=. Units 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 ts=2 time=10000 temp=300 rt=20 fs=0.048377 K=1/0.6333328d-5
% const relax=0 nitf=100 xtol=1d-4 gtol=1d-4 hess=T step=0.01 nkill=100
% if dyn==1
MD:{dyn},{ts/fs},{temp/K},{rt/fs},{time/fs}
% elseif relax>0
MSTAT:{relax},{hess},{xtol},{gtol},{step},{nkill} NIT={nitf}
% endif
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 ASA program lmdos.f (q.v.) 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:
...
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:
...
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):
...
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 the DOS is to be contructed using
lmdos --dos:tbdos:[other switches ..] extusing the program lmdos from the ASA suite. In the spin polarised tight binding, use
lmdos --dos:tbu:[other switches ..] extWhen 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. The program pldos.f in the distribution FPLOT.*.gz may used like
pldos -fplot dos.extto produce a file dosp.dat containing columns of data that may be used in a third party plotting package.
ccomp -uLM -dTBE lm.f tbe.fIt must then be linked as lm.f against the archive tb/subs.a. The latter is not yet automatically made by configure but is easily created as follows.
% cd $HOME/<lmto distribution directory> (the one containing configure) % cd tb % ln -s ../Make.inc . % ../startup/Maksmakfile > Makefile % make subs.a
The following are command line options specific to tbe.
--mxq do charge mixing --lc assume monopole moments only in the mixing --diagn assume a diagonal density matrix in TB+U --point point charge self consistency only (saves time in Ewald summation) --efield=#,#,# (MOL=T only) apply an electric field (atomic Rydberg units)