\chapter{Computing the transport function}
\label{chap:computing_phi}
The computation of the transport parameters of real materials as described in
the previous chapters has been implemented in a computer program which will be
described in this chapter.  The program is quite large, over 7000 lines, and
is written in Fortran 90 which allows for greater flexibility and modularity
than the standard Fortran 77.  This means that almost all arrays used in the
program are {\it dynamically allocated} and related variables are grouped
together in {\it MODULES} which serve a purpose similar to the common blocks
in Fortran 77 but are easier to use and less likely to cause errors.  

The program relies heavily on the {\it LAPACK} program library of linear
algebra subroutines which on Dec Alpha computers is included in the
pre-installed {\it DXML} library.  We have however written an interface to all
the routines that are called from the LAPACK library.  This is useful since
many of the LAPACK routines need work-arrays and job parameters and these are
now be allocated in the interface routines thereby minimizing the hassle and
errors associated with using the library routines.  Also if the program needs
to be compiled with a different linear algebra library only the interface
needs to be changed but actual calls to the interface routines will remain the
same.

\section{What is computed and unit conventions}
The main output of the program is the transport function,
\begin{equation}
  \label{eq:whatIsComputed}
  \phi^{xx}(\epsilon) = \frac{1}{V_{cell}}\sum_{k}\trace
  \left\{v^{x}_{k}\rho_{k}(\epsilon)v^{x}_{k}\rho_{k}(\epsilon)\right\},
\end{equation}
or more specifically a symmetrized version of this expression.  Here
$V_{cell}$ is the volume of the unit cell.  The units used in the program are
electron volt $[eV]$, Bohr radius $[a_0]$ and $[\hbar]$ so the units of the
spectral density $\rho$ are $[1/eV]$ and the units of the velocity are $[eV
a_0/\hbar]$, since $k$ is measured in units of $[1/a_0]$ and the velocity is
the momentum derivative of the Hamiltonian.  Hence the units of the transport
function are $[1/a_0\hbar^2]$.  The kinetic coefficients, $A_n$, are now given
by:
\begin{equation}
  \label{eq:MultiOrbitalAnFunction}
  A_n = N_{spin}\pi\hbar\int_{-\infty}^{\infty}d\epsilon
  \phi^{xx}(\epsilon)f(\epsilon)f(-\epsilon)(\beta\epsilon)^{n}
\end{equation}
and thus the units of $A_n$ are $[eV/a_0\hbar]$.  Note that if there
is some interaction between spins in the Hamiltonian then all spin sectors
must be included in the Hamiltonian and the factor of $N_{spin}$ is omitted.
The factor of $\hbar$ comes from the final division by $\omega$ in the
expression for the kinetic coefficients (see equation
(\ref{eq:omegaZeroLimit})).  The reason why we get a factor of $\hbar$ here is
due to the fact that this $\omega$ comes from the relation between the
electric field and the vector potential and there of course $\omega$ is not
measured in units of energy. The transport parameters and their units are now
given by:
\begin{alignat}{3}
\rho &= \frac{\kb T}{e^2}\frac{1}{A_{0}}  & \qquad
S &= \frac{-\kb}{e}{A_{1} \over A_{0}} & \qquad 
\kappa &= k_{B}\left(A_{2}-{A_{1}^2 \over A_{0}}\right)
\nonumber\\
&\left[\frac{a_0\hbar}{e^2}\right]  & \qquad
&\left[\frac{\kb}{e}\right]  & \qquad
&\left[\frac{\kb eV}{a_0\hbar}\right].
\end{alignat}
The numerical values of these units are shown in table \ref{tab:units}. 
\begin{table}[htb]
\label{tab:units}
\centering
\begin{tabular}{|c|c|}\hline
Quantity & Natural unit \\\hline\hline
electrical resistivity & 0.2173985 $\mu\Omega m$ \\\hline
thermopower            & 86.173857 $\mu V /K$  \\\hline
thermal conductivity   & 396.38665 $W/(mK)$    \\\hline
\end{tabular}
\caption{The natural units for the transport parameters.}
\end{table} 

