# Tutorial 2: FeSe and optimization of the structure

## Tutorial 2.a. : Preparing the job

This tutorial is on FeSe, one of the Fe superconductors. We will first give tutorial on computing a single DFT+EDMFT calculation. In the later sections we will show how to optimize the structure.

FeSe cristalizes in P4/nmm with lattice parameter a=7.121433 a.u.. In our structure file, two iron atoms are at position (1/4,3/4,0) and at (3/4,1/4,0), while Se atoms are at (3/4,3/4,1-z) and (1/4,1/4,z), respectively. According to experiment, the internal parameter z is close to 0.27 (see PRB 79, 014522 (2009)) while in DFT is 0.24 (GGA) or 0.23 (LDA). In this tutorial we want to see how much dynamic correlations on iron change this parameter.

The correlations are on Fe ion, which contains around 6 electrons. These electrons are forming fluctuating spin, which is at finite temperature relatively poorly screened in the solid. While the charge fluctuations are very well screened and behave just as in a Cu metal, the orbital fluctuations are somewhat slower and result in enhanced orbital susceptibility. However, the spin degrees of freedom are really slow due to Hund's coupling, which penalizes low spin configurations, therefore spins and very poorly screened, hence the Fermi liquid scales can become very small, and the material is a poor metal. This physics was called Hund's metal physics and was first discovered here.

We start by performing DFT calculation to get good charge density. We create a directory with the name FeSe and copy the FeSe.struct file into it. [Note that Wien2k requires to call directory with the same name as the structure file]. We choose to start with z=0.26 (which is slightly smaller than experimental value), in order to be slightly away from the theoretical minimum.

We initialized the DFT by
```  x init_lapw
```
and choose the default options, except we will use LDA instead of GGA because LDA+DMFT gives better energies than GGA+DMFT. This is because in LDA the energies and the band structure are more consistent than in GGA. One way to understand this is to realize that both LDA and GGA give too broad bands and quasiparticle-renormalization Z=1. Therefore both should overbind, and consequently give too small lattice parameters in moderately correlated systems. As is well known, LDA indeed gives too small volumes, while GGA tends to give too large volumes in weakly to moderately correlated materials. As DMFT is expected to reduce binding and increase volumes in most material (as Z is always < 1), the LDA+DMFT is expected to be universally applicable, while GGA+DMFT is definitely wrong in moderately correlated systems, and experience shows that it does not get better in strongly correlated systems. In order to correct energies by DMFT, it is thus better to start with a consistent theory, such as LDA.
In version 14, we will choose the following options:
```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
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 : 500
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
```
Here Ctrl-X-C stands for breaking out of your editor. Ctrl-XC is used in emacs and similar editors, but in vi, one should replace Ctrl-X-C by :q!.

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.] In DMFT calculation, one needs converged frequency dependent local Green's function, hence more k-points is typically needed. We will increase the number of k-points to 2000 in this tutorial. To do that, run
`x kgen`
and specify 2000 k-points with shifted mesh. It is also advisable to increase the reciprocal cutoff of the DFT basis. Namely, in LAPW we form the basis with the basis functions that have envelope form of plane waves. The most oscillating plane wave in the basis is chosen according to criteria K*RMT < RKmax, where RMT is the smallest muffin-thin radius, and RKmax is a parameter in FeSe.in1, by default 7. We will change it to 9, so that the basis becomes much more precise (but calculation is slower). If you are limited with computer power, you might want to keep RKmax=7.

Next, rerun wien2k on this k-mesh, i.e., run_lapw. (This step is not essential as the DFT charge will not be very precise anyway for our DFT+DMFT calculation).

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 4 atoms in the unit cell:
1 Fe  2+
2 Fe  2+
3 Se  2-
4 Se  2-
Specify correlated atoms (ex: 1-4,7,8): 1-2
You have chosen the following atoms to be correlated:
1 Fe  2+
2 Fe  2+
Do you want to continue; or edit again? (c/e): c

For each atom, specify correlated orbital(s) (ex: d,f):
1 Fe  2+: d
2 Fe  2+: d
You have chosen to apply correlations to the following orbitals:
1  Fe  2+-1 d
2  Fe  2+-2 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  Fe  2+-1 d: 2
2  Fe  2+-2 d: 2
Do you want to continue; or edit again? (c/e): c

Specify projector type (default = 2):
------  ------------------------------------------------------------
Projector  Drscription
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
Each set of equivalent correlated problems are listed below:
1   (Fe  2+1 d) (Fe  2+2 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): i

Ctrl-X-C

Is this a spin-orbit run? (y/n): n

