by Pasquale Pavone & Caterina Cocchi for exciting nitrogen
Purpose: In this tutorial you will learn how to set up and execute exciting calculations for graphene. Examples of electronic band-structure calculations and full structure optimization will be also shown. Finally, the investigation of the loss function of graphene will be presented as an example for the calculation of excited-states properties.
|
Table of Contents
|
0. Preliminary notes
Read the following paragraphs before getting started!
Be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts that are relevant for this tutorial, together with a short description.
- SETUP-elastic-strain.py: Python script for generating strained structures.
- SETUP-planar.py: Python script for setting calculations for two-dimensional materials.
- SETUP-graphene-along-c.py: Python script for generating graphene structures in which an atom is displaced along the c axis.
- EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
- EXECUTE-elastic-strain.sh: (Bash) shell script for running a series of exciting calculations.
- EXECUTE-diamond-phonon-g.sh: (Bash) shell script for running a series of exciting calculations.
- PLOT-birch.py: Python script for fitting energy-vs-strain curves using the Birch-Murnaghan equation of state.
- PLOT-multiple-dir.py: Python visualization tool for multiple files with the same name in different directories.
- PLOT-optimized-geometry.py: Python visualization tool for relaxed coordinates of atoms in the unit cell.
- PLOT-loss-function.py: Python script for generating plots of the loss function.
From now on the symbol $ will indicate the shell prompt.
Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.
1. Basic background on graphene
Graphene is a two-dimensional crystal that is built out of carbon atoms arranged in a hexagonal structure (honeycomb lattice), as shown in the left panel of the following figure (taken from Pereira et al, EPL 92, 67001 (2010)).
This structure can be seen as a triangular lattice with a basis of two atoms per unit cell. The lattice vectors can be written as
(1)where $a$ = 2.46 Å is the experimental lattice constant. The reciprocal-lattice vectors are given by
(2)Of particular importance for the physics of graphene are the K points at the corners of the graphene Brillouin zone. These are named Dirac points for reasons that are related to a particular feature of the electronic band structure. This feature is illustrated in the next figure where a three-dimensional plot of the uppermost occupied and the lowermost empty bands are shown (taken from Castro Neto et al, Rev. Mod. Phys. 81, 109 (2009) ).
2. Using exciting for calculating properties of graphene
i) Two dimensional structures are treated as special three-dimensional periodic systems
In order to perform the calculation, we have to let exciting know that we like to treat a two-dimensional system. However, all calculations performed by exciting assume a three-dimensional periodicity along the directions defined by the lattice vectors. Therefore, when considering an isolated monolayer, one has to prepare a special three-dimensional crystal where the distances between a monolayer and its replica in the neighboring unit cells are large enough so that each monolayer behaves as it would be isolated.
The procedure one has to follow in this case, is illustrated in the next picture in a view parallel to the graphene monolayer.
In this case, a three dimensional crystal is created by replicating the monolayer along the out-of-plane axis at a given interlayer distance, corresponding to the amount of "vacuum" that should be added to the system. Its size has to be chosen such to prevent interlayer interactions.
ii) Preparation of the input file
The first step before starting an exciting calculation, is to create a directory for the example that you want to investigate. Thus, we create a directory /home/tutorials/graphene, and we move inside it.
$ mkdir /home/tutorials/graphene
$ cd /home/tutorials/graphene
The first example of an exciting input file (input.xml) for graphene is given in the following.
<input> <title>Graphene</title> <structure speciespath="$EXCITINGROOT/species"> <crystal scale="4.65"> <basevect> 0.5 0.8660254040 0.0 </basevect> <basevect> -0.5 0.8660254040 0.0 </basevect> <basevect> 0.0 0.0000000000 6.0 </basevect> </crystal> <species speciesfile="C.xml" rmt="1.20"> <atom coord="0.00000000 0.00000000 0.0"/> <atom coord="0.33333333 0.33333333 0.0"/> </species> </structure> <groundstate do="fromscratch" xctype="GGA_PBE" ngridk="8 8 1" rgkmax="4" swidth="0.01"> </groundstate> </input>
Input files for exciting are written in the XML format and are typically named input.xml. The XML format allows your data to be written in a structured way. Figuratively speaking, an exciting input is pretty much like an article with its sections and subsections. In case of XML data, sections and subsections are called elements.
Let us examine this example bit-by-bit. The first thing to be said is that an XML is not sensitive to line indentation. However, for the sake of clarity, line indentation is used in all examples of this tutorial to illustrate the hierarchy among elements. Now, let's come back to the input.xml shown above. The first and the last line indicate the beginning and the end of the input.
The element title contains some freely chosen text simply to describe the calculation. Keywords <title> and </title> indicate the beginning and the end of the element (note in the second keyword the presence and position of the slash /).
The next element, structure, describes the geometry and the chemical composition of a studied system. Note that the declaration of the structure section contains as additional information the parameter speciespath.
<structure speciespath="$EXCITINGROOT/species/">
Such parameters in the XML language are called attributes, and their values are always given in quotes regardless whether it is numerical, symbolic, or boolean information. In particular, the attribute speciespath defines the location, where the files with the data about chemical elements are stored.
Please, note that exciting is always using atomic units (Hartree, Bohr, etc.)! Thus, you can notice that in the example above the value of the experimental lattice constant of graphene has been considered (in units of bohr!).
... <crystal scale="4.65"> ...
Let us now investigate in more detail the other instructions given in the input file. The lattice vectors are chosen according to Eq. (1):
... <basevect> 0.5 0.8660254040 0.0 </basevect> <basevect> -0.5 0.8660254040 0.0 </basevect> <basevect> 0.0 0.0000000000 6.0 </basevect> ...
Note also that the amount of "vacuum" along the out-of-plane direction is set to 6 times the value of the in-plane lattice constant.
The position of the two atoms of carbon inside the unit cell of graphene is chosen according to first figure of this tutorial.
... <species speciesfile="C.xml" rmt="1.20"> <atom coord="0.00000000 0.00000000 0.0"/> <atom coord="0.33333333 0.33333333 0.0"/> </species> ...
In the last part of the input file, it appears a list of the parameters which are set for this calculation:
... <groundstate do="fromscratch" xctype="GGA_PBE" ngridk="8 8 1" rgkmax="4" swidth="0.01"> </groundstate> ...
In particular, according to this part:
- The GGA_PBE functional is chosen for the exchange-correlation energy.
- Brillouin-zone integration are performed by summing contributions from a 8$\times$8$\times$1 k-points grid. Notice that, due to the fact that one assumes that no interaction is present between the monolayer and its periodic replicas, the number of k-points along the out-of-plane direction is set to 1.
- The number of APW basis fuctions used in the calculation is set by assigning a value to the attribute rgkmax.
- A Gaussian smearing width of 0.01 Ha is chosen for icncluding the zero-gap metallic effects on the k-points grid.
iii) Execute calculations
To perform the first example calculation, copy the complete input listed above inside a file input.xml which you should create in the current directory.
N.B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT using the command
$ SETUP-excitingroot.sh
You can now run the script EXECUTE-single.sh in the directory where the input file input.xml has been created. If you want to save the results of this run in the subdirectory first-example you should run as in the following.
$ EXECUTE-single.sh first-example
===> Output directory is "first-example" <===
Running exciting for file input.xml -------------------------------------
Elapsed time = 0m13s
Run completed for file input.xml ----------------------------------------
$
This (preliminary) calculation solve self-consistently the Kohn-Sham equations and allows to determine the ground-state total enery of graphene in the configuration specified in the input file. More information can be found inside the file first-example/INFO.OUT.
iv) Calculate density of states
After you have completed the ground-state 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 — about how many electronic states there are at any given energy. To calculate it, move inside the directory first-example.
$ cd first-example
Then, you need to do the following simple modifications in the file input.xml inside this directory (for more details, see Input Reference):
- modify the value of the attribute do to "skip" to the element groundstate;
- add the element properties after the groundstate element;
- insert the subelement dos into the element properties;
- add the attribute nsmdos = "3" to the element dos.
The corresponding part of the input.xml should now look like this:
... <groundstate do="skip" ... </groundstate> <properties> <dos nsmdos="3"/> </properties> ...
All other instructions inside input.xml remain unchanged.
Now, execute exciting again by typing on the command line:
$ time excitingser
The run should take only a few seconds. Then, to visualize the DOS, execute
$ DOS Graphene_dos.agr
This produces the file Graphene_dos.agr which can be visualized using the program xmgrace. Open it with the command
$ xmgrace Graphene_dos.agr >/dev/null 2>&1 &
Question: Does the DOS reveal the presence of a gap? Why?
v) Calculate electronic 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 graphene, insert the subelement bandstructure in the element properties with the following specifications:
... <properties> <bandstructure> <plot1d> <path steps="100"> <point coord="0.0 0.0 0.0" label="GAMMA"/> <point coord="0.33333333 -0.33333333 0.0" label="K" /> <point coord="0.5 -0.5 0.0" label="M" /> <point coord="1.0 0.0 0.0" label="GAMMA"/> </path> </plot1d> </bandstructure> </properties> ...
Now execute excitingser again on the command line:
$ time excitingser
This makes the code produce the band structure, which is written to bandstructure.xml. To visualize the band structure, execute the following command:
$ BAND
Now, you have produced the xmgrace file Graphene_bandstructure.agr, which you can open, visualize and manipulate with xmgrace:
$ xmgrace Graphene_bandstructure.agr >/dev/null 2>&1 &
3. Full structure optimization
Up to here, all calculations have been performed considering graphene in its experimental equilibrium configuration. However, exciting is also able to predict the equilibrium configuration, i.e., the equlilibrium lattice constant of graphene. In this section, it will be shown how to proceed in order to find the equlilibrium lattice constant of two dimensional systems. We take as a starting point the input file used in the last section. We create a new working directory /home/tutorials/graphene/structure-optimization and we move inside it.
$ mkdir /home/tutorials/graphene/structure-optimization
$ cd /home/tutorials/graphene/structure-optimization
Copy inside this directory the input file you used in the first example above as input.xml.
$ cp /home/tutorials/graphene/input.xml ./input.xml
Modify the input file including the relax element just after groundstate (see Input Reference for the meaning of the different attributes).
... <relax method="harmonic" taunewton="1.0"/> </input>
The goal is to find the lattice constant which minimizes the total energy of the system. In our special case of a graphene monolayer, we should be careful in changing the lattice parameter in the monolayer without changing the distance between the monolayers along the out-of-plane direction. The effect of changing the amount of "vacuum" will be addressed in the next section.
A series of distorted structures corresponding to an in-plane homogeneous strain are generated and minimized with respect to the internal degrees of freedom (if applicable). To this purpose we use the script SETUP-elastic-strain.py which can be executed as follows
$ SETUP-elastic-strain.py optimization
Enter maximum Lagrangian strain [smax] >>>> 0.10
Enter the number of strain values in [-smax,smax] >>>> 21
------------------------------------------------------------------------
List of deformation codes for strains in Voigt notation
------------------------------------------------------------------------
0 => ( eta, eta, eta, 0, 0, 0) | volume strain
1 => ( eta, 0, 0, 0, 0, 0) | linear strain along x
2 => ( 0, eta, 0, 0, 0, 0) | linear strain along y
3 => ( 0, 0, eta, 0, 0, 0) | linear strain along z
4 => ( 0, 0, 0, eta, 0, 0) | yz shear strain
5 => ( 0, 0, 0, 0, eta, 0) | xz shear strain
6 => ( 0, 0, 0, 0, 0, eta) | xy shear strain
7 => ( 0, 0, 0, eta, eta, eta) | shear strain along (111)
8 => ( eta, eta, 0, 0, 0, 0) | xy in-plane strain
9 => ( eta, -eta, 0, 0, 0, 0) | xy in-plane shear strain
10 => ( eta, eta, eta, eta, eta, eta) | global strain
11 => ( eta, 0, 0, eta, 0, 0) | mixed strain
12 => ( eta, 0, 0, 0, eta, 0) | mixed strain
13 => ( eta, 0, 0, 0, 0, eta) | mixed strain
14 => ( eta, eta, 0, eta, 0, 0) | mixed strain
------------------------------------------------------------------------
Enter deformation code >>>> 8
$
In this example, (on screen) input entries are preceded by the symbol ">>>>". Entry values must be typed on the screen when requested. The first entry (in our example 0.10) represents the absolute value of the maximum strain for which we want to perform the calculation. The second entry (21) is the number of deformed structures equally spaced in strain, which are generated between the maximum negative strain and the maximum positive one. The third (last) entry (8) is a label indicating the type of deformation we have chosen: an in-plane homogeneous strain. Strain tensors are given in the Voigt notation.
After running the script, a directory called optimization is created, which contains input files for different strain values. To execute the series of calculation with input files created by SETUP-elastic-strain.py you have to run the script EXECUTE-elastic-strain.sh.
$ EXECUTE-elastic-strain.sh optimization
===> Output directory is "optimization" <===
Running exciting for file input-01.xml ----------------------------------
...
Run completed for file input-21.xml -------------------------------------
$
After the complete run, move to the working directory optimization.
$ cd optimization
Inside this directory, results of the calculation for the input file input-i.xml are contained in the subdirectory rundir-i where i is running from 01 to the total number of strain values. The data for energy-vs-strain curves are contained in the file energy-vs-strain.
For extracting the value of the theoretical equilibrium lattice parameter for the monolayer, run first the script SETUP-planar.py.
$ SETUP-planar.py
Then, use the script PLOT-birch.py as indicated below.
$ PLOT-birch.py
Input file is "energy-vs-strain".
Modified version for strained planar systems!
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A2 A3 lagrangian strain at minimum log(chi)
885.627 -10393.683 0.00011 ( alat0 = 4.6505 ) -4.07
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
$
As can be seen from this output an equilibrium lattice constant (alat0) of 4.6505 Bohr is obtained. The comparison with between this theoretical value and the experimental lattice constant of 4.65 Bohr is very good. However, the accuracy of this result should be checked for improved values of the numerical parameters of the calculation, as explained in the next section. Furthermore, together with the screen output the script PLOT-birch.py generates a PNG output file PLOT.png, which can be visualized using standard Linux tools.
The presence of internal relaxation can be verified by running
$ PLOT-optimized-geometry.py
This script also generate an output file PLOT.png, which can be visualized. If every plotted curve is flat, there is no internal relaxation.
Question: Do you find any relaxation of the relative position of the two atoms in the unit cell of the graphene monolayer? Why? Are there deformations for which relaxation of the internal degrees of freedom is important?
4. Convergence tests
In order to be able to rely on your calculation, you need to understand how the choice of the computational parameters change of the values of the physical quantities which are relevant for you. For most of the calculations performed in this tutorial, the most relevant computational parameters are listed below:
- The mesh of k-points (groundstate attribute: ngridk).
- The size of the basis set for expanding the wave function (groundstate attribute: rgkmax).
- The amount of "vacuum" considered in the out-of-plane direction between two adjacent monolayers.
Try to get a feeling on how the variation of these parameters changes the physical properties which are relevant in this tutorial.
5. Is graphene really planar?
In the rest of this tutorial, we are assuming that the planar configuration of graphene is the most stable one. In this section, we will try to say if this is always true. Create a new working directory /home/tutorials/graphene/planar and we move inside it.
$ mkdir /home/tutorials/graphene/planar
$ cd /home/tutorials/graphene/planar
In order to check if the planar configuration is the most favorable one, we perform a series of calculations in which one of the two atoms in the unit cell of the graphene monolayer is displaced along the out-of-plane direction. The calculation can be repeated for different values of the applied strain. The starting point is the following input file, which should be copied as input.xml inside /home/tutorials/graphene/planar.
<input> <title>Graphene</title> <structure speciespath="$EXCITINGROOT/species"> <crystal scale="4.6505"> <basevect> 0.5 0.8660254040 0.0 </basevect> <basevect> -0.5 0.8660254040 0.0 </basevect> <basevect> 0.0 0.0000000000 6.0 </basevect> </crystal> <species speciesfile="C.xml" rmt="1.03"> <atom coord="0.00000000 0.00000000 0.0"/> <atom coord="0.33333333 0.33333333 0.0"/> </species> </structure> <groundstate do="fromscratch" xctype="GGA_PBE" ngridk="8 8 1" rgkmax="4" swidth="0.01"> </groundstate> </input>
Notice that the value for the muffin-tin radius of the carbon atoms has ben decreased to avoid the overlap of muffin-tin spheres at large negative strains
... <species speciesfile="C.xml" rmt="1.03"> ...
As usual, use
$ SETUP-excitingroot.sh
to assign the correct value to the variable $EXCITINGROOT. Now you can produce a series of input files for each displacement using the script SETUP-graphene-along-c.py. As an example for an applied compression strain of -0.10, one has
$ SETUP-graphene-along-c.py str-0.10
Enter lagrangian strain >>>> -0.10
Enter maximum displacement umax [c/a] >>>> 0.2
Enter the number of displacements in [0,umax] >>>> 11
$
where one of the two carbon atoms is displaced along the out-of-plane direction up to a value of 20% of the length of the out-of-plane lattice vector.
Results are stored in the str-0.10 subdirectory. To esplicitly execute the calculation, you should run the script EXECUTE-diamond-phonon-g.sh:
$ EXECUTE-diamond-phonon-g.sh str-0.10
The values of the energy for the different displacements can be obtained using
$ PLOT-multiple-dir.py energy-vs-displacement str-0.10
and visualizing the corresponding output file PLOT.png.
The calculation can be repeated for different value of the applied strain. In particular, it is very intersting to investigate what happens at relatively high compression strains. Repeat the above procedure for the following list of applied strains
strain directory
-0.12 str-0.12
-0.14 str-0.14
-0.16 str-0.16
-0.18 str-0.18
-0.20 str-0.20
and look at the results using
$ PLOT-multiple-dir.py energy-vs-displacement str-0.10 str-0.12 str-0.14 str-0.16 str-0.18 str-0.20
Question: What do you conclude from the previous results? Is graphene really planar? Always?
6. Excited-states properties: The loss function
In this section, it will be shown how to calculate the loss function of graphene using time-dependent density-functional theory (TDDFT) within the random-phase approximation (RPA). The loss function is the quantity which is directly measurable using electron-energy loss spectroscopy (EELS).
Create a new working directory /home/tutorials/graphene/loss-function and we move inside it.
$ mkdir /home/tutorials/graphene/loss-function
$ cd /home/tutorials/graphene/loss-function
Copy inside the current directory the input.xml given in Section 2 (do not forget to update the value of the attribute speciespath). In order to calculate excited-states properties, make the following changes in input.xml.
- Modify the inplane lattice vector to the optimized value of 4.6505 Bohr.
- Include the element xs just after groundstate as in the next example (see Input Reference and Excited states from TDDFT for more information about the meaning of the different attributes).
... <xs xstype ="TDDFT" rgkmax="4" ngridk="8 8 1" vkloff="0.097 0.273 0.493" nempty="80" gqmax="1.0" broad="0.02" tevout="true"> <energywindow intv="0.0 1.0" points="500" /> <tddft fxctype="RPA"/> <qpointset> <qpoint> 0.0 0.0 0.0 </qpoint> </qpointset> </xs> ...
Then, execute exciting using the script EXECUTE-single.sh:
$ EXECUTE-single.sh test-1
Wher the run is completed, move to the directory test-1 which includes the result of the calculation.
$ cd test-1
Now, the xx component of the loss function can be visualized with the help of the script PLOT-loss-function.py.
$ PLOT-loss-function.py LOSS_FXCRPA_OC11_QMT001.OUT.xml
As usual, the visualization script generates the output file PLOT.png. You can compare your result with the one displayed in Figure 3 of the paper Trevisanutto et al, Phys. Rev. Lett.101, 226405 (2008).
Question: Which are the main differences between your result and the one of Trevisanutto et al.? How can you improve the accuracy of your calculations? Which are the most important parameters to be optimized?