7. Programs for running an SCF cycle

In sections 7.1-7.9 we describe the main programs to run an SCF cycle as illustrated in figure 4.1.

1. LAPW0 (generates potential)

lapw0 computes the total potential $V_{tot}$ as the sum of the Coulomb $V_c$ and the exchange-correlation potential $V_{xc}$ using the total electron (spin) density as input. It generates the spherical part (l=0) as case.vsp and the non-spherical part as case.vns. For spin-polarized systems, the spin-densities case.clmup and case.clmdn lead to two pairs of potential files. These files are called: case.vspup, case.vnsup and case.vspdn, case.vnsdn.

The Coulomb potential is calculated by the multipolar Fourier expansion introduced by Weinert (81). Utilizing the spatial partitioning of the unit cell and the dual representation of the charge density [equ. 2.10], firstly the multipole moments inside the spheres are calculated (Q-sp). The Fourier series of the charge density in the interstitial also represent SOME density inside the spheres, but certainly NOT the correct density there. Nevertheless, the multipole moments of this artificial plane-wave density inside each sphere are also calculated (Q-pw). By subtracting Q-pw from Q-sp one obtains pseudo-multipole moments Q. Next a new plane-wave series is generated which has two properties, namely zero density in the interstitial region and a charge distribution inside the spheres that reproduces the pseudo-multipole moments Q. This series is added to the original interstitial Fourier series for the density to form a new series which has two desirable properties: it simultaneously represents the interstitial charge density AND it has the same multipole moments inside the spheres as the actual density. Using this Fourier series the interstitial Coulomb potential follows immediately by dividing the Fourier coefficients by $K^2$ (up to a constant).

Inside the spheres the Coulomb potential is obtained by a straightforward classical Green's function method for the solution of the boundary value problem.

The exchange-correlation potential is computed numerically on a grid. Inside the atomic spheres a least squares procedure is used to reproduce the potential using a lattice harmonics representation (the linear equations are solved with modified LINPACK routines). In the interstitial region a 3-dimensional fast Fourier transformation (FFT) is used.

The total potential $V$ is obtained by summation of the Coulomb $V_C$ and exchange-correlation potentials $V_{xc}$.

In order to find the contribution from the plane wave representation to the Hamilton matrix elements we reanalyze the Fourier series in such a way that the new series represents a potential which is zero inside the spheres but keeps the original value in the interstitial region.

The contribution to the total energy which involves integrals of the form $\rho*V$ is calculated according to the formalism of Weinert et al (82).

The Hellmann-Feynman force contribution to the total force is also calculated (Yu et al 91).

Finally, the electric field gradient (EFG) is calculated in case you have an L=2 term in the density expansion. The EFG tensor is given in both, the ``local-rotation-matrix'' coordinate system, and then diagonalized. The resulting eigenvectors of this rotation are given by columns.

For surface calculations the total and electrostatic potential at z=0 and z=0.5 is calculated and can be used as energy-zero for the determination of the workfunction. (It is assumed that the middle of your vacuum region is either at z=0 or z=0.5).

1. Execution

The program lapw0 is executed by invoking the command:

lapw0 lapw0.def or x lapw0

2. Dimensioning parameters

