Orbital construction and conversion

The first step for a DMFT calculation is to provide the necessary input based on a DFT calculation. We will not review how to do the DFT calculation here in this documentation, but refer the user to the documentation and tutorials that come with the actual DFT package. Here, we will describe how to use output created by Wien2k, as well as how to use the light-weight general interface.

Interface with Wien2k

We assume that the user has obtained a self-consistent solution of the Kohn-Sham equations. We further have to require that the user is familiar with the main in/output files of Wien2k, and how to run the DFT code.

Conversion for the DMFT self-consistency cycle

First, we have to write the necessary quantities into a file that can be processed further by invoking in a shell the command

x lapw2 -almd

We note that any other flag for lapw2, such as -c or -so (for spin-orbit coupling) has to be added also to this line. This creates some files that we need for the Wannier orbital construction.

The orbital construction itself is done by the Fortran program dmftproj. For an extensive manual to this program see TutorialDmftproj.pdf. Here we will only describe the basic steps.

Let us take the compound SrVO3, a commonly used example for DFT+DMFT calculations. The input file for dmftproj looks like

3                ! Nsort
1 1 3            ! Mult(Nsort)
3                ! lmax
complex          ! choice of angular harmonics
1 0 0 0          ! l included for each sort
0 0 0 0          ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0)
cubic            ! choice of angular harmonics
1 1 2 0          ! l included for each sort
0 0 2 0          ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0)
01               !
0                ! SO flag
complex          ! choice of angular harmonics
1 1 0 0          ! l included for each sort
0 0 0 0          ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0)
-0.11 0.14

The first three lines give the number of inequivalent sites, their multiplicity (to be in accordance with the Wien2k struct file) and the maximum orbital quantum number \(l_{max}\). In our case our struct file contains the atoms in the order Sr, V, O.

Next we have to specify for each of the inequivalent sites, whether we want to treat their orbitals as correlated or not. This information is given by the following 3 to 5 lines:

  1. We specify which basis set is used (complex or cubic harmonics).
  2. The four numbers refer to s, p, d, and f electrons, resp. Putting 0 means doing nothing, putting 1 will calculate unnormalized projectors in compliance with the Wien2k definition. The important flag is 2, this means to include these electrons as correlated electrons, and calculate normalized Wannier functions for them. In the example above, you see that only for the vanadium d we set the flag to 2. If you want to do simply a DMFT calculation, then set everything to 0, except one flag 2 for the correlated electrons.
  3. In case you have a irrep splitting of the correlated shell, you can specify here how many irreps you have. You see that we put 2, since eg and t2g symmetries are irreps in this cubic case. If you don’t want to use this splitting, just put 0.
  4. (optional) If you specifies a number different from 0 in above line, you have to tell now, which of the irreps you want to be treated correlated. We want to t2g, and not the eg, so we set 0 for eg and 1 for t2g. Note that the example above is what you need in 99% of the cases when you want to treat only t2g electrons. For eg’s only (e.g. nickelates), you set 10 and 01 in this line.
  5. (optional) If you have specified a correlated shell for this atom, you have to tell if spin-orbit coupling should be taken into account. 0 means no, 1 is yes.

These lines have to be repeated for each inequivalent atom.

The last line gives the energy window, relative to the Fermi energy, that is used for the projective Wannier functions. Note that, in accordance with Wien2k, we give energies in Rydberg units!

After setting up this input file, you run:

dmftproj

Again, adding possible flags like -so for spin-orbit coupling. This program produces the following files (in the following, take case as the standard Wien2k place holder, to be replaced by the actual working directory name):

  • case.ctqmcout and case.symqmc containing projector operators and symmetry operations for orthonormalized Wannier orbitals, respectively.
  • case.parproj and case.sympar containing projector operators and symmetry operations for uncorrelated states, respectively. These files are needed for projected density-of-states or spectral-function calculations in post-processing only.
  • case.oubwin needed for the charge density recalculation in the case of fully self-consistent DFT+DMFT run (see below).