Ctrl-X-C
```

The init_dmft.py script generates two necessary input files, case.indmfl and case.indmfi. They connect the solid and the impurity with the DMFT engine, respectively. Let us take a closer look at the header of the FeSe.indmfl file:

```19 44 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)
2                                     # number of correlated atoms
1     1   0                           # iatom, nL, locrot
2   2   1                           # L, qsplit, cix
2     1   0                           # iatom, nL, locrot
2   2   2                           # L, qsplit, cix
```
The first line specifies the hybridization window. We set it to "(-10eV,10eV)" around EF, which it turns out, contains bands numbered from 19-44. For each k-point we will use the same bands, i.e., we will follow the bands with the projector.
The same line also specified the projector type, which is 5.
The second line contains a switch for the real-axis or imaginary axis calculation. If we require the local green's function or density of states on real axis, we will set the switch "matsubara" to 0. Here, we will work on imaginary axis, hence "matsubara" is set to 1. The second line also contains information on small lorentzian broadening in k-point summation. The last three numbers are used for plotting the spectral functions, to which we will return latter.
The third line specifies the number of correlated atoms, here two irons.
The forth line specifies which atom (from case.struct) is correlated. These are the first two Fe. Note that each atom in the structure file counts, even when several atoms in structure files are equivalent.
The forth line also specifies how many orbital indeces (different L's) we will use for projection, and calculation of local green's function and partial dos. The last index "locrot" is used for rotation of local axis, and is not needed here (set to 0).
The fifth line specifies orbital momentum "l=2", "qsplit=2" which we chose during initialization. Each correlated orbital-set gets a unique index during initialization. The first and second Fe atoms were given index "cix=1" and "cix=2", respectively. For each correlated set, a more detailed specification is needed, which is given in the remaining of "case.indmfl" file.

The specification of correlated set is given in the second part of case.indmfl:
```#================ # Siginds and crystal-field transformations for correlated orbitals ================
2     5   5       # Number of independent kcix blocks, max dimension, max num-independent-components
1     5   5       # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'z^2' 'x^2-y^2' '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.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.70710679 0.00000000   0.00000000 0.00000000   0.00000000 0.00000000   0.00000000 0.00000000   0.70710679 0.00000000
```

The first line specifies the number of all correlated blocks, the dimension of the largest correlated block, and number of orbital components which differ. The second line specifies the dimension of the first correlated block (here "d" has dimension 5), and the number of orbital components which differ.

The transformation matrix from spheric harmonics to the DMFT basis is given next, and shows how these orbitals are defined: It is the expansion of each of these orbitals in spherical harmonics. The columns correspond to the following quantum numbers (m_l=-2,-1,0,1,2) and rows to ('z^2','x^2-y^2','xz','yz','xy'). For example, the 'z^2' orbital, which is expanded in the first line, is equal to the m=0 spherical harmonic, whereas the yz orbital on the forth line is an equal superposition of m=1 and m=-1, multiplied by i (note that two consecutive real numbers define one complex number). Note that the transformation matrix must be unitary.

The same information (except the first line above) is repeated once more for the second Fe atom.
This file FeSe.indmfl can be edited in text editor and can be adjusted if necessary.

The second file, created at the initialization, is case.indmfi and contains the following lines:
```1   # number of sigind blocks
5   # dimension of this sigind block
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```
Notice that it contains only one impurity problem, while FeSe.indmfl contains two correlated blocks. Due to the symmetry, the index table ("Sigind") relates both correlated blocks with the impurity problem, hence we do not need to run the impurity solver twice.

Now we want to start DMFT 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. These include the .struct file as well as the various Wien2K input files such as .inm, .in1, etc. The charge density file, .clmsum, is also copied. This is all we need from the initial Wien2K run, and now we proceed to create the rest of the files, needed by the DMFT calculation.

Copy the params.dat file. This file contains information about the impurity solver and additional parameteres, which appear in DFT+EDMFT. Its content is:
```solver         =  'CTQMC'   # impurity solver
DCs            =  'nominal' # 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             =  2e-6      # the charge density precision to stop the LDA+DMFT run
ec             =  2e-6      # 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.

