8 Programs for analysis, calculation of properties, and geometry optimization

1 TETRA (density of states)

This program calculates total and partial density of states (DOS) by means of the modified tetrahedron method (Blöchl et al 1994). Please note, the tetrahedron method will not work with just one k-point and tetra will automatically switch to a Gaussian broadening scheme (with default broadening of 0.01 Ry). It can also be selected by input (see below), but is not recommended for small unit cells.

It uses the partial charges in case.qtl generated by lapw2 (switch QTL) and generates the DOS in states/Ry (files case.dos1/2/3/...) and in states/eV (with respect to the Fermi energy; files case.dos1/2/3ev). In spin-polarized calculations the DOS is given in states/Ry/spin (or states/eV/spin).

Please note: The total DOS is equal to the sum over the atoms of the total-atomic DOS (inside spheres) and the interstitial-DOS. (Thus in the total-atomic DOS the ``multiplicity'' of an atom is considered). On the other hand, in the partial (lm-like) DOS the multiplicity is not considered and one obtains the total-atomic DOS as a sum over all partial DOS times the multiplicity.

The ``m-decomposed'' DOS (e.g. $p_z, p_y, p_x$) is given with respect to the local coordinate system for each atom as defined by the local rotation matrix (see Appendix A), unless you have used x qtl to generate the case.qtl and specified a specific coordinate system in case.inq (see Chapter 8.2).

Using the switches -rxes E1 E2 it is possible to generate a ``weight-file'', where each k-point is weighted according to its contribution to the DOS in the energy range E1-E2. This weight-file case.rxes can be used using the switch -rexs to calculate the DOS with these weights. This option might be useful to simulate the E-dependency of RXES spectra, or in general calculate a ``DOS'' of regions around selected k-points only.

The density of states in files case.dos1/2/3/... or case.dos1/2/3/...ev can be plotted by dosplot2_lapw (see 5.7.4).

It is strongly recommended that you use ``Run Programs $\rightarrow$ Tasks $\rightarrow$ Density of States'' from w2web.

1 Execution

The program tetra is executed by invoking the command:

tetra tetra.def or x tetra [-up|dn -rxes -rxesw E1 E2]

2 Dimensioning parameters

The following parameters are listed in file param.inc:

MG	max. number of DOS cases
LXDOS	usually 1, except for ``cross-DOS'' when using TELNES.2 = 3 (not needed anymore for TELNES3)

3 Input

The required input file case.int can optionally be created using the w2webinterface or the configure_int_lapw script (see 5.2.7).

An example is given below:

------------------ top of file: case.int ------------------
TiO2                         		# Title  
 -1.000     0.00250   1.200 0.003 	# EMIN, DE, EMAX for DOS, GAUSS-Broad
    7    N   0.000                   	# NUMBER OF DOS-CASES, G/L/B broadening
    0    1   tot             		# jatom, doscase, description
    1    2   Ti-s   
    1    3   Ti-p   
    1    4   Ti-px
    1    5   Ti-py
    1    6   Ti-pz
    2    1   O-tot
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

free format
title

line 2:

free format
emin, delta, emax, broad

