Home What is? Install Overview CTQMC MnO FeSe SrVO3 LaVO3 Cerium Sr2IrO4 Ta2NiSe5 phonon

Tutorial 7: Ta2NiSe5, irreducible representations, Fermi surface, non-local Hartree-Fock, and optical conductivity

This tutorial is on Ta2NiSe5, which was proposed to be excitonic insulator, as it undergoes q=0 structural phase transition, opening a semiconducting gap of unusual shape.

We will first compute single-site DMFT in both the high-temperature orthorombic phase, and the low temperature monoclinic phase. We will than show that the monoclinic distortion opens a small semiconducting gap within eDMFT (notice that Ta2NiSe5 is metallic both phases within LDA), which is however smaller than in experiment. We will use structural relaxation to prove that the monoclinic distortion is unstable within single-site eDMFT, and the system relaxes into high temperature orthorombic phase during relaxation. We will than identify an electronic symmetry breaking of B3g type, which enhances the gap size, and makes the monoclinic structural distortion stable. This electronic instability is here treated by non-local Hartree-Fock approximation between the trimer of Ta-Ni-Ta atoms. We note that in principle such trimerization could be be treated from ab-initio by cluster-DMFT calculation, where Ta-Ni-Ta trimer would be the quantum impurity. In this case, we would likely need only the on-site Coulomb repulsion, and the non-local Ta-Ni-Ta self-energy would be generated automatically. However, such quantum calculation are very expensive, and currently suffer from fermionic sign problem. We will therefore add the non-local self-energy using exchange interaction only (Hartree-Fock), which is a reasonable approximation in the ordered state. However, in this treatment we than require non-local Coulomb repulsion between Ta and Ni, here estimated to be half of the on-site Coulomb repulsion.

Crystal structures and LDA

The high-temperature orthorombic and the low-temperature monoclinic structure, with space group 63 and 15, respectively, have the following cif files:
Ta2NiSe5O.cif Ta2NiSe5M.cif.

We could use cif2struct tool in wien2k to produce the corresponding structure files. However, for later convenience, we will use slightly modified structure files, i.e.,
Ta2NiSe5O.struct Ta2NiSe5O.struct
which were produced with cif2struct, and slightly modified. In this structure files Ta follows Ni in the structure file (was manually reordered), and most importantly, the group operations are reordered, starting with identify operation first, followed by inversion. This is convenient later for buliding the cluster-trimer out of the Ni-Ta atoms, where the precise positions of the atoms in real space are essential, not just their wyckoff position. We mention that this changes are not essential for the calcullation, but we encourage user to use the same convention for easier comparison of the results.

Next we run LDA in wien2k. Create directory with name Ta2NiSe5O and Ta2NiSe5M, copy corresponding structure files into that directory, and execute
init_lapw
and answer the following questions:
Enter reduction in %:   0
Use old or new scheme (o/N) : N
Do you want to accept these radii; .... (a/d/r) : a
specify nn-bondlength factor: 2
Ctrl-X-C
continue with sgroup : c
Ctrl-X-C
continue with symmetry : c
Ctrl-X-C
continue with lstart : c
Eventually specify switches for instgen_lapw (or press ENTER): ENTER
SELECT XCPOT: 5
SELECT ENERGY : -6
Ctrl-X-C
continue with kgen : c
Ctrl-X-C
Ctrl-X-C
NUMBER OF K-POINTS : 100
Shift of k-mesh allowed. Do you want to shift: (0=no, 1=shift): 1
Ctrl-X-C
continue with dstart or execute kgen again or exit (c/e/x) : c
Ctrl-X-C
do you want to perform a spinpolarized calculation : n
Next, run the DFT calculation using the command
run_lapw
[Of course Wien2k can also be run in parallel, and the Wien2k manual discusses this option in detail. Alternatively, one can just increase number of cores using export OMP_NUM_THREADS=10]

For DMFT calculation, we will increase the number of k-poinst to 500. To do that, run
x kgen
and specify 500 k-points with shifted mesh.

Next, rerun wien2k on this k-mesh, i.e.,
run_lapw -NI

eDMFT initialization

Once the DFT calculation is done, we have a good charge density to start a DMFT calculation. The next step is to initialize the DMFT calculation. For this, use the script
init_dmft.py
and choose the following options:
There are 16 atoms in the unit cell:
  1 Ni1
  2 Ni1
  3 Ta1
  4 Ta1
  5 Ta1
  6 Ta1
  7 Se1
  8 Se1
  9 Se2
 10 Se2
 11 Se2
 12 Se2
 13 Se3
 14 Se3
 15 Se3
 16 Se3
Specify correlated atoms (ex: 1-4,7,8): 1-6
[We correlate both Ta and Ni atoms]
You have chosen the following atoms to be correlated:
  1 Ni1
  2 Ni1
  3 Ta1
  4 Ta1
  5 Ta1
  6 Ta1
Do you want to continue; or edit again? (c/e): c
For each atom, specify correlated orbital(s) (ex: d,f):
  1 Ni1: d
  2 Ni1: d
  3 Ta1: d
  4 Ta1: d
  5 Ta1: d
  6 Ta1: d
[correlations of Ni is on 3d and Ta on 5d, hence all are d]
  
You have chosen to apply correlations to the following orbitals:
  1  Ni1-1 d
  2  Ni1-2 d
  3  Ta1-3 d
  4  Ta1-4 d
  5  Ta1-5 d
  6  Ta1-6 d
Do you want to continue; or edit again? (c/e): c
Specify qsplit for each correlated orbital (default = 0):
    Qsplit  Description
------  ------------------------------------------------------------
     2  real harmonics basis, no symmetry, except spin (up=dn)
------  ------------------------------------------------------------
  1  Ni1-1 d: 2
  2  Ni1-2 d: 2
  3  Ta1-3 d: 2
  4  Ta1-4 d: 2
  5  Ta1-5 d: 2
  6  Ta1-6 d: 2
[For now we will just choose cubic spherical harmonics, and we will later optimize them with different tool]
  
You have chosen the following qsplits:
  1  Ni1-1 d: 2
  2  Ni1-2 d: 2
  3  Ta1-3 d: 2
  4  Ta1-4 d: 2
  5  Ta1-5 d: 2
  6  Ta1-6 d: 2
Do you want to continue; or edit again? (c/e): 
Specify projector type (default = 5):
  Projector  Description
------  ------------------------------------------------------------
     5  fixed projector, which is written to projectorw.dat. You can
        generate projectorw.dat with the tool wavef.py
------  ------------------------------------------------------------
> 5

Do you want to continue; or edit again? (c/e): c

Do you want to group any of these orbitals into cluster-DMFT problems? (y/n): n

Enter the correlated problems forming each unique correlated
problem, separated by spaces (ex: 1,3 2,4 5-8): 1-2 3-6
[ This specifies that atom 1 and 2, which are Ni atoms, are equivalent, while 3-6 are also an equivalent group of Ta atoms]
Each set of equivalent correlated problems are listed below:
  1   (Ni11 d) (Ni12 d) are equivalent.
  2   (Ta13 d) (Ta14 d) (Ta15 d) (Ta16 d) are equivalent.

Do you want to continue; or edit again? (c/e): c

Range (in eV) of hybridization taken into account in impurity
problems; default -10.0, 10.0: <ENTER>

Perform calculation on real; or imaginary axis? (r/i): r

Ctrl-X-C
Ctrl-X-C

Local coordinate axis attached to Ni and Ta atoms

In this part of the tutorial, we will optimize the orbitals by changing Ta2NiSe5O.indmfl file. If you want to compare your results with the final form of the file, you can download it here Ta2NiSe5O.indmfl.

First we want to orient the local coordinate axis on Ta and Ni, so that the orbitals on these sites produce mostly diagonal density matrix. To do that, we execute
localaxes.py
and select first Ni (1) and obtain 3x3 rotation matrix:
  0.70710678   0.00000000   0.70710678 
  0.00000000   1.00000000  -0.00000000 
 -0.70710678   0.00000000   0.70710678 
For convenience, we will flip z and y axis, and change y->-y to keep right coordinate system. In this choice the local z-axis on Ni points along the b axis. [ This choice is not essential, but just convenient].
  
  0.70710678   0.00000000   0.70710678  # local x-axis
  0.70710678   0.00000000  -0.70710678  # local y-axis
  0.00000000   1.00000000  -0.00000000  # local z-axis
Next we execute localaxes.py once more, and select atom 3, i.e., first Ta atom, and we get:
  0.70710678  -0.64655917  -0.28628873 
 -0.70710678  -0.64655917  -0.28628873 
 -0.00000000   0.40487341  -0.91437275 
Now we use these rotations to edit Ta2NiSe5O.indmfl file. We identify the specification for the first Ni atom at the line 4, and we set locrot from 0 to -1. The "locrot" switch is used to specify the local rotation axis for each atom. We will use the cartesian 3x3 unitary matrix of rotation to specify local x, y, and z axis as rows of the matrix. Then we paste the 3x3 rotation matrix one line below this line, and we repeat identical procedure for both Ni atoms. We do the corresponding procedure also for all La atoms, so that the header file of Ta2NiSe5O.indmfl will looks like
111 198 1 5                           # hybridization band index nemin and nemax, renormalize for interstitials, projection type
0 0.025 0.025 200 -3.000000 1.000000  # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)
6                                     # number of correlated atoms
1     1  -1                           # iatom, nL, locrot [THIS IS WHERE WE MADE THE FIRST CHANGE]
  2   2   1                           # L, qsplit, cix
  0.70710678   0.00000000   0.70710678  # local x-axis
  0.70710678   0.00000000  -0.70710678  # local y-axis
  0.00000000   1.00000000  -0.00000000  # local z-axis
2     1  -1                           # iatom, nL, locrot
  2   2   2                           # L, qsplit, cix
  0.70710678   0.00000000   0.70710678  # local x-axis
  0.70710678   0.00000000  -0.70710678  # local y-axis
  0.00000000   1.00000000  -0.00000000  # local z-axis
3     1  -1                           # iatom, nL, locrot
  2   2   3                           # L, qsplit, cix
  0.70710678  -0.64655917  -0.28628873
 -0.70710678  -0.64655917  -0.28628873
 -0.00000000   0.40487341  -0.91437275
4     1  -1                           # iatom, nL, locrot
  2   2   4                           # L, qsplit, cix
  0.70710678  -0.64655917  -0.28628873
 -0.70710678  -0.64655917  -0.28628873
 -0.00000000   0.40487341  -0.91437275
5     1  -1                           # iatom, nL, locrot
  2   2   5                           # L, qsplit, cix
  0.70710678  -0.64655917  -0.28628873
 -0.70710678  -0.64655917  -0.28628873
 -0.00000000   0.40487341  -0.91437275
6     1  -1                           # iatom, nL, locrot
  2   2   6                           # L, qsplit, cix
  0.70710678  -0.64655917  -0.28628873
 -0.70710678  -0.64655917  -0.28628873
 -0.00000000   0.40487341  -0.91437275
For detailed description of this file, please refer to earlier tutorials.

Choice of local orbitals for diagonal hybridization

Next we construct self-energy, which will contain only zeros, by executing
szero.py
Notice that we get self-energy on the real axis. Namely, we will perform calculation of the green's function on the real axis, and further optimize orbitals, to make hybridization as diagonal as possible. We execute
x_dmft.py dmft1
to produce density of states, the local Green's function, and density matrix. The latter is found in Ta2NiSe5O.outputdmf1. Towards the end of this file, we see
  
 Full matrix of impurity levels follows
 icix=           1
   -1.43999270     0.00000000       0.00000000     0.00000000       0.00000000     0.00000000      -0.00000000    -0.00000000       0.07727972     0.00000000   
    0.00000000    -0.00000000      -1.48498476     0.00000000      -0.00000000    -0.00000000      -0.00000000    -0.00000000      -0.00000000    -0.00000000   
    0.00000000    -0.00000000      -0.00000000     0.00000000      -1.39615433     0.00000000       0.05311168     0.00000000       0.00000000    -0.00000000   
   -0.00000000     0.00000000      -0.00000000     0.00000000       0.05311168    -0.00000000      -1.39615433    -0.00000000      -0.00000000     0.00000000   
    0.07727972    -0.00000000       0.00000000     0.00000000       0.00000000     0.00000000      -0.00000000    -0.00000000      -1.37291948     0.00000000   
 
 icix=           1 at omega=0
   -2.43920361     0.00000000      -0.00000000     0.00000000      -0.00000000     0.00000000       0.00000000    -0.00000000       0.31811819    -0.00000000   
   -0.00000000    -0.00000000      -4.18346033     0.00000000      -0.00000000    -0.00000000      -0.00000000    -0.00000000       0.00000000    -0.00000000   
   -0.00000000    -0.00000000      -0.00000000     0.00000000      -2.03348922     0.00000000       0.78590852     0.00000000      -0.00000000     0.00000000   
    0.00000000     0.00000000      -0.00000000     0.00000000       0.78590852    -0.00000000      -2.03348922    -0.00000000       0.00000000    -0.00000000   
    0.31811819     0.00000000       0.00000000     0.00000000      -0.00000000    -0.00000000       0.00000000     0.00000000      -1.62348389     0.00000000   