# Impurity problem number 0
iparams0={"exe"                : ["ctqmc"              , "# Name of the executable"],
"U"                  : [5.0                  , "# Coulomb repulsion (F0)"],
"J"                  : [0.8                  , "# Coulomb repulsion (F0)"],
"CoulombF"           : ["'Ising'"            , "# Can be set to 'Full'"],
"beta"               : [50                   , "# Inverse temperature"],
"svd_lmax"           : [25                   , "# 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"                : [200                  , "# 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"                : [6.0                  , "# Nominal occupancy nd for double-counting"],
}```
Note that the order of parameters is not important. Note also that this file is written in python syntax, hence it could be executed as python script. This could be used to check for possible errors by executing python params.dat. If no error is reported, the file is syntactically correct.

Below we explain the meaning of the variables
• `solver = 'CTQMC'`
The impurity solver will be Continous Time Quantum Monte Carlo, as implemented here.
• `DCs            =  'nominal'  # nominal double-counting `
In order to relax the structure, we choose nominal double-counting here, because it is much faster and quite accurate. For example, using DC='nominal' in relaxed structure, we will get nd=6.28 while using DC='exact' we get nd=6.25.
• ```max_dmft_iterations = 1
max_lda_iterations = 100
finish = 50```
Impurity solver is slower than DFT part, hence we will solve one impurity, and converge DFT in each global loop. The maximum number of global loops is set to 50. In this picture

we are choosing a single iteration for the red part, maximum 100 in the grey, and 50 iterations of both.
• `ntail         =  300`
The number of Matsubara points, distributed in logarithmic mesh, to describe the tail of Greens functions.
• ```cc             =  2e-6      # the charge density precision to stop the LDA+DMFT run
ec             =  2e-6      # the energy precision to stop the LDA+DMFT run```
The DFT loop will stop after it reaches the above given limit (max_lda_iterations). The loop will also stop when both the charge and energy are sufficiently converged (difference of charge less than cc, and energy less then ec).
• `recomputeEF = 1`
For metals, we must always re-calculate the chemical potential, so that the solid is charge neutral.
• ```wbroad         =  0.0      # broadening of sigma on the imaginary axis
This is important only when svd_lmax=0 and we sample directly in Matsubara basis. See Tutorial on MnO for more explanation.
• `"iparams0=...."`
For each impurity problem (in this case we have a single impurity) we need a python dictionary, containing parameters for the impurity solver. When more than single impurity is needed, one needs to specify "iparams1", "iparams2", ...
• `"U" : [5.0]`
The value of the on-site Coulomb repulsion (Slater F0). Note that in the charge self-consistent DMFT, one needs larger value of U than it is required for downfolded DMFT calculations, because part of the on-site screening (hybridization screening by many itinerant degrees of freedom) is explicitly taken into account in the such DMFT calculation. We will calculate the value of U in later tutorial.
• `"J"                  : [0.8]`
The Hund's coupling. Internally, we use Slater F2 and F4, but here we rather give more commonly used parameter "J", which is related to F2 = 14./1.625 J, and F4 = 0.625 * 14./1.625 * J. Notice that "J" can also be a list of two numbers, in which case F2=14./1.625 J[0] and F4 = 0.625 * 14./1.625 * J[1]. For more in-depth discussion, see the notes on the Coulomb repulsion.
• `"CoulombF" : ["'Ising'"]`
Here we use the density-density form of the Coulomb repulsion, which speeds up the calculation considerably. However, the Fermi liquid scale is different in the density-density and rotationally invariant calculation. We will adress this latter.
• `"beta"               : [50. ]`
Inverse temperature in eV. beta=50.0/eV means that temperature about 232 Kelvins.
• `"svd_lmax"           : [25 ],`
Instead of sampling self-energy in Matsubara basis, we will project it to the SVD basis, and determine 25 coefficients of expansion. This should give precision better than 1e-7.
• `"M"                  : [10e6]`
Number of Monte Carlo steps carried out on each processor. Notice that the total number of Monte Carlo steps in M*Ncores.
• `"mode"               : ["SH" ]`
We will calculate self-energy from the two particle Green's function, as explained here. The tail will be computed by the Hubbard I approximation.
• `"nom"                : [200]`
We calculate a finite number of Matsubara frequency points with sampling, while the rest are approximated by an analytic tail, which is much more stable. Typically, we need around nom ~ 3-4*beta points, which corresponds to an energy scale between 6*pi-8*pi ~ 22eV scale, outside which we can expect the analytic tail, computed up to the second order, to be an excellent approximation.
• `"nf0"               : [6.0 ]`
The nominal occupancy of the Fe-d shell, which is needed for double-counting. It is evaluated according to the formula:
and

• `"Ncout"              : [1000000]`
How often to print information of ctqmc solver to imp.?/nohup_imp.out.
• ```"GlobalFlip"         : [1000000              , "# How often to try a global flip"]
"warmup"             : [3e5                  , "# Warmup number of QMC steps"],
"tsample"            : [400                  , "# How often to record measurements"],```
Every "GlobalFlip"=1000000 MC steps we try a global flip, where we try to flip all up spin into down spins. We skip first "warmup" MC steps before we start measuring observables. Observables are being measured every "tsample" MC steps. Note that it is a good idea to look into imp.0/nohup_imp.out.000 to see how much time is spend in measuring the self-energy as compared to accepting the step. It is also a good idea to compare tsample with the Acceptance ratio printed in the same file.

There is only one more file that we need: A blank self energy (sig.inp) file. Generate it by the command
`szero.py`

## Tutorial 2.b. : Running the job

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. An example of creating such command is given by
echo "mpirun -np \$NSLOTS" > mpi_prefix.dat
or
echo "mpiexec -port \$port -np \$NSLOTS -machinefile \$TMPDIR/machines" > mpi_prefix.dat

but the precise form of the command depends on the supercomputer and its environment.
If many more CPU-cores are available than the number of k-points, we can optimize the execution further. Namely, a single k-point can be executed on many cores using open_mp instructions. To use this feature, one needs to specify "mpi_prefix.dat2" in addition to "mpi_prefix.dat" file. The former is than used in the dmft1, and dmft2 part of the loop, while the latter is used by the impurity solver. Clearly, "mpi_prefix.dat2" should specify the mpi command, which is started on a subset of machines. Therefore "-np" should be smaller than "\$NSLOTS" and "OMP_NUM_THREADS" should be larger than 1.
We provide a script that attempts to do that, hence you might want to insert the following line between "....> mpi_prefix.dat" and "run_dmft.py":

`\$WIEN_DMFT_ROOT/createDMFprefix.py \$JOBNAME.klist mpi_prefix.dat > mpi_prefix.dat2`
where "\$JOBNAME" should be set to wien2k case. Unfortunately, mpi execution is environment dependent, and therefore this script might not work on each system.

Before submitting the job, make sure that the computing notes have the following environmental variable specified: \$WIENROOT, \$WIEN_DMFT_ROOT, \$SCRATCH. Also make sure that \$WIENROOT and \$WIEN_DMFT_ROOT are in \$PATH as well as \$WIEN_DMFT_ROOT is in \$PYTHONPATH. Typically, the .bashrc will contain the following lines:
export WIENROOT=<wien-instalation-dir>
export WIEN_DMFT_ROOT=<dir-containing-dmft-executables>
export SCRATCH="."
export PATH=\$PATH:\$WIENROOT:\$WIEN_DMFT_ROOT
export PYTHONPATH=\$PYTHONPATH:\$WIEN_DMFT_ROOT
It might be necessary to repeat these commands (seeting variables) also in the submit script. This of ocurse depends on the system.

Here we paste an example of the submit script for SUN Grid Engine , but different system will require different script.
```#!/bin/bash
########################################################################
# SUN Grid Engine job wrapper
########################################################################
#\$ -N FeSe
#\$ -pe mpi2  250
#\$ -q <your_queue>
#\$ -j y
#\$ -M <email@your_institution>
#\$ -m e
########################################################################
source \$TMPDIR/sge_init.sh
########################################################################
export WIENROOT=<your_w2k_root>
export WIEN_DMFT_ROOT=<your_dmft_root>
export PATH=\$PATH:\$WIENROOT:\$WIEN_DMFT_ROOT
export PYTHONPATH=\$PYTHONPATH:.:\$WIEN_DMFT_ROOT
export SCRATCH='.'

JOBNAME="FeSe"   # must be equal to case name

mpi_prefix="mpiexec  -np \$NSLOTS  -machinefile \$TMPDIR/machines -port \$port -env OMP_NUM_THREADS 1 -envlist SCRATCH,WIEN_DMFT_ROOT,PYTHONPATH,WIENROOT"
echo \$mpi_prefix > mpi_prefix.dat

\$WIEN_DMFT_ROOT/createDMFprefix.py \$JOBNAME.klist mpi_prefix.dat > mpi_prefix.dat2

\$WIEN_DMFT_ROOT/run_dmft.py > nohup.dat 2>&1

rm -f \$JOBNAME.vector*
find . -size 0|xargs rm
```
The metallic FeSe is a bit harder to converge than the insulating MnO, as it typically takes a few more dmft steps. The impurity and local occupancy could be plotted from info.iterate file (column 9 and 10). By executing analizeInfo.py, we get a printout of only dmft steps, and if we plot the 9-th and 10-th column, we should see something like that:
After 30 or so dmft iterations (300 dft iterations), the occupancy converged to roughly 6.3.

If we let the job run for another 20 or so iterations (finish=50), we can get much better estimation of the total energy and free energy. By averaging over last 15 steps, i.e., steps 35-50, we can estimate the free energy of the system, which should be around -14802.540 (for RKmax=7) or -14802.570 (for RKmax=8) and statistical noise 1e-4 Ry (roughly 1meV) using 250 core computer.
Note that this is not the entire free energy yet. To compute the entire free energy, we need to subtract the entropy of the impurity, multiplied by temperature. We will learn how to compute entropy of the impurity in a subsequent tutorial, but for now it suffices to say that this term is very small at normal temperature as entropy is limited to a small fraction of log(2) in metallic systems and hence the term T*Simp is smaller than a meV.

## Tutorial 2.c.: Computing the Forces on atoms

We want to optimize the position of Se atom, therefore we want to see how large is force on Se atom. In order to compute forces on all atoms, we need to change only a single variable in FeSe.in2. [We can change that at the beginning of the run, or, during the run, before it converges, so that we have some statistics for the forces.] The first line in FeSe.in2 file should start with FOR rather than TOT, i.e,
```# FeSe.in2, line 1 :
FOR             (TOT,FOR,QTL,EFG,FERMI)```
After this change, FeSe.scf should contain a printout of the total force, for example:
```       TOTAL FORCE IN mRy/a.u. = |F|     Fx             Fy             Fz     with/without FOR in case.in2
:FOR001:   1.ATOM       0.000000       0.000000       0.000000       0.000000 total forces
:FOR002:   2.ATOM      19.938644       0.000000       0.000000     -19.938644 total forces
TOTAL FORCE WITH RESPECT TO GLOBAL CARTESIAN COORDINATES:
:FCA001:   1.ATOM                      0.000000       0.000000       0.000000 total forces
:FCA002:   2.ATOM                      0.000000       0.000000     -19.938644 total forces
TOTAL FORCE WITH RESPECT TO THE GLOBAL COORDINATE SYSTEM:
:FGL001:   1.ATOM                 0.000000000     0.000000000     0.000000000 total forces
:FGL002:   2.ATOM                 0.000000000     0.000000000   -19.938644357 total forces
```
This format is identical to the format in Wien2k, and it says that the Se atom has a force in z-direction of -19.94mRy/a.u. (for different RKmax, it might be somewhat different value). We can grep the force at each step by
`grep ':FGL002' FeSe.scf`
and its plot should look like
It shows the force on Se atom at each dft iteration. It is remarkable how stable is the force. It converges to 19.9 mRy/a.u. as soon as the occupancy is converged, and its noise is only around 0.1mRy/a.u. (when using 240 cores).

## Tutorial 2.d.: Structural relaxation

Once the calculation is converged, we can start structural relaxation. We can either resubmit the job again, or, just change a few variables while running the code.

1. First, we need to make sure that FeSe.in2 file starts with FOR, as described above, so that the forces are computed inside dmft2 step.

2. Next, make sure that max_lda_iterations in params.dat file is large enough, for example 100, so that dft will successfully relax structure at each dmft step.
3. Finally, we will change the FeSe.inm file, so that it starts with mixing method MSR1a, i.e.,
```# FeSe.inm, line 1:
MSR1a   0.0   YES  (BROYD/PRATT, extra charge (+1 for additional e), norm)  ```
In the next iteration, we should start seeing the following lines in FeSe.scf file:
```:APOS001  X=0.25000000 Y=0.75000000 Z=0.00000000 : DX=   0.00 DY=   0.00 DZ=   0.00 TOT=   0.00 mau
:APOS002  X=0.75000000 Y=0.75000000 Z=0.72948443 : DX=   0.00 DY=   0.00 DZ= -21.28 TOT=  21.28 mau

:FRMSA:  (mRyd/au)   73.175  103.485 RMS (au) 1.50E-02 MAX 2.13E-02 :F-condition (mRyd/au)    2.000 F
```
It prints the current positions of the two inequivalent atoms, and the attempted move of each. The third line shows the current condition to stop the relaxation of crystal structure. When that happens, the code continues to converge the charge at fixed crystal structure, so that at the end charge is very well converged.
The condition for stopping relaxation is the last number (2.000) and is called TRUST in wien2k. The code will attempt structural relaxation as long as forces are above TRUST. If all forces are smaller than TRUST for three consecutive dft iterations, and at the same time the charge and energy convergence are better than force_Rho_diff and force_Ene_diff, the structural relaxation is switched off.

We can change TRUST in FeSe.inM file (first line, second number). Let's reduce it to 0.5, because in such simples cases as FeSe we can achieve much better convergence.

We could also reduce force_Rho_diff, which is by default 0.01, or force_Ene_diff, which is 0.001. Both parameters can be added in params.dat file. Let's add parameter
`force_Rho_diff = 1e-4`
to params.dat file.

After some time passed, we can check how is the structural relaxation progressing:
`grep ':FRM' FeSe.scf`
should look like:
```:FRMSA:  (mRyd/au)  658.782  931.658 RMS (au) 2.47E-03 MAX 3.50E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)  622.141  879.841 RMS (au) 2.41E-03 MAX 3.41E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)  599.243  847.458 RMS (au) 5.76E-02 MAX 8.15E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)   73.175  103.485 RMS (au) 1.50E-02 MAX 2.13E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)   22.986   32.507 RMS (au) 3.19E-03 MAX 4.50E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)   14.839   20.985 RMS (au) 2.63E-03 MAX 3.72E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    4.120    5.826 RMS (au) 2.07E-03 MAX 2.92E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    1.355    1.916 RMS (au) 9.79E-03 MAX 1.38E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    9.552   13.509 RMS (au) 1.60E-02 MAX 2.26E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    5.092    7.201 RMS (au) 1.24E-02 MAX 1.75E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    7.584   10.725 RMS (au) 2.91E-03 MAX 4.11E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    3.919    5.542 RMS (au) 4.41E-03 MAX 6.23E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.782    1.105 RMS (au) 4.38E-03 MAX 6.20E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.759    1.073 RMS (au) 6.89E-03 MAX 9.75E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.727    1.028 RMS (au) 6.85E-03 MAX 9.68E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.089    0.126 RMS (au) 5.45E-03 MAX 7.71E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.649    0.918 RMS (au) 1.87E-03 MAX 2.65E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    1.281    1.811 RMS (au) 2.39E-03 MAX 3.38E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.475    0.671 RMS (au) 2.51E-03 MAX 3.55E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.991    1.402 RMS (au) 1.37E-03 MAX 1.94E-03 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.056    0.079 RMS (au) 1.14E-04 MAX 1.61E-04 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.046    0.066 RMS (au) 9.53E-05 MAX 1.35E-04 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.005    0.007 RMS (au) 1.51E-05 MAX 2.14E-05 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.425    0.600 RMS (au) 1.06E-05 MAX 1.50E-05 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.153    0.216 RMS (au) 1.14E-05 MAX 1.61E-05 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    3.312    4.684 RMS (au) 2.50E-04 MAX 3.53E-04 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)   13.764   19.465 RMS (au) 1.90E-04 MAX 2.68E-04 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    7.423   10.498 RMS (au) 2.81E-03 MAX 3.98E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.956    1.352 RMS (au) 3.48E-05 MAX 4.92E-05 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    1.352    1.912 RMS (au) 2.41E-03 MAX 3.41E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    1.693    2.394 RMS (au) 3.51E-03 MAX 4.96E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    2.032    2.874 RMS (au) 4.66E-03 MAX 6.60E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.223    0.315 RMS (au) 7.32E-03 MAX 1.04E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    1.224    1.731 RMS (au) 8.78E-03 MAX 1.24E-02 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.098    0.139 RMS (au) 2.79E-03 MAX 3.94E-03 :F-condition (mRyd/au)    2.000 F
:FRMSA:  (mRyd/au)    0.804    1.136 RMS (au) 7.93E-04 MAX 1.12E-03 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.001    0.001 RMS (au) 1.27E-04 MAX 1.79E-04 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.149    0.211 RMS (au) 1.80E-04 MAX 2.55E-04 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.001    0.001 RMS (au) 4.60E-05 MAX 6.50E-05 :F-condition (mRyd/au)    2.000 T
:FRMSA:  (mRyd/au)    0.010    0.014 RMS (au) 1.59E-05 MAX 2.25E-05 :F-condition (mRyd/au)    2.000 T```
We have TRUST at 2.0, and force_Rho_diff = 1e-4. The forces are very large at the beginning, but they quickly drop, and somewhere in the middle of the above output, the stop_condition is satisfied 6 times (line ends with "T"). Consequently, the relaxation was stopped at that point, and the charge density was converged to higher precision. But once the charge is converged, the impurity solver was called again, and it updated the self-energy. At this point, the structural relaxation is attempted again, and the force is again large (but considerably smaller than in the first step).
We can also print the current crystal structure parameter by
`grep ':APOS002' FeSe.scf`
As the printed parameter is 1-z, we want to print z directly, which can be achieved with a bit of transformation
`grep ':APOS002' FeSe.scf| sed 's/Z=/ /'|strans '\$j++' '1-#4'`
The plot of this z-parameter might look like
We started off with z=0.26, somewhat too small parameter. After a few hundred dft steps (30 dmft steps), we confidently converged to a value close to z=0.27, which is in much better agreement with experiment then GGA result (0.24). Of course the precise value of z depends on the value of Hund's coupling and also a bit on RKmax, and it could be slightly increased by increasing J or decreased by decreasing J. For physically relevant Hund's couplings (0.7-0.9eV), the z parameter is in the range 0.265-0.28.

## 2.e. Spectral function plots

Next we wanted to plot the spectral function. First we need self-energy on the real axis, therefore we will use once more the maximum entropy method on auxiliary Green's function (as explained in a previous tutorial), to analytically continue the self-energy. To reduce the MC noise, we average the self energy (sig.inp.*) files from the last few steps by
`saverage.py saverage.py sig.inp.2[5-9].1`
The output is written to sig.inpx. [We could append -o option, if we wanted to change the sig.inpx name.] Next, create a new directory (lets call it maxent), and copy the averaged self energy sig.inpx into it. Also copy the maxent_params.dat file, which contains the parameters for the analytical continuation process:
```params={'statistics': 'fermi', # fermi/bose
'Ntau'      : 300,     # Number of time points
'L'         : 20.0,    # cutoff frequency on real axis
'x0'        : 0.005,   # low energy cut-off
'bwdth'     : 0.004,   # smoothing width
'Nw'        : 300,     # number of frequency points on real axis
'gwidth'    : 2*15.0,  # width of gaussian
'idg'       : 1,       # error scheme: idg=1 -> sigma=deltag ; idg=0 -> sigma=deltag*G(tau)
'deltag'    : 0.005,   # error
'Asteps'    : 4000,    # anealing steps
'alpha0'    : 1000,    # starting alpha
'min_ratio' : 0.001,   # condition to finish, what should be the ratio
'iflat'     : 1,       # iflat=0 : constant model, iflat=1 : gaussian of width gwidth, iflat=2 : input using file model\$
'Nitt'      : 1000,    # maximum number of outside iterations
'Nr'        : 0,       # number of smoothing runs
'Nf'        : 40,      # to perform inverse Fourier, high frequency limit is computed from the last Nf points
}```
Perhaps the most important is mesh on the real axis, which is given by the upper cutoff L, the low energy cutoff x0 and Nw points. The mesh is non-uniform and more dense at low energy (tan-mesh). The MC error is here set to 0.004, and can be estimated from the variation in sig.inp.? files. Notice that we do not continue the self-energy, which would have very large error, and would not be stable to continue analytically. We will analytically continue an auxiliary function, as explained in a previous tutorial.
The script which does that is called maxent_run.py. We will thus run
`maxent_run.py sig.inpx`
which uses the maximum entropy method to do analytical continuation. This creates the self energy Sig.out on the real axis. It should be similar to this plot:
Notice that all orbitals have Fermi-liquid like behaviour at low energy, with substantial mass enhancement. The most correlated is the "xy" orbital. Notice also that the rotationally invariant calculation (CoulombF="Full") would give somewhat more coherent self-energy at low frequency, and should be preferred method of calculating spectra. We will show this below.

Now we need to make one last dmft1 calculation in order to obtain the Green's function and DOS on the real axis. Create a new directory inside your output directory, and while it the new directory, use
`dmft_copy.py ../`
to copy necessary files from the output of the DMFT run to the new directory. Also copy the analytically continued self energy Sig.out to sig.inp in this new directory:
`cp ../maxent/Sig.out sig.inp`
Since we need to run the code on the real axis, we need to change the NiO.indmfl file. The first number on the second line of NiO.indmfl file determines whether the code is run on the real axis (0) or the imaginary axis (1). Change it to 0. The header of FeSe.indmfl file should now look like this
```19 42 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)
2                                     # number of correlated atoms```
Then, inside the new directory run
```x lapw0 -f  FeSe
x_dmft.py lapw1
x_dmft.py dmft1
```
Before you execute these lines, you might want to set OMP_NUM_THREADS to a large number (but smaller then the number of cores on your node), such as 4, to parallelize some loops with open_mp.

