The first step is to make a DFT calculation to get a good charge density. Copy the Ce_alpha.struct file into directory with name Ce_alpha, and run

init_lapw

initso

run_lapw -so

run_lapw -so

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 using

init_dmft.py

There are 1 atoms in the unit cell: 1 Ce Specify correlated atoms (ex: 1-4,7,8): 1

For each atom, specify correlated orbital(s) (ex: d,f): 1 Ce: f

Specify qsplit for each correlated orbital (default = 0): Qsplit Description ------ ------------------------------------------------------------ 0 average GF, non-correlated 1 |j,mj> basis, no symmetry, except time reversal (-jz=jz) -1 |j,mj> basis, no symmetry, not even time reversal (-jz=jz) 2 real harmonics basis, no symmetry, except spin (up=dn) -2 real harmonics basis, no symmetry, not even spin (up=dn) 3 t2g orbitals -3 eg orbitals 4 |j,mj>, only l-1/2 and l+1/2 5 axial symmetry in real harmonics 6 hexagonal symmetry in real harmonics 7 cubic symmetry in real harmonics 8 axial symmetry in real harmonics, up different than down 9 hexagonal symmetry in real harmonics, up different than down 10 cubic symmetry in real harmonics, up different then down 11 |j,mj> basis, non-zero off diagonal elements 12 real harmonics, non-zero off diagonal elements 13 J_eff=1/2 basis for 5d ions, non-magnetic with symmetry 14 J_eff=1/2 basis for 5d ions, no symmetry ------ ------------------------------------------------------------ 1 Ce-1 f: 4

Specify projector type (default = 2): Projector Description ------ ------------------------------------------------------------ 1 projection to the solution of Dirac equation (to the head) 2 projection to the Dirac solution, its energy derivative, LO orbital, as described by P2 in PRB 81, 195107 (2010) 4 similar to projector-2, but takes fixed number of bands in some energy range, even when chemical potential and MT-zero moves (folows band with certain index) 5 fixed projector, which is written to projectorw.dat. You can generate projectorw.dat with the tool wavef.py ------ ------------------------------------------------------------ > 5

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

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

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

Next create a new folder, and when in that folder, use the command

dmft_copy.py <lda_results_directory>

Next, copy the params.dat file, which has the content

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. # mix_delta = 0.5 # uncomment, if experience instability # Impurity problem number 0 iparams0={"exe" : ["ctqmc" , "# Name of the executable"], "U" : [6.0 , "# Coulomb repulsion (F0)"], "J" : [0.7 , "# Coulomb repulsion (F0)"], "nc" : [[0,1,2,3] , "# Impurity occupancies"], "beta" : [100 , "# Inverse temperature"], "svd_lmax" : [25 , "# We will use SVD basis to expand G, with this cutoff"], "M" : [10e6 , "# Total number of Monte Carlo steps"], "mode" : ["GH" , "# We will use self-energy sampling, and Hubbard I tail"], "nom" : [200 , "# Number of Matsubara frequency points sampled"], "tsample" : [300 , "# How often to record measurements"], "GlobalFlip" : [1000000 , "# How often to try a global flip"], "warmup" : [3e5 , "# Warmup number of QMC steps"], "atom_Nmax" : [100 , "# Maximun size of the block in generatic atomic states"], "nf0" : [1.0 , "# Nominal occupancy nd for double-counting"], }

The other important difference compared to previous calculations on
*d* systems is in two additional
parameters, which both control the exact diagonalization of the
atom. These parameters are specific to *f*-systems, and have no effect
in *d* systems:

nc=[0,1,2,3]specifies that we will take into account only a finite number of valences, namely Ce f

atom_Nmax=100specifies that if a block of atomic eigenstates has dimension more than 100, it should be cut to size 100. This does not have any effect in Ce, as the largest block in f

To proceed, one should create a blank self-energy on imaginary axis by a command

szero.py

Moreover, to allow the impurity to restart from the previous calculation,
one can also copy *imp.0/status.xxx* files, which contain the
last configuration of the impurity solver. We also provide
*status_1.tgz*, *status_2.tgz*, and *status_3.tgz*
files, which contain status files from previous steps, from one, two
and three steps before. This is useful when we want to restart from a
status from a few steps ago, in case the calculation goes in the wrong
direction.

Now we have all files prepared, and we can submit the job. The
submit-script should invoke python
script *"$WIEN_DMFT_ROOT/run_dmft.py"*. Do not forget to
prepare *mpi_prefix.dat* file before invoking the python
script.

Once the DFT+DMFT is running, monitor its status by checking the
log files

- at the to level:
*info.iterate*,*:log*,*case.scf*,*case.dayfile*,*dmft_info.out*. - dmft1 step can be monitored by checking
*dmft1_info.out*and*Ce_alpha.outputdmf1* - impurity can be monitored by checking
*imp.0/ctqmc.log.xxx*,*imp.0/nohup_imp.out.000**imp.0/Sig.out.xxx*and*imp.0/Gf.out.xxx* - dmft2 step is best monitored through
*dmft2_info.out*and*Ce_alpha.outputdmf2*

To check what is the distribution of perturbation orders, we can plot

We can also check what is the probability for each atomic state, i.e., probability that a Ce-f electron is in any of the atomic states. This information is written in

Here is the plot of probability for first 106 states:

In the next plots, we show the self-energy on imaginary axis for both alpha and gamma phase or Cerium. We can monitor self-energy during the self-consistent run. After 3 DMFT iterations, the self-energy of alpha-Cerium is almost converged, and we see that it does not change much between iteration 5-9. The slope of imaginary part of the the 5/2 self-energy is dSigma/dw~4, which gives mass enhancement of m*/m~5, i.e., the quasiparticle peak is roughly 5-times narrower than in DFT.

The gamma-phase of Cerium is in local moment regime, where the scattering rate at low energy is very large (bad metal). Here we see the imaginary part of 5/2 self-energy to reach 1.5eV.

To obtain the self-energy on the real axis, we need to do analytical continuation from Matsubara to real axis. We will use maximum-entropy method on auxiliary quantity G

saverage.py sig.inp.[5-9].1

maxent_run.py sig.inpx

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, and use

dmft_copy.py <dmft_results>

x lapw0 -f alpha-Ce x_dmft.py lapw1 x_dmft.py lapwso x_dmft.py dmft1

Finally, let us resolve the spectra in momentum space. For that we first need to prepare k-list. We can take one from Wien2k templated for fcc crystal structure. We then run DFT by

x lapw1 -f alpha-Ce -band x lapwso -f alpha-Ce -band

We also need to edit alpha-Ce.indmfl file to set the range in frequency and make sure that matsubara is set to 0. We will use 400 frequency points on real axis in the energy window -4eV to 5eV, hence the header of alpha-Ce.indmfl will look like that

9 34 1 5 # hybridization band index nemin and nemax, renormalize for interstitials, projection type 1 0.025 0.025 400 -4.000000 5.000000 # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)

Next, we compute DMFT eigenvalues by executing

x_dmft.py dmftp

wakplot.py

x_dmft.py dmftp 0.02