In addition to the transport function, the program computes a few other
auxiliary quantities which will be described in the following sections.

\section{Input and overview}
We will now give a detailed prescription of what the input for the calculation
is and an overview of the different components of the computational program
used to implement the calculation.  In figure \ref{fig:ProgFlow} we have
displayed a schematic overview of the program that computes the transport
function.  As we can see there are quite a lot of data files needed for the
computation.  We will now describe each of the data files and the variables
they contain.
\begin{figure}[b]
\centering
\epsfig{file=Figures-6/flowchart.ps,angle=-90,width=0.99\linewidth}
\caption{A schematic overview of the computation of the transport function.
  The statements within the dashed line boxes are only executed in mode {\it
    bare}.  The statements within the dotted line boxes are only used for
  control purposes.  If the {\it dat} files connected to program statements
  with a dashed line do not exist default settings are used.}
\label{fig:ProgFlow}
\end{figure}
\clearpage                   

\subsection{Hamiltonian data (ham.dat), subroutine: ini\_ham}
This file contains the Hamiltonian matrix at each irreducible mesh point in
the Brillouin zone.  To facilitate and standardize the rotations of the
Hamiltonian within the Brillouin zone the Hamiltonian is always assumed to be
written in a spherical harmonics basis.  This means that if one wants to
consider a model Hamiltonian for only the $t_{2g}$ levels of a crystal field
split d-level one must read in the whole Hamiltonian for the d-level.  
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Number of irreducible k-points & nkp  \\\hline
  Dimension of Hamiltonian      & ndim \\\hline
  Hamiltonian                   & hamf (1:ndim,1:ndim,1:nkp) \\\hline
\end{MyProgTab}
\caption{Variables read in from file ham.dat}
\end{table}
All the computations within the code are done in electron volts, but the
default units on the Hamiltonian are Rydbergs and the unit conversion is done
within the program.  If the Hamiltonian is already in electron volts the
conversion can be suppressed by the command line switch {\bf elv}.  The program
checks if the Hamiltonian is Hermitian and will exit if that check fails.

\subsection{Overlap matrix data (olp.dat), subroutine: ini\_ham}
This file contains the overlap matrix at each irreducible mesh point in the
Brillouin zone.  The dimension of this matrix should be equal to the dimension
of the Hamiltonian and the number of irreducible k-points should also match
with the Hamiltonian data.  This is checked by the code and the program exits
if the check fails.  If one is working in an orthogonal basis the overlap
matrix is redundant and in that case the program can be run in mode {\bf orth}
in which case the file olp.dat is not read.  As with the Hamiltonian the
program checks for hermiticity.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Number of irreducible k-points & nkp  \\\hline
  Dimension of overlap matrix   & ndim \\\hline
  Overlap matrix                & olap (1:ndim,1:ndim,1:nkp) \\\hline
\end{MyProgTab}
\caption{Variables read in from file olp.dat}
\end{table}

\subsection{Correlation index data (cix.dat), subroutine: ini\_cix}
\label{ssec:cix.dat}
The file cix.dat contains the so called correlation index for each of the
orbitals in the Hamiltonian.  This index can be used when the self energy in
diagonal in the spherical harmonics basis.  Then we can simply read in the
different diagonal values of the self energy at each frequency point and apply
them to the appropriate orbital.  Thus the orbitals will be assigned an index
from zero to {\bf ncorr} where zero denotes an uncorrelated orbital and a
nonzero index corresponds to data in the self energy file.  If the file
cix.dat does not exist it is assumed that all orbitals are correlated with the
same correlation index, equal to 1.  If the self energy is not diagonal in the
spherical harmonics basis then we must of course read in the whole self energy
matrix (see \article\ref{sec:progSigma}).
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Correlation index & idnsu (1:ndim) \\\hline
\end{MyProgTab}
\caption{Variables read in from file cix.dat}
\end{table}

