Empirical tight-binding program TBE
Table of Contents
Simple unit cell
Producing the ctrl file
This information will be contained within the file ctrl.fe, where fe may be replaced by any string, and is used to name the file. The file describing the first, simplistic calculation would be as follows. This text may be copied directly into a ctrl file to run this simulation. Details of input categories may be found here.
VERS TB=10 LM=7 FP=7
IO SHOW=0 HELP=F VERBOS=31 WKP=F
%const amass=1.09716d-3 r0=5.423516
STRUC ALAT={r0} NBAS=2
PLAT= 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
SPEC ATOM=Fe Z=26 AMASS=55.845/{amass} NL=3 IDXDN=3 3 1
SITE
ATOM=Fe POS=0.0 0.0 0.0 RELAX= 1 1 1
ATOM=Fe POS=0.5 0.5 0.5 RELAX= 1 1 1
HAM NSPIN=1
BZ NKABC=4 4 4 METAL=3 TETRA=F EF0=0.643 DELEF=0.01 N=1 W=0.002
TB RMAXH=1.01*{r0} FORCES=1
ME 0
Fe Fe | 0 0 0 0 0 0 0 -0.6 0.4 -0.1
START CNTROL=T
ATOM=Fe P=4 4 3
Q=0 0.15 1
0 0.45 1
6 0.00 1
ITER
CONV=1d-4
CONVC=5d-7
NIT=100
MIX=A5,b=0.5
Constants declared by % const are variables parsed by the preprocessor. Alternatively they may be defined in the CONST category, the only difference being that they are then referenced without {} brackets.
The categories STRUC and SITE define the structure of the cell. The STRUC section includes information regarding the lattice structure such as lattice constant (ALAT=) and its primitive lattice translation vectors (PLAT=). The SITE section of the init file includes the basis information, including the atom species, coordinates (defined here in fractions of the lattice constant) and the directions along which relaxation is permitted (RELAX=).
SPEC will define the details of the species involved, such as the l-cutoff (NL) and the orbitals that will be included in the calculation (IDXDN), where a value of 1 leaves the orbital in the basis and a value of 3 excludes it. These iron atoms will have only d-orbitals, for simple sp-carbon atoms it would be:
ATOM=C Z=6 AMASS=12.0107/{amass} IDXDN=1 1 3
TB includes the Tight Binding specific switches, such as FORCES, which will define whether or not the forces are calculated, and RMAXH which sets the range for finding neighbours in the lattice.
ME defines the matrix elements for the atomic interactions. The value of zero after ME defines the hopping integrals as fixed during the calculation, where their values are set in the string of numbers following each pair of atoms. Here all the interactions have been set to zero apart from the d-orbital sigma, pi and delta integrals which are set to their canonical values.
START is primarily specific to the ASA, however the Q values are used by the TBE, while the P are not but must be set for consistency with the rest of the code. The Q values are used to set the number of electrons, the onsite energy and the Hubbard U for the s, p and d orbitals respectively. Here it defines 6 electrons in the d-channel, with an onsite energy of 0.00 and a Hubbard U equal to 1.
The ITER category is essentially identical to the rest of the LMTO, and so is well described within the documentation, however there are some minor differences. When TBE is called the mixing type may be defined in the command line; the default scheme is to mix elements of the Hamiltonian, but it may be more effective to use tbe --mxq to mix charges or tbe --mxq2 to mix charges and magnetic moments separately. If --mxq2 is being used with Broyden mixing, w=w1,w2 will still set the weights used for charge and magnetic moment mixing respectively, however if wa=1 is set it will carry out an alternate method called “clock mixing”, whereby it will mix the charges to self consistency before each iteration of magnetic mixing. Now CONV will set the charge rms difference requirement for self consistency and CONVC will set that for magnetic moments. This method can generally produce better converged magnetic moments than the alternatives.
Running the simulation
The TBE executable may be invoked (while in the same folder as the ctrl file) with charge mixing specified, using the command:
tbe --mxq ctrl.fe
First the structure will be read, the interaction parameters determined and the k-point map built, reducing the number of required k-points using the symmetry of the lattice unless --nosym is specified in the command line. At each iteration the eigenvalues at selected k-points will be printed:
SECMTB: kpt 1 of 10, k= 0.00000 0.00000 0.00000
-3.9667 -3.9667 -0.2333 -0.2333 0.1556 0.1556 0.1556 2.6444 2.6444
SECMTB: kpt 10 of 10, k= 0.50000 0.50000 0.50000
-1.4000 -1.4000 -1.4000 -1.4000 -1.4000 -1.4000 2.1000 2.1000 2.1000
The charges will then be mixed as follows:
QMIX mixing multipole moments:
Anderson iteration 1. 2 elements; mixed 0 of 0, beta=0.5, rms diff: 8.88178e-16
There are 2 elements, or the charge on each atom, and zero prior iterations to use in the next estimation of the charges as this is the first iteration. rms diff is the convergence criterion defined in MIX and so the defined limit is 1d-4.
The forces will be calculated once the charges have been mixed to self consistency, however they will obviously be zero here as this is a perfect lattice.
Site pos force
1 Fe 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
2 Fe 0.50000 0.50000 0.50000 0.00000 0.00000 0.00000
Finally it will print the results of the calculation, giving the trace of the density matrix Tr[rho] (or the total number of electrons), the trace of the density matrix with the Hamiltonian Tr[rho][H_0] (which will give the band structure energy), the second order correction to the energy, due to the charge transfer terms and so on, E_2, the internal Virial 3PV and the electrostatic energy emad. A description of the rest of these terms may be found in the source code file tbtote.f.
Note that due to the simplicity of this calculation most of these values are zero, but as the calculation becomes more complex with the inclusion of a pair potential, multiple spins, charge transfer, polarisation and so on, more values will be calculated.
Tr[rho] : 12.00000000
Tr[rho][H_0] : -11.12658182
band structure energy : -11.12658182
E_2 : -0.00000000
pair potential energy : 0.00000000
reference energy : 0.00000000
total energy : -11.12658182
free energy : -11.12676850
3PV pair : 0.00000000
c mmom=0 sev=-11.126582 erep=0 etot=-11.126582 emad=0 3pv=0 TS=.000187
Increasing complexity
Producing the ctrl file
Now attempt relaxation of a larger structure, including magnetism and a more complex description of the hopping integrals.
The DYN category can be used to relax a static structure as follows:
DYN
MSTAT[MODE=5 HESS=F XTOL=1d-3 GTOL=1d-3
STEP=0.01 NKILL=100] NIT=1000
MODE=5 sets the relaxation scheme to Fletcher-Powell and the tolerance for displacement of atoms and forces are set with XTOL and GTOL respectively.
A more interesting structure to relax may be a 54 atom super-cell of BCC iron with one atom removed to produce a vacancy:
STRUC NBAS=53 NSPEC=1 NL=3 ALAT=5.423516
PLAT=3 0 0 0 3 0 0 0 3
SITE
ATOM=Fe POS=0 0 0
ATOM=Fe POS=1.0 0 0
ATOM=Fe POS=0 1.0 0
ATOM=Fe POS=1.0 1.0 0
ATOM=Fe POS=2.0 0.0 0
ATOM=Fe POS=2.0 1.0 0
ATOM=Fe POS=2.0 2.0 0
ATOM=Fe POS=0 2.0 0
ATOM=Fe POS=1.0 2.0 0
ATOM=Fe POS=0.5 0.5 0.5
ATOM=Fe POS=1.5 0.5 0.5
ATOM=Fe POS=0.5 1.5 0.5
ATOM=Fe POS=1.5 1.5 0.5
ATOM=Fe POS=2.5 0.5 0.5
ATOM=Fe POS=2.5 1.5 0.5
ATOM=Fe POS=2.5 2.5 0.5
ATOM=Fe POS=1.5 2.5 0.5
ATOM=Fe POS=0.5 2.5 0.5
ATOM=Fe POS=0 0 1.0
ATOM=Fe POS=1.0 0 1.0
ATOM=Fe POS=0 1.0 1.0
ATOM=Fe POS=1.0 1.0 1.0
ATOM=Fe POS=2.0 0.0 1.0
ATOM=Fe POS=2.0 1.0 1.0
ATOM=Fe POS=2.0 2.0 1.0
ATOM=Fe POS=0 2.0 1.0
ATOM=Fe POS=1.0 2.0 1.0
ATOM=Fe POS=0.5 0.5 1.5
ATOM=Fe POS=1.5 0.5 1.5
ATOM=Fe POS=0.5 1.5 1.5
ATOM=Fe POS=1.5 1.5 1.5
ATOM=Fe POS=2.5 0.5 1.5
ATOM=Fe POS=2.5 1.5 1.5
# ATOM=Fe POS=2.5 2.5 1.5
ATOM=Fe POS=1.5 2.5 1.5
ATOM=Fe POS=0.5 2.5 1.5
ATOM=Fe POS=0 0 2.0
ATOM=Fe POS=1.0 0 2.0
ATOM=Fe POS=0 1.0 2.0
ATOM=Fe POS=1.0 1.0 2.0
ATOM=Fe POS=2.0 0.0 2.0
ATOM=Fe POS=2.0 1.0 2.0
ATOM=Fe POS=2.0 2.0 2.0
ATOM=Fe POS=0 2.0 2.0
ATOM=Fe POS=1.0 2.0 2.0
ATOM=Fe POS=0.5 0.5 2.5
ATOM=Fe POS=1.5 0.5 2.5
ATOM=Fe POS=0.5 1.5 2.5
ATOM=Fe POS=1.5 1.5 2.5
ATOM=Fe POS=2.5 0.5 2.5
ATOM=Fe POS=2.5 1.5 2.5
ATOM=Fe POS=2.5 2.5 2.5
ATOM=Fe POS=1.5 2.5 2.5
ATOM=Fe POS=0.5 2.5 2.5
The atom was removed by inserting # before it, which will cause the preprocessor to ignore this line.
Due to the more complex structure it is advisable to change the mixing parameters in order to ensure convergence:
ITER CONV=1d-4 CONVC=1d-3 NIT=1000 MIX=b13,b=0.008,bv=0.008,wc=-1,w=-1,-1,wa=1
Now TBE will use the Broyden mixing scheme, with a smaller beta and an increased number of prior iterations used in each estimate to improve stability. If --mxq2 is specified in the command line the weights for charge and spin mixing are set separately, with w=w1,w2, otherwise only wc will be used. Note that the weights are set to negative numbers, this causes a new weight to be calculated on each iteration using the previous RMS. As mentioned previously, setting wa=1 means that “clock mixing” will be used for the charges and spins. Obviously this mode may only be used with --mxq2.
HAM NSPIN=2 ELIND= -0.7 PWMODE=11 PWMAX=5 PWMIN=1
SPEC
ATOM=Fe Z=26 R=R I=0.05 A=0.025 AMASS=55.845/{amass}
IDU= 0 0 0 0 UH= 0 0 0 0 JH=0.05 0.05 0.05 0.05
COLOUR=0.1 0.1 0.1 RADIUS=0.5
IDXDN=3 3 1
START CNTROL=T
ATOM=Fe P= 4 4 3 4 4 3
Q= 0 0.15 1
0 0.45 1
(8.8)/2 0.0 1
0 0.15 1
0 0.45 1
(4.8)/2 0.0 1
In order to include magnetism the number of spins must be increased to 2 in HAM. The spin channels must now be specified separately using two sets of Q values.
ME 2
Fe Fe MEMODE=2 PPMODE=10 POLY=5 CUTMOD=2 CUTPP=4.88116 7.592922
| -0.35 0 0 0 -0.14067 0 0 -2.43826 1.99715735 -0.90723918
DECAY=0.3 1 1 1 0.3 1 1 0.9 0.9 0.9
CUT=5.9658 10.847 0 0 0 0 0 0 5.96586 10.847 0 0 0 0 4.88116 7.5929 4.88116 7.5929 4.88116 7.5929
@ 0.45 0 0 0 0.5 0 0 0 0 0
DECAY=0.3 1 1 1 0.3 1 1 0.9 0.9 0.9
CUT=5.9658 10.847 0 0 0 0 0 0 5.96586 10.847 0 0 0 0 4.88116 7.5929 4.88116 7.5929 4.88116 7.5929
! 683.1 0 1.5376 -459.5 0 1.4544 0 0 0
TB FORCES=1 EVDISC=T RMAXH=7.59292 TRH=T RHO=T 3PV=1
UL=1 IODEL=0 OVLP=0
SYMGRP find
BZ NKABC=8 TETRA=T METAL=3
EFMAX=5
The interaction between iron atoms has been set to an exponential form, where the first set of numbers (after “|”) define the values for the interactions at zero range, while DECAY sets the exponent of the exponential decay function for each interaction.
If the exponential interactions were left to decay normally an excessive number of neighbours would be included, and so the tail is replaced or scaled with some rapidly decaying function. The general form of this function is defined by POLY (the degree of the cut-off polynomial) and CUTMOD (whether the function is augmentative or multiplicative), while CUT defines the range at which the function starts and where it goes to zero for each interaction.
The values after “@” set the overlap, which has a similar form for its decay parameter as the interactions, however this is not relevant here as the d-orbitals are too localised for an overlap to be required and the s-orbitals have not been included. PPMODE defines the form of the pair potential, here set to a sum of power law multiplied by an exponential, the values of which are set after “!”.
Finally the value of RMAXH must be changed to reflect the new range of the interactions.
Running the simulation
This time, in order to mix the charges and spins separately, the simulation must be run with:
tbe --mxq2 ctrl.fe
First the charges will be mixed to a self-consistency of at least CONV:
Broyden iteration 62: 53 charges; mixed 13 of 13, rms diff: 2.61389e-5 (tol: 1e-4)
Before the moments are mixed once:
Broyden iteration 63: 53 spins; mixed 1 of 13, rms diff: 2.8871 (tol: 1e-3)
This will continue until the spins reach a consistency of CONVC, then the forces on each atom will be calculated:
Forces on atom 1 Species: Fe
Coordinates : 0.00000000 0.00000000 0.00000000
From bands : 0.00007160 0.00007160 0.00000000
From e'stx : 0.00000201 0.00000201 0.00000000
From pairs : -0.00000000 -0.00000000 0.00000000
Total : 0.00007361 0.00007361 0.00000000
And the maximum force:
Maximum force= 0.013959 on atom 25 (Fe)
This value will be compared with GTOL to determine whether the structure is relaxed. If not, the process gradzr will update the site postions using the Fletcher-Powell scheme:
gradzr: begin F.P. xtol=1e-3 and gtol=1e-3 gtll=0.25 dxmx=0.01 isw=161
gradzr new line 1: g.h=-0.00271 g.(h-g)=0 max g=-0.0135 |grad|=0.0521
p= 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000
g= -7.3610e-5 -7.3610e-5 -6.367e-21 -6.832e-21 -5.1871e-4 6.9440e-18
h= 7.36102e-5 7.36102e-5 6.3673e-21 6.8323e-21 5.18712e-4 -6.944e-18
rfalsi: new start (c) xtol=0.074 dxmn=0.037 dxmx=0.74
rfalsi ir=-1: seek xn=0.73960125
RELAX line 1: new line minimization; max shift=0.01000 |g|=0.0521
And will print out the new coordinates:
Updated atom positions:
Site Class Position(relaxed)
1 Fe 0.00005444(T) 0.00005444(T) 0.00000000(T)
Before mixing begins again with the atoms in their new positions, this will continue until the structure is sufficently relaxed, then the final energies and coordinates will be printed:
Tr[rho] total : 360.40000000
moment : 153.91791280
Tr[rho][H_0] up : -0.93026638
down : -9.09863081
total : -10.02889719
band structure energy : -21.16528983
E_2 : -5.58942732
pair potential energy : -0.78950422
reference energy : 0.00000000
Stoner magnetic energy : -5.58999333
Magnetic moment : 153.91791280
total energy : -16.40782873
free energy : -16.40783056
3PV pair : 42.22590387
bands : -45.87326581
charges : -0.00061289
total : -3.64797484
This tutorial covers the very basics of running a TBE simulation, however a great many features have not been mentioned, such as the non-orthogonal TB, polarisable ion self-consistency and molecular dynamics. A description of these features may be found in the documentation and source code.