emin, delta, emax		specifies the energy mesh (in Ry) where the DOS is calculated. (emin should be set slightly below the lowest valence band; emax will be checked against the lowest energy of the highest band in case.qtl, and set to the minimum of these two values; delta is the energy increment.
broad		Gauss-broadening factor. Must be greater than delta to have any effect.

line 3:

free format
ndos, Bmethod, broadening

ndos		specifies the number of DOS cases to be calculated. It should be at least 1. The corresponding output is written in groups of 7 to respective case.dosX files
Bmethod		optional input (can be omitted) to select instead of the tetrahedron method:
	G	Gaussian broadening
	L	Lorentzian broadening
	B	both, Gauss and Lorentzian broadening
broadening		parameters in Ry, typically below/around 0.01 (optional, specify two numbers for B)

line 4:

(2i5,3x,a6)
jatom, jcol, description

jatom		specifies for which atom the DOS is calculated. 0 means total DOS, $jatom=nat+1$ means DOS in the interstitial, where $nat$ is the number of inequivalent atoms. When spin-orbit is included, $jatom=nat+1$ gives total spin-up/dn DOS in a spinpolarized SO calculation, but is meaningless in a non-spinpolarized SO case.
jcol		specifies the column to be used in the respective QTL-file. 1 means total, 2 ...s, 3 ...p, ...The further assignment depends on the value of ISPLIT set in case.struct (see sec. 4.3); the respective description can be found in the header of case.qtl.
description		text used for further identification.

$»>$:line 4

is repeated ``ndos`` times

2 QTL (calculates special partial charges and population matrices)

This program was contributed by:

qtl creates the input for calculating total and projected density of states of selected atoms and selected $l$-subshells. It thus provides similar data as lapw2 -qtl, but it allows for additional options. In particular it supports calculation of DOS projected on relativistic states , , , DOS projected on states in a rotated coordinate system and DOS projected on individual $f$ states. qtl also allows to calculate population matrix and energy resolved population matrix. Comparing to lapwdm population matrix, the matrix created by qtl may contain also the cross terms between different orbital and spin numbers and it can be energy resolved. Important option of the qtl is the symmetrization that makes the calculation longer, but must be switched on whenever the quantities, which are not invariant are calculated. Detailed description may be found in QTL - technical report by P. Novák. The calculation is based on the spectral decomposition of a density matrix on a given atomic site and its transformation to the required basis.

The output is written to case.qtl [up/dn]. For the DOS calculation the file case.qtltext [up/dn] is created in which the ordering of partial charges is given. Please note, that in contrast to case.qtl [up/dn] from x lapw2 -qtl the total partial charge of an atom is NOT multiplied with its ``multiplicity'' and contains only the sum of the requested l,m terms (eg. s,p,d) and thus not all contributions. Also the interstital charge will usually be NOT correct.

1 Execution

The program qtl is executed by invoking the command:

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

2 Input

A sample input (a default is created automatically during init_lapw for case.inq is given below.

------------------ top of file: case.inq --------------------
-7.  2.          Emin Emax
 2               number of selected atoms
 1  2  0  0      iatom1  qsplit1 symmetrize loro
 2  1  2         nL1 p  d
 3  3  1  1      iatom2  qsplit2 symmetrize loro
 4  0  1  2  3   nL2 s  p  d  f
 1. 1. 1.        new axis z
 ------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1: free format
emin,emax	energy window
line 2: free format
natom	number of atoms selected for calculation
line 3: free format
iatom, QSPLIT, symmetrize, loro
iatom	integer, index of atom
QSPLIT	integer, analog of ISPLIT in case.struct: see below
symmetrize	integer, =0 (no symmetrization), 1 (symmetrization)
loro	integer =0 original coord. system preserved
	=1 (new z axis)
	=2 (new z and x axes)
line 4: free format
Nl(iatom), (l(iatom,i),i=1,Nl(iatom))
Nl	number of orbital numbers selected for calculation
l	orbital numbers selected for calculation for atom iatom
line 5: free format
hz, kz, lz	real*8, direction of new axis z (if loro=1,2)

Lines starting from line 3 are repeated for each selected atom. Line 5 only appears when calculation in new coordinate system is required (loro 0). Axis z in this system is along hz,kz,lz (in units of the lattice vectors, need not be normalized). If not only the z axis, but also the x axis need to be specified, then loro must be equal to 2 and additional line

hx, kx, lx (real*8)

giving the direction of the new axis x, perpendicular to the new axis z must appear. For relativistic splitting (QSPLIT=0,-1) this rotation is ignored and z points along the direction of magnetization as defined in case.inso.

Indices of selected atoms, as well as the orbital numbers, must form an ascending sequence.

Table 8.5: Possible values of QSPLIT and their interpretation

QSPLIT	meaning
-2	DOS in basis according to ISPLIT from case.struct
-1	DOS in relativistic basis
0	DOS in relativistic basis, summed over
1	DOS in basis (no symmetry)
2	DOS in basis of real orbitals (no symmetry)
3	axial symmetry
4	hexagonal symmetry
5	cubic symmetry
6	user written unitary transformation
88	population matrix, no crossterms corresponds to ISPLIT=88
99	full population matrix including crossterms (as ISPLIT=99)

For QSPLIT=6 (unitary transformation prepared by user) the unitary matrices are read as in WIEN2k_07 qtl: For the i-th atom selected for qtl calculation, they are stored in case.cf$i and ordered according to increasing $l$. The unitary transformation matrix must rotate from the standard ${}_{lms}$-basis to the desired one. A few examples (e.g. ${}_{jj_z}$, ${}_{lms}$, or $e_g-t_{2g}$) are supplied with the code in $WIENROOT/SRC_templates/case.cf_* and must be copied to case.cf$i . For less common cases these must be generated by hand.

3 Output

The results in file case.qtl[up/dn] are written in the same format as lapw2 file case.qtl[up/dn] and thus they may be directly used by tetra.

The data for the interstital DOS correspond to ($nat$ is number of atom types). The ordering of densities for all selected atoms is summarized in the file case.qtltext[up/dn]. The qtltext file that corresponds to the input data given above is:

Ordering of DOS in QTL file for: HoMnO3 (Munoz)                                              

 atom   1 ordering of projected DOS 
p,px,py,pz, real basis
d,dz2,d(x2-y2),dxy,dxz,dyz, real basis
 
 atom   3 ordering of projected DOS 
 s                                                                
p,pxy,pz, axial basis                                        
d,dz2,d(x2-y2),d(yz+xz),dxy,  axial basis                              
f,A2,[x(T1)+y(T1)],z(T1),[ksi(T2)+eta(T2)],zeta(T2),   axial basis 
      A2=xyz  x(T1)=x(x2-3r2/5)  y(T1)=y(y2-3r2/5)  z(T1)=z(z2-3r2/5)
             ksi(T2)=x(y2-z2)   eta(T2)=y(z2-y2)   zeta(T2)=z(x2-y2)

 Data for interstital DOS correspond to atom index    8

The output for the population matrix integrated over energy is written to case.dmat [up/dn] that has the same format as analogous file calculated by lapwdm.

3 SPAGHETTI (energy bandstructure plots)

This program generates an energy bandstructure plot (postscript file case.spaghetti_ps and xmgrace file case.bands.agr) using the eigenvalues printed in case.output1 or case.outputso (with switch -so). Using the SCF potentials one runs x lapw1 -band with a special k-mesh (case.klist_band) along some high-symmetry lines (some sample inputs can be found in SRC_templates/*.klist or you create your own k-mesh using Xcrysden). As an option, one can emphasize the character of the bands by additionally supplying corresponding partial charges (file case.qtl which can be obtained using x lapw2 -qtl -band , see 7.5). This will be called ``band-character plotting`` below, in which each energy is drawn by a circle whose radius is proportional to the specified character of that state. It allows to analyze the character of bands (see also figures 3.12 and 3.13).

The file case.bands.agr can be opened directly with xmgrace. Within xmgrace, all features of the plot, such as the plot range, the plot size, line properties (style, thickness and color), axis properties, labels, etc.  can easily be changed by either using the menu (submenus of the "Plot" menu) or double-klicking on the corresponding part of the figure. The size of the characters for a ``band-character plot`` can be changed in the menu "Plot/Graph appearance/Z normalization". The figures can directly be printed or exported in eps, jpg, png and other formats, via the menus "File/Print setup" and "File/Print".

C.Persson has modified this program and it allows now also to draw connected lines. For this purpose it uses the irreducible representations (from file case.irrep produced by program irrep together with a table of ``compatibility relations'' to decide which points should be connected (non-crossing rule !). (Note: This option will NOT work on the surface of the BZ for non-symmorphic spacegroups, because the corresponding group-theory has not been implemented.)

The presence of ``incompatible'' case.irrep or case.qtl files (from a previous run or qtls from a DOS calculation) may crash spaghetti. In such cases it is necessary to remove these files explicitly.

It is strongly recommended that you use ``Run Programs $\rightarrow$ Tasks $\rightarrow$ Bandstructure'' from w2web.

1 Execution

The program spaghetti is executed by invoking the command:

spaghetti spaghetti.def or x spaghetti [-up|dn] [-so] [-p]

The -p switch directs spaghetti to use the case.output1_* files of a k-point parallel lapw1.

2 Input

An example is given below:

----------------- top of file: case.insp -------------------
### Figure configuration
 5.0   3.0                         # paper offset of plot
10.0  15.0                         # xsize,ysize [cm]
 1.0   4                           # major ticks, minor ticks
 1.0   1                           # character height, font switch
 1.1   2    4                      # line width, line switch, color switch
### Data configuration            
-25.0 15.0  2                      # energy range, energy switch (1:Ry, 2:eV)
1      0.74250                     # Fermi switch,  Fermi-level (in Ry units)
1   999                            # number of bands for heavier plotting   1,1
0      1    0.02                   # jatom, jtype, size  of heavier plotting   
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

free format
test

test

line must start with '###'. Begin of figure description. This tests also if you use the new input (different from WIEN97 or early WIEN2k versions)

line 2:

free format
xoffset, yoffset

xoffset		x offset (in cm) of origin of plot
yoffset		y offset (in cm) of origin of plot

line 3:

free format
xsize,ysize

xsize		plotsize in x direction (cm)
ysize		plotsize in y direction (cm)

line 4:

free format
eincr, mtick

eincr		energy increment where y-axis labels are printed (major ticks)
mtick		number of minor ticks of y-axis

line 5:

free format
charh, font

charh		scaling factor for size of labels
font	0	no text
	1	Times and Symbol
	2	Times,Times-Italic and Symbol
	3	Helvetica, Symbol, and Helvetica-Italic
	4	include your own fonts in defins.f

line 6:

free format
linew, ilin, icol

linew		line width
ilin	0	dots or open circles
	1	lines
	2	lines and open circles
	3	lines and filled circles
icol	0	black
	1	one-color plot
	2	three-color plot
	3	multi-color plot
	4	multi-color plot,one color for each irred. representation

line 7:

free format
test

test

line must start with '###'. Begin of data description.

line 8:

free format
emin, emax, iunits

emin		energy minimum of plot
emax		energy maximum of plot
iunits
	1	energies in Ry (internal scale)
	2	energies in eV with respect to $E_f$

line 9:

free format
iferm, efermi

iferm	0	no line at EF
	1	solid line at EF
	2	dashed line at EF
	3	dotted line at EF
efermi		Fermi energy (Ry); can be found in the respective case.scf file. If set to 999., $E_f$ is not plotted (and iunits=2 cannot be used)

line 10:

free format

nband1, nband2

lower and upper band index for bands which should show ``band-character plotting`` (if case.qtl is present and the proper switch is set, see below). In addition the corresponding x and y coordinates are written to file case.spaghetti_ene (which can be used for plotting with an external xy-plotting program).

line 11:

free format
jatom, jcol, jsize

jatom		If a case.qtl file is present, jatom indicates the atom whose character (selected by jcol) is used for ``band-character plotting`` (dots are replaced by circles with radii proportional to the corresponding weight, requires ilin=0,2,3). If set to zero or if case.qtl is not present, ``band-character plotting`` does not occur.
jcol		specifies the column to be used in the respective QTL-file. 1 means total, 2 ...s, 3 ...p, ...The further assignment depends on the value of ISPLIT set in case.struct. (ignored for jatom=0). The description can be found in the header of case.qtl.
jsize		size factor for radii of circles used in ``band-character plotting''

if line 11

is repeated, one can average the QTLs for different atoms (but with identical jcol and jsize).

4 IRREP (Determine irreducible representations)

This program was contributed by:

This program determines the irreducible representation for each eigenvalue and all your k-points. It is in particular usefull to analyse energy bands and their connectivity.

You need a valid vector file, but no other input is required. The output can be found in case.outputir and case.irrep. For nonmagnetic SO calculations you must set IPR=1 in case.inso.

The output of this program is needed when you want to draw bandstructures with connected lines (instead of ``dots'').

It will not work in cases of non-symmorphic spacegroups AND k-points at the surface of the BZ. See also $WIENROOT/SRC_irrep/README.

1 Execution

The program irrep is executed by invoking the command:

irrep [up/dn]irrep.def or x irrep [-so -up/dn ]

2 Dimensioning parameters

The following parameters are listend in file param.inc:

LOMAX	max. no. of local orbital. should be consistent with lapw1 and lapwso
NLOAT	number of different types of LOs
MSTP	max. step to describe k as a fraction
MAXDG	max. no. of degenerate eigenfunctions
MAXIRDG	max. no. of degenerate irr. representations
FLMAX	size of flag (FL) array (should be 4)
MAXIR	max. no. of irreducible representations
NSYM	max. no. of symmetry operations
TOLDG	min. energy deviation of degenerate states, in units of Rydberg

5 LAPW3 (X-ray structure factors)

This program calculates X-ray structure factors from the charge density by Fourier transformation.

You have to specify interactively valence or total charge density (because of the different normalization of case.clmsum and case.clmval) and a maximum $sin \theta/\lambda$ value.

1 Execution

The program lapw3 is executed by invoking the command:

lapw3 lapw3.def or lapw3c lapw3.def or x lapw3 [-c ]

2 Dimensioning parameters

The following parameters are listend in file param.inc_r or param.inc_c :

LMAX2	highest L in in LM expansion of charge and potential
NCOM	number of LM terms in density
NRAD	number of radial mesh points

6 LAPW5 (electron density plots)

This program generates the charge density (or the potential) in a specified plane of the crystal on a two dimensional grid which can be used for plotting with an external contour line program of your choice. Depending on the input files one can generate valence (case.clmval) or difference densities (i.e. crystalline minus superposed atomic densities) using the additional file (case.sigma). In spinpolarized cases one can produce up-, dn- and total densities but also spin densities (difference up-dn). It is also possible to plot total densities (case.clmsum), Coulomb (case.vcoul), exchange-correlation (case.r2v) or total (case.vtotal) potentials, but in those cases the file lapw5.def has to be edited and you must replace case.clmval by the respective filename. The file case.rho contains in the first line

npx, npy, xlength, ylength;

and then the density (potential) written with:

      write(21,11) ((charge(i,j),j=1,npy),i=1,npx)
 11   format(5e16.8)

It is strongly recommended that you use ``Run Programs $\rightarrow$ Tasks $\rightarrow$ Electron density plots'' from w2web, see the TiC example in Fig.3.6 .

1 Execution

The program lapw5 is executed by invoking the command:

lapw5 lapw5.def or lapw5c lapw5.def or x lapw5 [-c -up|dn]

2 Dimensioning parameters

The following parameters are listend in file param.inc:

LMAX2	highest L in in LM expansion of charge and potential
NCOM	number of LM terms in density
NRAD	number of radial mesh points
NPT00	number of radial mesh points beyond RMT
NSYM	order of point group

3 Input

An example is given below. You may want to use XCRYSDEN by T.Kokalj to generate this file (see sect. 9.23.2).

---------------- top of file: case.in5 --------------------
0 0 0 1           # origin of plot: x,y,z,denominator
1 1 0 1           # x-end of plot
0 0 1 2           # y-end of plot
3 3 3             # x,y,z nshells (of unit cells)
100 100           # nx,ny
RHO               # RHO/DIFF/OVER; ADD/SUB or blank
ANG VAL NODEBUG   # ANG/ATU, VAL/TOT, DEBUG/NODEBUG
NONORTHO          # optional line: ORTHO|NONORTHO
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

free format

ix,iy,iz,idv

The plane and section of the plot is specified by three points in the unit cell, an origin of the plot, an x-end and an y-end. The first line specifies the coordinates of the origin, where x=ix/idv, ...in units of the lattice vectors (except fc, bc and c lattices, where the lattice vectors of the conventional cell are used)

line 2:

free format

ix,iy,iz,idv

coordinates of x-end

line 3:

free format

ix,iy,iz,idv

coordinates of y-end (The two directions x and y must be orthogonal to each other unless NONORTHO is selected). Since it is quite difficult to specify those 3 points for a rhombohedral lattice, an auxiliary program rhomb_in5 is provided, which creates those points when you specify 3 atomic positions which will define your plane. You can find this program using ``Run Programs $\rightarrow$ Other Goodies'' from w2web.

line 4:

free format

nxsh, nysh, nzsh

specifies the number of nearest neighbor cells (in x,y,z direction) where atomic positions are generated (needs to be increased for very large plot sections, otherwise some ``atoms'' are not found in the plot)

line 5:

free format

npx, npy

specifies number of grid points in plot. npy=1 produces a file case.rho_onedim containing the distance r (from the origin) and the respective density, which can be used in a standard x-y plotting program.

line 6:

format (2a4)
switch, addsub

switch	RHO	charge (or potential) plots, no atomic density is used (regular case)
	DIFF	difference density plot (crystalline - superposed atomic densities), needs file case.sigma (which is generated with LSTART, see section 6.4)
	OVER	superposition of atomic densities, needs file case.sigma
addsub	NO	(or blank field): use only the file from unit 9
	ADD	adds densities from units 9 and 11 (if present), e.g. to add spin-up and down densities.
	SUB	subtracts density of unit 11 (if present) from that of unit 9 (e.g. for the spin-density, which is the difference between spin-up and down densities).

line 7:

format (3a4)
iunits, cnorm, debug

iunits	ATU	density (potential) in atomic units e/$\mbox{a.u.}^3$ (or Ry)
	ANG	density in e/ ${\mbox{\AA}}^3$ (do not use this option for potentials)
cnorm		determines normalization factor
	VAL	used for files case.clmval, r2v, vcoul, vtotal
	TOT	used for files case.clmsum
debug	DEBU	debugging information is printed (large output)

line 8:

free format

noorth1	ORTHO	(default) enforces directions to be orthogonal
	NONORT	directions can be arbitrary; use this option only if your plotting program supports non orthogonal plots (e.g. for XCRYSDENS).

In order to plot total densities or potentials (see cnorm as above) you have to create lapw5.def using x lapw5 -d, then edit lapw5.def and insert proper filenames (case.clmval, case.r2v, case.vcoul, case.vtotal) for units 9 and 11, and finally run lapw5 lapw5.def.

7 AIM (atoms in molecules)

This program was contributed by:

This program analyses the topology of the electron density according to Bader's ``Atoms in molecules'' theory. For more information see Bader 2001 and Sofo and Fuhr 2001.

The original code has been significantly speeded-up by L.Marks (L-marks@northwestern.edu). There are some new optional keywords in the input (usually not needed, more for testing) and also more debugging output. All changes are described in $WIENROOT/SRC_aim/Notes.txt.

Basically it performs two different tasks, namely searching for ``critical points'' (CP) and/or determination of the atomic basins with an integration of the respective charge density.

It is important that you provide a ``good'' charge density, i.e. one which is well converged with respect to LMMAX in the CLM-expansion (you may have to increase the default LM-list to LM=8 or 10) and with as little ``core-leakage'' as possible (see lstart, sect. 6.4), otherwise discontinuities appear at the sphere boundary.

1 Execution

The program aim is executed by invoking the command:

aim aim.def or aimc aim.def or x aim [-c ]

2 Dimensioning parameters

The following parameters are listed in file param.inc:

LMAX2	highest L in in LM expansion of charge and potential
NRAD	number of radial mesh points
NSYM	order of point group

3 Input

The input file contains ``SWITCHES'', followed by the necessary parameters until an END-switch has been reached.

Examples for ``critical-point'' searches and ``charge-integration'' are given below:

---------------- top of file: case.inaim --------------------
CRIT
1                 # index of the atom (counting multiplicity)
ALL               # TWO/THRE/ALL /FOUR	
3 3 3             # x,y,z nshells (of unit cells)
END
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

A4

CRIT		Keyword to calculate critical points

line 2:

free format

iatom

index of the atom from where the search should be started. This count includes the multiplicity, i.e. if the first atom has MULT=2, the ``second atom'' has iatom=3 (Do not use simply the atom-numbers from case.struct)

line 3:

A4

KEY		TWO, THRE, ALL, or FOUR
		defines the starting point for the CP search to be in the middle of 2, 3 or 4 atoms. ALL combines option TWO and THRE together.

line 4:

free format
nxsh, nysh, nzsh

specifies the number of nearest neighbor cells (in x,y,z direction) where atomic positions are generated.

lines 1-4

can be repeated with different atoms or KEYs

line 5:

A4

END

specifies end of job.

In case.outputaim the critical points are marked with a label :PC

:PC a1 a2 a3 l1 l2 l3 c lap rho iat1 dist1 iat2 dist2

where a1,a2,a3 are the coordinates of the CP in lattice vectors; l1 l2 l3 are the eigenvalues of the Hessian at the CP; c is the character of the CP (-3, -1, 1 or 3); lap is the Laplacian of the density at the CP (lap=l1+l2+l3) and rho is the density at the CP (all in atomic units). In case of a bond critical point (c=-1) also the nearest distances (dist1, dist2) to the two nearest atoms (iat1, iat2) are given.

For convenience run extractaim_lapw case.outputaim (see 5.2.10) and get in the file critical_points_ang a comprehensive list of the CP (sorted and unique) with all values given in $\AA, e/\AA^3, ...$ (instead of bohr).

---------------- top of file: case.inaim --------------------
SURF
3                                atom in center of surface (including MULT)
40 0.0 3.1415926536              theta, 40 points, from zero to pi
40 -0.7853981634 2.3561944902    phi
0.07 0.8 2                       step along gradient line, rmin, check
1.65 0.1                         initial R for search, step (a.u)
3 3 3                            nshell
IRHO                             "INTEGRATE" rho
WEIT                             WEIT (surface weights from case.surf), NOWEIT 
30                               30 radial points outside min(RMIN,RMT)¨
END
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows:

line 1:

A4

SURF

Keyword to calculate the Bader surface.

line 2:

free format

iatom

index of the atom from where the search should be started. This count includes the multiplicity, i.e. if the first atom has MULT=2, the ``second atom'' has iatom=3 (Do not use simply the atom-numbers from case.struct)

line 3:

free format
ntheta, thmin, thmax

	ntheta	number of theta directions for the surface determination. This (and nphi) determines the accuracy (and computing time).
	thmin	starting angle for theta
	thmax	ending angle for theta. If you have higher symmetry, you can change the angles thmin=0, thmax=$\pi$ and use only the ``irreducible'' part, i.e. when you have a mirror plane normal to z (see case.outputs), restrict thmax to $\pi/2$.

line 4:

free format
nphi, phimin, phimax

	nphi	number of phi directions for the surface determination
	phimin	starting angle
	phimax	ending angle. (see comments for theta to reduce phi from the full $0-2\pi$ integration).

line 5:

free format
h0, frmin, nstep

	h0	step in real space to follow the gradient ( 0.1)
	frmin	defines the radius, for which the routine assumes that the search path has entered an atom, given as ``rmin = frmin * rmt'' ( 0.8-1.0)
	nstep	number of steps between testing the position being inside or outside of the surface ( 2-8).

line 6:

free format
r0, dr0

	r0	initial radius for the search of the surface radius ( 1.5)
	dr0	step for the search of the surface radius( 0.1)

line 7:

free format
nxsh, nysh, nzsh

specifies the number of nearest neighbor cells (in x,y,z direction) where atomic positions are generated.

line 8:

A4

IRHO

integrate function on ``unit 9'' (usually case.clmsum) inside previously defined surface (stored in case.surf).

line 9:

A4

WEIT

specifies the use of weights in case.surf.

line 9:

free format

npt

specifies number of points for radial integration outside the MT ( 30)

line 8:

A4

END

specifies end of job.

8 LAPW7 (wave functions on grids / plotting)

This program was contributed by:

The program lapw7 generates wave function data on spatial grids for a given set of $k$-points and electronic bands. lapw7 uses the wave function information stored in case.vector (or in reduced (filtered) form in case.vectorf which can be obtained from case.vector by running the program filtvec). Depending on the options set in the input file case.in7(c) one can generate the real or imaginary part of the wave functions, it's modulus (absolute value) or argument, or the complex wave function itself. For scalar-relativistic calculations both the large and the small component of the wave functions can be generated (only one at a time). The wave functions are generated on a grid which is to be specified in the input file(s). The grid can either be any arbitrary list of points (to be specified free-formatted in a separate file case.grid) or any $n$-dimensional grid ($n=0...3$). The operating mode and grid parameters are specified in the input file case.in7(c). As output lapw7 writes the specified wave function data for further processing - e.g. for plotting the wave functions with some graphical tools such as gnuplot - in raw format to case.psink. For quick inspection, a subset of this data is echoed to the standard output file case.outputf (the amount of data can be controlled in the input). In case, lapw7 is called many times for one and the same wave function, program overhead can be reduced, by first storing the atomic augmentation coefficients $A_{lm}$, $B_{lm}$ (and $C_{lm}$) to a binary file case.abc. For the spin-polarized case two different calculations have to be performed using either the spin-up or the spin-down wave function data as input.

It should be easy to run lapw7 in parallel mode, and/or to apply it to wave function data obtained by a spin-orbit interaction calculation. None of these options have been implemented so far. Also, lapw7 has not yet been adapted for w2web.

Please note: lapw7 requires an LAPW basis set and does not work with APW+lo yet.

1 Execution

The program lapw7 is executed by invoking the command:

lapw7 lapw7.def or lapw7c lapw7.def or x lapw7 [-c] [-up|dn] [-sel]

With the -sel option lapw7 expects data from the reduced (filtered) wave function file case.vectorf, otherwise the standard wave function file case.vector is used. The reduced vector file case.vectorf is assumed to resist in the current working directory, while the standard vector file case.vector (which may become quite large) is looked for in the WIEN scratch directory. For details see lapw7.def.

2 Dimensioning parameters

The following parameters are listed in file param.inc_(r/c):

NRAD		number of radial mesh points
NSYM		order of point group
LMAX7		maximum L value used for plane wave augmentation
LOMAX		maximum L value used for local orbitals

The meaning of LMAX7 is the same as that of LMAX2 in lapw2 and that of LMAX-1 in lapw1. Rather than being an upper bound it directly defines the number of augmentation functions to be used. It may be set different to LMAX2 in lapw2 or LMAX-1 in lapw1, but it must not exceed the latter one. Note that, the degree of continuity of the wave functions across the boundary of the muffin tin sphere is quite sensitive to the choice of the parameter LMAX7. A value of 8 for LMAX7 turned out to be a good compromise.

3 Input

A sample input is given below. It shows how to plot a set of wave functions on a 2-dim. grid.

- - - - - - - - - - - - - - - - - top of file - - - - - - - - - - - - - - - - -
2D ORTHO		# mode O(RTHOGONAL)|N(ON-ORTHOGONAL)
0 0 0 2		# x, y, z, divisor of origin
3 3 0 2		# x, y, z, divisor of x-end
0 0 3 2		# x, y, z, divisor of y-end
141 101 35 25		# grid points and echo increments
NO		# DEP(HASING)|NO (POST-PROCESSING)
RE ANG LARGE		# switch ANG|ATU|AU LARGE|SMALL
1 0		# k-point, band index
- - - - - - - - - - - - - - - - - end of file - - - - - - - - - - - - - - - - -

Interpretive comments on this file are as follows.

line 1:	format(A3,A1)
	mode flag
	mode		the type of grid to be used
		ANY	An arbitrary list of grid points is used.
		0D, 1D, 2D, or 3D	An $n$-dim. grid of points is used. $n = 0, 1, 2, \mbox{or } 3$.
	flag		orthogonality checking flag (for $n$-dim. grids only)
		N	The axes of the $n$-dim. grid are allowed to be non-orthogonal.
		O or $\langle\mbox{blank}\rangle$	The axes of the $n$-dim. grid have to be mutual orthogonal.
line 2:	free format -- (for $n$-dim. grids only)
	ix iy iz idiv		Coordinates of origin of the grid, where x=ix/idv etc. in units of the conventional lattice vectors.
line 3:	free format -- (for $n$-dim. grids with $n>0$ only)
	ix iy iz idiv		Coordinates of the end points of each grid axis. This input line has to be repeated $n$-times.
line 4:	free format -- (not for 0-dim. grids)
	np ... npo ...		In case of an $n$-dim. grid, first the number of grid points along each axis, and then the increments for the output echo for each axis. Zero increments means that only the first and last point on each axis are taken. In case of an arbitrary list of grid points, the total number of grid points and the increment for the output echo. Again a zero increments means that only the first and last grid point are taken. Hence, for $n$-dim. grids, altogether, $2*n$ integers must be provided; for arbitrary lists of grid points two intergers are expected.
line 5:	format(A3)
	tool		post-processing of the wave functions
		DEP	Each wave function is multiplied by a complex phase factor to align it (as most as possible) along the real axis (the so-called DEP(hasing) option).
		NO	No post-processing is applied to the wave functions.
line 6:	format(A3,1X,A3,1X,A5)
	switch iunit whpsi
	switch		the type of wave function data to generate
		RE	The real part of the wave functions is evaluated.
		IM	The imaginary part of the wave functions is evaluated.
		ABS	The absolute value of the wave functions is evaluated.
		ARG	The argument the wave functions in the complex plane is evaluated.
		PSI	The complex wave functions are evaluated.
	iunit		the physical units for wave function output
		ANG	Å units are used for the wave functions.
		AU or ATU	Atomic units are used for the wave functions.
	whpsi		the relativistic component to be evaluated
		LARGE	The large relativistic component of wave function is evaluated.
		SMALL	The small relativistic component of wave function is evaluated.

line 7:	free format
	iskpt iseig
	iskpt		The $k$-points for which wave functions are to be evaluated. Even if the wave function information is read from case.vectorf, iskpt refers to the index of the $k$-point in the original case.vector file! If iskpt is set to zero, all $k$-points in case.vector(f) are considered.
	iseig		The band index for which wave functions are to be evaluated. Even if the wave function information is read from case.vectorf, iseig refers to the band index in the original case.vector file! If iseig is set to zero, all bands (for the selected $k$-point(s)) which can found in case.vector(f) are considered.
line 8:	format(A4) -- this line is optional
	handle		augmentation coefficient control flag
		SAVE or STOR(E)	Augmentation coefficients are stored in case.abc). No wave function data is generated in this case. This option is only allowed if a single wave function is selected in the previous input line.
		READ or REPL(OT)	Previously stored augmentation coefficients are read in (from case.abc). This option is only allowed if the same single wave function as the one who's augmentation coefficients are stored in case.abc is selected in the previous input line.
		anything else	Augmentation coefficients are generated from the wave function information in case.vector(f).

9 FILTVEC (wave function filter / reduction of case.vector)

This program was contributed by:

The program filtvec reduces the information stored in case.vector files by filtering out a user-specified selection of wave functions. Either a fixed set of band indices can be selected which is used for all selected $k$-points (global selection mode), or the band indices can be selected individually for each selected $k$-point (individual selection mode). The complete wave function and band structure information for the selected $k$-points and bands is transferred to case.vectorf. The information on all other wave functions in the original file is discarded. The structure of the generated case.vectorf file is identical to that of the original case.vector file. Hence, it should be possible to use case.vectorf as substitutes for case.vector anywhere in the WIEN program package. (This has only been tested for lapw7.and filtvec.) To filter vector files from spin-polarized calculations, filtvec has to be run separately for both the spin-up and the spin-down files.

filtvec has not yet been adapted for w2web.

1 Execution

The program filtvec is executed by invoking the command:

filtvec filtvec.def or filtvecc filtvec.def or x filtvec [-c] [-up|dn]

In accordance with the file handling for lapw1 and lapw7 the input vector file case.vector is assumed to be located in the WIEN scratch directory, while the reduced output vector file case.vectorf is written to the current working directory. See filtvec.def for details.

2 Dimensioning parameters

The following parameters are listed in file param.inc_(r/c):

NKPT		number of $k$-points
LMAX		maximum number of L values used (as in lapw1)
LOMAX		maximum L value used for local orbitals (as in lapw1)

The parameter LMAX and LOMAX must be set precisely as in lapw1; all other parameters must not be chosen smaller than the corresponding parameters in lapw1.

3 Input

Two examples are given below. The first uses global selection mode; the second individual selection mode.

I. Global Selection Mode

- - - - - - - - - - - - - - - - - top of file - - - - - - - - - - - - - - - - -
3	1 17 33		# number of k-points, k-points
2	11 -18		# number of bands, band indices
- - - - - - - - - - - - - - - - - end of file - - - - - - - - - - - - - - - - -

Interpretive comments on this file are as follows.

line 1:	free format
	kmax ik(1) ... ik(kmax)	Number of $k$-point list items, followed by the list items themselves. Positive list items mean selection of the $k$-point with the specified index; negative list items mean selection of a range of $k$-points with indices running from the previous list item to the absolute value of the current one. E.g. the sequence 2 -5 stands for 2, 3, 4, and 5.
line 2:	free format
	nmax $$ ie(1) ... ie(nmax)	Number of band index items, followed by the list items themselves. Again, positive list items mean selection of a single band index; negative list items mean selection of a range of band indices.

II. Individual Selection Mode

- - - - - - - - - - - - - - - - - top of file - - - - - - - - - - - - - - - - -
2 :				# number of k-points
17	4	11 13 15 17		# k-point, number of bands, band indices
33	3	11 -14 18		# k-point, number of bands, band indices
- - - - - - - - - - - - - - - - - end of file - - - - - - - - - - - - - - - - -

Interpretive comments on this file are as follows.

line 1:	free format
	kmax	the number of individual $k$-points to be selected. This number must be followed by any text, e.g. 'SELECTIONS' or simply ':', to indicate individual selection mode.
line 2:	free format
	ik nmax ie(1) ... ie(nmax)	First the index of the selected $k$-point, then the number of band index items, followed by the list items for the current $k$-point themselves. Positive list items mean selection of the band with the specified index; negative list items mean selection of a range of band indices running from the previous list item to the absolute value of the current one. E.g. the sequence 3 -7 stands for 3, 4, 5, and 7. This input line has to be repeated kmax-times.

10 XSPEC (calculation of X-ray Spectra)

This program calculates near edge structure of x-ray absorption or emission spectra according to the formalism described by Neckel et al.75, Schwarz et al. 79 and 80. For a brief introduction see below. It uses the partial charges in case.qtl. This file must be generated separately using lapw2. Partial densities of states in case.dos1ev are generated using the tetra program. Spectra are calculated for the dipole allowed transitions, generating matrix elements, which are multiplied with a radial transition probability and the partial densities of states. Unbroadened spectra are found in the file case.txspec, broadened spectra in the file case.xspec. Other generated files are: case.m1 (matrix element for the selection rule L+1) and case.m2 (matrix element for the selection rule L-1) and case.corewfx (radial function of the core state). The calculation is done with several individual programs (initxspec, tetra, txspec, and lorentz). which are linked together with the c-shell script xspec.

It is strongly recommended that you use ``Run Programs $\rightarrow$ Tasks $\rightarrow$ X-ray spectra'' from w2web.

1 Execution

1 Execution of the shell script xspec

The program xspec is executed by invoking the command:

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

2 Sequential execution of the programs

Besides calculating the X-ray spectra in one run using the xspec script, calculations can be done ``by hand``, i.e. step by step, for the sake of flexibility.

initxspec

This program generates the appropriate input file case.int, according to the dipole selection rule, for the subsequent execution of the tetra program.

The program initxspec is executed by invoking the command:

initxspec xspec.def or x initxspec [-up|-dn]

tetra

The appropriate densities of states for (L+1) and (L-1) states respectively are generated by execution of the tetra program.

The program tetra is executed by invoking the command:

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

txspec

This program calculates energy dependent dipole matrix elements. Theoretical X-ray spectra are generated using the partial densities of states (in the case.dos1ev file) and multiplying them with the corresponding dipole matrix elements.

The program txspec is executed by invoking the command:

txspec xspec.def or x txspec [-up|-dn]

lorentz

The calculated spectra must be convoluted to account for lifetime broadening and for a finite resolution of the spectrometer before they can be compared with experimental spectra. In the lorentz program a Lorentzian is used to achieve this broadening.

The program lorentz is executed by invoking the command:

lorentz xspec.def or x lorentz [-up|-dn]

2 Dimensioning parameters

The following dimensioning parameters are collected in the files param.inc of SRC_txspec and SRC_lorentz:

IEMAX0	maximum number of energy steps in the spectrum (SRC_lorentz)
NRAD	number of radial mesh points
LMAX	highest l+1 in basis function inside sphere (consistent with input in case.in1)

3 Input

Two examples are given below; one for emission spectra and one for absorption spectra:

1 Input for Emission Spectra:

---------------- top of file: case.inxs --------------------
NbC: C K                (Title)
2               (number of inequivalent atom)
1               (n core)
0               (l core)
0,0.5,0.5       (split, int1, int2)
-20,0.1,3       (EMIN,DE,EMAX   in eV)
EMIS            (type of spectrum, EMIS or ABS)
0.35            (S)
0.25            (gamma0)
0.3             (W)
AUTO            (generate band ranges AUTOmatically or MANually
-7.21           (E0 in eV)
-10.04          (E1 in eV)
-13.37          (E2 in eV)
------------------- bottom of file ------------------------

2 Input for Absorption Spectra:

---------------- top of file: case.inxs --------------------
NbC: C K        (Title)
2               (number of inequivalent atom)
1               (n core)
0               (l core)
0,0.5,0.5       (split, int1, int2)
-2,0.1,30       (EMIN,DE,EMAX   in eV)
ABS             (type of spectrum)
0.5             (S)
0.25            (gamma0)
------------------- bottom of file ------------------------

Interpretive comments on these files are as follows.

line 1:

free format

TITLE

Title

line 2:

free format

NATO

Number of the selected atom (in case.struct file)

line 3:

free format

NC

principle quantum number of the core state

line 4:

free format

LC

azimuthal quantum number of the core state

The table below lists the most commonly used spectra:

Table 8.50: Quantum numbers of the core state involved in the x-ray spectra

Spectrum	n	l
$K$	1	0
$L_{II, III}$	2	1
$M_V$	3	2

line 5

free format

SPLIT, INT1, INT2

split in eV between e.g. $L_{II}$ and $L_{III}$ spectrum (compare with the respective core eigenvalues), INT1 and INT2 specifies the relative intensity between these spectra. Values of 0, 0.5, 0.5 give unshifted spectra.

line 6:

free format

EMIN, DE, EMAX		minimum energy, energy increment for spectrum, maximum energy; all energies are in eV and with respect to the Fermi level
		EMIN and EMAX are only used as limits if the energy range created by the lapw2 calculation (using the QTL switch) is greater than the selected range.

line 7:

Format A4

TYPE	EMIS	X-ray emission spectrum
	ABS	X-ray absorption spectrum (default)

line 8:

free format

S

broadening parameter for the spectrometer broadening. For absorption spectra S includes both experimental and core broadening. Set S to zero for no broadening.

line 9:

free format

GAMMA0

broadening parameter for the life-time broadening of the core states. Set GAMMA0 to zero to avoid lifetime broadening of the core states.

line 10:

free format

W

broadening parameter for the life-time broadening of valence states. Set W to zero to avoid lifetime broadening of the valence states.

line 11:

format A4

BANDRA
	AUTO	band ranges are determined AUTOmatically (default)
	MAN	band ranges have to be entered MANually

line 12:

free format

E0		Emission spectra: onset energy for broadening, E0, generated automatically if AUTO was set in line 10
		Absorption spectra: not used

line 13:

free format

E1		Emission spectra: onset energy for broadening, E1, generated automatically if AUTO was set in line 10
		Absorption spectra: not used

line 14:

free format

E2		Emission spectra: onset energy for broadening, E2, generated automatically if AUTO was set in line 10
		Absorption spectra: not used

11 TELNES3 (calculation of energy loss near edge structure)

This program was contributed by:

The TELNES3 program calculates the double differential scattering cross section (DDSCS) on a grid of energy loss values and impulse transfer vectors. This double differential cross section is integrated to yield a differential cross section, which is written to file. The differential cross section is either a function of energy (ELNES integrated over impulse transfer q); or a function of impulse transfer (ELNES integrated over energy loss E), which shows the angular behavior of scattering.

The DDSCS is calculated as described in a forthcoming publication by K. Jorissen, C. Hebert, and J. Luitz. (The Ph.D. thesis of K. Jorissen (http://www.wien2k.at/reg_user/textbooks/) also describes the formalism onto which TELNES3 is built in great detail.) This formalism allows calculation of relativistic EELS including transitions of ‘arbitrary’ order (i.e., non-dipole transitions). It takes into account the relative orientation between sample and beam. If this is not necessary (because the crystal is isotropic, or the sample is polycrystalline), the formula may be integrated over 4$\pi$, simplifying the calculation. Both scenarios are implemented in TELNES3.

A note to our faithful fans from the early days: it used to be necessary to play such tricks as recompiling lapw2 with lxdos=3 ; to create k-meshes without symmetry ; and to edit case.struct and set ‘ISPLIT’ to 99. This is no longer necessary. Just sit back, relax, and press the buttons in w2web. The integration with the package qtl will do the job.

1 Execution

1 Execution

The program telnes3 is executed by invoking the command:

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

2 Input

TELNES3 requires one input file - case.innes. We recommend using InnesGen $^{\mbox{\textsc{TM}}}$ of w2web to create this input file in a clear and intuitive way. If you wish to manually edit the file, please refer to the following description. Please note that input files created for TELNES2 may or may not work with TELNES3, depending on which optional keywords were used. There isn't a shred of compatibility with the old TELNES program.

The file case.innes consists of two parts: a first block with required input, and a second block with optional input. In fact, the second part may be omitted altogether. The simplest input file looks like this:

Graphite C K edge of first atom.
1            (atom)
1, 0         (n, l core)
285          (E-Loss of 1st edge in eV)
300          (energy of the incident electrons in keV)
0.0 20.0 0.1 (the energy mesh)
5.0  1.87    (collection semiangle, convergence semiangle, both in mrad)
10 1         (NR, NT, defining the integration mesh in the detector plane)
0.8          (spectrometer broadening in eV)
END

This first part of the file is not formatted and contains the following information:

line	value	explanation
1	`Graphite ...'	Title (of no consequence for the calculation)
2	1	Atom number as given in case.struct (the index which numbers inequivalent atoms)
3	1 0	main and orbital quantum number n and l of the core state; eg. 1 0 stands for 1 s
4	285	energy of the edge onset in eV (here for the C K edge)
5	300	beam energy in keV
6	0.0 20.0 0.1	energy mesh given as ; all values in eV. 0.0 is the edge threshold.
7	5.0 1.87	detector collection semiangle and microscope convergence semiangle in mrad
8	10 1	parameters NR and NT which determine the mesh used for sampling the distribution of Q-vectors allowed by collection and convergence angles
9	0.8	spectrometer broadening FWHM in eV
10	END	keyword telling the program that there is no more input to read. Optional keywords and values must be inserted before this line!

There are many other parameters that control the calculation, most of which are set to reasonable default values. To use these advanced parameters, add corresponding keywords before the END keyword. We recommend using InnesGen $^{\mbox{\textsc{TM}}}$ of w2web to create this input file.

Currently, the keywords listed below may be used. Although only the first four characters of each keyword are read, we recommend using the full keyword for clarity.

VERBOSITY
n               eg. : 1

Specifies how much output you'll get. n must be 0 (only basic output; default), 1 (medium output) or 2 (full output, including more technical information).

ATOMS
n1 n2            eg. : 1 3    (default : 1 0 == 1 mult(natom) )

The atom number on line 2 (see above) corresponds to a class of equivalent atoms in case.struct. Equivalent positions n1 to n2 will contribute to the spectrum (default : sum over all atoms in the equivalency class). Since all equivalent atoms have identical electronic structure up to a symmetry operation, this will simply yield a prefactor (n2-n1+1) for the orientation averaged spectrum, but as each equivalent atom has a different orientation with respect to the beam, this setting will influence the shape of an orientation sensitive spectrum.

DETECTOR POSITION
theta_x theta_y            eg. : 0.5 0.5   (default : 0 0)

By default, the detector is aligned with the incoming beam - i.e., source, sample, and detector are connected by a straight line. This card shifts the detector in a plane perpendicular to the incoming beam. The shift is expressed as an angle in mrad. If one draws a line between source and sample, and another line from the sample to the center of the detector aperture, these 2 lines will form an angle of mrad.

MODUS
m                eg. : angles  (default is energy)

The output is a spectrum as a function of energy if m=energy. The output is a spectrum as a function of impulse transfer/scattering angle if m=angle.

SPLIT
splitting energy eg. : 2.7

If the initial state has an orbital quantum number larger than 0, it will generate two superposed edges: one corresponding to $j=l-1/2$, and one corresponding to $j=l+1/2$ (eg., for the 2p initial state we have a L3 and a L2 edge). The splitting energy sets the energy separation of the two edges and should be given in eV (here, L3 is at the energy specified in the beginning of case.innes, and L2 is 2.7 eV higher). By default (keyword omitted), the splitting energy is calculated by the program. It is generally quite accurate.

BRANCHING RATIO
branching ratio  eg. : 1.4

The branching ratio is a scaling factor (eg., here the ratio of intensities L3/L2 would be set to 1.4). By default (keyword omitted), the branching ratio is set to its statistical value of $(2l+2)/2l$.

NONRELATIVISTIC

This key tells the program not to use the relativistic corrections to the scattering cross section. This option generates spectra identical to output of the old TELNES program. This produces incorrect results in many cases. By default, relativistic calculations are done.

INITIALIZATION
make_dos      write_dos      eg. N N (default : Y Y)
make_rot.mat. write_rot.mat  eg. Y N (default : Y Y)

TELNES3 needs many ingredients for its calculations, and this key defines how it gets two of them: the density of states, and the rotation matrices (used for transforming q-vectors from one atom to an equivalent atom). The first entry says whether or not the ingredient has to be calculated (Y : calculate; N : read from file), and the second entry says whether or not the ingredient has to be written to file (Y : write; N : don't write). If make_dos=Y, a file case.qtl must be present from which the dos will be calculated. If make_dos=N, then either a file case.dos or a file case.xdos containing the (x)dos must exist. If make_rot.mat=N, a file case.rotij containg the rotation matrices must exist. If write_rot.mat=Y, a file case.rotij is written. If write_dos=Y, a file case.dos or case.xdos is written. The calculation of the rotation matrices is computationally negligible, but it is recommended to write the xdos to file and not calculate it over and over again.

QGRID
qmodus           eg. L    (U by default)
theta_0               eg. 0.05 (no default value)  )

A collection angle $\alpha$ and convergence angle $\beta$ allow scattering angles up to $\alpha+\beta$ and a corresponding set of Q-vectors. This set (a disk of radius $\alpha+\beta$) is sampled using a discrete mesh. Three types of meshes are implemented :

U

a uniform grid, where each Q-vector samples an equally large part of the disk. Sampling is set up by drawing NR equidistant circles inside the big circle, and choosing $(2i-1)NT$ points on circle $i$, giving $NR^2*NT$ points in total.

L

a logarithmic grid with $NR$ circles. The distance between circles increases exponentially. There are $(2i-1)NT$ points on circle $i$, and $NR^2NT$ points in total. Circle $i$ is at radius theta_0 , where $dx$ depends on $NR$, $\alpha$ and $\beta$.

1

a one dimensional logarithmic mesh; there are $NR$ circles at exponential positions, and only one point on each circle (so $NR$ points in total). This means we sample a line in the detector*beam plane. An economic way of getting spectra as a function of scattering angle in cases with symmetric scattering.

The line specifying theta_0 is to be omitted for the U grid.

ORIENTATION SENSITIVE
g1 g2 g3        (eg. 0.0  40.0  0.0)  (no default value)

This key tells the program not to average over sample to beam orientations, but to use the particular sample to beam orientation defined by the three Euler angles (to be given in degrees). If the ORIENTATION SENSITIVE key is not set, the program will average over all orientations (default).

SELECTION RULE
type            (eg. : q)  (default : n)

The formula for the DDSCS contains an exponential factor in q, which we expand using the Rayleigh expansion. We identify each term in the expansion by the order lambda of the spherical Bessel function it contains. This key keeps some terms and discards others. This can be useful to eliminate unwanted transitions ; to study a spectrum in greater detail ; or simply to speed up the calculation significantly. Possible settings for `type' are :

m   :  use lambda = 0 only
d   :  use lambda = 1 only 
q   :  use lambda = 2 only
o   :  use lambda = 3 only
n   :  no selection rule, calculate all transitions
0-3 :  all transitions up to lambda (eg., 1 means lambda = 0 and 1)

Be aware that the availability of the DOS limits the possible transitions (WIEN2k gives us the DOS only up to l=3). In the nonrelativistic limit, the SELECTION RULE and LSELECTION RULE coincide – i.e., the terms correspond to dipole transitions etc. This is no longer true in the relativistic case.

LSELECTION RULE  
type            (eg. : q)  (default : d)

Whereas the previous key selects transitions by the order of the interaction potential, this key selects them by the L-character of the final states. Possible settings for `type' are (the orbital momentum of the initial state being denoted with l):

m   :  L=l
d   :  L=l +/- 1
q   :  L=l +/- 2
o   :  L=l +/- 3
n   :  no selection rule, calculate all transitions
0-3 :  |L-l| <= type

Be aware that the availability of the DOS limits the possible transitions (WIEN2k gives us the DOS only up to l=3). In the nonrelativistic limit, the SELECTION RULE and LSELECTION RULE coincide – i.e., the terms correspond to dipole transitions etc. This is no longer true in the relativistic case.

EXTEND POTENTIALS
Rmax  sampling  lmax  refine   (e.g.:  3.0  15  0  1.0)  (no defaults)

Calculate matrix elements beyond the muffin tin radius up to r = rmax (in Bohr units). Refine the radial grid by a factor ‘refine’ (1 means default sampling density). This is done by evaluating the potential as given in case.vtotal, which must be present for this type of calculation, and reexpanding it in spherical harmonics, using an angular grid with step of ‘sampling’ degrees, and expanding up to l=lmax. Currently, users should keep lmax to 0 and almost certainly refine to 1.0 . However, advanced users can play around with the software and tweak it to do interesting things if they wish. TELNES3 only requires the spherical potential l=0.

FERMI ENERGY
Ef     (e.g. 0.75)

Manually set the Fermi energy to Ef (needs to be given in Rydberg units). (The default behavior is to get Ef from the header of case.qtl.)

CORE WAVEFUNCTION
filename     (e.g.  case.cwf)

Read the wave function of the initial state from file. (Default behavior is to calculate it instead.)

FINAL STATE WAVEFUNCTION
filename     (e.g.  case.finalwf)

Read the radial wave functions of the final state from file. (Default behavior is to calculate it instead.)

RELATIVISTIC
Itype   (e.g. 1)

Determines which flavor of relativity to use : 0 means nonrelativistic (as in TELNES), 1 means fully relativistic (default), 2 means using the contracted q-vector (only valid for dipole transitions ; as in TELNES2).

NOHEADERS

Don't put headers in output files. This can be helpful if your plotting program doesn't like the headers. (Gnuplot doesn't mind them.)

DOSONLY

Don't calculate the EELS spectrum – halt the program after the calculation of the density of states is finished.

NBTOT
nb    (e.g.    200 )

Arrays for the DOS are first allocated at some initial size, and then reallocated at larger size if necessary. Unfortunately, these reallocation routines appear unstable in some circumstances. This card allows the user to set an array size manually and avoid the need to reallocate (nb is the number of bands). However, very large systems may lead the system to run out of memory and cause a crash.

The following cards are not yet activated (placeholders): TABULATE, SPIN

The following cards are no longer active and must be removed or renamed: XQTL, WRONG.

3 Practical considerations

A typical ELNES calculation consists of the following steps:

• initialize (init_lapw) and converge a SCF calculation (run_lapw)
• provide a suitable case.innes file
• if more excited states are needed than given by the SCF calculation, raise the upper energy limit in case.in1 and run x lapw1
• create the case.qtl file using x qtl -telnes
• calculate the EELS spectrum using x telnes3. It is generally a smart move to make the program calculate the DOS on the biggest energy grid you wiill ever need, save this to file, and simply read it from file for all future calculations (INITIALIZATION key). The same should be done for calculations using EXTEND POTENTIAL (use CORE WAVEFUNCTION key to save to file). This saves time. (In case of disk space problems, once the case.qtl file has been created, the case.vector files can be deleted. Similar, the case.qtl file can be deleted or compressed once the case.dos file exists.)
• add broadening to the spectrum using x broadening. If you wish, editing the case.inb file allows tweaking of the broadening.
• study the output (case.elnes or case.broadspec are the place to start).
• if you wish to do more calculations, save the current results using save_eels -d calculation1 . Edit case.innes and run x telnes3 again.

This sequence can conveniently be executed using w2web by simply clicking one button after the other.

4 Files

TELNES3 uses a lot of files. Many output files are only written if VERBOSITY is set to a high level. Many input files are required only for certain input settings in case.innes. We list here all files possibly used by TELNES3 (and listed in telnes3.def). Each filename is followed by “I� or “O� (input/output), a short description of the file content, and a comment on when the file is used.

• case.innes (I). Defines the ELNES calculation. Always read.
• case.struct (I). Defines the crystal. Always read.
• case.vsp (I). Spherical component of the crystal potential. Read unless core and final state wavefunctions are read from file.
• case.vtotal (I). Total crystal potential (can be generated by lapw0). Read if EXTEND POTENTIAL is used.
• case.rotij (I). Rotation matrices that transform q-vectors between equivalent atoms. Read if INITIALIZATION tells the program to do so.
• case.dos (IO). l-resolved density of states. Read or written depending on INITIALIZATION settings.
• case.xdos (IO). lm,l'm'-resolved density of states. Read or written depending on INITIALIZATION settings; only if the calculation is orientation resolved.
• case.qtl (I). contains partial ‘charge’ components and Fermi energy. Read if DOS needs to be calculated (INITIALIZATION) or if Fermi energy is not specified using FERMI.
• case.inc (I). Specifies core states. Only read if core states are calculated.
• case.kgen (I). contains k-mesh to sample the Brillouin Zone. Read if DOS needs to be calculated.
• case.outputelnes (O). Main log file. Always written. Content depends on VERBOSITY.
• case.elnes (O). Total spectrum. Always written.
• case.sdlm (O). Partial (l,m) spectra. Written if verbosity $>$ 0.
• case.ctr (O). (l,m,l'm') crossterms. Written if verbosity $>$ 0 and calculation is orientation sensitive.
• case.corewavef (O). Contains core wavefunctions. Written if core wavefunctions were calculated and verbosity $>$ 1.
• case.final (O). Contains APW radial basis functions for final states at selected energies. Written if verbosity $>$ 1.
• case.ortho (O). Contains scalar products of initial and final states. Written if verbosity $>$ 1.
• case.matrix (O). Proportionality between partial DOS and spectrum for each l-value. Written if verbosity $>$ 0 and MODUS is energy.
• case.cdos (O). Selected (l,m,l'm') cross-DOS terms. Written if calculation is orientation sensitive and verbosity $>$ 1 or INITIALIZATION causes DOS to be written to file.
• case.sp2 (O). Integrated cross sections as a function of collection angle for all l-values. Written if calculation is orientation sensitive, MODUS is set to angle and verbosity $>$ 1.
• case.angular (O). Differential cross section as a function of scattering angle for all l-values. Written if calculation is orientation sensitive, MODUS is set to angle and verbosity $>$ 1.
• case.inb (O). Settings for the broadening program. Always written.
• case.eelstable (O). Placeholder. Not currently used.
• telnes3.def (I). List of files used by TELNES3. Always read.
• telnes3.error (O). Error file containing current error message; empty after successful calculation. Always written.

12 BROADENING (apply broadening to calculated spectra)

This program was contributed by:

The broadening program can be used in conjunction with the TELNES3 or the xspec program to broaden theoretical spectra by applying a lorentzian broadening for core and valence life times and a gaussian broadening for spectrometer broadening.

1 Execution

1 Execution

The program broadening is executed by invoking the command:

broadening broadening.def or x broadening

2 Input

broadening needs one input file - case.inb. When running TELNES3 this input file is automatically created from settings given in case.innes.

GaN                                                                            
ELNES 
1   1   0 
0.0   1.0    0.0 
0.116  0.116  
1   2.15000000000000 
0.6   
dummy   
0.0 
0.0 
0.0

line	value	explanation
1	`GaN ...'	Title (of no consequence for the calculation)
2	ELNES | ABS | EMIS	Type of input spectrum
3	NC C1 C2	specification of input file: NC number of columns to read, C1 and C2 column to broaden (only in ``ELNES'' mode)
4	SPLIT XINT1 XINT2	split energy, XINT1|2 relative intensities of spectra in C1 and C2
5	GA GB	core hole lifetime of the two edges
6	W WSHIFT	W: type of valence broadening (1: linear with $E/10$, 2: Muller like $E^2$), edge offset
7	S	Spectrometer broadening FWHM in eV
8	dummy	dummy keyword for compatibility with lorentz
9-11	E0, E1, E2	quadratic energy dependent broadening (only used for type ELNES and EMIS when selecting valence broadening type W=2)

13 OPTIMIZE (Volume, c/a or 2-4 dimensional lattice parameter optimization)

This program generates a series of new struct files corresponding to different volumes, c/a ratios, or otherwise different lattice parameters (depending on your input choice) from an existing struct file (either case_initial.struct or case.struct). (When case_initial.struct is not present, it will be generated from the original case.struct.

Furthermore it produces a shell script optimize.job. You may modify this script and execute it. Further analysis of the results (at present only equilibrium volume or c/a ratio are supported in w2web) allows to find the corresponding equillibrium parameters (see Sec.5.3.1).

1 Execution

The program optimize is executed by invoking the command:

optimize optimize.def or x optimize

2 Input

You have to specify interactively which task should be performed (volume, c/a, b/a optimization, or full optimization for tetragonal, orthorhombic or monoclinic structure), how many cases you want to do and how large the change (+/- xx %) should be for each case.

14 ELAST (Elastic constants for cubic cases)

This program was contributed by:

This package calculates elastic constants for cubic crystals. It is described in detail by the author in Charpin 2001. Please note, at http://www.wien2k.at/reg_user/unsupported you can find a package Hex-elastic by Morteza Jamal, which should calculate elastic constants also for hexagonal symmetry.

1 Execution

The package is driven by three scripts:

• init_elast:
  It prepares the whole calculation and should be run in a directory with a valid case.struct and case.inst file. It creates the necessary subdirectories elast, elast/eos, elast/tetra, elast/rhomb, elast/result, the templates for tetragonal and rhombohedral distortion and initializes the calculations using init_lapw.

• elast_setup:
  It should be run in the elast directory, generates the distorted struct-files and eos.job, rhomb.job and tetra.job. These scripts must be adapted to your needs (spin-polarization, convergence,...) and run. elast_setup can be run several times (for different distortions,...).

• ana_elast:
  Once all calculations are done, change into elastresult and run this script. The final results are stored in elastresultoutputs.

• genetempl, setelast, anaelast:
  These three small programs are called by the above scripts.

15 MINI (Geometry minimization)

This program is usually called from the script min_lapw and performs movements of the atomic positions according to the calculated forces (please read Sec. 5.3.2). It generates a new case.struct file which can be used in the next geometry/time step. Depending on the input options, mini helps to find the equilibrium positions of the atoms or performs a molecular dynamics simulation (which might take very long time).

For finding the equilibrium positions different methods are available. We recommend PORT, a ``reverse-communication trust-region Quasi-Newton method'' from the Port library (http://www.bell-labs.com/project/PORT/doc/port3doc.tar.gz, Gay 1983), which was implemented by L.D.Marks (L-marks@northwestern.edu, http://www.numis.northwestern.edu). It minimizes the total energy and NOT the forces (using the forces as derivative of E vs. atomic positions). In cases when energy and forces are not "compatible", eg. because of numerical noise due to limited scf convergence, small RKmax or crude k-mesh, PORT may fail. An interesting alternative is a sophisticated modified steepest-descent method (NEW1), which minimizes the forces (does not use the total energy). Eventually a damped Newton dynamics is also available.

The forces are read from a file case.finM, while the ``history'' of the geometry optimization or MD is stored in case.tmpM

One can constrain individual positions in case.inM or define linear constrains for several positions using case.constraint (thanks to B.Yanchitsky (Kiev, yan@imag.kiev.ua); for details see comments in the SRC_templates/case.constraint file). In case of calculations with linear constrains one should use NEW1 (in case.inM). When constraining individual positions and using PORT, one should after modifications in case.inM rerun x pairhess -copy (which copies .minpair to .minrestart and .min_hess).

1 Execution

The program mini is executed by invoking the command:

mini mini.def or x mini

2 Dimensioning parameters

The following dimensioning parameters are collected in the file param.inc:

MAXIT	maximum number of geometry steps
NRAD	number of radial mesh points
NCOM	number of LM terms in density
NNN	number of neighboring atoms for nn
NSYM	order of pointgroup

3 Input

Two examples are given below; one for a PORT geometry optimization, and one for molecular dynamics using a NOSE thermostat:

1 Input for geometry optimization:

---------------- top of file: xxx.inM --------------------
PORT 2.0  0.25         (PORT/NEWT   tolf step0 (a4,2f5.2))
1.0 1.0 1.0 3.0        ( 1..3:delta, 4:BO/eta(1=friction zero))
1.0 1.0 1.0 6.0        ( 1..3=0 constraint)
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows.

line 1:

format(a4,2f5.2)

MINMOD		Modus of the calculation
	PORT	Geometry optimization with reverse-communication trust-region Quasi-Newton routine from the Port library. Recommended option.
	NEW1	Performs geometry optimization with "sophisticated" steepest-descent method with automatic adaptation of stepsize (still experimental, but when PORT fails, an interesting alternative)
	NEWT	Performs geometry optimization with damped Newton scheme according to
		$R^{\tau+1}_m = R^\tau_m + \eta_m (R^\tau_m - R^{\tau-1}_m) + \delta_m F^\tau_m$
		where $R^\tau_m$ and $F^\tau_m$ are the coordinate and force at time step $\tau$. When the force has changed its direction from the last to the present timestep (or is within the tolerance TOLF), $\eta_m$ will be set to $1-\eta_m$. Please see also the comments in Sect. 5.3.2
	BFGS	Performs geometry optimization with the variable metric method of BFGS. This option works only when a quadratic approximation is a good approximation to the specific potential surface. Obsolete.
TOLF		Force tolerance, geometry optimization will stop when all forces are below TOLF.
STEP0		Initial "Trust-region radius". Determines size of first geometry step.

line 2:

free format

DELTA(1-3)		For PORT (and BFGS): Precondition parameters: rescales the gradient and thus determines the size of the geometry steps
		For NEWT/NEW1: x,y,z-delta parameters. Determines speed of motion. Good values must be found for each individual system. They depend on the atomic mass, the vibrational frequencies and the starting point (see Sect. 5.3.2).
		DELTA(i) = 0 constrains the corresponding i-th coordinate (for PORT: after setting a DELTA(i)=0, also rerun pairhess to set a proper Hessian). The delta-x,y,z correspond to the global coordinates (the same as the positions in case.struct and the forces :FGL from case.scf).
		Whenever you change these DELTA(i) you must remove file case.tmpM !
ETA		For NEWT: damping (friction) parameter. ETA=1 means no friction, ETA=0 means no speed from previous time steps
		PORT: changes the strength of the bonds when running pairhess and ZWEIGHT is negative (see the pairhess description), otherwise not used
		NEW1: ETA is not used

$»>$ line 2:

must be repeated for every atom

2 Input for Molecular dynamics:

---------------- top of file: nbc.inM --------------------
NOSE                      (NOSE/MOLD (a4))
58.9332 400. 1273. 5.0    (Masse, delta t, T, nose-frequency) 
58.9332 400. 1273. 5.0    
58.9332 400. 1273. 5.0    
58.9332 400. 1273. 5.0
58.9332 400. 1273. 5.0
58.9332 400. 1273. 5.0
------------------- bottom of file ------------------------

Interpretive comments on this file are as follows.

line 1:

format(a4,f5.2)

MINMOD		Modus of the calculation
	MOLD	Performs next molecular dynamics timestep
	NOSE	Performs next molecular dynamics timestep using a NOSE thermostat

line 2:

free format

MASS		Atomic mass of $i^{th}$ atom
TIMESTEP		Time step of MD (in atomic units, depends on highest vibrational frequencies)

TEMP

Simulation Temperature (K)

NOSF

Nose-frequency

$»>$line 2:

must be repeated for every atom

16 OPTIC (calculating optical properties)

This program was contributed by:

The theoretical background is described in detail in Ref. Abt 1994 and Ambrosch-Draxl 06 (Please cite the latter when publishing optics results!). The calculation of optical properties requires a dense mesh of eigenvalues and the corresponding eigenvectors. For that purpose start kgen and generate a fine k-mesh (with many k-points). Run lapw1 and then lapw2 with the option FERMI (Note: You must also put TETRA / with value=101. for metallic systems case.in2) in order to generate the weight-file. After the vector-file has been generated by lapw1 run optic in order to produce the momentum matrix elements. Then the program joint carries out the BZ integration and computes the imaginary part of the complex dielectric tensor. In order to obtain the real part of the dielectric tensor kram may be executed which uses the Kramers-Kronig relations.

The program optic generates the symmetrized squared momentum matrix elements

${\bf M}_i = <n^{\prime}{\vec k} \vert {\vec p}.{\vec e}_i \vert n{\vec k}>^2$

between all band combinations for each k-point given in the vector-file and stores them in case.symmat. For the orthogonal lattices the squared diagonal components can be found in the file case.mat_diag. For non-orthogonal systems all 6 components $(M_j)^* M_k$ can be calculated according to the symmetry of the crystal. In systems without inversion symmetry the complex version opticc must be executed.

The matrix elements (and the imaginary part of the dielectric tensor) are given per spin in case of the spin-polarized calculation and as a sum of both spin directions if the calculation is non-spinpolarized.

Due to spin-orbit coupling imaginary parts of the nondiagonal elements may occur in spinpolarized cases. Thus in general, up to 9 components can be calculated at the same time.

Since version WIEN2k_11.1 an option for the calculation of XMCD (X-ray magnetic circular dichroism) has been added by Lorenzo Pardini (loren.pard@gmail.com). In the case of the XMCD calculation, the momentum matrix elements in the dipole approximation between the selected core state and conduction states are stored in case.symmat1up (higher energy core state, eg. L) and case.symmat2up (lower energy core state, eg. L) for each k-point and every band. For K, L, and M edges, only case.symmat1up is written, since in these cases there is only one edge, whereas both case.symmat1up and case.symmat2up are written for the remaining cases.
XMCD calculation can be only performed for system with spin-polarized AND spin-orbit set up.
In order to calculate XMCD and x-ray absorption spectra, eigenvalues must be evaluated over a mesh in the whole Brillouin zone; for that porpouse, the following procedure should be followed:

• copy case.struct to case.ksym (cp case.struct case.ksym) and remove all the symmetry operations but the identity;
• generate a k-mesh in the whole Brilouin zone (x kgen -so);
• change TOT to FERMI in case.in2c;
• set IPRINT=1 in case.inc to activate core-wavefunction output;
• for metallic systems, put TETRA with value 101;
• execute runsp_lapw -so -s lapw1 -e lcore;
• run optic: x optic -c -so -up;
• run joint: x joint -up.

You must not use p-1/2 ``relativistic'' LOs in LAPWSO, since this basis is not supported on OPTICS yet.

1 Execution

The program optic is executed by invoking the command:

optic(c) optic.def or x optic [-c -up|dn -so -p]

Recommended procedure for spin-orbit coupling:

In order to get the correct matrix elements, the files case.vectorso[up|dn] have to be used. For that purpose the following procedure is recommended:

• run SCF cycle: run[sp]_lapw -so
• generate a fine k-mesh for the optics part: x kgen [-so (if case.ksym has been created by symmetso) ]
• change TOT to FERMI in case.in2c
• execute run[sp]_lapw -so -s lapw1 -e lcore with this fine k-mesh
• run optic: x opticc -so [-up]
• run joint: x joint [-up]
• run kram: x kram [-up]

In cases of non-spinpolarized spin-orbit calculations WITHOUT inversion symmetry one must do some tricks and ``mimick'' a spinpolarized calculation:

• cp case.vsp case.vspup
• cp case.vsp case.vspdn
• cp case.vectorso case.vectorsoup

• x lapw2 -fermi -so -c
• cp case.weight case.weightup
• cp case.weight case.weightdn

• x optic -so -up
• x joint -up

Due to the ``paramagnetic'' weight files (which are normalized to 2 electrons per band instead of one) all your results (joint/sigma...) must be divided by a factor of two.

Note: In spin-polarized cases with spin-orbit only one call to optic, joint and/or kram (either up or down) is necessary, since the spins are not independent any more and both vector-files are used at the same time.

2 Dimensioning parameters

The following dimensioning parameters (listed in param.inc_r and param.inc_c) are used:

LMAX	highest l+1 in basis function inside sphere (reducing LMAX to 4 or 5 may dramatically speed-up optics for large cases, but of course the matrix elements will be truncated and do not have full precision)
LOMAX	highest l for local orbital basis (consistent with input in case.in1)
NRAD	number of radial mesh points
NSYM	order of point group

3 Input

An example is given below:

---------------- top of file: case.inop --------------------
99999 1          : NKMAX, NKFIRST
-5.0 2.0  18     : EMIN, EMAX, NBvalMAX
XMCD  1  L23      : optional line: for  XMCD of 1st atom and L23 spectrum
2                : number of choices (columns in *symmat)
1                : Re xx
3                : Re zz
OFF              : ON/OFF   writes MME to unit 4
------------------- bottom of file -------------------------

Interpretive comments on this file are as follows:

line 1:

free format

nkmax, nkfirst

maximal number of k-points , number of k-point to start calculation

line 2:

free format

emin, emax		absolute energy range (Ry) for which matrix elements should be calculated
nbvalmax		optional input. Setting this to the number of occupied bands (see case.output2) will reduce cpu-time of optics (for large cases, MM only between occupied and empty bands)

line 3:

optional line, must be omitted for ``normal'' optic; free format

XMCD		fixed keyword to indicate XMCD calculation. You should also use NCOL=6
natom		atom number (from case.struct file) for which XMCD should be calculated
edge		specify the edge: must be K, L1, L23, M1, M23 or M45

line 3+:

free format

ncol

number of choices (columns in case.symmat)

line 4+:

free format

icol		column to select. Choices are:
		1 ...Re $<x><x>$
		2 ...Re $<y><y>$
		3 ...Re $<z><z>$
		4 ...Re $<x><y>$
		5 ...Re $<x><z>$
		6 ...Re $<y><z>$
		7 ...Im $<x><y>$
		8 ...Im $<x><z>$
		9 ...Im $<y><z>$
		Options 7-9 apply only in presence of SO, options 4-6 only in non-orthogonal cases.

line 5:

free format
IMME, NATOMS (optional input)

IMME		OFF/ON; optionally prints unsquared momentum matrix elements to unit 4
NATOMS		number of atoms for which the opt. matrix elements should be calculated (The index of the atoms is read in the next line). Please note, that since we need the squared matrix elements, the sum of $\epsilon_2$ using atom ``1'' and atom ``2'' separately is NOT the same as using atom ``1 and 2'' together, since we miss crossterms. Nevertheless this can be a useful option to analyze the origin of certain peaks in $\epsilon_2$. I recommend to repeat this analysis for all possible combinations, and also for a list of ``all'' atoms, since this shows the effect of the interstitial (and crossterms involving the interstitial).

line 6: (optional)

free format

IATOMS

List of NMODUS atoms for which the opt. matrix elements should be calculated (see above).

17 JOINT (Joint Density of States)

This program was contributed by:

This program carries out the BZ integration using the momentum matrix elements case.symmat calculated before by optic. The interband or the intraband contributions to the imaginary part of the dielectric tensor ($\epsilon_2$) can be computed. Alternatively, the DOS or the joint DOS can be derived.

The output in case.joint can be plotted with any xy-plotting package or opticplot_lapw or Curve_lapw.

Warning: Negative values for $\epsilon_2$ may occur due to negative weights in Blöchl's tetrahedron method.

For optional XMCD calculations (see OPTICS) an integration of the Brillouin zone is carried out using the momentum matrix elements from case.symmat1up and case.symmat2up (if both edges are present, otherwise only from case.symmat1up). The broadened and unbroadened spectra are written in files case.xmcd and case.rawxmcd: in these files, the first coloumn is the energy mesh, the second and third coloumns the left and right polarized absorption spectra, the fourth column the XMCD and the last is the XAS. For L, M, and M edges, the broadened and unbroadened spectra for the single edges (useful for the application of Carra's and Thole's sum rules) are stored in case.broad1 and case.broad2 and case.raw1 and case.raw2, respectively, where "1" and "2" are refererred to the higher and lower energy core state.

1 Execution

The program joint is executed by invoking the command:

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

2 Dimensioning parameters

The following parameter is listend in files param.inc:

NSYM	order of point group
MG0	number of columns (usually 9)

3 Input

An example is given below:

---------------- top of file: case.injoint -----------------------
    1  9999  8                : LOWER,UPPER,upper-valence BANDINDEX
  -0.0000    0.00100   2.0000 : EMIN DE EMAX FOR ENERGYGRID IN ryd 
eV                            : output units  eV / ryd  
XMCD                            : omitt these 4 lines for non-XMCD
 -49.88 -50.80                  : core energies in Ry (grep :2P case.scfc)
   1.6    0.6                   : core-hole broadening (eV) for both core states
   0.1                          : spectrometer broadening (eV)
    4                         : SWITCH 
    2                         : NUMBER OF COLUMNS
   0.1  0.1  0.3              : BROADENING (FOR DRUDE MODEL - switch 6,7)
------------------- bottom of file -------------------------------

Interpretive comments on this file are as follows:

line 1:

free format

b1, b2, b3

lower, upper and (optional) upper-valence band-index (Setting b3 may allow for additional analysis (restricting the occupied bands from b1-b3) and in big cases it will reduce memory requirements. Otherwise set b3 equal b2)

line 2:

free format

emin, de, emax

Energy window and increment in Ry (emin must not be negative)

line 3:

free format

units	eV	output in units of eV
	Ry	output in units of Ry

line 4:

optional line for XMCD, must be omitted for ``normal'' optic; free format

XMCD

keyword for XMCD calculation, requires 3 more lines

line 4xmcd:

must be omitted for ``normal'' optic; free format

E_core1, E_core2

lower and higher core energies (in Ry, get them using eg. ``grep :2P case.scf'')

line 4xmcd:

must be omitted for ``normal'' optic; free format

broad_core1, broad_core2

lifetime broadening (eV) of lower and higher core state

line 4xmcd:

must be omitted for ``normal'' optic; free format

broad

spectrometer (Gaussian) broadening (eV)

line 4+:

free format

switch	0	joint DOS for each band combination
	1	joint DOS as sum over all band combinations
	2	DOS for each band
	3	DOS as sum over all bands
	4	imaginary part of the dielectric tensor ($\epsilon_2$)
	5	imaginary part of the dielectric tensor for each band combination
	6	intraband contributions: number of ``free`` electrons per unit cell assuming bare electron mass (calculated around $E_F \pm 10*de$ as defined in input line 4), plasma-frequency
	7	in addition to switch 6 the contributions from different bands to the plasma frequency are analyzed.

line 5:

free format

ncol

number of columns

line 6:

free format
broadening

x,y,z

broadening parameters (in units defined in line 3) for Drude-model

The band analysis for all options (switches 0, 2, 5, and 7) has been improved: For each tensor component additional files are created, where each column contains the contributions from a single band or band combination. The file names are e.g. .Im_eps_xx_1, .Im_eps_xx_2, or .jdos_1 etc. where the number of files depend on the number of bands/band combinations.

Warning: The number of band combinations might be quite large!

18 KRAM (Kramers-Kronig transformation)

This program was contributed by:

The Kramers-Kronig analysis is carried out for the actual number of columns contained in the case.joint[up|dn] file. For each real component its imaginary counterpart is created and vice versa. All dielectric tensor components can be found in file case.epsilon[up|dn]. The real and imaginary parts of the optical conductivity (in $10^{15}/s$) are written to file case.sigmak[up|dn]. In addition, file case.absorp contains the real parts of the optical conductivity (in $1/(\Omega cm)$ and the absorption coefficients. The loss function is also calculated (case.eloss), where for the previously calculated Plasma-frequency the intraband contributions can be added.

Please note, that for spin-polarized calculations, the Kramers-Kronig analysis is NOT really additive, i.e. most quantities (like ) cannot be obtained by simply adding the spin-up and dn results to get the total contribution (see equations in Ambrosch 06). Thus, one should add up both spin contributions of $\epsilon_2$ (in case.jointup and case.jointdn) using addjoint-updn_lapw (this will produce case.joint) before calling (non-spinpolarized) x kram.

The 3 sumrules are also checked and written to case.sumrules.

The output in case.epsilon[up|dn] and case.sigmak[up|dn] can be plotted with any xy-plotting package, opticplot_lapw or the "OPTIC"-task in w2web.

1 Execution

The program kram is executed by invoking the command:

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

2 Dimensioning parameters

The following parameters are listed in files param.inc:

MAXDE	maximum number of points in energy mesh
MPOL	fixed at 6

3 Input

An example is given below:

---------------- top of file: case.inkram -----------------------
 0.0    gamma for Lorentz broadening (in units selected in joint)
 0.0    energy shift (scissors operator) (in units selected in joint)
  1      add intraband contributions? yes/no: 1/0
 12.60   plasma frequencies (for each ``column'' in case.injoint) 
  0.20   Gammas for Drude terms (for each ``column'' in case.injoint) 
------------------- bottom of file -------------------------------

Interpretive comments on this file are as follows:

line 1:

free format

EGAMM

Lorentz broadening (in energy units selected in joint)

line 2:

free format

ESHIFT

Energy shift (scissors operator) (in energy units selected in joint)

line 3:

free format

INTRA	0	Intraband contributions are not added
	1	Intraband contributions are added (requires plasma-frequencies calculated by joint using switch ``6'') )

line 4:

free format

EPL

Plasma-frequencies (calculated by joint using SWITCH=6 for all columns)

line 5:

free format

EDRU

Broadening for Drude terms (for all columns)

19 DIPAN (Dipolar anisotropies)

This program was contributed by:

This program calculates the magnetic dipolar hyperfine field and the dipolar magnetocrystalline anisotropy by a direct lattice summation over the magnetic moments of all sites.

According to Wikipedia

(13)

where .
is direction of magnetization.
is permeability of free space; H/m.
is the dipolar field in T.
is magnetic dipolar moment in Am = J/T, assumed to be parallel to .
$r$ is in m.
We want to express $\mu$ in Bohr magnetons =9.274078.10 J/T and
$r$ in atomic units for length (Bohr radius) =5.2917706.10 m.
Inserting in (1) gives

(14)

Total dipolar field acting on atom $i$ is given by the lattice sum

(15)

Dipolar anisotropy energy is given by the sum

(16)

when the sum is over atoms in the unit cell, $V$ is the unit cell volume, Factor 1/2 appears because of the double summation.

Expressing in T, in and $V$ in (a.u.) gives

(17)

1 Execution

The program dipan is executed by invoking the command:

dipan dipan.def or x dipan

2 Dimensioning parameters

The following parameters are listed in files dipan.f:

NATO	number of inequivalent atoms in unit cell
NDIF	total number of atoms in unit cell

3 Input

An example is given below:

---------------- top of file: case.indipan -----------------------
 160.   0                      Rmax (a.u.),  ipr (printing option)
 -0.26                         Magnetic moment of 1s atom (Y) in mu_B
  1.525                        Magnetic moment of 2nd atom  (Co(2c))
  1.529                        Magnetic moments of 3rd atom (Co(3g)) in mu_B
 1381.                         Volume in a.u.**(-3)
 2                             ndir: numder of magnetization directions    
 0. 0. 1.                      first direction for the magnetization
 1. 1. 0.                      second direction
------------------- bottom of file -------------------------------

Interpretive comments on this file are as follows:

line 1:

free format
Rmax, IPR

Rmax		max distance (bohr) for lattice summation. Vary it for convergence check.
IPR		Print switch. IPR=2 produces very large files case.outputdipan and case.nn_dipan

line 2:

free format

mm

Magnetic moment () of first atom

line 2

must be repeated for every non-equivalent atom in the unit cell

line 3:

free format

VOLUME

Unit cell volume in bohr**3 (grep :VOL case.scf)

line 4:

free format

NDIR

number of magnetization directions for which the dipolar contributions will be calculated. For the differences are also calculated.

line 5:

free format

h,k,l

direction of magnetization

line 5

must be repeated NDIR times

20 FSGEN (Fermi-surface generation)

Unfortunately there is no really versatile tool for Fermi surface generation or analyzing FS properties. We have collected here a series of small programs together with some description on how to proceed to generate 2D-Fermisurfaces within WIEN.

• As usually, you have to run an scf cycle and determine a good Fermi-energy. "Good" means here a Fermi-energy coming from a calculation with a dense k-mesh.

• You should than create a mesh within a plane of the BZ, where you want to plot the FS. Some utility programs like sc_fs_mesh, (fcc, bcc, cxz_mon and hex are also available) may help you here, but only some planes of the BZ have been implemented so far. Please check these simple programs and modify them according to your needs. Copy the generated k-mesh fort.2 to case.klist.

• Run lapw1 with this k-mesh.

• Run spaghetti with input-options such that it prints the bands which intersect EF to
  case.spaghetti_ene (line 10, see sec. 8.3)

• Edit case.spaghetti_ene and insert a line at the top:
  NX, NY, x-len, y-len, NXinter, NYinter, Invers, Flip
  where
  NX, NY are the number of points in the two directions
  x-len, y-len are the length of the two directions of the plane (in bohr$^{-1}$, you can find this in case.spaghetti_ene)
  NXinter, NYinter: interpolated mesh, e.g. 2*NX-1
  Invers: 0/1: mirrors x,y
  FLIP: 0/1: flips x,y to y,x

• Run spagh2rho < case.spaghetti_ene to convert from this format into a format which is compatible with the case.rho file used for charge density plotting. It generates files fort.11, fort.12, ... (for each band separately) and you should use your favorite plotting program to generate a contourplot of the FS (by using a contourlevel = 0). Alternatively you can use for plotting:

• Run fsgen_lapw 11 xx save_filename, which is a small shell script that can plot all fermi surfaces using the data-files fort.11, fort.12, ... fort.xx generated in the previous steps. It requires the public domain package pgplot and the contour-plot program plotgenc. (The latter can be obtained from http://www.wien2k.at/reg_user/unsupported/, but you must have installed the pgplot library before.)