Tutorial 1: MnO

Tutorial 1.a. : Preparing the input files

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: 13
SELECT ENERGY : -6
Ctrl-X-C
continue with kgen : c
Ctrl-X-C
Ctrl-X-C
NUMBER OF K-POINTS : 500
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

Specify correlated atoms (ex: 1-4,7,8): 1
Do you want to continue; or edit again? (c/e): c
For each atom, specify correlated orbital(s) (ex: d,f):
 1 Mn: d
Do you want to continue; or edit again? (c/e): c
Specify qsplit for each correlated orbital (default = 0):
 1 Mn-1 d: 7
Do you want to continue; or edit again? (c/e): c
Specify projector type (default = 5): 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
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? (i/r) (default=i) : i
Ctrl-X-C
Is this a spin-orbit run? (y/n): n
Ctrl-X-C

• At the beginning, the code shows all atoms in the unit cell:

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

  and asks which one should be treated dynamically. We choose only the first atom, Mn, to be treated as a correlated atom in the DMFT step.
• The code next asks which orbital quantum numbers should be treated dynamically:

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

  The partially filled d shell of the Mn ion is to be treated in DMFT.
• Next we should specify what type of local symmetry we have on the given atom, and if we possibly want to correlated only a sub-shell

  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 Mn-1 d: 7
  

  This chooses the local basis for the DMFT calculation. Notice that current implementation of the DMFT in combination with CTQMC is basis dependent. This is because by default we ignore off-diagonal components of the hybridization function. (The code also supports off-diagonal hybridization, but its use requires more advanced options, and leads to a minus sign problem in the solution of the impurity.)
  In MnO, we ignored the spin-orbit coupling, and the crystal symmetry is cubic, hence hybridization is exactly diagonal in the cubic harmonics basis. We thus want to use this basis, which is the option 7. We could also use option 2, but this would lead to a slight worse statistics as impurity solver now averages over 3 t2g (or 2 eg) equivalent orbitals.
  

• Specify projector type (default = 5):
   Projector Drscription
  ------ ------------------------------------------------------------
   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
  

  As discussed above, DMFT requires a local basis, i.e., projector to local Green's function. If we could treat by DMFT all degrees of freedom on a given atom, the solution would not depend on the choice of the projector (it would depend on the range though). Fortunately, the orbitals on Mn ion are quite nicely separated in energy, and the only narrow states close to EF on Mn ion are the 3d states, hence they become correlated. Note that the Mn 4s orbital is not far above EF, but the 4s state is very spatially extended, hence weakly correlated.
  All choices of the projector project to spherical/cubic harmonics inside a sphere around an atom, but the radial dependence of the function is slightly different. "Projector 1" chooses solution of the Dirac equation for the radial part, "Projector 2" chooses a combination of the Dirac equation solution and its energy derivative, "Projector 4" has the same radial component as Projector 2, but when choosing an energy range for the projector, it will use equal number of bands at each k-point (projector will follow bands). Finally, projector 5 fixes the projector to the solution of the Dirac equation on the LDA level.
  If we are interested in the total free energy of the system, we should use projector 5, because in this case the projector is fixed during DMFT self-consistent run, which guarantees that the DFT+Embedded-DMFT method is stationary.
• Cluster versus single-site DMFT option:

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

  We do not attempt to do a cluster-DMFT calculation at this stage.
• For efficiency reasons, we want to solve impurity only for those atoms which are different (and not related by symmetry). While we could figure this out from the structure file, we allow the user to set this here independently of the symmetry in the structure file. This allows one to treat spontaneously broken states (such as magnetically ordered states) in DMFT without reducing the crystal symmetry.

  Enter the correlated problems forming each unique correlated
  problem, separated by spaces (ex: 1,3 2,4 5-8): 1
  

  There is only one impurity in this system anyway.
• Choice for the energy range in the Kohn-Sham basis

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

  This sets the energy range around the Fermi level, where the hybridizations will be taken into account. Usually a value of the order of U is a good starting guess, but one should check the resultant hybridization to see that it is enough. "Projector 1" and "Projector 2" will use all bands in this energy range. "Projector 4" and "5" will use the same set of bands at each k-point, and to faithfully represent the states in this energy range, it will automatically increase the energy range in some k-points.

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

  The CTQMC is an imaginary axis method.

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

  The spin-orbit coupling is ignored.

5 15 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)
1 # number of correlated atoms
1 1 0 # iatom, nL, locrot
 2 7 1 # L, qsplit, cix

#================ # Siginds and crystal-field transformations for correlated orbitals ================
1 5 2 # Number of independent kcix blocks, max dimension, max num-independent-components
1 5 2 # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'eg' 't2g' 
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 1 0 0 0
0 0 2 0 0
0 0 0 2 0
0 0 0 0 2
#---------------- # 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 

