Purpose: This tutorial gives a basic introduction into electronic-structure calculations. It explains how to set up and execute a simple exciting calculation, using elemental Ag as example. It is described how to prepare the input, how to run the calculation, and how to analyze the output. In addition, it is shown how basic properties like the density of states, and the band structure.
|
Table of Contents
|
0. General Preparations
Units in exciting: By default, all quantities in the exciting code are given atomic units: Energies in Hartree, lengths in Bohr, etc. In case other units are desirable, they can be converted using templates as a post-processing to exciting's standard output.
Very important: Before starting, the following shell variable must be set by the user:
EXCITINGROOT = Directory where exciting has been downloaded, e.g.: /home/user/exciting .
The setting of this variable can be done in a bash shell by typing (from now on the symbol $ will indicate the shell prompt):
$ export EXCITINGROOT=/local_path_to_exciting_root
Please note: In the input-files shown on this page, the placeholder EXCITINGROOT always needs to be replaced by the value of the $EXCITINGROOT variable.
1. Electronic structure of silver
1.1 Groundstate calculation
The first step of any density-functional calculation is the determination of the groundstate total energy and electron density.
The starting point of a groundstate calculation is the crystal structure, only. At the beginning of a groundstate calculation, an initial electron density is generated, which is obtained from a superposition of atomic densities. Thus, this initial electron density lacks electron-electron and electron-ion interactions between atoms and is normally a rather crude approximation of the density.
Then, the calculation iteratively goes through the following steps:
- Determine the potential from the electron density.
- Solve of the Kohn-Sham (KS) equations to get the eigenfunctions and eigenvalues as well as the total energy.
- Calculate the electron density from the KS eigenfunctions.
- Create a new charge density, mixing the electron density from the current iteration with the ones of previous iteration (to ensure a good convergence behavior).
- Start again with (1).
Such a sequence of steps is usually called an iteration. The code will repeat such iterations, until the potential (or total energy, or charge density, …) obtained at the end of the last iteration is consistent with the one of the previous iteration. Thus, this kind of calculations is often called self-consistent field (SCF) calculation, and an iteration is often referred to as an scf cycle.
To prepare your calculation, create a new, empty directory named Ag/ somewhere on your filesystem. In this directory, save the following lines
as input.xml (do not forget to replace EXCITINGROOT by the value of the $EXCITINGROOT variable!):
<?xml version="1.0" encoding="UTF-8"?> <?xml-stylesheet href="http://xml.exciting-code.org/inputfileconverter/inputtohtml.xsl" type="text/xsl"?> <input xsi:noNamespaceSchemaLocation="EXCITINGROOT/xml/excitinginput.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsltpath="EXCITINGROOT/xml/"> <title>Ag</title> <structure speciespath="EXCITINGROOT/species"> <crystal scale="7.7201"> <basevect>0.5 0.5 0.0</basevect> <basevect>0.5 0.0 0.5</basevect> <basevect>0.0 0.5 0.5</basevect> </crystal> <species speciesfile="Ag.xml"> <atom coord="0.0 0.0 0.0" /> </species> </structure> <groundstate ngridk="8 8 8"></groundstate> </input>
If the visualization program XCrySDen is set up appropriately (find here how to do this: XcrysdenExcitingSetup) you can visualize the structure in an exciting input.xml executing
$ xcrysden --exciting input.xml
After this, start the groundstate calculation by executing the following command in the Ag/ directory:
$ $EXCITINGROOT/bin/excitingser
The calculation should roughly take 1 minute. During the calculation, output files are created, which contain all kind of information on your material system and on the calculation. Some of the output files are already created at the beginning of the calculation and will not be changed anymore during the run. The most important among them are:
| filename | description |
|---|---|
| LATTICE.OUT | Information on the lattice: Primitive lattice vectors, unit cell volume, reciprocal lattice vectors, etc. |
| SYMCRYS.OUT | Information on the symmetry operations of the crystal; more symmetry information is found in the files SYMT2.OUT, SYMSITE.OUT, SYMMULT_TABLE.OUT, SYMMULT.OUT, SYMLAT.OUT, SYMINV.OUT, and SYMGENR.OUT. |
| KPOINTS.OUT | List of k-points, their coordinates (in units of the reciprocal lattice vectors), weights, matrix size. |
| IADIST.OUT | Interatomic distances; useful to check the correctness of an input file. |
| EQATOMS.OUT | Information on equivalency of atoms due to the crystal symmetry. |
| geometry.OUT.xml | Structural information on your system. This will often be identical to the element structure in your input file, but may differ in case you set the groundstate attribute primcell = "true". |
Other files are updated or extended in each iteration, and contain information about the scf calculation. Here are the most important ones:
| filename | description |
|---|---|
| info.xml | "Master" output file in XML format, containing the essential information on the material system, on the parameters of the calculation, on the results (total energy, energy contributions, charge contributions, Fermi energy …) of each iteration, and more … |
| INFO.OUT | Extended version of info.xml, displayed as plain text file. |
| TOTENERGY.OUT | Total energy (Hartree); each line corresponds to one iteration. |
| EFERMI.OUT | Fermi energy (Hartree). |
| FERMIDOS.OUT | Density of states at the Fermi level (states/unit cell/Hartree); each line corresponds to one iteration. |
| EIGVAL.OUT | Eigenvalues (energies) of the valence bands, for each k-point and band. |
| EVALCORE.OUT | Eigenvalues (=energy levels) of the core states. |
| LINENGY.OUT | Linearization energies as fixed in the species files (if searchE="false" for the corresponding linearisation energy in the species.xml file) or determined by exciting (if searchE="true" for the corresponding linearisation energy in the species.xml file). |
Exercises
With the table above, find out the following properties of your calculation — to do so, find out first in which output files they are contained:
- How many iterations did the calculation go through?
- What is the total energy for the first iteration (which started from the superposition of atomic electron densities), and of the converged calculation (last iteration)?
- What was the change in total energy between the
- first two iterations?
- last two iterations?
- What is the Fermi energy of the system?
- How many occupied valence bands are there in the system?
- How much charge is there inside the muffin-tin sphere, and how much is found in the interstitial?
1.2 Density of states
After you have completed the groundstate run and have obtained the corresponding total energy, you can go for more properties of the system. One of the most fundamental ones is the density of states (DOS). The DOS gives you information on the energy levels in your system, or — more precisely — how many electronic states there are at any given energy.
To calculate it, you need to do the following to simple modifications in input.xml:
- add the attribute do="skip" to the xml-element <groundstate>;
- add the element <properties> … </properties> after the <groundstate> element;
- insert the subelement <dos> … </dos> into the element <properties>;
- add the attribute nsmdos="1" to the element <dos> (see Input Reference for more details).
The corresponding part of the input.xml should now look like this:
... <groundstate ngridk="8 8 8" do="skip"></groundstate> <properties> <dos nsmdos="1"> </dos> </properties> ...
Then execute excitingser again on the command line:
$ $EXCITINGROOT/bin/excitingser
This time, the program will produce the following files:
- TDOS.OUT … total DOS
- PDOS_S01_A0001.OUT … partial density of states for each atom, l and m quantum number.
- IDOS.OUT … interstitial DOS
- dos.xml … all densities of states assembled in XML format.
Let's have a look at the total DOS in TDOS.OUT. To visualize it, execute
$ xmgrace TDOS.OUT
This will open the plotting tool xmgrace and display the total density of states. The units in this plot are
- Energy (Hartree) for the x axis
- DOS (states/unit cell/Hartree) for the y axis.
Inside xmgrace, you can change the plot appearance in any way you want, zoom in to the see details of the DOS, or produce a figure in any format you like (ps, jpg, png, etc. …). See Xmgrace Quickstart for more info.
Exercises
- We were using the attributes do="skip" for the element <groundstate> for generating the DOS after the groundstate SCF run. Find out why, by searching for the element groundstate in Input Reference, and then going down to its attribute do.
- Consider the electronic configuration (e.g., see http://www.webelements.com/silver/atoms.html -> Ground state electron configuration) and the electron binding energies (http://www.webelements.com/silver/orbital_properties.html -> Electron binding energies) of atomic silver:
- Which states do you expect to form the valence band in the Ag compound?
- Which states are there with energy greater than 3 Hartree? (The conversion factor between Hartree and eV is found in Input Reference.)
- What about the silver 4p states? Try to find out where they are located in the Ag compound. To do so, you will need to use the attribute winddos of the element <dos> — check out Input Reference to see why, and how this works.
- Do the following modifications inside xmgrace (see Xmgrace Quickstart for help …):
- change the x-region of the plot to -0.4 — 0.25
- autoscale x axis
- add axes labels
- add a title
- increase the line-width of the DOS curve to 3
- change the colour of the DOS curve to red
- add a legend
- generate the file Ag_dos.png.
1.3 Band structure
Now we are ready for a more detailed view on the electronic structure: The band structure. In addition to the energy of each state, the band structure shows the dependence of the energy eigenvalues on the coordinates in k-space.
To calculate the band structure of silver, insert the subelement <bandstructure> in the element <properties> with the following specifications:
... <properties> <bandstructure> <plot1d> <path steps="100"> <point coord="0.75000 0.50000 0.25000" label="W" ></point> <point coord="0.50000 0.50000 0.50000" label="L" ></point> <point coord="0.00000 0.00000 0.00000" label="GAMMA"></point> <point coord="0.50000 0.50000 0.00000" label="X" ></point> <point coord="0.75000 0.50000 0.25000" label="W" ></point> <point coord="0.75000 0.37500 0.37500" label="K" ></point> </path> </plot1d> </bandstructure> </properties> ...
As you may have realized, we have removed the subelement <dos> now.
Now execute excitingser again on the command line:
$ $EXCITINGROOT/bin/excitingser
This makes the code produce the band structure, which is written to bandstructure.xml. To visualize the band structure, execute the following command:
$ xsltproc $EXCITINGROOT/xml/visualizationtemplates/xmlband2agr.xsl bandstructure.xml
Now you have produced the xmgrace file Ag_bandstructure.agr, which you can open, visualize and manipulate with xmgrace:
$ xmgrace Ag_bandstructure.agr
The result should look like this:
Exercises
- Use http://www.cryst.ehu.es/ -> Space Groups Retrieval Tools -> KVEC to find out about the location of the special k-points within the Brillouin zone. The spacegroup of Ag is 225, Fm-3m. Select Choose to choose the corresponding spacegroup, and then click [ Brillouin zone ] to see the Brillouin zone with the special k-points.
- The band structure of our plot is plotted along the path W -> L -> Γ (Gamma) -> X -> W -> K. In which parts of this path do you find strong or weak dispersion, respectively? Try to explain why by using the Brillouin zone plot from above.
- In xmgrace: Use the autoscale button to see all bands. Look at the dispersion of the bands for low energies and high energies:
- What trend do you see relating the band width to the energy?
- How can you explain this trend (think about how the "dispersion" of an isolated atom would look like …)?
1.4 Convergence issues — crucial parameters in a calculation
A groundstate calculation using a DFT code fundamentally depends on 2 parameters:
- the mesh of k-points (<groundstate> attribute: ngridk);
- the size of the basis set for expanding the wave function (<groundstate> attribute: rgkmax).
In order to be able to rely on your calculation, you need to understand what these two parameters mean and make a convergence tests for each new system you are dealing with.
Let us first have a closer look at the meaning of these two parameters:
The k-point mesh
From basic considerations of solid-state theory, it can be shown that the Schrödinger (or, in our case, Kohn-Sham) equation for a periodic system has a bunch of solutions, where each solution is characterized by a vector k in reciprocal space. This k-vector is essentially related to the periodicity of the corresponding solution of the Schrödinger equation: Roughly speaking, the solution will have the same value at coordinates in real space that, along the direction of k, have a distance of 2π/|k|.
Many properties of a solid, including the total energy, are represented as integrals, performed over all possible k-vectors. Obviously, the direct calculation of these integrals is very demanding (one should calculate the solution of the Schrödinger equation for a huge number of k-points). Therefore, the integrals are approximated by sums on a set of k-points distributed on a finite grid. The spacing of the points on the grid is a measure of the accuracy of the calculation of the integrals which also depends on how fast the integrand (the solution of the Schrödinger equation) varies by changing k. The challenge is to find a good number of k-points: Large enough to capture the physical properties of your system well, but small enough to keep calculations as fast as possible.
In fact, the number of k-points in a given direction expresses, over how many unit cells electrons feel each other in this direction. This has two important consequences:
- The larger the unit cell, the smaller is the required k-point mesh. As a rool of thumb, for a given kind of material, two calculations have the same quality with respect to the k mesh, if the product (no. of atoms)*(no. of k-points) is the same for both of them. For example: If you create a supercell of your starting structure by doubling it in z direction, it is enough to use only half of the k-points along z.
- Systems with longer-range interactions need larger k-point meshes.
In metals you need a large number of k-points, since the conduction electrons are delocalized over the whole system, i.e. the interactions are very long-ranged. Moreover, for a good description of the system it is important to know exactly where the conduction bands cross the Fermi level, which also requires a dense k mesh.
In contrast, semiconductors and insulators usually have much more localized electronic states, and a gap between valence and conduction bands. Thus, the number of k-points required for a good calculation is much smaller.
While for a metal with 1 atom per unit cell, up to ngridk="20 20 20" or larger may be necessary to get reliable physical properties, for a semiconductor with the same unit-cell size ngridk="4 4 4" may already do a good job. If you want to produce high-quality results for publication, a careful convergence test is indispensable. Please note, that, for metals, this test must go hand in hand with a test of the swidth. Check out the Input Reference for the <groundstate> attributes stype and swidth to find out more about this topic.
The basis-set size
Most of the DFT codes solve the Kohn-Sham equation in reciprocal space, expanding the wavefunctions in terms of suitable periodic basis functions. In our case, the basis functions are the linearized augmented plane waves (LAPW). The larger the size of the basis set, the more accurate is the calculation. The corresponding «groundstate» attribute in the exciting code is rgkmax (see Input Reference). The default value of rgkmax is 7.0. In many cases (e.g., for a good description of forces or equations of state), larger values (8.0 — 9.5) may be necessary to get reliable results. Therefore, for this parameter it is very important to test how large rgkmax must be in order to get the physical property under investigation with the desired accuracy.
Terminology: When speaking about the number of the basis functions, you may hear the expressions "rgkmax", "basis-set size", "energy cut-off", or "plane-wave cut-off" — all of them are a measure (in different units) of the same quantity.
Scaling of the calculation time…
The calculation time scales
- linearly with respect to the number of k-points;
- as rgkmax to the power of 9.
Exercises
- Convergence of the k-mesh: Repeat the calculation of groundstate and DOS using for ngridk the values "3 3 3" ,"6 6 6", …, "21 21 21" in steps of 3. The best way to do so is to (1) create a subdirectory for each k-mesh (e.g. dir_nkx_3/ for ngridk="3 3 3"), (2) copy the input file there and modify the attribute ngridk appropriately, (3) start the calculation in each subdirectory by executing $EXCITINGROOT/excitingser on the commmand line.(Please note: Since both the groundstate and the DOS should be calculated at once, now, set the attribute do="fromscratch" for the <groundstate> element).
- Extract the total energy from the last iteration of each groundstate calculation and plot total energy vs. nkx, where nkx is the number of k-points in x direction (e.g. 3 in the case of ngridk="3 3 3"). Taking the total energy of the calculation with ngridk="21 21 21" as a reference: What k-mesh is needed to get a total energy within 10-3 Hartree from the converged result?
- Plot the total DOS (file TDOS.OUT) for all k-meshes (You can do so by invoking xmgrace with all TDOS.OUT files at once, see Xmgrace Quickstart. Sect. 6 ). For which value of ngridk can the DOS be considerd to be converged?
- Convergence of the basis-set size: Do calculations with a fixed value of ngridk="12 12 12", but varying rgkmax from 5 to 10 in steps of 1. Make the same checks as in exercise about the k-mesh, and try to find out for which values of rgkmax the total energy and DOS, respectively, can be considered as converged.