pwscfLDA+U机理和操作的详细介绍

LSDA+U practical session
In this session we will learn how to run a LSDA+U calculation and how to compute the U from a linear response approach. We will study three examples: Fe, FeO and NiO.
Fe is a metal. For this material the DFT+U approach is not particularly useful. However we will study this system to learn how to compute the U.
FeO and NiO are, instead, “Mott” and “charge-transfer” insulators (correlated systems). For these systems the use of DFT+U can yield much more accurate results than “standard” approxi-mations to the DFT exchange-correlation energy functionals.
FeO
In FeO_ldau you will find three scripts to run: run_FeO_NO_U is to run a regular LDA calcula-tion, run_FeO_U and run_FeO_U2 for LDA+U calculations.
run_FeO_NO_U is the input file for the antiferromagnetic (AFM) phase of this system. The “system” namelist of this input file is:
&system
ibrav= 0, celldm(1)=8.19, nat= 4, ntyp= 3,
ecutwfc = 30.0, ecutrho = 240.0,
occupations='smearing', smearing='mp', degauss=0.02,
nspin=2,
starting_magnetization(2)= 0.5, starting_magnetization(3)=-0.5,
lda_plus_u = .true.,
Hubbard_U(2) = 1.d-10,
Hubbard_U(3) = 1.d-10,
/
lda_plus_u is set to .true. (necessary for DFT+U calculations) but the U on Fe atoms is very small (practically 0): 1.d-10. U won’t have any effect on the final result of the calculations; how-ever its finite value is needed for the code to print the occupation matrices in the output file. The occupation matrix is written as follows:
atom 4 Tr[ns(na)]= 6.6962905 -----> Tr n^I: total occupation of the d states of atom 4
atom 4 spin 1
eigenvalues: 0.3032823 0.3032823 0.3982914 0.4540635 0.4540635----> occupation numbers. in this case two eigenvectors are almost half-filled, the others have lower occupations.
eigenvectors ---> for d states the basis set is: z^2, xz, yz, x^2-y^2, xy. Eigenvectors are given in the same order as eigenvalues. So, for example, the state: 1/sqrt(3)(xz+yz+xy), third eigenvector of the list, (which points along the [111] direction) has an occupation of 0.3982914, third eigen-value of the list.
1 -0.2491645 -0.1536785 0.1132056 -0.9486023 -0.0404729
2 -0.948602
3 -0.0419923 -0.112093
4 0.249164
5 -0.1540857
3 0.0000000 0.5773503 0.5773503 0.0000000 -0.5773503
4 0.1726237 -0.0309771 -0.6775086 -0.0909491 -0.7084856
5 -0.0909491 0.8002041 -0.4269290 -0.1726237 0.3732751
occupations ---> this is the 5x5 occupation matrix.
0.309 -0.012 -0.012 0.000 -0.024
-0.012 0.432 -0.017 -0.020 0.017
-0.012 -0.017 0.432 0.020 0.017
0.000 -0.020 0.020 0.309 0.000
-0.024 0.017 0.017 0.000 0.432
This example relates to a DFT calculation. Since the material is antiferromagnetic the occupation eigenvalues of the the two Fe ions are like follows:
atom 3 spin 1
eigenvalues: 0.9879041 0.9879041 0.9976713 0.9979445 0.9979445
atom 3 spin 2
eigenvalues: 0.3032823 0.3032823 0.3982914 0.4540635 0.4540635
atom 4 spin 1
eigenvalues: 0.3032823 0.3032823 0.3982912 0.4540633 0.4540633
atom 4 spin 2
eigenvalues: 0.9879041 0.9879041 0.9976713 0.9979445 0.9979445
Please notice the splitting of the 5 d states into two doubly degenerate groups (E_g symmetry) and 1 single state (A1g).
If a finite (4.3 eV) U is applied to the Fe ions (run_FeO_U) the eigenvalues of the occupation matrices are as follows:
atom 3 spin 1
eigenvalues: 0.9955893 0.9955893 0.9999277 1.0001248 1.0001248
atom 3 spin 2
eigenvalues: 0.1067696 0.2378999 0.2378999 0.5950639 0.5950639
atom 4 Tr[ns(na)]= 6.7643928
atom 4 spin 1
eigenvalues: 0.1067751 0.2379167 0.2379167 0.5952144 0.5952144
atom 4 spin 2
eigenvalues: 0.9955899 0.9955899 0.9999272 1.0001243 1.0001243
So the most relevant change wrt the U = 0 case (LDA) consists in the fact that the single state has now the lowest occupation (minority spin). This situation leaves the most occupied doublet states half filled and results in a metallic ground state (see DOS) contradicting the observed insulating character of FeO. A more “natural” (or, at least intuitive) choice for the occupation of the minor-ity spin manifold (where, nominally, only 1 electron is present) would be to fill the single state (A1g). This can be obtained by “suggesting” the system to do so from the input file (run_FeO_U2):
&system
ibrav= 0, celldm(1)=8.19, nat= 4, ntyp= 3,
ecutwfc = 30.0, ecutrho = 240.0,
occupations='smearing', smearing='mp', degauss=0.02,
nspin=2,
starting_magnetization(2)= 0.5, starting_magnetization(3)=-0.5,
lda_plus_u=.true.
Hubbard_U(2)=4.3, Hubbard_U(3)=4.3,
starting_ns_eigenvalue(3,2,2)=1.d0, starting_ns_eigenvalue(3,1,3)=1.d0
/
through starting_ns_eigenvalue the occupation matrix is forced to have a certain occupation on a specific state. The order of indexes is as follows: state, spin, atom kind. How do we know what state we want to force the occupation of? We need to look at the occupations of the first iteration of the run_FeO_U output and decide from there (in fact, the code will enforce the specified oc-cupation just after the first iteration):
atom 3 spin 2
eigenvalues: 0.1701462 0.1701462 0.2159735 0.2623305 0.2623305
So the “lone” state is the third one for atom 3 (type 2) spin 2. This is the state we would like to put the electron on.
Result from the calculation with the new input:
atom 3 spin 2
eigenvalues: 0.0979003 0.0979003 0.2586863 0.2586863 0.9927447
So the single state has become the most occupied one. Also, the corresponding ground state has a lower energy (please check this point).
This teaches us that, when DFT+U is used, the system may acquire some metastable solutions (corresponding to relatively close local minima of the energy). To reach the ground state of the system (global minimum) some insight on how localized electrons occupy atomic levels is needed
Computing the DOS you should be able to observe that in the latter case a gap opens between valence and conduction states as expected from experiments.
Fe
Although DFT+U is not very useful for bulk Fe, we will use this system to get familiar with the calculation of the Hubbard U.
Preliminarily (probably there won’t be time to) you can study what the effect of U is on bulk Fe just comparing the results (e.g., DOS) you get with DFT and DFT+U. To this purpose you can use the scripts: run_Fe_NOU and run_Fe_U.
Calculation of U.
U is computed as the second derivative of the energy with respect to the on-site occupation num-bers. A non-interacting term (due to re-hybridization of orbitals) is subtracted:
U =d 2E GGA d (n I )2−d 2E 0GGA d (n I )2
The second derivatives appearing in this expressions are computed with a linear-response ap-proach studying the response of the system (the variation of on-site occupations) to a perturba-tion on the localized states of a “Hubbard” atoms:
ΔV =
αϕ
m I ϕm
I m ∑If we define J I IJ
d n d αχ00=J I
IJ d n d αχ=we can write U as:
U =−d αI dn I +d αI dn 0I =χ0−1−χ−1()II So we aim to compute the two response matrices X and X 0. We proceed through the following steps:
1) run an scf calculation with no perturbation
2) starting from the potential saved from the scf calculation we perform a series of calculations for few values of alpha around 0.
3) studying the (linear) response of all the on-site occupations we construct the response matri-ces.
4) we invert them and take the difference. The diagonal elements will be the U parameters on different atoms.
The perturbation we apply should be isolated: we don’t want that perturbations on periodically repeated atoms interfere with each other. So we normally use supercells of increasing size (up to when convergence of U is reached). Data obtained from smaller cells can also be extrapolated to bigger ones in the hypothesis that further shells of neighbors do not participate significantly to the response.
Let’s start from the smallest supercell possible for BCC Fe: the simple cubic cell with 2 atoms/ cell. Let’s inspect the script. The main part of input file for the scf calculation is as follows:
&system
ibrav= 1, celldm(1)=5.42, nat= 2, ntyp= 2,
ecutwfc = 30.0, ecutrho = 360.0,
occupations='smearing', smearing='mp', degauss=0.01,
nbnd = 14,
nspin=2,
starting_magnetization(1)= 0.6
starting_magnetization(2)= 0.6
lda_plus_u = .true.
Hubbard_U(1)= 1.d-20
Hubbard_U(2)= 1.d-20
/
&electrons
mixing_beta = 0.7
conv_thr = 1.0d-9,
/
ATOMIC_SPECIES
Fe1 1. Fe.pbe-nd-rrkjus.UPF
Fe 1. Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
Fe1 0.0 0.0 0.0
Fe 0.5 0.5 0.5
We notice that it is a DFT+U calculation (lda_plus_u = .true.) with a very small U (1.d-20). As explained earlier this small U has no effect on the results but is necessary to have the code print the occupation matrices we will need to evaluate the response of the system to the external per-turbation. It is important to notice that the two atoms in the unit cell are described as of different kind: in fact one of them will be subject to the perturbative potential (if we had them described as of the same kind the perturbation would be added on both). In general the atom on which the per-turbation is applied needs to be described as of different type from the others (even though with the same pseudopotential) to which it would be equivalent otherwise.
After the scf run is complete we are ready to start the perturbative calculation:
ethr=`grep ethr results_Fe_sc1_ucalc_pbe/fe.scf.out |tail -1 |awk '{print $3}'`
## the above instruction extracts the threshold of the iterative diagonalization of the Hamiltonian ## from the last iteration of the scf run. This is used as the initial diagonalization threshold for ## the first iteration of the perturbed run (the one used to compute X0). See diago_thr_init below
rm -rf $temp1/*
mv $temp/* $temp1/ # we need to save the potential file from the unperturbed scf run on
# another directory
for a in 0.0 -0.05 0.05 -0.1 0.1 # this is the series of alpha values around 0
do
rm -rf $temp/* # we always start from the same potential (of the unperturbed run)
cp -r $temp1/* $temp/
# perturbed calculations (finite alpha)
rm -f fe.scf.in
cat > fe.scf.in << EOF
&control
pseudo_dir = '$pseudo',
outdir='$temp'
restart_mode='from_scratch'
verbosity = 'high',
prefix='fe',
/
&system
ibrav= 1, celldm(1)=5.42, nat= 2, ntyp= 2,
ecutwfc = 30.0, ecutrho = 360.0,
occupations='smearing', smearing='mp', degauss=0.01,
nbnd = 14,
nspin=2,
starting_magnetization(1)= 0.6
starting_magnetization(2)= 0.6
lda_plus_u = .true.
Hubbard_U(1)= 1.d-20
Hubbard_U(2)= 1.d-20
Hubbard_alpha(1)= $a # this is alpha, only applied on atom 1 /
&electrons
startingpot = 'file' # we start from potential and wfcs saved from the scf
startingwfc = 'file' # (unperturbed) run
diago_thr_init = $ethr # see above
mixing_beta = 0.7
conv_thr = 1.0d-9,
/
ATOMIC_SPECIES
Fe1 1. Fe.pbe-nd-rrkjus.UPF
Fe 1. Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
Fe1 0.0 0.0 0.0
Fe 0.5 0.5 0.5
K_POINTS automatic
8 8 8 0 0 0
EOF
pw.x < fe.scf.in > results_Fe_sc1_ucalc_pbe/fe.pert_$a.out
done
Obviously, since the two atoms are equivalent, we don’t need to perturb them separately. We can construct the response to a perturbation on the second one by symmetry.
If we look inside the file fe.scf.out and we search “Tr” we will find the total occupation of each atom. At the last iteration (self-consistency) we have:
atom 1 Tr[ns(na)]= 7.5336622
atom 2 Tr[ns(na)]= 7.5336535
The slight difference between these numbers (they should be strictly equal) is within the numeri-cal precision of our calculations.
Now we open the output file of the perturbed run with alpha = -0.1 (fe.pert_-0.1.out) and we look for the on-site occupation to study the change du to the perturbation. After the first iteration we have:
atom 1 Tr[ns(na)]= 7.6039784
atom 2 Tr[ns(na)]= 7.4771671
When self-consistency is reached we obtain:
atom 1 Tr[ns(na)]= 7.5605331
atom 2 Tr[ns(na)]= 7.5217903
So the first atom acquires more electrons than the second (the alpha is negative) but this effect is screened by electron-electron interactions (readjusting the charge density during the scf cycle), which is what we want to measure.
After the whole series of calculation is complete we are ready to compute U. For this we can “cd Ucalc” where we will use different scripts and codes.
first: issue the command “./grepnalfa_sc1_pbe”. This script will collect all the occupations for all the atoms from the output of the perturbed run and write them in files dn.at1.da.at2.dat where at1 is the atom whose response we are measuring, at2 is the atom where the perturbation was ap-plied. The file named is dn0.at1.da.at2.dat contains the response at the first iteration (to compute X0). The script grepnalfa_sc1_pbe will also create file_sc1 where the names of the dn*.dat files are collected.
Using xmgrace try to plot a dn.at1.da.at2.dat file and the corresponding dn0.at1.da.at2.dat. You should observe a linear behavior of n wrt alpha for both cases. Also the two lines should cross each other at alpha = 0.
The code that computes the response matrices and invert them is resp_mat.f90. Its input file resp_mat.in is constructed as follows:
&input_mat
ntyp = 1 # number of types of atoms
na(1) = 2 # number of atoms per type
nalfa = 5 # number of perturbations (alpha)
filepos = 'pos_sc1' # file containing the atomic positions and the vectors of the unit cell
back = 'neutral' # to add a neutralizing “background”. see PRB 71 35105 (2005)
filednda = 'file_sc1' # file containing the names of the files dn*dat
n1 = 4
n2 = 4 # number of unit cells in each direction to construct the supercell for
n3 = 4 # extrapolation (in this case a 4x4x4 supercell of the SC one)
&end
issue “./r.x < resp_mat.in”. a Umat.out file should be created containing the response matrices and the differences between their inverses that define U. In the file Umat.out search for the string: “ CHI0^-1 - CHI^-1 Matrix”. The interaction matrix is given in the format:
typ(at1) typ(at2) spin(1) spin(2) at1 at2 distance(at1,at2) U
The script ucalc_sc1_pbe.j does all these operations necessary to compute U automatically.
The whole procedure can now be repeated starting from calculations in bigger supercells. Go back to the parent directory and launch run_Fe_bcc8_ucalc_pbe (probably there won’t be time to do anything else....) that performs the perturbative calculation on a 2x2x2 BCC supercell with 8 atoms. Atomic species and positions are as follows:
ATOMIC_SPECIES
Fe1 1. Fe.pbe-nd-rrkjus.UPF
Fe 1. Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
Fe1 0.00000 0.00000 0.00000
Fe 0.50000 0.00000 0.00000
Fe 0.00000 0.50000 0.00000
Fe 0.00000 0.00000 0.50000
Fe 0.50000 0.50000 0.00000
Fe 0.00000 0.50000 0.50000
Fe 0.50000 0.00000 0.50000
Fe 0.50000 0.50000 0.50000
Please notice that the atom where alpha is applied is always described as of different type from the others although with the same PP and always placed at the origin (0,0,0) of the unit cell. In this way symmetry is preserved as much as possible and calculations complete in shorter time. After repeating the same series of calculations starting also from a 2x2x2 SC supercell with 16 atoms we find the following results from the convergence to bigger and bigger unit cells:
The black line describes the extrapolation of data from the 2 atoms SC unit cell, the red from the 8 atoms BCC cell, the blue from a 16 atoms SC super cell. The plateau of the blue line is the converged value of U that can be approximated by extrapolating the 8 atoms BCC cell results. NiO
Move to NiO_ldau. Even though there won’t be time to run an actual calculation it is instructive to study the input file (run_NiO_r16_ucalc) to compute the U in this case. For this system we want to evaluate U both for Ni and O. Things to notice:
1) Formally there are two kinds of Ni atoms; however they only differ for having opposite mag-netization, so only one of them will be perturbed. The response to a perturbation on the second is obtained using the symmetry between the two spin.
2) When computing the response of the perturbation on oxygen a new scf calculation is per-formed with the perturbed oxygen translated to the origin of the unit cell. This is done to pre-serve the largest possible number of symmetries.
3) the atom where the perturbation is applied is of different kind from the others described by the same PP.
In Ucalc/ the input file resp_mat.in looks as follows:
&input_mat
ntyp = 2
na(1) = 8 # Ni atoms The order of species should be the same as in the pos_nio_r16
na(2) = 8 # O atoms file
nalfa = 5
magn = .true. # because we want to specify the spin of Ni (see file pos_nio_r16)
filepos = 'pos_nio_r16'
back = 'no'
filednda = 'file.nio.r16'
n1 = 5
n2 = 5
n3 = 5
&end
the pos_nio_r16 should be as follows:
10.0 10.0 0.0
10.0 0.0 10.0
0.0 10.0 10.0
0.d0 0.d0 0.d0 1
0.5d0 0.5d0 0.d0 1 # spin up Ni
0.5d0 0.d0 0.5d0 1
0.d0 0.5d0 0.5d0 1
0.5d0 0.d0 0.d0 -1
0.d0 0.5d0 0.d0 -1 # spin down Ni
0.d0 0.d0 0.5d0 -1
0.5d0 0.5d0 0.5d0 -1
-0.25d0 0.25d0 0.25d0 0
0.25d0 0.25d0 0.25d0 0
-0.25d0 0.25d0 0.75d0 0 # O: non spin polarized
-0.25d0 0.75d0 0.25d0 0
0.25d0 0.75d0 0.25d0 0
0.25d0 0.25d0 0.75d0 0
-0.25d0 0.75d0 0.75d0 0
0.25d0 0.75d0 0.75d0 0
Notice that the O atoms have 0 spin (they are not magnetic) while Ni atoms have 1 or -1
In the file Umat.out search for the string: “ CHI0^-1 - CHI^-1 Matrix”. The interaction matrix is given in the format:
typ(at1) typ(at2) spin(1) spin(2) at1 at2 distance(at1,at2) U。