Now we convert these files into an hdf5 file that can be used for the DMFT calculations. For this purpose we use the python module Wien2kConverter. It is initialized as:

from pytriqs.applications.dft.converters.wien2k_converter import *
Converter = Wien2kConverter(filename = case)

The only necessary parameter to this construction is the parameter filename. It has to be the root of the files produces by dmftproj. For our example, the Wien2k naming convention is that all files are called the same, for instance SrVO3.*, so you would give filename = “SrVO3”. The constructor opens an hdf5 archive, named case.h5, where all the data is stored. For other parameters of the constructor please visit the Converters section of the reference manual.

After initializing the interface module, we can now convert the input text files to the hdf5 archive by:

Converter.convert_dft_input()

This reads all the data, and stores it in the file case.h5. In this step, the files case.ctqmcout and case.symqmc have to be present in the working directory.

After this step, all the necessary information for the DMFT loop is stored in the hdf5 archive, where the string variable Converter.hdf_filename gives the file name of the archive.

At this point you should use the method dos_wannier_basis contained in the module SumkDFTTools to check the density of states of the Wannier orbitals (see Tools for analysis).

You have now everything for performing a DMFT calculation, and you can proceed with the section on single-shot DFT+DMFT calculations.

Data for post-processing

In case you want to do post-processing of your data using the module SumkDFTTools, some more files have to be converted to the hdf5 archive. For instance, for calculating the partial density of states or partial charges consistent with the definition of Wien2k, you have to invoke:

Converter.convert_parproj_input()

This reads and converts the files case.parproj and case.sympar.

If you want to plot band structures, one has to do the following. First, one has to do the Wien2k calculation on the given \(\mathbf{k}\)-path, and run dmftproj on that path:

x lapw1 -band
x lapw2 -band -almd
dmftproj -band

Again, maybe with the optional additional extra flags according to Wien2k. Now we use a routine of the converter module allows to read and convert the input for SumkDFTTools:

Converter.convert_bands_input()

After having converted this input, you can further proceed with the Tools for analysis. For more options on the converter module, please have a look at the Converters section of the reference manual.

Data for transport calculations

For the transport calculations, the situation is a bit more involved, since we need also the optics package of Wien2k. Please look at the section on Transport calculations to see how to do the necessary steps, including the conversion.

A general H(k)

In addition to the more complicated Wien2k converter, DFTTools contains also a light converter. It takes only one inputfile, and creates the necessary hdf outputfile for the DMFT calculation. The header of this input file has a defined format, an example is the following:

64          ! number of k-points
1.0         ! Electron density
2           ! number of total atomic shells
1 1 2 5     ! iatom, isort, l, dimension
2 2 1 3     ! iatom, isort, l, dimension
1           ! number of correlated shells
1 1 2 5 0 0 ! iatom, isort, l, dimension, SO, irep
1 5         ! # of ireps, dimension of irep

