Graphene: from the Ground State to Excitations

by Pasquale Pavone, Caterina Cocchi, Wahib Aggoune, & Ronaldo Rodrigues Pela for exciting oxygen

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 loss function of is presented as an example for the calculation of excited-states properties.

### 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-files.py: Python visualization tool for multiple files with the same name in different directories (for more details, see The python script "PLOT-files.py").
• PLOT-optimized-geometry.py: Python visualization tool for relaxed coordinates of atoms in the unit cell.
• PLOT-band-structure.py: Python visualization script for plotting and comparing energy bands(for more details, see The python script "PLOT-band-structure.py") .
• PLOT-dos.py: Python visualization script for plotting and comparing density of states (for more details, see The python script "PLOT-dos.py").

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 material composed by carbon atoms in a hexagonal structure (honeycomb lattice), as shown in the left panel of the following figure. This structure can be described as an hexagonal lattice with a basis of two atoms. One choice for the lattice vectors is (1) \begin{align} {\bf a}_1= a \left(-\frac{1}{2}, \frac{\sqrt{3}}{2}\right) \qquad {\rm and} \qquad {\bf a}_2= a \left(\frac{1}{2}, \frac{\sqrt{3}}{2}\right)\;, \end{align} wherea= 2.46 Å is the experimental lattice constant. The reciprocal-lattice vectors are (2) \begin{align} {\bf b}_1= \frac{2\pi}{a} \left(-1, \frac{\sqrt{3}}{3}\right) \qquad {\rm and} \qquad {\bf b}_2= \frac{2\pi}{a} \left(1, \frac{\sqrt{3}}{3}\right)\;. \end{align} Of particular importance for the physics of graphene are the K points at the corners of the first Brillouin zone. These are named Dirac points due to the particular feature of the electronic band structure, illustrated in the next figure with a three-dimensional plot of the uppermost occupied (blue) and the lowermost empty (red) bands. ### 2. Using exciting for calculating properties of graphene ##### 2.1 Two dimensional structures are treated as special three-dimensional periodic systems Since exciting assumes three-dimensional periodicity along the directions defined by the lattice vectors, to simulate the two-dimensional system, some care needs to be taken. The solution is to consider a special three-dimensional crystal where the distances between a monolayer and its replica are large enough so that each monolayer behaves as it would be isolated. This illustrated in the next figure. As seen, a three dimensional structure is created by replicating the monolayer along the out-of-plane axis at a given interlayer distance, corresponding to the amount of "vacuum" added to the system. Its size has to be selected to prevent interlayer interactions. ##### 2.2 Preparation of the input file Initially, create a directory for this tutorial, graphene, and enter it.  mkdir -p graphene
$cd graphene  The 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>


##### 2.3 Execute calculations

To perform the first calculation for graphene, create the input.xml (as provided above) inside the current directory (graphene/).

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 current directory (where input.xml is). To save the results in the subdirectory first-example: $ EXECUTE-single.sh first-example

===> Output directory is "first-example" <===

Running exciting for file input.xml -------------------------------------

Elapsed time = 0m16s

Run completed for file input.xml ----------------------------------------

$  This calculation solves self-consistently the Kohn-Sham equations and allows to determine the ground-state total energy. For more information, check the file first-example/INFO.OUT. ##### 2.4 Density of states After the ground-state is completed, it is time to investigate more properties. One of the most fundamental ones is the density of states (DOS): it provides information on how many electronic states there are at a given energy. To calculate it, move inside the directory first-example. $ cd first-example/


Then, modify the file input.xml inside this directory (for more details, see Input Reference) as described below:

1. in the element groundstate, change the attribute do to "skip";
2. add the element properties after groundstate;
3. insert the sub-element dos into the element properties;
4. add the attribute nsmdos = "3" to the element dos.

The input.xml file should look like this:

...
<groundstate
do="skip"
...
</groundstate>

<properties>

<dos nsmdos="3"/>

</properties>
...


Now, execute exciting with:

