First steps
In this tutorial we will install Fleur and run a few basic calculations to show you several basic workflows when using Fleur.
1. Installation
This exercise will be performed together with a presenter who demonstrates each step. In general information about the Fleur installation process can be obtained at the respective page.
2. A first calculation
For the first calculation we choose a perfect Si crystal in diamond structure as example system. The inpgen input for such a system is:
Si bulk &lattice latsys='cF', a0=1.8897269, a=5.43 / 2 14 0.125 0.125 0.125 14 -0.125 -0.125 -0.125
The first line in this input is just a comment line. The '&lattice' line defines the lattice. With 'latsys='cF' we here specify a cubic close-packed lattice. 'a=5.43' specifies the lattice constant. inpgen and fleur in general assume atomic units but in this case we provide the lattice constant in Angstrom. Therefore we need a conversion factor from Angstrom to atomic units. This is specified by 'a0=1.8897269'.
The next line contains the number of atoms in the unit cell. It is followed by a list containing for each atom the atomic number (14 for Si) and the relative position in the unit cell.
A documentation of the general layout of inpgen inputs is provided at the respective page.
Create a directory for this calculation and in it a text file 'inpSi.txt' (or similar) with the content above. To generate the Fleur input invoke in the new directory:
pathToInpgen/inpgen < inpSi.txt
Several files are created. Of these the input files for the Fleur calculation are 'inp.xml' for the full parametrization of the calculation and 'sym.out' for a list of symmetry operations present in the crystal. The file 'struct.xsf' is an XCrysDen structure file that can be used to visualize the unit cell. 'out' is the general text output of inpgen (or fleur after executing fleur in the directory). If something went wrong in the generation of the Fleur input it is a good idea to have a look in that File to see whether there is a hint what was wrong. 'FleurInputSchema.xsd' is not relevant to the user. It is a general specification of the inp.xml file format in terms of an XML Schema definition.
Have a closer look at 'inp.xml' and identify in it the different cutoff parameters and other numerical parameters used for the FLAPW method and for DFT in general. We will discuss the contents. The proper choice of the parameters is not a topic of this tutorial. Here we will mostly rely on the default values. If you perform real research with Fleur make sure that you obtain results that are also converged with respect to the parameters. In this older tutorial the different parameters are investigated in more detail (Note that it needs some adaptions to the here used Fleur version).
We invoke fleur in the directory with the 'inp.xml' and 'sym.out' files:
pathToFleur/fleur
Note that we use 'fleur' as an acronym for the Fleur executable. Typically this either refers to 'fleur' if the program is compiled without MPI support or 'fleur_MPI' if MPI support is available. Note: MPI is an acronym for Message Passing Interface. It relates to a certain type of parallelization scheme which is especially needed for a parallelization over multiple compute nodes. Beyond MPI and even without MPI support Fleur also features an OpenMP parallelization that can be used for the parallelization on each compute node. The degree of OpenMP parallelization is typically controlled by external environment variables, e.g., 'OMP_NUM_THREADS'. In practice it turns out that the most efficient parallelization uses a mixture of OpenMP and MPI on each compute node.
In practice DFT is implemented as an iterative algorithm that starts with a first guess for the electron density and ends after several iterations with a self-consistent density. In Fleur by default up to 15 iterations of the self-consistency loop are performed (can be changed in inp.xml - /fleurInput/calculationSetup/scfLoop/@itmax parameter). You can observe the development of the distance between the input and output densities of each iteration in the terminal output (not in this tutorial because the used compiler redirects this output to the out file). Alternatively this can also be obtained after the calculation by invoking 'grep dist out' to find the respective entries in the generated out file. You will find that actually less than 15 iterations were performed because the convergence criterion in /fleurInput/calculationSetup/scfLoop/@minDistance was fulfilled. In other calculations that are not converged after the default number of iterations we can invoke Fleur again to perform more iterations or adjust itmax prior to the Fleur run.
The output of the fleur calculation is available in the 'out' file and also in the 'out.xml' file. Investigate these files and find out what quantities are written out. Note that certain calculations or unit cell setups may trigger additional output. We will discuss the contents of the files. Note that whenever Fleur is executed an old out file that may be present in the directory is overwritten.
3. The Si lattice constant
To calculate the equilibrium lattice constant for a material typically the total energies for several test lattice constants are calculated. The equilibrium lattice constant minimizes the total energy. This minimum can be determined by fitting the Birch-Murnaghan equation of state to the obtained total energies. In this tutorial we make a less sophisticated approach and just guess the minimum from a plot of the total energies.
To set up the calculations for the different test lattice constants we start with the Si inp.xml file and modify it slightly to allow for the calculation of the different lattice constants. Then we adapt it to each of the lattice constants. Note: If calculations with differing calculational parameters (e.g. different lattice constants) are to be compared it is important to base the different calculations on the same inp.xml file and not on the same inpgen input. If parameters are varied in the inpgen input this may lead to inp.xml files that differ by more than only the considered parameter. In such a case one may observe differences in the results of the Fleur calculations that are due to unintended variations in the Fleur input.
We create a new directory for the determination of the Si lattice constant and copy the inp.xml and sym.out file from the previous calculation to it. We want to test lattice constants in 1% steps from 97% to 103% of the lattice constant we initially provided to inpgen. It has to be ensured that the input works with each of the intended calculations. In this case this implies that we have to adapt the Si MT radius. Since the MT spheres of the atoms in FLAPW calculations nearly touch a change in the lattice constant might lead to overlapping MT spheres: This would be an invalid setup. We use the fact that we can use mathematical expressions within the inp.xml file wherever a floating point number is required. Replace the entry in atomSpecies/species/mtSphere/@radius for the Si-1 species by '0.97*2.17000000'. Furthermore we don't want to invoke Fleur multiple times for each lattice constant. Therefore we set the number of iterations in calculationSetup/scfLoop/@itmax to 40. This is a save value that will never be reached since the convergence criterion in calculationSetup/scfLoop/@minDistance will cause a stop of the Fleur calculation as soon as convergence is reached.
Generate the subdirectories latt-0-97 to latt-1-03 and copy inp.xml and sym.out into each of them. In each directory we can easily adjust a scaling factor in the respective inp.xml file to obtain the desired lattice constant. This scaling factor is found in the XML attribute cell/bulkLattice/@scale. We perform the 7 calculations by calling fleur in each of the subdirectories.
Check whether all calculations are converged, e.g., by looking at the direct Fleur output or by calling "grep dist out". The total energies can be found for each iteration in the out file (grep "total energy" out) or in the out.xml file (grep totalEnergy out.xml). We collect them in a file totalEnergies.txt with two columns: 1. scaling factor, 2. converged total energy in Htr.
Plot totalEnergies.txt with some plotting tool like gnuplot, xmgrace, or Matplotlib and guess the equilibrium lattice constant on the basis of these plots.
Note: Within gnuplot an adequate plotting command is
plot 'totalEnergies.txt' u 1:2 w lines
4. The Si band structure and density of states
Consider your self-consistent calculation nearest to the equlibrium lattice constant and create the two subdirectories 'bands' and 'dos'. Copy all files from the self-consistent calculation to these subdirectories. In 'bands' we will perform a band structure calculation, in 'dos' we will calculate the density of states.
For the band structure calculations a special k point set has to be chosen that is distributed along the desired k point path. For common k point paths this is automatically done by fleur as long as the k point set is provided by a kPointCount XML element and not as a mesh or as a list of k points. The current version of inpgen provides such a k point set specification in the optional XML element /fleurInput/calculationSetup/bzIntegration/altKPointSet with the purpose "bands". This alternative k point set is used by fleur whenever a band structure is calculated. If this optional XML element is not present the standard k point set is used for the band structure calculation and the user has to make sure that this k point set is adapted to the needs of the calculation. For more information of how to specify a k point set consult the respective section in the inp.xml file reference.
To activate the band structure calculation mode in fleur the XML attribute /fleurInput/output/@band has to be set to T.
Do this and run fleur. You will obtain several files:
| day | topic |
|---|---|
| band.1 | This text file contains the eigenvalues (column 2) together with the x coordinate (column 1) of the band structure plot. |
| band.gnu | This file contains a gnuplot script that can be used for the visualization of the band.1 file. |
Note that in spin-polarized calculations also a band.2 file is generated. In such a case the two band.X files contain the data for the two different spins.
Use the script with
gnuplot < band.gnu > bands.ps
You can either visualize the resulting postscript file with a viewer like gv (ghostscript) or convert it into a pdf file with
ps2pdf bands.ps
The outcome of this conversion can then be visualized with a pdf viewer.
Note that only basic assumptions are made for the definition of the Fermi energy in the band structure files. If you have a fleur version that is compiled with HDF5 support the Fermi energy of the last iteration before the band structure calculation is used, otherwise the Fermi energy of the band structure calculation. Neither of them is always correct since the physically highest occupied state may either be found on the high-symmetry k point path or not. Especially for metals it is not always obvious which Fermi energy is the better choice. In such a case the k point set of the self-consistency calculation should be converged until nothing changes anymore. You can adjust a wrong Fermi energy in the band.gnu plotting script.
For the density of states (DOS) calculation the k point set should sample the Brillouin zone very well. Typically this means that a DOS calculation should not be performed with the same k point set as the band structure calculation. In practice you can take the k point set from the self-consistency calculation (or even refine it slightly).
To activate a DOS calculation you have to set /fleurInput/output/@dos to T and /fleurInput/output/densityOfStates/@ndir to -1. See the page on the integrated DOS program for details on the ndir switch.
Perform the calculation. You now have a file DOS.1 (and for spin-polarized calculations also a file DOS.2 for the other spin). In this file you finde several columns with numbers. The first column is an energy mesh in eV relative to the Fermi energy. In the other columns for each energy you find the total DOS, the DOS in the interstitial region, the vacuum DOS for vacuum 1 and 2 (filled with 0.0 if the unit cell does not represent a film system), multiple columns for the DOS at each atom type (atom group in the inp.xml file), and finally multiple columns for the s,p,d,f projected DOS at each atom type (grouped for each atom type). Note that the selection of the correct columns for a plot is error-prone. Be careful!
Use gnuplot, or another plotting tool, to plot for the given energy mesh the total DOS, the interstitial region DOS, the DOS at the Si atom, and the s-, and p-projected DOS at the Si atom.
You will observe that the DOS is strongly smeared out. The smearing is required because the k point sampling of the Brillouin zone is finite. But we may try out a slightly finer smearing. This is controlled by the /fleurInput/output/densityOfStates/@sigma parameter. Divide it by 3, rerun the calculation, and investigate the differences in the plotted DOS. You may notice that the DOS of the latter calculation features some small spikes. In a third DOS calculation we use the refined sigma parameter but also increase the k point sampling. Set /fleurInput/calculationSetup/bzIntegration/kPointCount/@count to 100 and rerun the calculation. How does the DOS change? Why does it change?
5. Relaxing the atom positions in a single MoS2 layer
MoS2 is a transition metal dichalcogenide that crystallizes in a layered structure in which each layer features a honeycomb lattice with a 3-atomic basis. The interlayer-bonding has a Van-der-Waals nature which tranlates into a large interlayer distance. In this exercise we consider a single MoS2 layer as a thin film system. The following inpgen input describes the respective unit cell:
MoS2 monolayer &input film=T cartesian=f symor=t/ 0.5 -0.5 0. 0.5 0.5 0. 0. 0. 1. 5.952636864 1. -3. 30.00 3 42 0.333333 0.666667 0.00 16 0.666667 0.333333 2.85 16 0.666667 0.333333 -2.85 &factor 1. 1. 1. / &end /
The '&input' namelist here specifies that this is a film system ('film=T'). We are using internal coordinates for the atom positions ('cartesian=f') and only want to use symmorphic symmetry operations ('symor=t'). The latter switch is used because non-symmorphic symmetry operations are not compatible to all lattices in Fleur. If you set up own input files and obtain a Fleur error message about this you now know what to do.
Next comes the Bravais matrix which is scaled according to the two lines following it. The first line is a general scaling parameter for the Bravais matrix, the second line specifies scaling parameters for each column of the Bravais matrix separately. Note that a scaling of '-3' actually means a scaling with . The last entry ('30.00') is just a large number which only ensures that the film fits into the here-specified unit cell. It is automatically reduced to a reasonable value by the program.
For film systems the last atom coordinate is always cartesian, independent of what was declared in the '&input' namelist. This is why we have larger numbers in comparison to the other two coordinates which are provided in terms of Bravais lattice vectors. There are more scaling parameters for the atom positions below the last atom position entry. The atom positions are divided by these numbers. Here these scaling parameters are .
Note that the x and y atom positions are supposed to be at and . We provide these numbers with a precision better than . With such a precision they are automatically turned into the exact values. A less precise specification of the atom positions will lead to a reduced set of symmetry operations being detected. A consequence of this is the generation of a k-point set that breaks the non-detected symmetries resulting, e.g., in the calculaton of wrong forces that are not compatible with all symmetries that are actually present in the system. Of course, such a problem can be overcome by using a very fine k-point mesh, but it is more efficient and precise to make sure that all relevant symmetry operations are included in the calculation. You can also directly provide the exact fractions by using the scaling factor. With this the atom positions would look like:
42 1.0 2.0 0.0 16 2.0 1.0 2.85 16 2.0 1.0 -2.85 &factor 3.0 3.0 1.0 /
The z-position of the sulfur atoms provided in this inpgen input is only a guess. In this exercise we will optimize this position with repect to the forces on the atoms and compare the band structures with and without the relaxation.
Use the inpgen input to generate a Fleur input for the MoS2 monolayer.
The calculation of forces is based on the presence of self-consistent densities. Therefore we will first adjust a few parameters in the inp.xml file to allow possible movements of the atoms, then perform a self-consistency calculation and finally calculate the forces and suggested atom displacements. The repetition of several self-consistency calculations (for every new atom displacement), each followed by a force calculation step is then automatized by a small script.
Note that formally Fleur does not require the presence of a self-consistent density to calculate the forces. Using Fleur in such a way will result in an automatic self-consistency loop until the forces seem to be converged. But this force-convergence may be misleading or inaccurate. As a result the number of force steps needed to reach equilibrium atom positions may increase.
Perform the following changes in the inp.xml file:
- Set the number of iterations (itmax) to 50. The convergence here is not problematic. Therefore this just leads to the convergence criterion ending the calculation.
- Multiply the MT radii for the atoms by .
Run Fleur to reach self-consistency. After this has succeeded copy the files to two subdirectories 'bands' and 'relax'.
Perform a band structure calculation in the 'bands' folder.
In the 'relax' folder we now modify the inp.xml file to perform the force calculations and generate the atom displacements. For this set the /fleurInput/calculationSetup/geometryOptimization/@l_f switch to "T".
In the geometryOptimization XML element you can also change the relaxation algorithm to something different (we keep BFGS), change the mixing parameter, or adjust the criteria for a successful relaxation. Note that you can limit the force calculations for each atom type to certain orientations. This is done in the /fleurInput/atomGroups/atomGroup/force/@relaxXYZ XML attribute. Each logical in this variable stands for one direction. But for this exercise we don't use this switch.
Invoke Fleur to calculate the forces for this configuration.
You will obtain a file 'relax.xml'. In every Fleur run this file is automatically included in the inp.xml file whenever it is present. Investigate the file. You will first find a list of displacements for each atom group. This is used to modify the atom positions in the inp.xml file in Fleur calculations. The displacements are updated by every new force calculation. The other big block in the relax.xml file is the relaxation history. Here you find for every relaxation step and atom type a list of atom positions in cartesian coordinates followed by the forces on the representative atom of this atom type in Hartree per . You can use this list as an output, e.g., look up the final atom positions in it, but it is also an input for Fleur that is used by the relaxation algorithm to determine the next atom positions. An experienced user may modify this file whenever a relaxation process seems to be problematic. For example in certain situations it may be a good idea to delete a part of the relaxation history.
We now have to repeat the self-consistency step followed by the force step several times to reach the equilibrium atom positions. Since we don't want to do this manually we write a small shell script that performs these steps. Note that a systematic approach to automatization would be the use of Aiida workflows. This is something we will address in one of the following days. Note also that the Fleur version we are using for this tutorial is compiled with HDF5 support. Without HDF5 support the output files would be slightly different and the user would have to delete the 'stars' and 'wkf' files after every force calculation step. These files depend on the atom positions.
Copy the inp.xml file to a file 'inp-scf.xml' and to a file 'inp-force.xml'. Adjust these files such they only differ in the l_f switch which is set to F in the inp-scf.xml file and to T in the inp-force.xml. Copy the following script to a file runScript.sh:
cp inp-scf.xml inp.xml Fleur cp inp-force.xml inp.xml Fleur cp inp-scf.xml inp.xml Fleur cp inp-force.xml inp.xml Fleur cp inp-scf.xml inp.xml Fleur cp inp-force.xml inp.xml Fleur cp inp-scf.xml inp.xml Fleur cp inp-force.xml inp.xml Fleur cp inp-scf.xml inp.xml Fleur cp inp-force.xml inp.xml Fleur
Replace in the script 'Fleur' by the Fleur call you perform on your system. Make the script executable by running the command 'chmod a+x runScript.sh'.
Start the relaxation by executing './runScript.sh' and wait until it is finished (it will take a few minutes). Note that Fleur indicates in the ending message of the force step whether more force steps are required (Structural relaxation: new displacements generated) or the atoms have reached their equilibrium positions (Structural relaxation: Done).
Find the sulfur z-positions of the last force step in the relax.xml file and compare them to the original positions. How much did they change? What is the energy difference between the two configurations?
Create a subdirectory bands in the relax folder and copy everything from the relax folder into it.
Note: The cdn.hdf file has become quite large because many densities are stored in it. you can infestigate it by calling fleur with the "-info" command line switch. You will obtain a list of the densities together with time stamps and the distance to the previous density. Of course, there are ways to get rid of certain densities. One way is to perform the Fleur calculations with the command line switch "-last_extra". This will write another file out in every iteration that only contains the last charge density. Another way is to use Fleur command line options to actually delete densities from the cdn.hdf file. A list of all Fleur command line options can be obtained with the command line option "-h". But in this exercise we will not manipulate the cdn.hdf file in such a way.
Copy the inp-scf.xml file to inp.xml and adjust the inp.xml file such that a band structure calculation can be performed. Calculate the band structure.
Compare the band structures for the guessed atom positions and the relaxed atom positions. Where are the valence band maxima and conduction band minima in each case? How large are the band gaps for each configuration. If everything worked correctly you can observe that the small change in the atom position has a rather large effect on these quantities.