Once the run is done, we now have the densities of states on the real axis. Both the total (column 2) and partial Fe-d (column 3 and 4) are stored in FeSe.cdos:

Next we create a desired k-path in the Brillouin zone. One can do that with the xcrysden tool, or, by editing a script klist3_generate.py. Let's create a path from Gamma, to X and to M and back to Gamma with 240 points in total with the name FeSe.klist_band.
Next, adjust the frequency range in FeSe.indmfl file. Suppose that we want to plot spectra in the window between -1 eV to 1 eV with 200 frequency points. The header file of FeSe.indmfl should then be corrected to
```19 42 1 5                             # hybridization band index nemin and nemax, renormalize for interstitials, projection type
0 0.025 0.025 200 -1.000000 1.000000  # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)
```
Notice that we changed only the last three numbers of the second row.

Next, in directory with real axis self-energy, run lapw on these k-points by executing
```x_dmft.py lapw1 --band
x_dmft.py dmftp
```
The first line computes the Kohn-Sham eigenstates on the new k-points from FeSe.klist_band. The second line computes the frequency dependent eigenvalues, and stores them to eigvals.dat. Now we can display the spectra by executing
```wakplot.py 1.0
```
Note that you might need to adjust brightness by reducing the argument to a number between 0.0-1.0, so that the excitations are clearly seen in the picture.
The spectra is extremely incoherent here, so that all excitations are washed out. This is because we used quite high temperature, and partially also because we used density-density form of the interaction.