The lines of this header define

  1. Number of \(\mathbf{k}\)-points used in the calculation

  2. Electron density for setting the chemical potential

  3. Number of total atomic shells in the hamiltonian matrix. In short, this gives the number of lines described in the following. IN the example file give above this number is 2.

  4. The next line(s) contain four numbers each: index of the atom, index of the equivalent shell, \(l\) quantum number, dimension of this shell. Repeat this line for each atomic shell, the number of the shells is given in the previous line.

    In the example input file given above, we have two inequivalent atomic shells, one on atom number 1 with a full d-shell (dimension 5), and one on atom number 2 with one p-shell (dimension 3).

    Other examples for these lines are:

    1. Full d-shell in a material with only one correlated atom in the unit cell (e.g. SrVO3). One line is sufficient and the numbers are 1 1 2 5.
    2. Full d-shell in a material with two equivalent atoms in the unit cell (e.g. FeSe): You need two lines, one for each equivalent atom. First line is 1 1 2 5, and the second line is 2 1 2 5. The only difference is the first number, which tells on which atom the shell is located. The second number is the same in both lines, meaning that both atoms are equivalent.
    3. t2g orbitals on two non-equivalent atoms in the unit cell: Two lines again. First line is 1 1 2 3, second line 2 2 2 3. The difference to the case above is that now also the second number differs. Therefore, the two shells are treated independently in the calculation.
    4. d-p Hamiltonian in a system with two equivalent atoms each in the unit cell (e.g. FeSe has two Fe and two Se in the unit cell). You need for lines. First line 1 1 2 5, second line 2 1 2 5. These two lines specify Fe as in the case above. For the p orbitals you need line three as 3 2 1 3 and line four as 4 2 1 3. We have 4 atoms, since the first number runs from 1 to 4, but only two inequivalent atoms, since the second number runs only form 1 to 2.

    Note that the total dimension of the hamiltonian matrices that are read in is the sum of all shell dimensions that you specified. For example number 4 given above we have a dimension of 5+5+3+3=16. It is important that the order of the shells that you give here must be the same as the order of the orbitals in the hamiltonian matrix. In the last example case above the code assumes that matrix index 1 to 5 belongs to the first d shell, 6 to 10 to the second, 11 to 13 to the first p shell, and 14 to 16 the second p shell.

  5. Number of correlated shells in the hamiltonian matrix, in the same spirit as line 3.

  6. The next line(s) contain six numbers: index of the atom, index of the equivalent shell, \(l\) quantum number, dimension of the correlated shells, a spin-orbit parameter, and another parameter defining interactions. Note that the latter two parameters are not used at the moment in the code, and only kept for compatibility reasons. In our example file we use only the d-shell as correlated, that is why we have only one line here.

  7. The last line contains several numbers: the number of irreducible representations, and then the dimensions of the irreps. One possibility is as the example above, another one would be 2 2 3. This would mean, 2 irreps (eg and t2g), of dimension 2 and 3, resp.

After these header lines, the file has to contain the Hamiltonian matrix in orbital space. The standard convention is that you give for each \(\mathbf{k}\)-point first the matrix of the real part, then the matrix of the imaginary part, and then move on to the next \(\mathbf{k}\)-point.

The converter itself is used as:

from pytriqs.applications.dft.converters.hk_converter import *
Converter = HkConverter(filename = hkinputfile)
Converter.convert_dft_input()

where hkinputfile is the name of the input file described above. This produces the hdf file that you need for a DMFT calculation.

For more options of this converter, have a look at the Converters section of the reference manual.

Wannier90 Converter

Using this converter it is possible to convert the output of wannier90 Maximally Localized Wannier Functions (MLWF) and create a HDF5 archive suitable for one-shot DMFT calculations with the SumkDFT class.

The user must supply two files in order to run the Wannier90 Converter:

  1. The file seedname_hr.dat, which contains the DFT Hamiltonian in the MLWF basis calculated through wannier90 with hr_plot = true (please refer to the wannier90 documentation).
  2. A file named seedname.inp, which contains the required information about the \(\mathbf{k}\)-point mesh, the electron density, the correlated shell structure, ... (see below).

Here and in the following, the keyword seedname should always be intended as a placeholder for the actual prefix chosen by the user when creating the input for wannier90. Once these two files are available, one can use the converter as follows:

from pytriqs.applications.dft.converters import Wannier90Converter
Converter = Wannier90Converter(seedname='seedname')
Converter.convert_dft_input()

The converter input seedname.inp is a simple text file with the following format:

 0  6 4 6
8.0
 4
 0  0  2  3  0  0
 1  0  2  3  0  0
 2  0  2  3  0  0
 3  0  2  3  0  0

