Band-resolved local density matrices

The DOS mode densityMatrix stores, for every eigenstate, k point and selected atom, the full local density matrix of the state in the muffin-tin sphere in the banddos.hdf file. From these matrices a few lines of Python give projections that are not available as a built-in DOS mode, e.g. weights for arbitrary combinations of spherical harmonics, t2g/eg-resolved DOS, weights in a rotated local frame, j-resolved weights, or orbital occupation matrices and spin moments of selected states. What is stored, the conventions and the frames are described in The band-resolved local density matrix.

Prerequisites:

• A self-consistent density, as for any DOS or band structure calculation.
• FLEUR compiled with HDF5. Without HDF5, only the traces are written (to text files).
• An inp.xml with fleurInputVersion="0.39" or newer.

Setting up the input

The mode is an addition to a DOS or band-structure calculation. It is switched on by the element densityMatrix in bandDOS.

<output dos="T">
   <bandDOS minEnergy="-0.5" maxEnergy="0.5" sigma="0.015">
      <densityMatrix l="spd" lCross="none" symmetrize="auto"/>
   </bandDOS>
</output>

Three choices determine what is stored:

• l: the $l$ channels, any subset of spdf. Usually only the channels that are relevant for the analysis, e.g. l="d" for a transition metal.
• lCross: blocks between different $l$, e.g. lCross="pd" for p-d hybridization. They are only needed for quantities that mix $l$ channels, such as weights of hybrid orbitals. The default none stores no such blocks.
• symmetrize: with the default auto, the matrices are symmetrized for a DOS on a k-point mesh and left unsymmetrized for a band structure on a path. Use symmetrize="F" if you need weights of single orbitals that are mixed by the site symmetry (e.g. $d_{xy}$ vs. $d_{yz}$ at a cubic site). Their sums over the irreducible representations (t2g, eg) are the same in both cases. See the theory page.

The selection of atoms and frames is shared with the other DOS modes:

Running the calculation

Run FLEUR once on the converged density:

fleur_MPI

As for any DOS run, FLEUR stops after the postprocessing. The matrices are calculated together with the standard $l$-resolved DOS. In non-collinear calculations without l_mperp, both spins are computed in a single pass when this mode is active. The results are unchanged, but the mtCharges of both spins then appear in a single out.xml element.

Output

Two lines in the out file report the memory needed for the matrices and a consistency check of the traces against the $l$-resolved weights of the standard DOS:

grep "DensityMatrix DOS" out

 DensityMatrix DOS: storage of the matrices      0.003 GB
 DensityMatrix DOS: max. deviation of the traces from the l-resolved DOS weights:   1.332E-15

The deviation should be at the level of numerical precision.

The data are written to the group DensityMatrix of banddos.hdf:

• DensityMatrix/DOS, DensityMatrix/EV and DensityMatrix/BS hold the traces of the diagonal blocks as ordinary weights DM:<atom><l> (non-collinear: also DM:<atom><l>mx, my, mz), as for every other DOS mode. They can be plotted with the masci-tools and are also written to the text files DensityMatrix.1 (and DensityMatrix.2), with the same layout as Local.1.
• DensityMatrix/matrix/elem_<n>/pair_<l><l'> holds the matrices, one group per stored atom (elem_<n>) and one dataset per block, e.g. pair_dd or pair_pd. DensityMatrix/matrix/eigenvalues holds the eigenvalues in Hartree, not shifted by the Fermi energy.

In Python (h5py) a dataset pair_<l><l'> has the shape (nkpt, nbands, nSpinBlocks, 2l'+1, 2l+1, 2). The last index separates real and imaginary part, the two before it are $m'+l'$ (row) and $m+l$ (column). A complete description is in the reference.

import h5py, numpy as np

f = h5py.File("banddos.hdf", "r")
m = f["DensityMatrix/matrix"]
eig = np.asarray(m["eigenvalues"])                 # (spin, k, band) in Htr
a = np.asarray(m["elem_1/pair_dd"])
rho = a[..., 0] + 1j*a[..., 1]                     # (k, band, spin block, m', m)

# d charge of each state and spin (non-collinear: blocks 0 and 1 are 11 and 22)
nd = np.trace(rho, axis1=-2, axis2=-1).real        # (k, band, spin block)

# real d orbitals from the complex harmonics m=-2..2, same order as the orbcomp mode
h = np.sqrt(0.5)
T = np.array([[h, 0, 0, 0, -h],      # xy
              [0, h, 0, h, 0],       # yz
              [0, h, 0, -h, 0],      # zx
              [h, 0, 0, 0, h],       # x2-y2
              [0, 0, 1, 0, 0]])      # z2
w = np.einsum("am,kbmn,an->kba", T, rho[:, :, 0], np.conj(T)).real   # (k, band, orbital), spin block 0
t2g = w[..., :3].sum(-1)
eg = w[..., 3:].sum(-1)

Any other unitary transformation T (cubic f harmonics, crystal-field states, $|j m_j\rangle$ states built from all four spin blocks with Clebsch-Gordan coefficients) is used in the same way.

Example: bcc Fe

A ferromagnetic bcc Fe calculation shows the mode on a cubic site. The input for inpgen:

bcc Fe, ferromagnetic
&lattice latsys='bcc' a=2.87 a0=1.88973 /
1
26 0.0 0.0 0.0
&kpt div1=10 div2=10 div3=10 /

After inpgen -f and the self-consistency cycle (fleur_MPI), add the element to bandDOS in inp.xml, set output/@dos="T" and run FLEUR again:

<bandDOS minEnergy="-.5*Htr" maxEnergy="0.5*Htr" sigma="0.015" storeEVData="T">
   <densityMatrix l="spd"/>
</bandDOS>

The k-point mesh is used, so the matrices are symmetrized. Summing the matrices of the occupied states with the k-point weights (/kpts/weights) gives the d occupation matrix of each spin. Transformed to the real d orbitals with T from above, it is diagonal, as required by the cubic site symmetry. The occupations of t2g ($d_{xy}, d_{yz}, d_{zx}$) and eg ($d_{x^2-y^2}, d_{z^2}$) orbitals are

The traces agree with the d charges in the mtCharges of out.xml (4.265 and 1.822) up to the small difference from occupying the states with a step function at the Fermi energy instead of the Fermi smearing of FLEUR.

The regression tests basic/CuDM (all blocks of s, p, d) and noco/FeBccNocoDM (non-collinear, unsymmetrized, all four spin blocks) are further examples.

Hints

• Size. Memory and file size grow with the number of spin blocks, bands, k points, blocks and atoms. Check the storage line in out before using large k-point sets, and restrict l, lCross and the atoms to what you need.
• Per-orbital weights. Symmetrized matrices average the weights of orbitals within an irreducible representation of the site symmetry state by state. Use symmetrize="F" for single orbitals, and T (or auto on a mesh) for DOS-like sums over the Brillouin zone.
• Non-collinear calculations with symmetry. Prefer symmetrize="F", see the limitations.