.....
The first matrix prints the impurity levels, and the second the impurity levels with the real part of the hybridization at zero frequency, i.e., Eimp+Delta(omega=0). For metals and small gap semiconductors, the latter should be mostly diagonal for successful solution of the impurity problem by the current Monte Carlo solver. For large gap insulators, the low energy hybridization will dissapear, hence it is better to optimize Eimp only. In most cases, both can be reasonably diagonal with the same set of orbitals.
The order of orbitals is compatible with the Ta2NiSe5O.indmfl input, namely, z2, x2-y2, xz, yz, xy. We notice that xz and yz orbitals mix, and since they are degenerate, their sum and difference will give a diagonal matrix, hence we will change the Ta2NiSe5O.indmfl file, and input xz+yz, followed by xz-yz orbital. For the Ni-atoms, we will hence change the Ta2NiSe5O.indmfl such that we use xz+yz and xz-yz orbital. In addition, there is small mixing of z2 and xy orbital, which can be diagonalized with taking 0.94565163*z2-0.32518147*xy as the new z2 orbital, and 0.94565163*xy+0.32518147*z2 as the new xy orbital. We thus get:
  
#---------------- # Independent components are --------------
'z^2' 'x^2-y^2' 'xz+yz' 'xz-yz' 'xy'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
 0.00000000  0.22993802    0.00000000  0.00000000    0.94565163  0.00000000    0.00000000  0.00000000    0.00000000 -0.22993802  
 0.70710679  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.70710679  0.00000000  
 0.00000000  0.00000000    0.50000000  0.50000000    0.00000000  0.00000000   -0.50000000  0.50000000    0.00000000  0.00000000  
 0.00000000  0.00000000    0.50000000 -0.50000000    0.00000000  0.00000000   -0.50000000 -0.50000000    0.00000000  0.00000000  
 0.00000000 -0.66867669    0.00000000  0.00000000    0.32518147  0.00000000    0.00000000  0.00000000    0.00000000  0.66867669  
We than rerun
x_dmft.py dmft1
and see that indeed the impurity level matrix is basically diagonal, i.e.,
 Full matrix of impurity levels follows
 icix=           1
   -1.48042850     0.00000000       0.00000000     0.00000000      -0.00000000     0.00000000       0.00000000     0.00000000       0.04031057     0.00000000   
    0.00000000    -0.00000000      -1.48498476     0.00000000      -0.00000000     0.00000000      -0.00000000    -0.00000000       0.00000000    -0.00000000   
   -0.00000000    -0.00000000      -0.00000000    -0.00000000      -1.34304266    -0.00000000       0.00000000     0.00000000       0.00000000     0.00000000   
    0.00000000    -0.00000000      -0.00000000     0.00000000       0.00000000    -0.00000000      -1.44926601    -0.00000000       0.00000000    -0.00000000   
    0.04031057    -0.00000000       0.00000000     0.00000000      -0.00000000    -0.00000000       0.00000000     0.00000000      -1.33248368    -0.00000000   
 
 icix=           1 at omega=0
   -2.54859499     0.00000000       0.00000000    -0.00000000      -0.00000000     0.00000000      -0.00000000     0.00000000      -0.00000000    -0.00000000   
    0.00000000     0.00000000      -4.18346033     0.00000000      -0.00000000     0.00000000      -0.00000000    -0.00000000       0.00000000    -0.00000000   
   -0.00000000     0.00000000      -0.00000000    -0.00000000      -1.24758070    -0.00000000       0.00000000     0.00000000       0.00000000    -0.00000000   
   -0.00000000    -0.00000000      -0.00000000     0.00000000       0.00000000    -0.00000000      -2.81939774    -0.00000000      -0.00000000     0.00000000   
   -0.00000000     0.00000000       0.00000000     0.00000000       0.00000000    -0.00000000      -0.00000000    -0.00000000      -1.51409251    -0.00000000   
We can ckeck the partial density of states for all orbitals by plotting Ta2NiSe5O.gc1. The density of states should look like
Of course the number of k-points is not large enough to see a smooth curve, but we will increase the k-points once we have converged self-energy. We notice that eg orbitals do not contribute at low energy, while t2g's do. For convenience, we will reorder the orbitals, so that the orbital that contributes most at the low energy is entered first, followed by the next important orbital, etc. Specifically, we will choose the order xz+yz, xy, xz-yz, z2, x2-y2. This reordering is easy, as each row of the matrix corresponds to an orbital, hence we just exchange the order of rows of the Transformation matrix. We emphasize that this reordering is not necessary for the calculation. It is merely for convenience, so that when running DMFT on imaginary axis, we will know that the first orbital is most important, followed by the second, etc. We then repeat the same combination of orbitals also for the second Ni atom. We arrive at the following entries for the first two atoms
  
1     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'xz+yz' 'xy' 'xz-yz' 'z^2' 'x^2-y^2'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
 0.00000000  0.00000000    0.50000000  0.50000000    0.00000000  0.00000000   -0.50000000  0.50000000    0.00000000  0.00000000
 0.00000000 -0.66867669    0.00000000  0.00000000    0.32518147  0.00000000    0.00000000  0.00000000    0.00000000  0.66867669
 0.00000000  0.00000000    0.50000000 -0.50000000    0.00000000  0.00000000   -0.50000000 -0.50000000    0.00000000  0.00000000
 0.00000000  0.22993802    0.00000000  0.00000000    0.94565163  0.00000000    0.00000000  0.00000000    0.00000000 -0.22993802
 0.70710679  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.70710679  0.00000000
2     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'xz+yz' 'xy' 'xz-yz' 'z^2' 'x^2-y^2'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
 0.00000000  0.00000000    0.50000000  0.50000000    0.00000000  0.00000000   -0.50000000  0.50000000    0.00000000  0.00000000
 0.00000000 -0.66867669    0.00000000  0.00000000    0.32518147  0.00000000    0.00000000  0.00000000    0.00000000  0.66867669
 0.00000000  0.00000000    0.50000000 -0.50000000    0.00000000  0.00000000   -0.50000000 -0.50000000    0.00000000  0.00000000
 0.00000000  0.22993802    0.00000000  0.00000000    0.94565163  0.00000000    0.00000000  0.00000000    0.00000000 -0.22993802
 0.70710679  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.70710679  0.00000000
Next we need to diagonalize also Ta hybridization function. Currently, the hybridization function, corresponding to Ta, in the Ta2NiSe5O.outputdmf1 is
  
 icix=           3
    1.90366052    -0.00000000       0.00000000    -0.00000000      -0.00703653     0.00000000      -0.00703653    -0.00000000       0.03657682    -0.00000000   
    0.00000000     0.00000000       1.99571598     0.00000000      -0.02338245    -0.00000000       0.02338245     0.00000000      -0.00000000     0.00000000   
   -0.00703653    -0.00000000      -0.02338245     0.00000000       1.46810935     0.00000000      -0.03545005     0.00000000       0.02419652    -0.00000000   
   -0.00703653     0.00000000       0.02338245    -0.00000000      -0.03545005    -0.00000000       1.46810935    -0.00000000       0.02419652    -0.00000000   
    0.03657682     0.00000000       0.00000000    -0.00000000       0.02419652     0.00000000       0.02419652     0.00000000       1.30833154     0.00000000   
 
 icix=           3 at omega=0
   10.74050286    -0.00000000      -0.00000000    -0.00000000       1.11134247     0.00000000       1.11134247    -0.00000000       0.40444646    -0.00000000   
   -0.00000000     0.00000000      18.37514373     0.00000000      -1.53761686    -0.00000000       1.53761686     0.00000000       0.00000000     0.00000000   
    1.11134247    -0.00000000      -1.53761686     0.00000000       4.01205491     0.00000000      -1.49417656     0.00000000       0.23770269    -0.00000000   
    1.11134247     0.00000000       1.53761686    -0.00000000      -1.49417656    -0.00000000       4.01205491    -0.00000000       0.23770269    -0.00000000   
    0.40444646     0.00000000       0.00000000    -0.00000000       0.23770269     0.00000000       0.23770269     0.00000000       1.64248086     0.00000000   
and has a fair number of off-diagonal components. We will diagonalize it with find3dRotation.py, described already in LaNiO3 tutorial. We create a small text file (we will call it rot.dat) of python form, and enter the hybridization matrix (at zero frequency) from Ta2NiSe5O.outputdmf1, as well as the current orbital set from Ta2NiSe5O.indmfl. It should have the form
strHc="""
   10.74050286    -0.00000000      -0.00000000    -0.00000000       1.11134247     0.00000000       1.11134247    -0.00000000       0.40444646    -0.00000000
   -0.00000000     0.00000000      18.37514373     0.00000000      -1.53761686    -0.00000000       1.53761686     0.00000000       0.00000000     0.00000000
    1.11134247    -0.00000000      -1.53761686     0.00000000       4.01205491     0.00000000      -1.49417656     0.00000000       0.23770269    -0.00000000
    1.11134247     0.00000000       1.53761686    -0.00000000      -1.49417656    -0.00000000       4.01205491    -0.00000000       0.23770269    -0.00000000
    0.40444646     0.00000000       0.00000000    -0.00000000       0.23770269     0.00000000       0.23770269     0.00000000       1.64248086     0.00000000
"""
strT2C="""
 0.00000000  0.00000000    0.00000000  0.00000000    1.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000
 0.70710679  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.70710679  0.00000000
 0.00000000  0.00000000    0.70710679  0.00000000    0.00000000  0.00000000   -0.70710679  0.00000000    0.00000000  0.00000000
 0.00000000  0.00000000    0.00000000  0.70710679    0.00000000  0.00000000    0.00000000  0.70710679    0.00000000  0.00000000
-0.00000000 -0.70710679    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.70710679
"""
Clearly strHc has been copied from Ta2NiSe5O.outputdmf1 and strT2C from Ta2NiSe5O.indmfl. We then execute
find3dRotation.py rot.dat
and obtain the following set of orbitals:
  
 0.00000000 -0.66323970   -0.17313727 -0.17313727    0.01789324  0.00000000    0.17313727 -0.17313727    0.00000000  0.66323970  # xy
-0.00000000  0.24274639   -0.46008182 -0.46008182    0.18827449  0.00000000    0.46008182 -0.46008182   -0.00000000 -0.24274639  # -xz-yz
 0.11471345 -0.00000000    0.49337655 -0.49337655    0.00000000  0.00000000   -0.49337655 -0.49337655    0.11471345  0.00000000  # xz-yz
 0.00000000 -0.03445729    0.09136855  0.09136855    0.98195343  0.00000000   -0.09136855  0.09136855    0.00000000  0.03445729  # z^2
 0.69773981  0.00000000   -0.08111466  0.08111466   -0.00000000  0.00000000    0.08111466  0.08111466    0.69773981 -0.00000000  # x^2-y^2
We added comment at the end of the lines to emphasize the dominant orbital. In this case, the orbitals are already order in desired order, i.e, from most important to least important, and we will keep this order. We will change Ta2NiSe5O.indmfl file for each Ta-atom to have the appropriate form
  
#---------------- # Independent components are --------------
'xy' 'xz+yz' 'xz-yz' 'z2' 'x2-y2'
#---------------- # Sigind follows --------------------------
 6  0  0  0  0
 0  7  0  0  0
 0  0  8  0  0
 0  0  0  9  0
 0  0  0  0 10
#---------------- # Transformation matrix follows -----------
 0.00000000 -0.66323970   -0.17313727 -0.17313727    0.01789324  0.00000000    0.17313727 -0.17313727    0.00000000  0.66323970  
 0.00000000 -0.24274639    0.46008182  0.46008182   -0.18827449  0.00000000   -0.46008182  0.46008182    0.00000000  0.24274639  
 0.11471345 -0.00000000    0.49337655 -0.49337655    0.00000000  0.00000000   -0.49337655 -0.49337655    0.11471345  0.00000000  
 0.00000000 -0.03445729    0.09136855  0.09136855    0.98195343  0.00000000   -0.09136855  0.09136855    0.00000000  0.03445729  
 0.69773981  0.00000000   -0.08111466  0.08111466   -0.00000000  0.00000000    0.08111466  0.08111466    0.69773981 -0.00000000  
