Commit 01e372f4 authored by Daniel Wortmann's avatar Daniel Wortmann

Added more of v26 doco

parent b7d59964
# Relaxation
Here, we describe how to do a structural optimization.
Suppose, you have converged a charge density, e.g. of a 3 layer Cu film as described in the example inp-file. Now, our interest is to optimize the position of the topmost layer (the position of the central layer is of course fixed by symmetry). What you have to do:
## How to get forces
* change ` force =F ` to ` force =T ` in the inp-file, for those atoms, which you want to relax.
* change ` l_f=F ` to ` l_f=T ` in the inp-file, to allow the generation of Pulay-forces
* edit the line ` relax 000 001 ` at the end of the inp-file to allow relaxations in specific directions (here, ` 000` means no relaxation of the first atom, ` 001 ` means only in z-direction for the second).
* then run a few iterations, until the program stops with ` GEO: new inp created !`. This will happen, if the forces of two subsequent iterations do not differ more than 0.00001 htr/a.u. (This parameter should not be changed to ensure good convergence of the forces. It can, however, in cases of emergency be changed by creating an "eps_force"-file, which contains a different convergence parameter.)
When this is finished, you will notice that you have two new files in your working directory:
* `inp_new` containing a new guess for the atomic positions and
* `forces.dat` which stores the positions and forces of this optimization step (and, if you keep it, also of the next steps).
* `out` contains the different Pulay force terms (cf. [![Symbol - externer Link][2]Rici Yu et al., PRB 43, 6411 (1991)][2] and the files `force*` of the Fleur source code).
## How to start with a new geometry
Now there are different ways how to proceed. In any case, you should not forget to remove a couple of files, that depend on the geometry:
* wkf2, the file containing the step-function (which cuts out the muffin-tins from the interstitial), and
* broyd and broyd.7, containing the (in/out)-density history of this self-consistency cycle.
To be on the safe side, also the file "stars" can be removed, along with all "cdnXX", "pot[tot,coul]", "eig" and "tma[t,s]" files. Keep your directories clean.
One possibility is to reuse just a minimum of files and create a new starting density and proceed as usual:
* ` mkdir ../GEO2 `
* ` cp enpara kpts sym.out forces.dat ../GEO2 `
* ` mv inp_new ../GEO2/inp `
* ` cd ../GEO2 `
* set ` strho=T `
* run ` fleur.x ` etc.
In other cases, it might be useful to reuse the old charge density with the new inp-file:
* ` mkdir ../GEO2 `
* ` cp enpara kpts sym.out forces.dat cdn1 ../GEO2 `
* ` mv inp_new ../GEO2/inp `
* ` cd ../GEO2 `
* ` echo " F" > qfix `
* run ` fleur.x ` etc.
In this case, you did not create a new starting density, but instead -- via the "F" in the file "qfix" indicate to the program that this is an old charge density, which is reused for shifted atomic positions (and therefore, the program renormalized the charge density a bit differently than normal). After the first iteration you will find that the "F" in this file changed to a "T".
## Finishing the structural optimization
This procedure can be repeated over and over again:
1. converge a calculation
2. obtain forces and an "inp_new" file
3. set up a new geometry and go to 1
For once, you will notice that the "forces.dat" file gets longer and longer, containing the history of your structural optimization (like the broyd-files store your convergence of the electronic degrees of freedom). Again, a quasi-Newton scheme is used to minimize your forces. This has advantages, but can also lead to problems:
* if you are already in a parabolic region of your total-energy landscape, this scheme converges fast. Of course, this will depend on the number of degrees of freedom you have.
* if not, you might run into trouble. ` STOP 'bfgs: <p|y><0 ` is an indication of this case, but also very large relaxations (overlapping muffin-tins, atoms in the vacuum) can result. Then you should delete the "forces.dat" file.
If you delete the "forces.dat" file after every step, you will effectively do a steepest-descent minimization. Since the relation between force and displacement is in general unknown, you can enter a Debye-temperature in the inp-file ("thetad= 330.000" near the end of the file).
As in every good mixing method, you can also specify a mixing parameter for your quasi-Newton scheme (called "xa" in the inp-file, typically set to 2%).
When you encounter a message "GEO: Des woas" (Austrian for: "that's it"), you either converged your system very well (either all forces are smaller than the parameter "epsforce" of the inp-file, or the displacements in the last step were smaller than "epsdisp), or all forces vanished for some other reason: Maybe, you set
* ` force =F ` for all atoms where forces can occur or
* ` relax 000 ` for directions where forces can be expected.
In any case, you should check. If your forces are smaller than 0.001 htr/a.u., for most cases the structure can be considered relaxed.
[]: http://link.aps.org/doi/10.1103/PhysRevB.43.6411
# HeisenbergInteractionParametersJij
## Heisenberg interaction parameters Jij
Calculation of Heisenberg interaction parameters:
1) Collinear calculation
You will use as a starting point the charge density obtained in a
collinear calculation. You may start from a ferro- or antiferromagnetic
state. In the case of an ideal Heisenberg system it wouldn't matter
which one of the two you chose, but it is almost certain that
your system is not an ideal Heisenberg system.
Therefore, it is wise to chose the starting collinear state
which is closer to the ground state of your system.
In the initial collinear calculation take care that you don't use
symmetries which group magnetic atoms of the basis into the same atom-type;
all the magnetic atoms in the input file should be inequivalent!
Whichever symmetries are left, use them!
The files you should keep after the collinear calculation are:
inp, enpara, cdn1, sym.out (in any case, DO DELETE fl7para and kpts file!).
Rename cdn1 to rhomat_inp.
2) In the inp-file:
(a) l_noco=T,l_J=T
(b) As for any non-collinear calculation, ctail=F
(c) Define the number of k-points (nkpt) you want to use; [kpts][2] file will be
generated with k-points from the whole Brillouin zone.
(d) Define in the last line the number of spin-spiral vectors (nqpt)
you want to use; [qpts][3] file will be generated with spin-spiral vectors ("qss")
from the irreducible wedge of the Brillouin zone.
3) Create a nocoinp file (use the switches for Jij calculation)
4) You are now ready to run!
CAUTION:
-- In Jij calculation force theorem is used. This usually means that your
k-points set should be considerably larger than in a self-consistent calculation.
It is advisable to check the convergence of a magnon dispersion
curve with respect to the number of k-points you use before you
proceed with the calculation of Jij. This means: take a direction in the
reciprocal space (it is convenient to take a high-symmetry direction) and
calculate self-consistently a dispersion curve on a set of few qss
along this direction. Save this set in the qpts file (take care it has the right format!).
For this calculation use the same cone angle you will use for the
Jij calculation. This angle shouldn't be too big - usually 30 degrees works fine.
Now see which is the smallest set of k-points you have to use to reproduce this
curve with the Force theorem. To make this test easier, you can use the l_disp
switch in the nocoinp file. When set to true (along with the changes (a) and (b)),
this switch tells the program to use the Force theorem to calculate the dispersion
curve on a predefined set of qss (the qpts set you previously generated).
A too small k-points set can often give you a wrong ground state of the system!
--The accuracy with which Jij(r) are calculated for a given qss set depends on the distance
of the atom i and the shell in which atom j lies, since the real-space Jij(r) are calculated
from their Fourier transforms Jij(qss). If you want Jij for a larger number of shells,
make sure the number of qss is big enough.
More about the theory of these calculations, the convergence tests and the applications examples you can read in [Lezaic_Marjana.pdf][6] (6.9 MB)
[6]: http://darwin.bth.rwth-aachen.de/opus/volltexte/2006/1359/pdf/Lezaic_Marjana.pdf
# Calculation a Bandstructure
Starting with FLEUR version .24 (used for the NIC winter school 2006) there is a convenient new run mode that further simplifies the creation of band structure plots.
For a proper band structure plot, also the correct Fermi energy is desired. This leads to two options:
## The quick way.
In this way no special care is taken for the Fermi energy, it is obtained from the eigenvalues of our band structure. Due to the very selection of the k-points, this can deviate from the true Fermi energy you obtained in your self-consistency. In the case of semiconductors, your band structure as well as your self-consistency k-point set might contain the k-point comprising the highest occupied state, then there should not be any difference. In the case of metals, one will expect a larger shift. The new mode tries to guess a reasonable path in the Brillouine zone for the k-points. For many simple cases it will give a suitable "default bandstructure". However, by supplying a band_inp file you can also choose the k-points yourself.
In the new mode you tell FLEUR to create a new k-point set and calculate the eigenvalues. Therefore you either need to remove all k-point dependent files,
> rm broyd* fl7para kpts pot* wkf2
or you continue the calculation in a new directory:
> mkdir band && cp inp cdn1 enpara sym.out band && cd band
As a next step, you modify the switch `band` in the inp file, changing the line from
iplot=F,score=F,plpot=F,band=F
to
iplot=F,score=F,plpot=F,band=T
That's it, now you can run FLEUR,
> fleur.x
DOS OK
and see newly created files `bands.1` and `band.gnu`. The latter one is a GNUplot script that can be used to easily create a postscript file: (Note that you have to adjust the Fermi-energy as showed below.
> gnuplot < band.gnu > band.ps
that can be viewed e.g. by
> gv band.ps
For the Silicon setup used in this guide, it looks like this:
![][3]
By the way, if you watch your `kpts` file,
> head -1 kpts
301 1.0000000000 you see that the number of k-points used is 301, since `nkpt= 300` was specified in the `inp` file.
## The proper way.
In the case of silicon, the Gamma point (comprising the largest occupied eigenvalue) is contained in our band structure path but not in the Monkhorst pack set used for self-consistency, therefore we expect a small change in the Fermi energy. To make matter clearer, the numbers for the Fermi energy below are taken from a cupper bulk calculation. Nevertheless, the procedure itself is general.
First remember the correct Fermi energy, so do
> grep -i Fermi out
fermi energy and band-weighting factors: FERMIE: FERHIS: Fermi-Energy by histogram:
EF_NEWTON: Adjust Fermi-Energy by Newton-Method.
--> new fermi energy : 0.281417 htr
after your self-consistency (This is only relevant for metals, for insulators or semiconductors the Fermi level is not printed explicitely in the out file). Then remove only the two files,
> rm fl7para kpts
and modify the switches `dos` and `ndir` in the first line of the inp file (same as above) as well as the `pot8` switch, changing it from
strho=F,film=F,dos=F,isec1=99,ndir= 0,secvar=F
:
:
vchk=F,cdinf=F,pot8=F,gw=0,numbands= 0
:
:
iplot=F,score=F,plpot=F,band=F to
strho=F,film=F,dos=T,isec1=99,ndir= 1,secvar=F
:
:
vchk=F,cdinf=F,pot8=T,gw=0,numbands= 0
:
:
iplot=F,score=F,plpot=F,band=T Now run FLEUR again
> fleur.x
DOS OK and grep for the new Fermi energy:
> grep -i Fermi out
FERMIE: FERHIS: Fermi-Energy by histogram:
EF_NEWTON: Adjust Fermi-Energy by Newton-Method.
--> new fermi energy : 0.271041 htr
There is a difference of `0.281417 htr - 0.271041 htr = 0.010376 htr = 0.282227 eV` that the correct Fermi energy is higher, i.e. with respect to the Fermi level the newly calculated eigenvalues should be lowered by this number. For this, edit the last line of the GNUplot script `band.gnu` from
"< awk '{print $1,$2+0.00}' bands.1" w p pt 7 ps 0.5
to
"< awk '{print $1,$2-0.282}' bands.1" w p pt 7 ps 0.5
Now you can continue as above,
> gnuplot < band.gnu > band.ps && gv band.ps
[3]: http://www.flapw.de/pm/datapool/sbsg/band/siband_2.jpg
# SymmetryOfBands
When doing band structure plots, FLEUR identifies the symmetry of the eigenvectors. Though not implemented for all point-group symmetries, the common cases are present. This shall be illustrated with a band structure of fcc Nickel.
# Band structure
fcc Nickel has space group [225][2] or Fm-3m. Its band structure of majority spin along the Delta line (i.e. Gamma-X) is discussed here. Starting from the file dosinp given below, the files bnd1.dat and bnd2.dat were obtained with the command
band3xmgr 15 0.314191
and bnd1.dat is shown below. The k-points Gamma and X have point group symmetry [0h][3] and [D4h][4] comprising 10 characters each, while the points on the Gamma line have point group symmetry [C4v][5] with five characters.
![][6]
The Fermi energy is marked with a dashed red line.
# FLEUR output
The symmetry data is placed in two files: fort.444 and dosinp.
When performing a band structure calculation, the file fort.444 is generated (visit the link at the bottom of the page) which contains the symmetry informations of each k-point, but only if it differs from the previous k-point. Therefore we get three tables, for Gamma, for the first point on the Delta line, and for X. Each section begins with the text "Character table for k: ". You can compare this to the tables linked above. These table also contains the degeneracy of eigenvalues.
In the file dosinp, you have have for each k-point a section listing all eigenvalues followed by its symmetry kind. (About the other information in that line, you can read here: [Calculating the Density of States][6].) For instance for the dosinp file given below starting at line 823 (refering to point number 20),
0.377277849596E+00 0.000000000000E+00 0.000000000000E+00 20 0.196078431373E-01
0.07324 1 0 0.00000
0.43875 0.08257 0.02693 0.00002
0.21743 4 0 0.00000
0.00000 0.00000 0.92980 0.00284
0.26096 5 0 0.00000
0.00000 0.00567 0.95576 0.00226
0.26096 5 0 0.00000
0.00000 0.00567 0.95576 0.00226
0.26498 1 0 0.00000
0.01314 0.00646 0.92652 0.00301
0.29293 3 0 0.00000
0.00000 0.00000 0.98334 0.00085
you get the symmetry kinds 1, 4, 5, 5, 1, 3 referring to the character listed in tables in file fort.444.
With this data, you can add some information to the plot. Below the same band structure is shown, with multiplicity and symmetry info for the eigenvalues at the Gamma point and an arbitrarily chosen k-point on the Delta line marked with the vertical red line (to which the cited data a few lines above refer to). The kind of band was read from the dosinp file, name and multiplicity from fort.444.
![][7]
This provides the proper way to detect band crossings, though it has not been done in the plot above.
# Footnote
Example files used:
* [![Symbol - externer Link][9]dosinp][9]
* [![Symbol - externer Link][10]fort.444][10]
Comments:
* The contents in fort.444 are printed twice. (For each spin, but of course the information is the same. You may ignore this.)
* In the given fort.444, the character table is not printed again for Gamma point with spin down (small bug).
* In the given fort.444, the band names at the X point are not yet implemented but named 'unknown'.
[2]: http://www.webelements.com/webelements/elements/text/Ni/xtal.html
[3]: http://www.cryst.ehu.es/cgi-bin/rep/programs/sam/point.py?sg=221
[4]: http://www.cryst.ehu.es/cgi-bin/rep/programs/sam/point.py?sg=123
[5]: http://www.cryst.ehu.es/cgi-bin/rep/programs/sam/point.py?sg=99
# SymmetryOfBands
When doing band structure plots, FLEUR identifies the symmetry of the eigenvectors. Though not implemented for all point-group symmetries, the common cases are present. This shall be illustrated with a band structure of fcc Nickel.
# Band structure
fcc Nickel has space group [![Symbol - externer Link][2]225][2] or Fm-3m. Its band structure of majority spin along the Delta line (i.e. Gamma-X) is discussed here. Starting from the file dosinp given below, the files bnd1.dat and bnd2.dat were obtained with the command
band3xmgr 15 0.314191
and bnd1.dat is shown below. The k-points Gamma and X have point group symmetry [![Symbol - externer Link][3]0h][3] and [![Symbol - externer Link][4]D4h][4] comprising 10 characters each, while the points on the Gamma line have point group symmetry [![Symbol - externer Link][5]C4v][5] with five characters.
![][5]
The Fermi energy is marked with a dashed red line.
# FLEUR output
The symmetry data is placed in two files: fort.444 and dosinp.
When performing a band structure calculation, the file fort.444 is generated (visit the link at the bottom of the page) which contains the symmetry informations of each k-point, but only if it differs from the previous k-point. Therefore we get three tables, for Gamma, for the first point on the Delta line, and for X. Each section begins with the text "Character table for k: ". You can compare this to the tables linked above. These table also contains the degeneracy of eigenvalues.
In the file dosinp, you have have for each k-point a section listing all eigenvalues followed by its symmetry kind. (About the other information in that line, you can read here: [Calculating the Density of States][6].) For instance for the dosinp file given below starting at line 823 (refering to point number 20),
0.377277849596E+00 0.000000000000E+00 0.000000000000E+00 20 0.196078431373E-01
0.07324 1 0 0.00000
0.43875 0.08257 0.02693 0.00002
0.21743 4 0 0.00000
0.00000 0.00000 0.92980 0.00284
0.26096 5 0 0.00000
0.00000 0.00567 0.95576 0.00226
0.26096 5 0 0.00000
0.00000 0.00567 0.95576 0.00226
0.26498 1 0 0.00000
0.01314 0.00646 0.92652 0.00301
0.29293 3 0 0.00000
0.00000 0.00000 0.98334 0.00085
you get the symmetry kinds 1, 4, 5, 5, 1, 3 referring to the character listed in tables in file fort.444.
With this data, you can add some information to the plot. Below the same band structure is shown, with multiplicity and symmetry info for the eigenvalues at the Gamma point and an arbitrarily chosen k-point on the Delta line marked with the vertical red line (to which the cited data a few lines above refer to). The kind of band was read from the dosinp file, name and multiplicity from fort.444.
![][7]
This provides the proper way to detect band crossings, though it has not been done in the plot above.
Comments:
* The contents in fort.444 are printed twice. (For each spin, but of course the information is the same. You may ignore this.)
* In the given fort.444, the character table is not printed again for Gamma point with spin down (small bug).
* In the given fort.444, the band names at the X point are not yet implemented but named 'unknown'.
[2]: http://www.webelements.com/webelements/elements/text/Ni/xtal.html
[3]: http://www.cryst.ehu.es/cgi-bin/rep/programs/sam/point.py?sg=221
[4]: http://www.cryst.ehu.es/cgi-bin/rep/programs/sam/point.py?sg=123
[5]: http://www.cryst.ehu.es/cgi-bin/rep/programs/sam/point.py?sg=99
[6]: http://www.flapw.de/pm/datapool/sympsi/fleursym_ni_up_gamma_x_up_1.png ""
[7]: http://www.flapw.de/pm/datapool/sympsi/fleursym_ni_up_gamma_x_up_2.png ""
# Density Of States
To calculate the density of States different switches and modes are present in the inp-file:
* Most important, `dos` hast to bet set to .true.
strho=F,film=F,dos=T
* The energy window, where the DOS is calculates is specified at the end of the inp-file:
emin_dos= -0.50000,emax_dos= 0.50000,sig_dos= 0.00050
Note, that this is in Hartree units and not referred to some Fermi-level. The DOS is broadened by a Gaussian with the width ` sig_dos ` (also in htr.).
* A little bit nicer DOS' can be produced with the tetrahedron (3D) or triangular (2D) method. For films this is used whenever the k-point set allows it (corner points are included). For bulk, you have to create a k-point set with ` tria=T ` in the inp-file.
* Finally, there is the switch ` ndir ` in the first line of the inp-file. it is used trigger the following:
* Producing raw-data for a DOS in the a dosinp-file
* Integrated DOS programm gives DOS.x files
* Orbital decomposition of the DOS
Note that you have to delete the fl7para and kpts files if you want to calculate the DOS.
# Producing Raw-data For A DOS
1. Setting dos=t and ndir=0 in the inp-file switches into the DOS-mode of FLEUR. In this mode a dosinp-file is created which contains the eigenvalues and thier relative weights for all kpts.
2. Similarly, one obtains a vacDOS file containing the LDOS in the vacuum by setting vacdos=t.
3. By setting `dos=t` and `ndir>0` one obtains a dosinp-file with additional symmetry-information which is used for band-structure plotting
# IntegratedDOSProgramm
## total & partial DOS
By setting dos=t and ndir<0 the integrated DOS program is used. Three additional parameters in the inp file have to be specified in this case. The two energies emin\_dos and emax\_dos specify the energy-range in which the DOS is evaluated. The parameter sig_dos specifies the broadening, as shown below (where the last few lines of the input file are given.
vacdos=T,layers= 1,integ=F,star=T,nstars= 2
50
iplot=F,score=F,plpot=F
0 0.000000 0.000000,nnne= 0,pallst=F
xa= 2.00000,thetad= 300.00000,epsdisp= 0.00010,epsforce= 0.00010
relax 000
emin_dos= -1.000,emax_dos= 0.00000,sig_dos= 0.00050
Three different negative values for ndir are implemented at the moment:
* ndir=-1: In this mode FLEUR calculates the DOS-output and stops without creating a new charge density. This is analogous to the old DOS-mode
* ndir=-2: In this mode FLEUR calculates the DOS-output and creates a new charge density afterwards. This is done only if it=itmax, so that one can use this setting to create a DOS during self-consistency iterations.
* ndir=-3: In this mode the Orbital decomposition of the DOS is calculated
For information about the output see the documentation of the DOS.x file.
## vacuum DOS
If one wishes to calculate the DOS in selected parts of the vacuum region, the parameter ` vacdos=T ` has to be set. In combination with `dos=T` and `ndir=-1`, this will generally produce a file VACDOS.[1,2] which has a format as describec in the documentation of the DOS.x file.
This is done by integrating over the 2-D unit cell and in z-direction one can either integrate (`integ=T`) up several layers (defined by `layers=` in the inp-file), or choose planes, where the DOS is evaluated.
* If `integ=F`, the line below should provide as many integer numbers, as the number of `layers` that has been defined. The z-position of a layer with number N is given by:
dvac N
z = ---- + -- (a.u.)
2 10
i.e. in the sample inp the layer # 50 is at 11 a.u. or about 3.7 Angstroem above the topmost Cu atom. The format is (in FORTRAN) ` (20(i3,1x))`, i.e. a maximum of 20 planes can be defined. By default, the largest number (N) that can be entered is 250.
Notice, that there is a parameter ` layerd ` in the fl7para file. Typically, this should be increased to the number of layers.
* If integration is choosen, then the next line in the inp-file defines pairs of z\_low and z\_up for integration. Again the pairsare given in integers as described above. The format here is ` (10(2(i3,1x),1x)) `, i.e. up to 10 regions can be integrated.
* If you are interested in the lateral modulation of the local DOS in the vacuum, you can plot the (symmetrized) Fourier-components of the DOS by setting `star=T`, as long as you are within 10 a.u. from the vacuum matching plane (i.e. dvac/2). The number of Fourier-components you want so see in your VACDOS.[1,2] file is defined by `nstars=`.
* Many more features are built in, but usually the ones described above should suffice. In any case they are mainly tested within collinear calculations (without SOC), so be careful.
# Orbital Decomposition Of The DOS
In cases it might be useful to obtain also the orbital decomposed density of states of a certain atom. To get a DOS.x file with this information, you have to specify
dos=T
and
ndir=-3
in the inp file. Additionally, a file named `orbcomp` has to be present in the working directory, where the number of an atom (counted in the order as it appears in the inp file) is specified. Then, the DOS-file contains information about the selected atom (not atom-type!) in the following form (the energy is in eV):
energy total-DOS selected atom-type
s p p p d d d d d f f f f ...
x y z xy yz zx x²-y² z² x³ y³ z³ x²y
(Note: you have to provide k-points in the whole Brillouin-zone to get a correct DOS here! Proper symmetrization will be implemented in future releases.) The other f-components can be found in the file ek_orco.[1,2], where also the orbital composition is given for each k-point. To get this file,
cdinf=T
should be set in the inp-file.
### rotated coordinate frames
Since the decomposition in spherical harmonics is per default performed in the global coordinate frame, it can be sometimes useful to have a tool to perform the decomposition in a rotated coordinate frame. For this, another file `orbcomprot` has to be specified that contains three lines, specifying the Euler-angles \alpha, \beta and \gamma for the rotation.
## layer resolved DOS
If a zero is specified in the `orbcomp` file, or no such file has been provided, then the DOS.x files will contain the information about the layer resolved DOS. This is especially meaningful for a film-calculation, where the film can be divided into layers; the boundaries of the layers are defined as (x,y)-planes that intersects halfway between two neighbouring atoms of different height in the film. The results are listed in the following form:
energy total-DOS layer interstitial muffin-tin
contribution to layer
1, 2, ... n 1, 2, ... n 1, 2, ... n
This diff is collapsed.
This diff is collapsed.
# LDAu
To introduce a Hubbard U for the description of strongly correlated electrons (e.g. in transition-metal oxides or f-electron systems) it is possible to use the LDA+U method (see e.g. [Anisimov et al. J.Phys.: Condens. Matter 9 (1997) p. 767-808][2]). It was implemented similar to [Shick et al., Phys. Rev. B 60, 10763 (1999)][3] but without the pseudo-perturbative treatment described there. For the input of the U and J parameters for an atom the second line (previously containing the local symmetry) is used:
***************************
Gd 64 ...
&ldaU l=3,u=6.7,j=0.7 /
This is a namelist input, the format is free. The 3 parameters, l, u and j specify the orbital (l=0 for s, 1 for p, 2 for d and 3 for f), the Hubbard U and exchange interaction parameter J (given in eV). The latter two parameters can be taken from literature or be calculated [Solovyev and Dederichs, Phys. Rev. B 49, 6736 (1994)][4]. The density matrix is stored on a file [n\_mmp\_mat] and should be kept with the charge density.
For vi-users, a convenient way to add LDA+U to all atoms of a specific type is (in the above example)
:g/Gd 64/+1 s/ /\&ldaU l=3,u=6.7,j=0.7 \//
Since the double-counting correction of the LDA+U method can be formulated in different ways, the atomic and the "around mean-field" limes, a switch has been introduced in the line
&ldaU l=2,u=4.0,j=2.0,l_amf=T/
If l_amf=T is set, the around mean-field limes is used, otherwise the atomic limes is taken. For practical purposes the difference is, that in the atomic limes (default in the version 0.22) states that are more than half occupied are lowered in energy, while in the around mean-field limes states that are more occupied than the average are lowered in energy. For a comparson of the two methods and another variant (not included) see e.g. [Petukhov et al., PRB 67, 153106 (2003)][6].
**Using LDA+U:**
(1) *Setting up the density matrix:* When you have converged a charge density with LDA or GGA, run a single iteration with `&ldaU ...` in the inp-file. A line in the output will inform you that LDA+U has been skipped, since no density-matrix was found, but together with the charge-density a `n_mmp_mat` file will be created, that contains the density-matrix to start with. Then delete the files `broyd` and `broyd.*`.
(2) *Choosing the right executable:* Note, that if the density matrix contains off-diagonal elements you will get a complex-hermitian Hamiltonian, e.g. if a orbital moment arises. Then, the program should be compiled without inversion, even if inversion-symmetry is present in the lattice.
(3) *Converging the calculation:* If you now run the executable, the charge-density and density-matrix have to be converged again. By default, the density-matrix is converged separately using straight-mixing with the mixing parameter given in the last line of the `n_mmp_mat`-file. If you delete this line, the density-matrix is mixed together with the charge-density in the `broyd`-files. If you switch from straight- to broyden-mixing, delete these files.
(4) *Finishing:* Monitor the convergence of the charge-density and density-matrix by e.g. `grep dist out`. Note that both densities have to be converged, which is sometimes not trivial. Also note, that the final result may depend on the starting point, e.g. different orbital momenta can be stabilized if you calculate an isolated atom starting from different density-matrices.
**Limitations and Problems**
* Only one set of orbitals (one `l`-value) can be specified per atom (anyway, you should use LDA+U with care and only for those cases, where it is really needed).
* If local orbitals are present with the same `l`-value as the one where the U acts on, the method will not discriminate between the local orbital and the normal LAPW's.
* Forces calculated for an atom, where the LDA+U method is applied, are currently inaccurate (of course this depends on the variation of the occupation matrix with the atomic displacement). See the discussion in [Tran et al., Comp. Phys. Comm. **179**, 784 (2008)][7].
* In principle, there exists also a spin-offdiagonal part of the density-matrix, which is not implemented so far (P. Novak, A. Shick et al.)
[2]: http://dx.doi.org/10.1088/0953-8984/9/4/002
[3]: http://link.aps.org/abstract/PRB/v60/p10763
[4]: http://link.aps.org/abstract/PRB/v49/p6736
[6]: http://link.aps.org/abstract/PRB/v67/e153106
[7]: http://dx.doi.org/10.1016/j.cpc.2008.06.015
Guide to Local Orbitals (v26)
==========================
## When do I need LO's?
When the accurate description of an atom requires two 'valence' bands with equal symmetry (e.g. 3p and 4p in Ti compounds) then the use of a 'second window' (energy panel) or of LO's is recommended. There are three categories of problems that might require additional local orbitals:
### (a) Atoms that tend to give raise to 'ghost-band' problems
Some materials with high-lying core states (say, 5p in Cs) tend to produce ghost-band, i.e. the unwanted appearance of these high-lying states in the valence window. Also in more robust atoms, the use of small muffin-tin radii (that makes the basis more flexible) can give such problems. A typical precursor of a ghost-band is the output 'n eigenvalue(s) below the energy e\_low', with n >= 1 and a reasonably chosen lower energy bound e\_low. When the ghost-band appears in the valence window, i.e. its energy is higher than e\_low, the charge-density distance 'jumps' to a large value and one of the energy parameters drops to a value near e\_low. The band associated with this energy parameters should then be described with LO's.
### (b) Low-lying bonding orbitals
Sometimes an atom binds with low- and high lying states of equal symmetry (e.g. in Ti oxides, the 3p binds to the oxygen while the 4p is involved in metallic states). Then one of these states should be described with LO's.
### (c) Extreme spin-splitting of f-bands
An extremely spin-split valence band, e.g. the 4f of Gd, requires two energy parameters to be described correctly. Here, LO's can be used for one spin channel.
## How do I calculate with LO's?
Having identified the atom and state that should be calculated with a LO, two steps are necessary to switch from a normal to a LO calculation:
### (a) Specify the number and l-character of the LO's (nlo and llo) and subtract them from the core-levels and adjust the number of electrons in the valence band.
First, Identify the atom in the input file. Say, an iron has a very small muffin-tin radius and gives 3p ghost-bands:
```
| ************************* |
| Fe 26 7 8 461 1.700000 0.023000 |
| |
| 1,force =F,nlo= 0,llo= |
```
* Now, we want a single LO: set nlo= 1
* it should be a p-orbital: set llo= 1 (use the l-quantum number)
* subtract two core-levels: 7 -> 5 (the 3p1/2 and 3p3/2)
This gives us:
```
| ************************* |
| Fe 26 5 8 461 1.700000 0.023000 |
| |
| 1,force =F,nlo= 1,llo= 1 |
```
Finally, add the electrons to the valence band, in our example change from:
```
| Window # 1 |
| -0.50000 0.80000 8.00000
```
to:
```
| Window # 1 |
| -2.00000 0.80000 14.00000 |
```
(the 3p contains 6 electrons, therefore 8+6=14. Here, we also had to
change e_low from -0.5 to -2.0, so that the 3p fits in the valence
window.)
### (b) Choose the energy parameter and fix these values or skip the states in the calculation of the energy parameters.
Now set the energy parameters in enpara. (If you start a calculation, this should be done automatically, but check, whether it did, what you intended!) Taking the Fe-example from above, your enpara file probably contained energy parameters like:
```
| energy parameters for window 1 spin 1 |
| atom s p d f |
| --> 1 0.13301 0.21323 0.24802 0.24094 change: TTTT skiplo: 0 |
```
and you expect the 3p somewhere near -1.7 htr., so you add:
```
| --> 1 0.13301 0.21323 0.24802 0.24094 change: TTTT skiplo: 3 |
| --> lo -1.70000 |
| --> change T |
```
Here, in the last line we specified 'change T ', i.e this LO energy parameter should be adjusted to the center of gravity of the p-band. But, since we have two p-bands, we also want two energy parameters ! Therefore we specified in the first line 'skiplo: 3', that means: skip the first three bands for the evaluation of the 4p energy parameter. Alternatively, we could have set 'change: TFTT skiplo: 0', then the first p-energy parameter (0.21323) would have kept fixed at this value. In general, you have to avoid two equal energy parameters. Otherwise, the basis gets overcomplete and the program crashes.
### Summary:
* set number and type of LO's
* reduce the number of core levels
* add valence electrons
* adjust lower energy bound of window
* specify energy parameter
* choose, which energy parameters to change