$time exciting_smp  This should take only a few seconds. Then, to visualize the DOS, execute $ PLOT-dos.py


The script PLOT-dos.py (discussed in details in The python script "PLOT-dos.py") produces the PNG file PLOT.png. You can visualize this file with standard tools.

Questions:

• Does the DOS reveal the presence of a gap? Why?
• What would you expect due to the existence of the Dirac cones?
• Is the K point included in the ngridk mesh?
• What does it change if you use ngridk = "9 9 1"? Why?

##### 2.5 Electronic band-structure

The electronic band-structure details the dependence of the energy eigenvalues on the wavevectors k. To obtain it, proceed as follows:

• Comment or delete the sub-element dos inside properties.
      <!--dos nsmdos="3"/-->

• Insert the sub-element bandstructure inside 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 exciting with the command:

$time exciting_smp  To visualize the electronic band-structure (which is written inside the file BAND.OUT), you can use the script PLOT-band-structure.py (discussed in details in The python script "PLOT-band-structure.py"). $ PLOT-band-structure.py


This script produces the PNG file PLOT.png which can be visualized with standard tools.

### 3. Full structure optimization

Up to here, all calculations assumed graphene in its experimental equilibrium configuration. However, exciting is also able to predict this equilibrium configuration, i.e., the equlilibrium lattice constant of graphene. This is addressed in this section.

Move to previous directory, create a new folder structure-optimization, and move inside it.

$cd ..$ mkdir -p structure-optimization/
$cd structure-optimization/  We take as a starting point the input file used in the last section. Copy it using: $ cp ../input.xml ./


Modify the input file including the relax element just after groundstate (see Input Reference for the meaning of the different attributes).

...
</groundstate>