When we run
x_dmft.py dmft1
again, we see in Ta2NiSe5O.outputdmf1 that the hybridizations on Ta are now mostly diagonal, i.e.,
 icix=           6
    1.30255268    -0.00000000      -0.02320125    -0.00000000       0.00000000     0.00000000       0.04495231     0.00000000       0.00000000     0.00000000   
   -0.02320125     0.00000000       1.45504129     0.00000000      -0.00000000     0.00000000      -0.08213064     0.00000000      -0.00000000    -0.00000000   
    0.00000000    -0.00000000      -0.00000000    -0.00000000       1.50592516    -0.00000000       0.00000000    -0.00000000       0.04745737     0.00000000   
    0.04495231    -0.00000000      -0.08213064    -0.00000000       0.00000000     0.00000000       1.88705739     0.00000000      -0.00000000     0.00000000   
    0.00000000    -0.00000000      -0.00000000     0.00000000       0.04745737    -0.00000000      -0.00000000     0.00000000       1.99335023    -0.00000000   
 
 icix=           6 at omega=0
    1.52609284    -0.00000000       0.00000000    -0.00000000       0.00000000     0.00000000       0.00000002     0.00000000      -0.00000000     0.00000000   
    0.00000000     0.00000000       2.32171379     0.00000000      -0.00000000     0.00000000      -0.00000003    -0.00000000      -0.00000000    -0.00000000   
    0.00000000    -0.00000000      -0.00000000    -0.00000000       5.14872492    -0.00000000       0.00000000    -0.00000000      -0.00000009     0.00000000   
    0.00000002    -0.00000000      -0.00000003     0.00000000       0.00000000     0.00000000      11.05305544     0.00000000      -0.00000000     0.00000000   
   -0.00000000    -0.00000000      -0.00000000     0.00000000      -0.00000009    -0.00000000      -0.00000000    -0.00000000      18.73265027    -0.00000000   
We mention in passing that the entire optimization of the orbitals can now be achieved with a single command
findRot.py
as discussed in LaVO3 tutorial. Newertheless, this step by step algorithm is also useful to understand. The second file that was created during initialization (dmft_init.py) is Ta2NiSe5O.indmfi. It contains two impurity problems, which will executed sequentially. No change is needed in this file.

Other eDMFT parameters

Both Ni and Ta harbor correlations. Estimates of U (by constrained DMFT) gives of the order of 6eV for Ta and around 8eV for Ni. We notice in passing that the eDMFT results do not depend strongly on U, and even U=9eV on both Ta and Ni gives quite similar results as U=6eV and U=8eV. We only need params.dat file, and blank self-energy on imaginary axis. We will use the following params.dat file
solver         =  'CTQMC'   # impurity solver
DCs            =  'exact'  # double counting scheme

max_dmft_iterations = 1     # number of iteration of the dmft-loop only
max_lda_iterations  = 100   # number of iteration of the LDA-loop only
finish         =  50        # number of iterations of full charge loop (1 = no charge self-consistency)

ntail          =  300       # on imaginary axis, number of points in the tail of the logarithmic mesh

cc             =  1e-5      # the charge density precision to stop the LDA+DMFT run
ec             =  1e-5      # the energy precision to stop the LDA+DMFT run

recomputeEF    = 1          # Recompute EF in dmft2 step. If recomputeEF = 2, it tries to find an insulating gap.

wbroad         =  0.0       # broadening of sigma on the imaginary axis
kbroad         =  0.0       # broadening of sigma on the imaginary axis

# Impurity problem number 0
iparams0={"exe"                : ["ctqmc"              , "# Name of the executable"],
          "U"                  : [8.0                  , "# Coulomb repulsion (F0)"],
          "J"                  : [0.8                  , "# Coulomb repulsion (F0)"],
          "CoulombF"           : ["'Ising'"            , "# Can be set to 'Full'"],
          "beta"               : [100                  , "# Inverse temperature"],
          "svd_lmax"           : [30                   , "# We will use SVD basis to expand G, with this cutoff"],
          "M"                  : [5e6                  , "# Total number of Monte Carlo steps"],
          "mode"               : ["SH"                 , "# We will use self-energy sampling, and Hubbard I tail"],
          "nom"                : [350                  , "# Number of Matsubara frequency points sampled"],
          "tsample"            : [400                  , "# How often to record measurements"],
          "GlobalFlip"         : [1000000              , "# How often to try a global flip"],
          "warmup"             : [3e5                  , "# Warmup number of QMC steps"],
          "nf0"                : [8.0                  , "# Nominal occupancy nd for double-counting"],
          }

# Impurity problem number 0
iparams1={"exe"                : ["ctqmc"              , "# Name of the executable"],
          "U"                  : [6.0                  , "# Coulomb repulsion (F0)"],
          "J"                  : [0.8                  , "# Coulomb repulsion (F0)"],
          "CoulombF"           : ["'Ising'"            , "# Can be set to 'Full'"],
          "beta"               : [100                  , "# Inverse temperature"],
          "svd_lmax"           : [30                   , "# We will use SVD basis to expand G, with this cutoff"],
          "M"                  : [5e6                  , "# Total number of Monte Carlo steps"],
          "mode"               : ["SH"                 , "# We will use self-energy sampling, and Hubbard I tail"],
          "nom"                : [350                  , "# Number of Matsubara frequency points sampled"],
          "tsample"            : [400                  , "# How often to record measurements"],
          "GlobalFlip"         : [1000000              , "# How often to try a global flip"],
          "warmup"             : [3e5                  , "# Warmup number of QMC steps"],
          "nf0"                : [2.0                  , "# Nominal occupancy nd for double-counting"],
          }
A few comments are in order. The impurity solve is CTQMC. We will use the exact double-counting, which is here essential, as we have two atoms, which have unusual valencies (far from nominal). We have two impurity problems, hence we need to specify iparams0 and iparams1 dictionaries. They differ only in the value of U, and nf0. The latter is important only for initialization of the self-energy. Namely, we need to start with a reasonable value of the double-counting, which determines the position of the impurity levels. If we start too far from the final occupancy, we will need very many iterations to converge. However, we have never found multiple solutions, i.e., possibility of getting different solution depending on what is choosen for nf0. Hence nf0 is here only to speed up the convergence.

Finally, we will slightly edit Ta2NiSe5O.indmfl file, to switch from real axis to imaginary axis, and than we will create a blank self-energy on imaginary axis. Open
Ta2NiSe5O.indmfl
and edit the second line, so that it starts with 1, i.e, matsubara=1
111 198 1 5                           # hybridization band index nemin and nemax, renormalize for interstitials, projection type
1 0.025 0.025 200 -3.000000 1.000000  # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)
Finally, remove the existing sig.inp file, and execute
szero.py
This should create new sig.inp file, which should start like this
  
# s_oo= [57.2, 57.2, 57.2, 57.2, 57.2, 8.6, 8.6, 8.6, 8.6, 8.6]
#  Edc= [57.2, 57.2, 57.2, 57.2, 57.2, 8.6, 8.6, 8.6, 8.6, 8.6]
   0.031415926535898  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
   0.094247779607694  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
   0.157079632679490  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
   0.219911485751286  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
....
Clearly, the self-energy at infinity and the double-counting correspond to our choosen nf0, U and J, using nominal formula U*(nf-1/2)-1/2*J*(nf-1). When computing the Green's function and hybridization, only s_oo-Edc is used, which is zero at the first iteration, hence we start with no self-energy. However, the impurity solver gets impurity level from Edc, and impurity level is rougly -Edc. Therefore if we choose correct nf0, the valence of the impurity will be close to converged valence from the very beginning.

The low temperature monoclinic phase

We repeat the above steps also for the low temperature monoclinic phase Ta2NiSe5O.struct. The dft-initialization by init_lapw as well as the dmft initialization by init_dmft.py are identical as above. Notice however that a monoclinic structure in wien2k needs the third angle gamma to be different from 90 degrees, which is properly done by the new version of the cif2struct converter. Consequently the space group convention used for this monoclinic phase has a, b, and c flipped compared to the the orthorombic structure. More precisely, the monoclinic aM=-a, bM=-c, and cM=-b. We want to choose the local axis on Mn and Ta atoms to be equal in both phases, hence we need to take into account this change of the real-space global axis.

The difference between the high and low-T calculation appears in optimizing the orbitals for this low-symmetry structure. The local axis is now slightly different, namely for Ni the script
localaxes.py
gives
 -0.71715269  -0.69691608   0.00000000  # local x-axis
 -0.69691608   0.71715269   0.00000000  # local y-axis
  0.00000000   0.00000000  -1.00000000  # local z-axis
and for Ta, we get
 -0.70889896  -0.28456479   0.64535660  # local x-axis
  0.70529973  -0.29095784   0.64644862  # local y-axis
  0.00381505   0.91343659   0.40696318  # local z-axis
The optimization of the orbitals is done in the same way as above, and we obtain the following set of orbitals for Ni
#---------------- # Independent components are --------------
'xz+yz' 'xy' 'xz-yz' 'z^2' 'x^2-y^2'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
 0.00000000  0.00000000    0.59348956  0.38440882    0.00000000  0.00000000   -0.59348956  0.38440882    0.00000000  0.00000000  
-0.00000000 -0.70710679    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.70710679
 0.00000000  0.00000000    0.38440882 -0.59348956    0.00000000  0.00000000   -0.38440882 -0.59348956    0.00000000  0.00000000  
 0.00000000  0.00000000    0.00000000  0.00000000    1.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000
 0.70710679  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.70710679  0.00000000
and the following set for Ta atom
#---------------- # Independent components are --------------
'xy' 'xz+yz' 'xz-yz' 'z^2' 'x^2-y^2'
#---------------- # Sigind follows --------------------------
 6  0  0  0  0
 0  7  0  0  0
 0  0  8  0  0
 0  0  0  9  0
 0  0  0  0 10
#---------------- # Transformation matrix follows -----------
-0.02919978 -0.48390653   -0.26387152  0.41063810    0.23121432  0.00000000    0.26387152  0.41063810   -0.02919978  0.48390653  
-0.01036665  0.13702680    0.51778210  0.44563043    0.16989067  0.00000000   -0.51778210  0.44563043   -0.01036665 -0.13702680  
-0.04911212 -0.49412774    0.40030280 -0.28446927   -0.15659305  0.00000000   -0.40030280 -0.28446927   -0.04911212  0.49412774  
-0.13888519  0.02205362    0.03219037 -0.22398487    0.92630351  0.00000000   -0.03219037 -0.22398487   -0.13888519 -0.02205362  
 0.69089711 -0.04908719    0.03154328 -0.04120560    0.18739669  0.00000000   -0.03154328 -0.04120560    0.69089711  0.04908719  
To compare, you might want to donwload Ta2NiSe5M.indmfl file. Note that params.dat file for the low-temperature and high temperature phase are identical.

Submitting and running the job

Now we want to start eDMFT calculatuion with the minimal number of necessary files. To do this, create a new folder (with any name), and when in that folder, use the command
dmft_copy.py <dft_results_directory>
where <dft_results_directory> is the directory where the output of the DFT calculation with DMFT initialization is. This command copies the necessary files to start a DMFT calculation. Now we are ready to submit the job. The executable script is
run_dmft.py
which calls all other executable. To make the code parallel, we need to prepare at least one file mpi_prefix.dat, but it is better two prepare also mpi_prefix.dat2. The first contains the parallel command used by ctqmc code, and the second is used by the dmft1 and dmft2 step. Please refer to tutorial on FeSe for these details.

The Ta2NiSe5 job is a bit more expensive, because there are two impurity problems to solve, and in particular Ta is very itinerant, hence showing very high perturbation order in ctqmc simulation. The histogram shows perturbation order over 1100 for Ta and 600 for Ni at 100K:
To speed up the calculation, one could increase temperature to reduce perturbation order. Newertheless, the convergence is quite rapid.

Since the symmetry is relatively low, and off-diagonal terms of hybridization tends to be generate during the run, epsecially on Ta, it is advisable to reorthogonalize the hybridization function a few times during the run. Just prepare rot.dat by taking data from Ta2NiSe5M.outputdmf1 and Ta2NiSe5M.indmfl, and execute find3dRotation.py rot.dat as before, and change Transformation matrix during the run. Make sure that you keep the same order of orbitals in the Transformation, otherwise the self-energy from previous iteration will not be compatible with the new choice of orbitals.
The converged hybridization at zero and large frequency (from Ta2NiSe5O.outputdmf1) on Ni and Ta in orthorombic phase should be similar to:
 Full matrix of impurity levels follows
 icix=           1
   -1.93346562     0.00000000       0.00000000    -0.00000000       0.00000000     0.00000000       0.00000000    -0.00000000       0.00000000    -0.00000000   
    0.00000000     0.00000000      -1.94286171    -0.00000000      -0.00000000     0.00000000       0.05339190    -0.00000000      -0.00000000    -0.00000000   
    0.00000000    -0.00000000      -0.00000000    -0.00000000      -2.06397297     0.00000000       0.00000000    -0.00000000       0.00000000    -0.00000000   
    0.00000000     0.00000000       0.05339190     0.00000000       0.00000000     0.00000000      -2.08102735    -0.00000000       0.00000000    -0.00000000   
    0.00000000     0.00000000      -0.00000000    -0.00000000       0.00000000    -0.00000000       0.00000000     0.00000000      -2.10262888     0.00000000   
 
 icix=           1 at omega=0
   -1.67969507     0.00000000       0.00000000    -0.00000000       0.00000000     0.00000000       0.00000000    -0.00000000       0.00000000     0.00000000   
    0.00000000     0.00000000      -2.23193952    -0.00000000      -0.00000000     0.00000000       0.00052338    -0.00000000       0.00000000     0.00000000   
    0.00000000    -0.00000000      -0.00000000    -0.00000000      -3.65710751     0.00000000       0.00000000    -0.00000000      -0.00000000    -0.00000000   
    0.00000000     0.00000000       0.00052338     0.00000000       0.00000000     0.00000000      -3.11368414    -0.00000000      -0.00000000    -0.00000000   
    0.00000000    -0.00000000       0.00000000    -0.00000000      -0.00000000     0.00000000      -0.00000000     0.00000000      -4.48062986     0.00000000   