### 2.e.a. Low temperature spectra and Orbitally resolved spectra

Next, we will reduce the temperature a bit, to recover the coherence of the spectral function. We will also use the rotationally invariant Coulomb interaction, which should give better results (less scattering at low energy, but higher effective masses). Note that the rotationally invariant Coulomb interaction never leads to spin-freezing at zero temperature. Only when the Ising interaction is used, one can get spin-freezing down to lowest temperatures.

To proceed, we will create a new directory, and while in that directory, copy results from previous dmft run into this new directory, i.e.,

`dmft_copy.py <dmft_previous_run>`
We will then change the following parameters in params.dat file:
```          "CoulombF"           : ["'Full'"             , "# Can be set to 'Full'"],
"beta"               : [100                  , "# Inverse temperature"],
"nom"                : [400                  , "# Number of Matsubara frequency points sampled"],
"M"                  : [20e6                 , "# Total number of Monte Carlo steps"],
"mode"               : ["GH"                 , "# We will use self-energy sampling, and Hubbard I tail"],
"PChangeOrder"       : [0.97                 , "# How often to add/rm versus move a kink"],
"warmup"             : [1e6                  , "# Warmup number of QMC steps"],
```
and resubmit the job. The crucial parameters to change are CoulombF and beta, however, to converge easier, we have adjusted several other parameters. We increased nom, because we increased beta, and we want to keep similar frequency cutoff for sampled self-energy (w=(2*nom+1)*pi/beta). We also increased the number of MC steps, because the acceptance rate is much smaller when Full interaction is used. We also changed mode to sample the Green's function rather than the self-energy. The latter still has some issues when the full Coulomb repulsion is used. In addition, we notice that using the default value for PChangeOrder=0.9, the acceptance ratio for adding/removing kinks is very small (less than 1%) while the acceptance for moving a kink is above 20%. For better statistics, it is a good idea to more often try to add/rm a kink. We also notice that the trial time is very short compared to acceptance time, hence we will not increase much the trial time by modifying PChangeOrder. We also increased the warmup, because the first DMFT step will produce very poor self-energy as it starts without a good starting configuration. In the next DMFT step, we can reduce warmup for an order of magnitude.