1 # number of sigind blocks
5 # dimension of this sigind block
1 0 0 0 0
0 1 0 0 0
0 0 2 0 0
0 0 0 2 0
0 0 0 0 2

solver = 'CTQMC' # impurity solver
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 = 10 # number of iterations of full charge loop (1 = no charge self-consistency).
 # You should probably use 30 or so, but for testing 10 is OK.
ntail = 300 # on imaginary axis, number of points in the tail of the logarithmic mesh
cc = 5e-6 # the charge density precision to stop the LDA+DMFT run
ec = 5e-6 # the energy precision to stop the LDA+DMFT run
recomputeEF = 0 # Recompute EF in dmft2 step. If recomputeEF = 0, EF is fixed. Good for insulators.
DCs = 'exactd' # exact DC with dielectric model 
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" : [9.0 , "# Coulomb repulsion (F0)"],
 "J" : [1.14 , "# Coulomb J"],
 "CoulombF" : ["'Ising'" , "# Density-density form. Can be changed to 'Full' "],
 "beta" : [38.68 , "# Inverse temperature"],
 "svd_lmax" : [25 , "# We will use SVD basis to expand G, with this cutoff"],
 "M" : [5e6 , "# Total number of Monte Carlo steps per core"],
 "mode" : ["SH" , "# We will use Greens function sampling, and Hubbard I tail"],
 "nom" : [100 , "# Number of Matsubara frequency points sampled"],
 "tsample" : [30 , "# How often to record measurements"],
 "GlobalFlip" : [500000 , "# How often to try a global flip"],
 "warmup" : [1e5 , "# Warmup number of QMC steps"],
 "nf0" : [5.0 , "# Nominal valence we expect. Used only for starting DC."],
 }

• solver = 'CTQMC'

  The impurity solver to be used in the calculation. CTQMC stands for Continous Time Quantum Monte Carlo. Only this solver will be used in these tutorials. Alternative solver is 'OCA' (one crossing approximation).

• max_dmft_iterations = 1

  First recall the global computational scheme from the introduction:
  
  There is the dmft loop in red, the DFT loop in grey, and the global loop that connects both loops. This variable sets the maximum number of steps inside the dmft (ref) loop. As we set the variable to 1, a single DMFT step (solution of the impurity problem) will be performed every time the global loop is executed. In this example the impurity solver is much slower than the DFT part, hence we will solve the impurity problem only once inside the charge DFT+Embedded-DMFT loop.

• max_lda_iterations = 100

  After each DMFT step (impurity solver), the code performs up-to 100 DFT steps (the grey loop in figure). Once we are close to the solution, the charge will converge must faster than in 100 steps, hence we will completely converge charge at each self-energy update. Notice that broyden is very efficiently used here (when self-energy is fixed).

• finish = 10

  Maximum total number of DFT+Embedde-DMFT global steps, each of which include single DMFT step and up to 100 DFT steps. Default value is 50.

• ntail = 300

  The self-energy on imaginary axis contains very many Matsubara points (20,000 is quite common), hence we do not want to compute hybridization function for all these points. Fortunately, high frequency tail on imaginary axis contains very little information, hence it can be very well described by a logarithmic mesh of Matsubara points. In computing hybridization function and charge density, we compute first "nom" points exactly, and the rest are approximated with ntail=300 logarithmically distributed Matsubara points. Internally, we either use spline or second order interpolation, to get functions on all Matsubara points.

• cc = 5e-6 # the charge density precision to stop the LDA+DMFT run
  ec = 5e-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 = 0

  The chemical potential will not be re-calculated at each step. For metals, we should always re-calculate the chemical potential. For Mott insulators, the continuous update of the chemical potential does not work well, because the position of the Fermi level is not well defined and its position will jump inside the gap (due to numerics).
  
  One way is to first switch on the search for the chemical potential (recomputeEF=1), and then monitor the progress of the calculation. Once it starts jumping inside the gap, one can fix it roughly into the middle of the gap, to allow nice convergence of the DMFT loop. There are several ways to estimate what is the "middle of the gap". Maybe the easiest is to look at the impurity Green's function (imp.0/Gf.out.?.1). When the gap opens, the imaginary part of G (for all orbitals) will go to zero at zero frequency. The real part of G will be small when G appears roughly particle-hole symmetric at low energy, which means that chemical potential is roughly in the middle of the gap. When we see that for some iterations Im(G(0))=0 and Re(G(0)) small, we should look into info.iterate what was the chemical potential at that step, and set this value into EF.dat.
  
  Far more simple way is to just fix the chemical potential to the LDA value (fix it at the very beginning of the run), and than wait to see if the gap in the Green's function develops after a few iterations. If it does, the chemical potential is most likely in the gap. To make sure that this is not a "fake" gap, we need to scroll to the end of dmft2_info.out, and check the line which says someting like:

  Ratio to renormalize= 1.000120921355 rho-rho_expected= -0.0022972280

  The same information is also printed in MnO.scf file with the key-word :DRHO, hence

  grep ':DRHO' MnO.scf

  will print the entire history of charge mismatch. This number shows that the charge density is very close to what we expected. Currently it has 2.3e-3 electron missing, which is definitely due to finite precision in computing the electronic charge. Notice that if we hit a "fake" gap, the charge density will be wrong for an integer number, i.e., "rho-rho_expected" will be close to +1.0 or -1.. Note also that at the first few iterations the density might not be so precise and "rho-rho_expected" might be finite. But once the Green's function stops changing, the density should be correct.