.....
 icix=           3
   -0.32161541    -0.00000000       0.00375517     0.00000000      -0.00000000    -0.00000000       0.05617608     0.00000000       0.00000000     0.00000000   
    0.00375517    -0.00000000      -0.13557610    -0.00000000       0.00000000     0.00000000      -0.06914885    -0.00000000      -0.00000000    -0.00000000   
   -0.00000000    -0.00000000       0.00000000    -0.00000000      -0.09401932     0.00000000      -0.00000000    -0.00000000       0.02745149    -0.00000000   
    0.05617608    -0.00000000      -0.06914885     0.00000000      -0.00000000    -0.00000000       0.27253730     0.00000000      -0.00000000    -0.00000000   
    0.00000000     0.00000000      -0.00000000     0.00000000       0.02745149     0.00000000      -0.00000000     0.00000000       0.38293218    -0.00000000   
 
 icix=           3 at omega=0
   -0.77999416    -0.00000000      -0.00026256    -0.00000000       0.00000000    -0.00000000      -0.00093024     0.00000000      -0.00000000     0.00000000   
   -0.00026256     0.00000000       0.39222606    -0.00000000      -0.00000000     0.00000000       0.00098537    -0.00000000       0.00000000    -0.00000000   
    0.00000000     0.00000000      -0.00000000    -0.00000000       2.02177792     0.00000000      -0.00000000    -0.00000000       0.00261617     0.00000000   
   -0.00093024    -0.00000000       0.00098537     0.00000000      -0.00000000     0.00000000       6.15932483     0.00000000      -0.00000000     0.00000000   
   -0.00000000    -0.00000000       0.00000000     0.00000000       0.00261617    -0.00000000      -0.00000000    -0.00000000      11.07735301    -0.00000000   
while in monoclinic phase (from Ta2NiSe5M.outputdmf1) should be
 Full matrix of impurity levels follows
 icix=           1
   -1.83410635    -0.00000000      -0.00000000    -0.00000000      -0.02357286     0.00000000      -0.00000000    -0.00000000      -0.00000000     0.00000000   
   -0.00000000     0.00000000      -1.84971870    -0.00000000      -0.00000000     0.00000000       0.05626837     0.00000000      -0.00047353     0.00000000   
   -0.02357286    -0.00000000      -0.00000000    -0.00000000      -1.98027763     0.00000000      -0.00000000    -0.00000000       0.00000000    -0.00000000   
   -0.00000000     0.00000000       0.05626837    -0.00000000      -0.00000000     0.00000000      -1.99984875    -0.00000000       0.00910335    -0.00000000   
   -0.00000000    -0.00000000      -0.00047353    -0.00000000       0.00000000     0.00000000       0.00910335     0.00000000      -2.02272977    -0.00000000   
 
 icix=           1 at omega=0
   -1.70219389    -0.00000000       0.00000000     0.00000000       0.08654367    -0.00000000      -0.00000000     0.00000000      -0.00000000    -0.00000000   
   -0.00000000    -0.00000000      -2.29364202    -0.00000000      -0.00000000    -0.00000000      -0.00841268     0.00000000      -0.04648391     0.00000000   
    0.08654367     0.00000000      -0.00000000     0.00000000      -4.09259410     0.00000000      -0.00000000    -0.00000000      -0.00000000     0.00000000   
   -0.00000000    -0.00000000      -0.00841268    -0.00000000      -0.00000000     0.00000000      -3.11398196    -0.00000000      -0.00023317    -0.00000000   
   -0.00000000     0.00000000      -0.04648391    -0.00000000      -0.00000000    -0.00000000      -0.00023317     0.00000000      -4.64904718    -0.00000000   
....
 icix=           3
   -0.28410623     0.00000000      -0.00413198     0.00000000       0.03102557    -0.00000000       0.04128152    -0.00000000      -0.00071236    -0.00000000   
   -0.00413198    -0.00000000      -0.09877598    -0.00000000       0.01175201     0.00000000       0.06864815     0.00000000       0.00183100    -0.00000000   
    0.03102557     0.00000000       0.01175201    -0.00000000      -0.03689735    -0.00000000       0.06654623    -0.00000000      -0.04505486    -0.00000000   
    0.04128152     0.00000000       0.06864815    -0.00000000       0.06654623     0.00000000       0.31033130    -0.00000000       0.01039805    -0.00000000   
   -0.00071236     0.00000000       0.00183100     0.00000000      -0.04505486     0.00000000       0.01039805     0.00000000       0.41883044     0.00000000   
 
 icix=           3 at omega=0
   -0.60250154     0.00000000      -0.00277749    -0.00000000      -0.00140154     0.00000000       0.00179854     0.00000000      -0.00126328    -0.00000000   
   -0.00277749     0.00000000       0.28428764    -0.00000000      -0.00046112     0.00000000      -0.00238284     0.00000000      -0.00021410    -0.00000000   
   -0.00140154    -0.00000000      -0.00046112    -0.00000000       2.71506856    -0.00000000      -0.00626211    -0.00000000       0.00679348     0.00000000   
    0.00179854    -0.00000000      -0.00238284    -0.00000000      -0.00626211     0.00000000       6.01793755    -0.00000000      -0.00479451     0.00000000   
   -0.00126328     0.00000000      -0.00021410     0.00000000       0.00679348    -0.00000000      -0.00479451    -0.00000000      10.20280268     0.00000000   
The converged Ta2NiSe5O.indmfl file, which specifies the orbitals in the local coordinate axis of orthorombic run, is
....
2     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'xz+yz' 'xy' 'xz-yz' 'z^2' 'x^2-y^2'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
 0.00000000  0.00000000    0.50000000  0.50000000    0.00000000  0.00000000   -0.50000000  0.50000000    0.00000000  0.00000000
 0.00000000 -0.66293346    0.00000000  0.00000000    0.34790583  0.00000000    0.00000000  0.00000000    0.00000000  0.66293346
 0.00000000  0.00000000    0.50000000 -0.50000000    0.00000000  0.00000000   -0.50000000 -0.50000000    0.00000000  0.00000000
 0.00000000  0.24600657    0.00000000  0.00000000    0.93752948  0.00000000    0.00000000  0.00000000    0.00000000 -0.24600657
 0.70710679  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.00000000  0.00000000    0.70710679  0.00000000
3     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'xy' 'xz+yz' 'xz-yz' 'z2' 'x2-y2'
#---------------- # Sigind follows --------------------------
 6  0  0  0  0
 0  7  0  0  0
 0  0  8  0  0
 0  0  0  9  0
 0  0  0  0 10
#---------------- # Transformation matrix follows -----------
 0.00000000 -0.70194732   -0.05758947 -0.05758947    0.03569007  0.00000000    0.05758947 -0.05758947    0.00000000  0.70194732
 0.00000000 -0.08442539    0.49015525  0.49015525   -0.15727709  0.00000000   -0.49015525  0.49015525    0.00000000  0.08442539
 0.08506837  0.00000000    0.49636851 -0.49636851    0.00000000  0.00000000   -0.49636851 -0.49636851    0.08506837  0.00000000
 0.00000000  0.01193054    0.08019538  0.08019538    0.98690938  0.00000000   -0.08019538  0.08019538    0.00000000 -0.01193054
 0.70197107  0.00000000   -0.06015242  0.06015242    0.00000000  0.00000000    0.06015242  0.06015242    0.70197107  0.00000000
....
and for the monoclinic phase (Ta2NiSe5M.indmfl) is
....
2     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'xz+yz' 'xy' 'xz-yz' 'z^2' 'x^2-y^2' 
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
 0.00000000  0.00000000    0.54510879  0.45039586    0.00000000  0.00000000   -0.54510879  0.45039586    0.00000000  0.00000000
-0.00402422 -0.65017840    0.00000000  0.00000000    0.39306964  0.00000000    0.00000000  0.00000000   -0.00402422  0.65017840
 0.00000000  0.00000000    0.45039586 -0.54510879    0.00000000  0.00000000   -0.45039586 -0.54510879    0.00000000  0.00000000
 0.07091418  0.27617429    0.00000000  0.00000000    0.91509447  0.00000000    0.00000000  0.00000000    0.07091418 -0.27617429
 0.70353038 -0.03155676    0.00000000  0.00000000   -0.08999096  0.00000000    0.00000000  0.00000000    0.70353038  0.03155676
3     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'xy' 'xz+yz' 'xz-yz' 'z^2' 'x^2-y^2'
#---------------- # Sigind follows --------------------------
 6  0  0  0  0
 0  7  0  0  0
 0  0  8  0  0
 0  0  0  9  0
 0  0  0  0 10
#---------------- # Transformation matrix follows -----------
-0.00489259 -0.69416147    0.13411482 -0.00335834    0.01535367  0.00000000   -0.13411482 -0.00335834   -0.00489259  0.69416147
-0.01110550  0.08356195    0.43666523  0.53858106    0.15587173  0.00000000   -0.43666523  0.53858106   -0.01110550 -0.08356195
-0.10146142  0.10139414    0.49645295 -0.45391163    0.23204965  0.00000000   -0.49645295 -0.45391163   -0.10146142 -0.10139414
 0.07153571 -0.02612921   -0.18700502  0.01851820    0.95800429  0.00000000    0.18700502  0.01851820    0.07153571  0.02612921
 0.69601793  0.01391990    0.09950018 -0.05950198   -0.06204045  0.00000000   -0.09950018 -0.05950198    0.69601793 -0.01391990
...
We can plot lattice and impurity occupations as a number of iterations. This can be found in info.iterate, for Ni we compare columns 9,10 and for Ta columns 13,14. The plot should looks like
Notice that the impurity occupations starts around 8.1 and 2.08 for Ni and Ta, respectively, which is the result of our nf0 choice. With iterations, they converge to 8.46 and 2.25 on Ni and Ta in orthorombic phase, and to 8.48 and 2.23 in monoclinic phase, respectively. Across the phase transition, the occupation changes a bit, and Ni gains and Ta looses a small fraction (0.02) of an electron. Notice that this valence is determined by the exact double-counting, and different form of the double-counting (nominal or FLL) would change this occupancy. Notice also that in the exact double-counting the value of the double-counting potential (from sig.inp) is now orbitally dependent, and is
Edc=[59.50744806514058, 59.497945473363345, 59.46111542446767, 59.460057277212776, 59.490055307869966]
for Ni and
[8.421196770934927, 8.404531442236912, 8.361844515541126, 8.32570382132488, 8.288944700100567]
for Ta in monoclinic phase. Usually the orbital dependence of double-counting is very weak, but on Ta this is not the case, as the dominant xy orbital is pushed lower for 0.13eV as compared to the x2-y2 orbitals. This might help to open the gap in the low-temperature phase.

Plotting the orbitally resolved spectral functions

Next we want to plot results on the real axis. We can either take the last self-energy on imaginary axis, or average over a few last steps, if enough iterations is available (using tool saverage.py sig.inp.4[5-9].1). Than we create a new directory, and copy the resulting self-energy (sig.inpx) to the new directory. We also need maxent_params.dat, for which we use the same set of parameters as in other tutorials. We than execute
maxent_run.py sig.inpx
We next create yet another directory, in which the spectra will be plotted, and copy necessary files from the output directory to this new directory with
dmft_copy.py ../
and we take the real-axis self-energy produced by maxent_run.py and overwrite current sig.inp in this new directory
cp ../maxent/Sig.out sig.inp 
We also change matsubara flag to
0
in Ta2NiSe5O.indmfl file, so that the calculation will be performed on the real axis. We might want to parallelize the job from now on by specifying "mpi_prefix.dat", or, just increase number of cores by OMP_NUM_THREADS (but do not oversubscribe the machine by using both). We next calculate the density of states by executing
x lapw0  -f Ta2NiSe5O
x_dmft.py lapw1
x_dmft.py dmft1
we can check the content of Ta2NiSe5O.cdos, Ta2NiSe5O.gc1 and Ta2NiSe5O.cgc3 and plot density of state and orbitally resolved partial dos, which should look like
Notice that for smoother dos, we should increase number of k-points before running lapw1 and dmft1 above.