We will also change the mixing method back to MSR1 in FeSe.inm file, because the structure will not change appreciably anymore. Once we converge the solution with the full Coulomb repulsion, using previously determined optimized structure, the force should still be appreciably smaller than our previous cutoff (0.5mRy/a.u.).

The first iteration is not very good, because we are starting without status.xxx files, therefore the impurity problem starts from atomic configuration. But nevertheless we achieve rapid convergence of the self-energy. After a few dmft steps, we can also increase slightly the number of MC steps, to have better statistics:

`     "M"                  : [30e6                 , "# Total number of Monte Carlo steps"]`

The results are converged within 5-10 DMFT steps. To have some confidence in convergence, we will let the code run for a few more iterations. Here is the plot of total energy and free energy as a function of dmft iterations (up to 30 dmft iterations):
We have some MC noise, but the results seem converged after just two iterations. To obtain this plot (printing only the dmft steps), remember to run analizeInfo.py script.

However, to be confident in the convergence, it is important to also check how well is the self-energy converged. Here is the plot obtained on 250 cores:
Clearly, the convergence was not good until iteration 5 or so. But after iteration 10, the results are converged within the statistical noise, in particular at low energy, where most of the physics is.

Next we will perform the analytic continuation of the self-energy. First average self-energy over a few iterations to improve statistics:
`saverage.py sig.inp.2?.1 sig.inp.30.1`
The averaged self-energy sig.inpx should look like that:
Then copy sig.inpx into a new directory for analytic continuation, and copy maxent_params.dat from previous analytic continuation. Finally, run
`maxent_run.py sig.inpx`
to obtain Sig.out, which should look like that:
On this scale, the self-energy does not look very different from the previous self-energy at higher temperature (and Ising-type interaction), however, the imaginary part at zero goes to zero now, so that there is very little scattering left at low frequency. Note that this scattering rate can be read off from the imaginary axis data (sig.inpx) above.