\subsection{Group operations (grp.dat and wig.dat), subroutine: ini\_grp}
The file grp.dat contains information about the space group operations of the
crystal.  Here we read in both the number of atoms and the number of group
operations that exist.  In principle the maximum number of symmetry operations
is 48 but in low symmetry crystals they could be less and if one is using a
reducible (i.e. not minimal) unit cell they could be more.  To each of the
symmetry operation corresponds a group operation element but since the maximum
number of group elements is 48, two or more symmetry operations could
correspond to a given group element.  In addition to knowing the
correspondence between the symmetry operations and the group elements we
must know how the atoms move under the action of the symmetry operation if
there are more than one atom of a given type in the unit cell.  There is also
the possibility that the symmetry operation takes the atom out of the unit
cell and then it must be translated back by a lattice vector back to the
original unit cell.  Since we are using the Hamiltonian in the Brillouin zone
these translations must be known because they introduce phase factors into
the equation for the transformation of the Hamiltonian. 
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Number of symmetry operations & ngroup  \\\hline
  Number of atoms in unit cell  & natom \\\hline
  Symmetry operation to group element pointer & ikov (1:ngroup) \\\hline
  Atom to atom mapping for each symmetry operation 
  & igatom (1:natom,1:ngroup) \\\hline
  Back to cell translation vectors for each atom and symmetry operation 
  & rdc (1:3,1:natom,1:ngroup) \\\hline
\end{MyProgTab}
\caption{Variables read in from file grp.dat}
\end{table}

The file wig.dat contains Wigner matrices for spherical harmonics angular
momenta up $l = 3$.  This is needed to both rotate the Hamiltonian and the
overlap matrix and to create the vector rotation matrices for each group
operation.  Here we read in the number of group operations (should always
be 48) and the maximum angular momentum (should be 3).  The Wigner matrix
is stored as a vector in the computer and the mapping from the angular 
momentum labels, $(l,m_1,m_2)$, to the vector component, $i$, is given by:
$i = l(2l-1)(2l+1)/3+(2l+1)(l+m_2)+(l+m_1+1)$
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Number of group operations & ng $\equiv$ 48  \\\hline
  Maximum angular momentum  & lmax $\equiv$ 3 \\\hline
  Wigner matrices  & wigmat (1:imax,1:48) \\\hline
\end{MyProgTab}
\caption{Variables read in from file wig.dat}
\end{table}