Next we create a k-point path along the high symmetry lines. In this step we could use wien2k tools to construct the path, or run xcrysden to plot the path. Alternatively, there is a klist3_generate.py, which can also generate a k-point path, but we need to slightly modify the script. Near the top there is a section specifying the path, which we will modify:
Nt = 200                                  # number of all points along the path
legend=['$Z$','$\Gamma$','$X/2$']         # name of  the points
BS = [[0, 0, 0.5], [0, 0, 0], [0.25,0,0]] # coordinates of the path. CAREFUL: Has to be specified in conventional brillouine zone, not primitive
BR1=[[0.94919,   0.00000,   0.00000],     # reciprocal vectors from case.rotlm of the conventional BZ.
     [0.00000,   0.25835,   0.00000],     # absolute value is not important, only ratios are important
     [0.00000,   0.00000,   0.21209]]     # Needed to find the length of each interval, and redistribute points along the path
Let's choose 200 points (variable Nt), and the path from Z=[0,0,1/2], Gamma, and X=[0.5373,0,0] direction (variables BS and legend), but lets only go to approximately half way, because the rest of the path contains no low energy bands. We also need to specify the reciprocal vectors of the conventional unit cell (BR1), which is written in Ta2NiSe5O.rotlm. We than execute
python klist3_generate.py > Ta2NiSe5O.klist_band
We will change the frequency range for the plot in Ta2NiSe5O.indmfl to -0.5 to 0.5eV, with 200 points. We edit the second line of this file to be:
0 0.025 0.025 200 -0.500000 0.500000  # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)
and than we execute
x_dmft.py lapw1 --band
x_dmft.py dmftp
Next we display the results by
wakplot.py 0.2
We can play with the intensity (0.2) to make plot brighter or dimmer, so that it best displays the bands. We should get a plot like
Occasionally, the plot might show unphysical jumps of spectra like that
which comes from the fact that projector introduces cuts (at -10eV and 10eV) where bands might be highly topologically entangled. In this case, usually increasing the number of bands for just one, will eliminate the jumps (first line in Ta2NiSe5O.indmfl).

To orbitally resolve the spectra, we use dmftgk tool, already described in FeSe tutorial. At this stage, it is important to print projector into a single file "Udmft.0", which can be achieved when mpi parallelization is switched off. We hence eliminate "mpi_prefix.dat" and print the projector by executing
x_dmft.py dmftu -g --band
Next we prepare a file dmftgke.in, which should look like:
e                    # mode [g/e]: we use mode to compute eigenvalues and eigenvectors
BasicArrays.dat      # filename for projector
0                    # matsubara
Ta2NiSe5O.energy     # LDA-energy-file, case.energy(so)(updn)
Ta2NiSe5O.klist_band # k-list
Ta2NiSe5O.rotlm      # for reciprocal vectors
Udmft.0              # filename for projector
0.0025               # gamma for non-correlated
0.0025               # gammac
sig.inp1_band  sig.inp2_band sig.inp3_band sig.inp4_band sig.inp5_band sig.inp6_band # self-energy filenames
eigenvalues.dat     # eigenvalues
UR.dat_             # right eigenvector in orbital basis
UL.dat_             # left eigenvector in orbital basis
-1.                 # emin for printed eigenvalues
 1.                 # emax for printed eigenvalues
Next we execute
dmftgk dmftgke.in
to print the projectors as well as the eigenvalues to UR.dat_, UL.dat_ and eigenvalues.dat. We then copy wakplot_sophisticated.py script, and slightly edit it.
if __name__ == '__main__':
    aspect_scale=1.0
    ....
    ....
    COHERENCE=True
    orb_plot=[0,1,3,3,3,3,3,3,3,3,2,3,3,3,3]  # should have integers 0,1,2 for red,green,blue
Here we changed "aspect_scale=1.0" from 0.3 to 1.0, and we edited "orb_plot". In the projectors "UL.dat_", we have 6x5 orbitals printed, in which first ten correspond to Ni, and the other 20 to Ta. For the red/green/blue color, we will choose red==Ni:xz+yz, green==Ni:xy, blue:Ta:xy. Here Ni:xz+yz and Ni:xy appear as the 0-th and 1-th component, while Ta:xy is the 10-th component, because there are two Ni-atoms, each with 5 orbitals, printed before. The components that we do not want to hghlight (they will appear black) get number 3. We can also skip the rest of the orbitals, which we are not interested in. We than display the result by
python wakplot_sophisticated.py eigenvalues.dat 
which should look like:


Next we repeat the same procedure for the low-temperature monoclinic phase, which reveals a small semiconducting gap.
We can convince ourselves that it is semiconducting gap by looking into the self-energy file sig.inp, which shows very small scattering rate for all orbitals. The orbitally resolved spectra of the monoclinic phase than looks like
which has an unusual type of shape, and lead some to argue that bands are formed by excitonic condensation.

Finding the irreducible representation of the quasiparticle bands

The spectra which is very sharp has well defined quasiparticles, which can be assigned irreducible representation, just like in the band-theory case. We will use wien2k irrep tool to do that. This tool needs a vector file, which contains quasiparticle bands. To obtain such bands, we will run dmftp tool with the following command
x_dmft.py dmftp --hermitianw
which forces the (Hk+Sigma(0)) to be symmetrized (Hermitian) and writes a new vector file Ta2NiSe5M.vector_dmft to the disc. Notice that here we take self-energy at zero-frequency, hence the bands will look somewhat different. Notice that this part has to be execute in serial, and one should remove "mpi_prefix.dat" from the current directory.

To see the quasiparticle-bands, which will be analized by irrep, we need to display energies from Ta2NiSe5M.outputdmfp. We could use wien2k tool spaghetti for that purpose, but it is a bit tedious to substitute Ta2NiSe5M.output1 by Ta2NiSe5M. outputdmfp, and also resulting data is hard to analize in poscript form. We therefore use alternative tool with similar capability x_spaght.py. We will execute:
x_spaght.py Ta2NiSe5M.outputdmfp -y-0.5:0.5 -x: -g -Ry
to display bands between -0.5 and 0.5eV:
We choose to plot the bands in Ry units (-Ry) rather than eV, because it will be easier to compare with data in case.outputir. We also choose empty x-axis (-x:), which adds integer labels to each k-point, which will be convenient in analyzing case.outputir. We also show the grid (-g). This plot is convenient to find which k-point is being analized by irrep. Note that the mouse will show energy and momentum index of each point.
Next we produce irrep.def file for our analisys of dmft-vector file:
x irrep -f Ta2NiSe5M -d
This produces irrep.def. We now need to replace Ta2NiSe5M.vector with Ta2NiSe5M.vector_dmft in irrep.def file, i.e, irrep def should have a line:
10,'./Ta2NiSe5M.vector_dmft',      'old',    'unformatted',9000
Finally, we execute
$WIENROOT/SRC_irrep/irrep irrep.def
and we should get Ta2NiSe5M.irrep and Ta2NiSe5M.outputir. We can then browse through Ta2NiSe5M.outputir, and with cursor on matplolib plots analyze irreps. For the monoclinic phase we get:
We than repeat the same procedure for orhorombic phase, and the quasiparticle plot is
while the annotated spectral function is:

Plotting the Fermi surface of the orthorombic phase

To plot the Fermi surface, we will create another directory named FS, and while in that directory we will copy crucial files from previous real-axis work, i.e.,
dmft_copy.py ../onreal/
cp  ../onreal/SOlapm.dat .
Since we will not run dmft1 again, we also need the renormalization factors for the projector, stored in SOlapm.dat. We then construct a very large k-point mesh for 3D plotting of the Fermi surface.
x kgen -f Ta2NiSe5O
and choose 5000 k-points, and we do not shift the mesh. We then compute eigenvalues on this mesh by
x lapw0 -f Ta2NiSe5O
x_dmft.py lapw1
note that here we can use parallel version of the codes, hence "mpi_prefix.dat" can be prepared, and written to the disc in the current directory.

Next we construct self-energy with a single frequency in sig.inp file, namely the zero frequency. We eliminate all other frequencies.

We than run dmft0 tool, which produces energies for all k-points. We can either execute
x_dmft.py dmft0
or in metallic phase we can also assume that at the Fermi level the scattering rate vanishes, therefore the eigenproblem is hermitian
x_dmft.py dmft0 --hermitian
Finally, we can use FSconvert.py tool, which produces input file for xcrysden, i.e.,
python FSconvert.py 
We should get Ta2NiSe5O.bxsf file, which can be opened by
xcrysden --bxsf Ta2NiSe5O.bxsf
The resulting Fermi surface, looking from the top down, is

Constructing the cluster: Ti-Ni-Ti trimer

In the monoclinic phase, the chains of Ta atoms move parallel to the chain direction, as shown below with orange arrows. As a result, Ni-Ta distances, which were all equal in orthorombic phase, disproportionate and hence each Ni atom comes closer with two Ta octahedra (red arrows), while the other two octahedra move away (navy blue arrows). As a result, we have a trimer formation, indicated by red arrows in the figure.
There are two non-overlaping trimes in the unit cell, with positions specified in the above figure. We notice that Ni1 from structure file forms bond with Ta1 and Ta3, while Ni2 forms bond with Ta2 and Ta4. The positions of Ni1, Ta1, and Ta3 from the structure file are (0,3/4,z0), (1-x,y,z), and (x,3/2-y,z). However, it turns out that in the calculation a different set of atoms is actually used. Even more confusing is the fact that positions of the atoms used in the calculation depends on the order of the symmetry operations specified in the structure file. If we reorder the symmetry operations, we will be using a different (but equivalent) set of atoms in the actual calculations. Within local approximations, this is not an issue, however, for cluster constructions the bond between Ni1:Ta(1-x,y,z) is not the same as the bond between Ni1:Ta(-x,y,z). Hence, we need to carefully select correct atoms which form clusters.

The file Ta2NiSe5M.outputdmf1 prints the actual positions of atoms, as used in the calculations. The positions are:
Actual position of atom  1 is:  0.000000    0.750000    0.201210
Actual position of atom  2 is:  0.000000   -0.750000   -0.201210
Actual position of atom  3 is:  0.512048    0.610655    0.221366
Actual position of atom  4 is: -0.512048   -0.610655   -0.221366
Actual position of atom  5 is: -0.512048   -0.110655    0.221366
Actual position of atom  6 is:  0.512048    1.110655   -0.221366
which are not exactly the same as those specified in the structure file. We notice that Ni1 and Ni2 are already positions exactly as indicated in the above figure, namely at Ni1=(0,3/4,z0) and Ni2=(0,-3/4,-z0). However, Ta1=atom(3) is at (-x+1,y,z) and Ta3=atom(5) is at (x-1,3/2-y-1,z). We thus need to add a shift of (-1,0,0) to Ta1 and (1,1,0) to Ta3. Furthermove, Ta2=atom(4) is at (x-1,-y,-z) and Ta4=atom(6) is at (-x+1,y-3/2+2,-z). We thus require shift of (1,0,0) and (-1,-2,0) for Ta2 and Ta4, respectively.

Next we will modify Ta2NiSe5M.indmfl file, so that it will correspond to the cluster (cellular-DMFT type) calculation. We create a new directory, say nlHF, and copy previous real-axis data into this directory by
dmft_copy.py ../onreal/
To indicate that a set of atoms form a cluster in real space, we must give all such atom the same index cix in the Ta2NiSe5M.indmfl file. Each atom contains a line like
  2   2   1                           # L, qsplit, cix
where the last number specifies cix. Notice that currently every atom has its own cix-index, i.e.,
  2   2   1                           # L, qsplit, cix for Ni1
  2   2   2                           # L, qsplit, cix for Ni2
  2   2   3                           # L, qsplit, cix for Ta1
  2   2   4                           # L, qsplit, cix for Ta2
  2   2   5                           # L, qsplit, cix for Ta3
  2   2   6                           # L, qsplit, cix for Ta4
This indicates that each atom is its own cluster, i.e., corresponding to single-site DMFT approximation. Now we want to modify these lines, to indicate the correct grouping, i.e.,
  2   2   1                           # L, qsplit, cix for Ni1
  2   2   2                           # L, qsplit, cix for Ni2
  2   2   1                           # L, qsplit, cix for Ta1
  2   2   2                           # L, qsplit, cix for Ta2
  2   2   1                           # L, qsplit, cix for Ta3
  2   2   2                           # L, qsplit, cix for Ta4
This means that Ta1 and Ta3 form a cluster with Ni1, and Ta2 and Ta4 form a cluster with Ni2.