<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 must be generated. The script SETUP-elastic-strain.py serves to this purpose and 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 ">>>>" and must be typed on the screen when requested. The first entry (in our example, 0.10) represents the absolute value of the maximum strain considered in the calculations. The second entry (21) is the number of deformed structures equally spaced in strain generated between the maximum negative strain and the maximum positive one. The third entry (8) is a label indicating the type of deformation: an in-plane homogeneous strain, in our example. 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 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, 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.566 -10393.618 0.00010 ( 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, 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 PLOT.png, overwritting the existing one. If every plotted curve is flat, there is no internal relaxation. Questions: • Do you find any relaxation of the relative position of the two atoms in the unit cell of graphene? 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 the values of the physical quantities relevant for you. For most of the calculations performed in this tutorial, the most important 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. Now, perform some calculation to verify how changes in these parameters modify the physical properties relevant in this tutorial. ### 5. Is graphene really planar? In the rest of this tutorial, we assumed that the planar configuration of graphene is the most stable one. In this section, we can check if this is always true. Create a new folde planar inside the directory graphene and we move into it. $ cd ..
$mkdir -p planar$ cd planar


To analyze if the planar configuration is the most favorable one, we perform a series of calculations displacing one of the two atoms in the unit cell 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, to be copied as input.xml inside the current folder 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 muffin-tin radius of carbon has been decreased to avoid the overlap of muffin-tin spheres at large negative strains ... <species speciesfile="C.xml" rmt="1.03"> ...  As usual, employ $ 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 compressive 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. The generated input files for each displacement are stored in the str-0.10 subdirectory. To carry out the calculations, run the script EXECUTE-diamond-phonon-g.sh: $ EXECUTE-diamond-phonon-g.sh str-0.10


The calculated energy for the different displacements can be obtained using

$PLOT-files.py -f energy-vs-displacement -d str-0.10 -pm -ly 'Energy [Ha]' -lx '% of out-of-plane displacement [c/a]'  Inspect now the corresponding output file PLOT.png. The same calculation can be repeated for different applied strains. It is specially interesting to investigate what happens at relatively high compression strains. Repeat the above procedure for the following 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-files.py -f energy-vs-displacement  -d str-0.10 str-0.12 str-0.14 str-0.16 str-0.18 str-0.20  -pm  -ly 'Energy [Ha]'  -lx '% of out-of-plane displacement [c/a]'


Questions:

• At a given compressive strain, which percentage of out-of-plane displacement minimizes the total energy?
• What do you conclude from the previous results?
• Is graphene really planar? Always?

### 6. Excited-states properties: The loss function

In this section, we calculate the loss function of graphene using time-dependent density-functional theory (TDDFT) within the random-phase approximation (RPA). The loss function is a quantity directly measurable using electron-energy loss spectroscopy (EELS).

Go to the directory graphene, create a new folder loss-function, and enter it.

$cd ../../$ mkdir -p loss-function
$cd loss-function  Copy the input.xml given in Section 2 (if this is the case, update the attribute speciespath). You can do this with $ cp ../input.xml ./


To calculate excited-states properties, make the following changes in input.xml.

• Modify the in-plane lattice vector to the optimized value of 4.6505 Bohr.
• Change inside groundstate the value of the attribute ngridk to "9 9 1".
• Include the element xs just after groundstate as given below (see Input Reference and Excited states from TDDFT for more information about the meaning of the different attributes).
...
</groundstate>

<xs
xstype ="TDDFT"
rgkmax="4"
ngridk="9 9 1"
vkloff="0.097 0.273 0.493"
nempty="80"
gqmax="1.0"
tevout="true">

<energywindow
intv="0.0 1.0"
points="500" />

<tddft
fxctype="RPA"/>

<qpointset>
<qpoint> 0.0 0.0 0.0 </qpoint>
</qpointset>

</xs>

</input>


Then, execute exciting using the script EXECUTE-single.sh:

$EXECUTE-single.sh rpa-k09  When the calculation finishes, the results are stored in the directory rpa-k09 where the suffix k09 is representative for the choice ngridk = "9 9 1" inside both the groundstate and xs elements. Now, the xx component of the loss function can be copied to the current directory as follows $ cp rpa-k09/LOSS_FXCRPA_OC11_QMT001.OUT loss-k09


You can visualize the loss function employing the script PLOT-files.py (for a detailed description of the script arguments see The python script "PLOT-dos.py") as follows:

$PLOT-files.py -f loss-k09 -rc -lx 'Energy [eV]' -ly 'Loss function'  This generates the file PLOT.png, which should look as follows. ##### 6.1 Converging the loss spectrum w.r.t. the mesh in reciprocal space The results shown above are obtainedd using ngridk = "9 9 1" for both the ground-state (groundstate) and excited-states (xs) calculations. The loss spectrum is not fully converged at this mesh. Further calculations should be performed with finer meshes until the desired convergence is reached. In order to do so, repeat the calculations above for "18 18 1" --> rpa-k18 --> loss-k18 "27 27 1" --> rpa-k27 --> loss-k27 "36 36 1" --> rpa-k36 --> loss-k36 "45 45 1" --> rpa-k45 --> loss-k45 ...  Notice that at each step the complexity of the calculation is increased in a non linear way, resulting soon in quite large execution times. Hint: You can compare two or more loss spectra by using the script PLOT-files.py similarly to what follows $ PLOT-files.py -f loss-k09 loss-k18 loss-k27 -rc  -lx 'Energy [eV]'  -ly 'Loss function'


Questions:

• At which value of ngridk can you consider the loss spectrum as converged?
• Consider separately the features in the loss spectrum below and above 10 eV. Do they converge at the same value of ngridk?

##### 6.2 Converging the loss spectrum w.r.t. the number of basis functions

Another very inportant parameter for the convergence of the loss spectrum is the attribute rgkmax of the element groundstate. This parameter, together with the total number of local orbitals, is related to the number of basis function used for the expansion of the Kohn-Sham wave-functions.

The loss spectrum is not fully converged for rgkmax = "4.0", which is the value we used above. Further calculations should be performed with larger values (e.g., "4.5", "5.0", "5.5", etc.) until the desired convergence is reached.

Questions:

• How can you further improve the precision of your calculations?
• Which further parameters are relevant to be optimized?