• DCs = 'exactd' # exact double-counting with dielectric constant 

  The double counting scheme (DCs): In this calculation we are using the exact double-counting assuming dielectric screening.
  

  There are several possible choices for the double-counting.

  • The most stable is the nominal choice
    DCs='nominal'
    in which case we need to specify the expected nominal occupancy of the ionic state. This is done through parameter nf0 in iparams0 (see below). We then evaluate
    
    

    As explained in PRL 115, 196403, this choice is quite close to the exact DC, which we explain below.

  The other choices are called exact, as they are derivable from the common part of the DFT and DMFT functional (PRL 115, 196403). However, they differ in their choice for the real space form of the Coulomb repulsion.
  • For many metallic systems the Yukawa form of the screening is quite good, i.e.,
    DCs='exacty'
    which means that in real space, the DMFT Coulomb repulsion is assumed to have the form
    and the screening length is determined such that
  • In some insulators, the dielectric screening is better, which can be used by selecting the option
    DCs='exactd'
    which assumes
    
  • Most often the best choice is somewhere between the two described options (yukawa and dielectric), and we implemented a combination of the two, which we call
    DCs='exact'
    and assumes the form
    Of course, the equation still needs to be satisfied.
  • We also allow the so-called fully-localized-limit
    DCs='FLL'
    double-counting, but we do not recommend it.
  Note that in MnO the results of exactd, nominal and exact double-countings are similar, for example converged nd=5.03 with exactd and nd=4.97 with exact. However, when using exact, the chemical potential needs to be moved to higher energy. As explained above, this can be achieved by first selecting recomputeEF=1, and after a few steps we switch it back to recomputeEF=0 to better converge the equations.

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

  This is relevant only for sampling in Matsubara frequency (svd_lmax=0), but not when svd_lmax>0, in which case we expand the Green's function in the SVD basis. It can be omitted for svd_lmax>0. The first "nom" Matsubara points of the self-energy are sampled by CTQMC, while the rest of it is replaced by a tail, which however is also computed by CTQMC. Due to Monte Carlo noise, we sometimes see substantial scattering of data at intermediate frequency. Sharp spikes in self-energy can lead to numerical instabilities and we want to avoid them. If we use Matsubara sampling, it is advisable to use finite broadening, for example wbroad=0.02 and kbroad=0.05. so that the self-energy is broadened with gaussian broadening, which has width = wbroad + kbroad*omega. Therefore the low energy self-energy is not affected, while the high frequency is smoothened out. To see the effect of this broadening, one can compare the function imp.0/Sig.outb (data) and imp.0/Sig.out (broadened).

• "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" : [9.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 many itinerant states) is explicitly taken into account in the such DMFT calculation. We will calculate the value of U in later tutorial.

• "J" : [1.14]

  The Hund's coupling. It can be estimated by running

  RCoulombU.py -U 9.

  which will give something like

  lambda= 0.975327482315 epsilon= 1.0
  Fk= [8.999999999999998, 9.807040121910346, 6.655370991435106]
  U,J= [8.999999999999998, 1.1383171570074508, 1.2359974698379483]

  This uses Yukawa screening to estimate the relation between U and J. When U is 9eV, J is around 1.14eV. Notice that unscreened J is 1.3eV for our very localized orbital (you can get it by running RCoulombU.py with no arguments). Internally, we use Slater F2 and F4, but here we rather give more commonly used parameter "J", which gives 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].

• "CoulombF" : ["'Ising'"]

  Here we use the density-density form of the Coulomb repulsion, which speeds up the calculation considerably. However, one should be careful with use of this form of Coulomb repulsion, as it can sometimes lead to considerably different results than the rotationally invariant form. To use rotationally invariant Slater form, set

  "CoulombF" : ["'Full'"]

• "beta" : [38.68 ]

  Inverse temperature in eV. beta=38.68 means that temperature is 3868 meV, which is about 300 Kelvins.

• "svd_lmax" : [25 ]