Next we want to take care of the necessary shifts of atoms. Currently, we are specifying local rotation for each atom, but now we will require also shift, in addition to rotation. For that, we will change locrot, which is currently -1 (indicating nontrivial local rotation) to -4. [ Notice that if local rotation was not necessary, we could use flag -3, and omit local rotation. ] When locrot is -4 it means that we will specify below 3x3 rotation as well as shift of an atom. For example, for Ta(1), which currently contains this block
3     1  -1                           # iatom, nL, locrot
  2   2   3                           # L, qsplit, cix
-0.70889896 -0.28456479  0.64535660   # new x-axis
 0.70529973 -0.29095784  0.64644862   # new y-axis
 0.00381505  0.91343659  0.40696318   # new z-axis
we will modify to
3     1  -4                           # iatom, nL, locrot
  2   2   1                           # L, qsplit, cix
-0.70889896 -0.28456479  0.64535660   # new x-axis
 0.70529973 -0.29095784  0.64644862   # new y-axis
 0.00381505  0.91343659  0.40696318   # new z-axis
-1 0 0
We have done three things: i) changed cix from 3 to 1, changed locrot from -1 to -4, and we specified in the last line the necessary shift of the Ta(1) atom.
We now repeat this procedure for all atoms. The resulting hader of the Ta2NiSe5M.indmfl file will look like
111 198 1 5                           # hybridization band index nemin and nemax, renormalize for interstitials, $
0 0.025 0.025 200 -0.500000 0.500000  # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_$
6                                     # number of correlated atoms
1     1  -1                           # iatom, nL, locrot
  2   2   1                           # L, qsplit, cix
-0.71715269 -0.69691608  0.00000000   # new x-axis
-0.69691608  0.71715269  0.00000000   # new y-axis
 0.00000000  0.00000000 -1.00000000   # new z-axis
2     1  -1                           # iatom, nL, locrot
  2   2   2                           # L, qsplit, cix
-0.71715269 -0.69691608  0.00000000   # new x-axis
-0.69691608  0.71715269  0.00000000   # new y-axis
 0.00000000  0.00000000 -1.00000000   # new z-axis
3     1  -4                           # iatom, nL, locrot
  2   2   1                           # L, qsplit, cix
-0.70889896 -0.28456479  0.64535660   # new x-axis
 0.70529973 -0.29095784  0.64644862   # new y-axis
 0.00381505  0.91343659  0.40696318   # new z-axis
-1 0 0
4     1  -4                           # iatom, nL, locrot
  2   2   2                           # L, qsplit, cix
-0.70889896 -0.28456479  0.64535660   # new x-axis
 0.70529973 -0.29095784  0.64644862   # new y-axis
 0.00381505  0.91343659  0.40696318   # new z-axis
1 0 0
5     1  -4                           # iatom, nL, locrot
  2   2   1                           # L, qsplit, cix
-0.70889896 -0.28456479  0.64535660   # new x-axis
 0.70529973 -0.29095784  0.64644862   # new y-axis
 0.00381505  0.91343659  0.40696318   # new z-axis
1 1 0
6     1  -4                           # iatom, nL, locrot
  2   2   2                           # L, qsplit, cix
-0.70889896 -0.28456479  0.64535660   # new x-axis
 0.70529973 -0.29095784  0.64644862   # new y-axis
 0.00381505  0.91343659  0.40696318   # new z-axis
-1 -2 0
Now that we specified the two clusters, we also need to modified the definition of the orbitals. The cluster will require 15 orbitals rather than 5, and each orbital will have 15 dimensional vector, i.e, the transformation must be 15x15 rather than 5x5. But there will be only two such entries, corresponding to the two clusters, rather than previously 6 entries.

The next lines specifies the number of cix-blocks and maximum dimension, which used to be
#================ # Siginds and crystal-field transformations for correlated orbitals ================
6     5   5       # Number of independent kcix blocks, max dimension, max num-independent-components
1     5   5       # cix-num, dimension, num-independent-components
and should be changed to
#================ # Siginds and crystal-field transformations for correlated orbitals ================
2    15  12       # Number of independent kcix blocks, max dimension, max num-independent-components
1    15  12       # cix-num, dimension, num-independent-components
In addition to 10 different diagonal components of the self-energy, which we have in single site DMFT, we can here add any number of off-diagonal components. For simplicity, we will concentrate only on the low energy orbitals, here Ni(xz+yz) and Ta(xy). Since there are two Ta atoms, we need to add at least two such components, hence we will have 10+2=12 components in total.

Next we need to specify the order of the orbitals (their names), the Sigind block, and the orbital Transformation. Lets just add Ni:d orbitals first, followed by Ta1, and Ta3 orbitals. Sigind is than just a combination of Ni and Ta1+Ta3 components.
#---------------- # Independent components are --------------
'Ni(xz+yz)' 'Ni(xy)' 'Ni(xz-yz)' 'Ni(z^2)' 'Ni(x^2-y^2)' 'Ta1(xy)' 'Ta1(xz+yz)' 'Ta1(xz-yz)' 'Ta1(z^2)' 'Ta1(x^2-y^2)' 'Ta2(xy)' 'Ta2(xz+yz)' 'Ta2(xz-yz)' 'Ta2(z^2)' 'Ta2(x^2-y^2)'
#---------------- # Sigind follows --------------------------
 1  0  0  0  0 -11 0  0  0  0 -12 0  0  0  0
 0  2  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  3  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  4  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  5  0  0  0  0  0  0  0  0  0  0
-11 0  0  0  0  6  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  7  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  8  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  9  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0 10  0  0  0  0  0
-12 0  0  0  0  0  0  0  0  0  6  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  7  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  8  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  9  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0 10
Here we also added the off-diagonal components between the low energy orbitals, i.e., between 1-st and 6-th and 1-st and 11-th orbital. Notice that 6-th orbital corresponds to Ta1(xy) and 11-th orbital is Ta3(xy) orbital, while 1-st is Ni(xz+yz).
We could use positive numbers (11 and 12) for these off-diagonal components. However, we than must add four more columns into the self-energy (sig.inp) file, i.e., two new components of the self-energy with real and imaginary part. If, however, a component has a negative index, it means that it is frequency independent, and just needs to be specified at the top of sig.inp file by modifying s_oo, which is the self-energy at infinity. In this case, the dynamic part of the self-energy does not need any extra columns. This is convenient when additing Hartree-Fock off-diagonal self-energy, which is static.

Finally, the Transformation needs to reflect the chosen order of the orbitals, i.e.,
 0.00000000  0.00000000    0.54510879  0.45039586    0.00000000  0.00000000   -0.54510879  0.45039586    0.00000000  0.00000000   0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0
-0.00402422 -0.65017840    0.00000000  0.00000000    0.39306964  0.00000000    0.00000000  0.00000000   -0.00402422  0.65017840   0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0
 0.00000000  0.00000000    0.45039586 -0.54510879    0.00000000  0.00000000   -0.45039586 -0.54510879    0.00000000  0.00000000   0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0
 0.07091418  0.27617429    0.00000000  0.00000000    0.91509447  0.00000000    0.00000000  0.00000000    0.07091418 -0.27617429   0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0
 0.70353038 -0.03155676    0.00000000  0.00000000   -0.08999096  0.00000000    0.00000000  0.00000000    0.70353038  0.03155676   0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 -0.00489259 -0.69416147    0.13411482 -0.00335834    0.01535367  0.00000000   -0.13411482 -0.00335834   -0.00489259  0.69416147   0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 -0.01110550  0.08356195    0.43666523  0.53858106    0.15587173  0.00000000   -0.43666523  0.53858106   -0.01110550 -0.08356195   0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 -0.10146142  0.10139414    0.49645295 -0.45391163    0.23204965  0.00000000   -0.49645295 -0.45391163   -0.10146142 -0.10139414   0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0  0.07153571 -0.02612921   -0.18700502  0.01851820    0.95800429  0.00000000    0.18700502  0.01851820    0.07153571  0.02612921   0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0  0.69601793  0.01391990    0.09950018 -0.05950198   -0.06204045  0.00000000   -0.09950018 -0.05950198    0.69601793 -0.01391990   0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0   -0.00489259 -0.69416147    0.13411482 -0.00335834    0.01535367  0.00000000   -0.13411482 -0.00335834   -0.00489259  0.69416147
 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0   -0.01110550  0.08356195    0.43666523  0.53858106    0.15587173  0.00000000   -0.43666523  0.53858106   -0.01110550 -0.08356195
 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0   -0.10146142  0.10139414    0.49645295 -0.45391163    0.23204965  0.00000000   -0.49645295 -0.45391163   -0.10146142 -0.10139414
 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0    0.07153571 -0.02612921   -0.18700502  0.01851820    0.95800429  0.00000000    0.18700502  0.01851820    0.07153571  0.02612921
 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0    0.69601793  0.01391990    0.09950018 -0.05950198   -0.06204045  0.00000000   -0.09950018 -0.05950198    0.69601793 -0.01391990
Finally, this structure of Sigind and transformation has to be repeated twice, as we have two clusters. Notice that the two clusters here are equivalent, and hence these entries are identical, except the first line needs to reflect cix-num=2:
2    15   12       # cix-num, dimension, num-independent-components
Finally, we will modify the self-energy sig.inp file. Currently, the top has the following entries
# s_oo= [58.34111987610038, 57.99864526007574, 58.66792046452633, 58.15374868968896, 58.68284274853787, 11.52900772267362, 11.456232499198762, 11.35548455310819, 11.192604230189128, 11.111987546730667]
# Edc= [59.842921188735275, 59.83356456338613, 59.79961530249243, 59.79826650028635, 59.826546613721476, 8.319863999593021, 8.299004410493579, 8.302522649084414, 8.236828534352366, 8.22545468607832]
We will add two extra zeros at the end for our 11-th and 12-th self-energy component. We modify both s_oo and Edc.
  
# s_oo= [58.34111987610038, 57.99864526007574, 58.66792046452633, 58.15374868968896, 58.68284274853787, 11.52900772267362, 11.456232499198762, 11.35548455310819, 11.192604230189128, 11.111987546730667, 0, 0]
# Edc= [59.842921188735275, 59.83356456338613, 59.79961530249243, 59.79826650028635, 59.826546613721476, 8.319863999593021, 8.299004410493579, 8.302522649084414, 8.236828534352366, 8.22545468607832, 0, 0]
Notice that the double-counting for off-diagonal components should always vanish, as the LDA/GGA does not contain any correlations between sites in the LDA/GGA Luttinger-Ward functional. The self-energy s_oo will eventually be modified within Hartree-Fock method.

To check that Ta2NiSe5M.indmfl was properly modified, we will plot density of states, and compare with previous results. Since the off-diagonal self-energy is set to zero here, the results for density of states and all green's functions should be identical. Notice however, that the hybridizations can be very different, because the hybridization contains the inverse of the green's function, which is now a much larger (15x15) matrix. The hybridization now correponds to the cluster-DMFT hybridization for the trimer.
x lapw0 -f Ta2NiSe5M
x_dmft.py lapw1
x_dmft.py dmft1
Note that these steps can be done in parallel. Here is the complete Ta2NiSe5M.indmfl file, which we just created.
The local green's function from previous real-axis run and this one should match, which means that columns 1 to 11 of previous ../onreal/Ta2NiSe5M.gc1 and current Ta2NiSe5M.gc1 should match. In addition, the columns of 12 to 21 of Ta2NiSe5M.gc1 should match columns 2 to 11 of previous ../onreal/Ta2NiSe5M.gc3. For example, a plot of
plot -g -u1:3,1:13 -x-3:3 Ta2NiSe5M.gc1 -u1:3 ../onreal/Ta2NiSe5M.gc1 ../onreal/Ta2NiSe5M.gc1 
should result in plot like that:

We also need to check the symmetry of the off-diagonal green's function by displaying 23-th and 25-th column (11-th and 12-th component of the off-diagonal green's function).
We notice that in this monoclinic phase, the symmetry requires Ni-Ta1 Green's function to have opposite sign to Ni-Ta2 bond. We could add an extra pi phase to the Ta2 orbitals, and than the two components would be identical. We note that the self-energy and the Green's function have the same symmetry, hence if we redefine the Ta2 orbitals with opposite sign, we would need only 11 component self-energy rather than 12 component. But we will not do that now.

Non-local Hartree-Fock

Now that we have the trimer cluster, we can perform Hartree-Fock calculation for the non-local components. In particular, we can consider the effects on non-local components of the screened Coulomb repulsion, i.e., .
The Hartree part of the Coulomb interaction has been already correctly taken into account by LDA/GGA, hence we should not double-count the Hartree-part. However, the exchange and correlation between sites is completely missing in LDA/GGA or DMFT, hence at the lowest order, we need to consider the exchange contribution of the form
To consider the self-energy between Ni(xz+yz) and Ta(xy) orbital, we thus need the corresponding entry in the density matrix. We will thus run dmft1 on imaginary axis, and check the density matrix in Ta2NiSe5M.outdmft1.