We will now repeat several steps from above, to obtain the density of states on the real axis:
• create a new directory, for example on_real_axis
• while in the new director, copy converged results into the new directory by
`dmft_copy.py ../`
• in the new directory, overwrite the imaginary axis self-energy with the real axis counterpart
`cp ../maxent/Sig.out sig.inp`
• Change the second line of FeSe.indmfl file, so that the first flag (matsubara) is set to 0. Also change the frequency range for plotting to -1.0 to 1.0, which we will need in the next step. The second line should look like:
`0 0.025 0.025 200 -1.000000 1.000000  # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)`
• Obtain DOS by running lapw0, lapw1, and dmft1 steps:
```x lapw0 -f  FeSe
x_dmft.py lapw1
x_dmft.py dmft1```
Not surprisingly, the DOS looks quite similar to previous results,
except that the quasiparticle peak is considerably sharper.

Next we copy FeSe.klist_band from previous calculation into the current directory, and re-run Kohn-Sham diagonalization on the selected k-path by
`x_dmft.py lapw1 --band`
followed by computation of the DMFT eigenvalues
`x_dmft.py dmftp`

Finally display results by
`wakplot.py 0.1`
Notice that the intensity (0.1) should be tuned until a nice plot like that
is obtained.

Now we want to resolve the spectra in its orbital components, i.e., xz, yz, and xy orbital. To do that, we will use another executable, which can print out both the eigenvalues and the eigenvectors of the Dyson equation.
First, we will print DMFT projector, which relates the Fe-orbitals with the Kohn-Sham orbitals.
`x_dmft.py dmftu -g --band`
Note that -g stands here for printing only the specified k-points in the klist, but not their star (equivalent k-points). Moreover --band switch takes FeSe.klist_band rather than FeSe.klist.
The self-energy on real axis, with the desired frequency range, should normally be prepared here. However, since we already run x_dmft.py dmftp above, we should have files like sig.inp1_band already present.

