WIEN2k User Guide for Crystal Properties
WIEN2k User Guide for Crystal Properties
Peter Blaha
Karlheinz Schwarz
Georg Madsen
Dieter Kvasnicka
Joachim Luitz
ISBN 3-9501031-1-2
Contents
1 Introduction 1
2 Basic concepts 7
2.1 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 The APW Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 The LAPW Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 The APW+lo Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.3 General considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3 Quick Start 13
3.1 Naming conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Starting the server . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Connecting to the w2web server . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Creating a new session . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5 Creating a new case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.6 Creating the struct file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.7 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.8 The SCF calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.9 The [Link] file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.10 Saving a calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.11 Calculating properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.11.1 Electron density plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.11.2 Density of States (DOS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.11.3 X-ray spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.11.4 Bandstructure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.11.5 Bandstructure with band character plotting / full lines . . . . . . . . . . . . . 27
3.11.6 Volume Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.12 Setting up a new case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.12.1 Manual setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.12.2 Setting up a new case using w2web . . . . . . . . . . . . . . . . . . . . . . . . 29
3
II Detailed description of the files and programs of the WIEN2k package 31
5 Shell scripts 51
5.1 Job control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.1.1 Main execution script (x lapw) . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.1.2 Job control for initialization (init lapw) . . . . . . . . . . . . . . . . . . . . . . 52
5.1.3 Job control for iteration (run lapw or runsp lapw) . . . . . . . . . . . . . . . . 53
5.2 Utility scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2.1 Save a calculation (save lapw) . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2.2 Restoring a calculation (restore lapw) . . . . . . . . . . . . . . . . . . . . . . . 56
5.2.3 Remove unnecessary files (clean lapw) . . . . . . . . . . . . . . . . . . . . . . 56
5.2.4 Migrate a case to/from a remote computer (migrate lapw) . . . . . . . . . . . 56
5.2.5 Generate [Link] (instgen lapw) . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.6 Set R-MT values in your [Link] file (setrmt lapw) . . . . . . . . . . . . . . 57
5.2.7 Create [Link] file (for DOS) (configure int lapw) . . . . . . . . . . . . . . . . 57
5.2.8 Check for running WIEN jobs (check lapw) . . . . . . . . . . . . . . . . . . . . 58
5.2.9 Cancel (kill) running WIEN jobs (cancel lapw) . . . . . . . . . . . . . . . . . . 58
5.2.10 Extract critical points from a Bader analysis (extractaim lapw) . . . . . . . . . 58
5.2.11 scfmonitor lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.2.12 analyse lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.2.13 Check parallel execution (testpara lapw) . . . . . . . . . . . . . . . . . . . . . 59
5.2.14 Check parallel execution of lapw1 (testpara1 lapw) . . . . . . . . . . . . . . . 59
5.2.15 Check parallel execution of lapw2 (testpara2 lapw) . . . . . . . . . . . . . . . 59
5.2.16 grepline lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2.17 initso lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2.18 vec2old lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2.19 clmextrapol lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3 Structure optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3.1 Lattice parameters (Volume, c/a, lattice parameters) . . . . . . . . . . . . . . . 61
5.3.2 Minimization of internal parameters (min lapw) . . . . . . . . . . . . . . . . . 63
5.4 Phonon calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.4.1 init phonon lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.4.2 analyse phonon lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.5 Parallel Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.5.1 k-Point Parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.5.2 MPI parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.5.3 How to use WIEN2k as a parallel program . . . . . . . . . . . . . . . . . . . . 68
5.5.4 The .machines file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.5.5 How the list of k-points is split . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.5.6 Flow chart of the parallel scripts . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.5.7 On the fine grained parallelization . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.6 Getting on-line help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.7 Interface scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.7.1 eplot lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.7.2 parabolfit lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.7.3 dosplot lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.7.4 dosplot2 lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.7.5 Curve lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.7.6 specplot lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.7.7 rhoplot lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.7.8 opticplot lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.7.9 addjoint-updn lapw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6 Initialization 75
6.1 NN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.1.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.2 SGROUP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.2.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3 SYMMETRY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4 LSTART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.5 KGEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.5.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.5.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.6 DSTART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.6.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.6.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7 SCF cycle 83
7.1 LAPW0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.1.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.1.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.1.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.2 ORB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
7.2.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.2.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.2.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.3 LAPW1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
7.3.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.3.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.3.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.4 LAPWSO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.4.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.4.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.4.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
7.5 LAPW2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.5.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.5.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.5.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.6 SUMPARA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.6.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.6.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.7 LAPWDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.7.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.7.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.7.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.8 LCORE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.8.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.8.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7.8.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7.9 MIXER . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
7.9.1 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
7.9.2 Dimensioning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
7.9.3 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
13 References 195
IV Appendix 199
11
List of Figures
2.1 Partitioning of the unit cell into atomic spheres (I) and an interstitial region (II) . . . 8
3.1 TiC in the sodium chloride structure. This plot was generated using BALSAC (see
9.23.1). Interface programs between WIEN2k and BALSAC are available. . . . . . . 14
3.2 Startup screen of w2web . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Main window of w2web . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4 StructGen of w2web . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5 List of input files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.6 Task “Electron Density Plots” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.7 Electron density of TiC in (100) plane using Xcrysden . . . . . . . . . . . . . . . . . . 23
3.8 Electron density of TiC in (100) plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.9 Density of states of TiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.10 Density of states of TiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.11 Ti LIII spectrum of TiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.12 Bandstructure of TiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.13 Bandstructure of TiC, showing t2g-character bands of Ti in character plotting mode . 28
3.14 Energy vs. volume curve for TiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.1 Data flow during a SCF cycle ([Link], [Link], [Link], [Link] and
optional files are omitted) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.2 Program flow in WIEN2k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
13
Licence conditions of WIEN2k
P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka and J. Luitz
DEFINITIONS:
In the following, the term “the authors”, refers to P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvas-
nicka and J. Luitz at the above address. “Program” shall mean that copyrighted APW+LO code
(in source and object form) comprising the computer programs known as WIEN2k or the graphical
user interface w2web.
1. All title, ownership and rights to the program or to copies of it remain with the authors,
irrespective of the ownership of the media on which the program resides.
2. I will not supply a copy of the code to anyone for any reason whatsoever. This in no way
limits my making copies of the code for backup purposes, or for running on more than one
computer system at my institution (it is a site license for the registered group). I will refer
any request for copies of the program to the authors.
3. I will not incorporate any part of WIEN2k or w2web into any other program system, without
prior written permission of the authors.
4. I will keep intact all copyright notices.
5. I understand that the authors supply WIEN2k and w2web and its documentation on an “as
is” basis without any warranty, and thus with no additional responsibility or liability. I agree
to report any difficulties encountered in the use of WIEN2k or w2web to the authors.
6. In any publication in the scientific literature I will reference the program as follows:
P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka and J. Luitz, WIEN2k, An Aug-
mented Plane Wave + Local Orbitals Program for Calculating Crystal Properties
(Karlheinz Schwarz, Techn. Universität Wien, Austria), 2001. ISBN 3-9501031-1-2
Please enter your publications with WIEN2k on our web-page for “papers”, so that we can
easily include them in the list of WIEN-publications. In addition we like to receive a copy
(ps-, pdf-file or reprint), especially for less common journals. Please send it to the second
author, K. Schwarz.
7. It is understood that modifications of the WIEN2k or the w2web code can lead to problems
where the authors may not be able to help. Please report useful modifications or major ex-
tensions to the authors.
8. I understand that support for running the program can not be provided in general, except on
the basis of a joint project between the authors and the research partner.
ii
1 Introduction
The Linearized Augmented Plane Wave (LAPW) method has proven to be one of the most accurate
methods for the computation of the electronic structure of solids within density functional theory.
A full-potential LAPW-code for crystalline solids has been developed over a period of more than
twenty years. A first copyrighted version was called WIEN and it was published by
In the following years significantly improved and updated UNIX versions of the original WIEN-
code were developed, which were called WIEN93, WIEN95 and WIEN97. Now a new version,
WIEN2k, is available, which is based on an alternative basis set. This allows a significant improve-
ment, especially in terms of speed, universality, user-friendliness and new features.
WIEN2k is written in FORTRAN 90 and requires a UNIX operating system since the programs are
linked together via C-shell scripts. It has been implemented successfully on the following computer
systems: Pentium systems running under Linux, IBM RS6000, HP , SGI , Compac DEC Alpha, and
SUN. It is expected to run on any modern UNIX (LINUX) system.
Hardware requirements will change from case to case (small cases with 10 atoms per unit cell can
be run on any Pentium PC with 128 Mb under Linux), but generally we recommend a powerful PC
or workstation with at least 256 Mb (better 512 Mb or more) memory and 1 Gb (better a few Gb)
of disk space. For coarse grain parallization on the k-point level, a cluster of PCs with a 100 Mb/s
network is sufficient. Faster communication is recommended for the fine grain (single k-point)
parallel version.
In order to use all options and features (such as the new graphical user interface w2web or some
of its plotting tools) the following public domain program packages in addition to a F90 compiler
must be installed:
Usually these packages should be available on modern systems. If one of these packages is not
available, it can either be installed from public domain sources (see Chapt. 11) or the corresponding
configuration may be changed (e.g. using vi instead of emacs). None of the principal components
of WIEN2k requires these packages, only for advanced features or w2web they are needed.
WIEN2k has the following features that are new with respect to WIEN97:
1
2 CHAPTER 1. INTRODUCTION
I due to the new APW+lo basis set it is significantly faster (up to an order of magnitude).
Optimizations in the most time consuming parts of LAPW1 and LAPW2 have been made.
I iterative diagonalization (for cases with large matrices and few eigenvalues)
I beside the k-point parallelization (including heterogeneous workstation clusters) a fine grain
parallelization based on MPI is also available.
I A new web-based graphical user interface w2web has been developed. It does NOT require
an X-environment and thus WIEN2k can be controlled from (but not run on !) any Windows-
PC. This should particularly help the novice to get acquainted with WIEN2k but it should be
useful for the regular user as well.
I support for AFM and FSM calculations
I spin-orbit coupling, including a new p1/2 -LO for higher accuracy
I wavefunction plotting
I determination of irreducible representations
I elastic constants (cubic cases only)
I Topological analysis based on Bader’s “atoms in molecules” concept
I LDA+U, orbital polarization (OP), magnetic and electric fields
I Exact-exchange and Hybrid functionals inside spheres
I new PKZB and TPSS meta-GGA functionals
The development of WIEN2k was made possible by support from many sources. We try to give
credit to all who have contributed. We hope not to have forgotten anyone who made an important
contribution for the development or the improvement of the WIEN2k code. If we did, please let us
know (we apologize and will correct it). The main developers in addition to the authors are the
following groups:
We want to thank those WIEN97 users, who reported bugs or made suggestions and thus con-
tributed to new versions as well as persons who have made major contributions in the develop-
ment of previous versions of the code:
pblaha@[Link]
kschwarz@[Link]
4 CHAPTER 1. INTRODUCTION
Part I
5
2 The basic concepts of the present
band theory approach
with EN N the repulsive Coulomb energy of the fixed nuclei and the electronic contributions, la-
belled conventionally as, respectively, the kinetic energy (of the non-interacting particles), the
electron-electron repulsion, nuclear-electron attraction, and exchange-correlation energies. Two
approximations comprise the LSDA, i), the assumption that Exc can be written in terms of a local
exchange-correlation energy density µxc times the total (spin-up plus spin-down) electron density
as
Z
Exc = µxc (ρ↑ , ρ↓ ) ∗ [ρ↑ + ρ ↓]dr (2.1)
and ii), the particular form chosen for that µxc . Several forms exist in literature, we use the most
recent and accurate fit to the Monte-Carlo simulations of Ceperly and Alder by Perdew and Wang
92. Etot has a variational equivalent with the familiar Rayleigh-Ritz principle. The most effective
way known to minimize Etot by means of the variational principle is to introduce orbitals χσik
constrained to construct the spin densities as
X
ρσ (r) = ρσik |χσik (r)|2 (2.2)
i,k
Here, the ρσik are occupation numbers such that 0 ≤ ρσik ≤ 1/wk , where wk is the symmetry-required
weight of point k. Then variation of Etot gives the Kohn-Sham equations (in Ry atomic units),
which must be solved and thus constitute the primary computational task. This Kohn-Sham equa-
tions must be solved self-consistently in an iterative process, since finding the Kohn-Sham orbitals
requires the knowledge of the potentials which themselves depend on the (spin-) density and thus
on the orbitals again.
7
8 CHAPTER 2. BASIC CONCEPTS
Recent progress has been made going beyond the LSDA by adding gradient terms of the electron
density to the exchange-correlation energy or its corresponding potential. This has led to the gen-
eralized gradient approximation (GGA) in various parameterizations, e.g. the one by Perdew et al
92 or Perdew, Burke and Ernzerhof (PBE) 96, which is the recommended option.
A recent version called meta-GGA by Perdew et al (1999) and Tao et al. (2003) employes for the
evaluation of the exchange-correlation energy not only the gradient of the density, but also the
kinetic energy density τ (r). Unfortunately, such schemes are not yet self-consistent.
Recently, the development of the Augmented Plane Wave (APW) methods from Slater’s APW, to
LAPW and the new APW+lo was described by Schwarz et al. 2001.
The linearized augmented plane wave (LAPW) method is among the most accurate methods for
performing electronic structure calculations for crystals. It is based on the density functional theory
for the treatment of exchange and correlation and uses e.g. the local spin density approximation
(LSDA). Several forms of LSDA potentials exist in the literature , but recent improvements using
the generalized gradient approximation (GGA) are available too (see sec. 2.1). For valence states
relativistic effects can be included either in a scalar relativistic treatment (Koelling and Harmon 77)
or with the second variational method including spin-orbit coupling (Macdonald 80, Novák 97).
Core states are treated fully relativistically (Desclaux 69).
A description of this method to linearize Slater’s old APW method (i.e. the LAPW formalism) and
further programming hints are found in many references: Andersen 73, 75, Koelling 72, Koelling
and Arbman 75, Wimmer et al. 81, Weinert 81, Weinert et al. 82, Blaha and Schwarz 83, Blaha et al.
85, Wei et al. 85, Mattheiss and Hamann 86, Jansen and Freeman 84, Schwarz and Blaha 96). An
excellent book by D. Singh (Singh 94) describes all the details of the LAPW method and is highly
recommended to the interested reader. Here only the basic ideas are summarized; details are left
to those references.
Like most “energy-band methods“, the LAPW method is a procedure for solving the Kohn-Sham
equations for the ground state density, total energy, and (Kohn-Sham) eigenvalues (energy bands)
of a many-electron system (here a crystal) by introducing a basis set which is especially adapted to
the problem.
I
I
II
Figure 2.1: Partitioning of the unit cell into atomic spheres (I) and an interstitial region (II)
This adaptation is achieved by dividing the unit cell into (I) non-overlapping atomic spheres (cen-
tered at the atomic sites) and (II) an interstitial region. In the two types of regions different basis
sets are used:
2.2. THE APW METHODS 9
1. (I) inside atomic sphere t, of radius Rt , a linear combination of radial functions times spheri-
cal harmonics XYlm (r) is used (we omit the index t when it is clear from the context)
φkn = [Alm,kn ul (r, El ) + Blm,kn u̇l (r, El )]Ylm (r̂) (2.4)
lm
where ul (r, El ) is the (at the origin) regular solution of the radial Schroedinger equation for
energy El (chosen normally at the center of the corresponding band with l-like character)
and the spherical part of the potential inside sphere t; u̇l (r, El ) is the energy derivative of
ul evaluated at the same energy El . A linear combination of these two functions constitute
the linearization of the radial function; the coefficients Alm and Blm are functions of kn (see
below) determined by requiring that this basis function matches (in value and slope) each
plane wave (PW) the corresponding basis function of the interstitial region; ul and u̇l are
obtained by numerical integration of the radial Schroedinger equation on a radial mesh
inside the sphere.
The solutions to the Kohn-Sham equations are expanded in this combined basis set of LAPW’s
according to the linear variation method
X
ψk = cn φkn (2.6)
n
and the coefficients cn are determined by the Rayleigh-Ritz variational principle. The convergence
of this basis set is controlled by a cutoff parameter Rmt Kmax = 6 - 9, where Rmt is the smallest
atomic sphere radius in the unit cell and Kmax is the magnitude of the largest K vector in equation
(2.6).
In order to improve upon the linearization (i.e. to increase the flexibility of the basis) and to make
possible a consistent treatment of semicore and valence states in one energy window (to ensure
orthogonality) additional (kn independent) basis functions can be added. They are called “local
orbitals (LO)“ (Singh 91) and consist of a linear combination of 2 radial functions at 2 different
energies (e.g. at the 3s and 4s energy) and one energy derivative (at one of these energies):
φLO
lm = [Alm ul (r, E1,l ) + Blm u̇l (r, E1,l ) + Clm ul (r, E2,l )]Ylm (r̂) (2.7)
The coefficients Alm , Blm and Clm are determined by the requirements that φLO should be normal-
ized and has zero value and slope at the sphere boundary.
φlo
lm = [Alm ul (r, E1,l ) + Blm u̇l (r, E1,l )]Ylm (r̂) (2.9)
This new lo (denoted with lower case to distinguish it from the LO given in equ. 2.7) looks almost
like the old “LAPW”-basis set, but here the Alm and Blm do not depend on kn and are determined
by the requirement that the lo is zero at the sphere boundary and normalized.
Thus we construct basis functions that have “kinks” at the sphere boundary, which makes it nec-
essary to include surface terms in the kinetic energy part of the Hamiltonian. Note, however, that
the total wavefunction is of course smooth and differentiable.
As shown by Madsen et al. (2001) this new scheme converges practically to identical results as the
LAPW method, but allows to reduce “RKmax” by about one, leading to significantly smaller basis
sets (up to 50 %) and thus the corresponding computational time is drastically reduced (up to an
order of magnitude). Within one calculation a mixed “LAPW and APW+lo” basis can be used for
different atoms and even different l-values for the same atom (Madsen et al. 2001). In general one
describes by APW+lo those orbitals which converge most slowly with the number of PWs (such
as TM 3d states) or the atoms with a small sphere size, but the rest with ordinary LAPWs. One
can also add a second LO at a different energy so that both, semicore and valence states, can be
described simultaneously.
P
VLM (r)YLM (r̂) inside sphere
LM
V (r) = (2.10)
VK eiK·r
P
outside sphere
K
and the charge densities analogously. Thus no shape approximations are made, a procedure fre-
quently called a “full-potential“ method.
The “muffin-tin“ approximation used in early band calculations corresponds to retaining only the
l = 0 component in the first expression of equ. 2.10 and only the K = 0 component in the second.
This (much older) procedure corresponds to taking the spherical average inside the spheres and
the volume average in the interstitial region.
The total energy is computed according to Weinert et al. 82.
Rydberg atomic units are used except internally in the atomic-like programs (LSTART and LCORE)
or in subroutine outwin (LAPW1, LAPW2), where Hartree units are used. The output is always
given in Rydberg units.
The forces at the atoms are calculated according to Yu et al (91). For the implementation of this
formalism in WIEN see Kohler et al (96) and Madsen et al. 2001. An alternative formulation by
Soler and Williams (89) has also been tested and found to be equivalent, both in computationally
efficiency and numerical accuracy (Krimmel et al 94).
The Fermi energy and the weights of each band state can be calculated using a modified tetrahe-
dron method (Blöchl et al. 94), a Gaussian or a temperature broadening scheme.
Spin-orbit interactions can be considered via a second variational step using the scalar-relativistic
eigenfunctions as basis (see Macdonald 80, Singh 94 and Novák 97). In order to overcome the prob-
lems due to the missing p1/2 radial basis function in the scalar-relativistic basis (which corresponds
to p3/2 ), we have recently extended the standard LAPW basis by an additional “p1/2 -local orbital”,
i.e. a LO with a p1/2 basis function, which is added in the second-variational SO calculation (Kuneš
et al. 2001).
2.2. THE APW METHODS 11
It is well known that for localized electrons (like the 4f states in lanthanides or 3d states in some
TM-oxides) the LDA (GGA) method is not accurate enough for a proper description. Thus we have
implemented various forms of the LDA+U method as well as the “Orbital polarization method”
(OP) (see Novák 2001 and references therein). In addition you can also calculate exact-exchange
inside the spheres and apply various hybrid functionals (see Tran et al. 2006 for details).
One can also consider interactions with an external magnetic (see Novák 2001) or electric field (via
a supercell approach, see Stahn et al. 2000).
PROPERTIES:
The density of states (DOS) can be calculated using the modified tetrahedron method of Blöchl et
al. 94.
X-ray absorption and emission spectra are determined using Fermi’s golden rule and dipole matrix
elements (between a core and valence or conduction band state respectively). (Neckel et al. 75,
Schwarz et al 79,80)
X-ray structure factors are obtained by Fourier Transformation of the charge density.
Optical properties are obtained using the “Joint density of states” modified with the respective
dipole matrix elements according to Ambrosch et al. 95, Abt et al. 94, Abt 97. and in particular
Ambrosch 06. A Kramers-Kronig transformation is also possible.
An analysis of the electron density according to Bader’s “atoms in molecules” theory can be made
using a program by J. Sofo and J. Fuhr (2001)
12 CHAPTER 2. BASIC CONCEPTS
3 Quick Start
Contents
3.1 Naming conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Starting the server . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Connecting to the w2web server . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Creating a new session . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5 Creating a new case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.6 Creating the struct file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.7 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.8 The SCF calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.9 The [Link] file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.10 Saving a calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.11 Calculating properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.12 Setting up a new case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
We assume that WIEN2k is properly installed and configured for your site and that you ran
userconfig lapw to adjust your path and environment. (For a detailed description of the in-
stallation see chapter 11.
This chapter is intended to guide the novice user in the handling of the program package. We
will use the example of TiC in the sodium chloride structure to show which steps are necessary to
initialize a calculation and run a self consistent field cycle. We also demonstrate how to calculate
various physical properties from these SCF data. Along the way we will give all important infor-
mation in a very abridged form, so that the novice user is not flooded with information, and the
experienced user will be directed to more complete information.
In this chapter we will also show, how the new graphical user interface w2web can be utilized to
setup and run the calculations.
Before we begin with our introductory example, we describe the naming conventions, to which we
will adhere throughout this user’s guide.
On UNIX systems the files are specified by [Link] and it is required that all files reside in a
subdirectory ./case. Here and in the following sections and in the shell scripts which run the
package themselves, we follow a simple, systematic convention for file labeling.
For the general discussion (when no specific crystal is involved), we use case, while for a specific
case, e.g. TiC, we use the following notation:
13
14 CHAPTER 3. QUICK START
Figure 3.1: TiC in the sodium chloride structure. This plot was generated using BALSAC (see
9.23.1). Interface programs between WIEN2k and BALSAC are available.
case=TiC
The filetype “type” always describes the content of the file (e.g.,
[Link]
Start the user interface w2web on the computer where you want to execute WIEN2k(you may have
to telnet, ssh,.. to this machine) with the command
If the default port (7890) used to serve the interface is already in use by some other process,
you will get the error message w2web failed to bind port 7890 - port already in
use!. Then you will have to choose a different port number (between 1024 and 65536) . Please
remember this port number, you need it when connecting to the w2web server.
Note: Only user root can specify port numbers below 1024!
At the first startup of this server, you will also be asked to setup a username and password, which
is required to connect to this server.
3.3. CONNECTING TO THE W2WEB SERVER 15
Use your favorite WWW-browser to connect to w2web, specifying the correct portnumber, e.g.
(If you do not remember the portnumber, you can find it by using “ps -ef | grep w2web” on the
computer where w2web is running.) You should see a screen as in Fig.3.2.
The user interface w2web uses sessions to distinguish between different working environments
and to quickly change between different calculations. First you have to create a new session (or
select an old one). Enter “TiC” and click the “Create” button.
Note: Creating a session does not automatically create a new directory!
You will be placed in your home directory if no working directory was designated to this session
previously (or if the directory does not exist any more).
Title TiC
Lattice F (for face centered)
a 4.328 Å(make sure the Ang button is selected)
b 4.328 Å
c 4.328 Å
α, β, γ 90
Atom Ti, enter position (0,0,0)
Atom C, enter position (.5,.5,.5)
Click “Save Structure” (Z will be updated automatically) and “set automatically RMT and con-
tinue editing ”:
3.6. CREATING THE STRUCT FILE 17
This will compute the nearest neigbor distances using the program nn and setrmt lapw will then
determine the optimal RMT values (muffin-tin radius, atomic sphere radius). To learn more about
the philosophy of setting RMTs see http : //[Link]/regu ser/f aq. Since it is essential to keep
RMTs constant within a series of calculations (eg. when you do a Volume-optimization, see 3.11.6
), you should already now decide whether you want to do just one single calculation with fixed
structural parameters, or whether you intend a relaxation of internal parameters (using forces and
min lapw) or a volume optimization, which would required reduced RMT values.
Choose a reduction of 3 % so that we can later optimize the lattice parameter.
When you are done, exit the StructGen with “save file and clean up”. This will generate the file
[Link] (shown now in view-only mode with a different background color), which is the
master input file for all subsequent programs.
A few other hints on StructGen:
18 CHAPTER 3. QUICK START
You have to click on Save Structure after every modifications you make in the white fields.
Add/remove a position/atom only if you have made no other changes before.
In a face-centered (body-centered) spacegroup you have to enter just one atom (not the ones in
(.5,.5,0),. . . ).
StructGen offers a built in calculator: Each position of equivalent atoms can be entered as a num-
ber, a fraction (e.g. 1/3) or a simple expression (e.g. 0.21 + 1/3). The first position defines the
variables x, y and z, which can be using in expression defining the other positions (e.g. −y, x,
−z + 1/2).
When you now choose “Files o show all files”, you will see, that [Link] has been created.
For a detailed description of these files consult sections 4.3 and 6.4.3.
After the two basic input files have been created, initalization of the calculation is done by “Execu-
tion o initialize calc.”. This will guide you through the steps necessary to initialite the calculation.
Simply follow the steps that are highlighted in green and follow the instructions.
The initialization process is described in detail in section 5.1.2.
Alternatively you could run the script init lapw from the command line. All actions of this script
are logged in short in :log and in detail in the file [Link], which can easily be accessed
by Utils. o show dayfile.
Initializing the calculation will run several steps automatically, where x is the script to start WIEN2k
programs (see section: 5.1.1).
x nn calculates the nearest neighbors up to a specified distance and thus helps to determine the
atomic sphere radii (you must specify a distance factor f, e.g. 2, and all distances up to f *
NN-dist. are calculated)
view [Link] : check for overlapping spheres, coordination numbers and nearest neighbor
distances, (e.g. in the sodium chloride structure there must 6 nearest and 12 next nearest
neighbors). Using these distances and coordinations you can check whether you put the
proper positions into your struct file or if you made a mistake. nn also checks whether
your equivalent atoms are really crystallographically equivalent and eventually writes a new
struct-file which you may or may not accept. If you have not done so at the very begin-
ning, go back to StructGen and choose proper RMT values. You can save a lot of CPU-time by
changing RMT to almost touching spheres. See Sec.4.3
x sgroup calculates the point and spacegroups for the given structure
view [Link] : Now you can either accept the [Link] file generated by sgroup
(if you want to use the spacegroup information or a different cell has been found by sgroup)
or keep your original file (default).
x symmetry generates from a raw [Link] file the space group symmetry operations, de-
termines the point group of the individual atomic sites, generates the LM expansion for the
lattice harmonics (in case.in2 st) and local rotation matrices (in [Link] st).
view [Link] : check the symmetry operations (they have been written to or compared with
already available ones in [Link] by the program symmetry) and the point group sym-
metry of the atoms (You may compare them with the “International Tables for X-Ray Crys-
tallography“). If the output does not match your expectations from the “Tables”, you might
have made an error in specifying the positions. The [Link] file will be updated with
symmetry operations, positive or negativ atomic counter (for “cubic” point group symme-
tries) and the local rotation matrix.
3.7. INITIALIZATION 19
instgen lapw : You are requested to generate an input file [Link] and can define the spin-
polarization of each atom. While this is not important for TiC, it is very important for spin-
polarized calculations and in particular for anti-ferromagnetic cases, where you should “flip”
the spin of the AFM atoms and/or set the spin of the “non-magnetic” atoms (eg. oxygen in
NiO) to zero.
x lstart generates atomic densities (see section 6.4) and determines how the orbitals are treated in
the band structure calculations (i.e. as core or band states, with or without local orbitals, . . . ).
You are requested to specify the desired exchange correlation potential and an energy that
separates valence from core states. For TiC select the recommended potential option “GGA
of Perdew-Burke-Ernzerhof 96” and a separation energy of -6.0 Ry.
edit [Link] : check the output (did you specify a proper atomic configuration, did lstart
converge, are the core electrons confined to the atomic sphere?). Warnings for the radial
mesh can usually be neglected since it affects only the atomic total energy. lstart generates
TiC.in0 st, in1 st, in2 st, inc st and inm st. For Ti it selects automatically 1s, 2s,
and 2p as core states, 3s and 3p will be treated with local orbitals together with 3d, 4s and 4p
valence states.
edit TiC.in1 st : As mentioned, the input files are generated automatically with some default val-
ues which should be a reasonable choice for most cases. Nevertheless we highly recommend
that you go through these inputs and become familiar with them. The most important param-
eter here is RKMAX, which determines the number of basis functions (size of the matrices).
Values between 5-9 (APW) and 6-10 (LAPW) are usually reasonable. You may change here
the usage of APW or LAPW (set 1 or 0 after the CONT/STOP switch), since often APW is
necessary only for orbitals more difficult to converge (3d, 4f). Here we will just change EMAX
of the energy window from 1.5 to 2.0 Ry in order to be able to calculate the unoccupied DOS to higher
energies.
edit TiC.in2 st : Here you may limit the LM expansion (for some speedup), change the value of
GMAX (in cases with small spheres (e.g. systems with H-atoms) values of 15-24 are recom-
mended) or specify a different BZ-integration method to determine the Fermi energy. For
this example you should not change anything so that you can compare your results with the
test run.
Copy all generated inputs (from [Link]∗ st to [Link]*). In cases without inversion sym-
metry the files case.in1c, in2c are produced.
x kgen generates a k-mesh in the Brillouin zone (BZ). You must specify the number of k-points
in the whole BZ (use 1000 for comparison with the provided output, a “good” calculations
needs many more). For details see section 6.5.
view [Link] : check the number of k-points in the irreducible wedge of the BZ (IBZ) and the
energy interval specified for the first k-point. You can now either rerun kgen (and generate
a different k-mesh) or continue.
x dstart generates a starting density for the SCF cycle by superposition of atomic densities gener-
ated in lstart. For details see section 6.6.
view [Link] (check if gmax >gmin)
Now you are asked , whether or not you want to run a spin-polarized calculation (in such a case
case dstart is re-run to generate spin-densities). For TiC say No.
Alternatively, w2web provides an “expert-mode”, where some inputs can be specified right at the
beginning and then the whole initialization runs at once. Please check carefully the STDOUT-
listing and some output-files for possible errors or warnings!!
Initialization of a calculation (running init lapw) will create all inputs for the subsequent SCF
calculation choosing some default options and values. You can find a list of input files using “Files
o input files” ( 3.5).
20 CHAPTER 3. QUICK START
After the case has been set up, a link to “run SCF” is added, (“Run Programs o run SCF” and you
should invoke the self-consistency cycle (SCF). This runs the script run lapw with the desired
options.
The SCF cycle consists of the following steps:
After selecting “run SCF” from the “Execution” menu, the SCF-window will open, and you can
now specify additional parameters. For this example we select charge convergence to 0.0001: Spec-
ify “charge” to be used as convergence criterion, and select a value of 0.0001 (-cc 0.0001).
To run the SCF cycle, click on “Run!”
Since this might take a long time for larger systems; you can specify the “Execution type” to be batch
or submit (if your system is configured with a queuing system and w2web has been properly set
up, see section 11.3).
While the calculation is running (as indicated by the status frame in the top right corner of the
window), you can monitor several quantities (see section 3.9).
Once the calculation is finished (11 iterations), view [Link] for timing and errors and
compare your results with the files in the provided example (TiC/case scf).
3.9. THE [Link] FILE 21
For magnetic systems you would run a spin-polarized calculation with the script runsp lapw.
The program flow of such a calculation is described in section 4.5.2 and the script itself in section
5.1.3.
Then just execute another run lapw using “Execution o run SCF”.
I The total charge density includes the Ti 3s and 3p states and the resulting density
around Ti would be very large and dominated by these semicore states. To get
a “meaningful” picture of the chemical bonding effects one must remove these
states. Inspection of TiC.scf1 and TiC.scf2 should allow you to select an
EMIN value to eliminate the Ti 3s and 3p semicore states.
22 CHAPTER 3. QUICK START
ORTHO
Compare the result with the electron density plotted in the (100) plane (see figure 3.8). The pro-
gram gnuplot (public domain) must be installed on your computer. For more advanced graphics
use your favorite plotting package or specify other options in gnuplot (see rhoplot lapw how
gnuplot is called).
24 CHAPTER 3. QUICK START
Select “Density of States (DOS)” from the “Tasks” menu and click on the buttons one by one:
I Calculate partial charges (x lapw2 -qtl). (This is only possible, when you still
have a valid [Link] file on a tetrahedral mesh.)
I Create [Link], either using “configure [Link]” or/and by “editing” the offered
template input file. Select: total DOS, Ti-d, Ti-deg , Ti-dt2g , C-s and C-p-like DOS.
TiC
-0.50 0.00200 1.500 0.003 EMIN, DE, EMAX, Gauss-broadening
6 NUMBER OF DOS-CASES
0 1 tot (atom,case,description)
1 4 Ti d
1 5 Ti eg
1 6 Ti t2g
2 2 C s
2 3 C p
For a detailed description of input options consult section 8.1.3
I Calculate DOS (x tetra).
I Preview output using “dosplot”
If you want to use the supplied plotting interface dosplot2 to preview the results, the program
gnuplot (public domain) must be installed on your computer.
The calculated DOS can be compared with figures 3.9 and 3.10. Together with the electron density
the partial DOS allows you to analyse the chemical bonding (covalency between T i−deg and C −p,
non-bonding T i − dt2g , charge transfer estimates,....)
3.11. CALCULATING PROPERTIES 25
Select “X-Ray Spectra” from the “Tasks menu” and click on the buttons one by one:
I Calculate partial charges (x lapw2 -qtl). This is only possible, when you still
have a valid [Link] file on a tetrahedral mesh. To reproduce this figure you
will have to increase the EMAX value in your TiC.in1 to 2.5 Ry and rerun x lapw1
I Edit [Link]; choose the offered template. This template will calculate the LIII -
spectrum of the first atom (Ti in this example) in the energy range between -2 and
15 eV. For a detailed description of the contents of this input file refer to section
8.10.3.
I Calculate spectra
I Preview spectra
If you want to use the supplied plotting interface specplot to preview the results, the public domain
program gnuplot must be installed on your computer. The calculated TiC Ti-LIII -spectrum can be
compared with figure 3.11.
3.11.4 Bandstructure
Select “Bandstructure” from the “Tasks” menu and click on the buttons one by one:
I Edit [Link]: insert the correct Fermi energy (which can be found in the saved
scf-file) and specify plotting parameters. For comparison with figure 3.12 select
an energy-range from -13 to 8 eV.
I Calculate Bandstructure (x spaghetti).
I Preview Bandstructure (needs ghostscript installed).
If you want to preview the bandstructure, the program ghostview (public domain) must be in-
stalled on your computer. You can compare your calculated bandstructure with figure 3.12.
tic atom 0 size 0.40
8.0
7.0
6.0
5.0
4.0
3.0
2.0
1.0
0.0 EF
-1.0
Energy (eV)
-2.0
-3.0
-4.0
-5.0
-6.0
-7.0
-8.0
-9.0
-10.0
-11.0
-12.0
-13.0
W L Λ Γ ∆ X Z W K
8.0
7.0
6.0
5.0
4.0
3.0
2.0
1.0
0.0 EF
-1.0
Energy (eV)
-2.0
-3.0
-4.0
-5.0
-6.0
-7.0
-8.0
-9.0
-10.0
-11.0
-12.0
-13.0
W L Λ Γ ∆ X Z W K
Figure 3.13: Bandstructure of TiC, showing t2g-character bands of Ti in character plotting mode
Select “Optimize (V,c/a)” from the “Execution” menu. Setup the shell script [Link]
script using x optimize and volume variations of -10, -5, 0, +5 and +10%. Then run the
[Link]. When the job has finished, you should click on Plot and then preview the en-
ergy curve.
You should get an energy curve as in figure 3.14. On the screen you will find the fitting parameters
for the “equation of states” (Murnaghan, Birch-Murnaghan and the EOS2 equation, see sec. 9.12).
This information is also written to [Link].
Usually this information can be collected from the “International Tables of Crystallography” once
you know the space group, the Wyckoff position and the internal free coordinates.
Usually for a new “case“ the input is not created from scratch, but one uses the struct file from a
similar case as pattern. Change into the lapw subdirectory and proceed as follows:
Now edit case [Link] (see section 4.3) as necessary (Note: this is a fixed formatted file,
so all values must remain at their proper columns). Afterwards generate case [Link] using
instgen lapw.
Use the menu Session Mgmt. o change session of w2web to create a new session (enter the name of
the new session and click on “Create”). Then you should also create a new directory and “select”
it..
When you select “Execution o StructGen”, you have several choices:
You can just specify the number of non-equivalent atoms and a template file will be created. In
StructGen you simply specify the lattice (type or spacegroup), cell parameters and name and po-
sitions of atoms. When you “save file and clean up” the new [Link] file and the [Link]
file are created automatically.
Alternatively, you can use cif2struct or xyz2struct to convert a “cif”, “txt” or “xyz” file
into the WIEN2k [Link] file. Check page 168 for more info on the specific file formats.
For more information on the StructGen refer to page 169.
30 CHAPTER 3. QUICK START
Part II
31
4 File structure and program flow
Contents
4.1 Flow of input and output files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.2 Input/Output files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 The [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.4 The [Link] file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.5 Flow of programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
Each program is started with (at least) one command line argument, e.g.
programX [Link]
in which the arguments specifies a filename, in which FORTRAN I/O units are connected to unix
filenames. (See examples at specific programs). These “def“-files are generated automatically
when the standard WIEN2k scripts x, init lapw or run lapw are used, but may be tailored by
hand for special applications. Using the option
x program -d
a def-file can be created without running the program. In addition each program reads/writes the
following files:
The programs of the SCF cycle (see figure 4.1) write the following files:
[Link] a file containing only the most significant output (see description below).
[Link] error report file, should be empty after successful completion of a program (see
chapter 6)
33
34 CHAPTER 4. FILES AND PROGRAM FLOW
clmsum
broyd2
broyd1
(new)
(new)
scfm
MIXER
scf
MERGE
clmval
help3*
clmsc
scfm
(old)
scf2
LAPW2
clmcor
vector
kgen
scf1
scfc
LAPW1
LCORE
scf0
vsp
vns
LAPW0
necessary
optional
clmsum
(old)
Figure 4.1: Data flow during a SCF cycle ([Link], [Link], [Link], [Link] and
optional files are omitted)
4.1. FLOW OF INPUT AND OUTPUT FILES 35
The following tables describe input and output files for the initialization programs nn, sgroup,
symmetry, lstart, kgen, dstart (table 4.1), the utility programs tetra, irrep, spaghetti,
aim, lapw7, elnes, lapw3, lapw5, xspec, optic, joint, kram, optimize and mini (table
4.2) as well as for a SCF cycle of a non-spin-polarized case (table 4.2). Optional input and output
files are used only if present in the respective case subdirectory or requested/generated by an
input switch. The connection between FORTRAN units and filenames are defined in the respective
[Link] files. The data flow is illustrated in Fig. 4.1.
[Link]
[Link]
[Link]
[Link]
[Link] diag
KRAM [Link] [Link] [Link]
[Link] [Link] [Link]
[Link]
OPTIMIZE [Link] case [Link] [Link] case vol [Link]
case c/a [Link]
MINI [Link] [Link] mini [Link] [Link] inter
[Link] [Link] case.tmpM1
[Link] [Link] case.struct1
[Link] [Link] [Link] mini1
[Link] .min hess .minrestart
IRREP [Link] [Link]
[Link] [Link]
AIM [Link] [Link] [Link]
[Link] [Link]
[Link]
LAPW7 [Link] case.output7 [Link]
[Link] [Link]
case.in7 [Link]
[Link]
QTL [Link] [Link]
[Link] [Link]
[Link]
[Link]
case.scf2
LCORE [Link] [Link] [Link] [Link]
[Link] [Link]
[Link] [Link]
[Link] [Link]
After LCORE the [Link] files are appended to [Link] and the
[Link] file is renamed to [Link] old (see run lapw)
MIXER [Link] [Link] old [Link] [Link]*
[Link] [Link] [Link]
[Link] [Link] [Link]
[Link] [Link] [Link]
case.broyd1
case.broyd2
After MIXER the file [Link] is appended to [Link], so that after an iteration is
completed, the two essential files are [Link] and [Link].
Table 4.3: Input and output files of main programs in an SCF cycle
[Link] Contains the Alm , Blm , Clm coefficients of the wavefunctions (generated optional by
lapw2).
[Link] Contains the charge density of previous iterations if you use Broyden’s method for
mixing. They are removed when using save lapw. They should be removed by hand when
calculational parameters (RKMAX, kmesh, . . . ) have been changed, or the calculation crashed
due to a too large mixing and are restarted by using a new density generated by dstart.
[Link] Contains the core charge density (as σ(r) = 4πr2 ρ(r) and has only a spherical part).
In spin-polarized calculations two files [Link] and [Link] are used instead.
[Link] Contains the semi-core charge density in a 2-window calculation, which is no longer
recommended. In spin-polarized calculations two files are used instead: [Link] and
[Link].
[Link] Contains the total charge density in the lattice harmonics representation and as
Fourier coefficients. (The LM=0,0 term is given as σ(r) = 4πr2 ρ(r), the others as r2 ρLM (r);
suitable for generating electron density plots using lapw5 when the TOT-switch is set,
(see section 8.6). In spin-polarized calculations two additional files [Link] and
[Link] contain the spin densities. Generated by dstart or mixer.
[Link] Contains the valence charge density as r2 ρLM (r); suitable for generating valence elec-
tron density plots using lapw5 when the VAL-switch is set, (see 8.6). In spin-polarized cal-
culations two files [Link] and [Link] are used instead.
[Link]/dn Contains the density matrix generated by lapwdm for LDA+U, OP or Hybrid-
DFT calculations.
[Link] Contains the density of states (states/Ry) and corresponding energy (in Ry at the inter-
nal energy scale) generated by tetra. X can be 1-3. Additional files [Link] contain
the DOS in (states/eV) and the energy in eV with respect to EF.
case.help03X Contains eigenvalues and partial charges for atom number X.
[Link] This file contains the indices of the tetrahedra in terms of the list of k-points. It is used
in lapw2 (if EFMOD switch in case.in2 is set to TETRA, see 7.5.3) and in tetra.
[Link] This file contains a list of k-points in the first BZ and represents a tetrahedral (special
point) mesh. It is generated in kgen and can either be inserted into the case.in1 file or
used directly in kgen.
[Link] Contains eigenvalues and corresponding partial charges (bandwise) in a form suitable for
tetra and band structure plots with “band character”. The decomposition of these charges
is controlled by ISPLIT in [Link].
[Link] Contains the radial basis functions inside spheres (generated optional by lapw2).
[Link] Contains the electron densities on a grid in a specified plane generated by lapw5. This
file can be used as input for your favorite contour or 3D plotting program.
38 CHAPTER 4. FILES AND PROGRAM FLOW
[Link] Contains the atomic densities generated by lstart. They are used by dstart to gener-
ate a first crystalline density ([Link]).
case.r2v Contains the exchange potential (in the lattice harmonics representation as r2 ∗ VLM (r)
and as Fourier coefficients) in a form suitable for plotting with lapw5.
[Link] mini Contains the last scf-iteration of each individual time (geometry) step during a struc-
tural minimization using mini. Thus this file contains a complete history of properties (en-
ergy, forces, positions) during a structural minimization.
[Link] Contains the atomic densities for those states with a “P” in [Link]. Generated in
lstart and used for difference densities in lapw5.
[Link] ps A ps file with the energy bandstructure plot generated by spaghetti.
[Link] A xmgrace file with the energy bandstructure plot generated by spaghetti.
[Link] Contains the Coulomb potential (in the lattice harmonics representation as r2 ∗ VLM (r)
and as Fourier coefficients) in a form suitable for plotting with lapw5.
[Link] Contains the orbital potential (in Ry) generated by orb for LDA+U or hybrid-DFT cal-
culations in form of a (2l+1,2l+1) matrix.
[Link] Contains the total potential (in the lattice harmonics representation as r2 ∗ VLM (r) and
as Fourier coefficients) in a form suitable for plotting with lapw5.
[Link] Binary file, contains the eigenvalues and eigenvectors of all k-points calculated in
lapw1. In spin-polarized calculations two files [Link] and [Link] are
used instead. lapwso generates [Link].
[Link] Contains the eigenvalues of all k-points calculated in lapw1. In spin-polarized calcu-
lations two files [Link] and [Link] are used instead. lapwso generates
[Link].
[Link] Contains the non-spherical part of the total potential V. Inside the sphere the radial co-
efficients of the lattice harmonics representation are listed (for L greater than 0), while for
the interstitial region the reanalyzed Fourier coefficients are given (see equ. (2.10)). In spin-
polarized calculations two files [Link] and [Link] are used instead.
[Link]/dn Contains the orbital dependent part of the potential in LDA+U, OP or Hybrid-
DFT calculations. Generated in orb, used in lapw1.
[Link] Contains the spherical part of the total potential V stored as r ∗ V (thus the first val-
ues should be close to −2 ∗ Z). In spin-polarized calculations two files [Link] and
[Link] are used instead.
If you are using the “Struct Generator” from the graphical user interface w2web, you don’t have to
bother with this file directly! However, the description of the fields of the input mask can be found
here.
Note: If you are changing this file manually, please note that this is a formatted file and the proper column
positions of the characters are important! Use REPLACE instead of DELETE and INSERT during edit!
We start the description of this file with an abridged example for rutile TiO2 (adding line numbers):
MULT= 2 ISPLIT= 8 6
ATOM -1: X= 0.5000000 Y= 0.5000000 Z= 0.5000000
Titanium NPT= 781 R0=.000022391 RMT=2.00000000 Z:22.0 7
LOCAL ROT MATRIX: -.7071068 0.7071068 0.0000000 8
0.7071068 0.7071068 0.0000000 9
0.0000000 0.0000000 1.0000000 10
ATOM -2: X= 0.3050000 Y= 0.3050000 Z= 0.0000000
MULT= 4 ISPLIT= 8
ATOM -2: X= 0.6950000 Y= 0.6950000 Z= 0.0000000
ATOM -2: X= 0.8050000 Y= 0.1950000 Z= 0.5000000
ATOM -2: X= 0.1950000 Y= 0.8050000 Z= 0.5000000
Oxygen NPT= 781 R0=.000017913 RMT=1.60000000 Z: 8.0
LOCAL ROT MATRIX: 0.0000000 -.7071068 0.7071068
0.0000000 0.7071068 0.7071068
1.0000000 0.0000000 0.0000000
16 SYMMETRY OPERATIONS: 11
1 0 0 0.00 12
0 1 0 0.00 13
0 0 1 0.00 14
1 15
1 0 0 0.00
0 1 0 0.00
0 0-1 0.00
2
........
15
0 1 0 0.50
-1 0 0 0.50
0 0 1 0.50
16
------------------ bottom of file ---------------------------
P all primitive lattices except hexagonal [a sin(γ) sin(β), a cos(γ) sin(β), cos(β)], [0, b sin(α), b
cos(α)], [0, 0, c]
F face-centered [a/2, b/2, 0], [a/2, 0, c/2], [0, b/2, c/2]
B body-centered [a/2, -b/2, c/2],[a/2, b/2, -c/2], [-a/2, b/2, c/2]
CXY C-base-centered (orthorhombic only) [a/2, -b/2, 0], [a/2, b/2, 0], [0, 0, c]
CYZ A-base-centered (orthorhombic only) [a, 0, 0], [0, -b/2, c/2], [0, b/2, c/2]
CXZ B-base-centered (orthorh. and monoclinic [a sin(γ)/2, a cos(γ)/2, -c/2], [0, b, 0], [a sin(γ)/2, a
symmetry) cos(γ)/2,
√ c/2] √ √
R rhombohedral [a/
√ 3/2, -a/2, c/3],[a/ 3/2, a/2, c/3],[-a/ 3, 0, c/3]
H hexagonal [ 3a/2, -a/2, 0],[0, a, 0],[0, 0, c]
Table 4.4: Lattice type, description and bravais matrix used in WIEN2k
lattice type as defined in table 4.4. For definitions of the triclinic lattice see
SRC nn/dirlat.f
NAT number of inequivalent atoms in the unit cell
a, b, c unit cell parameters (in a.u., 1 a.u. = 0.529177 Å). In face- or body-centered
structures the non-primitive (cubic) lattice constant, for rhombohedral (R) lat-
tices the hexagonal lattice constants must be specified. (The following may help
you to convert between hexagonal and rhombohedral specifications:
π−αrhomb
ahex = 2cos(
q 2
)arhomb
chex = 3 a2rhomb − 13 a2hex
√
and (for fcc-like lattices) arhomb = acubic / 2
◦
α, β, γ angles between unit axis (if omitted, 90 is set as default). Set it only for P and
CXZ lattices
>>>: line 5 must now be repeated MULT-1 times for the other positions of each equivalent atom according
to the Wyckoff position in the “International Tables of Crystallography”.
line 7: format (A10,5X,I5,5X,F10.8,5X,F10.5,5X,F5.2)
name of atom, NPT, R0, RMT, Z
name of Use the chemical symbol. Positions 3-10 for further labeling of nonequivalent
atom atoms (use a number in position 3)
4.4. THE [Link] FILE 41
NPT number of radial mesh points (381 gives a good mesh for LDA calculations,
but for GGA twice as many points are recommended; always use an odd number
of mesh points!) the radial mesh is given on a logarithmic scale: r(n) = R0 ∗
e[(n−1)∗DX]
R0 first radial mesh point (typically between 0.0005 and 0.00005, smaller for heavy
elements, bigger for light ones; a struct-file generated by w2web will have
proper R0 values.)
RMT atomic sphere radius (muffin-tin radius), can easily be estimated after running
nn (see 6.1) and are set automatically with setrmt lapw see 5.2.6). The follow-
ing guidelines will be given here: Choose spheres as large as possible as this
will save MUCH computer time. But: Use identical radii within a series of cal-
culations (i.e. when you want to compare total energies) — therefore consider
first how close the atoms may possibly come later on (volume or geometry op-
timization); do NOT make the spheres too different (even when the geometry
would permit it), instead use the largest spheres for f-electron atoms, 10-20 %
smaller ones for d-elements and again 10-20 % smaller for sp-elements; H is a
special case, you may choose it much smaller (e.g. 0.6 and 1.2 for H and C) and
systems containing H need a much smaller RKMAX value (3-5) in case.in1.
Z atomic number
One can monitor the energy eigenvalues (listed for the first k-point only), the Fermi-energy or
the total energy. Often the electronic charges per atom reflect the convergence. Charge transfer
between the various atomic spheres is a typical process during the SCF cycles: large oscillations
should be avoided by using a smaller mixing parameter; monotonic changes in one direction sug-
gest a larger mixing parameter.
In spin-polarized calculations the magnetic moment per atomic site is an additional crucial quan-
tity which could be used as convergence criterion.
If a system has electric field gradients and one is interested in that quantity, one should monitor
the EFGs, because these are very sensitive quantities.
It is best to monitor several quantities, because often one quantity is converged, while another still
changes from iteration to iteration. The script run lapw has three different convergence criteria
built in, namely the total energy, the atomic forces and the charge distance (see 5.1.2, 5.1.3).
4.5. FLOW OF PROGRAMS 43
NN a program which lists the nearest neighbor distances up to a specified limit (defined by a
distance factor f) and thus helps to determine the atomic sphere radii. In addition it is a
very usefull additional check of your [Link] file (equivalency of atoms)
SGROUP determines the spacegroup of the structure defined in your [Link] file.
SYMMETRY generates from a raw [Link] file the space group symmetry operations, deter-
mines the point group of the individual atomic sites, generates the LM expansion for the
lattice harmonics and determines the local rotation matrices.
LSTART generates free atomic densities and determines how the different orbitals are treated in
the band structure calculations (i.e. as core or band states, with or without local orbitals,. . . ).
KGEN generates a k-mesh in the irreducible part of the BZ.
DSTART generates a starting density for the scf cycle by a superposition of atomic densities
generated in LSTART.
Then a self-consistency cycle is initiated and repeated until convergence criteria are met (see 3.8
and 5.1.3). This cycle can be invoked with a script run lapw, and consists of the following steps:
In many cases it is desirable to distinguish three types of electronic states, namely core, semi-core
and valence states. For example titanium has core (1s, 2s, 2p), semi-core (3s, 3p) and valence (3d,
4s, 4p) states. In our definition core states are only those whose charge is entirely confined inside
the corresponding atomic sphere. They are deep in energy, e.g., more than 7-10 Ry below the Fermi
energy. Semi-core states lie high enough in energy (between about 1 and 7 Ry below the Fermi
energy), so that their charge is no longer completely confined inside the atomic sphere, but has a
few percent outside the sphere. Valence states are energetically the highest (occupied) states and
always have a significant amount of charge outside the spheres.
44 CHAPTER 4. FILES AND PROGRAM FLOW
NN LSTART
check for atomic calculation
overlap. spheres H ψ nl = E nl ψ nl
DSTART
atomic densities
superposition of
SGROUP SYMMETRY input files
atomic densities
struct files struct files
input files ρ
KGEN
k−mesh
generation
LAPW0
2
VC = −8 πρ Poisson
VXC ( ρ) LDA
ORB
LDA+U, OP potentials V= VC + VXC
V VMT
LAPW1 LCORE
2 atomic calculation
− +V ψk = Ek ψk
H ψ nl = E nl ψ nl
Ek ψk ρcore Ecore
LAPWSO
add spin−orbit interaction
LAPW2
ρval = Σ ψ∗k ψ k
Ek < EF
ρval ρold
MIXER
LAPWDM
ρnew = ρold ( ρval + ρcore)
calculates density matrix
ρnew
yes
STOP converged ?
no
The energy cut-off specified in lstart during init lapw (usually -6.0 Ry) defines the separation
into core- and band-states (the latter contain both, semicore and valence). If a system has atoms
with semi-core states, then the best way to treat them is with “local orbitals“, an extension of the
usual LAPW basis. An input for such a basis set will be generated automatically. (Additional LOs
can also be used for valence states which have a strong variation of their radial wavefunctions with
energy (e.g. d states in TM compounds) to improve the quality of the basis set, i.e. to go beyond
the simple linearization).
For magnetic systems spin-polarized calculations can be performed. In such a case some steps are
done for spin-up and spin-down electrons separately and the script runsp lapw consists of the
following steps:
The use of spin-polarized calculations is illustrated for fcc Ni (section 10.2), one of the test cases
provided in the WIEN2k package.
Using the script runfsm lapw -m XX it is possible to constrain the total spin magnetic moment
per unit cell to a fixed value XX and thus force a particular ferromagnetic solution (which may
not correspond to the equillibrium). This is particularly useful for systems with several metastable
(non-) magnetic solutions, where conventional spin-polarized calculation would not converge or
the solution may depend on the starting density. Additional SO-interaction is not supported.
Please note, that once runfsm lapw has finished, only [Link] is ok, but
[Link] is NOT the proper up-spin vector and MUST NOT be used for the calculations
of QTLs (and DOS). It must be regenerated by x lapw1 -up (see also the comments for iterative
diagonalization in section 5.2.18).
Several considerations are necessary, when you want to perform an AFM calculation. Please have
also a look into $WIENROOT/SRC afminput/afminput test.
I You must construct a unit cell which allows for the desired AF ordering. For example for
bcc Cr you must select a “P” lattice and specify both atoms, Cr1 at (0,0,0) and Cr2 at (.5,.5,.5),
corresponding to a CsCl structure. Note, that it is important to label the two Cr atoms with
“Cr1” and “Cr2”, since only then the symmetry programs can detect that those atoms should
be different (although they have the same Z). If sgroup has interchanged some axis, try to undo
these changes, since afminput may not properly find the correct symmetry operations in such a case.
46 CHAPTER 4. FILES AND PROGRAM FLOW
I When you generate [Link] you must specify the correct magnetic order and flip the
spin of the AF atoms (i.e. invert the spin up and dn occupation numbers). In addition you
should set a zero moment (identical spin up and dn occupations) for all “non-magnetic”
atoms. This can be done conveniently using instgen lapw -ask or during “initialization”
using w2web.
I Now you can run either a “normal” spinpolarized initialization (without AFM option) and
runsp lapw or:
I Create a struct file of the non-magnetic (or ferro-magnetic) supergroup (run init lapw up to
lstart). Name it [Link] supergroup. (For example for bcc Cr, this would be a struct
file with the ordinary cubic lattice parameters, “B” type lattice and just one Cr at (0,0,0).)
I Run init lapw. At the end AFMINPUT creates an input file for the program CLMCOPY.
Depending on the presence of [Link] supergroup and the specific symmetry it
may/may not ask you to supply a symmetry operation/nonprimitive translation (see Sect.
9.5 .
I Run runafm lapw. This script calls LAPW1 and LAPW2 only for spin-up but the corre-
sponding spin-dn density is created by CLMCOPY according to the rules defined during
initialization. This reduces the required cpu time by a factor of 2 (and in addition the scf
cycle is much more stable).
I It is highly recommended that you save your work (save lapw) and check the results by
continuing with a regular runsp lapw. If nothing changes (E-tot and other properties), then
you are ok, otherwise make sure the scf calculation is well converged (-cc 0.0001 or better).
Eventually the system may not want to be antiferromagnetic (but for instance it is ferrimag-
netic!).
runafm lapw saves you more than a factor of 2 in in computer time, since only spin-up is cal-
culated and in addition the scf-convergence may be MUCH faster. It works also with LDA+U
([Link]/dn are also copied), but does NOT work with Hybrid-DFT nor spin-orbit cou-
pling, since this requires the presence of both vector files in the LAPWSO step.
You can add spin-orbit interaction in LAPWSO (called directly after LAPW1) using a second-
variational method with the scalar-relativistic orbitals (from LAPW1) as basis. The number of
eigenvalues will double since SO couples spin-up and dn states, so they are no longer separable.
In addition, LOs with a “p1/2 ” radial basis can be added. (Kunes et al. 2001)
To assist with the generation of the necessary input files and possible changes in symmetry, a script
initso lapw exists. For non-spinpolarized cases nothing particular must be taken into account
and SO can be easily applied by running run lapw -so. It will automatically use the complex
version of LAPW2.
However, for spin-polarized cases, the SO interaction may change (lower) the symmetry depend-
ing on how you choose the direction of magnetization and care must be taken to get a proper setup.
initso lapw together with symmetso generates the proper symmetry.
Just a few hints what can happen:
I Suppose you have a cubic system and put the magnetization along [001]. This will create
a tetragonal symmetry (and you can temporarely tell this to the initialization programs by
changing the respective lattice parameter c to a tetragonal system).
I If you put the magnetization along [111], this creates most likely a rhombohedral (or hexag-
onal) symmetry. (Try to visualize this for a fcc lattice, XCRYSDEN is very usefull for this
purpose).
4.5. FLOW OF PROGRAMS 47
I Symmetry operations can be classified into operations which invert the magnetization,others
which leave it unchanged and some which do some arbitrary rotation. The program
symmetso (part of initso lapw) sorts these operations in the proper way.
I If you don’t have inversion symmetry in the original structure, you must not “add inversion”
in KGEN.
The recommended way to include SO in the calculations is to run a regular scf calculation first,
save the results, initialize SO and run another scf cycle including SO:
I run[sp] lapw
I save lapw case nrel
I initso lapw
I run[sp] lapw -so
For spin-polarized systems you may want to add the “-dm” switch to calculate also the orbital
magnetic moment.
If you want to force a non-magnetic solution you can constrain the spin-polarization to zero using
runsp c lapw.
Without SO, [Link]/dn will be considered in LAPW1(c). With SO, it will be applied in
LAPWSO (and allows coupling of nondiagonal spin-terms).
I LDA-Hartree-Fock
Functional 5 in case.in0. mode = EECE and fraction = 1 in [Link].
LDA-HF [ρ] = E LDA [ρ] + E HF [Ψ LDA
Exc xc x corr ] − Exc [ρcorr ]
48 CHAPTER 4. FILES AND PROGRAM FLOW
I LDA-Fock-α
Functional 5 in case.in0. mode = HYBR and fraction = α in [Link].
LDA-Fock-α [ρ] = E LDA [ρ] + α E HF [Ψ LDA [ρ
Exc xc x corr ] − Ex corr ]
I PBE-Fock-α
Functional 13 in case.in0. mode = HYBR and fraction = α in [Link].
PBE-Fock-α [ρ] = E PBE [ρ] + α E HF [Ψ
Exc ] − E PBE [ρ ]
xc x corr x corr
I WC-Fock-α
Functional 11 in case.in0. mode = HYBR and fraction = α in [Link].
WC-Fock-α [ρ] = E WC [ρ] + α E HF [Ψ WC
Exc xc x corr ] − Ex [ρcorr ]
I TPSS-H-Fock-α
Functional 27 in case.in0. mode = HYBR and fraction = α in [Link].
TPSS-H-Fock-α [ρ] = E TPSS [ρ] + α E HF [Ψ TPSS [ρ
Exc xc x corr ] − Ex corr ]
In addition to the input files which are necessary for an usual LDA or GGA calculation, the
input file [Link] is necessary to start a calculation. You may copy a template from
SRC templates, but must modify it according to your needs. In particular you must select for
which atoms and which orbitals (usually d-Orbitals of late transition metal atoms or f-orbitals for
4f/5f atoms) you want to add such a potential and which type of functional you want to use.
A sample input for calculations with exact exchange is given below.
emin lower energy cutoff, to be selected so that the energy of correlated states
is larger than emin
natorb number of atoms for which the exact exchange is calculated
As with LDA+U , hybrid functionals can be used only for spin-polarized calculations (runsp lapw
with the switch -eece). runsp lapw will internally call runeece lapw, which will cre-
ate all necessary additional input files (it requires a case.in0 file including the op-
tional IFFT line as generated by init lapw): [Link] ([Link]), [Link],
case.in0eece, case.in2eece (case.in2ceece) and once this is done, calculates in a se-
ries of lapw2/lapwdm/lapw0/orb calculations the corresponding orbital dependend potentials.
The modified Becke-Johnson exchange potential + LDA-correlation (Tran and Blaha 2009) allows
the calculation of band gaps with an accuracy similar to very expensive GW calculations. It is a
local approximation to an atomic “exact-exchange”-potential and a screening term. This is just a
XC-potential, not a XC-energy functional, thus Exc is taken from LSDA and the forces cannot be
used with this option.
We recommend the following steps to perform such calculations:
I cp case.in0 case.in0 grr and change indxc in case.in0 grr to 50. This option will
calculate the average of ∇ρ/ρ over the unit cell. (The presence of case.in0 grr will
be detected during the scf-cycle and lapw0 will be called twice, first with the input file
case.in0 grr, then with case.in0.)
I edit [Link] and change to PRATT mixing.
I run another scf cycle.
As mentioned above, you should use PRATT mixing since in most cases MSEC1 mixing will lead
to convergence problems of the scf cycle (V also depends on the kinetic energy density and this is
not mixed in mixer). PRATT mixing can be slow, or can start to oszillate and [Link], first
using a smaller mixing factor (eg. 0.2), later increasing it to about 0.50 to make sure that you do not
stop at false convergence.
The TB-mBJ potential uses an average of ∇ρ/ρ over the unit cell. This does not make sense for
surfaces or molecules. In such cases, run a similar bulk structure first, then cp case [Link] to
[Link] and remove case.in0 grr. This runs mBJ with a fixed value of “c”.
5 Shell scripts for running programs
Contents
5.1 Job control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.2 Utility scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.3 Structure optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4 Phonon calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.5 Parallel Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.6 Getting on-line help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.7 Interface scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
In order to run WIEN2k several c-shell scripts are provided which link the individual programs to
specific tasks.
All available (user-callable) commands have the ending lapw so you can easily get a list of all
commands using
ls $WIENROOT/∗ lapw
in the directory of the WIEN2k executables. (Note: all of the more important commands have a link to a
short name omitting “ lapw”.) All these commands have at least one option, -h, which will print a
small help indicating purpose and usage of this command.
The main WIEN2kscript, x lapw or x, executes a single WIEN2kprogram. First it creates the corre-
sponding [Link]-file, where the connection between Fortran I/O-units and filenames are
defined. One can modify its functionality with several switches, modifying file definitions in case
of spin-polarized or complex calculations or tailoring special behaviour. All options are listed with
the help switch
x -h or x lapw -h
With some of the options the corresponding input files may be changed temporarely, but are set
back to the original state upon completion.
51
52 CHAPTER 5. SHELL SCRIPTS
FLAGS:
-f FILEHEAD -> FILEHEAD for path of struct & input-files
-t/-T -> suppress output of running time
-h/-H -> help
-d -> create only the def-file
-up -> runs up-spin
-dn -> runs dn-spin
-du -> runs up/dn-crossterm
-sc -> runs semicore calculation
-c -> complex calculation (no inversion symmetry present)
-p -> run lapw1/2/so in parallel (needs .machines file)
-orb -> runs lapw1 with LDA+U/OP or B-ext correction
-it -> runs lapw1 with iterative diagonalization
-noHinv -> runs lapw1 with iterative diag. without Hinv
-noHinv0 -> runs lapw1 with iterative diag. writing new Hinv
-nohns-> runs lapw1 without HNS
-nmat_only-> runs lapw1 and yields only the matrixsize
-in1orig -> runs lapw2 but does not modify case.in1
-emin X -> runs lapw2 with EMIN=X (in bin9_blaha.in2)
-all X Y -> runs lapw2 with ALL and E-window X-Y (in bin9_blaha.in2)
-qtl -> calculates QTL in lapw2
-alm -> calculates ALM,BLM in lapw2
-almd -> calculates ALM,BLM in lapw2 for DMFT (Aichhorn/Georges/Biermann)
-qdmft -> calculates charges including DMFT (Aichhorn/Georges/Biermann)
-help_files -> creates [Link] files in lapw2
-vresp-> creates [Link] (for TAU/meta-GGA) in lapw2
-eece -> for hybrid-functionals (lapw0,lapw2,mixer,orb,sumpara)
-grr -> lapw0 for grad rho/rho (using SRC.in0_grr with ixc=50)
-band -> for bandstructures: unit 4 to 5 (in1), sets QTL and ROOT (in2)
-fermi-> calculates Fermi energy and weights in lapw2
-efg -> calculates lapw2 with EFG switch
-so -> runs lapw2 with def-file for spin-orbit calculation
-fbz -> runs kgen and generates a full mesh in the BZ
-fft -> runs dstart only up to case.in0_std creation
-super-> runs dstart and creates new_super.clmsum (and not [Link])
-sel -> use reduced vector file in lapw7
-settol 0.000x -> run sgroup with different tolerance
-sigma-> run lstart with case.inst_sigma (autogenerated) for [Link].
-rxes-> run tetra using [Link] weight file for RXES-spectroscopy.
-rxesw E1 E2-> run tetra and create [Link] file for RXES for energies E1-E2
-delta-> run arrows program with difference between two structures
-lcore-> runs dstart with [Link] (produces [Link])
-copy -> runs pairhess and copies .minpair to .minrestart and .minhess
-telnes -> run qtl after generating [Link] based on [Link]
USE: x -h PROGRAMNAME for valid flags for a specific program
Note: To make use of a scratch file system, you may specify such a filesystem in the environment variable
SCRATCH (it may already have been set by your system administrator). However, you have to make sure
that there is enough disk-space in the SCRATCH directory to hold your [Link]* and [Link]*
files.
In order to start a new calculation, one should make a new subdirectory and run all calculations
from there. At the beginning one must provide at least one file (see 3), namely [Link]
(see 4.3) ([Link] can be created automatically on the “fly”, see 6.4.3), then one runs a series of
programs using init lapw. This script is described briefly in chapter 4.5) and in detail in “Getting
started” for the example TiC (see chapter 3). You can get help with switch -h. All actions of this
script are logged in short in :log and in detail in the file [Link], which also gives you a
“restart” option when problems occurred. In order to run init lapw starting from a specific point
on, specify -s PROGRAM.
5.1. JOB CONTROL 53
Ignoring ERRORS and in many cases also WARNINGS during the execution of this script, most
likely will lead to errors at a later stage. Neglecting warnings about core-leakage creates .lcore,
which directs the scf-cycle to peform a superposition of core densities.
init lapw supports switch -b, a “batch” mode (non-interactive) for trivial cases AND experi-
enced users. You can supply various options and specify spin-polarization, XC-potential, RKmax,
k-mesh or mixing. See init lapw -h for more details. Changes to [Link] by nn will be
accepted, but by sgroup will be neglected. Please check the terminal output for ERRORS and WARN-
INGS !!!
In order to perform a complete SCF calculation, several types of scripts are provided with the
distribution. For the specific flow of programs see chapter 4.5.
Cases with/without inversion symmetry and with/without semicore or core states are handled au-
tomatically by these scripts. All activities of these scripts are logged in short in :log (appended)
and in detail together with convergence information in [Link] (overwriting the old “day-
file“). You can always get help on its usage by invoking these scripts with the -h flag.
run lapw -h
PROGRAM: /zeus/lapw/WIEN2k/bin/run_lapw
OPTIONS:
-cc LIMIT -> charge convergence LIMIT (0.0001 e)
-ec LIMIT -> energy convergence LIMIT (0.0001 Ry)
-fc LIMIT -> force convergence LIMIT (1.0 mRy/a.u.)
default is -ec 0.0001; multiple convergence tests possible
-e PROGRAM -> exit after PROGRAM ()
-i NUMBER -> max. NUMBER (40) of iterations
-s PROGRAM -> start with PROGRAM ()
-r NUMBER -> restart after NUMBER (99) iterations (rm *.broyd*)
-nohns NUMBER ->do not use HNS for NUMBER iterations
-in1new N -> create "new" in1 file after N iter (write_in1 using scf2 info)
-ql LIMIT -> select LIMIT (0.05) as [Link] for E-L setting in new in1
-qdmft NP -> including DMFT from Aichhorn/Georges/Biermann running on NP proc
FLAGS:
-h/-H -> help
-I -> with initialization of in2-files to "TOT"
-NI -> does NOT remove [Link]* (default: rm *.broyd* after 60 sec)
-p -> run k-points in parallel (needs .machine file [speed:name])
-it -> use iterative diagonalizations
-it1 -> use iterative diag. with recreating H_inv (after basis change)
-it2 -> use iterative diag. with reinitialization (after basis change)
-noHinv -> use iterative diag. without H_inv
-vec2pratt -> use vec2pratt instead of vec2old for iterative diag.
-so -> run SCF including spin-orbit coupling
-renorm-> start with mixer and renormalize density
-in1orig-> if present, use case.in1_orig file; do not modify case.in1
54 CHAPTER 5. SHELL SCRIPTS
CONTROL FILES:
.lcore runs core density superposition producing [Link]
.stop stop after SCF cycle
.fulldiag force full diagonalization
.noHinv remove [Link] files
case.inm_vresp activates calculation of vresp files for meta-GGAs
case.in0_grr activates a second call of lapw0 (mBJ pot., or E_xc analysis)
ENVIRONMENT VARIBLES:
SCRATCH directory where vectors and help files should go
Additional flags valid only for magnetic cases (runsp lapw) include:
-dm -> calculate the density matrix (when -so is set, but -orb is not)
-eece -> use "exact exchange+hybrid" methods
-orb -> use LDA+U, OP or B-ext correction
-orbc -> use LDA+U correction, but with constant V-matrix
Calling run lapw (after init lapw) from the subdirectory case will perform up to 40 iterations
(or what you specified with switch -i) unless convergence has been reached earlier. You can choose
from three convergence criteria, -ec (the total energy convergence is the default and is set to 0.0001
Ry for at least 3 iterations), -fc (magnitude of force convergence for 3 iterations, ONLY if your
system has “free” structural parameters!) or -cc (charge convergence, just the last iteration), and
any combination can also be specified. Be careful with these criteria, different systems will require
quite different limits (e.g. fcc Li can be converged to µRy, a large unit cell with heavy magnetic
atoms only to 0.1 mRy). You can stop the scf iterations after the current cycle by generating an
empty file .stop (use eg. touch .stop in the respective case-directory).
The scf-cycle creates [Link]* files which contain the ”charge-history”. Once run lapw has
finished, you should usually ”save lapw” (see below) the results. When you continue with an-
other run lapw without ”save lapw” (because the previous run did not fulfill the convergence
criteria or you want to specify a more strict criterium) the ”broyden-files” will be deleted unless
you specify -NI.
With -e PROGRAM you can run only part of one scf cycle (e.g. run lapw0, lapw1 and lapw2),
with -s PROGRAM you can start at an arbitrary point in the scf cycle (e.g. after a previous cycle
has crashed and you want to continue after fixing the problem) and continue to self-consistency.
Before mixer is invoked, [Link] is copied to [Link] old, and the final “important“
files of the scf calculation are [Link] and [Link].
Invoking
will first set in case.in2 the TOT-switch (if FOR was set) to save cpu time, then run up to 30 scf cycles
till the force criterion of 0.5 mRy/a.u. is met (for 3 consecutive iterations). Then the calculation of
all terms of the forces is activated (setting FOR in case.in2) for a final iteration.
By default the file case.in1 is updated after lapw2 and the current Fermi-energy is inserted.
This will force lapw1 to use instead of the default energy parameters (0.30) an energy “EF − 0.2”.
The switch -in1orig can be used to keep the present case.in1 file unmodified (or to copy
case.in1 orig back after -in1new).
The switch -in1new N preserves for N iteration the current case.in1 file. After the first N
iterations write in1 lapw is called and a new case.in1 file is generated, where the energy pa-
rameters are set according to the :EPLxx and :EPHxx values of the last scf iteration and the -ql
value (see sections 4.4 and 7.3). In this way you may select in some cases better energy-parameters
and also additional LOs to improve the linearization may be generated automatically. Note, how-
ever, that this option is potentially unsave and dangerous, since it may set energy-parameters of
LOs and APW+lo too close (leading to ghostbands) or in cases where you have a “bad” last iter-
ation (or large changes from one scf iteration to the next. The original case.in1 file is saved in
case.in1 orig and is used as template for all further scf-cycles.
5.2. UTILITY SCRIPTS 55
saves [Link], [Link], [Link], [Link] and [Link] under the new
name and removes the [Link]* files. Now you are ready to modify structural parameters
or input switches and rerun run lapw, or calculate properties like charge densities (lapw5), total
and partial DOS (tetra) or energy bandstructures (spaghetti).
For more complicated situations, where many parameters will be changed, we have extended
save lapw so that calculations can not only be saved under the head of save filename but
also a directory can be specified. If you use any of the possible switches (-a, -f, -d, -s) all input files
will be saved as well (and can be restored using restore lapw).
Options to save lapw can be seen with
save lapw -h
To restore a calculation the script restore lapw can be used. This script restores the struct,
clmsum, vorb and dmat files as well as all input files. Note: The input files will only be restored
when save lapw -d was used, i.e. when you have saved a calculation in an extra directory.
After restore lapw you can continue and either run an scf cycle (run lapw) or recreate the
scf-potential (x lapw0) and the corresponding eigenvectors (x lapw1) for further tasks (DOS,
electron density,...).
Options to restore lapw are:
-h help
-f force restore lapw to overwrite previous files
-d directory restore calculation from directory specified
-s silent operation (no output)
-t only test which files would be restored
Once a case has been completed you can clean up the directory with this command. Only the most
important files (scf, clmsum, struct, input and some output files) are kept. It is very important to
use this command when you have finished a case, since otherwise the large vector and helpXX files
will quickly fill up all your disk space.
This script migrates a case to a remote computer (to be called within the case-dir). Needs working
ssh/scp without password; local and remote case-dir must have the same name.
Call it within the desired case-dir as:
migrate lapw [FLAGS OPTIONS] [user@]host:path/case-dir
with the following options:
FLAGS:
-h -> help
-clean -> a clean_lapw is issued before copying
-r -> files in source directory are removed after copying
-R -> source directory (and all files) are removed after copying
-s -> do it silent (in batch mode)
-z -> gzip files before scp (slow network)
5.2. UTILITY SCRIPTS 57
which will prepare [Link] (eg. for the TiC example) with:
tic #Title
-1.000 0.00250 1.200 0.003 #Emin, DE, Emax, Gauss-Broad
8 #Number of DOS
0 1 total-DOS
1 1 tot-Ti
1 4 d-Ti
1 5 d-eg-Ti
1 6 d-t2g-Ti
2 1 tot-C
2 2 s-C
2 3 p-C
58 CHAPTER 5. SHELL SCRIPTS
This script searches for .running.* files within the current directory (or the directory specified
with “-d full path directory”) and then performs a ps command for these processes. If the specified
process has not been found, it removes the corresponding .running.* file after confirmation
(default) or immediately (when “-f” has been specified).
This script searches for .running.* files within the current directory (or the directory specified
with “-d full path directory”) and then kills the corresponding process after confirmation (default)
or immediately (when “-f” has been specified). It is particular usefull for killing “k-point parallel”
jobs.
This script extracts the critical points (CP) after a Bader analysis (x aim (-c)) from
[Link]. It sorts them (according to the density), removes duplicate CPs, converts units
into Å, e/Å3 , ... and produces critical points ang.
It is used with: extractaim lapw [Link]
Hartmut Enkisch
Institute of Physics E1b
University of Dortmund
Dortmund, Germany
enkisch@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
-h help switch
-i n show only the last n iterations
-f scf-file use "scf-file" instead of the default "[Link]"
-p produces file "[Link]" instead of X-window plot
arg1,... arguments to monitor (like ":ENE" or ":DIS" , see analyse_lapw )
5.2. UTILITY SCRIPTS 59
The scfmonitor can also be called directly from w2web using the ”Analyse” tool.
In order to have a reasonable behavior of scfmonitor the GNUPLOT window should stay in back-
ground. This can be achieved by putting a line into your .Xdefaults file like:
gnuplot*raise: off
Note: It does not make sense to start scfmonitor before the first cycle has finished because no [Link]
exists at this point.
The script analyse lapw is usually called from scfmonitor lapw. It ”greps” from an scf-file
the specified arguments and produces [Link].
analyse lapw is called using:
analyse lapw [-h] scf-file arg1 [arg2 arg3 arg4 arg5 arg6]
-h help switch
scf-file "scf-file" to analyse (there’s no default "[Link]" !)
arg1,... arguments to analyse:
atom independend: :ENE :DIS :FER :MMT
atom iii dependend: :CTOiii :CUPiii :CDNiii :NTOiii :NUPiii :NDNiii
:DTOiii :DUPiii :DDNiii :RTOiii :EFGiii :HFFiii
:MMIiii
vector quantities: :FORiii[x/y/z] :POSiii[x/y/z] :FGLiii[x/y/z]
where magnitude z z is the defaul
For vector quantities like :FGLiii or :POSiii (usefull with [Link] mini) one can specify the
respective coordinate by adding x/y/z to the corresponding labels.
testpara lapw is a small script which helps you to determine an optimal selection for the file
.machines for parallel calculations (see sec. 5.5).
testpara1 lapw is a small script which determines how far the execution of lapw1para has
proceeded.
testpara2 lapw is a small script which determines how far the execution of lapw2para has
proceeded.
60 CHAPTER 5. SHELL SCRIPTS
Using
you can get a list of a quantity “:label” (e.g. :ENE for the total energy) from several scf files at
once.
initso lapw helps you to initialize the calculations for spin-orbit coupling. It helps together
with make inso lapw (contributed by Morteza Jamal, m jamal57@[Link]) to create/modify
all required input files ([Link], case.in1, case.in2c). In a spinpolarized case SO may
reduce symmetry or equivalent atoms may become non-equivalent, and the script calls symmetso
and will help you to find proper symmetries and setup the respective input files. It is called using
initso lapw or
initso
and you should carefully follow the instructions and explanations of the script and the explanations
for [Link] given in section 7.4.
clmextrapol lapw extrapolates the charge density ([Link]/up/dn) from old to new po-
sitions (or from old to new lattice parameters). It takes the density from the old positions (copied
into [Link]) and subtracts an atomic superposition density (new [Link]) fom the
old positions and adds an atomic superposition density fom the new ones (generated by dstart).
5.3. STRUCTURE OPTIMIZATION 61
If new [Link] (generated automatically by init lapw) is not present, it will be generated
and for the next geometry step an extrapolation will take place.
It is usually called from “min lapw ” after a geometry step has finished and a new struct file has
been generated.
It can significantly reduce the number of scf-cycles for the new geometry step.
The auxilliary program optimize (x optimize) generates from an existing [Link] (or
case [Link], which is generated at the first call of optimize) a series of struct files
with various volumes (or c/a ratios, or other modified parameters) (depending on your input):
#!/bin/csh -f
foreach i ( \
tic_vol_-10.0 \
tic_vol__-5.0 \
tic_vol___0.0 \
tic_vol___5.0 \
tic_vol__10.0 \
)
cp $[Link] [Link]
# cp $[Link] [Link]
# x dstart
# run_lapw -ec 0.0001 -in1new 3 -renorm
run_lapw -ec 0.0001
set stat = $status
if ($stat) then
echo "ERROR status in" $i
exit 1
endif
save_lapw ${i}
# save_lapw -f -d XXX $i
end
You may modify this script according to your needs: use runsp lapw or even min lapw, or spec-
ify different convergence parameters; modify the save lapw command and change the save-name
or save into a directory to separate e.g. “gga” and “lda” results. Eventually you may activate the
line “ cp $[Link] [Link]” to use a previously saved clmsum file, e.g. from a calcula-
tion with smaller RKmax, ... and deactivate the ”clmextrapol lapw” lines, but usually the latter is
so efficient that this is no longer recommended.
Note: You must have a [Link] file (either from init lapw or from a previous scf calculation) in
order to run [Link].
After execution of this script you should have a series of scf-files with energies corresponding to the
modified parameters, which should allow you to find the corresponding equillibrium parameters.
For the volume optimization an analysis tool is available, other tools are under development).
62 CHAPTER 5. SHELL SCRIPTS
Using the script grepline (or the “Analysis o Analyze multiple SCF-files” menu of w2web) you
get a summary of the total energy vs. volume (c/a). The file [Link] can be used in
eplot lapw to find the minimum total energy and the equilibrium volume (c/a). Supported
equation of states include the EOS2, Murnaghan and Birch-Murnaghan EOS.
Using such strategies also higher-dimensional optimizations (e.g. c/a ratio and volume) are possi-
ble in combination with the -d option of save lapw.
For optimization of more degrees of freedom (2-4 lattice parameters), you can use the correspond-
ing option and for analysis of the data the script parabolfit lapw together with the program
eosfit6. It performs a non-linear least squares fit, using a parabolic fit-function in your variables
and get an analytic description of your energy surface. Please note, this is only a harmonic fit (no
odd or higher terms) and the description may not be very good if your parameter range is large
and/or the function is quite anharmonic, or you suffer from numerical noise.
For the determination of elastic constants see the description of ELAST in sec 8.14.
Package 2Doptimize
Morteza Jamal
email: m jamal57@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
This package performs a convenient 2D structure optimization (Volume and c/a, i.e. tetragonal
or hexagonal spacegroups). After initialization of a case, one generates a set of structures and a
job-file [Link] using the command
set2D lapw
This calls setup2D and you have to specify the “order of the fit” (stored in .fordfitcoa) and
the changes in volume and c/a. The resulting [Link] script should be adapted (eg. use
min lapw instead of run lapw; insert switches,...) and executed. Finally
ana2D lapw
can be executed and will analyze the results. It uses a set of [Link]* files (produced by
[Link] and stored also in subdirectory Vconst) and the numbvcoa file. You can mod-
ify the order of the fits in .fordfitcoa and rerun ana2D lapw to check the sensitivity of the
results with the fitting. Note: Fits of high order (and few “data points”) may lead to artificial
results due to unphysical oszillations of the fit.
You can see results for
- energy vs. c/a for each volume,
- energy vs. volume (with optimized c/a) and
- c/a vs. volume.
5.3. STRUCTURE OPTIMIZATION 63
Optionally you can specify more cases by rerunning set2D lapw. Specify also your ‘‘old’’
volume and c/a points again (or leave them out on purpose in case they were very bad
(eg. very far from the minimum). The old results will then be taken automatically into account
without recalculation (unless you modify [Link], see the comments at the top of this
file). Thus a “good” strategy is to use only 3x3 points (order of fit = 3) and in a second step you
add points where they are needed.
When you want to rerun such an optimization with different parameters (RKmax, k-mesh, XC-
potentials) modify the top of [Link] and set answscf=no and a new savename (eg.
” pbe rk8 1000k”).
Most of the more complicated structures have free internal structural parameters, which can either
be taken from experiment or optimized using the calculated forces on the nuclei.
Starting with WIEN2k 11.1 there are two possibilities to determine the equilibrium position of all
individual atoms automatically (obeying the symmetry constraints of a certain space group). One
can use either
I the shell script min lapw, together with the program mini, which will run a scf-cycle, update
the positions using the calculated forces and restarts a new scf cycle. This continues until
forces drop below a certain value;
I or use the normal scf-scripts run lapw with a special input in [Link] such that the charge
density and the positions are simultaneously optimized during the scf-cycle.
At present we recommend the second (new) option, although there are cases where this scheme
can be slower or may even fail to converge.
A typical sequence of commands for an optimization of the internal positions would look like:
I or edit [Link] and put MSR1a (or MSEC1a) as “mixing method”. Then continue with
run lapw -fc 0.5 -ec 0.0001 -cc 0.001 [-it]. It will run x pairhess (unless
[Link] is already present) and then run (several hundreds) scf-cycles, simultaneously
updating positions and charge densities. Once the forces seem to be smaller than the limit
defined in [Link] it will switch to “mixing method” MSR1 and finalize the scf-cycle with
fixed positions. Because of this, the final forces may not be as small as desired and eventually
you have to restart this step using MSR1a again.
When using the second method we recommend you read carefully $WIEN-
ROOT/SRC mixer/Mixer README [Link]. Overall the method is very good for semiconductors
(or well behaved metals), and allows “tricks” like small k-mesh or small RKMax at the beginning
of the minimization and using higher accuracy only towards the end.
The following text refers (mainly) to the first method using min lapw:
When [Link] is not present, an scf-cycle will be performed first, otherwise the corresponding
forces are extracted into [Link] and the program mini generates a new [Link] with
modified atomic positions. The previous step is saved under case 1/2/3.... Then a new scf-
cycle is executed and this loop continues until convergence (default: forces below 2mRy/bohr) is
reached.
The last iteration of each geometry step is appended to [Link] mini, so that this file contains
the complete history of the minimization and can be used to monitor the progress (grep :ENE *mini;
or :FORxxx ...).
By default (unless switch -noex is specified), min will call the script clmextrapol lapw after the
first geometry step and try to extrapolate the charge density to the new positions. This procedure
usually significantly reduces the number of scf-cycles and is thus highly recommended.
mini requires an input file [Link] (see Sec. 8.15) which is created automatically and MUST
NOT be changed while min lapw is running (except the force tolerance, which terminates the
optimization).
We recommend the PORT minimization method, a reverse-communication trust-region Quasi-
Newton method from the Port library, which seems to be stable, efficient and does not depend too
much on the users input (DELTAs, see below with NEWT). The PORT option also uses/produces a
file .min hess, which contains the (approximate) Hessian matrix (lower-triangle Cholesky factor)
If you restart a minimization with different k-points, RMT, RKmax, ... or do a similar calcula-
tion (eg. for a different volume, ...) it will be copied to .minrestart (unless -nohess is speci-
fied), so that you start with a reasonable approximation for the Hessian. The program pairhess,
which calculates the first Hessian, also prints out the average Hessian eigenvalue for the symmet-
ric, symmetry-preserving modes in mRyd/au2 as well as the minimum and maximum, and also
the vibration frequencies. A list of these is given at the end of [Link]. Note that these
are not all possible modes, but only the symmetry preserving ones. Therefore if you have prior
information about the vibrations of the system you can adjust the rescaling term so the average
vibration frequency is about correct. (see the description of pairhess in 9.2). (In addition there is a
program eigenhess, which will analyze the Hessian after the minimization has been completed.
It also prints vibrational frequencies and may give you hints about dynamical instability of your
system. Some more description is given in $WIENROOT/SRC pairhess/README and at the top
of the output file [Link].
When using PORT you may also want to check its progress using
where :LABEL is :ENE (should decrease), :GRAD (should also go down, but could sometimes also
go up for some time as long as the energy still decreases), :MIN (provides a condensed summary
of the progress), :WARN may indicate a problem), :DD (provides information about the step sizes
and mode used). Some general explanations are:
5.3. STRUCTURE OPTIMIZATION 65
1) The algorithm takes steps along what it considers are good directions (using some internal logic),
provided that these steps are smaller than what is called the trust-region radius. After a good step
(e.g. large energy decrease) it expands the trust-region; after a bad one it reduces it. Sometimes it
will try too large a step then have to reduce it, so the energy does not always go down. You can see
this by using ”:DD” and “:MIN” .
2) A grep on :MIN gives a condensed progress output, in which the most significant terms are
E (energy in some rescaled units), RELDF (last energy reduction), PRELDF (what the algorithm
predicted for the step), RELDX (RMS change in positions in Angstroms) and NPRELDF (predicted
change in next cycle). Near the solution RELDF and RELDX should both become small. However,
sometimes you can have soft modes in your structure in which case RELDX will take a long time
before it becomes small.
3) A warning that the step was reduced due to overlapping spheres if it happens only once (or
twice) is not important; the algorithm tested too large a step. However, if it occurs many times it
may indicate that the RMT’s are too big.
4) A warning ”CURVATURE CONDITION FAILED” indicates that you are still some distance from
the minimum, and the Hessian is changing a lot. If you see many of these, it may be that the forces
and energy are not consistent.
Sometimes PORT gets ”stuck” (often because of inconsistencies of energy and forces due to in-
sufficient scf convergence or a very non-harmonic potential energy surface). A good alternative is
NEW1, which is a ”sophisticated” steepest-descent method with optimized step size. It can be very
efficient in certain cases, but can also be rather slow when the potential energy surface is rather flat
in one, but steep in another direction (eg. a weakly bound molecule on a surface, but constraining
the sensitive parameters, like the bond distance of the molecule, may help).
Another alternative is NEWT, where one must set proper ”DELTAs” and a ”FRICTION” for each
atom. Unfortunately, these DELTAs determine crucially how the minimization performs. Too small
values lead to many (unnecessary) ”geometry steps”, while too large DELTAs can even lead to
divergence (and finally to a crash). Thus you MUST control how the minimization performs. We
recommend the following sequence after 2-3 geometry steps:
Good, since the force (only a force along z is present here) is decreasing reasonably fast towards
zero. You must check this for every atom in your structure.
When you detect oszillations or too small changes of the forces during geometry optimization, you
will have to decrease/increase the DELTAs in [Link] and rm [Link]. (NOTE: You must
not continue with modified DELTAs but keeping [Link].) Alternatively, stop the minimiza-
tion (touch .minstop and wait until the last step has finished), change [Link] and restart.
You can get help on its usage with:
PROGRAM: min
OPTIONS:
-j JOB -> job-file JOB (default: run_lapw -I -fc 1. -i 40 )
-noex -> does not extrapolate the density for next geometry step
-p -> adds -p (parallel) switch to run_lapw
-it -> adds -it (iterative diag.) switch to run_lapw
-it1 -> adds -it1 ([Link]. with recreating H_inv) switch to $job
-it2 -> adds -it2 ([Link]. with reinitialization) switch to $job
-noHinv -> adds -it -noHinv ([Link]. without H_inv) switch to $job
-sp -> uses runsp_lapw instead of run_lapw
-nohess -> removes .minrestart (initial Hessian) from previous minimization
-m -> extract force-input and execute mini (without JOB) and exit
-mo -> like -m but without copying of case.tmpM1 to [Link]
-h/-H -> help
-NI -> without initialization of minimization (eg. continue after a crash)
-i NUMBER -> max. NUMBER (50) of structure changes
-s NUMBER -> save_lapw after NUMBER of structure changes
CONTROL FILES:
.minstop stop after next structure change
For instance for a spin-polarized case, which converges more difficultly, you would use:
init phonon lapw uses case.d45 from PHONON and creates subdirectories case XX and
case [Link] files for all required displacements. It allows you to define globally RMT values
for the different atoms and
- initializes every case individually (batch option of init lapw is now supported) or
- initializes every second case (useful for pos. and neg. displacements, which have the same sym-
metry and thus only one initialization is necessary), or
- initializes only the first case and copies the files from the first case to all others. This is most
convenient in low symmetry cases with P1 symmetry for all cases and thus just one init lapw needs
to be executed (while for higher symmetry a separate initialization is required (but computational
effort is reduced).
Please use mainly nn to reduce equivalent atoms. sgroup might change the unitcell and than the
collection of forces into the original supercell is not possible (or quite difficult).
A script run phonon has been created. Modify it according to your needs (parallelization,....) and
run all cases to selfconsistency.
5.5. PARALLEL EXECUTION 67
Note that good force convergence is essential (at least 0.1 mRy/bohr) and if your structure has
free parameters, either very good equillibrium positions must have been found before, or even
better, use both, positive and negative displacements to average out any resulting error from non-
equillibrium positions.
analyse phonon lapw uses the resulting scf files and generates the “Hellmann-Feynman”-file
required by PHONON. When you have positive and negative displacements an automatic averag-
ing will be performed. The resulting [Link] and [Link] filse should be transfered back to
MS-Windows and imported into PHONON.
Parts of the code are executed in parallel, namely LAPW1, LAPWSO, LAPW2, LAPWDM, and
OPTIC. These are the numerically intensive parts of most calculations.
Parallelization is achieved on the k-point level by distributing subsets of the k-mesh to different
processors and subsequent summation of the results. The implemented strategy can be used both
on a multiprocessor architecture and on a heterogeneous (even multiplatform) network.
To make use of the k-point parallelization, make sure that your system meets the following require-
ments:
NFS: All files for the calculation must be accessible under the same name and path. Therefore you
should set up your NFS mounts in such a way, that on all machines the path names are the
same.
68 CHAPTER 5. SHELL SCRIPTS
Remote login: rlogin or ssh to all machines must be possible without specifying a password.
Therefore you must either edit your .rhosts file to include all machines you intend to use
(not necessary for a shared memory machine), or correctly specify public/private keys for
ssh. This can be done by running “ssh-keygen -t rsa” and copying the id [Link]
key into /̃.ssh/authorized keys at the remote sites.
The command for launching a remote shell is platform dependent, and usually can be ’ssh’,
’rsh’ or ’remsh’. It should be specified during installation when siteconfig lapw is exe-
cuted (see chapter 11).
Fine grained MPI parallel versions are available for the programs lapw0, lapw1, and lapw2.
This parallelization method is based on parallelization libraries, including MPI, ScaLapack, PBlas
and FFTW 2.1.5 (lapw0). The required libraries are not included with WIEN2k. On parallel comput-
ers, however, they are usually installed. Otherwise, free versions of these libraries are available1 .
The parallelization affects the naming scheme of the executable programs: the fine grained par-
allel versions of lapw0/1/2 are called lapw0 mpi, lapw1[c] mpi, and lapw2[c] mpi. These
programs are executed by calls to the local execution environments, as in the sequential case, by
the scripts x, lapw0para, lapw1para, and lapw2para. On most computers this is done by
calling mpirun and should also be configured using siteconfig lapw.
To start the calculation in parallel, a switch must be set and an input file has to be prepared by the
user.
I The switch -p switches on the parallelization in the scripts x and run lapw.
I In addition to this switch the file .machines has to be present in the current working direc-
tory. In this file the machine names on which the parallel processes should be launched, and
their respective relative speeds must be specified.
If the .machines file does not exist, or if the -p switch is omitted, the serial versions of the pro-
grams are executed.
Generation of all necessary files, starting of the processes and summation of the results is done
by the appropriate scripts lapw1para, lapwsopara,lapwdmpara and lapw2para (when using
-p), and parallel programs lapw0 mpi, lapw1 mpi, and lapw2 mpi (when using fine grained
parallelization has been selected in the .machines file).
The following .machines file describes a simple example. We assume to have 5 computers, (al-
pha, ... epsilon), where epsilon has 4, and delta and gamma 2 cpus. In addition, gamma, delta and
epsilon are 3 times faster than alpha and beta.:
3:gamma:2 delta
3:delta:1 epsilon:4
residue:delta:2
lapw0:gamma:2 delta:2 epsilon:4
To each set of processors, defined by a single line in this file, a certain number of k-points is as-
signed, which are computed in parallel. In each line the weight (relative speed) and computers are
specified in the following form:
where weight is an integer (e.g. a three times more powerful machine should have a three times
higher weight). The name of the computer is machine name[1/2/...], and the number of pro-
cessors to be used on these computers are number[1/2/...]. If there is only one processor on a
given computer, the :1 may be omitted. Empty lines are skipped, comment lines start with #.
Assuming there are 8 k-points to be distributed in the above example, they are distributed as fol-
lows. The computers alpha and beta get 1 each. Two processors of computer gamma and one pro-
cessor of computer delta cooperate in a fine grained parallelization on the solution of 3 k-points,
and one processor of computer delta plus four processors of computer epsilon cooperate on
the solution of 3 k-points. If there were additional k-points, they would be calculated by the first
processor (or set of processors) becoming available. With higher numbers of k-points, this method
ensures dynamic load balancing. If a processor is busy doing other (e. g., interactive) work, the
overall calculation will not stall, but most of its work will be done by other processors (or sets of
processors using MPI). This is, however, not an implementation for fail safety: if a process does
not terminate (e. g., due to shutdown of a computer) the calculation will never terminate. It is up
to the user to handle with such hardware failures by modifying the .machines file and restarting
the calculation at the appropriate point.
During the run of lapw1para the file .processes is generated. This file is used by lapw2para
to determine which [Link]* to read.
By default lapw1para will generate approximately 3 vector-files per processor, if enough k-points
are available for distribution. The factor 3 is called “granularity” and allows for some load bal-
ancing in heterogeneous environments. If you can be sure that load balancing is not an issue (eg.
because you use a queuing-system and can be sure that you will get 100% of the cpus for your jobs)
it is recommended to set
granularity:1
residue:machine name:number
extrafine:1
70 CHAPTER 5. SHELL SCRIPTS
defines the computers used for running lapw0 mpi. In this example the 6 processors of the com-
puters gamma, delta, and epsilon run lapw0 mpi in parallel.
If fine grained parallelization is used, each set of processors defined in the .machines file is con-
verted to a single file .machine[1/2/...], which is used in a call to mpirun (or another parallel
execution environment).
When using a queuing system (like PBS, LoadLeveler or SUN-Gridengine) one can only request
the NUMBER of processors, but does not know on which nodes the job will run. Thus a “static”
.machines file is not possible. On can write a simple shell script, which will generate this file on
the fly once the job has been started and the nodes are assigned to this job. Examples can be found
at our web-site “[Link] users/faq”.
In the setup of the k-point parallel version of LAPW1 the list of k-points in [Link] (Note, that
the k-list from case.in1 cannot be used for parallel calculations) is split into subsets according to the
weights specified in the .machines file:
$ %
weighti ∗ klist
newweighti = P
granularity ∗ j weightj
lapw0 mpi is parallelized over the number of atoms. This method leads to good scalability as long
as there are more atoms than processors. For very many processors, however, the speedup is
limited, which is usually not at all critical, since the overall computing time of lapw0 mpi is
quite small.
lapw1 mpi uses a two-dimensional processor setup to distribute the Hamilton and overlap matri-
ces. For higher numbers of processors two-dimensional communication patterns are clearly
preferable to one-dimensional communication patterns.
Let us assume, for example, 64 processors. In a given processing step, one of these processors
has to communicate with the other 63 processors if a one-dimensional setup was chosen. In
the case of a two-dimensional processor setup it is usually sufficient to communicate with
the processors of the same processor row (7) or the same processor column (7), i. e. with 14
processors. p
In general
the processor array
P × Q is chosen as follows: P = number of processors ,
number of processors
Q= P . Because of SCALAPACK, often P × P arrays (i.e. 4, 9, 16,...
processors) give best performance and of course it is not recommended to use eg. 17 proces-
sors.
lapw2 mpi is parallelized in two main parts: (i) The density inside the spheres is parallelized over
atoms, and (ii) the fast Fourier transforms are done in parallel.
In addition the density calculation for each atom can be further parallelized by distributing
the eigenvector on a certain subset of processors (ususlly 2-4). This is not so efficient, but
most usefull if the memory requirement is too big otherwise. You set it in .machines using
lapw2 vector split:2
If more than one k-point is distributed at once to lapw1 mpi or lapw2 mpi, these will be treated
consecutively.
Depending on the parallel computer system and the problem size, speedups will vary to some
extend. Matrix setup in lapw1 should scale nearly perfect, while diagonalization (using SCALA-
PACK) will not. Usually, “iterative” scales better than “full” diagonalization and is preferred for
large scale computations. Scalability over atoms will be very good if processor and atom numbers
are compatible. Running the fine grained parallelization over a 100 Mbit/s Ethernet network is not
recommended, even for large problem sizes.
– Update information: When you think the program has an error, please check wether
newer versions are available, which might have fixed the problem you encounter.
– A mailing list:
Please check the ”digest”! In many cases your questions may have been answered be-
fore.
Locate your problem: If a calculation crashes, please locate the problem. Check
the content of files like [Link], *.error, [Link], [Link],
[Link] where X specifies the program which crashed.
Posting questions: Please provide enough information so that somebody can help you.
A question like: “My calculation crashed. Please help me!” will most likely not be
answered.
The script eplot lapw plots total energy vs. volume or total energy vs. c/a-ratio using the file
[Link]. The latter should have been created with grepline (using :VOL and :ENE
labels) or the “Analysis o Analyze multiple SCF-files” menu of w2web and the file names must be
generated (or compatible) with “[Link]”.
For a description of how to use the script for batch like execution call the script using
eplot lapw -h
The script parabolfit lapw is an interface for a harmonic fitting of E vs. 2-4-dim lattice param-
eters by a non-linear least squares fit (eosfit6) using PORT routines. Once you have several scf
calculations at different lattice parameters (usually generated with [Link]) it generates
the required [Link] and [Link] from your scf files. Using
you can optionally specify the dimensionality of the fit or the specific scf-filenames.
The script dosplot lapw plots total or partial Density of States depending on the input used
by [Link] and the interactive input. A more advanced plotting interface is provided by
dosplot2 lapw, see below.
For a description of how to use the script for batch like execution call the script using
dosplot lapw -h
74 CHAPTER 5. SHELL SCRIPTS
The script dosplot2 lapw plots total or partial Density of States depending on the input used by
[Link] and the interactive input. It can plot up to 4 DOS curves into one plot, and simultane-
ously plot spin-up/dn DOS.
It was provided by Morteza Jamal (m jamal57@[Link]).
For a description of how to use the script for batch like execution call the script using
dosplot2 lapw -h
The script Curve lapw plots x,y data from a file specified interactively. It asks for additional
interactive input. It can plot up to 4 curves into one plot and is a simple gnuplot interface.
It was provided by Morteza Jamal (m jamal57@[Link]).
specplot lapw provides an interface for plotting X-ray spectra from the output of the xspec or
txspec program.
For a description of how to use the script for batch like execution call the script using
specplot lapw -h
The script rhoplot lapw produces a surface plot of the electron density from the file [Link]
created by lapw5.
Note: To use this script you must have installed the C-program reformat supplied in SRC reformat.
The script opticplot lapw produces XY plots from the output files of the optics package us-
ing the [Link], [Link], [Link], [Link] or [Link]. For a
description of how to use the script for batch like execution call the script using
opticplot lapw -h
The script addjoint-updn lapw adds the files [Link] and [Link] together
and produces [Link]. It uses internally the program add columns. It should be called for
spin-polarized optics calculations after x joint -up and x joint -dn, because the Kramers-
Kronig transformation to the real part of the dielectric function (1 ) is not a simple additive quan-
tity concerning the spin (see Ambrosch-Draxl 06). The KK transformation should then be done
non-spinpolarized (x kram) resulting in files: [Link], [Link], [Link]
or [Link].
6 Programs for the initialization
Contents
6.1 NN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.2 SGROUP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3 SYMMETRY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4 LSTART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.5 KGEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.6 DSTART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
In sections (6.1-6.6) we describe the initial utility programs. These programs are used to set up a
calculation.
75
76 CHAPTER 6. INITIALIZATION
6.1.1 Execution
nn [Link] or x nn
6.2 SGROUP
6.2.1 Execution
6.3 SYMMETRY
This program uses information from [Link] (lattice type, atomic positions). If NSYM was
set to zero it generates the space group symmetry operations and writes them to [Link] st
to complete this file. Otherwise (NSYM > 0) it compares the generated symmetry operations
with the already present ones. If they disagree a warning is given in the output. In addition
the point group of each atomic site is determined and the respective symmetry operations and
LM values of the lattice harmonics representation are printed. The latter information is written
into case.in2 sy, while the local rotation matrix, the positive or negative IATNR values and the
proper ISPLIT parameter are written to [Link] st. (See appendix A and Sec. 4.3).
6.4. LSTART 77
6.3.1 Execution
NSTOP= .....
in [Link], this means, that a proper solution for at least one orbital could not be ob-
tained. In such a case the input must be changed and one should provide different occupation
numbers for these states (e.g. Cu can not be started with 3d10 4s1 , but it works with 3d9 4s2 ).
The program produces “WARNINGS” if R0 is too big or core-density leaks out of RMT.
6.4.1 Execution
The files [Link](up|dn) are generated and contain the atomic (spin) densities, which will be
used by DSTART later on.
Using -sigma generates [Link] sigma with modified input to generate [Link] used
for difference densities (see below).
The following parameters are defined in file [Link] (static and not allocatable arrays):
NPT total number of radial mesh points, must be gt.(NRAD+NPT00), where NRAD is
the number of mesh-points up to RMT specfied in [Link].
NPT00 max. number of radial mesh points beyond RMT
RMAX0 max. distance of radial mesh
78 CHAPTER 6. INITIALIZATION
6.4.3 Input
When running lstart you will first be asked interactively to specify an XC-potential switch. Cur-
rently 5 (LSDA, Perdew and Wang 92) as well as 11, 13 and 19 (three GGAs, Wu,Cohen 06; the
standard “PBE” Perdew et al. 96, as well as “PBEsol”, Perdew et al. 08; respectively) are officially
supported, 13 is the “standard PBE-GGA”.
In addition the program asks for an energy cut-off, separating core from valence states. Usually
-6.0 Ry is a good choice, but you should check for each atom how much core charge leaks out of the
sphere (WARNINGS in [Link]). If this is the case one should lower this energy cut-off
and thus include these low lying states into the valence region. Alternatively you can also select
a “charge localization” criterium (usually between 0.97 and 0.9999). This allows a more localized
state (like a 4f of 5d elements) to be core, while a more delocalized state at lower energy (like the
5p states of 5d elements) to be semi-core.
The rest of the input is described in the sample input below.
Note: Only the data at the beginning of the line are read whereas the comment describes the respective
orbitals. This file can be generated automatically in w2web during “Initialize calc. or using “Sin-
glePrograms o instgen lapw” or with the script instgen lapw. To edit this file by hand choose
“View/Edit o Input Files” and choose [Link].
line 1: format(a4,a6)
title, keyword
title
keyword The keyword Watson enables a stabilization of negative ions using a
“Watson”-sphere of radius R-wat with charge Q-wat, which must be
given in the next line when this keyword is specified.
The keyword PRATT enables a scf mixing using standard PRATT
scheme. It might be usefull if a certain atomic configuration does not
converge with the standard mixing scheme and requires a (usually
quite small) mixing factor, which must be given in the next line when
this keyword is specified.
config specifies the core state configuration by an inert gas (He, Ne, Ar, Kr,
Xe, Rn) and the number of (valence) orbitals (without spin). (In the
example given above one could also use Ar 3 and omit the 3s and 3p
states.) The atomic configurations are listed in the appendix and can
also be found online using periodic table, a shell script which dis-
plays SRC/[Link] with ghostview)
line 3: format(i1,1x,i2,1x,f5.3,a1)
n, kappa, occup, plot
>>>:line 3 is repeated for the other spin and for all orbitals specified above by config.
>>>: the last two lines must be
****
****
optional inserted as line 2 when “Watson” has been specified in line 1: free format
R-wat, Q-wat
The quantum numbers are defined as follows (see e.g. Liberman et al 65):
Spin quantum number: s = +1 or s = −1
Orbital quantum number j = l + s/2
Relativistic quantum number κ = −s(j + 1/2)
lattice vectors) for a mesh yourself. If inversion symmetry is not present, it will be added auto-
matically unless you specified the “-so” switch (for magnetic cases with spin-orbit coupling). The
k-mesh is then created with this additional symmetry. If symmetry permits, it further asks whether
or not the k-mesh should be shifted away from high symmetry directions. The file [Link] is
used in lapw1 and [Link] is used in tetra and lapw2, if the EF switch is set to TETRA, i.e.
the tetrahedron method for the k-space integration is used. For the format of the [Link] see
page 94.
6.5.1 Execution
With the switch -so it uses a file [Link] (usually generated by symmetso) instead of
[Link] and does not add inversion symmetry.
6.6.1 Execution
The program dstart is executed by invoking the command:
With the switch -fft dstart will terminate after case.in0 std has been created. The switch
-super will produce new [Link] instead of [Link], which is necessary for
charge extrapolation (clmextrapol lapw). -lcore produces [Link] from the radial core
densities [Link].
NPT number of r-mesh points in atomic density (should be the same as in LSTART)
LMAX2 max l in LM expansion
NCOM number of LM terms in density
82 CHAPTER 6. INITIALIZATION
7 Programs for running an SCF cycle
Contents
7.1 LAPW0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.2 ORB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
7.3 LAPW1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
7.4 LAPWSO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.5 LAPW2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.6 SUMPARA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.7 LAPWDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.8 LCORE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.9 MIXER . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
In sections 7.1-7.9 we describe the main programs to run an SCF cycle as illustrated in figure 4.1.
83
84 CHAPTER 7. SCF CYCLE
(the linear equations are solved with modified LINPACK routines). In the interstitial region a 3-
dimensional fast Fourier transformation (FFT) is used.
The total potential V is obtained by summation of the Coulomb VC and exchange-correlation po-
tentials Vxc .
In order to find the contribution from the plane wave representation to the Hamilton matrix el-
ements we reanalyze the Fourier series in such a way that the new series represents a potential
which is zero inside the spheres but keeps the original value in the interstitial region and this series
is put into [Link].
The contribution to the total energy which involves integrals of the form ρ ∗ V is calculated accord-
ing to the formalism of Weinert et al (82).
The Hellmann-Feynman force contribution to the total force is also calculated (Yu et al 91).
Finally, the electric field gradient (EFG) is calculated in case you have an L=2 term in the density
expansion. The EFG tensor is given in both, the “local-rotation-matrix” coordinate system, and
then diagonalized. The resulting eigenvectors of this rotation are given by columns.
For surface calculations the total and electrostatic potential at z=0 and z=0.5 is calculated and can
be used as energy-zero for the determination of the workfunction. (It is assumed that the middle
of your vacuum region is either at z=0 or z=0.5).
7.1.1 Execution
The following parameters are used (they are collected in file [Link], but usually need not to
be changed:
7.1.3 Input
The input is very simple. It is generated automatically by init lapw, and needs to be changed
only if a different exchange-correlation potential should be used:
switch
TOT total energy contributions and total potential calculated
KXC total energy contributions and total potential calculated. In addition the
kinetic energy contribution as well as the XC-energy will be printed.
POT total potential is calculated, but not the total energy
MULT multipole moments calculated only
COUL Coulomb potential calculated only
EXCH exchange correlation potential calculated only
NOTE: MULT, COUL, and EXCH are for testing only, whereas POT,
saves some CPU time if total energy is not needed
indxc index to specify type of exchange and correlation potential. Supported
options (for more options see the SRC lapw0/vxclm2.f subroutine) in-
clude:
5 Perdew and Wang 92, reparameterization of Ceperly-Alder data, the
recommended LDA option
13 Generalized Gradient approximation (PBE, Perdew-Burke-Ernzerhof
96)
11 Generalized Gradient approximation (Wu-Cohen 2006, Tran et al. 2007)
19 Generalized Gradient approximation (PBEsol, Perdew 2008)
20 Generalized Gradient approximation (AM05, Mattsson 2008)
12 Meta-GGA PKZB (Perdew et al. 1999). In order to generate the
requiered [Link]* files, you need [Link] vresp (cp
$WIENROOT/SRC templates/[Link] vresp [Link] vresp
and run one scf cycle with PBE (indxc=13) after creation of
[Link] vresp. Only afterwards change indxc to 12. In addi-
tion you must use very large IFFT parameters, otherwise it might be
numerically unstable.
27 Meta-GGA TPSS (Tao et al. 2003). At present the “best” meta-GGA.
(See also the option above.) Note: Only EXC is obtained that way, VXC
of standard PBE is used.
28 modified Becke-Johnson (mBJ-LDA) potential VXC (Tran and Blaha
2009). Uses a mBJ exchange potential + LDA correlation potential and
yields gaps in very good agreement with experiment. The xc-energy
EXC is from LDA. For detailed usage see chapter about mBJ calcula-
tions (4.5.8).
50 calculates the average of ∇ρ(r)/ρ(r) over the unit cell. This is used for
the mBJ potential mentioned above.
FFTopt IFFT optional keyword, which lets you define the IFFTx mesh and an en-
hancement factor in the next line (necessary for runeece lapw)
LUSE optional l-max value for the angular grid used in xcpot1. For stan-
dard LDA/GGA the recommended value is max L value of LM-list in
case.in2 + 2; for EECE one should use a better, antialiased grid, thus a
large negative LUSE-value is recommended (and set automatically by
runeece lapw)
line 3: free format (must be omitted when IFFT is not specified above)
IFFTx, IFFTy, IFFTz, IFFTfactor, iprint
The following line is optional and can be omitted. It is used to introduce an electric field via
a zig-zag potential (see [Link] et al. 2001):
line 4: free format
IFIELD, EFIELD, WFIELD
[Link]́k
Inst. of Physics, [Link], Prague, Czeck Republic
email: novakp@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
orb calculates the orbital dependent potentials, i.e. potentials which are nonzero in the atomic
spheres only and depend on the orbital state numbers l, m. In the present version the potential is
assumed to be independent of the radius vector and needs the density matrix calculated in lapwdm.
Four different potentials are implemented in this package:
7.2. ORB 87
I LDA+U. There are three variants of this method, two of them are discussed in Novák et al.
2001
1. LDA+U(SIC) - introduced by Anisimov et al. 1993, with an approximate correction for
the self-interaction correction. This is probably best suited for strongly correlated sys-
tems and for a full potential method we recommend to use an “effective” Uef f = U − J;
setting J = 0.
2. LDA+U(AMF) - introduced by Czyzyk and Sawatzky 1994 as ’Around the Mean Field’
method. (In Novák et al. 2001 it is denoted as LDA+U(DFT)). This version is (probably)
more suitable for metallic or less strongly correlated systems.
3. LDA+U(HMF) - in addition the Hubbard model in the mean field approximation, as
introduced by Anisimov et al. 1991 is also implemented. Note, however, that it is to be
used with the LDA (not LSDA) exchange-correlation potential in spin polarized calcu-
lations!
All variants are implemented in the rotationally invariant way (Liechtenstein et al. 1995). If
LDA+U is used in an unrestricted, general way, it introduces an orbital field in the calcu-
lation (in analogy to the exchange field in spin-polarized calculations, but it interacts with
the orbital, instead of spin momentum). The presence of such an orbital field may lower the
symmetry. In particular the complex version of LAPW1 must be used. Care is needed when
dealing with the LDA+U orbital field. It may be quite large, and without specifying its direc-
tion it may fluctuate, leading to oscillations of scf procedure or/and to false solutions. It is
therefore necessary to use it in combination with the spin-orbit coupling, preferably running
first LSDA+(s-o) and then slowly switching on the LDA+U orbital field. If the LDA+U orbital
polarization is not needed, it is sufficient to run real version of LAPW1, which then automat-
ically puts the orbital field equal to zero. For systems without the center of inversion, when
LAPW1 must be complex, an extra averaging of the LDA+U potential is necessary.
I Orbital polarization. The additional potential has the form (Brooks 1985, Eriksson et al. 1989):
VOP = cOP < Lz > lz (7.1)
where cOP is the orbital polarization parameter, < Lz > is projection of the orbital momen-
tum on the magnetization direction and lz is single electron orbital momentum component z
parallel to M~.
I Exact exchange and Hybrid methods: see Tran et al. 2006 and 4.5.7
I Interaction with the external magnetic field. In this case the additional potential has a simple
form:
VBext = µB B~ ext (~l + 2~s). (7.2)
The interaction with the electronic spin is taken into account by shifting the spin up and spin
down exchange correlation potentials in LAPW0 by the energy +µB Bext − µB Bext , respec-
tively. The interaction of Bext with spin could be as well calculated using the ’Fixed spin
moment’ method. For an interaction with the orbital momentum it is necessary to specify the
atoms and angular momentum numbers for which this interaction will be considered. Cau-
tion is needed when considering interaction of the orbital momentum with Bext in metallic
or metallic-like systems. For the analysis see the paper by Hirst 1997
In all cases the resulting potential for a given atom and orbital number l is a Hermitian, (2l +
1)x(2l + 1) matrix. In general this matrix is complex, but in special cases it may be real.
For more information see also section 4.5.6.
7.2.1 Execution
7.2.3 Input
Since this program can handle three different cases, examples and descriptions of [Link] for
all cases are given below:
line 2: (A5,f8.2)
mixmod,amix
mixmod PRATT or BROYD (should not be changed, see MIXER for more infor-
mation)
amix coefficient for the Pratt mixing of Vorb
This option is now only used for testing. The mixing should be set to
PRATT, 1.0
U(li,i), Coulomb and exchange parameters, U and J, for LDA+U in Ry for atom
J(li,i) type i and orbital number li. We recommend to use Uef f only.
7.2. ORB 89
Example of the input file for NiO (LDA+U included for two inequivalent Ni atoms that have in-
dexes 1 and 2 in the structure file):
Ncalc(i)
1 [Link]. parameters are calculated ab-initio
0 [Link]. parameters are read from input
pop(li,i) OP parameter in Ry
~
Example of the input file for NiO (total < Lz > used in (1), OP parameters calculated ab-initio, M
along [001]):
lapw1 sets up the Hamiltonian and the overlap matrix (Koelling and Arbman 75) and finds by
diagonalization eigenvalues and eigenvectors which are written to [Link]. Besides the
standard LAPW basis set, also the APW+lo method (see Sjöstedt et al 2000, Madsen et al. 2001) is
supported and the basis sets can be mixed for maximal efficiency. If the file [Link] exists (i.e.
non-spherical terms in the potential), a full-potential calculation is performed.
For structures without inversion symmetry, where the hamilton and overlap matrix elements are
complex numbers, the corresponding program version lapw1c must be used in connection with
lapw2c.
Since usually the diagonalization is the most time consuming part of the calculations, two options
exist here. In WIEN2k we include highly optimized modifications of LAPACK routines. We call all
these routines “full diagonalization”, but we also provide an option to do an “iterative diagonal-
ization” using a new preconditioning of a block-Davidson method (see Singh 89 and Blaha et al.
09). The scheme uses an old eigenvector from the previous scf-iteration, and produces approximate
(but usually still highly accurate) eigenvalues/vectors. The preconditioner (inverse of (H − λS)
can be calculated at the first iterative step (which will therefore take longer than subsequent iter-
ative steps), stored on disk ([Link]) and reused in all subsequent scf-iterations (until
the next “full” diagonalization or when it is recreated (x lapw1 -it -noHinv0)). Usually this is
the fastest scheme, but storage of [Link] can be large (and slow when you have a
slow network) and when the Hamiltonian changes too much, performance may degrade. Alter-
natively, the preconditioner can be recalculated all the time (x lapw1 -it -noHinv). Depending on
the ratio of matrix size to number of eigenvalues (cpu time increases linearly with the number
of eigenvalues, but a sufficiently large number is necessary to ensure convergence) a significant
speedup compared to “full” diagonalization (LAPACK) can be reached. Iterative diagonalization
is activated with the -it switch in x lapw1 -it or run lapw -it. Often the preconditioner
is so efficient, that it does not need to be recalculated, even within a structural optimization and
one can use min lapw -j ‘‘run lapw -I -fc 1 -it’’. In some cases it is preferable to use
min lapw -j ‘‘run lapw -I -fc 1 -it1’’, which will recreate [Link], or do
not store these files at all using min lapw -j ‘‘run lapw -I -fc 1 -it -noHinv ’’
Parallel execution (fine grain and on the k-point level) is also possible and is described in detail in
Sec. 5.5.
The switch -nohns skips the calculation of the nonspherical matrix elements inside the sphere.
This may be used to save computer time during the first scf cycles.
7.3. LAPW1 91
7.3.1 Execution
In cases without inversion symmetry, the default input filename is case.in1c. For 2-window (not
recommended) semi-core calculations the [Link] file uses a case.in1s file and creates the
files case.output1s and [Link]. For the spin-polarized case lapw1 is called twice with
[Link] and [Link]. To all relevant files the keywords “up“ or “dn“ are appended
(see the fcc Ni test case in the WIEN2k package).
LMAX highest l+1 in basis function inside sphere (consistent with input in case.in1)
LMMX number of LM terms in potential (should be at least NCOM-1)
LOMAX highest l for local orbital basis (consistent with input in case.in1)
NGAU number of Gaunt coefficients for the non-spherical contributions to the matrix
elements
NMATMAX maximum size of H,S-matrix (defines size of program, should be chosen accord-
ing to the memory of your hardware, see chapter 11.2.2!)
NRAD number of radial mesh points
NSLMAX highest l+1 in basis functions for non-muffin-tin matrix elements (consistent with
input in case.in1).If set larger than 5, parameter MAXDIM (modules.F) and LO-
MAX=8, P(10,10) (gaunt2.f) must also be increased.
NSYM order of point group
NUME maximum number of energy eigenvalues per k-point
NVEC1 defines the largest triple of integers which define reciprocal
NVEC2 K-vectors when multiplied with the reciprocal Bravais matrix
NVEC3
NLOAT max number of LOs for one l-quantum number
restrict output for mpi-jobs, limits the number of case.output1 X proc XXX files to “re-
strict output”
7.3.3 Input
Below a sample input is shown for T iO2 (rutile), one of the test cases provided in the WIEN2k
package. The input file is written automatically by LSTART, but was modified to set APW only for
Ti-3d and O-2p orbitals.
switch WFFIL standard option, writes wave functions to file [Link] (needed
in lapw2)
SUPWF suppresses wave function calculation (faster for testing eigenvalues
only)
WFPRI prints eigenvectors to case.output1 and writes [Link] (pro-
duces long outputs!)
EF optional input. If “EF=” key is present, lapw1 reads EF and sets all
default energy parameters (0.3) to “EF-0.2” Ry.
rkmax Rmt ∗ Kmax determines matrix size (convergence), where Kmax is the
plane wave cut-off, Rmt is the smallest of all atomic sphere radii. Usu-
ally this value should be between 5 and 9 (APW+lo) or 6 - 10. (LAPW-
2
basis) (Kmax would be the plane wave cut-off parameter in Ry used
in pseudopotential calculations.) Note that d (f) wavefunctions con-
verge slower than s and p. For systems including hydrogen with short
bondlength and thus a very small Rmt (e.g. 0.7 a.u.), RKmax = 3 might
already be reasonable, but convergence must be checked for a new type
of system.
Note, that the actual matrix size is written on case.scf1. It is determined
by whatever is smaller, the plane wave cut-off (specified with RKmax)
or the maximum matrix dimension NMATMAX, (see previous section).
lmax maximum l value for partial waves used inside atomic spheres (should
be between 8 and 12)
lnsmax maximum l value for partial waves used in the computation of non-
muffin-tin matrix elements (lnsmax=4 is quite good)
Etrial default energy used for all El to obtain ul (r, El ) as regular solution of
the radial Schrödinger equation [used in equ.2.4,2.7] (see figure 7.1).
ndiff number of exceptions (specified in the next ndiff lines)
Napw 0 ... use LAPW basis, 1 ... use APW-basis for all “global” l values of this
atom. We recommend to use LAPW here.
line 4: format(I2,2F10.5,A4)
l, El, de, switch, NAPWL
l l of partial wave
El El for L=l
de energy increment
7.3. LAPW1 93
u l(r,E) E
E bottom
E top
El El
E top E bottom
r
RMT DOS
>>>:line 4 is repeated ndiff times (see line 3) for each exception. If the same l value is specified
twice, local orbitals are added to the (L)APW basis. The first energy (E1 ) is used for the usual
LAPW’s and the second energy (E2 ) for the LOs, which are formed according to (see equ.
2.7): uE1 + u̇E1 + uE2 .
Note: The default energy parameters (0.30) are replaced by an energy “EF − 0.2” if the EF-switch
was read before. Please read also the comments about run lapw in section 5.1.3. In addition, you may
want to change the automatically created input and add d- or f-local orbitals to reduce the linearization
error (e.g. in late transition metals you could put E3d at 0.0 and 1.0 Ry) or s, p, d, and/or f-LOs at
very high energy (e.g. 2.0 - 3.0 Ry) to better describe unoccupied states.
>>>:lines 3 and 4 are repeated for each non equivalent atom
line 5: format (20x,i1,2f10.1,i6)
unit-number, Emin, Emax nband
94 CHAPTER 7. SCF CYCLE
unit- file number from which the k-vectors in the irreducible wedge of the
number Brillouin zone are read. Default is 4, for which the corresponding infor-
mation is contained in [Link] (generated by KGEN). Should not
be changed.
EMIN, energy window in which eigenvalues shall be searched (overrides set-
EMAX ting in [Link]. A small window saves computer time, but it also
limits the energy range for the DOS calculation of unoccupied states.
nband number of eigenvalues calculated with iterative diagonalization. Set
automatically to nband = ne ∗ 2.0 + 5 in lstart. Larger values will
lead to more cpu-time. (Optional input)
line 6: free format; optional input line, but necessary if k-vectors are read from unit 5
spro limit, lambda iter
spro limit limit for detection of linear dependency for iterative diagonalization
(optional input), typical around 1.d-15)
lambda iter optional λ value for preconditioner of iterative diagonalization (see
above). By default we use λ = 0, but in some cases convergence can
be improved by a small (around 1.0) positive or negative λ
line 7: format (A10,4I10,3F5.2); (only when unit-number=5, not recommended, use unit 4 and
[Link])
name, ix,iy,iz, idv, weight
>>>: line 7 is repeated for each k-vector in the IBZ. The utility program kgen (see section 6.5)
provides a list of such vectors (on a tetrahedral mesh) in [Link].
>>>: the last line must be END
7.4. LAPWSO 95
I [Link]
I [Link]
I [Link]
I case.in2c
I case.in1
We recommend to use initso (see Sec.5.2.17) which helps you together with symmetso (see
Sec.9.1) to setup spinorbit calculations.
Note: SO eigenvectors are complex and thus lapw2c must be used in a selfconsistent calculation.
7.4.1 Execution
where here -up indicates a spin-polarized calculation (no “-dn” is needed, since spin-orbit will mix
spin-up and dn states in one calculation).
FLMAX constant = 3
LMAX highest l of wave function inside sphere (consistent with lapw1)
LABC highest l of wave function inside sphere where SO is considered
LOMAX max l for local orbital basis
NRAD number of radial mesh points
NLOAT number of local orbitals
96 CHAPTER 7. SCF CYCLE
7.4.3 Input
A sample input for lapwso is given below. It will be generated automatically by initso
line 1: format(A5)
switch
Emin minimum energy for which the output eigenvectors and eigenenergies
will be printed (Ry)
Emax maximum energy
h,k,l vector describing the direction of magnetization. For R lattice use h,k,l
in rhombohedral coordinates (not in hexagonal)
STOP calculation stops if not both Etop and Ebottom are found (especially use-
ful for semi-core states)
>>>: line 6 must be repeated “nlr” times (or should be omitted if nlr=0).
line 7: free format
noff, (iatoff(i),i=1,noff)
noff number of atoms for which SO is switched off (for “light” elements,
saves time)
iatoff atom-numbers
7.5.1 Execution
x lapw2 [-c -up|dn -p -so -qtl -fermi -efg -band -eece -vresp
-help files -emin X -all X Y] or
lapw2 [Link] [proc#] or lapw2c [Link] [proc#]
where proc# is the i-th processor number in case of parallel execution (see Fig. 5.2). The -so switch
sets -c automatically.
For complex calculations case.in2c is used. For a spin-polarized case see the fcc Ni test case in
the WIEN2k package.
7.5.3 Input
A sample input for lapw2 is listed below, it is generated automatically by the programs lstart
and symmetry.
line 1: format(2A5)
switch, EECE
switch TOT total valence charge density expansion inside and outside spheres
FOR same as TOT, but in addition a “Pulay” force contribution is calculated
(this option costs extra computing time and thus should be performed
only at the final scf cycles, see run lapw script in sec. 5.1.3)
QTL partial charges only (generates file [Link] for DOS calculations), set
automatically by switch -qtl
EFG computes decomposition of electric field gradient (EFG), contributions
from inside spheres (the total EFG is computed in lapw0), set automat-
ically by switch -efg.
ALM this generates two files, [Link] and [Link], where the
radial wavefunctions and the Alm , Blm , Clm coefficients of the wave-
function inside spheres are listed. The file [Link] can get very
big.
CLM CLM charge density coefficients only
FERMI Fermi energy only, this produces weight files for parallel execution
and for the optics and lapwdm package, set automatically by switch
-fermi.
>>>: TOT and FOR are the standard options, QTL is used for density of states
(or energy bandstructure) calculations, EFG for analysis, while FOURI,
CLM are for testing only.
EECE if set to “EECE”, calculates the density for specified atoms and angu-
lar momentum only. Used for exact-exchange or hybrid-calculations,
usually set automatically by runsp lapw -eece
emin lower energy cut-off for defining the range of occupied states, can be
set termporarely to “X” by switch -emin X or -all X Y
7.5. LAPW2 99
line 3: format(a5,f10.5)
efmod, eval
nat rho number of atoms for which a specific density should be calculated
jatom rho index of atom for which a specific density should be calculated
100 CHAPTER 7. SCF CYCLE
l rho angular momentum l-value for which a specific density should be cal-
culated
Symmetry LM combinations
23 0 0, 4 0, 4 4, 6 0, 6 4,-3 2, 6 2, 6 6,-7 2,-7 6, 8 0, 8 4, 8 8,-9 2,-9 6,-9 4,-9 8,10 0, 10 4,10 8, 10 2, 10 6, 10 10
M3 0 0, 4 0, 4 4, 6 0, 6 4, 6 2, 6 6, 8 0, 8 4, 8 8,10 0, 10 4,10 8, 10 2, 10 6, 10 10
432 0 0, 4 0, 4 4, 6 0, 6 4, 8 0, 8 4, 8 8,-9 4,-9 8,10 0, 10 4,10 8
-43M 0 0, 4 0, 4 4, 6 0, 6 4,-3 2,-7 2,-7 6, 8 0, 8 4, 8 8,-9 2,-9 6,10 0, 10 4,10 8
M3M 0 0, 4 0, 4 4, 6 0, 6 4, 8 0, 8 4, 8 8,10 0, 10 4,10 8
Table 7.41: LM combinations of “Cubic groups” (3k(111)) direction, requires “positive atomic in-
dex” in [Link]. Terms that should be combined (Kara and Kurki-Suonio 81) must follow one
another.
Table 7.42: LM combination and local coordinate system of “non-cubic groups” (requires “negative
atomic index” in [Link])
line 6: A4
reclist FILE writes list of K-vectors into file [Link] or reuses this list if
the file exists. The saved list will be recalculated whenever GMAX, or a
lattice parameter has been changed.
7.6.1 Execution
where # of proc is the numbers of parallel processors used. It is usually called automatically
from lapw2para or x lapw2 -p.
The following parameters are listend in file [Link], but usually they need not to be modified:
NCOM number of LM terms in density
NRAD number of radial mesh points
NSYM order of point group
102 CHAPTER 7. SCF CYCLE
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
lapwdm calculates the density matrix needed for the orbital dependent potentials generated in orb.
Optionally it also provides orbital moments, orbital and dipolar contributions to the hyperfine field
(only for the specified atoms AND orbitals). It calculates the average value of the operator X which
behaves in the same way as the spin-orbit coupling operator: it must be nonzero only within the
atomic spheres and can be written as a product of two operators - radial and angular:
X = Xr (r) ∗ Xls (~l, ~s)
Xr (r) and Xls (~l, ~s) are determined by RINDEX and LSINDEX in the input as described below:
I RINDEX=0 LSINDEX=0: the density matrix is calculated (this is needed for LDA+U calcula-
tions)
I RINDEX=1 LSINDEX=1: <X> is number of electrons inside the atomic sphere (for test)
I RINDEX=2 LSINDEX=1: <X> is the < 1/r3 > expectation value inside the atomic sphere
I RINDEX=1 LSINDEX=2: <X> is the projection of the electronic spin inside the atomic sphere
(must be multiplied by g=2 to get the spin moment)
I RINDEX=1 LSINDEX=3: <X> is the projection of the orbital moment inside the atomic
sphere (in case of SO-calculations WITHOUT LDA+U)
I RINDEX=3 LSINDEX=3: <X> is the orbital part of the hyperfine field at the nucleus (for a
converged calculation at the very end)
I RINDEX=3 LSINDEX=5: <X> is the spin dipolar part of the hyperfine field at the nucleus
(for a converged calculation at the very end)
To use the different operators, set the appropriate input. More information and extentions to op-
erators of similar behavior may be obtained directly from P. Novák (2006). (RINDEX=3 includes
now an approximation to the relativistic mass enhancement and LSINDEX=5 includes nondiagonal
terms)
lapwdm needs the occupation numbers, which are calculated in lapw2. Note: You must not use
ROOT in case.in2 for that purpose.
7.7.1 Execution
The program lapwdm is executed by invoking the command:
FLMAX constant = 3
LMAX highest l of wave function inside sphere (consistent with lapw1)
LABC highest l of wave function inside sphere where SO is considered
LOMAX max l for local orbital basis
NRAD number of radial mesh points
7.7.3 Input
iatom index of atom for which the density matrix should be calculated
nl number of l-values for which the density matrix should be calculated
l l-values for which the density matrix should be calculated
7.8.1 Execution
The program lcore is executed by invoking the command:
104 CHAPTER 7. SCF CYCLE
For the spin-polarized case see fcc Ni on the distribution tape. If [Link] and [Link]
are present for spin-polarized calculations, different core-occupation (“open core” approximation
for 4f states or spin-polarized core-holes) for both spins are possible.
7.8.3 Input
>>>: line 2 is repeated for each orbital (nrorb times; see line 1)
>>>: line 1 and 2 are repeated for each inequivalent atom. Atoms without core states (e.g. H or
Li) must still include a 1s orbital, but with occupation zero.
line 3: free format
7.9. MIXER 105
2. a Multi-Secant mixing scheme contributed by L. Marks (see Marks and Luke 2008), in which
all the expansion coefficients of the density from several preceding iterations (usually 6-10)
are utilized to calculate an optimal mixing fraction for each coefficient in each iteration. It is
very robust and stable (works nicely also for magnetic systems with 3d or 4f states at EF, only
for ill-conditioned single-atom calculations you can break it) and usually converges at least
30 % faster than the old BROYD scheme.
3. Two new variants on the Multi-Secant method including a rank-one update (see Marks 2011)
which appear to be 5-10% faster and equally robust.
At the outset of a new calculation (for any changed computational parameter such as k-mesh, ma-
trix size, lattice constant etc.), any existing [Link] files must be deleted (since the iterative
history which they contain refers to a “different“ incompatible calculation).
If the file [Link] old can not be found by mixer, a “PRATT-mixing“ with mixing factor
(greed) 1.0 is done.
Note: a [Link] file must always be present, since the LM values and the K-vectors are read from
this file.
The total energy and the atomic forces are computed in mixer by reading the [Link] file and
adding the various contributions computed in preceding steps of the last iteration. Therefore
[Link] must not contain a certain “iteration-number” more than once and the number of it-
erations in the scf file must not be greater than 999.
For LDA+U calculations [Link]/dn and for hybrid-DFT (switch -eece) [Link]/dn
files will be included in the mixing procedure.
With the new mode MSR1a (or MSECa) (contributed by L. Marks) atomic positions will also be
mixed and thus optimized. This scheme can (unfortunately not in all cases) be a facter or 2-3 faster
then the traditional optimization using min lapw.
7.9.1 Execution
The program mixer is executed by invoking the command:
7.9.3 Input
Below a sample input (written automatically by lstart) is provided for T iO2 (rutile), one of the
test cases provided with the WIEN2k package.
line 1: (A5,*)
switch, bgch, norm
bgch Background charge for charged cells (+1 for additional electron, -1 for
core hole, if not neutralized by additional valence electron)
greed mixing greed Q. Essential for Pratt, rather less important for MSEC1. In
the first iteration using Broyden’s scheme: Q is automatically reduced
by the program depending on the average charge distance :DIS andthe
relative improvement in the last cycle. In case that the scf cycle fails
due to large charge fluctuations, this can be further reduced but this
can lead to stagnation. One should rarely reduce this below 0.05.)
nuse For all Multisecant methods: Only nuse steps are used during modified
broyden (this value has some influence on the optimal convergence.
Usually 6-10 seems reasonable and 8 is recommended).
108 CHAPTER 7. SCF CYCLE
8 Programs for analysis, calculation of
properties, and geometry
optimization
Contents
8.1 TETRA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8.2 QTL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
8.3 SPAGHETTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
8.4 IRREP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
8.5 LAPW3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.6 LAPW5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.7 AIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
8.8 LAPW7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
8.9 FILTVEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.10 XSPEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
8.11 TELNES3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
8.12 BROADENING . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
8.13 OPTIMIZE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.14 ELAST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
8.15 MINI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
8.16 OPTIC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
8.17 JOINT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
8.18 KRAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
8.19 DIPAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
8.20 FSGEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
109
110 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
8.1.1 Execution
8.1.3 Input
The required input file [Link] can optionally be created using the w2webinterface or the
configure int lapw script (see 5.2.7).
An example is given below:
emin, specifies the energy mesh (in Ry) where the DOS is calculated. (emin
delta, should be set slightly below the lowest valence band; emax will be
emax checked against the lowest energy of the highest band in [Link],
and set to the minimum of these two values; delta is the energy incre-
ment.
broad Gauss-broadening factor. Must be greater than delta to have any effect.
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 [Link] (see sec. 4.3); the respective description can be
found in the header of [Link].
description text used for further identification.
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
qtl creates the input for calculating total and projected density of states of selected atoms and se-
lected 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 p1/2 , p3/2 , d3/2 , d5/2 ,
f5/2 , f7/2 , 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 when-
ever 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 [Link] [up/dn]. For the DOS calculation the file [Link] [up/dn] is
created in which the ordering of partial charges is given. Please note, that in contrast to [Link]
[up/dn] from x lapw2 -qtl the total partial charge of an atom is NOT multiplied with its “mul-
tiplicity” and contains only the sum of the requested l,m terms (eg. s,p,d) and thus not all contri-
butions. Also the interstital charge will usually be NOT correct.
8.2.1 Execution
8.2.2 Input
A sample input (a default is created automatically during init lapw for [Link] is given be-
low.
Lines starting from line 3 are repeated for each selected atom. Line 5 only appears when calculation
in new coordinate system is required (loro 6= 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 [Link].
Indices of selected atoms, as well as the orbital numbers, must form an ascending sequence.
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 [Link]$i and
ordered according to increasing l. The unitary transformation matrix must rotate from the stan-
dard lms -basis to the desired one. A few examples (e.g. jjz , lms , or eg − t2g ) are supplied with the
code in $WIENROOT/SRC templates/[Link] * and must be copied to [Link]$i . For less
114 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
8.2.3 Output
The results in file [Link][up/dn] are written in the same format as lapw2 file [Link][up/dn] and
thus they may be directly used by tetra.
The data for the interstital DOS correspond to n = nat + 1 (nat is number of atom types). The
ordering of densities for all selected atoms is summarized in the file [Link][up/dn]. The qtltext
file that corresponds to the input data given above is:
The output for the population matrix integrated over energy is written to [Link] [up/dn] that
has the same format as analogous file calculated by lapwdm.
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 space-
groups, because the corresponding group-theory has not been implemented.)
The presence of “incompatible” [Link] or [Link] 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 o Tasks o Bandstructure” from w2web.
8.3.1 Execution
The -p switch directs spaghetti to use the case.output1 * files of a k-point parallel lapw1.
8.3.2 Input
test line must start with ’###’. Begin of figure description. This tests also if
you use the new input (different from WIEN97 or early WIEN2k ver-
sions)
eincr, mtick
eincr energy increment where y-axis labels are printed (major ticks)
mtick number of minor ticks of y-axis
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 [Link] file. If set
to 999., Ef is not plotted (and iunits=2 cannot be used)
nband1, lower and upper band index for bands which should show “band-
nband2 character plotting“ (if [Link] is present and the proper switch is
set, see below). In addition the corresponding x and y coordinates are
written to file [Link] ene (which can be used for plotting
with an external xy-plotting program).
jatom If a [Link] file is present, jatom indicates the atom whose charac-
ter (selected by jcol) is used for “band-character plotting“ (dots are re-
placed by circles with radii proportional to the corresponding weight,
requires ilin=0,2,3). If set to zero or if [Link] 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 [Link]. (ignored for jatom=0). The description can be
found in the header of [Link].
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).
Clas Persson
Condensed Matter Theory Group,Department of Physics,
University of Uppsala, Sweden
email: [Link]@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
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
[Link] and [Link]. For nonmagnetic SO calculations you must set IPR=1 in
[Link].
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.
8.4.1 Execution
The program irrep is executed by invoking the command:
118 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
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
8.5.1 Execution
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 ([Link]) or difference den-
sities (i.e. crystalline minus superposed atomic densities) using the additional file ([Link]).
In spinpolarized cases one can produce up-, dn- and total densities but also spin densities (differ-
ence up-dn). It is also possible to plot total densities ([Link]), Coulomb ([Link]),
exchange-correlation (case.r2v) or total ([Link]) potentials, but in those cases the file
[Link] has to be edited and you must replace [Link] by the respective filename. The
file [Link] contains in the first line
8.6. LAPW5 119
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 o Tasks o Electron density plots” from
w2web, see the TiC example in Fig.3.6 .
8.6.1 Execution
8.6.3 Input
An example is given below. You may want to use XCRYSDEN by [Link] to generate this file (see
sect. 9.23.2).
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 con-
ventional cell are used)
nxsh, specifies the number of nearest neighbor cells (in x,y,z direction) where
nysh, atomic positions are generated (needs to be increased for very large plot
nzsh sections, otherwise some “atoms” are not found in the plot)
switch RHO charge (or potential) plots, no atomic density is used (regular case)
DIFF difference density plot (crystalline - superposed atomic densities),
needs file [Link] (which is generated with LSTART, see section
6.4)
OVER superposition of atomic densities, needs file [Link]
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 densi-
ties).
In order to plot total densities or potentials (see cnorm as above) you have to create [Link] us-
ing x lapw5 -d, then edit [Link] and insert proper filenames ([Link], case.r2v,
[Link], [Link]) for units 9 and 11, and finally run lapw5 [Link].
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
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 [Link] (L-marks@[Link]).
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/[Link].
Basically it performs two different tasks, namely searching for “critical points” (CP) and/or deter-
mination 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.
8.7.1 Execution
8.7.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:
line 1: A4
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 “sec-
ond atom” has iatom=3 (Do not use simply the atom-numbers from
[Link])
line 3: A4
specifies the number of nearest neighbor cells (in x,y,z direction) where
atomic positions are generated.
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.
8.7. AIM 123
For convenience run extractaim lapw [Link] (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 Å, e/Å3 , ... (instead of bohr).
line 1: A4
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 “sec-
ond atom” has iatom=3 (Do not use simply the atom-numbers from
[Link])
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=π and use only the “irreducible” part, i.e.
when you have a mirror plane normal to z (see [Link]), restrict
thmax to π/2.
specifies the number of nearest neighbor cells (in x,y,z direction) where
atomic positions are generated.
line 8: A4
line 9: A4
npt specifies number of points for radial integration outside the MT ( 30)
line 8: A4
Uwe Birkenheuer
Max-Planck-Institut für Physik komplexer Systeme
Nöthnitzer Str. 38, D-01187 Dresden, Germany
email: birken@[Link]
and
Birgit Adolph,
University of Toronto, T.O., Canada
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
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 [Link] (or in re-
duced (filtered) form in [Link] which can be obtained from [Link] 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 argu-
ment, or the complex wave function itself. For scalar-relativistic calculations both the large and
8.8. LAPW7 125
the small component of the wave functions can be generated (only one at a time). The wave func-
tions 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 [Link]) 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 process-
ing – e.g. for plotting the wave functions with some graphical tools such as gnuplot – in raw
format to [Link]. For quick inspection, a subset of this data is echoed to the standard out-
put file [Link] (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 Alm , Blm (and Clm ) to a binary file [Link]. 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.
8.8.1 Execution
With the -sel option lapw7 expects data from the reduced (filtered) wave function file
[Link], otherwise the standard wave function file [Link] is used. The reduced
vector file [Link] is assumed to resist in the current working directory, while the stan-
dard vector file [Link] (which may become quite large) is looked for in the WIEN scratch
directory. For details see [Link].
8.8.3 Input
A sample input is given below. It shows how to plot a set of wave functions on a 2-dim. grid.
126 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
- - - - - - - - - - - - - - - - - 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 - - - - - - - - - - - - - - - - -
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, 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 hblanki The axes of the n-dim. grid have to be mutual or-
thogonal.
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 in-
crement for the output echo. Again a zero incre-
ments 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 func-
tions.
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 eval-
uated.
ABS The absolute value of the wave functions is eval-
uated.
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.
128 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
Uwe Birkenheuer
Max-Planck-Institut für Physik komplexer Systeme
Nöthnitzer Str. 38, D-01187 Dresden, Germany
email: birken@[Link]
and
Birgit Adolph
University of Toronto, T.O., Canada
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
The program filtvec reduces the information stored in [Link] 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 individ-
ually 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 [Link]. The
information on all other wave functions in the original file is discarded. The structure of the gener-
ated [Link] file is identical to that of the original [Link] file. Hence, it should be
possible to use [Link] as substitutes for [Link] anywhere in the WIEN program
package. (This has only been tested for [Link] 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.
8.9.1 Execution
In accordance with the file handling for lapw1 and lapw7 the input vector file [Link]
is assumed to be located in the WIEN scratch directory, while the reduced output vector file
[Link] is written to the current working directory. See [Link] for details.
8.9.3 Input
Two examples are given below. The first uses global selection mode; the second individual 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 - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - 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 - - - - - - - - - - - - - - - - -
8.10. XSPEC 131
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 be-
low. It uses the partial charges in [Link]. 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 [Link], broadened spectra in the file [Link]. 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 [Link] (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 o Tasks o X-ray spectra” from w2web.
8.10.1 Execution
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.
132 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
initxspec This program generates the appropriate input file [Link], according to the dipole
selection rule, for the subsequent execution of the tetra program.
The program initxspec is executed by invoking the command:
initxspec [Link] 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 [Link] or x tetra [-up|-dn]
txspec This program calculates energy dependent dipole matrix elements. Theoretical X-ray spec-
tra are generated using the partial densities of states (in the case.dos1ev file) and multi-
plying them with the corresponding dipole matrix elements.
The program txspec is executed by invoking the command:
txspec [Link] 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 [Link] or x lorentz [-up|-dn]
The following dimensioning parameters are collected in the files [Link] of SRC txspec and
SRC lorentz:
8.10.3 Input
Two examples are given below; one for emission spectra and one for absorption spectra:
TITLE Title
Spectrum n l
K 1 0
LII,III 2 1
MV 3 2
Table 8.50: Quantum numbers of the core state involved in the x-ray spectra
SPLIT, split in eV between e.g. LII and LIII spectrum (compare with the re-
INT1, spective core eigenvalues), INT1 and INT2 specifies the relative inten-
INT2 sity between these spectra. Values of 0, 0.5, 0.5 give unshifted spectra.
EMIN, minimum energy, energy increment for spectrum, maximum energy; all
DE, energies are in eV and with respect to the Fermi level
EMAX
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
GAMMA0 broadening parameter for the life-time broadening of the core states.
Set GAMMA0 to zero to avoid lifetime broadening of the core states.
BANDRA
AUTO band ranges are determined AUTOmatically (default)
MAN band ranges have to be entered MANually
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
8.11. TELNES3 135
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 ([Link] user/textbooks/) also de-
scribes the formalism onto which TELNES3 is built in great detail.) This formalism allows calcu-
lation 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 (be-
cause the crystal is isotropic, or the sample is polycrystalline), the formula may be integrated over
4π, 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 [Link]
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.
8.11.1 Execution
Execution
8.11.2 Input
TELNES3 requires one input file - [Link]. We recommend using InnesGenTM 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 [Link] 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:
This first part of the file is not formatted and contains the following information:
136 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
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 [Link].
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
q line from the sample to the center of the detector aperture, these 2 lines will form an
angle of theta2x + theta2y mrad.
MODUS
m eg. : angles (default is energy)
8.11. TELNES3 137
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 [Link], 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 [Link] must be present from which
the dos will be calculated. If make dos=N, then either a file [Link] or a file [Link] con-
taining the (x)dos must exist. If make [Link]=N, a file [Link] containg the rotation matrices
must exist. If write [Link]=Y, a file [Link] is written. If write dos=Y, a file [Link] or
[Link] 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 α and convergence angle β allow scattering angles up to α +β and a correspond-
ing set of Q-vectors. This set (a disk of radius α + β) 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)N T points
on circle i, giving N R2 ∗ N T points in total.
138 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
L a logarithmic grid with N R circles. The distance between circles increases exponentially. There
are (2i − 1)N T points on circle i, and N R2 N T points in total. Circle i is at radius theta 0
e((i−1)dx) , where dx depends on N R, α and β.
1 a one dimensional logarithmic mesh; there are N R circles at exponential positions, and only one
point on each circle (so N R 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.
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 partic-
ular 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 spheri-
cal Bessel function jλ (q) 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 :
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 λ = 1 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 se-
lects them by the L-character of the final states. Possible settings for ‘type’ are (the orbital momen-
tum 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
8.11. TELNES3 139
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 λ = 1 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 [Link], which must be present for this type of calculation, and reexpand-
ing 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 [Link].)
CORE WAVEFUNCTION
filename (e.g. [Link])
Read the wave function of the initial state from file. (Default behavior is to calculate it instead.)
Read the radial wave functions of the final state from file. (Default behavior is to calculate it in-
stead.)
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 )
140 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
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.
This sequence can conveniently be executed using w2web by simply clicking one button after the
other.
8.11.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 [Link]. We list here all files
possibly used by TELNES3 (and listed in [Link]). 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.
Joachim Luitz
IAST Austria
wien2k@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
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.
8.12.1 Execution
Execution
8.12.2 Input
broadening needs one input file - [Link]. When running TELNES3 this input file is automat-
ically created from settings given in [Link].
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
8.13.1 Execution
The program optimize is executed by invoking the command:
8.14. ELAST 143
8.13.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.
Thomas Charpin
Lab. Geomateriaux de l’IPGP, Paris, France
(In September 2001 we received the sad notice that Thomas Charpin died in a
car accident).
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
This package calculates elastic constants for cubic crystals. It is described in detail by the author
in Charpin 2001. Please note, at [Link] user/unsupported you can
find a package Hex-elastic by Morteza Jamal, which should calculate elastic constants also for
hexagonal symmetry.
8.14.1 Execution
I init elast:
It prepares the whole calculation and should be run in a directory with a valid [Link]
and [Link] file. It creates the necessary subdirectories elast, elast/eos,
elast/tetra, elast/rhomb, elast/result, the templates for tetragonal and rhombo-
hedral distortion and initializes the calculations using init lapw.
I elast setup:
It should be run in the elast directory, generates the distorted struct-files and [Link],
[Link] and [Link]. These scripts must be adapted to your needs (spin-
polarization, convergence,...) and run. elast setup can be run several times (for different
distortions,...).
I ana elast:
Once all calculations are done, change into elastresult and run this script. The final
results are stored in elastresultoutputs.
I genetempl, setelast, anaelast:
These three small programs are called by the above scripts.
144 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
8.15.1 Execution
8.15.3 Input
Two examples are given below; one for a PORT geometry optimization, and one for molecular
dynamics using a NOSE thermostat:
line 1: format(a4,2f5.2)
DELTA(1- For PORT (and BFGS): Precondition parameters: rescales the gradient
3) and thus determines the size of the geometry steps
For NEWT/NEW1: x,y,z-delta parameters. Determines speed of mo-
tion. Good values must be found for each individual system. They de-
pend 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 [Link] and the forces :FGL from [Link]).
Whenever you change these DELTA(i) you must remove file [Link] !
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 1: format(a4,f5.2)
NOSF Nose-frequency
Claudia Ambrosch-Draxl
Atomistic Modelling and Design of Materials
University Leoben
A-8700 Leoben, AUSTRIA
email: cad@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
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.
8.16. OPTIC 147
The program optic generates the symmetrized squared momentum matrix elements
Mi =< n0~k|~
p.~ei | n~k >2
between all band combinations for each k-point given in the vector-file and stores them in
[Link]. For the orthogonal lattices the squared diagonal components can be found in the
file [Link] diag. For non-orthogonal systems all 6 components (Mj )∗ Mk 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 spinpolar-
ized 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 dichro-
ism) has been added by Lorenzo Pardini ([Link]@[Link]). Please cite Pardini et al. 2011
when using XMCD and check the paper for further details. In the case of the XMCD calcu-
lation, 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. L3 ) and
case.symmat2up (lower energy core state, eg. L1 ) for each k-point and every band. For K, L1 , and
M1 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:
I copy [Link] to [Link] (cp [Link] [Link]) and remove all the symmetry
operations but the identity;
I generate a k-mesh in the whole Brilouin zone (x kgen -so);
I change TOT to FERMI in case.in2c;
I set IPRINT=1 in [Link] to activate core-wavefunction output;
I for metallic systems, put TETRA with value 101;
I execute runsp lapw -so -s lapw1 -e lcore;
I run optic: x optic -c -so -up;
I 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.
8.16.1 Execution
I generate a fine k-mesh for the optics part: x kgen [-so (if [Link] has been
created by symmetso) ]
I change TOT to FERMI in case.in2c
I execute run[sp] lapw -so -s lapw1 -e lcore with this fine k-mesh
I run optic: x opticc -so [-up]
I run joint: x joint [-up]
I run kram: x kram [-up]
I cp [Link] [Link]
I cp [Link] [Link]
I cp [Link] [Link]
I x lapw2 -fermi -so -c
I cp [Link] [Link]
I cp [Link] [Link]
I x optic -so -up
I 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.
The following dimensioning parameters (listed in [Link] r and [Link] c) are used:
LMAX highest l+1 in basis function inside sphere (reducing LMAX to 4 or 5 may dra-
matically 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
8.16.3 Input
An example is given below:
---------------- top of file: [Link] --------------------
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 -------------------------
emin, absolute energy range (Ry) for which matrix elements should be calcu-
emax lated
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 [Link] file) for which XMCD should be
calculated
edge specify the edge: must be K, L1, L23, M1, M23 or M45
IATOMS List of NMODUS atoms for which the opt. matrix elements should be
calculated (see above).
8.17. JOINT 151
Claudia Ambrosch-Draxl
Atomistic Modelling and Design of Materials
University Leoben
A-8700 Leoben, AUSTRIA
email: cad@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
This program carries out the BZ integration using the momentum matrix elements [Link]
calculated before by optic. The interband or the intraband contributions to the imaginary part of
the dielectric tensor (2 ) can be computed. Alternatively, the DOS or the joint DOS can be derived.
The output in [Link] can be plotted with any xy-plotting package or opticplot lapw or
Curve lapw.
Warning: Negative values for 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 [Link] and [Link]: 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 L2,3 , M2,3 , and M4,5 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.
8.17.1 Execution
8.17.3 Input
b1, b2, lower, upper and (optional) upper-valence band-index (Setting b3 may
b3 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 4: optional line for XMCD, must be omitted for ‘‘normal’’ optic; free
format
E core1, lower and higher core energies (in Ry, get them using eg. “grep :2P
E core2 [Link]”)
broad core1, lifetime broadening (eV) of lower and higher core state
broad core2
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!
Claudia Ambrosch-Draxl
Atomistic Modelling and Design of Materials
University Leoben
A-8700 Leoben, AUSTRIA
email: cad@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
The Kramers-Kronig analysis is carried out for the actual number of columns contained in the
[Link][up|dn] file. For each real component its imaginary counterpart is created and vice
versa. All dielectric tensor components can be found in file [Link][up|dn]. The real
and imaginary parts of the optical conductivity (in 1015 /s) are written to file
[Link][up|dn]. In addition, file [Link] contains the real parts of the optical
conductivity (in 1/(Ωcm) and the absorption coefficients. The loss function is also calculated
([Link]), 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 1 ) 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 2 (in [Link] and [Link]) using addjoint-updn lapw
(this will produce [Link]) before calling (non-spinpolarized) x kram.
The 3 sumrules are also checked and written to [Link].
The output in [Link][up|dn] and [Link][up|dn] can be plotted with any
xy-plotting package, opticplot lapw or the ”OPTIC”-task in w2web.
154 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
8.18.1 Execution
8.18.3 Input
ESHIFT Energy shift (scissors operator) (in energy units selected in joint)
P. Novák
Inst. of Physics, [Link], Prague, Czeck Republic
email: novakp@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary,
we will communicate the problem to the authors.
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
~ = µ0 µ [3(~nr̂)r̂ − ~n]
B (8.1)
4πr3
where r̂ = ~r/r.
~n = M ~ /M is direction of magnetization.
µ0 is permeability of free space; µ0 = 4π10−7 H/m.
~ is the dipolar field in T.
B
~ is magnetic dipolar moment in Am2 = J/T, assumed to be parallel to ~n.
µ
r is in m.
We want to express µ in Bohr magnetons µB =9.274078.10−24 J/T and
r in atomic units for length a0 (Bohr radius) a0 =5.2917706.10−11 m.
Inserting in (1) gives
1 X~
Ean = − Bj µ
~j (8.4)
2V j
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 Bj in T, µj in µB and V in (a.u.)3 gives
3.129232.107 X ~
Ean (J/m3 ) = − Bj (T)~
µj (µB ) (8.5)
V (a.u.)3 j
156 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
8.19.1 Execution
8.19.3 Input
Rmax max distance (bohr) for lattice summation. Vary it for convergence
check.
IPR Print switch. IPR=2 produces very large files [Link] and
[Link] dipan
line 2 must be repeated for every non-equivalent atom in the unit cell
line 3: free format
I 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.
I 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 [Link].
I Run lapw1 with this k-mesh.
I Run spaghetti with input-options such that it prints the bands which intersect EF to
[Link] ene (line 10, see sec. 8.3)
I Edit [Link] 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
[Link] 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
I Run spagh2rho < [Link] ene to convert from this format into a format
which is compatible with the [Link] 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:
I 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, ... [Link] generated in the
previous steps. It requires the public domain package pgplot and the contour-plot
program plotgenc. (The latter can be obtained from
[Link] user/unsupported/, but you must have installed the
pgplot library before.)
158 CHAPTER 8. ANALYSIS, PROPERTIES AND OPTIMIZATION
9 Utility Programs
Contents
9.1 symmetso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
9.2 pairhess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
9.3 eigenhess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
9.4 patchsymm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
9.5 afminput . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
9.6 clmcopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
9.7 reformat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
9.8 hex2rhomb and rhomb in5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
9.9 plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
9.10 add columns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
9.11 clminter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
9.12 eosfit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
9.13 eosfit6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
9.14 spacegroup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
9.15 join vectorfiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
9.16 arrows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
9.17 xyz2struct . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
9.18 cif2struct . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
9.19 struct2cif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
9.20 StructGen of w2web . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
9.21 supercell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
9.22 structeditor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
9.23 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172
9.24 Unsupported software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
9.1 symmetso
This program helps to setup spin-orbit calculations in magnetic systems. Since SO may break
symmetry in certain spacegroups, it classifies your symmetry operations into operations A, which
do not invert the magnetization (identity, inversion, rotations with the rotation axis parallel to
magnetization), B, which invert it (mirror planes) and C, which change the magnetization in some
other way. (Note: magnetization is a result of a circular current, or equivalently, an axial vector
resulting from a vector product ẑ ∼ x̂ × ŷ). symmetso will keep all A-type and throw away all
C-type symmetry operations. Depending on the presence of inversion symmetry it will keep
(inversion is present) or remove the B-type operations. Finally, symmetso uses the remaining
159
160 CHAPTER 9. UTILITY PROGRAMS
symmetry operations to check/generate equivalent atomic positions (it can happen that some
equivalent atoms become non-equivalent after inclusion of SO interaction).
In essence, it reads your [Link] and [Link] (for the direction of magnetization) files
and creates an ordered [Link] orb file with proper symmetry and equivalent atoms. It
also generates a file [Link], which is a struct file with valid operations to generate a proper
k-mesh using ’’x kgen -so’’. In addition proper input files case.in1, case.in2,
[Link], [Link]/dn, [Link]/dn, [Link], [Link]/dn are
generated, so that you can continue with runsp -so without any further changes.
9.1.1 Execution
The program symmetso is executed by invoking the command:
Usually it is called from the script initso lapw and thus needs not to be invoked manually.
9.2 pairhess
This program was contributed by:
This program creates an approximate hessian matrix (in .minpair) for structure minimization
using the PORT option. It uses a harmonic model with exponentially decaying bond strenght and
in many cases reduces the number of geometry steps during min lapw significantly. It is
described in detail in Rondinelli et al. 2006.
For its usage see the comments in sect. 5.3.2.
9.2.1 Execution
The program pairhess is executed by invoking the command:
The switch -copy copies .minpair to .minrestart and .min hess, which are needed in
min lapw.
9.2.3 Input
pairhess uses an optional input file [Link], which is needed only for an experienced
user for better tailoring of certain default parameters.
An example is given below:
RESCALE A scaling term to multiply the pairwise hessian by. This number is
rather important; 0.25 appears to be best for a system with soft modes,
0.35 for a stiffer system. You can save substantial time by adjusting
RESCALE so it is approximately correct using a .min hess from a previ-
ous run (adjust until numbers for similar multiplicities are similar), or
by adjusting the frequencies (see also eigenhess).
CUTOFF When the weighting (via an exponential decay) becomes smaller than
this number the pairwise bonds are ignored.
DIAG The value to multiply a unitary matrix by, this is added to the hessian
estimate
9.3 eigenhess
Laurence Marks
Dept. Materials Science and Engineering
Northwestern University
Evanston, USA
l-marks@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
This program analyses / manipulates .min hess, which was created by a structural
minimization using min lapw and the “PORT” option. In particular, such an analysis can yield
approximate vibrational frequencies and corresponding eigenmodes, which eventually can give a
hint about a dynamically unstable structure (imaginary frequencies). Some more description is
given in $WIENROOT/SRC pairhess/README.
The program eigenhess is executed by invoking the command:
x eigenhess
9.4 patchsymm
This program performs a symmetry check on the positions and produces a new struct file
[Link] new. It is useful in case something went wrong during min lapw (rounding
errors of positions) or the cif/amc file did not have enough digits (eg. “1/3” was prepresented by
“0.33333” only). The file [Link] gives information on how parameters changed.
9.4.1 Execution
9.5 afminput
This program creates the inputfile [Link] st for the program clmcopy, which copies
spin-up densities of atom i to spin-down densities of the related antiferromagnetic atom j and vice
versa in an anti-ferromagnetic system. It uses a symmetry operation to find out how and which
atomic densities must be interchanged and how the Fourier coefficients of the density transform.
It is based on the ideas of Manuel Perez-Mato (Bilbao, Spain).
See $WIENROOT/SRC afminput/afminput test for several examples.
The best way is to supply a file [Link] supergroup, which is the struct file of the
nonmagnetic supergroup. If the two spacegroups are “TRANSLATIONENGLEICH”, it will find
out automatically the proper symmetry operation. Please note, this automatic way works only when
the coordinate system remains identical. In some cases sgroup may interchange eg. the y and z axis. In such
cases reverse this change, both, for the lattice parameters as well as for all positions, set NSYM=0 and run
init lapw again (ignoring any suggestion of sgroup).
If the two spacegroups are “KLASSENGLEICH” (i.e. have the same number of symmetry
operations), you will be asked to supply a translation which transforms the AF atoms into each
other. A typical example would be bcc Cr: the bcc supergroup and the AF subgroup (simple
cubic) have both 48 symmetry operations and the proper translation is (0.5,0.5,0.5).
Finally, if you don’t give [Link] supergroup, you have to supply a symmetry operation
(rotation + non-primitive translation) as input. For bcc Cr or the famous NiO-AFII structure this
would be simply
1.0 0.0 0.0 0.5
0.0 1.0 0.0 0.5
0.0 0.0 1.0 0.5
Please see the comments in sect. 4.5.4 on how to proceed in detail for AFM calculations and find
further examples in SRC afminput.
9.5.1 Execution
9.6 clmcopy
This program generates the spin-dn density ([Link]) from a given spin-up density
([Link]) according to rules and symmetry operations in [Link] (generated
earlier by afminput) for an AFM calculation.
Please see the comments in sect. 4.5.4 on how to proceed in detail for AFM calculations.
164 CHAPTER 9. UTILITY PROGRAMS
9.6.1 Execution
9.6.3 Input
NATOM Number of atoms for which rules for copying the density will be de-
fined
NLM Number of LM values, for which you have to change the sign when
swapping up and dn-densities
L1,M1,L2,M2,Fac NLM pairs of L1,M1 (spin-up), which change into L2,M2 (spin-dn) and
the respecting CLMs are multiplied by Fac
9.7 reformat
To produce a surface plot of the electron density using rhoplot lapw (which is an interface to
gnuplot), data from the file [Link] created by lapw5 must be converted using reformat
The sources of the program reformat.c are supplied in SRC reformat.
9.9 plane
plane helps to generate case.in5 for density plots with lapw5 (for orthogonal and hex lattices
only). The plane will be specified by 3 atoms and you need an auxiliary file [Link], which
contains:
9.11 clminter
clminter interpolates the density in [Link]/up/dn to a new radial mesh as defined in
[Link] new. This utility is usefull when you run a structural minimization (min lapw),
some atoms start to overlap and you have to reduce RMT (the size of the atomic spheres) of
certain atoms. In such a case:
Note: Please be aware the the total energy will change with modified RMT (by some constant)
and you must not compare energies comming from different RMTs (but most likely you can
determine the constant shift by repeating (at least) ONE calculation with identical structure but
different RMTs).
The source of this program is supplied in SRC trig.
9.12 eosfit
Small program to calculate the Equation of States (EOS; Equilibrium volume V0 , Bulk modulus B0
and it’s derivative B00 . The Murnaghan (1944), the Birch-Murnaghan and the EOS2 equation of
states are supported. It relies on the file [Link] (containing lines with ”volume, E-tot”, usually
created from w2web using ”Volume optimization”), or alternatively is called from eplot lapw
using [Link] (see 5.7.1 and 5.3.1).
The sources are supplied in SRC eosfit.
9.13 eosfit6
Nonlinear least squares fit (using PORT routines) for a parabolic fit of the energy vs. 2-4 dim.
lattice parameters. It requires [Link] and [Link], usually generated by
parabolfit lapw. It can optionally produce [Link], which contains energies on a
9.14. SPACEGROUP 167
specified grid for plotting purposes (in 2D same format as [Link], which can be used in
contourplot programs). (See 5.3.1).
The sources are supplied in SRC eosfit6.
9.14 spacegroup
This program was contributed by:
Vaclav Petricek
Institute of Physics
Academy of Sciences of the Czech Republic
Na Slovance 2
182 21 Praha (Prague) 8
Czech Republic
petricek@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
Interactive program to generate equivalent positions for a given spacegroup and lattice. The
program is also used internally from w2web to generate positions when selecting spacegroups in
the StructGen.
Phillipp Wissgott
Institute of Solid State Physics
TU Vienna
wissgott@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
Interactive program to combine parallel vector and energy files ([Link] xx and
[Link] xx) into single files ([Link] and [Link]).
Executed by:
9.16 arrows
This program was contributed by:
168 CHAPTER 9. UTILITY PROGRAMS
Evgeniya Kabliman
Institute for MaterialsChemistry
TU Vienna
evgeniya@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
Small program which together with Xcrysden allows to display the “forces acting on all atoms”
or the “differences between two structures” using arrows which indicate the movement of the
atoms. The recommended sequence to visualize forces is:
I Prepare (copy) a struct and scf file with the initial structure using the names
case [Link] and case [Link].
I View case [Link] in Xcrysden and “File/Save as xsf-structure” with the
name case [Link].
I x arrows
I View the resulting case [Link] using: xcrysden --xsf case [Link].
Switch on “Display/Forces” and adjust the length of the arrows in “Modify/Force-settings”.
while differences between the inital and relaxed structure can be viewed by:
I Prepare (copy) two struct files with the initial and the relaxed structure using the names:
case [Link] and case [Link].
I View case [Link] in Xcrysden and “File/Save as xsf-structure” with the
name case [Link].
I x -delta arrows
I View the resulting case [Link] using: xcrysden --xsf case [Link]. Switch
on “Display/Forces” and adjust the length of the arrows in “Modify/Force-settings”.
9.17 xyz2struct
ang
7.47700 7.47700 7.47700
B 4.98466667 1.24616667 0.00000000
C 6.23083333 2.49233333 0.00000000
.....
or
BOHR F L
17.47700
0.470724637 0.492808141 0.000000000
-0.471118220 0.493012774 0.000000000
0.000000000 0.000000000 0.680559876
B 0.00000000 0.00000000 0.00000000
C 0.14300000 0.14300000 0.25000000
.....
9.18 cif2struct
cif2struct reads structural data in cif-format from [Link] and writes them into
[Link]. It is executed using:
cif2struct [Link]
The required cif files can be for example be obtained from Cystallographic databases (e.g. the
Inorganic Crystal Structure DataBase ICSD) or from other programs.
Alternatively, cif2struct can work with [Link], which contains the following data:
a # a..Ang, b..Bohr
0.0 0.0 0.0 # shift of origin
4.7554 4.7554 12.991 90. 90. 120. # a,b,c,angles
’R-3c’ # spacegroup-symbol (see \STRUCTGEN{})
’Al’ # atom-name
0.0000000 0.0000000 0.3520000 # atomic position
’O’ # ...
0.3063000 0.0000000 0.2500000 # ...
...
9.19 struct2cif
The new StructGen helps to generate the master input file [Link]. It has the following
additional features:
I built in calculator to carry out simple arithmetic operations to specify the position pameters
(of the equivalent atoms). Each position of equivalent atoms can be entered as a number, a
fraction (e.g. 1/3) or a simple expression (e.g. 0.21 + 1/3). The first position defines the
variables x, y and z, which can be using in expression defining the other positions (e.g. −y,
x, −z + 1/2).
9.21 supercell
9.21.1 Execution
supercell or x supercell
9.22 structeditor
Robert Laskowski
email: rolask@[Link]
Please make comments or report problems with this program to the WIEN-mailinglist. If necessary, we will
communicate the problem to the authors.
9.22. STRUCTEDITOR 171
This package helps to manipulate structures. Usually one would start from an appropriate
(simple) [Link] file, and this tool allows to add or manipulate atoms (with or without
symmetry considerations), or generate arbitrary supercells or surfaces. It is commandline driven
and targeted for the more experienced user, who “knows what he wants to do” and is just looking
for a convenient tool.
It consists of a couple of octave (mathlab) routines and some fortran code, thus it requires octave
(the free mathlab version) and for visualization the opendx package ([Link]
and [Link] ).
A full documentation and some examples can be found in
$WIENROOT/SRC structeditor/doc, but the main commands are:
9.22.1 Execution
octave
s=loadstruct(”[Link]”)
# make an orthorhombic supercell and visualize it
a=[1 0 0; 1 1 0; 0 0 2]
sout=makesupercell (s,a);
showstruct(sout);
# save it as [Link]
savestruct (sout,”[Link]”);
# get help on all commands
helpstruct
172 CHAPTER 9. UTILITY PROGRAMS
9.23 Visualization
9.23.1 BALSAC
balsac (Build and Analyze Lattices, Surfaces and Clusters) was written by Klaus Hermann
(Fritz-Haber Institut, Berlin). It provides high quality postscript files. In SRC balsac-utils we
provide the following interface programs to convert from WIEN2k to balsac:
For an example see figure 3.1 For scientific questions concerning BALSAC please contact Klaus
Hermann at hermann@[Link]
Balsac is available from:
9.23.2 XCrysDen
XCrysDen (Kokalj 1999) is a render and analysis package. It has the following features (see also
[Link]
Tone Kokalj
Jozef Stefan Institute, Dept. of Physical and Organic Chemistry
Jamova 39, SI-1000 Ljubljana, Slovenia
Tel.: +386 61 177 3520, Fax: +386 61 177 3811
[Link]@[Link]
[Link]
9.24. UNSUPPORTED SOFTWARE 173
Three test cases are provided in the WIEN2k package. They contain the two starting files
[Link] and [Link] and all the output so that you can compare your results with
them.
The test cases are the following (where the names correspond to what was called CASE in the rest
of this User’s Guide)
TiC
Fccni
TiO2
We recommend to run these test cases (in a different directory) and compare the output to the
provided one. All test cases are setup such that the CPU-time remains small (seconds). For real
production runs the value of RKMAX in case.in1 must be increased and a better (denser)
k-mesh should be used.
In addition we provide a subdirectory example struct files were various more complicated
struct files can be found.
10.1 TiC
Ferromagnetic Nickel is a test case for a spin-polarized calculation. Ni has the atomic
configuration 1s2 , 2s2 , 2p6 , 3s2 , 3p6 , 3d8 , 4s2 or [Ar] 3d8 , 4s2 . We treat the 1s, 2s, 2p and 3s as core
states, and 3p (as local orbital), 3d, 4s and 4p are handled as valence states. In a spin-polarized
calculation the file structure and the sequence of programs is different from the
non-spin-polarized case (see 4.5.2).
Create a new session and its corresponding directory. Generate the structure with the following
data (we can use a large sphere as you will see from the output of nn):
175
176 CHAPTER 10. EXAMPLES
Title fcc Ni
Lattice F
a 6.7 bohr
b 6.7 bohr
c 6.7 bohr
α, β, γ 90
Atom Ni, enter position (0,0,0) and RMT = 2.3
Initialize the calculation using the default RKmax and use 3000 k-points (a ferromagnetic metal
needs many k-points to yield reasonably converged magnetic moments). Allow for
spin-polarization.
Start the scf cycle (runsp lapw) with ”-cc 0.0001” (in particular for magnetic systems charge
convergence is often the best choice). At the bottom of the converged scf-file ([Link]) you
find the magnetic moments in the interstital region, inside the sphere and the total moment per
cell (only the latter is an “observable”, the others depend on the sphere size).
Title TiO2
Spacegroup P 42 /mnm (136)
a 8.682 bohr
b 8.682 bohr
c 5.592 bohr
α, β, γ 90
Atom Ti, enter position (0,0,0) and RMT = 2.0
Atom O, enter position (0.3,0.3,0) and RMT = 1.6
cp $WIENROOT/SRC_templates/.machines .
and edit this file. You would start the scf-cycle (in background) simply by typing
During the scf-cycle monitor [Link] and check convergence (:ENE, :DIS, :FGL002),
either using “Utils/Analysis” in w2web, or ‘‘grep :ENE [Link]’’. You should see some
convergence of :FGL002 and then a big jump in the final cycle, when the valence-force corrections
are added. Only the last force (including this correction) is valid.
Since this force is quite large, you can now optimize the position of the O-atom:
Start the structure minimization in w2web using “Execution o [Link]”. This will generate
[Link], and you can try option PORT with tolf=1.0 (instead of 2.0), otherwise stay with the
default parameters. Repeat “Execution o [Link]” and start the minimization.
Alternatively you can use
min_lapw -p
This will create [Link] automatically, call the program min, which generates a new struct file
using the calculated forces, and continues with the next scf cycle. It will continue until the forces
are below 1 mRy/bohr ([Link]) and the final results are not “saved” automatically but can be
found in the “current” calculation.
You should watch the minimization (:ENE, :FGL002, :POS002) using the file [Link] mini,
which contains the final iteration of each geometry step (see also Sec.5.3.2). If the forces in this file
oscillate from plus to minus and seem to diverge, or if they change very little, you can edit
[Link] (change the method, reduce or increase the stepsize), and remove [Link]
(contains the “history” of the minimization and is used to calculate the velocities of the moving
atoms). (This should not be neceaasry for the rutile example, but may occur in more complex
minimizations. See comments in Sec. 5.3.2).
The final structural parameter of the O-atom should be close to x=0.304, which compares well
with the experimental x=0.305.
mkdir super
cd super
cp ../TiC/[Link] .
x supercell
Specify “[Link]”, a “2x2x2” supercell, “F” lattice (this will create a cell with 16 atoms, you can
also create 32 or 64 atom cells using B or P lattice type. Note: surfaces require a P supercell).
cp TiC_super.struct [Link]
and edit this file to make some changes. You could eg.
178 CHAPTER 10. EXAMPLES
Note: it is important to make at least one of these chages. Otherwise the initialization will restore the
original unit cell (or the calculations will fail later on because symmetry is most likely not correct)
Run init lapw. You will see that nn complains and finds a new set of equivalent atoms
(originally all atoms were non-equivalent, but nn finds that some atoms have identical neighbors,
thus should be in an equivalent set). Accept the automatically generated struct file and continue.
Remember, supercells normally require less k-points than the original small cell.
After the complete initialization you may in principle restore the original struct file (eg. without a
displacement) in case you want to “repeat” the undistorted structure in supercell geometry.
For a “core-hole” calculation you would now edit [Link] and remove one core electron from
the desired atom and state (1s or 2p, ...). In addition you should add the missing electron either in
[Link] (background charge) or super.in2 (add it to the valence electrons). In the latter
case, you should remove this extra electron AFTER scf and BEFORE calculation of the spectra.
Once this has been done, you could start a scf-cycle (for impurities, vacancies,.. you should most
likely also optimize the internal positions).
179
11 Installation and Dimensioning
Contents
11.1 Requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
11.2 Installation of WIEN2k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
11.3 w2web . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
11.4 Environment Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
11.1 Requirements
WIEN2k is written in FORTRAN 90 and requires a UNIX operating system since the programs are
linked together via C-shell scripts. It has been implemented successfully on the following
computer systems: Intel and AMD based PCs running under Linux, IBM RS6000, HP, SGI, DEC
Alpha and SUN. Hardware requirements will change from case to case (small cases with 10 atoms
per unit cell can be run on a Pentium-II PC with 128 MB under Linux), but we recommend a more
powerful PC or workstation with at least 256 MB (better 512 MM or more) memory and plenty of
disk space (a few Gb).
For coarse grain parallization on the k-point level, a cluster of PCs with a 100 Mb/s (better 1Bg/s)
network is sufficient. Faster communication (Infiniband) is recommended for the fine grain
(single k-point) parallel version.
For Intel (AMD) based systems we recommend the Intel ifort compiler and the Intel
mkl library (which includes blas, lapack and Scalapack) (see [Link] If you
have installed ifort yourself on your local PC, don’t forget to configure your environment
properly. Add some thing like:
181
182 CHAPTER 11. INSTALLATION AND DIMENSIONING
I pdf-reader, acroread
I Tcl/Tk-Toolkit (for Xcrysden only)
I MPI+SCALAPACK (for fine grain parallelization only)
I FFTW-2.1.5 (mpi-version for fine grain parallelization only)
Usually these packages should be available on modern systems. If one of these packages is not
available it can either be installed from public domain sources (ask your computing center, use the
WWW to search for the nearest location of these packages) or the corresponding configuration
may be changed (e.g. using vi instead of emacs). Brief installation instructions for mpich and fftw
are given below. None of the principal components of WIEN2k requires these packages, only
w2web needs them.
Optionally, one can also use in the sequential (non-mpi) version of lapw0 and lapw2 the fftw
routines. However, there is some speedup only when you use the MKL-fft routines, not the
self-compiled fftw-binaries. The mkl-interface to fftw is not active by default, but you may have to
compile it yourself. To do so (syntax ifort12):
I cd $MKLROOT/interfaces/fftw2xf
I make libintel64 (or make libia32)
I in case you do not have icc installed, but use GNU-C (gcc) you must:
– edit makefile, and remove -D GNU from the line “ CC=gcc -D GNU” (to remove
additional from the object names)
– make libintel64 compiler=gnu
I add -DFFTW2 to FOPT in the Makefile of lapw0 and lapw2
I add -lfftw2xf or -lfftw2xf gnu to R LIBS in the Makefile of lapw0 and lapw2
(e.g. ./WIEN2k). In addition you can download three examples, namely [Link],
[Link] and [Link].
Uncompress and expand all files using:
tar -xvf wien2k [Link] (skip this if you downloaded files separately)
gunzip *.gz
chmod +x ./expand lapw
./expand lapw
./SRC
SRC_2Doptimize
SRC_afminput
SRC_aim
SRC_arrows
SRC_balsac-utils
SRC_broadening
SRC_cif2struct
SRC_clmaddsub
SRC_clmcopy
SRC_dipan
SRC_dstart
SRC_elast
SRC_eosfit
SRC_eosfit6
SRC_filtvec
SRC_fsgen
SRC_initxspec
SRC_irrep
SRC_joint
SRC_kgen
SRC_kram
SRC_lapw0
SRC_lapw1
SRC_lapw2
SRC_lapw3
SRC_lapw5
SRC_lapw7
SRC_lapwdm
SRC_lapwso
SRC_lcore
SRC_lib
SRC_lorentz
SRC_lstart
SRC_mini
SRC_mixer
SRC_nn
SRC_optic
SRC_optimize
SRC_orb
SRC_pairhess
SRC_phonon
SRC_qtl
SRC_reformat
SRC_sgroup
SRC_spacegroup
SRC_spaghetti
SRC_structeditor
SRC_sumpara
SRC_supercell
SRC_symmetry
SRC_symmetso
SRC_telnes3
SRC_templates
SRC_tetra
SRC_trig
SRC_txspec
SRC_usersguide_html
SRC_vecpratt
example_struct_files
TiC
TiO2
fccni
184 CHAPTER 11. INSTALLATION AND DIMENSIONING
Thus, each program has its source code (split into several files) in its own subdirectory. All
programs are written in FORTRAN90 (except SRC sgroup and SRC reformat, which are in C).
/SRC contains the users guide (in form of a postscript file [Link] and as pdf-file
[Link]), all c-shell scripts and some auxiliary files.
/SRC usersguide html contains the html version of the UG.
/Fccni, /TiC and /TiO2 contain three example inputs and the respective outputs.
/example struct files contains a collection of various struct files, which could be of use
especially for the less experienced user.
/SRC templates contains various input templates.
In addition to the expansion of the tar-files ./expand lapw copies also all csh-shell scripts from
/SRC to the current directory and creates links for some abbreviated commands.
At the end of expand lapw you will be prompted to start the script
./siteconfig lapw
When you start this script for the first time (file INSTALLDATE not present), you will be guided
through the setup process.
Later on you can use siteconfig lapw to redimension parameters, update individual packages
and recompile the respective programs.
During the first run, you will be asked to specify:
I your system; at this point system specific files (e.g. cputim.f will be installed. If your
system is not listed, use the system generic, which should compile on any machine.
I your FORTRAN90 and C compilers;
I your compiler and linker options as well as the place for LAPACK and BLAS libraries.
Depending on the system you selected, we have included some recommended compiler and
linker options, which are known to work on our systems (use generic when you have
problems here; see also sec. 11.2.4). On some systems it is required to specify LAPACK and
BLAS libraries twice (i.e. R LIBS:-llapack lapw -lblas lapw -llapack lapw -lblas lapw). This
generates Makefiles from the corresponding [Link] in all subdirectories.
I configuration of parallel execution will ask whether your system is shared memory, so that
default parameters can be set accordingly ( $WIENROOT/parallel options is the file
where this information is stored).
I to configure parallel execution for distributed systems, specify the command to open a
remote shell, which on most systems is rsh or ssh.
I You will then be asked wether you want to run fine-grained parallel. This is only possible if
FFTW, MPI and SCALAPACK (included in newer mkl-versions) are installed on your
system and requires a fast network (100Mb/s is not enough) or a shared memory machine.
It pays off only for bigger cases (matrixsize > 5000).
I You should define NMATMAX, i.e. the maximum matrixsize (number of basisfunctions).
This value should be adjusted according to the memory of your hardware. Rough estimates
are:
NMATMAX= 5000 ==> 256MB (real, i.e. with inversion symmetry)
NMATMAX=10000 ==> 1GB (real) (==> cells with about 80-150 atoms/unitcell)
If you choose it too large, lapw1 will start to “page” leading
√ to inacceptable performance or
a crash. NMATMAX will √ be automatically recuced (by 2) for complex (without inversion)
cases and increased by N P E for mpi-parallel cases.
11.2. INSTALLATION OF WIEN2K 185
I Now you are prompted to compile all programs (this will be done using make) and the
executables are copied to the $WIENROOT directory. Compilation might take quite some
time.
I During compilation watch for error messages on the screen. If there are errors, you may
need to change into the corresponding SRC * directory and examine file [Link] for
details.
Common errors are wrong specification of compiler, linking or library options. In such
cases, adopt the Makefile in this directory and recompile using make. Once you have proper
options, correct them globally in siteconfig lapw and recompile.
Later on you can use siteconfig lapw to change parameters, options or to update a package.
to your ˜/.cshrc or˜/.bashrc file. Eventually you should also edit these files and set the
$LD LIBRARY PATH variable (path where compiler-libs or blas-libraries are located).
Note: This will work only when the csh, tcsh or bash-shell is your login shell. Depending on your settings
you may have to add similar lines also in your .login file. If you are using a different login-shell, edit
your startup files manually.
I Find out which compiler options (man f90) make these programs run faster. You could
specify a higher optimization (-O3), or specify a particular processor architecture
(-qarch=pwr5 or -R10000, ....).
I Good performance depends on highly optimized BLAS (and much less on LAPACK)
libraries. Whenever it is possible, replace the supplied libraries (SRC lib/blas lapw
SRC lib/lapack lapw), by routines from your vendor (mkl for Intel or AMD processors,
aclm for AMD, essl for IBM, sunperf for SUN, [Link] on SGI, ...). If such libraries
are not available use the GOTO-library
([Link] which is a very competitive
alternative. (Eventually you may try to optimize them yourself using the “ATLAS” system
(see [Link] but this is no longer recommended. We provide an
“old” ATLAS-BLAS for a Pentium3 with WIEN2k.
WIEN2k is written in Fortran 90 and all important arrays are allocated dynamically. The only
important parameters left are NMATMAX and NUME, specifying the maximum matrixsize
(should be adjusted to the memory of your hardware, see above) and the maximum number of
eigenvalues (must be increased for unitcells with large number of electrons)
Some less important parameters are still present and described in chapter “dimensioning
parameters“ of the respective section in chapter 6.
We recommend to use siteconfig lapw for redimensioning and recompilation. In order to
work properly, the parameter XXXX in the respective [Link] files must obey the following
syntax:
PARAMETER(XXXX= ....)
w2web requires perl, which should be available on most systems. (If not contact your system
administrator or install it yourself from the WWW)
When you start w2web for the first time on the computer where you want to execute WIEN2k
(you may have to telnet, ssh,.. to this machine) with the command w2web [-p xxxx], you will
be asked for a username/password (I recommend you use the same as for your UNIX login).
You must also specify a “port” number (which can be changed the next time you start w2web). If
the default port (7890) used to serve the interface is already in use by some other process, you will
get the error message w2web failed to bind port 7890 - port already in use!.
Then you will have to choose a different port number (between 1024 and 65536) . Please
remember this port number, you need it when connecting to the w2web server.
Note: Only user root can specify port numbers below 1024!
Once w2webhas been started, use your favorite WWW-browser to connect to w2web, specifying
the correct portnumber, e.g. firefox [Link] where w2web runs:7890
On certain sites a firewall may block all high ports and one cannot connect to this machine. In
these cases you can create a ssh-tunnel using the following commands:
11.3. W2WEB 187
At your “local host” (the PC in front of you) connect to the “w2web host” (where you started
w2web) using
On your local host use a web browser and connect with: firefox http:[Link]:2000.
Using “Configuration ” you can further tailor the behaviour according to your wishes. In particular
you can define new “execution types” to adjust to your queuing system.
For example the line
batch=batch < %f
defines an execution type “batch” using the UNIX batch command. (w2web collects its
commands in a temporary script and you can access it using %f).
If you run on a machine with a queuing system (like loadleveler, sun-grid-engine, or pbs) you
may define an “execution type”
The following scripts may serve as templates: qsub-wienjob lapw in $WIENROOT needs a
master-job-template qsub-job0 lapw and examples for loadleveler and SGE are provided in
$WIENROOT (you may need to adapt them ! Other examples you can find on our FAQ-page on the
web). Of course, with some small modifications you can define several “execution types” with eg.
different number of processors or mpi vs. k-point parallel runs,....
w2web saves several variables in startup files which are in the (˜/.w2web) directory.
w2web acts like a normal web-server - except that it runs on a ”user level port” instead of the
default http-port 80. It serves html-files and executes perl-scripts or executes system or user
commands on the server host.
w2web creates on the first start of w2web on host “hostname” the directory .w2web/hostname
in your home directory with the following content:
I .w2web/hostname/conf
I .w2web/hostname/logs
I .w2web/hostname/sessions
In this file various configuration parameters are stored by w2web. To restrict the access to certain
IP addresses you can add lines like:
deny=*.*.*.*
allow=128.130.134.* [Link]
188 CHAPTER 11. INSTALLATION AND DIMENSIONING
In order to use the https-protocol the perl-library Net::SSLeay in addition to the OpenSSL package
must be installed on your system. Both are freely available.
Then you must include a line with ssl=1 in [Link].
If you run w2web-server in ssl-mode you need a site certificate for your server. You may use the
supplied certificate in $WIENROOT/SRC w2web/bin/[Link] (copy this file to your
conf-directory and set the keyfile= /.w2web/¡hostname¿/conf/[Link] line in your
[Link]).
This certificate will not expire until 2015, but usually browsers will complain that they do not
know the Certificate Authority who issued this certificate - if you don’t like this message, you
must buy a certificate from VeriSign, Thawte or a similar CA.
Of course you must connect to https: instead of http:, i.e. use:
netscape [Link] where w2web runs:7890.
In this chapter hints are given for solving some difficulties that have occurred frequently. This
chapter is by no means complete and the authors would appreciate further suggestions which
might be useful for other users. Beside the printed version of the users guide, an online pdf
version is available using help lapw. You can search for a specific keyword (use ∧ f keyword)
and hopefully find some information.
If an error occurs in one of the SCF programs, a file [Link] is created and an error message
is printed into these files. The run lapw script checks for these files and will automatically stop if
a non-empty error file occurs.
Check the files [Link] (which is written by init lapw and run lapw), :log (where a
history of all commands using x is given) and *.error for possible explanations.
In several places the dimensions are checked. The programs print a descriptive error message and
stop.
[Link]: This file gives error messages if the atomic spheres overlap. But it should also be
used to check the distances between the atoms and the coordination number (same
distance). If inconsistencies exists, your [Link] file may contain an error. A check
for overlapping spheres is also included in mixer and lapw1.
[Link]: Possible stops or warnings are:
“NO SYMMETRY OPERATION FOUND IN ROTDEF“: This indicates that in your
[Link] file either the positions of equivalent atoms are not specified correctly (only
positive coordinates allowed!!) or the symmetry operations are wrong.
case.output1: Possible stops or warnings are:
“NO ENERGY LIMITS FOUND IN SELECT“: This indicates that Etop or Ebottom could
not be found for some ul (r, El ). Check your input if it happens in the zeroth iteration.
Later, (usually in the second to sixth iteration) it may indicate that in your SCF cycle
something went wrong and you are using a crazy potential. Usually it means that
mixing of the charge densities was diverging and large charge fluctuations occured.
Check previous charges for being physically reasonable (grep for labels :NTOxx
:CTOxx :DIS :NEC01). Usually this happens when your input is not ok, or for very ill
conditioned problems (very rare), or more likely, when “Ghostbands” appeared (or
189
190 CHAPTER 12. TROUBLE SHOOTING
some states were missing) because of bad energy parameters in case.in1. You will
probably have to delete [Link]* and [Link], rerun x dstart and then
change some calculational parameters. These could be: fixing some energy parameter
(modify both, case.in1 and case.in1 orig or try the -in1orig switch if you have
used -in1new); switch to a broadening method (TEMP with eg. 0.010 mRy); or increase
the k-mesh (magnetic metals); or reduce the mixing parameter in [Link] slightly
(eg. to 0.1). In very difficult (magnetic) cases a PRATT mixing with eg. 0.01 mixing
might be helpful at the beginning of the scf cycle (but later switch to MSEC1 again) !
“STOP RDC 22“: This indicates that the overlap matrix is not positive definite. This
usually happens if your [Link] file has some error in the structure or if you
have an (almost) linear dependent basis, which can happen for large RKMAX values
and/or if you are using very different (extremely small and large) sphere radii RM T .
“X EIGENVALUES BELOW THE ENERGY emin“: This indicates that X eigenvalues were
found below emin. Emin is set in case.in1 (see sec. 7.3.3) or in [Link]
generated by KGEN, see 6.3, 6.5). It may indicate that your value of emin is too high or
the possibility of ghostbands, but it can also be intentional to truncate some of the low
lying eigenvalues.
If you don’t find enough eigenvalues (e.g.: in a cell with 4 oxygens you expect 4 oxygen s
bands at roughly -1 Ry) check the energy window (given at the end of the first k-point
in case.in1 or in [Link]) and make sure your energy parameters are found by
subroutine SELECT or set them by hand at a reasonable value.
case.output2: Possible stops or warnings are:
“CANNOT BE FOUND“: This warning, which could produce a very long output file,
indicates that some reciprocal K-vector would be requested (through the k-vector list of
lapw1), but was not present in the list of the K generated in lapw2. You may have to
increase the NWAV, and/or KMAXx parameters in lapw2 or increase GMAX in
case.in2. The problems could also arise from wrong symmetry operations or a
wrong structure in [Link].
“QTL-B VALUE“: If larger than a few percent, this indicates the appearance of ghost bands,
which are discussed below in section 12.1.
The few percent message (e.g up to 10 %) does not indicate a ghost band, but can
happen e.g. in narrow d-bands, where the linearization reaches its limits. In these cases
one can add a local orbital to improve the flexibility of the basis set. (Put one energy
near the top and the other near the bottom of the valence band, see section 7.3.3).
FERMI LEVEL not converged (or similar messages). This can have several reasons: i) Try a
different Fermi-Method (change TETRA to GAUSS or TEMP in case.in2). ii) Count
the number of eigenvalues in case.output1 and compare it with the number of
valence electrons. If there are too few eigenvalues, either increase EMAX in
[Link] (from 1.5 to e.g. 2.5) or check if your scf cycle had large charge
oszillations (see SELECT error above)
If the SCF cycle stops somewhere (especially in the first few iterations), it is quite possible, that the
source of the error is actually in a previous part of the cycle or even in a previous (e.g. the zeroth)
iteration. Check in the [Link] file previous charges, eigenvalues, . . . whether they are
physically reasonable (see SELECT error above).
The first case may occur in a system which has atoms with very different atomic sphere radii.
Suppose you calculate a hydroxide with very short O-H bonds so that you select small RM T radii
for O and H such as e.g. 1.0 and 0.6 a.u., respectively. The cation may be large and thus you could
choose a large RM T of e.g. 2.4 a.u. However, this gives a four time larger effective RKmax for the
cation than for H, (e.g. 16.0 when you select RKmax=4.0 in case.in1). This enormous difference
in the convergence may lead to unphysical eigenvalues. In such cases choose lmax=12 in
case.in1 (in order to get a very good re-expansion of the plane waves) and reduce RM T for the
cation to e.g. 1.8 a.u.
The second case can occur when you don’t use a proper set of local orbitals. In this situation the
energy region of interest (valence bands) falls about midway between two states with different
principle quantum numbers, but with the same l-value (for one atom).
Take for example Ti with its 3p states being occupied as (semi-core) states, while the 4p states
remain mostly unoccupied. In the valence band region neither of those two states (Ti 3p, 4p)
should appear. If one uses 0.2 Ry for the expansion energy E(1) for the p states of Ti, then Ti-p
states do appear as ghost bands. Such a run is shown below for T iO2 (rutile).
The lowest six eigenvalues at GAMMA fall between about -1.30 and -1.28 Ry. They are ghost
bands derived from fictitious Ti-p states. The next four eigenvalues between -0.94 and -0.78 Ry
correspond to states derived from O 2s states, which are ok, since there are four O’s per unit cell,
four states are found.
The occurrence of such unphysical (indeed, unchemical!) ghostbands is the first warning that
something went wrong. A more definite warning comes upon running LAPW2, where the
corresponding charge densities are calculated. If the contribution to the charge density from the
energy derivative of the basis function [the Blm coefficient in equ. 2.4,2.7] is significant (i.e. much
more than 5 per cent) then a warning is issued in LAPW2.
In the present example it reads:
Next we show tio2.output2 for the first of the ghost bands at -1.297 Ry. One sees that it
corresponds mainly to a p-like charge, which originates from the energy derivative part Q(UE) of
the Kohn-Sham orbital. Q(UE) contributes 40.1% compared with 8.5% from the main component
Q(U). Q(UE) greater than Q(U) is a good indication for a ghost band.
Another file in which the same information can be found is tio2.help031, since the ghost band
is caused by a bad choice for the Ti-p energy parameter:
Note again for L=1 the percentage of charge associated with the primary (APW) basis functions ul
(8.5%) versus that coming from the energy derivative component (40.1%).
12.1. GHOST BANDS 193
If a ghost band appears, one should first analyze its origin as indicated above, then use
appropriate local orbitals to improve the calculation and get rid of these unphysical states.
Do not perform calculations with “ghost-bands”, even when the calculation converges.
Good luck !
194 CHAPTER 12. TROUBLE SHOOTING
13 References
195
196 CHAPTER 13. REFERENCES
— 1975 Comp. Phys. Commun. 9, 31; this paper contains much of the Dirac-Fock treatment
used in Desclaux’s relativistic LSDA code.
Feldman J.L., Mehl M.J., and Krakauer H. 1987 Phys. Rev. B 35, 6395
Gay David M., ”ALGORITHM 611 – Subroutines for Unconstrained Minimization Using a
Model/Trust-Region Approach”, ACM Trans. Math. Software 9 (1983), pp. 503-524.
Hébert-Souche C., Louf P.-H., Blaha P., M. Nelhiebel, Luitz J., Schattschneider P., Schwarz K.
and Jouffrey B.; The orientation dependent simulation of ELNES, Ultramicroscopy, 83, 9
(2000)
“International Tables for X-Ray Crystallography“ 1964 Vol.1; The Kynoch Press,
Birmingham UK
Jansen H.J.F. and Freeman A.J. 1984 Phys. Rev. B 30, 561
Koelling D.D. and Arbman G.O. 1975 [Link]. F: Met. Phys. 5, 2041
Koelling D.D. and Harmon B.N. 1977 J. Phys. C: Sol. St. Phys. 10, 3107
Kohler B., Wilke S., Scheffler M., Kouba R. and Ambrosch-Draxl C. 1996
[Link]. 94, 31
Krimmel H.G., Ehmann J., Elsässer C., Fähnle M. and Soler J.M. 1994, [Link]. B50, 8846
Kuneš J, Novák P., Schmid R., Blaha P. and Schwarz K. 2001, Phys. Rev. B64, 153102
Liberman D., Waber J.T., and Cromer D.T. 1965, Phys. Rev. 137A, 27
Luitz J., Maier M., Hébert C., Schattschneider P., Blaha P., Schwarz K., Jouffrey B. 2001 Eur.
Phys J. B 21, 363-367
Madsen G. K. H., Blaha P, Schwarz K, Sjöstedt E and Nordström L 2001, Phys. Rev. B64,
195134
Mattheiss L.F. and Hamann D.R. 1986 Phys. Rev. B 33, 823
Mattsson A., Armiento R., Paier J., Kresse G., Wills J. and Mattsson T 2008 J. Chem. Phys.
128, 084714
197
Moruzzi V.L., Janak J.F., and Williams A.R. 1978 “Calculated Properties of Metals“
(Pergamon, NY)
Neckel A., Schwarz K., Eibler R. and Rastl P. 1975 [Link], Suppl.6, 257
Nelhiebel M., Louf P. H., Schattschneider P., Blaha P., Schwarz K. and Jouffrey B.; Theory of
orientation sensitive near-edge fine structure core-level spectroscopy, [Link]. B59, 12807
(1999)
Novák P. , Boucher F., Gressier P., Blaha P. and Schwarz K. 2001 Phys. Rev. B 63, 235114
Pardini L., Bellini V., Manghi F. and Ambrosch-Draxl C. 2011, [Link]. subm.
([Link]
Perdew J.P, Chevary J.A., Vosko S.H., Jackson K.A., Pederson M.R., Singh D.J., and Fiolhais
C. 1992 [Link].B46, 6671
Perdew J.P., Kurth S., Zupan J. and Blaha P. 1999, [Link]. 82, 2544
Ray A.K. and Trickey S.B. 1981 Phys. Rev. B24, 1751; erratum 1983, Phys. Rev. B28, 7352
Rondinelli JM, Beng Bin and Marks LD. 2007, Comp. Mater. Sci. 40, 345-353 (also: Los
Alamos archive, physics/0608160 ([Link]
Schwarz K. and Blaha P.: Lecture Notes in Chemistry 67, 139 (1996)
Singh D., Krakauer H., and Wang C.-S. 1986 Phys. Rev. B34, 8391
Singh D. and Nordström L 2006, Plane waves, pseudopotentials and the LAPW method, 2nd
edition, Springer, New York
Tran F, Laskowski R, Blaha P and Schwarz K. 2007, Phys. Rev. B 75, 115131
Tran F and Blaha P 2009, Phys. Rev. Lett. 102, 226401
von Barth U. and Hedin L. 1972 J. Phys. C.: Sol. St. Phys. 5, 1629
Wei S.H., Krakauer H., and Weinert M. 1985 Phys. Rev. B 32, 7792
Weinert M. 1981 J. Math. Phys. 22, 2433
Weinert M., Wimmer E., and Freeman A.J. 1982 Phys. Rev. B26, 4571
Wimmer E., Krakauer H., Weinert M., and Freeman A.J. 1981 Phys. Rev. B24, 864
Appendix
199
A Local rotation matrices
Local rotation matrices are used to rotate the global coordinate system (given by the definition of
the Bravais matrix) to a local coordinate system for each atomic site. They are used in the program
for two reasons:
I to minimize the number of LM combinations in the lattice harmonics expansion (of potential
and charge density according to equ. 2.10). For example for point group mm2 one needs for
L=1 just LM=1,0 if the coordinate system is chosen such that the z-axis coincides with the
2-fold rotation axis, while in an arbitrary coordinate system the three terms 1,0; 1,1 and -1,1
are needed (and so on for higher L).
I The interpretation e.g. of the partial charges requires a proper orientation of the coordinate
system. In the example given above, the p orbitals split into 2 irreducible representations,
but they can be attributed to pz and px , py only if the z-axis is the 2-fold rotation axis.
It is of course possible to perform calculations without “local rotation matrices“, but in such a case
the LM combinations given in Table 7.42 (and by SYMMETRY) may not be correct. (The LM
values for arbitrary orientations may be obtained from a procedure described in Singh 94.)
Fortunately, the “local rotation matrices“ are usually fairly simple and are now automatically
inserted into your [Link] file. Nevertheless we recommend to check them in order to be
sure.
The most common coordinate transformations are
Inspection of the output of SYMMETRY tells you if the local rotation matrix is the unit matrix or it
gives you a clear indication how to find the proper matrix.
The local rotation matrix R , which transforms the global coordinates r to the rotated ones r0 , is
defined by Rr = r0 .
There are two simple ways to check the local rotation matrixes together with the selected LM
combinations:
I charge density plots generated with LAPW5 must be continuous across the atomic sphere
boundary (especially valence or difference density plots are very sensitive, see 8.6)
I Perform a run of LAPW1 and LAPW2 using the GAMMA-point only (or a complete star of
another k point). In such a case, “wrong“ LM combinations must vanish. Note that the latter
is true only in this case. For a k mesh in the IBZ “wrong“ LM combinations do not vanish
and thus must be omitted!!
A first example for “local rotation matrices“ is given for the rutile TiO2, which has already been
described as an example in section 10.3. Also two other examples will be given (see below).
201
202 APPENDIX A. LOCAL ROTATION MATRICES
....
Titanium operation # 1 1
Titanium operation # 2 -1
Titanium operation # 5 2 || z
Titanium operation # 6 m n z
Titanium operation # 12 m n 110
Titanium operation # 13 m n -110
Titanium operation # 18 2 || 110
Titanium operation # 19 2 || -110
pointgroup is mmm (neg. iatnr!!)
axes should be: m n z, m n y, m n x
This output tells you, that for Ti a mirror plan normal to z is present, but the mirror planes normal
to x and y are missing. Instead, they are normal to the (110) plane and thus you need to rotate x, y
by 45◦ around the z axis. (The required choice of the coordinate system for mmm symmetry is
also given in Table 7.42)
....
Oxygen operation # 1 1
Oxygen operation # 6 m n z
Oxygen operation # 13 m n -110
Oxygen operation # 18 2 || 110
pointgroup is mm2 (neg. iatnr!!)
axes should be: 2 || z, m n y
For O the 2-fold symmetry axes points into the (110) direction instead of z. The appropriate
rotation matrices for Ti and O are:
√1 √1 −1 √1
0 0
√
2 2 2 2
−1 √1 √1 √1
√
2 2
0 0
2 2
0 0 1 1 0 0
A.2 Si Γ-phonon
Si possesses a face-centered cubic structure with two equivalent atoms per unit cell, at (±0.125,
±0.125, ±0.125). The site symmetry is -43m. For the Γ-phnon the two atoms are displaced in
opposite direction along the (111) direction and cubic symmetry is lost. The output of
SYMMETRY gives the following information:
Si operation # 1 1
Si operation # 13 m n -110
Si operation # 16 m n -101
Si operation # 17 m n 0-11
Si operation # 24 3 || 111
Si operation # 38 3 || 111
pointgroup is 3m (neg. iatnr!!)
axis should be: 3 || z, m n y
lm: 0 0 1 0 2 0 3 0 3 3 4 0 4 3 5 0 5 3 6 0 6 3 6
A.3. TRIGONAL SELENIUM 203
Therefore the required local rotation matrix should rotate z into the (111) direction and thus the
matrix in the “struct“ file should be:
√ √ √
6
0.4082483 -.7071068 0.5773503 √6
−√ 22 √3
3
6 2 3
0.4082483 0.7071068 0.5773503 6√ √2 √3
6 2 3
-.8164966 0.0000000 0.5773503 −2 6 2 3
H LATTICE,[Link]: 1
MODE OF CALC=RELA POINTGROUP:32
8.2500000 8.2500000 9.369000
ATOM= -1: X= .7746000 Y= .7746000 Z= 0.0000000
MULT= 3 ISPLIT= 8
ATOM= -1: X= .2254000 Y= .0000000 Z= 0.3333333
ATOM= -1: X= .0000000 Y= .2254000 Z= 0.6666667
Se NPT= 381 R0=.000100000 RMT=2.100000000 Z:34.0
LOCAL [Link]: 0.0 0.5000000 0.8660254
0.0000000 -.8660254 0.5000000
1.0000000 0.0000000 0.0
6 IORD OF GROUP G0
......
Se operation # 1 1
Se operation # 9 2 $|$$|$ 110
pointgroup is 2 (neg. iatnr!!)
axis should be: 2 || z
lm: 0 0 1 0 2 0 2 2 -2 2 3 0 3 2 -3 2 4 0 4 2 -4 2 ......
Point group 2 should have its 2-fold rotation axis along z, so the local rotation matrix can be
constructed in two steps: firstly interchange x and z (that leads to z k (011) ) and secondly rotate
from (011) into (001) (see the struct file given above). Since this is a hexagonal lattice, SYMMETRY
uses the hexagonal axes, but the local rotation matrix must be given in cartesian coordinates.
204 APPENDIX A. LOCAL ROTATION MATRICES
B Periodic Table
205
206
19K 20Ca 21Sc 22Ti 23V 24Cr 25Mn 26Fe 27Co 28Ni 29Cu 30Zn 31Ga 32Ge 33As 34Se 35Br 36Kr
2 1 2 2 2 3 2 5 1 5 2 6 2 7 2 8 2 10 1 10 2 10 2 10 2 2 10 2 3 10 2 4 10 2 5 10 2 6
Ar4s Ar4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s Ar3d 4s 4p Ar3d 4s 4p Ar3d 4s 4p Ar3d 4s 4p Ar3d 4s 4p Ar3d 4s 4p
37Rb 38Sr 39Y 40Zr 41Nb 42Mo 43Tc 44Ru 45Rh 46Pd 47Ag 48Cd 49In 50Sn 51Sb 52Te 53I 54Xe
2 1 2 2 2 3 2 5 1 5 2 6 2 7 2 8 2 10 1 10 2 10 2 10 2 2 10 2 3 10 2 4 10 2 5 10 2 6
Kr5s Kr5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s Kr4d 5s 5p Kr4d 5s 5p Kr4d 5s 5p Kr3d 5s 5p Kr4d 5s 5p Kr4d 5s 5p
55Cs 56Ba 57-71 72Hf 73Ta 74W 75Re 76Os 77Ir 78Pt 79Au 80Hg 81Tl 82Pb 83Bi 84Po 85At 86Rn
2 14 2 2 14 3 2 14 5 1 14 5 2 14 6 2 14 7 2 14 8 2 14 10 1 14 10 2 14 10 2 14 10 2 2 14 10 2 3 14 10 2 4 14 10 2 5 14 10 2 6
Xe6s Xe6s La-Lu Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s Xe4f 5d 6s 6p Xe4f 5d 6s 6p Xe4f 5d 6s 6p Xe4f 3d 6s 6p Xe4f 5d 6s 6p Xe4f 5d 6s 6p
57La 58Ce 59Pr 60Nd 61Pm 62Sm 63Eu 64Gd 65Tb 66Dy 67Ho 68Er 69Tm 70Yb 71Lu
2 2 2 3 2 4 2 5 2 6 2 7 2 7 2 9 2 10 2 11 2 12 2 13 2 14 2 14 2
Xe5d6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 5d6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 6s Xe4f 5d6s
89Ac 90Th 91Pa 92U 93Np 94Pu 95Am 96Cm 97Bk 98Cf 99Es 100Fm 101Md 102No 103Lr
2 2 2 2 1 2 3 1 2 4 1 2 6 2 7 2 7 2 9 2 10 2 11 2 12 2 13 2 14 2 14 2
Rn6d7s Rn6d 7s Rn5f 6d 7s Rn5f 6d 7s Rn5f 6d 7s Rn5f 7s Rn5f 7s Rn5f 6d7s Rn5f 7s Rn5f 7s Rn5f 7s Rn5f 7s Rn5f 7s Rn5f 7s Rn5f 6d7s
APPENDIX B. PERIODIC TABLE
WIEN2k
ISBN 3-9501031-1-2