We first save the real axis self-energy and replace it with the imaginary axis self-energy.
mv sig.inp sig.inp_real_axis
cp ../sig.inp .
We than add two more components to the header of the self-energy, just like before, i.e., the existing header
# s_oo= [58.33960607278515, 57.99799075799591, 58.66764718991543, 58.1528394846167, 58.68155194019003, 11.52723220068357, 11.45504198755027, 11.35418395795326, 11.19128337408098, 11.11088607272979]
# Edc= [59.841935471914276, 59.83260781272945, 59.79850099666909, 59.79715237803764, 59.82554966046209, 8.318434113332705, 8.297724288472981, 8.301180281321397, 8.23549025137547, 8.224172638351167]
is replaced by
# s_oo= [58.33960607278515, 57.99799075799591, 58.66764718991543, 58.1528394846167, 58.68155194019003, 11.52723220068357, 11.45504198755027, 11.35418395795326, 11.19128337408098, 11.11088607272979, 0, 0]
# Edc= [59.841935471914276, 59.83260781272945, 59.79850099666909, 59.79715237803764, 59.82554966046209, 8.318434113332705, 8.297724288472981, 8.301180281321397, 8.23549025137547, 8.224172638351167, 0, 0]
We change the flag matsubara to 1 in Ta2NiSe5M.indmfl to perform calculation on imaginary axis. And than we run
x_dmft.py dmft1
Than we look into Ta2NiSe5M.outdmft1 how the density matrix looks like
:NCOR   12.994665  1   15  # nf, icix, cixdm
   1.7166866   -0.0000000     -0.0000755   -0.0000000     -0.0040139   -0.0000000     -0.0000393   -0.0000000     -0.0000129    0.0000000     -0.0275268   -0.0000000     -0.0872614    0.0000000     -0.0468143   -0.0000000      0.0544655   -0.0000000     -0.0224366  
  -0.0000755    0.0000000      1.8119760    0.0000000      0.0000172   -0.0000000     -0.0039443   -0.0000000     -0.0048936   -0.0000000      0.0955963    0.0000000      0.0540922    0.0000000     -0.0030930    0.0000000      0.0239057   -0.0000000     -0.0313942  
  -0.0040139    0.0000000      0.0000172    0.0000000      1.6023202    0.0000000      0.0000040   -0.0000000     -0.0000009    0.0000000     -0.0500300    0.0000000      0.0658608   -0.0000000      0.1551905   -0.0000000     -0.0459625   -0.0000000     -0.0036675  
  -0.0000393    0.0000000     -0.0039443    0.0000000      0.0000040    0.0000000      1.7505693    0.0000000     -0.0135653    0.0000000      0.0110898   -0.0000000      0.0319096    0.0000000      0.1413475   -0.0000000     -0.0643905   -0.0000000      0.0137550  
  -0.0000129   -0.0000000     -0.0048936    0.0000000     -0.0000009   -0.0000000     -0.0135653   -0.0000000      1.6309636   -0.0000000     -0.0067397   -0.0000000     -0.1657144    0.0000000     -0.1909387   -0.0000000      0.1093709    0.0000000      0.0035428  
  -0.0275268    0.0000000      0.0955963   -0.0000000     -0.0500300   -0.0000000      0.0110898    0.0000000     -0.0067397    0.0000000      0.3507117   -0.0000000     -0.0196900    0.0000000      0.0080640   -0.0000000      0.0266478    0.0000000      0.0007574  
  -0.0872614   -0.0000000      0.0540922   -0.0000000      0.0658608    0.0000000      0.0319096   -0.0000000     -0.1657144   -0.0000000     -0.0196900   -0.0000000      0.3827676    0.0000000      0.0096011    0.0000000     -0.0031817   -0.0000000      0.0021048  
  -0.0468143    0.0000000     -0.0030930   -0.0000000      0.1551905    0.0000000      0.1413475    0.0000000     -0.1909387    0.0000000      0.0080640    0.0000000      0.0096011   -0.0000000      0.4679031   -0.0000000     -0.0091721    0.0000000      0.0022629  
   0.0544655    0.0000000      0.0239057    0.0000000     -0.0459625    0.0000000     -0.0643905    0.0000000      0.1093709   -0.0000000      0.0266478   -0.0000000     -0.0031817    0.0000000     -0.0091721   -0.0000000      0.4997402    0.0000000     -0.0008033  
  -0.0224366    0.0000000     -0.0313942    0.0000000     -0.0036675    0.0000000      0.0137550    0.0000000      0.0035428   -0.0000000      0.0007574    0.0000000      0.0021048   -0.0000000      0.0022629   -0.0000000     -0.0008033   -0.0000000      0.5399797  
   0.0275940   -0.0000000      0.0955899   -0.0000000      0.0500146    0.0000000      0.0110912   -0.0000000     -0.0067380    0.0000000      0.0042034   -0.0000000      0.0000261    0.0000000      0.0104573    0.0000000      0.0070014   -0.0000000     -0.0085243  
   0.0872932   -0.0000000      0.0540954   -0.0000000     -0.0658473   -0.0000000      0.0319100   -0.0000000     -0.1657210   -0.0000000     -0.0000452   -0.0000000     -0.0051224    0.0000000     -0.0226664   -0.0000000      0.0233843    0.0000000     -0.0035038  
   0.0468201    0.0000000     -0.0031192   -0.0000000     -0.1551889   -0.0000000      0.1413369    0.0000000     -0.1909445   -0.0000000      0.0104435   -0.0000000     -0.0226546   -0.0000000     -0.0497412   -0.0000000      0.0170714    0.0000000      0.0019080  
  -0.0544860   -0.0000000      0.0239183    0.0000000      0.0459665   -0.0000000     -0.0643833   -0.0000000      0.1093756    0.0000000      0.0069854    0.0000000      0.0233694   -0.0000000      0.0170699    0.0000000     -0.0061118   -0.0000000     -0.0023756  
   0.0224315    0.0000000     -0.0314047   -0.0000000      0.0036647    0.0000000      0.0137501   -0.0000000      0.0035420    0.0000000     -0.0085373   -0.0000000     -0.0034940   -0.0000000      0.0019112    0.0000000     -0.0023789    0.0000000      0.0010117  
We notice that all off-diagonal components are pretty small, and hence in the high-temperature phase the corrections to the band structure due to screened non-local contribution is small. However, in the symmetry broken phase, some components which are not allowed in the orthorombic phase, can become nonzero. The non-local Ni(xz+yz):Ta(xy), which is finite only within our chosen trimer, (non-zero on red bonds in the figure, and vanishes on the blue) breaks the symmetry of the high-temperature phase, and is thus not allowed in the orthorombic phase. But is allowed in the monoclinic phase. We will thus set it to some nonzero-value, and optimize it self-consistently for choosen value of non-local Coulomb interaction V. Let's consider V=3eV, which is half of the on-site interaction on Ta. We will set self-energy first to a value 0.35eV, i.e., header of the sig.inp becomes
# s_oo= [58.33960607278515, 57.99799075799591, 58.66764718991543, 58.1528394846167, 58.68155194019003, 11.52723220068357, 11.45504198755027, 11.35418395795326, 11.19128337408098, 11.11088607272979, -0.35,0.35]
# Edc= [59.841935471914276, 59.83260781272945, 59.79850099666909, 59.79715237803764, 59.82554966046209, 8.318434113332705, 8.297724288472981, 8.301180281321397, 8.23549025137547, 8.224172638351167, 0.0, 0.0]
and we will rerun
x_dmft.py dmft1
:NCOR   12.985283  1   15  # nf, icix, cixdm
   1.7063583   -0.0000000     -0.0000687    0.0000000      0.0235177   -0.0000000     -0.0000345    0.0000000     -0.0000082    0.0000000      0.1168346   -0.0000000     -0.0826159    0.0000000     -0.0411286   -0.0000000      0.0539939   -0.0000000     -0.0210488  
  -0.0000687   -0.0000000      1.8050129    0.0000000      0.0000301   -0.0000000     -0.0021135   -0.0000000     -0.0051496    0.0000000      0.0973357    0.0000000      0.0557467    0.0000000     -0.0075498   -0.0000000      0.0233133    0.0000000     -0.0325156  
   0.0235177    0.0000000      0.0000301    0.0000000      1.6058929    0.0000000      0.0000002   -0.0000000      0.0000003    0.0000000     -0.0455411    0.0000000      0.0676536    0.0000000      0.1479968   -0.0000000     -0.0424905   -0.0000000     -0.0063524  
  -0.0000345   -0.0000000     -0.0021135    0.0000000      0.0000002    0.0000000      1.7506449    0.0000000     -0.0135649   -0.0000000      0.0113038    0.0000000      0.0312383    0.0000000      0.1421482   -0.0000000     -0.0638839    0.0000000      0.0138251  
  -0.0000082   -0.0000000     -0.0051496   -0.0000000      0.0000003   -0.0000000     -0.0135649    0.0000000      1.6311759    0.0000000     -0.0069271   -0.0000000     -0.1662423    0.0000000     -0.1908940   -0.0000000      0.1094664    0.0000000      0.0033805  
   0.1168346    0.0000000      0.0973357   -0.0000000     -0.0455411   -0.0000000      0.0113038   -0.0000000     -0.0069271    0.0000000      0.3609223    0.0000000     -0.0230903    0.0000000     -0.0085177   -0.0000000      0.0382409    0.0000000     -0.0075248  
  -0.0826159   -0.0000000      0.0557467   -0.0000000      0.0676536   -0.0000000      0.0312383   -0.0000000     -0.1662423   -0.0000000     -0.0230903   -0.0000000      0.3809779    0.0000000      0.0116748    0.0000000     -0.0009808   -0.0000000      0.0025175  
  -0.0411286    0.0000000     -0.0075498    0.0000000      0.1479968    0.0000000      0.1421482    0.0000000     -0.1908940    0.0000000     -0.0085177    0.0000000      0.0116748   -0.0000000      0.4621919   -0.0000000     -0.0110672    0.0000000      0.0006904  
   0.0539939    0.0000000      0.0233133   -0.0000000     -0.0424905    0.0000000     -0.0638839   -0.0000000      0.1094664   -0.0000000      0.0382409   -0.0000000     -0.0009808    0.0000000     -0.0110672   -0.0000000      0.4997230    0.0000000     -0.0016861  
  -0.0210488    0.0000000     -0.0325156   -0.0000000     -0.0063524    0.0000000      0.0138251    0.0000000      0.0033805    0.0000000     -0.0075248    0.0000000      0.0025175   -0.0000000      0.0006904   -0.0000000     -0.0016861   -0.0000000      0.5391607  
  -0.1167733   -0.0000000      0.0972739    0.0000000      0.0454683   -0.0000000      0.0113219    0.0000000     -0.0069194   -0.0000000     -0.0116659   -0.0000000      0.0046546    0.0000000      0.0255121    0.0000000     -0.0044614   -0.0000000     -0.0014033  
   0.0826389   -0.0000000      0.0557496   -0.0000000     -0.0676305   -0.0000000      0.0312384   -0.0000000     -0.1662507   -0.0000000      0.0046007   -0.0000000     -0.0043375    0.0000000     -0.0221338    0.0000000      0.0214008   -0.0000000     -0.0034290  
   0.0411290   -0.0000000     -0.0075684    0.0000000     -0.1479888   -0.0000000      0.1421394   -0.0000000     -0.1908945    0.0000000      0.0254995   -0.0000000     -0.0221219   -0.0000000     -0.0500560    0.0000000      0.0175730    0.0000000      0.0020463  
  -0.0540128   -0.0000000      0.0233154    0.0000000      0.0424935    0.0000000     -0.0638766    0.0000000      0.1094696    0.0000000     -0.0044630    0.0000000      0.0213845    0.0000000      0.0175725    0.0000000     -0.0060172   -0.0000000     -0.0020753  
   0.0210495   -0.0000000     -0.0325247   -0.0000000      0.0063521   -0.0000000      0.0138209   -0.0000000      0.0033808    0.0000000     -0.0014127   -0.0000000     -0.0034195    0.0000000      0.0020470   -0.0000000     -0.0020772    0.0000000      0.0010964  
Now the corresponding density matrix entry at (1,6) and (1,11) as well as at (6,1) and and (11,1), corresponding to Ni(xz+yz):Ta(xy), becomes equal to 0.1168346. More precisely, the (1,6) and (6,1) component is 0.1168346, and (1,11) and (11,1) is -0.1168346. With V=3eV, this requires self-energy of 0.35eV, i.e., our input value for the self-energy.

