ABINIT tutorial, lesson on calculation of U and J with constrained Random Phase Approximation (cRPA) :

the case of SrVO3


This lesson aims at showing how to perform a calculation of U and J in Abinit using cRPA. This method is well adapted in particular to determine U and J as they can be used in DFT+DMFT. The paper describing the implementation is [Amadon2014].

It might be useful that you already know how to do PAW calculations using ABINIT but it is not mandatory (you can follow the two lessons on PAW in ABINIT (PAW1, PAW2)). The DFT+U tutorial in ABINIT (DFT+U) might be useful to know some basic variables about correlated orbitals.

The first GW tutorial in ABINIT (GW) is useful to learn how to compute the screening, and how to converge the relevant parameters (energy cutoffs and number of bands for the polarizability).

This lesson should take two hours to complete (you should have access to more than 8 processors).

Copyright (C) 2008-2017 ABINIT group (BAmadon)
This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

Content of the cRPA lesson


1. The cRPA method to compute effective interaction: summary and key parameters

The cRPA method aims at computing the effective interactions among correlated electrons. Generally, these highly correlated materials contain rare-earth metals or transition metals, which have partially filled d or f bands and thus localized electrons. cRPA relies on the fact that screening processes can be decomposed in two steps: Firstly, the bare Coulomb interaction is screened by non correlated electrons to produce the effective interaction Wr. Secondly, correlated electrons screened this interaction to produce the fully screening interaction W. (see [Aryasetiawan2004]). However, the second screening process is taken into account when one uses a method which describes accurately the interaction among correlated electrons (such as Quantum Monte Carlo within the DFT+DMFT method). So, to avoid a double counting of screening by correlated electrons, the DFT+DMFT methods needs as an input the effective interaction Wr. The goal of this tutorial is to present the implementation of this method using Projected Local Orbitals Wannier orbitals in ABINIT (The implementation of cRPA in ABINIT is described in [Amadon2014] and projected local orbitals Wannier functions are presented in [Amadon2008] ). The discussion about the localization of Wannier orbitals has some similarities with the beginning on the DMFT tutorial (see here and there)

Several parameters (both physical and technical) are important for the cRPA calculation:


2. Electronic Structure of SrVO3 in LDA

You might create a subdirectory of the ~abinit/tests/tutoparal directory, and use it for the tutorial. In what follows, the names of files will be mentioned as if you were in this subdirectory

Copy the files ../Input/tucrpa_1.in and ../Input/tucrpa_1.files in your Work directory, and run ABINIT: (as usual, the actual "abinit" command is something like ../../../../src/98_main/abinit):

abinit < tucrpa_1.files > log_1

This run should take some time. It is recommended that you use at least 10 processors (and 32 should be fast). It calculates the LDA ground state of SrVO3 and compute the band structure in a second step. The variable pawfatbnd allows to create files with "fatbands" (see description of the variable in the list of variables): the width of the line along each k-point path and for each band is proportional to the contribution of a given atomic orbital on this particular Kohn Sham Wavefunction. A low cutoff and a small number of k-points are used in order to speed up the calculation. During this time you can take a look at the input file. There are two datasets. The first one is a ground state calculations with nnsclo=3 and nline=3 in order to have well diagonalized eigenfunctions even for empty states. In practice, you have however to check that the residue of wavefunctions is small at the end of the calculation. In this calculation, we find 1.E-06, which is large (1.E-10 would be better, so nline and nnsclo should be increased, but it would take more time). When the calculation is finished, you can plot the fatbands for Vanadium and l=2 with

xmgrace tucrpa_O_DS2_FATBANDS_at0001_V_is1_l0002
The band structure is given in eV.

FatbandV

and the fatbands for one oxygen atom and l=1 with

 
xmgrace tucrpa_O_DS2_FATBANDS_at0003_O_is1_l0001.

FatbandV

In these plots, you recover the band structure of SrVO3 (see for comparison the band structure of Fig.3 of [Amadon2008]), and the main character of the bands. Bands 21 to 25 are mainly d and bands 12 to 20 are mainly oxygen p. However, we clearly see an important hybridization. The Fermi level (at 0 eV) is in the middle of bands 21-23.

One can easily check thats bands 21-23 are mainly d-t2g and bands 24-25 are mainly eg: just use pawfatbnd=2 in tucrpa_1.in and relaunch the calculations. Then the file tucrpa_O_DS2_FATBANDS_at0001_V_is1_l2_m-2, tucrpa_O_DS2_FATBANDS_at0001_V_is1_l2_m-1 and tucrpa_O_DS2_FATBANDS_at0001_V_is1_l2_m1 give you respectively the xy,yz and xz fatbands (ie d-t2g) and tucrpa_O_DS2_FATBANDS_at0001_V_is1_l2_m+0 and tucrpa_O_DS2_FATBANDS_at0001_V_is1_l2_m+2 give the z2 and z2-y2 fatbands (ie eg).

So in conclusion of this study, the Kohn Sham bands which are mainly t2g are the bands 21, 22 and 23.

Of course, it could have been anticipated from classical crystal field theory: the vanadium is in the center of an octahedron of oxygen atoms, so d orbitals are splitted in t2g and eg. As t2g orbitals are not directed toward oxygen atoms, t2g-like bands are lower in energy and filled with one electron, whereas eg-like bands are higher and empty.

In the next section, we will thus use the d-like bands to built Wannier functions and compute effective interactions for these orbitals.

3. Definition of screening and orbitals, input and log file

3.1. The constrained polarization calculation.

As discussed briefly in Appendix A of [Amadon2014] as well as in section III B [Vaugier2012] and Section II B of [Sakuma2013], one can use different schemes for the cRPA calculations. Let us discuss these different models namely the (d-d), (dp-dp), (d-dp(a)), and (d-dp(b)) models. In the notation (A-B), A and B refers respectively to some bands of A-like and B-like character. Moreover, A refers to the definition of screening and B refers to the definition of correlated orbitals. To clarify this definition, we give below some examples: Which way of computing interactions is the most relevant depends on the way the interactions will be used. Some aspects of it are discussed in Ref. [Vaugier2012]. Also, for Mott insulators, and if self-consistency over interactions is carried out, the choice of models is discussed in [Amadon2014].

3.2. The input file for cRPA calculation: correlated orbitals, Wannier functions

In this section, we will present the input variables and discuss how to extract useful information in the log file in the case of the d-d model. The input file for a typical cRPA calculation (tucrpa_1.in) contains four datasets (as usual GW calculations, see the GW tutorial): the first one is a well converged LDA calculation, the second is non self-consistent calculation to compute accurately full and empty states, the third computes the constrained non interacting polarizability, and the fourth computes effective interaction parameters U and J. We discuss these four datasets in the next four subsections.

Copy the files ../Input/tucrpa_2.in and ../Input/tucrpa_2.files in your Work directory. The input file tucrpa_2.in contains standard data to perform a LDA calculation on SrVO3. We focus in the next subsections on some peculiar input variables related to the fact that we perform a cRPA calculation. Before reading the following section, launch the abinit calculation:

abinit < tucrpa_2.files > log_2
3.2.1. The first DATASET: A converged LDA calculation
The first dataset is a simple calculation of the density using LDA using 50 bands. We do not use symmetry in this calculation, so nsym=1. If symmetry is used (nsym=0), then the calculation of U and J will be valid, but the full interaction matrix described in section 3.2.4 will not be correct.
3.2.2. The second DATASET: A converged LDA calculation and definition of Wannier functions

Before presenting the input variables for this dataset, we discuss two important physical parameters relevant to this dataset.

iscf2          -2    # Perform a non self-consistent calculation              
nbandkss2      -1    # Number of bands in KSS file (-1 means the maximum possible) 
kssform         3    # Format of the Wavefunction file (should be 3) 
nbdbuf2         4    # The last four bands will not be perfectly diagonalized 
tolwfr2   1.0d-18    # The criterion to stop diagonalization 

# == Compute Wannier functions 
usedmft2        1    # Mandatory to enable the calculation of Wannier functions as in DMFT. 
dmftbandi2     21    # Precise the definition of Wannier functions (also used for DMFT calculations) 
dmftbandf2     25    # Precise the definition of Wannier functions (also used for DMFT calculations) 
3.2.2. The third DATASET: Compute the constrained polarizability and dielectric function
The third DATASET drives the computation of the constrained polarizability and the dielectric function. Again, we reproduce below the input variables peculiar to this dataset.

We add some comments here on three most important topics

 optdriver3     3     # Keyword to launch the calculation of screening 
 gwcalctyp3     2     # The screening will be used later with gwcalctyp 2 in the next dataset 
 getwfk3       -1     # Obtain WFK file from previous dataset
 ecuteps3     5.0     # Cut-off energy of the planewave set to represent the dielectric matrix.
                              #It is important to converge effective interactions as a function of this parameter.

# -- Frequencies for dielectric matrix
 nfreqre3       1     # Number of  real frequencies 
 freqremax3    10 eV  # Maximal value of frequencies 
 freqremin3     0 eV  # Minimal value of frequencies 
 nfreqim3       0     # Number of  imaginary frequencies 

# -- Ucrpa: screening
 ucrpa_bands3  21 25  # Define the bands corresponding to the d contribution 

# -- Parallelism
 gwpara3        1
3.2.3. The fourth DATASET: Compute the effective interactions
The fourth DATASET drives the computation of the effective interactions. It is similar computationally to the computation of exchange in GW calculations. Again, we reproduce below the input variables peculiar to this dataset.

We add some comments on convergence properties

 optdriver4  4      # Self-Energy calculation
 gwcalctyp4  2      # activate HF or ucrpa
 getwfk4     1      # Obtain WFK file from dataset 1
 getscr4     2      # Obtain SCR file from previous dataset
 ecutsigx4  30.0    # Dimension of the G sum in Sigma_x.

# -- Frequencies for effective interactions
 nfreqsp4    1
 freqspmax4 10 eV
 freqspmin4  0 eV
 nkptgw4     0      # number of k-point where to calculate the GW correction: all BZ
 mqgrid4   300      # Reduced but fine at least for SrVO3
 mqgriddg4 300      # Reduced but fine at least for SrVO3

# -- Parallelism
 gwpara4 2          # do not change if nsppol=2
3.2.4 The cRPA calculation: the log file (for the d-d model)

We are now going to browse quickly the log file for this calculation.

4. Convergence studies

We give here the results of some convergence studies, than can be made by the readers. Some are computationally expensive. It is recommanded to use at least 32 processors. Input files are provided in tucrpa_3.in and tucrpa_3.files for the first case.

4.1 ecuteps

The number of plane waves used in the calculation of the cRPA polarizability is determined by ecuteps (see Notes on RPA calculations). The convergence can be studied simply by increasing the values of ecuteps and gives:
   ecuteps (Ha)   U (eV)     J (eV)
      3            3.22       0.43
      5            3.01       0.42
      7            2.99       0.40
      9            2.99       0.40
So, for nband=30, ecutsigx=30 and a 4x4x4 k-point grid, the effective interactions are converged with a precision of 0.02 eV for ecuteps=5.

4.2 nband

The RPA polarizability depends also of the number of Kohn-Sham bands as also discussed in the Notes on RPA calculations.
   nband          U (eV)     J (eV)
    30             3.01       0.42
    50             2.75       0.41
    70             2.71       0.41
    90             2.70       0.40
So, for ecuteps=5, ecutsigx=30 and a 4x4x4 k-point grid, the effective interactions are converged with a precision of 0.05 eV for nband=50.

4.3 ecutsigx

ecutsigx is the input variable determining the number of plane waves used in the calculation of the exchange part of the self-energy. The same variable is used here to determine the number of plane waves used in the calculation of the effective interaction. The study of the convergence gives:
     ecutsigx (Ha)            U_bare (eV)    J_bare (eV)      U (eV)       J (eV)
         30                    15.38          0.45            3.22         0.43
         50                    15.38          0.46            3.22         0.44
         70                    15.38          0.47            3.22         0.45
For ecuteps=3, nband=30 and a 4 4 4 k-point grid, effective interactions are converged with a precision of 0.02 eV for ecutsigx=30 Ha.

4.4 k-point

The dependence of effective interactions with the k-point grid is also important to check.
   kpoint grid     U_bare (eV)    J_bare (eV)     U (eV)       J (eV)
      4 4 4          15.38          0.45           3.22         0.43
      5 5 5          15.38          0.44           3.14         0.43
      6 6 6          15.36          0.44           3.19         0.43
      7 7 7          15.35          0.44           3.19         0.43
For ecuteps=3, nband=30 and ecutsigx=30, effective interactions are converged to 0.01 eV for a k-point grid of 6 6 6. The converged parameters are thus ecuteps=5,ecutsigx=30,ngkpt=6 6 6, nband=50 with an expected precision of 0.1 eV. We will use instead ecuteps=5,ecutsigx=30,ngkpt=4 4 4, nband=50 to lower the computational cost. In this case, we find values of U and J of 2.75 eV and 0.41 eV and U is probably underestimated by 0.1 eV.

5. Effective interactions for different models

In this section, we compute, using the converged values of parameters, the value of effective interactions for the models discussed in section and 3.1. The table below gives for each model, the values of dmftbandi, dmftbandf and ucrpa_bands that must be used, and it sums up the value of bare and effective interactions.

      model      ucrpa dmftbandi/dmftbandf ucrpa_bands |  Ubare    Ubare diag   Jbare |   U      Udiag     J 
      d - d        1        21/25            21 25     |  15.4      16.3      0.45  |  2.8     3.5     0.41
     t2g-t2g       1        21/25            21 25     |  15.3      16.0      0.48  |  2.8     3.4     0.44
      dp-dp        1        12/25            12 25     |  19.4      20.6      0.67  | 10.1    11.2     0.64
      d -dp (a)    1        12/25            21 25     |  19.4      20.6      0.67  |  3.4     4.4     0.62
      d -dp (b)    2        12/25            12 25     |  19.4      20.6      0.67  |  1.6     2.6     0.60
Even if our calculation is not perfectly converged with respect to k-points, the results obtained for the d-d, dp-dp, and d-dp (a) models are in agreement (within 0.1 or 0.2 eV) with results obtained in Table V of [Amadon2014] and references cited in this table.

To obtain the results with only the t2g orbitals, one must use a specific input file, which is tucrpa_4.in, which uses specific keywords, peculiar to this case (compare it with tucrpa_2.in). We now briefly comment the physics of the results.

Finally, fully screened interactions (not discussed here) could also be computed for each definition correlated orbitals by using the default values of ucrpa_bands .

6. Frequency dependent interactions

In this section, we are going to compute frequency dependent effective interactions. Just change the following input variables in order to compute effective interaction between 0 and 30 eV. We use nsym=0 in order to decrease the computational cost.
 # -- Frequencies for effective interactions
 nsym        0
 nfreqsp3   60
 freqspmax3 30 eV
 freqspmin3  0 eV
 nfreqsp4   60
 freqspmax4 30 eV
 freqspmin4  0 eV
An example of input file can be found in tucrpa_5.in. Note that we have decreased some parameters to speed-up the calculations. Importantly, however, we have increased the number of Kohn Sham bands, because calculation of screening at high frequency involves high energy transitions which requires high energy states (as well as semicore states). If the calculation is too time consuming, you can reduce the number of frequencies. The following figure has been plotted with 300 frequencies, but using 30 or 60 frequencies is sufficient to see the main tendencies. Extract the value of U for the 60 frequencies using:
grep "U(omega)" -A 60 tucrpa_5.out > Ufreq.dat
Remove the first line, and plot the data:

dynamical

We recover in the picture the main results found in the Fig. 3 of [Aryasetiawan2006] . Indeed, the frequency dependent effective interactions exhibits a plasmon excitation.

7. Conclusion

This simple tutorial showed how to compute effective interactions for SrVO3. It emphasized the importance of the definition of correlated orbitals and the definition of constrained polarizability.

References

[Aryasetiawan2004] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein Phys. Rev. B 70, 195104 (2004)
[Amadon2008] B. Amadon, F. Lechermann, A. Georges, F. Jollet, T. O. Wehling, and A. I. Lichtenstein Phys. Rev. B 77, 205112 (2008)
[Amadon2008a] B. Amadon, F. Jollet, and M. Torrent Phys. Rev. B 77, 155104 (2008)
[Vaugier2012] Loig Vaugier, Hong Jiang, and Silke Biermann Phys. Rev. B 86, 165105 (2012)
[Sakuma2013] R. Sakuma and F. Aryasetiawan Phys. Rev. B 87, 165118 (2013)
[Amadon2014] Bernard Amadon, Thomas Applencourt, and Fabien Bruneval Phys. Rev. B 89, 125110 (2014).
[Geneste2016] G. Geneste, B. Amadon, M. Torrent and G. Dezanneau (unpublished)
[Aryasetiawan2006] F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schönberger, Phys. Rev. B 74, 125106 (2006).