Calculation of U and J using the linear response approach¶
How to determine U for DFT+U in ABINIT? Cococcioni’s linear response approach.¶
This tutorial aims to show how you can determine U for further DFT+U calculations consistently and in a fast an easy way. You will learn how to prepare the input files for the determination and how to use the main parameters implemented for this aim. It is supposed that you already know how to run ABINIT with PAW (tutorial PAW1). Obviously, you should also read the tutorial DFT+U, and likely the tutorial PAW2, to generate PAW atomic data.
This tutorial should take about ½ hour.
Note
Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.
In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:
export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/ # Do not change this line: path to tests dir
export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/ # Do not change this line: path to pseudos dir
Examples in this tutorial use these shell variables: copy and paste
the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively,
source the set_abienv.sh script located in the ~abinit directory:
source ~abinit/set_abienv.sh
The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.
To execute the tutorials, create a working directory (Work*) and
copy there the input files of the lesson.
Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.
Summary of linear response method to determine U¶
The linear response method has been introduced by several authors [Cococcioni2002], [Cococcioni2005],[Dederichs1984],[Hybertsen1989], [Anisimov1991],[Pickett1998]. It is based on the fact that U corresponds to the energy to localize an additional electron on the same site: U=E[n+1]+E[n-1]-2E[n] [Hybertsen1989]. This can be reformulated as the response to an infinitesimal change of of occupation of the orbital by the electrons dn. Then U is the second derivative of the energy with respect to the occupation U=\frac{\delta^2 E}{\delta^2 n}. The first method fixed the occupation by cutting the hopping terms of localized orbitals. Later propositions constrained the occupation through Lagrange multipliers [Dederichs1984],[Anisimov1991]. The Lagrange multiplier \alpha corresponds to a local potential that has to be applied to increase or decrease the occupation by ±1 electron. Note that the occupation need not to vary by 1 electron, but the occupation shift can be infinitesimal.
It is recommended to read the following papers to understand the basic concepts of the linear response calculations to calculate U:
[1] “A LDA+U study of selected iron compounds “, M. Cococcioni, Ph.D. thesis, International School for Advanced Studies (SISSA), Trieste (2002) [Cococcioni2002]
[2] “Linear response approach to the calculation of the effective interaction parameters in the LDA + U method”, M. Cococcioni and S. de Gironcoli, Physical Review B 71, 035105 (2005) [Cococcioni2005]
Some further reading:
[3] “Ground States of Constrained Systems: Application to Cerium Impurities”, P. H. Dederichs, S. Blugel, R. Zeller, and H. Akai, Phys. Rev. Lett. 53, 2512 (1984) [Dederichs1984]
[4] “Calculation of Coulomb-interaction parameters for La2CuO4 using a constrained-density-functional approach”, M. S. Hybertsen, M. Schluter, and N. E. Christensen, Phys. Rev. B 39, 9028 (1989) [Hybertsen1989]
[5] “Density-functional calculation of effective Coulomb interactions in metals”, V. I. Anisimov and O. Gunnarsson, Phys. Rev. B42, 7570 (1991) [Anisimov1991]
[6] “Reformulation of the LDA+U method for a local-orbital basis”, W. E. Pickett, S. C. Erwin, and E. C. Ethridge, Phys. Rev. B58, 1201 (1998) [Pickett1998]
The implementation of the determination of U in ABINIT is briefly discussed in [Gonze2016].
How to determine U in ABINIT¶
Before continuing, you may consider to work in a different subdirectory as for the other tutorials. Why not Work_ucalc_lr?
Important
In what follows, the name of files are mentioned as if you were in this subdirectory. All the input files can be found in the $ABI_TESTS/tutorial/Input directory You can compare your results with reference output files located in $ABI_TESTS/tutorial/Refs directory (for the present tutorial they are named tucalc_lr*.abo).
The input file tucalc_lr_1.abi is an example of a file to prepare a wave function for further processing. The corresponding output file is ../Refs/tucalc_lr_1.abo).
Copy the files tucalc_lr_1.abi in your work directory, and run ABINIT:
cd $ABI_TESTS/tutorial/Input
mkdir Work_ucalc_lr
cd Work_calc_lr
cp ../tucalc_lr_1.abi .
abinit tucalc_lr_1.abi > log 2> err &
In the meantime, you can read the input file and see that this is a usual DFT+U calculation, with U=0.
################################################################# # Automatic test for ABINIT: # # Preliminary step for test v5#39 (macro_uj) # # and v5#40 (testirdden) # # # # Fe bcc 2 atomic supercell - ferromag.- PAW DJA 2010 & MT 2009 # ################################################################# #Unit cell acell 3*5.42 chkprim 0 # 0: do not check if uc primitive rprim 1.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 1.00 #Spin polarization nsppol 2 #1 unpolarized / 2 polarized spinat 0 0 2.843062 0 0 2.843062 #Definition of the atom types ntypat 1 znucl 26 #Definition of the atoms natom 2 typat 2*1 # atomic types xred 0.0 0.0 0.0 0.5 0.5 0.5 ecut 8 # Energy cutoff pawecutdg 20 # pawecutdg > 2*ecut nband 25 # Fe_2 minband=17 #Definition of the SCF procedure nstep 15 # max number SCF cycles tolvrs 10d-12 #Definition of the k-point grid kptopt 1 # 1: automatic generation of k points ngkpt 3 3 3 # n x n x n nshiftk 1 shiftk 0.5 0.5 0.5 #Smearing occopt 4 tsmear 0.05 eV #DFT+U usepawu 1 # 1 at lim dble cnt / 2 rnd m fld dle cnt lpawu 2 # ang moments corrrected #Save disk space & Miscelaneous prteig 0 prtden 1 # This is the default value #PAW atomic data outdata_prefix = "tucalc_lr_1.o" pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_paw_pbe_standard/Fe.xml" ixc 11 # Imposing use of internal PBE for portability purposes ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% test_chain = tucalc_lr_1.abi, tucalc_lr_2.abi, tucalc_lr_3.abi #%% [files] #%% files_to_test = #%% tucalc_lr_1.abo, tolnlines= 2, tolabs= 7.000e-09, tolrel= 1.7e-09 #%% output_file = "tucalc_lr_1.abo" #%% [paral_info] #%% max_nprocs = 1 #%% [extra_info] #%% authors = M. Torrent #%% keywords = PAW, DFTU #%% description = #%% Fe bcc 2 atomic supercell - ferromag.- PAW DJA 2010 & MT 2009 #%% Preliminary step for test v5#39 (macro_uj) and v5#40 (testirdden) #%%<END TEST_INFO>
This setting allows us to read the occupations of the Fe 3d orbitals (lpawu 2). The cell contains 2 atoms. This is the minimum to get reasonable response matrices. We converge the electronic structure to a high accuracy (tolvrs 10d-12), which usually allows one to determine occupations with a precision of 10d-10. The ecut is chosen very low, in order to accelerate calculations. We do not suppress the writing of the WFK file, because this is the input for the calculations of U.
Once this calculation has finished, run the second one: Copy the file tucalc_lr_2.abi in your work directory, and run ABINIT:
abinit tucalc_lr_2.abi > tucalc_lr_2.log
################################################################### ## Automatic test for ABINIT: ## ## determine U from change of occupation on atoms upon potential ## ## shift on atom 1 ## ## Fe bcc structure - ferromagnetic PAW DJA 2010 ## ################################################################### # 2 atomic supercell acell 3*5.42 chkprim 0 # 0: do not check if uc primitive rprim 1.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 1.00 nsppol 2 #Definition of the atom types ntypat 1 znucl 26 #Definition of the atoms natom 2 typat 2*1 # atomic types xred 0.0 0.0 0.0 0.5 0.5 0.5 ecut 8 # Energy cutoff pawecutdg 40 # pawecutdg > 2*ecut nband 25 # Fe_2 minband=17 #Definition of the k-point grid kptopt 1 # 1: automatic generation of k points ngkpt 3 3 3 # n x n x n nshiftk 1 shiftk 0.5 0.5 0.5 #Smearing occopt 4 tsmear 0.05 eV #DFT+U usepawu 1 # 1 at lim dble cnt / 2 rnd m fld dle cnt lpawu 2 # ang moments corrrected nsym 48 # nsym&symrel: break cubic symmetry of crystal: allow # individual ionicoccupations symrel 1 0 0 0 1 0 0 0 1 -1 0 0 0 -1 0 0 0 -1 -1 0 0 0 1 0 0 0 -1 1 0 0 0 -1 0 0 0 1 -1 0 0 0 -1 0 0 0 1 1 0 0 0 1 0 0 0 -1 1 0 0 0 -1 0 0 0 -1 -1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 0 1 0 -1 0 -1 0 0 0 0 -1 0 -1 0 1 0 0 0 0 -1 0 1 0 -1 0 0 0 0 1 0 -1 0 -1 0 0 0 0 1 0 1 0 1 0 0 0 0 -1 0 1 0 -1 0 0 0 0 -1 0 -1 0 1 0 0 0 0 1 0 0 1 1 0 0 0 1 0 0 0 -1 -1 0 0 0 -1 0 0 0 -1 1 0 0 0 -1 0 0 0 1 -1 0 0 0 1 0 0 0 -1 -1 0 0 0 1 0 0 0 1 1 0 0 0 -1 0 0 0 1 -1 0 0 0 -1 0 0 0 -1 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 -1 0 0 0 0 -1 0 -1 0 -1 0 0 0 0 1 0 -1 0 1 0 0 0 0 -1 0 1 0 -1 0 0 0 0 -1 0 1 0 1 0 0 0 0 1 0 -1 0 1 0 0 0 0 -1 0 -1 0 -1 0 0 0 0 1 0 1 0 0 1 0 0 0 1 1 0 0 0 -1 0 0 0 -1 -1 0 0 0 -1 0 0 0 1 -1 0 0 0 1 0 0 0 -1 1 0 0 0 -1 0 0 0 -1 1 0 0 0 1 0 0 0 1 -1 0 0 0 1 0 0 0 -1 -1 0 0 0 -1 0 0 0 1 1 0 0 0 0 1 0 1 0 1 0 0 0 0 -1 0 -1 0 -1 0 0 0 0 -1 0 1 0 -1 0 0 0 0 1 0 -1 0 1 0 0 0 0 -1 0 -1 0 1 0 0 0 0 1 0 1 0 -1 0 0 0 0 1 0 -1 0 -1 0 0 0 0 -1 0 1 0 1 0 0 pawujat 1 # default, the atom on which U is determined pawujv 0.1 eV # default, size of the potential shift macro_uj 1 # activate determination of U pawujrad 2.66866 # optional, radius ASA-sphere to which U should be extrapolated #Only to accelerate test irdwfk 1 # default for macro_uj = 1 # nline 2 # nnsclo 2 tolvrs 10d-9 # default for macro_uj = 1 #Save disk space prteig 0 prtwf 0 prtden 0 #PAW atomic data indata_prefix = "tucalc_lr_1.o" pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_paw_pbe_standard/Fe.xml" ixc 11 # Imposing use of internal PBE for portability purposes ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% test_chain = tucalc_lr_1.abi, tucalc_lr_2.abi, tucalc_lr_3.abi #%% [files] #%% files_to_test = #%% tucalc_lr_2.abo, tolnlines= 4, tolabs= 8.000e-03, tolrel= 6.500e-01, fld_options = -easy #%% output_file = "tucalc_lr_2.abo" #%% [paral_info] #%% max_nprocs = 1 #%% [extra_info] #%% authors = D.J. Adams #%% keywords = PAW, DFTU #%% description = #%% Fe bcc structure - ferromagnetic PAW #%% determine U from change of occupation on atoms upon potential #%% shift on atom 1 #%%<END TEST_INFO>
As you can see from the tucalc_lr_2.abi file, this run uses the tucalc_lr_1o_WFK as an input (as indata_prefix = “tucalc_lr_1.o”). In the tucalc_lr_2.abi all the symmetry relations are specified explicitly. In the tucalc_lr_2.log you can verify that none of the symmetries connects atoms 1 with atom 2:
symatm: atom number 1 is reached starting at atom
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
symatm: atom number 2 is reached starting at atom
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
This is important. Otherwise the occupation numbers have no freedom to evolve separately on the atoms surrounding the atom on which you apply the perturbation.
You can generate these symmetries, in a separate run, where you specify the atom where the perturbation is done as a different species. From the output you read the number of symmetries (nsym), the symmetry operations (symrel) and the non-symmorphic vectors (tnons). This is already done here and inserted in the tucalc_lr_2.abi file. Note that you can alternatively disable all symmetries with nsym = 1, or break specific symmetries by displacing the impurity atom in the preliminary run. However, for the determination of U, the positions should be the ideal positions and only the symmetry should be reduced.
For the rest, usually it is enough to set macro_uj 1 to run the calculation of U. Note also, that the irdwfk 1 and the tolvrs 1d-8 need not be set explicitly, because they are the defaults with macro_uj 1.
Once the calculation tucalc_lr_2 is converged, you can have look at the output. You can see, that the atomic shift (atvshift) is automatically set:
atvshift 0.00367 0.00367 0.00367 0.00367 0.00367
0.00367 0.00367 0.00367 0.00367 0.00367
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
This means, that all the 10 3d spin-spin orbitals on the first Fe atom where shifted by 0.1 eV (=0.00367 Ha). On the second atom no shift was applied. Self-consistency was reached twice: Once for a positive shift, once for the negative shift:
grep SCF tucalc_lr_2.abo
The lines starting with URES
URES ii nat r_max U(J)[eV] U_ASA[eV] U_inf[eV]
URES 1 2 4.69390 3.86321 3.10851 2.71778
URES 2 16 9.38770 7.28015 5.85793 5.12160
URES 3 54 14.08160 7.60761 6.12142 5.35197
URES 4 128 18.77540 7.67652 6.17686 5.40045
URES 5 250 23.46930 7.69879 6.19478 5.41611
contain U for different supercells. These values of U are computed using the extrapolation procedure proposed in [Cococcioni2005]. In this work, it is shown that using a two atom supercell for the DFT calculation and an extrapolation procedure allows a estimation of the value of U. More precise values can be obtained using a larger supercell for the DFT calculation and an extrapolation. For simplicity, in this tutorial, we restrict to the two atoms supercell for the DFT calculation.
The column “nat” indicates how many atoms were involved in the extrapolated supercell, r_max indicates the maximal distance of the impurity atoms in that supercell. The column U indicates the actual U you calculated and should use in your further calculations. U_ASA is an estimate of U for more extended projectors and U_inf is the estimate for a projector extended even further.
Although it is enough to set macro_uj 1, you can further tune your runs. As a standard, the potential shift to the 1st atom treated in DFT+U, with a potential shift of 0.1 eV. If you wish to determine U on the second atom you put pawujat 2. To change the size of the potential shift use e.g. pawujv 0.05 eV. Our tests show that 0.1 eV is the optimal value, but the linear response is linear in a wide range (1-0.001 eV).
The ujdet utility¶
In general the calculation of U with abinit as described above is sufficient. For some post-treatment that goes beyond the standard applications, a separate executable ujdet was created. The output of abinit is formatted so that you can easily “cut” the part with the ujdet input variables: you can generate the standard input file for the ujdet utility by typing:
sed -n "/MARK/,/MARK/p" tucalc_lr_2.abo > ujdet.in
Note that the input for the ujdet utility is always called ujdet.in
It contains the potential shifts applied vsh (there are 4 shifts: vsh1, vsh3 for non-selfconsistent calculations that allows to extract the contribution to U originating from a non-interacting electron gas, and vsh2, vsh4 for positive and negative potential shift). The same applies for the occupations occ[1-4].
We now calculate U for an even larger supercell: Uncomment the line scdim in ujdet.in and add
scdim 6 6 6
to specify a 6 6 6 supercell or
scdim 700 0 0
to specify the maximum total number of atoms in the supercell. Then, run ujdet (the executable is in the same directory as the abinit executable):
rm ujdet.[ol]* ; ujdet > ujdet.log
grep URES ujdet.out
URES ii nat r_max U(J)[eV] U_ASA[eV] U_inf[eV]
URES 1 2 4.69390 3.86321 3.10851 2.71778
URES 2 16 9.38770 7.28015 5.85793 5.12160
URES 3 54 14.08160 7.60761 6.12142 5.35197
URES 4 128 18.77540 7.67652 6.17686 5.40045
URES 5 250 23.46930 7.69879 6.19478 5.41611
URES 6 432 28.16310 7.70813 6.20230 5.42268
As you can see, U has now been extrapolated to a supercell containing 432 atoms.
The value of U depends strongly on the extension of the projectors used in the calculation. If you want to use U in LMTO-ASA calculations you can use the keyword pawujrad in the ujdet.in file to get grips of the U you want to use there. Just uncomment the line and add the ASA-radius of the specific atom e.g.
pawujrad 2.5
Running
rm ujdet.[ol]* ; ujdet > ujdet.log
gives now higher values in the column U_ASA than in the runs before (6.91 eV compared to 6.20 eV): For more localized projectors the U value has to be bigger.