The example shows the input for the perovskite crystal of LaVO3 in the room-temperature Pnma symmetry. The unit cell contains four symmetry-equivalent correlated sites (the V atoms) and the total number of electrons per unit cell is 8 (see second line). The first line specifies how to generate the \(\mathbf{k}\)-point mesh that will be used to obtain \(H(\mathbf{k})\) by Fourier transforming \(H(\mathbf{R})\). Currently implemented options are:

  • \(\Gamma\)-centered uniform grid with dimensions \(n_{k_x} \times n_{k_y} \times n_{k_z}\); specify 0 followed by the three grid dimensions, like in the example above
  • \(\Gamma\)-centered uniform grid with dimensions automatically determined by the converter (from the number of \(\mathbf{R}\) vectors found in seedname_hr.dat); just specify -1

Inside seedname.inp, it is crucial to correctly specify the correlated shell structure, which depends on the contents of the wannier90 output seedname_hr.dat and on the order of the MLWFs contained in it.

The number of MLWFs must be equal to, or greater than the total number of correlated orbitals (i.e., the sum of all dim in seedname.inp). If the converter finds fewer MLWFs inside seedname_hr.dat, then it stops with an error; if it finds more MLWFs, then it assumes that the additional MLWFs correspond to uncorrelated orbitals (e.g., the O-2p shells). When reading the hoppings \(\langle w_i | H(\mathbf{R}) | w_j \rangle\) (where \(w_i\) is the \(i\)-th MLWF), the converter also assumes that the first indices correspond to the correlated shells (in our example, the V-t2g shells). Therefore, the MLWFs corresponding to the uncorrelated shells (if present) must be listed after those of the correlated shells. With the wannier90 code, this can be achieved this by listing the projections for the uncorrelated shells after those for the correlated shells. In our Pnma-LaVO3 example, for instance, we could use:

Begin Projections
 V:l=2,mr=2,3,5:z=0,0,1:x=-1,1,0
 O:l=1:mr=1,2,3:z=0,0,1:x=-1,1,0
End Projections

where the x=-1,1,0 option indicates that the V–O bonds in the octahedra are rotated by (approximatively) 45 degrees with respect to the axes of the Pbnm cell.

The converter will analyze the matrix elements of the local Hamiltonian to find the symmetry matrices rot_mat needed for the global-to-local transformation of the basis set for correlated orbitals (see section hdf5 structure). The matrices are obtained by finding the unitary transformations that diagonalize \(\langle w_i | H_I(\mathbf{R}=0,0,0) | w_j \rangle\), where \(I\) runs over the correlated shells and i,j belong to the same shell (more details elsewhere...). If two correlated shells are defined as equivalent in seedname.inp, then the corresponding eigenvalues have to match within a threshold of 10-5, otherwise the converter will produce an error/warning. If this happens, please carefully check your data in seedname_hr.dat. This method might fail in non-trivial cases (i.e., more than one correlated shell is present) when there are some degenerate eigenvalues: so far tests have not shown any issue, but one must be careful in those cases (the converter will print a warning message).

The current implementation of the Wannier90 Converter has some limitations:

  • Since wannier90 does not make use of symmetries (symmetry-reduction of the \(\mathbf{k}\)-point grid is not possible), the converter always sets symm_op=0 (see the hdf5 structure section).
  • No charge self-consistency possible at the moment.
  • Calculations with spin-orbit (SO=1) are not supported.
  • The spin-polarized case (SP=1) is not yet tested.
  • The post-processing routines in the module SumkDFTTools were not tested with this converter.
  • proj_mat_all are not used, so there are no projectors onto the uncorrelated orbitals for now.

MPI issues

The interface packages are written such that all the file operations are done only on the master node. In general, the philosophy of the package is that whenever you read in something from the archive yourself, you have to manually broadcast it to the nodes. An exception to this rule is when you use routines from SumkDFT or SumkDFTTools, where the broadcasting is done for you.

Interfaces to other packages

Because of the modular structure, it is straight forward to extend the TRIQS package in order to work with other band-structure codes. The only necessary requirement is that the interface module produces an hdf5 archive, that stores all the data in the specified form. For the details of what data is stored in detail, see the hdf5 structure part of the reference manual.