Next, we will prepare an input 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
FeSe.energy         # LDA-energy-file, case.energy(so)(updn)
FeSe.klist_band     # k-list
FeSe.rotlm          # for reciprocal vectors
Udmft.0             # filename for projector
0.0025              # gamma for non-correlated
0.0025              # gammac
sig.inp1_band  sig.inp2_band   # self-energy name, sig.inp(x)
eigenvalues.dat     # eigenvalues
UR.dat              # right eigenvector in orbital basis
UL.dat              # left eigenvector in orbital basis
-2.                 # emin for printed eigenvalues
2.                 # emax for printed eigenvalues
```
Note that this input file is a normal text file read by fortran, and the order in which parameters appear is important here. Most of the parameters are self-explanatory. The files like BasicArrays.dat and Udmft.0 and FeSe.rotlm were obtained before when printing projector (x_dmftp.py dmftu). Further FeSe.energy file was obtained by running x_dmft.py lapw1. Furthermore sig.inp1_band and sig.inp2_band were obtained before when running x_dmft.py dmftp step. Note that the number of files in this line (line number 10) needs to match the Number of independent kcix blocks from FeSe.indmfl file. In the case of FeSe, we have two Fe-atoms, and for each we needs to specify the self-energy, although the two self-energies are equal. Finally, the output eigenvalues will be written to eigenvalues.dat and the left and the right eigenvectors will appear in UL.dat UR.dat, respectively. Note that we will print only the projection of the eigenvectors to the orbital basis.
The last two lines contain the range in which the eigenvalues/eigenvectors will be printed. Once you have dmftgke.in file, you can execute
`dmftgk dmftgke.in`
which prints out the eigenvalues and eigenvectors. You should see that eigenvalues.dat, UL.dat and UR.dat appeared.

Finally, copy the python script
`wakplot_sophisticated.py`
to your current directory, and edit it. Scroll down to the main part of the script, which should look something like this
```if __name__ == '__main__':
aspect_scale=0.3
small = 1e-5
DY = 0.0 #
intensity = [3.,3.,3.,1.]
COHERENCE=True
orb_plot=[1,1,2,2,1,1,0,0,0,0,0,0,0,0] #2,2,2,2,2,2]  # should have integers 0,1,2 for red,green,blue
#orb_plot=[]  # plotting unfolded band structure, but all orbits with hot map
```
To produce similar plot as before, tune the first number in the intensity list to similar number then you used above in wakplot.py. For example, setting
```intensity = [0.1,3.,3.,1.]
```
will make the plot brighter. If you have much larger number of k-points than frequency points, you might also need to tune aspect_scale, to make plot more/less elongated. For example
```aspect_scale=1.
```
will plot similar plot as wakplot.py above.

To display the plot, run
`./wakplot_sophisticated.py eigenvalues.dat`
Once we have a nice aspect ratio plot, we start adding orbital components.
First we will make two symbolic links to our eigenvectors
```ln -s UL.dat UL.dat_
ln -s UR.dat UR.dat_
```
This is because UL.dat_ might be in general different than UL.dat if we were attempting to unfold to a different size Brillouin zone. We will not do that here.
Next we will tune orb_plot (in wakplot_sophisticated.py) so that different orbitals are plotted with desired color. We have three color options, i.e., red, green, blue, which will be assigned numbers 0,1,2. We also have 10 orbitals in total, which corresponds to 5 orbitals of the first Fe atom and 5 of the second Fe atom, which appear in the following order [z^2, x^2-y^2, xz, yz, xy, z^2, x^2-y^2, xz, yz, xz]. This is of course specified in FeSe.indmfl file. We will assign red color to xy orbital, green to xz/yz orbital and blue to z^2 and x^2-y^2.
```#         z2, x2, xz, yz, xy, z2, x2, xz, yz, xy
orb_plot=[ 2,  2,  1,  1,  0,  2,  2,  1,  1,  0]
```
The intensity now stands for relative intensity of the three colors (RGB) and transparency (alpha). Lest set it to:
```intensity = [3.,3.,3.,1.]
```
to get
You might notice that the k-points are now displayed with Greek symbol Gamma, rather than the text GAMMA. You can achieve that by replacing word GAMMA in FeSe.klist_band with \$\Gamma\$, and M by \$M\$,....
As one can see, the inside hole pocket at Gamma is almost entirely of xz/yz character, the second pocket has a small admixture of z^2/x^2-y^2 to the xz/yz. The outside hole pocket at Gamma is almost entirely of xy character when it crosses the Fermi level, but it mixes with z^2/x^2-y^2 at x-point. The electron pockets are very yellow, hence equal mixture of xz/yz and xy character.

If we want to make the xy orbital a bit stronger, we might want to increase the first component of intensity, or better, decrease the second and third component of intensity. Furthermore, if we want to make the non-Fe orbitals more transparent, we might want to increase the last number for transparency. For example, the list
```#           red,green,blue,transparent
intensity = [3.,  2.,   2.,   3.]
```
gives the following plot: