Note that we do not include spin-orbit coupling here for pedagogical reasons. For the real material it is necessary to include also SOC.
DFT (Wien2k) and Wannier orbitals
=================================
DFT setup
---------
First, we do a DFT calculation, using the Wien2k package. As main input file we have to provide the so-called struct file :file:`Sr2MgOs6_noSOC.struct`. We use the following:
The DFT calculation is done as usual, for instance you can use for the initialisation
init -b -vxc 5 -numk 2000
This is setting up a non-magnetic calculation, using the LDA and 2000 k-points in the full Brillouin zone. As usual, we start the DFT self consistent cycle by the Wien2k script ::
Note that, due to the distortions in the crystal structure, we need to include all five d orbitals in the calculation (line 8 in the input file above). The projection window is set such that all d orbitals are included.
Before starting the DMFT calculation it is beneficial to look a bit more closely at the block structure of the problem. Eventually, we want to use a basis that is as diagonal as possible, and include only the partially filled orbitals in the correlated problem. All this can be done using the functionalities of the :class:`BlockStructure <dft.block_structure.BlockStructure>` class, see section :ref:`blockstructure`.
We first initialise the SumK class::
from triqs_dft_tools.sumk_dft import *
SK = SumkDFT(hdf_file='Sr2MgOsO6_noSOC.h5',use_dft_blocks=True)
The flag *use_dft_blocks=True* determines, as usual, the smallest possible blocks with non-zero entries, and initialises them as *solver* block structure. In order to disentangle the :math:`d_{x^2-y^2}` and the :math:`d_{xy}` orbitals, and pick out only the partially filled one, we do a transformation into a basis where the local Hamiltonian is diagonal::
mat = SK.calculate_diagonalization_matrix(prop_to_be_diagonal='eal',calc_in_solver_blocks=True)
We can look at the diagonalisation matrix, it is::
This transformation is already stored in the SK.block_structure class. The next step is actually not needed for a DMFT calculation, but it is good to do this check to see what the transformation does to the local Hamiltonian. We can calculate it before rotation, rotate it, and look at the 2x2 block with and without off-diagonals::
print eal[0]['up'][1:3,1:3] # prints the 2x2 block with offiagonals
[[ 0.391+0.j -0. -0.815j]
[-0. +0.815j 4.884-0.j ]]
print eal2['up_1'] # prints the 2x2 block after rotation
[[0.247-0.j 0. +0.j]
[0. -0.j 5.028+0.j]]
So the local Hamiltonian has been diagonalised. From the other entries we can see that the *up_0* block and the [1,1] entry of the *up_1* block correspond to :math:`e_g`-like orbitals, and the others are the
We see that we arrived at a structure with 3 orbitals per spin only, and blocks of size 1x1.
The interaction Hamiltonian
---------------------------
We now set up the interaction Hamiltonian. Since we want to rotate the interaction matrix into the local basis, we are using the Slater convention for it::
In the last line we use the Wien2k convention to write the U matrix in the cubic harmonics. Next, we want to set up a Hamiltonian and rotate it into the *solver* basis::
Note that we needed to set up the interaction Hamiltonian for the full set of five *d* orbitals. The :meth:`convert_operator` method then takes care of rotating and picking the relevant orbitals. In the last line above we made the Hamiltonian real, since we know it this case that there are only real numbers in the interaction Hamiltonian. Note that this is not generally the case!
S = Solver(beta=beta, gf_struct=SK.block_structure.gf_struct_solver_list[0])
# Solver parameters:
p = {}
# solver
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 200
p["n_warmup_cycles"] = 100000
p["n_cycles"] = 1000000
# tail fit
p["perform_tail_fit"] = True
p["fit_max_moment"] = 4
p["fit_min_n"] = 40
p["fit_max_n"] = 100
The DMFT loop with automatic basis rotations
--------------------------------------------
After these initialisation steps, the formal DMFT cycle is very similar to a calculation without these basis rotations, since these rotations are done automatically, once the :class:`BlockStructure` property *transformation* is set, see :ref:`basisrotation`.
The DMFT loop itself looks very much the same as in :ref:`SrVO3`::
# double counting correction:
dc_type = 0 # FLL
# DMFT loops:
n_loops = 1
for iteration_number in range(1,n_loops+1):
mpi.report("Iteration = %s"%iteration_number)
SK.symm_deg_gf(S.Sigma_iw) # symmetrizing Sigma
SK.set_Sigma([ S.Sigma_iw ]) # put Sigma into the SumK class
chemical_potential = SK.calc_mu( precision = 0.01 ) # find the chemical potential for given density
The only difference to the other example is in the initialisation of the real part of the self energy. We cannot just take an element of the *dc_imp* array, since this array is stored in the *sumk* structure. Therefore, we first need to transform this matrix into *solver* space, and then take the appropriate matrix element. After the first iteration (here done with 24e6 MC sweeps and using the real version of the CTMQC solver), you should get self energies like this:
A complete python script for this tutorial, including some more input/output, is available (:download:`Sr2MgOsO6_noSOC.py <images_scripts/Sr2MgOsO6_noSOC.py>`). When running the script, you will encounter warnings during the transformation from the *sumk* to the *solver* basis. These warnings just reflect that the off-diagonal elements of the full Greens function are not zero at all frequencies, although the local Hamiltonian is. In that sense, we still do an approximation when restricting ourselves to the :math:`t_{2g}`-like orbitals.