In addition to reading in the variables discussed above we also initialize
rotation matrices for vectors in cubic representation as well as a pointer
from a symmetry operation to its inverse.  This is needed to compute the
semirotated velocity (\ref{eq:semirotatedvelocity}) which is used in the
computation of the transport function. Here we also compute the displacement
term,
$[(\vec{R}^{g}_{\tau'})^j-(\vec{R}^{g}_{\tau})^j](R^{\gamma})^{-1}_{jx}$,
in the semirotated velocity.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Cubic rotation matrices & rotmat (1:3,1:3,1:ngroup) \\\hline
  Operation to inverse operation pointer  & groupinv (1:ngroup) \\\hline
  Velocity displacement term  & drmat (1:3,1:natom,1:natom,1:ngroup) \\\hline
\end{MyProgTab}
\caption{Additional variables initialized in symmetry group phase.}
\end{table}

\subsection{Brillouin zone data (bz.dat), subroutine: ini\_bz}
The file bz.dat contains the data on the Brillouin zone mesh.  Here we read in
various information about the lattice and the Brillouin zone such as the
lattice constant, the mesh divisions, the irreducible k-points and tetrahedra
as well as information about the group transformations that connect each
k-point to its irreducible representative.
\begin{table}[htb]
  \label{tab:progBZ}
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Lattice constant $[a_0]$ & latcon \\\hline
  Number of irreducible k-points & nkp \\\hline
  Number of irreducible tetrahedra & ntet \\\hline
  Total number of k-points & nktot \\\hline
  Number of k-points along each direction & ndiv (1:3) \\\hline
  Coordinate array for irreducible k-points & bk (1:3,1:nkp) \\\hline
  Integration weight array for irreducible k-points & wk (1:nkp) \\\hline
  Pointer from each irreducible tetrahedron to its multiplicity, 
  irreducible k-point corners and type & itt (1:6,1:ntet) \\\hline
  kcut matrices for all tetrahedron types (see figure \ref{fig:CubeSplit}) & 
  kcut (1:3,1:4,1:6) \\\hline
  Pointer from irreducible k-point to the four nearest irreducible neighbors
  in each direction & nbr (1:3,1:4,1:nkp) \\\hline
  Pointer from irreducible k-point to the group operation that rotates the
  irreducible neighbor to the real neighbor & grop (1:3,1:4,1:nkp) \\\hline
  Bravais lattice vectors, multiplied by the division along the 
  reciprocal direction in Brillouin zone (see eq. (\ref{eq:Amatrix}))& 
  rvec (1:3,1:3) \\\hline
  Brillouin zone minicell vectors 
  (i.e. $\vec{b}_j/N_j$ from eq. (\ref{eq:BZmeshpoints})) & 
  mvec (1:3,1:3) \\\hline
  Pointer from each mesh points to its irreducible representative & 
  kirrptr \mbox{(0:ndiv(1)-1,0:ndiv(2)-1,0:ndiv(3)-1)} \\\hline
  Pointer from each mesh point to the group operation that takes the
  irreducible representative to that point & 
  gropptr \mbox{(0:ndiv(1)-1,0:ndiv(2)-1,0:ndiv(3)-1)} \\\hline
\end{MyProgTab}
\caption{Variables read in from file bz.dat.}
\end{table}
\begin{figure}[t]
\centering
\epsfig{file=Figures-6/cubesplit.ps,angle=-90,width=0.95\linewidth}
\caption{One of four different ways of splitting the unit cube into six
  tetrahedra, here the common diagonal is chosen from $(0,1,0)$ to $(1,0,1)$.
  The coordinates of a given corner $i=1,\ldots,4$ of each tetrahedron
  $j=1,\ldots,6$ are given by the pointer kcut($\mu$,$i$,$j$)$\in\{0,1\}$
  where $\mu=1,2,3$.}
\label{fig:CubeSplit}
\end{figure}              
\begin{figure}[t]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-6/neighbor.eps,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-6/back2cell.eps,width=0.45\linewidth}} 
\end{tabular}
\caption{Left figure: An illustration of the information encoded in the 
  arrays {\bf nbr} and {\bf grop} along the 3rd ($k_z$) direction.  Right
  figure: An illustration of the back-to-cell vectors {\bf rdc} that take an
  atom that is rotated by the transformation igr out of the unit cell back to
  its equivalent position in the unit cell.} 
\label{fig:Neighbor}
\end{figure}       
In addition to the variables in table \ref{tab:progBZ} that are read in from
the file bz.dat there are two more variables that are initialized in this
phase of the program.  One of them is the total number of tetrahedra which is
computed from the variable {\it itt} but should of course be equal to $6\times
ndiv(1)\times ndiv(2)\times ndiv(3)$.  The other one is a pointer from each
corner of every tetrahedron to its irreducible representative.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Total number of tetrahedra & mtet \\\hline
  Pointer from each tetrahedron corner to its irreducible representative &
  coridx (1:4,1:6,1:ndiv(1),1:ndiv(2),1:ndiv(3)) \\\hline
\end{MyProgTab}
\caption{Additional variables initialized in Brillouin zone phase.}
\end{table}

\subsection{Auxiliary index data (idx.dat), subroutine: ini\_idx} 
The file idx.dat contains pointers from the computer representation of the
physical quantities (Hamiltonian etc.) as array to the corresponding atoms and
angular momenta.  Here we also get information about how many sorts of atoms
there are, a pointer from each atom to its sort and what is the maximum
angular momentum for a given atomic sort.  We also have a logical array that
indicates if some angular momenta are omitted from the Hamiltonian.  This is
useful if we want to consider simplified models where certain angular momenta
have been ignored, usually because they are of much to low energy to be of
interest.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Number of atoms in unit cell & natom \\\hline
  Number of atomic sorts & nsort \\\hline
  Pointer from each atom to its sort & is (1:natom) \\\hline
  Pointer from each sort to its maximum angular momentum 
  & ilsb (1:nsort) \\\hline
  Logical array indicating whether a given angular momentum is present for a
  given sort of atoms & ltpan (1:lmax,1:natom) \\\hline
  A pointer from each orbital state of every to its position in the
  Hamiltonian matrix & lmtind (1:(lmax+1)$^2$,1:natom) \\\hline
\end{MyProgTab}
\caption{Variables read in from file idx.dat.}
\end{table}
                                                                         
\subsection{Current combination data (xyz.dat), subroutine: ini\_xyz}
This data is actually now computed inside the subroutine ini\_xyz but it was
originally read from the file xyz.dat.  This data is used to properly
symmetrize the quantity that is being calculated.  In our case we want to
compute the transport function $\phi^{xx}$ but for instance in a cubic crystal
we would instead want compute $(\phi^{xx}+\phi^{yy}+\phi^{zz})/3$ and what
this routine does is to properly generate the symmetrization factors,
$S_{ij}$, for the case $\phi^{xx}$, i.e.
\begin{equation}
  \phi^{xx} \longrightarrow \sum_{ij}S_{ij}\phi^{ij}.
\end{equation}
If another case is needed, say $\phi^{xy}$ then the routine would have to be
changed appropriately.  In addition to the symmetrization parameter we also
create the "rotation" matrices for all the symmetry operations, expressed in a
cubic coordinate system.  These rotation matrices are used to create the
symmetrization parameters and for control purposes they are also printed to
the file {\bf rot.out}.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Cubic rotation matrices & rot\_car (1:3,1:3) \\\hline
  Transport function symmetrization parameters & curxyz (1:3,1:3) \\\hline
\end{MyProgTab}
\caption{Variables generated in subroutine ini\_xyz.}
\end{table}

\subsection{Bare function data (shft.dat), subroutine: ini\_shift}
The file shft.dat contains data that is used in the computation of the bare
transport function and the bare density of states (see
\article\ref{sec:smallscattering}).  When we compute the bare quantities from
the LMTO Hamiltonian we might want to apply a shift to some of the orbitals to
account for, say, a Hartree shift in certain orbitals.  In addition we can can
envision using a different value of the scattering when we compute the
lifetimes that go into the computation of the bare transport function.  This
corresponds to taking $\Sigma_0$ from equation (\ref{eq:smallscatteringSigma})
to be different from the identity, which is the default setting.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  The real part of the diagonal values of the self energy $\Sigma$ from
  equation (\ref{eq:smallscatteringSigma}) & levelshift (0:ncorr) \\\hline
  The imaginary part of the diagonal values of the $\Sigma_0$ matrix from
  equation (\ref{eq:smallscatteringSigma}) & gammaratio (0:ncorr) \\\hline
\end{MyProgTab}
\caption{Variables read in from file shft.dat.}
\end{table}


\subsection{External parameter data (par.dat), subroutine: ini\_par}
The file par.dat contains several parameter, the only one of which is used
here is the chemical potential.  The rest of the parameters that are read in
are not needed for the calculation as it is done now, they are remains of the
original DMFT code from the Russian group of V.~Anisimov.
\begin{table}[htb]
\centering
\begin{MyProgTab}\hline
  Quantity         &  Variable \\\hline\hline
  Chemical potential [eV] & efmu \\\hline
  Inverse temperature [1/eV] & beta \\\hline
  Hubbard U parameter [eV]   & upar \\\hline
  Density of LDA electrons   & nlda \\\hline
  Degeneracy of the correlated band & degen \\\hline
\end{MyProgTab}
\caption{Variables read in from file par.dat.}
\end{table}

\subsection{Self energy (STDIN):}
\label{sec:progSigma}
The retarded self energy is needed in a region around the Fermi surface.  A
rough estimate is that one needs to have data of about 20$\kb T$ around the
Fermi surface in order to get a good estimate of the transport parameters.
Depending on the structure of the self energy we have two choices as to how it
is read in from the standard input.  If the self energy is diagonal in the
spherical harmonics basis we can simply read in all the different diagonal
values and assign them to the proper orbitals (see
\article\ref{ssec:cix.dat}).  In that case we keep the self energy in tabular
form in the datafile that is read from standard input.  If the self energy is
not diagonal in the spherical harmonic basis then we need to read in the whole
matrix at each frequency and this is done by using the mode {\bf mat} and then
we read in the frequency first and then each matrix element separately.

It is necessary to make sure that the self energy is actually retarded and for
that we must make sure that the eigenvalues of the matrix have negative
imaginary part and this is checked in the program and the program will exit if
the check fails.  In addition we add some artificial broadening to the self
energy.  This is done in order to regularize the integral and decrease noise
from weakly correlated bands.  It should of course be emphasized that if the
bands around the Fermi-surface are very weakly correlated then this artificial
broadening is going to dominate the behavior of the transport function and the
bare transport function is probably a better choice for the computation of the 
transport parameters.

\section{Computational modes}
The program can be run in several different computational modes depending on
what information is desired.  The self energy is needed for all these modes
even though only the external frequency is used like in the bare mode.
\begin{enumerate}
\item {\bf bare:}  Program computes the bare dos and transport functions.
\item {\bf dos:} Program computes the full\footnote{Here {\it full} refers to
    quantities computed using a realistic self energy, as opposed to the
    simpler {\it bare} quantities (see equation
    (\ref{eq:BareTransportFunction})).} spectral density.
\item {\bf phi:}   Program computes the full transport function.
\item {\bf eig:}   Program computes eigenvalue characteristics (not useful).
\item {\bf orth:}  Program assumes orthonormal basis.
\item {\bf elv:}   Program assumes Hamiltonian is given in electron Volts.
\item {\bf mat:}   Program reads in the whole self energy matrix.
\end{enumerate}

\section{Data written by the program}
The program writes out several files that contain information about the
progress of the program and the results of computations.

\begin{enumerate}
\item {\bf phi.info} This file contains information about the progress of the
  program.  Most of the information is written out during the preparation
  process, i.e. before we start looping over the frequency data.  It contains
  such things as the dimensions of the Brillouin zone data, information about
  the Hamiltonian and the overlap matrix and information about the CPU time
  taken for various tasks.
\item {\bf bare.out} This file contains the bare quantities
  \
  \begin{itemize}
  \item Column \#1: Frequency
  \item Column \#2: Bare density of states
  \item Column \#3: Bare transport function for lifetime = $1/2\pi$
  \item Column \#4: Bare transport function for averaged lifetime
  \item Column \#5: Bare transport function for k-dependent lifetime
  \end{itemize}
\item   {\bf phi.out} This file contains the actual transport function data.
  \begin{itemize}
  \item Column \#1: Frequency
  \item Column \#2: Transport function
  \item Column \#3: Transport function, SSDD part
  \item Column \#4: Transport function, SSOD part
  \item Column \#5: Transport function, OSDD part
  \item Column \#6: Transport function, OSOD part
  \end{itemize}
\item {\bf dos.out} This file contains the density of states data.  
  \begin{itemize}
  \item Column \#1: Frequency
  \item Column \#2: Density of states
  \item Column \#n+2: Density of states for Hamiltonian orbital \#n (not
    necessarily positive)
  \end{itemize}
\item   {\bf rot.out} This file contains the vector rotation matrices for
  each group operation.
\item   {\bf wig.out} This file contains the Wigner matrices written out
  in nice form.
\item   {\bf drmat.out} This file contains the possible values of  the 
  direct lattice vector difference $R^{g}_{\tau'}-R^{g}_{\tau}$ used
  in the rotation of the Hamiltonian.
\end{enumerate}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 