The following parameters are used (they are collected in file param.inc, but usually need not to be changed:

NCOM	number of lm components in charge density and potential representation; it must satisfy the following condition: NCOM+3 .gt. $\{$[number of $l,m$ with $m=0$] + [2 * number of $l,m$ with $m>0$]$\}$
NRAD	number of radial mesh points
NSTD	defines the angular grid points used in fitting the xc-potential inside spheres
LMAX1D	defines the angular grid points used in fitting the xc-potential inside spheres
LMAX2	highest L in the LM expansion of charge and potential
NSYM	order of point group

3. Input

The input is very simple. It is generated automatically by init_lapw, and needs to be changed only if a different exchange-correlation potential should be used:

------------------ top of file: case.in0 --------------------
TOT   13              MULT/COUL/EXCH/POT /TOT ;  VXC-SWITCH
NR2V    (R2V / NR2V)
0 0.0     (#of FK in E-field expansion, EFELD (Ry)
0         (list of K-vectors)

------------------- bottom of file ---------------------------

Interpretive comments follow:

line 1:

format(A4,I4)
switch, indxc

switch
	TOT	total energy contributions and total potential calculated
	POT	total potential is calculated, but not the total energy
	MULT	multipole moments calculated only
	COUL	Coulomb potential calculated only
	EXCH	exchange correlation potential calculated only
		NOTE: MULT, COUL, and EXCH are for testing only, whereas POT, saves some CPU time if total energy is not needed
indxc		index to specify type of exchange and correlation potential. Supported options include:
	1	Moruzzi, Janak, Williams(78)
	5	Perdew and Wang 92, reparameterization of Ceperly-Alder data, the recommended LDA option
	13	Generalized Gradient approximation (Perdew-Bourke-Ernzerhof 96)
	14	Generalized Gradient approximation (Perdew-Wang 91)
	12	Meta GGA (Perdew et al. 2000). You must have a case.inm_vresp file without renormalization to generate the required input files for this option. In addition you must use a very large GMAX (24).

line 2:

format(A4)
option

	NR2V	no additional output
	R2V	Exchange-correlation (case.r2v), Coulomb and total potentials (case.vcoul) are written as ($r^2 V$) to a file for plotting with lapw5

The following 2 lines are optional and can be omitted. They are used to introduce an electric field via a zig-zag potential (see J.Stahn et al. 2000):

line 3:

free format
IEFIELD, EFIELD

	IEFIELD	number of Fourier coefficients to model the zig-zag potential
	EFIELD	value of the electric field (mRy)
		Please note: The zig-zag potential is only implemented for the LDA functionals.

line 4:

free format

KFIELD

Indices of K-vectors to be used in the zig-zag E-field description. I.e.: 2,3,4 means that you model this field by using the 2-4 th K-vector in case.clmsum, which are usually vectors like (0,0,1), (0,0,2),...)

2. LAPW1 (generates eigenvalues and eigenvectors)

lapw1 sets up the Hamiltonian and the overlap matrix (Koelling and Arbman 75) and finds by diagonalization eigenvalues and eigenvectors which are written to case.vector. Besides the standard LAPW basis set, also the APW+lo method (see Sjöstedt et al 2000, Madsen et al. 2001) is supported and the basis sets can be mixed for maximal efficiency. If the file case.vns exists (i.e. non-spherical terms in the potential), a full-potential calculation is performed.

For structures without inversion symmetry, where the hamilton and overlap matrix elements are complex numbers, the corresponding program version lapw1c must be used in connection with lapw2c.

Since usually the diagonalization is the most time consuming part of the calculations, several options exist here. In WIEN2k we include highly optimized modifications of LAPACK routines. We call all these routines ``full diagonalization'', but we also provide an option to do an ``iterative diagonalization'' using a block-Davidson method (see Singh 89). This scheme starts from an old eigenvector (previous scf-iteration), needs additional memory and produces only approximate eigenvalues/vectors, but can be significantly faster than LAPACK, in particular if the ratio of matrix size to number of relevant (e.g. occupied) eigenvalues is large. In any case, it is recommended that convergence is checked by another scf-iteration with ``full diagonalization''. Often the best performance can be obtained by first running to a crude convergence, next doing a ``full diagonalization'' and finally continuing iteratively to self consistency. (e. g. use first: run_lapw -it -fc 10 and then run_lapw -it -fc 1)

Parallel execution (fine grain and on the k-point level) is also possible and is described in detail in Sec. 5.4.

1. Execution

The program lapw1 is executed by invoking the command:

x lapw1 [-c -up|dn -it -p -nohns -orb] or

lapw1 lapw1.def or lapw1c lapw1.def

In cases without inversion symmetry, the default input filename is case.in1c. For semi-core calculations the lapw1s.def file uses a case.in1s file and creates the files case.output1s and case.vectors. For the spin-polarized case lapw1 is called twice with uplapw1.def and dnlapw1.def. To all relevant files the keywords ``up`` or ``dn`` are appended (see the fcc Ni test case in the WIEN2k package). The switch -nohns skips the calculation of the nonspherical matrix elements inside the sphere. This may be used to save computer time during the first scf cycles.

2. Dimensioning parameters

The following parameters (collected in file param.inc_r or param.inc_c) are used:

KMAX1START	a first guess for the largest h,k,l Fourier components of the potential
KMAX2START
KMAX3START
LMAX	highest l+1 in basis function inside sphere (consistent with input in case.in1)
LMMX	number of LM terms in potential (should be at least NCOM-1)
LOMAX	highest l for local orbital basis (consistent with input in case.in1)
NGAU	number of Gaunt coefficients for the non-spherical contributions to the matrix elements
NKPTSTART	a starting guess for the total number of k-points in irreducible wedge of Brillouin zone
NMATMAX	maximum size of H,S-matrix (defines size of program, should be chosen according to the memory of your hardware!)
NRAD	number of radial mesh points
NSLMAX	highest l+1 in basis functions for non-muffin-tin matrix elements (consistent with input in case.in1)
NSYM	order of point group
NUME	maximum number of energy eigenvalues per k-point
NVEC1	defines the largest triple of integers which define reciprocal
NVEC2	K-vectors when multiplied with the reciprocal Bravais matrix
NVEC3
NLOAT	max number of LOs for one $l$-quantum number

3. Input

Below a sample input is shown for $TiO_2$ (rutile), one of the test cases provided in the WIEN2k package. The input file is written automatically by LSTART, but was modified to set APW only for Ti-3d and O-2p orbitals. In addition UNIT was changed to 5 and k-points where inserted by hand for bandstructure plotting.

------------------ top of file: case.in1 --------------------
WFFIL       (WFPRI,WFFIL,SUPWF ; wave fct. print,file,suppress
  7.500   10    4  (R-mt*K-max; MAX l, max l for hns )
  0.30    5 0 (global energy parameter E(l), with 5 other choices,   LAPW)
 0   -3.00      0.020 CONT 0   ENERGY PARAMETER for s,               LAPW
 0    0.30      0.000 CONT 0   ENERGY PARAMETER for s-local orbital, LAPW-LO
 1   -1.90      0.020 CONT 0   ENERGY PARAMETER for p                LAPW
 1    0.30      0.000 CONT 0   ENERGY PARAMETER for p-local orbitals LAPW-LO
 2    0.20      0.020 CONT 1                                         APW   
  0.20    3 0 (global energy parameter E(l), with 1 other choice,    LAPW)
 0   -0.90      0.020 STOP 0                                         LAPW
 0    0.30      0.000 CONT 0                                         LAPW-LO
 1    0.30      0.000 CONT 1                                         APW
K-VECTORS FROM UNIT:5   -9.0      1.5          IUNIT, Emin, Emax
GAMMA         0    0    0    2 1.00 
              0    0    1   10 2.00
              0    0    2   10 2.00
              0    0    3   10 2.00
              0    0    4   10 2.00
Z             0    0    5   10 1.00
END
------------------- bottom of file ------------------------

Interpretive comments follow:

line 1:

format(A5)
switch

	WFFIL	standard option, writes wave functions to file case.vector (needed in lapw2)
	SUPWF	suppresses wave function calculation (faster for testing eigenvalues only)
	WFPRI	prints eigenvectors to case.output1 and writes case.vector (produces long outputs!)

line 2:

free format
rkmax, lmax, lnsmax

rkmax		$R_{mt} * K_{max}$ determines matrix size (convergence), where Kmax is the plane wave cut-off, Rmt is the smallest of all atomic sphere radii. Usually this value should be between 5 and 9 (APW+lo) or 6 - 10. (LAPW-basis) ($K_{max}^2$ would be the plane wave cut-off parameter in Ry used in pseudopotential calculations.) Note that d (f) wavefunctions converge slower than s and p. For systems including hydrogen with short bondlength and thus a very small $R_{mt}$ (e.g. 0.7 a.u.), RKmax = 3 might already be reasonable, but convergence must be checked for a new type of system.
		Note, that the actual matrix size is written on case.scf1. It is determined by whatever is smaller, the plane wave cut-off (specified with RKmax) or the maximum matrix dimension NMATMAX, (see previous section).
lmax		maximum l value for partial waves used inside atomic spheres (should be between 8 and 12)
lnsmax		maximum l value for partial waves used in the computation of non-muffin-tin matrix elements (lnsmax=4 is quite good)

line 3:

free format
Etrial, ndiff, Napw

Etrial		default energy used for all $E_l$ to obtain $u_l(r,E_l)$ as regular solution of the radial Schrödinger equation [used in equ.2.4,2.7] (see figure 7.1).
ndiff		number of exceptions (specified in the next ndiff lines)
Napw		0 ... use LAPW basis, 1 ... use APW-basis for all ``global'' $l$ values of this atom. We recommend to use LAPW here.

line 4:

format(I2,2F10.5,A4)
l, El, de, switch, NAPWL

l		l of partial wave
El		$E_l$ for L=l
de		energy increment
		de=0: this E(l) overwrites the default energy (from line 3)
		de$\ne$ 0: a search for a resonance energy using this increment is done. The radial function $u_l(r,E)$ up to the muffin-tin radius RMT varies with the energy. A typical case is schematically shown in Fig. 7.1.
		At the bottom of the energy bands u has a zero slope (bonding state), but it has a zero value (antibonding state) at the top of the bands. One can search up and down in energy starting with $E_l$ using the increment de to find where $u_l(R_{MT},E)$ changes sign in value to determine $E_{top}$ and in slope to specify $E_{bottom}$. If both are found $E_l$ is taken as the arithmetic mean and replaces the trial energy. Otherwise $E_l$ keeps the specified value. For $E_{top}$ and $E_{bottom}$ bounds of +1 and -10 Ry are defined respectively, and if they are not found, they remain at the initial value set to -200.
switch		used only if de.ne.0
	CONT	calculation continues, even if either $E_{top}$ or $E_{bottom}$ are not found
	STOP	calculation stops if not both $E_{top}$ and $E_{bottom}$ are found (especially useful for semi-core states)
NAPWL		0 ... use LAPW basis, 1 ... use APW-basis for this $l$ value of this atom. We recommend to use APW+lo when the corresponding wavefunction is ``localized'' and thus difficult to converge with standard LAPW (like 3d functions) and/or when the respective atomic sphere is small compared to the other spheres in the unit cell.

Figure 7.1: Schematic dependence of DOS and $u_l(r,E_l)$ on the energy

$>>>$:line 4

is repeated ndiff times (see line 3) for each exception. If the same l value is specified twice, local orbitals are added to the (L)APW basis. The first energy ($E_1$) is used for the usual LAPW's and the second energy ($E_2$) for the LOs, which are formed according to (see equ. 2.7): $u_{E_1} + \dot u_{E_1} + u_{E_2}$.
Note: You may change the automatically created input and add d- or f-local orbitals to reduce the linearization error (e.g. in late transition metals you could put $E_{3d}$ at 0.0 and 1.0 Ry) or s, p, d, and/or f-LOs at very high energy (e.g. 2.0 - 3.0 Ry) to better describe unoccupied states. This input is also adapted during scf when -in1new is specified, while the original in1 file is saved in case.in1_orig.

$>>>$:lines 3 and 4

are repeated for each non equivalent atom

line 5:

format (20x,i1,2f10.1)
unit-number, Emin, Emax

unit-number		file number from which the k-vectors in the irreducible wedge of the Brillouin zone are read. 5 specifies the input file itself (as shown in the example), default is 4, for which the corresponding information is contained in case.klist (generated by KGEN).
EMIN, EMAX		energy window in which eigenvalues shall be searched (overrides setting in case.klist. A small window saves computer time, but it also limits the energy range for the DOS calculation of unoccupied states.

line 6:

format (A10,4I5,3F5.2)
name, ix,iy,iz, idv, weight

name		name of k-vector (optional)
		$>>>$: the last line must be END !!
ix,iy,iz, idv		defines the k-vector in units of $2\pi/a$, $2\pi/b$, $2\pi/c$, where x= ix/idv etc.
weight		of k-vector (order of group of k)

$>>>$: line 6

is repeated for each k-vector in the IBZ, but Emin and Emax may be omitted after the first k point. The utility program kgen (see section 6.5) provides a list of such vectors (on a tetrahedral mesh) in case.klist.

$>>>$: the last line

must be END

3. LAPWSO (adds spin orbit coupling)

lapwso includes spin-orbit (SO) coupling in a second-variational procedure and computes eigenvalues and eigenvectors (stored in case.vectorso) using the scalar-relativistic wavefunctions from lapw1. For reference see Singh 94 and Novák 97. The SO coupling must be small, as it is diagonalized in the space of the scalar relativistic eigenstates. For large spin orbit effects it might be necessary to include many more eigenstates from lapw1 by increasing EMAX in case.in1 (up to 10 Ry!). We also provide an additional basisfunction, namely an LO with a $p_{1/2}$ radial wavefunction, which improves the basis and removes to a large degree the dependency of the results on EMAX and RMT (see Kunes et al. 2001). SO is considered only within the atomic spheres and thus the results may depend to some extent on the choice of atomic spheres radii. The nonspherical potential is neglected when calculating $dV\over{dr}$.

In spin-polarized calculations the presence of spin-orbit coupling can split equivalent atoms into non-equivalent ones or at least may lower the symmetry of the system. It is then necessary to consider a larger part of the Brillouin zone and the input for lapw2 should also be modified since the potential has lower symmetry than in the non-relativistic case. The following inputs must be changed:

• case.struct
• case.klist
• case.kgen
• case.in2c
• case.in1

We recommend to use initso (see Sec.5.2.10) which helps you to setup spinorbit calculations, but in some cases an even more complicated procedure might be necessary.

Note: SO eigenvectors are complex and thus lapw2c must be used in a selfconsistent calculation.

1. Execution

The program lapwso is executed by invoking the command:

x lapwso [ -up -p -c ] or
lapwso lapwso.def

where here -up indicates a spin-polarized calculation (no ``-dn'' is needed, since spin-orbit will mix spin-up and dn states in one calculation).

2. Dimensioning parameters

The following parameters are used (collected in file param.inc):

FLMAX	constant = 3
LMAX	highest l of wave function inside sphere (consistent with lapw1)
LABC	highest l of wave function inside sphere where SO is considered
LOMAX	max l for local orbital basis
NRAD	number of radial mesh points
NLOAT	number of local orbitals

3. Input

A sample input for lapwso is given below. It will be generated automatically by initso

------------------ top of file: case.inso --------------------
WFFIL
 4  0  0                      llmax,ipr,kpot 
 -10.0000   1.5000           Emin, Emax
     0 0 1                       h,k,l (direction of magnetization)
   1                             number of RLO
 1    -3.5      0.005 STOP       E-param for RLO
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

format(A5)
switch

WFFIL

wavefunctions will also be calculated for scf-calculation. Otherwise only eigenvalues are calculated.

line 2:

free format
LLMAX, IPR, KPOT

LLMAX		maximum l for wavefunctions
IPR		print switch, larger numbers give additional output
KPOT	0	V(dn) potential is used for $<dn\vert V\vert dn>$ elements, V(up) for $\mbox{$<up\vert V\vert up>$}$ and [V(dn)+V(up)]/2 for $<up\vert V\vert dn>$.
	1	averaged potential used for all matrix elements.

line 3:

free format
Emin, Emax

Emin		minimum energy for which the output eigenvectors and eigenenergies will be printed (Ry)
Emax		maximum energy

line 4:

free format

h,k,l

vector describing the direction of magnetization

line 5:

free format

nlr

number of $p_{1/2}$ LOs for this atom

line 6:

free format
l, El, de, switch

l		l of relativistic LO
El		$E_l$ for L=l
de		energy increment (see lapw1)
switch		used only if de.ne.0
	CONT	calculation continues, even if either $E_{top}$ or $E_{bottom}$ are not found
	STOP	calculation stops if not both $E_{top}$ and $E_{bottom}$ are found (especially useful for semi-core states)

$>>>$: line 5 and 6

are repeated for each atom in the cell.

4. LAPW2 (generates valence charge density expansions)

lapw2 uses the file case.vector and computes the Fermi-energy and the expansions of the electronic charge densities in a representation according to eqn. 2.10 for each occupied state and each k-vector; then the corresponding (partial) charges inside the atomic spheres are obtained by integration. The partial charges for each state (energy eigenvalue) and each k-vector are written to files case.help31, case.help32 etc., where the last digit gives the atomic index of inequivalent atoms. In addition ``Pulay-corrections`` to the forces at the nuclei are calculated here. For systems without inversion symmetry you have to use the program lapw2c (in connection with lapw1c).

1. Execution

The program lapw2 is executed by invoking the command:

x lapw2 [-c -up|dn -p -so -qtl -fermi] or
lapw2 lapw2.def [proc#] or lapw2c lapw2.def [proc#]

where proc# is the i-th processor number in case of parallel execution (see Fig. 5.2). The -so switch sets -c automatically.

For each inequivalent atom a file case.helpXX is defined, starting with XX=31. For complex calculations case.in2c is used, for semi-core calculations lapw2s.def differs from the regular file lapw2.def only in few points, namely case.in2 must be replaced by case.in2s and case.clmval by case.clmsc. For a spin-polarized case see the fcc Ni test case in the WIEN2k package.

2. Dimensioning parameters

The following parameters are used (collected in file param.inc):

IBLCK	Blocking parameter (32-255) in fourier.frc, optimize for best performance
IBLOCK	Blocking parameter (32-255) in l2main.frc, optimize for best performance
LMAX2	highest l in wave function inside sphere (smaller than in lapw1, at present must be .le. 6)
LOMAX	max l for local orbital basis
NCOM	number of LM terms in density
NGAU	max. number of Gaunt numbers
NRAD	number of radial mesh points

3. Input

A sample input for lapw2 is listed below, it is generated automatically by the programs lstart and symmetry.

------------------ top of file: case.in2 --------------------
TOT                         (TOT,FOR,QTL,EFG)
-1.2       32.000   0.5  0.05  (EMIN, # of electrons,ESEPERMIN, ESEPER0 )
TETRA       0.0      (EF-method (ROOT,TEMP,GAUSS,TETRA,ALL),value)
 0 0 2 0 2 2 4 0 4 2 4 4 
 0 0 1 0 2 0 2 2 3 0 3 2 4 0 4 2 4 4
14.0      (GMAX)
FILE      (NOFILE, optional)
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

format(A5)
switch

TOT		total valence charge density expansion inside and outside spheres
FOR		same as TOT, but in addition a ``Pulay'' force contribution is calculated (this option costs extra computing time and thus should be performed only at the final scf cycles, see run_lapw script in sec. 5.1.3)
QTL		partial charges only (generates file case.qtl for DOS calculations)
EFG		computes decomposition of electric field gradient (EFG), contributions from inside spheres (the total EFG is computed in lapw0).
CLM		CLM coefficients only
FERMI		Fermi energy only, this produces weight files for parallel execution and for the optics package.
$>>>$:		TOT and FOR are the standard options, QTL is used for density of states (or energy bandstructure) calculations, EFG for analysis, while FOURI, CLM are for testing only.

line 2:

free format
emin, ne, esepermin, eseper0

emin		lower energy cut-off for defining the range of occupied states
ne		number of electrons (per unit cell) in that energy range
esepermin		LAPW2 tries to find the ``mean'' energies for each $l$ channel, for both the valence and the semicore states. To define ``valence'' and ``semicore'' it starts at (EF - ``esepermin'') and searches for a ``gap'' with a width of at least ``eseper0'' and defines this as seperation energy of valence and semicore
eseper0		minimum gap width (see above). The values esepermin and eseper0 will only influence results if the option -in1new is used

line 3:

format(a5,f10.5)
efmod, eval

efmod		determines how E$_F$ is determined
	ROOT	E$_F$ is calculated and k space integration is done by root sampling (this can be used for insulators, but for metals poor convergence is found)
	TEMP	E$_F$ is calculated where each eigenvalue is temperature broadened using a Fermi function with a broadening parameter of eval Ry. (e.g. eval=0.002 gives good total energy convergence, but has no ``physical`` justification).

	GAUSS	E$_F$ is calculated as above but a Gaussian smearing method is used with a width of eval Ry. (e.g. eval=0.002 gives good total energy convergence, but has no ``physical`` justification).
	TETRA	E$_F$ is calculated and k space integration is done by the modified (if eval is .lt. 100) or standard (eval .gt. 100) tetrahedron-method (Blöchl 94). This ``standard'' scheme is recommended for optic. In this case the file case.kgen, consistent with the k-mesh used in lapw1, must be provided (see Sec. 7.2). This is the recommended option although convergence may be slower than with Gauss- or temperature-smearing.
	ALL	All states up to eval are used. This can be used to generate charge densities in a specified energy interval.
eval		when efmod is set to TEMP or GAUSS, eval specifies the width of the broadening (in Ry), if efmod is set to ALL, eval specifies the upper limit of the energy window, if efmod is set to TETRA, eval .gt. 100 specifies the use of the standard tetrahedron method instead of the modified one (see above).

line 4:

format (20I2)

L,M

LM values of lattice harmonics expansion (equ. 2.10), defined according to the point symmetry of the corresponding atom; generated in SYMMETRY, MUST be consistent with the local rotation matrix defined in case.struct (details can be found in Kara and Kurki-Suonio 81). CAUTION: additional LM terms which do not belong to the lattice harmonics will in general not vanish and thus such terms must be omitted. Automatic termination of the $lm$ series occurs when a second 0,0 pair appears within the list. When you change the $l,m$ list during an SCF calculation the Broyden-Mixing is restarted in MIXER.

$>>>$line 4:

must be repeated for each inequivalent atom

Table 7.26: LM combinations of ``Cubic groups'' (3$\parallel$(111)) direction, requires ``positive atomic index'' in case.struct. Terms that should be combined (Kara and Kurki-Suonio 81) must follow one another.

Symmetry	LM combinations
23	0 0, 4 0, 4 4, 6 0, 6 4,-3 2, 6 2, 6 6,-7 2,-7 6, 8 0, 8 4, 8 8,-9 2,-9 6,-9 4,-9 8,10 0, 10 4,10 8, 10 2, 10 6, 10 10
M3	0 0, 4 0, 4 4, 6 0, 6 4, 6 2, 6 6, 8 0, 8 4, 8 8,10 0, 10 4,10 8, 10 2, 10 6, 10 10
432	0 0, 4 0, 4 4, 6 0, 6 4, 8 0, 8 4, 8 8,-9 4,-9 8,10 0, 10 4,10 8
-43M	0 0, 4 0, 4 4, 6 0, 6 4,-3 2,-7 2,-7 6, 8 0, 8 4, 8 8,-9 2,-9 6,10 0, 10 4,10 8
M3M	0 0, 4 0, 4 4, 6 0, 6 4, 8 0, 8 4, 8 8,10 0, 10 4,10 8

Table 7.27: LM combination and local coordinate system of ``non-cubic groups'' (requires ``negative atomic index'' in case.struct)

Symmetry	Coordinate axes	Indices of $Y_{\pm LM}$	crystal system
1	any	ALL ($\pm$l,m)	triclinic
-1	any	($\pm$2l,m)
2	2$\parallel$ z	($\pm$l,2m)	monoclinic
M	m$\perp$z	($\pm$l,l-2m)
2/M	2$\parallel$z, m$\perp$z	($\pm$2l,2m)
222	2$\parallel$z, 2$\parallel$y, (2$\parallel$x)	(+2l,2m), (-2l+1,2m)	orthorhombic
MM2	2$\parallel$z, m$\perp$y, (2$\perp$x)	(+l,2m)
MMM	2$\perp$z, m$\perp$y, 2$\perp$x	(+2l,2m)
4	4$\parallel$z	($\pm$l,4m)	tetragonal
-4	-4$\parallel$z	($\pm$2l,4m), ($\pm$2l+1,4m+2)
4/M	4$\parallel$z, m$\perp$z	($\pm$2l,4m)
422	4$\parallel$z, 2$\parallel$y, (2$\parallel$x)	(+2l,4m), (-2l+1,4m)
4MM	4$\parallel$z, m$\perp$y, (2$\perp$x)	(+l,4m)
-42M	-4$\parallel$z, 2$\parallel$x (m=xy$\rightarrow$yx)	(+2l,4m), (-2l+1,4m+2)
4MMM	4$\parallel$z, m$\perp$z, m$\perp$x	(+2l,4m)
3	3$\parallel$z	($\pm$l,3m)	rhombohedral
-3	-3$\parallel$z	($\pm$2l,3m)
32	3$\parallel$z, 2$\parallel$y	(+2l,3m), (-2l+1,3m)
3M	3$\parallel$z, m$\perp$y	(+l,3m)
-3M	-3$\parallel$z, m$\perp$y	(+2l,3m)
6	6$\parallel$z	($\pm$l,6m)	hexagonal
-6	-6$\parallel$z	(+2l,6m), ($\pm$2l+1,6m+3)
6/M	6$\parallel$z, m$\perp$z	($\pm$2l,6m)
622	6$\parallel$zm, 2$\parallel$y, (2$\parallel$x)	(+2l,6m), (-2l+1,6m)
6MM	6$\parallel$z, m$\parallel$y, (m$\perp$x)	(+l,6m)
-62M	-6$\parallel$z, m$\perp$y, (2$\parallel$x)	(+2l,6m), (+2l+1,6m+3)
6MMM	6$\parallel$z, m$\perp$z, m$\perp$y, (m$\perp$x)	(+2l,6m)

line 5:

free format

GMAX

max. G (magnitude of largest vector) in charge density Fourier expansion. For systems with short H bonds larger values (e.g. GMAX up to 25) could be necessary. Calculations using GGA (potential option 13 or 14 in case.in0) should also employ a larger GMAX value (e.g. 14), since the gradients are calculated numerically on a mesh determined by GMAX. When you change GMAX during an scf calculation the Broyden-Mixing is restarted in mixer.

line 6:

A4

reclist

FILE

writes list of K-vectors into file case.recprlist or reuses this list if the file exists. The saved list will be recalculated whenever GMAX, or a lattice parameter has been changed.

NOFILE

always calculate new list of K-vectors

5. SUMPARA (summation of files from parallel execution)

sumpara is a small program which (in parallel execution of WIEN2k) sums up the densities (case.clmval_*) and quantities from the case.scf2_* files of the different parallel runs.

1. Execution

The program sumpara is executed by invoking the 2 commands as follows:

x sumpara -d [-updn] and then
sumpara sumpara.def #_of_proc

where #_of_proc is the numbers of parallel processors used. It is usually called automatically from lapw2para or x lapw2 -p.

2. Dimensioning parameters

The following parameters are listend in file param.inc, but usually they need not to be modified:

NCOM	number of LM terms in density
NRAD	number of radial mesh points
NSYM	order of point group

6. LAPWDM (calculate density matrix)

This program was contributed by:

lapwdm calculates the density matrix needed for the orbital dependent potentials generated in orb. Optionally it also provides orbital moments, orbital and dipolar contributions to the hyperfine field (only for the specified atoms AND orbitals). It calculates the average value of the operator X which behaves in the same way as the spin-orbit coupling operator: it must be nonzero only within the atomic spheres and can be written as a product of two operators - radial and angular:

$X =X_r(r)*X_{ls}(\vec{l},\vec{s})$

$X_r(r)$ and $X_{ls}(\vec{l},\vec{s})$ are determined by RINDEX and LSINDEX in the input as described below:

• RINDEX=1 LSINDEX=1: $<$X$>$ is number of electrons inside the atomic sphere (for test)
• RINDEX=2 LSINDEX=1: $<$X$>$ is the $<1/r^3>$ expectation value inside the atomic sphere
• RINDEX=1 LSINDEX=2: $<$X$>$ is the projection of the spin moment inside the atomic sphere
• RINDEX=1 LSINDEX=3: $<$X$>$ is the projection of the orbital moment inside the atomic sphere
• RINDEX=3 LSINDEX=3: $<$X$>$ is the orbital part of the hyperfine field at the nucleus
• RINDEX=3 LSINDEX=5: $<$X$>$ is the spin dipolar part of the hyperfine field at the nucleus

To use the different operators, set the appropriate input. More information and extentions to operators of similar behavior may be obtained directly from P. Novák (97)

1. Execution

The program lapwdm is executed by invoking the command:

x lapwdm [ -up/dn -p -c -so ] or
lapwdm lapwdm.def

2. Dimensioning parameters

The following parameters are used (collected in file param.inc):

FLMAX	constant = 3
LMAX	highest l of wave function inside sphere (consistent with lapw1)
LABC	highest l of wave function inside sphere where SO is considered
LOMAX	max l for local orbital basis
NRAD	number of radial mesh points

3. Input

A sample input for lapwdm is given below.

------------------ top of file: case.indm --------------------
-9.                      Emin cutoff energy
 1                       number of atoms for which density matrix is calculated
 1  1  2      index of 1st atom, number of L's, L1
 0 0          r-index, (l,s)-index
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

free format

emin

lower energy cutoff (usually set to very low number).

line 2:

free format

natom

number of atoms for which the density matrix is calculated

line 3:

free format
iatom, nl, l

	iatom	index of atom for which the density matrix should be calculated
	nl	number of l-values for which the density matrix should be calculated
	l	l-values for which the density matrix should be calculated

line 3 is repeated natom times

t

line 4:

free format, optional
RINDEX, LSINDEX

	RINDEX	0-3, as described in the introduction to lapwdm
	LSINDEX	0-5, as described in the introduction to lapwdm

7. ORB (Calculate orbital dependent potentials)

This program was contributed by:

orb calculates the orbital dependent potentials, i.e. potentials which are nonzero in the atomic spheres only and depend on the orbital state numbers $l, m$. In the present version the potential is assumed to be independent of the radius vector und needs the density matrix calculated in lapwdm. Three different potentials are implemented in this package:

• LDA+U. There are three variants of this method, two of them are discussed in Novák et al. 2001
  1. LDA+U(SIC) - introduced by Anisimov et al. 1993, with an approximate correction for the self-interaction correction.
  2. LDA+U(AMF) - introduced by Czyzyk and Sawatzky 1994 as 'Around the Mean Field' method. (In Novák et al. 2001 it is denoted as LDA+U(DFT)). This version is (probably) more suitable for full potential methods.
  3. LDA+U(HMF) - in addition the Hubbard model in the mean field approximation, as introduced by Anisimov et al. 1991 is also implemented. Note, however, that it is to be used with the LDA (not LSDA) exchange-correlation potential in spin polarized calculations!
  All variants are implemented in the rotationally invariant way (Liechtenstein et al. 1995). If LDA+U is used in an unrestricted, general way, it introduces an orbital field in the calculation (in analogy to the exchange field in spin-polarized calculations, but it interacts with the orbital, instead of spin momentum). The presence of the orbital field lowers the symmetry. In particular the complex version of LAPW1 must be used. Care is needed when dealing with the LDA+U orbital field. It may be quite large, and without specifying its direction it may fluctuate, leading to oscilations of scf procedure or/and to false solutions. It is therefore necessary to use it in combination with the spin-orbit coupling, preferably running first LSDA+(s-o) and then slowly switching on the LDA+U orbital field. If the LDA+U orbital polarization is not needed, it is sufficient to run real version of LAPW1, which then automatically puts the orbital field equal to zero. For systems without the centre of inversion, when LAPW1 must be complex, an extra averaging of the LDA+U potential is necessary.
• Orbital polarization. The additional potential has the form (Brooks 1985, Eriksson et al. 1989):
  

  $$V_{OP}= c_{OP}<L_z>l_z$$ (11)

  
  
  where $c_{OP}$ is the orbital polarization parameter, $<L_z>$ is projection of the orbital momentum on the magnetization direction and $l_z$ is single electron orbital momentum component $z$ parallel to $\vec{M}$.
• Interaction with the external magnetic field. In this case the additional potential has a simple form:
  

  $$V_{Bext}= \mu_B \vec{B}_{ext} (\vec{l} + 2\vec{s}) .$$ (12)

  
  
  The interaction with the electronic spin is taken into account by shifting the spin up and spin down exchange correlation potentials in LAPW0 by the energy $+\mu_B B_{ext} -\mu_B B_{ext}$, respectively. The interaction of $B_{ext}$ with spin could be as well calculated using the 'Fixed spin moment' method. For an interaction with the orbital momentum it is necessary to specify the atoms and angular momentum numbers for which this interaction will be considered. Caution is needed when considering interaction of the orbital momentum with $B_{ext}$ in metallic or metallic-like systems. For the analysis see the paper by Hirst 1997

In all cases the resulting potential for a given atom and orbital number $l$ is a Hermitean, $(2l+1)$x$(2l+1)$ matrix. In general this matrix is complex, but in special cases it may be real.

1. Execution

The program orb is executed by invoking the command:

x orb [ -up/dn ] or orb up/dnorb.def

2. Dimensioning parameters

The following parameters are used (collected in file param.inc):

LABC	highest l+1 value of orbital dependent potentials
NRAD	number of radial mesh points

3. Input

Since this program can handle three different cases, examples and descriptions for all cases are given below:

1. Input for all potentials

line 1:

free format
nmod,natorb,ipr

	nmod	defines the type of potential 1...LDA+U, 2...OP, 3...$B_{ext}$
	natorb	number of atoms for which orbital potential $V_{orb}$ is calculated
	ipr	printing option, the larger ipr, the longer the output

line 2:

(A5,f8.2)
mixmod,amix

	mixmod	BROYD or PRATT (see MIXER for more information)
	amix	coefficient for the Pratt mixing (and first BROYD iteration) of $V_{orb}$

line 3:

free format
iatom(i),nlorb(i),(lorb(li,i),li=1,nlorb(i))

	iatom	index of atom in struct file
	nlorb	number of orbital moments for which Vorb shall be applied
	lorb	orbital numbers (repeated nlorb-times)

3rd line repeated natorb-times

2. Input for LDA+U (nmod=1)

line 4:

free format

nsic		defines 'double counting correction'
	nsic=0	'AMF method' (Czyzyk et al. 1994)
	nsic=1	'SIC method' (Anisimov et al. 1993, Liechtenstein et al. 1995)
	nsic=2	'HMF method' (Anisimov et al. 1991)

line 5:

free format

U(li,i), J(li,i)

Coulomb and exchange parameters, U and J, for LDA+U in Ry for atom type i and orbital number li

5th line repeated natorb-times, for each natorb repeated nlorb-times

Example of the input file for NiO (LDA+U included for two inequivalent Ni atoms that have indexes 1 and 2 in the structure file):

1  2  0                        nmod, natorb, ipr
BROYD 0.3                      mixmod, amix
1 1 2                          iatom nlorb, lorb
2 1 2                          iatom nlorb, lorb
1                              nsic (LDA+U(SIC) used)
0.59 0.07                      U J
0.59 0.07                      U J

3. Input for Orbital Polarization (nmod=2)

line 4:

(free format)

nmodop		defines mode of 'OP'
	1	average $L_z$ taken separatelly for spin up, spin down
	0	average $L_z$ is the sum for spin up and spin down

line 5:

(free format)

Ncalc(i)
	1	Orb.pol. parameters are calculated ab-initio
	0	Orb.pol. parameters are read from input

this line is repeated natorb-times

line 6:

(free format) (only if Ncalc=0, then repeated nlorb-times)

pop(li,i)

OP parameter in Ry

line 7:

(free format)
xms(1), xms(2), xms(3)

direction of magnetization expressed in terms of lattice vectors

Example of the input file for NiO (total $<L_z>$ used in (1), OP parameters calculated ab-initio, $\vec{M}$ along [001]):

2  2  0                        nmod, natorb, ipr
BROYD  0.2                     mixmod, amix
1 1 2                          iatom nlorb, lorb
2 1 2                          iatom nlorb, lorb
0                              nmodop
1                              Ncalc
0. 0. 1.                       direction of M in terms of lattice vectors

4. Input for interaction with $B_{ext}$ (nmod=3)

line 4:

(free format)

$B_{ext}$

external field in Tesla

line 5:

(free format)
xms(1), xms(2), xms(3)

direction of magnetization expressed in terms of lattice vectors

Example of the input file for NiO, ($B_{ext}$= 4 T, along [001]):

3  2  0                        nmod, natorb, ipr
BROYD, 0.3                     mixmod, amix
1 1 2                          iatom nlorb, lorb
2 1 2                          iatom nlorb, lorb
4.                             Bext in T
0. 0. 1.                       direction of Bext in terms of lattice vectors

8. LCORE (generates core states)

lcore is a modified version of the Desclaux (69, 75) relativistic LSDA atomic code. It computes the core states (relativistically including SO, or non-relativistically if ``NREL'' is set in case.struct) for the current spherical part of the potential (case.vsp). It yields core eigenvalues, the file case.clmcor with the corresponding core densities, and the core contribution to the atomic forces.

1. Execution

The program lcore is executed by invoking the command:

lcore lcore.def or x lcore [-up|-dn]

For the spin-polarized case see fcc Ni on the distribution tape.

2. Dimensioning parameters

The following parameter is listend in file param.inc:

NRAD	number of radial mesh points

3. Input

Below is a sample input (written automatically by lstart)

for $TiO_2$ (rutile), one of the test cases provided with the WIEN2k

package.

------------------ top of file: case.inc --------------------
 4 0.0       # of orbitals,  shift of potential
 1,-1,2      n (principal quantum number), kappa, occup. number
 2,-1,2      2s
 2,-2,4      2p
 2, 1,2      2p*
 1   0.0          # of orbital of second atom
 1,-1,2      1s
 0           end switch
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

free format
nrorb, shift

nrorb		number of core orbitals
shift		shift of potential for ``positive'' eigenvalues (e.g. for 4f states as core states in lanthanides)

line 2:

free format
n, kappa, occup

n		principle quantum number
kappa		relativistic quantum number (see Table 6.4)
occup		occupation number (including spin), fractial occupations supported

$>>>$: line 2

is repeated for each orbital (nrorb times; see line 1)

$>>>$: line 1 and 2

are repeated for each inequivalent atom. Atoms without core states (e.g. H or Li) must still include a 1s orbital, but with occupation zero.

line 3:

free format

0

zero indicating end of job

9. MIXER (adding and mixing of charge densities)

In mixer the electron densities of core, semi-core, and valence states are added to yield the total new (output) density (in some calculations only one or two types will exist). Proper normalization of the densities is checked and enforced (by adding a constant charge density in the interstitial). As it is well known, simply taking the new densities leads to instabilities in the iterative SCF process. Therefore it is necessary to stabilize the SCF cycle. In WIEN2k this is done by mixing the output density with the (old) input density to obtain the new density to be used in the next iteration. Two mixing schemes are implemented:

1. straight mixing as originally proposed by Pratt (52) with a mixing factor Q
   
   $$\rho_{new}(r)=(1-Q)\rho_{old}(r)+Q\rho_{output}(r)$$
   
   

2. the Broyden-II mixing scheme (Singh et al., 86), in which all the expansion coefficients of the density from several preceding iterations are utilized to calculate an optimal mixing fraction for each coefficient in each iteration.

At the outset of a new calculation (for any changed computational parameter such as k-mesh, matrix size, lattice constant etc.), any existing case.broydX files must be deleted (since the iterative history which they contain refers to a ``different`` incompatible calculation). In addition, in some cases better convergence can be achieved, if these files are removed every 15-20 iterations. Usually the Broyden scheme is much better than Pratt's scheme and thus is recommended.

After modifications to the case.struct file (lattice parameters, atomic positions) a run with mixing factor 0.0 can be used to renormalize the case.clmsum_old file from the previous case. If the file case.clmsum_old can not be found by mixer, a ``PRATT-mixing`` with mixing factor 1.0 is done.

Note: a case.clmval file must always be present, since the LM values and the K-vectors are read from this file.

The total energy and the atomic forces are computed in mixer by reading the case.scf file and adding the various contributions computed in preceding steps of the last iteration. Therefore case.scf must not contain a certain ``iteration-number'' more than once and the number of iterations in the scf file must not be greater than 99.

1. Execution

The program mixer is executed by invoking the command:

mixer mixer.def or x mixer

A spin-polarized case will be detected automatically by x due to the presence of a case.clmvalup file. For an example see fccNi (sec. 10.2) in the WIEN2k package.

2. Dimensioning parameters

The following parameters are collected in file param.inc, :

NCOM	number of LM terms in density
NRAD	number of radial mesh points
NSYM	order of point group

3. Input

Below a sample input (written automatically by lstart) is provided for $TiO_2$ (rutile), one of the test cases provided with the WIEN2k package.

------------------ top of file: case.inm --------------------
BROYD 0.d0  YES  (PRATT/BROYD, background charge (+1 for additional e), NORM
  0.4            MIXING FACTOR
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

(A5,f8.2,A5)
switch, bgch, norm

switch	BROYD	Broyden's scheme
	PRATT	Pratt's scheme

bgch

Background charge for charged cells (+1 for additional electron, -1 for core hole, if not neutralized by additional valence electron)

norm	YES	Charge densities are normalized to sum of Z
	NO	Charge densities are not normalized

line 2:

free format

factor

mixing parameter Q (for Pratt and in the first iteration of Broyden). In the first iteration using Broyden's scheme: Q is automatically reduced by the program depending on the average charge distance :DIS and the number of non-equivalent TM (f)-elements. In case that the scf cycle fails due to large charge fluctuations, this factor must be further reduced (sometimes by an order of magnitude) before the calculations can be restarted (see sect.12