We note that such non-local Coulomb interaction is finite, but rather small, and the value of 3eV is somewhat too large to be realistic. Consider the Coulomb repulsion on Ta atom, which we estimated to be U=6eV and J=0.8eV. This translates to Yukawa screening of lambda=0.408/rB, and dielectric screening epslion=1.260. This can be simply obtained by running
RCoulombU.py -U 6. -J 0.8
which prints the value of the screening parameter lambda and dielectric constant epsilon for the choosen values of U and J on Ta. The corresponding Yukawa form interaction is
V(r) = 2Ry \frac{\exp(-\lambda\,r)}{\epsilon\, r}
and a simplified estimation of the non-local Coulomb repulsion at distance 5.22 bohr radius (corresponding to 2.76A trimer distance) is 0.5eV. Taking the value of U=8eV on Ni will make such non-local repulsion even smaller. However, such non-local correlations can be generated by purely local interaction, when considering higher order processes, or equivalently, cluster-DMFT. We should thus consider such non-local Hartree-Fock calculation as a qualitative estimation of the non-local effects, which are currently hard to compute by the available impurity solvers.

Now that we have required non-local self-energy, we switch back to real axis, and plot the spectral function. We first copy the saved self-energy
cp sig.inp_real_axis sig.inp
and correct its header, i.e., s_oo=[....,0.35,-0.35]. Next we change the flag matsubara to 0 in Ta2NiSe5M.indmfl file. Than we copy the k-path file from previous calculation to the current directory ../onreal/Ta2NiSe5M.klist_band. And finally we rerun
x_dmft.py lapw1 --band
x_dmft.py dmftp
and display results with
wakplot.py
We see that the gap has been substantially increased, and is now much closer to the experimental value of approximately 0.3eV.

We can now plot also the orbitally resolved spectra as above, using dmftgk tool, which gives the following plot:

Structural optimization with and withouth non-local Hartree Fock self-energy

We can relax internal parameters of the monoclinic crystal structure within eDMFT. It turns out that the monoclinic distortion is not stable, and the internal structural parameters closely resemble those of the high-temperature orthorhombic phase (even though the unit cell has monoclinic distortion).

Now we can also check if the broaken symmetry, which we guessed above with the non-local Hartee Fock, stabilizes the monoclinic phase. Note that the symmetry of the Ta-Ni self-energy corresponds to B3g irrep, which is indeed experimentally determined symmetry of the order parameter.
Indeed this order parameter stabilizes the monoclinic phase, and the Ta-Ni bonds disproportionate into long and short.

To do this type of calculation, we will create yet another directory relax, and copy data from previous directory nlHF the input files
dmft_copy.py  ../nlHF/
cp ../nlHF/sig.inp_imag_axis sig.inp 
We will also set flag matsubara to 1 in Ta2NiSe5M.indmfl.
We next edit Ta2NiSe5M.in2 file, and set TOT to FOR, i.e., the first line now looks like
FOR             (TOT,FOR,QTL,EFG,FERMI)
and we will also change the mixer to MSR1a, which mixes both the charge density and atomic coordinates for relaxation of the structure. The first line of the Ta2NiSe5M.inm will look like
MSR1a   0.0   YES  (BROYD/PRATT, BG charge (-1 for core hole), norm)
It is advisable to compute an approximate starting Hessian by wien2k program pairhess:
x pairhess -f Ta2NiSe5M
before submitting the job.

At the moment the run will still fail because case.indmfi and case.indmfl files are incompatible. We created cluster in case.indmfl, but we did not appropriately update case.indmfi. We will edit case.indmfi such that it will contain only one impurity problem (actually cluster problem now), and identical Sigind matrix as found in case.indmfl. Our new case.indmfi looks like:
1   # number of sigind blocks
15  # dimension of this sigind block
 1  0  0  0  0 -11 0  0  0  0 -12 0  0  0  0
 0  2  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  3  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  4  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  5  0  0  0  0  0  0  0  0  0  0
-11 0  0  0  0  6  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  7  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  8  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  9  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0 10  0  0  0  0  0
-12 0  0  0  0  0  0  0  0  0  6  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  7  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  8  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  9  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0 10
We will add an option
runIMP         =  False
into params.dat file, because we do not have cluster-DMFT solver to compute self-consistent cluster self-energy. We will keep the single-site self-energy + Hartree Fock value for now, and relax structure by optimizing charge density. We note that the band gap in this system is not opened by correlations, but it is a band gap, in which the DMFT self-energy is quite insensitive to small changes of the structural parameters, as seen by structural relaxation of the single-site DMFT calculation.

Now we are ready to submit the job.

After a few hundred iterations, the forces are reduced below 0.3, with converged charge, and relaxation confirms the disproportionate bonds, namely, the two inequivalent Ta-Ni distances are
theory:
Ni-Ta1 = 2.759 Å
Ni-Ta2 = 2.871 Å
experiment:
Ni-Ta1 = 2.763 Å
Ni-Ta2 = 2.848 Å
Moreover, all other Wickkoff positions are in reasonable agreement with experiment, i.e,
theory:
Ni = [0,3/4,0.19482]
Ta = [0.51826406,0.60978623,0.21383811]
Se1= [1/2,3/4,0.32546859]
Se2= [0.00089185,0.55309075,0.36032579]
Se3= [0.00502183,0.64416589,0.06569934]
experiment:
Ni = [0,3/4,0.20121]
Ta = [0.512048,0.610655,0.221366]
Se1= [1/2,3/4,0.327557]
Se2= [0.991386,0.548924,0.354601]
Se3= [0.991752,0.638047,0.079931]

Optical conductivity

For optics calculation we create anothe directory optics and copy the real-axis data into the new directory:
dmft_copy.py ../nlHF/
We will increase the k-mesh, as converged optics calculation needs very fine mesh. For finite plot below, we will use 5000 k-points. But when testing the input parameters, much smaller mesh should be first attempted.
x kgen -f Ta2NiSe5M
   5000,
for this purpose you can create mpi_prefix.dat for parallel execution.
x lapw0 -f Ta2NiSe5M
x_dmft.py lapw1
Calculate dos by
x_dmft.py dmft1
Move the chemical potential into the middle of the gap (by changing EF in EF.dat) and recalculate the density of states. Note that for (undoped) insulators, the fermi level should be close to the center of the gap at low temperature.
Next print projector, but first remove parallelization, because projector should be stored in a single file:
rm mpi_prefix.dat
x_dmft.py dmftu
We now have Udmft[updn].0 and BasicArrays.dat.

Next, prepare the input file for computing the optics matrix elements within wien2k. The input files should be Ta2NiSe5M.inop:
99999 1       number of k-points, first k-point
-5.0 3.0 9999 Emin, Emax for matrix elements, NBvalMAX
2             number of choices (columns in *outmat): 2: hex or tetrag. case
1             Re xx
3             Re zz
ON            ON/OFF   writes MME to unit 4
Note the last line has flag ON, which is different than in standard DFT calculations. Namely, we require not only the diagonal, but also the off-diagonal components of the momentum operator in band representation.

Next run
x optic[c] [-so] [-up]  -f Ta2NiSe5M
This computes the matrix elements of the momentum operator. They appear in case.mommat2 and case.symop file.

Next we prepare an input file for dmft-optics calculation dmftopt.in, which should look like
    0          # Temperature  (this is new in 2014)
    Ta2NiSe5M  # case
    0.002      # gamma  -- broadening of all bands / we recommend to keep nonzero value
    0.0        # gammac -- broadening of the correlated bands (in addition to the self-energy)
    4.0        # ommax  -- maximum frequency for the optics calculation
    1e-2       # delta  -- minimum separation of frequency for logarithmic mesh in frequency
    50         # Nd     -- number of points in each linear mesh, a subset of logarithmic mesh
    T          # Qsym: [F|T]. Do we need to symmetrize? (over all k-points or just irreducible.) If Qsym is F, it goes over irreducible, if Qsym is T, it goes over all.
    F          # InterbandOnly [F|T] (F -- all, T--interband)
    5          # dwindow -- Not all bands are used in computation of optics, but only those from [-omeg-dwindow, omega+dwindow]. You should increase this to make sure you are not cutting some bands which are needed.
    3          # Ndirection -- How many optics type do you want to compute: xx, yy, zz,...
    1.0 0.0 0.0    # alphaV(1,:)
    0.0 0.0 0.0    # alphaV(2,:)
    0.0 0.0 0.0    # alphaV(3,:)
    0.0 0.0 0.0    # alphaV(1,:)
    0.0 1.0 0.0    # alphaV(2,:)
    0.0 0.0 0.0    # alphaV(3,:)
    0.0 0.0 0.0    # alphaV(1,:)
    0.0 0.0 0.0    # alphaV(2,:)
    0.0 0.0 1.0    # alphaV(3,:)
If temperature is not very high (like in the Earth's core) we can safely set temperature to 0, and we will still get very accurate results, corresponding to the temperature at which DMFT was run, because we have scattering rate from DMFT calculations, which gives finite width to the Drude peak. The temperature in optics calculation is just the Fermi function temperature in calculating the bubble, and is unfortunately quite expensive to use, i.e., the step function needs to be replaced by dense mesh of points near the step, and for most cases does not make noticable difference.

Here we choose very small broadening gamma=0.002, which will produce very sharp features in the optical spectra, and requires very many k-points. In most calculations, we could increase gamma for factor of 5, and reduce number of k-points, and still resolve most of the features. We will broaden the spectra at the end, with the tool broad, which has very similar effect as increasing gamma. Note however, that we need finite gamma to not loose weight in the integration over frequency, which is done analytically.

The frequency mesh, which is used in the optics calculation for the internal frequency integration, is a hybrid of equidistant and logarithmic mesh, sketched in this figure:
For very low frequency, where we expect Drude peak in metals, we use very fine mesh, which starts with distance delta (here 1e-2). In many metals at low-T one needs to reduce this parameter, to resolve the Drude peak. Next, we use 2*Nd+1 (here 2*50+1=101) points in equidistant mesh, to resolve every frequency between 0 and Nd*delta with high accuracy. In the next step, we keep only every second point of the previous mesh, and extent it to 2*Nd*delta, using again 2*Nd+1 points. We keep increasing the frequency intervals until we reach ommax (here 4), which is reached when delta*Nd*2**N0==ommax. N0 here determines number of levels at which we double the extent of the mesh. The mesh with Nd=50 is very precise, and in many situations one could use Nd=10 or Nd=5, and still get very reasonable optics spectra. We suggest one to experiment with these parameters when number of k-points is not too large, and cutoff frequency ommax is not too large, and than increase number of k-points once a good set of parameters is found.

Qsym should be set to true whenever we need more than just the average optical conductivity. Only for the average conductivity (sigma_xx+sigma_yy+sigma_zz)/3. it is sufficient to consider the irreducible k-points only. While this symmetrization works quite well in most cases, we advice one to check it in any particular low symmetry calculation, by using calculation with all k-points in the first Brillouine zone (not just irreducible) and seeting Qsym to F. The symmetrization is done in similar way as in wien2k optics calculation.

InterbandOnly is normally set to F, because within DMFT it is hard to split the interband contributions from Drude peak, and even though the equation is straighforward, the result might be unexpected.

Finally, we can calculate several directions of optics simultaneously. Here we choose three directions, sigmaxx, sigmayy and sigmazz. Any other entry in the matrix is possible, in principle, but has not been well tested.

Now we can remove large vector file, which is hard to copy when submitting the job, and proceed only with the essential files. The optics calculation needs the following files: case.mommat2, case.symop, case.symmat, case.energy[so][updn], sig.inp?, Udmft.0, BasicArrays, and dmftopt.in.

We can run the executable dmftopt in parallel, for example
mpirun [-n nkp] dmftopt > nohup.dat
It is advisable to compare density of states, calculated by optdos, with the previously calculated dos (case.cdos). The comparison might look like:
Notice that because we use here very small broadening, the gap is very clearly developed, while in dmft1 the gap was somewhat less well developed due to finite broadening.

The calculated optical conductivity is found in optics.dat, and should look like:
Notice that the calculation with 2000 k-points and 5000 k-points is similar, newertheless the details of the peaks are somewhat different in the two calculations. The 2000 k-point calculation seems to have s small peak in the gap, which completely dissapears when using 5000 k-points. Finally, it is advisable to broaden the sharp features from the calculated optical conductivity, which has similar effect as increasing gamma in the calculation, except that in postprocessing we can apply broadening many times, until we find reasonably smooth curve. For example, with gaussian convolution of width 0.02eV
broad -w 0.02 optics.dat > b_optics.dat
we get the following curves:
The first sharp peak is clearly much better developed at 5000 k-points than with 2000 k-points. Finally, the other two directions have less low energy features, and are easier to converge:
Comparison with experiment (from PRB 95, 195144 (2017)) reveals very reasonable agreement of all three directions:
Note that the high-frequency peaks are somewhat red-shifted, which is quite common feature of dft-type calculations, as the high-frequency orbitals are somewhat less precisely predicted by dft, and are not corrected by dmft.