Physical Properties of Graphene

by Ronaldo Rodrigues Pela, Pasquale Pavone, Keith Gilmore, & Caterina Cocchi for exciting oxygen

Purpose: In this computational experiment, you will:

  • learn how to set up and execute exciting calculations for graphene;
  • calculate the lattice parameter of graphene and estimate the precision and accuracy of the calculation;
  • generate the loss function of graphene and comment on the precision and accuracy of your results.

Before starting the experiment, please read carefully the tutorial below.

All the experiments will be done remotely on the computers of the PC-Pool that belongs to the Physics Department of the Humboldt-Universität zu Berlin. You need a "Physik" account to have access to the PC-Pool.

If you have any problems, please write an email to this ed.nilreb-uh.kisyhp|eromligk#sserdda.

Practical details about the report, the grades, and the zoom link for the experiment are provided in this file.

1. Remote access to the PC-Pool

2. Preliminary notes

Read the following before starting with the laboratory exercise!

Requirements: Bash shell, Python numpy, lxml, matplotlib.pyplot and sys libraries.

Define relevant environment variables

Before starting, we need to set a few environment variables. We first define the root directory for the exciting code, called EXCITINGROOT, and then use it to define the path EXCITINGTOOLS that contains various plotting tools and other scripts that we will use. Issue the commands

export EXCITINGROOT=/tmp/exciting
export TIMEFORMAT="   Elapsed time = %0lR"
export WRITEMINMAX="1"

Since the exciting has been compiled with an Intel compiler, you need to load the specific libraries:

source /usr/global/intel-2019/bin/ intel64
Text editor

In this laboratory description, we assume that you are using nano as a text editor. Although it is quite intuitive, if you need a tutorial, please refer to this site or this one. Other possible choices are: vi, kwrite, and nedit. Independent of your preference, keep in mind that you need to use this text editor remotely when connected to the computers of the PC-Pool. All the editors mentioned are installed there.

Programs for making plots

As we will see below, exciting has built-in tools for making plots. However, you may choose to use alternate plotting tools. Gnuplot is a very simple option. Matplotlib is a more elaborate choice and requires some knowledge of Python Programming. If you want something more intuitive, you may prefer Scidavis, or Qtiplot (recommended by the PC-Pool - see here). In any case, you will need to install the program of your choice on your own computer to make final plots for your report.

3. Basic background on graphene

Graphene is a two-dimensional crystal that is built from carbon atoms arranged in a hexagonal structure (honeycomb lattice), as shown in the left panel of the following figure.


This structure can be seen as a triangular lattice with a basis of two atoms per unit cell. The lattice vectors are

\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}

where $a$ = 2.46 Ã… is the experimental lattice constant. The reciprocal-lattice vectors are given by

\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 points K and K' at the corners of the Brillouin zone. These are named Dirac points because the bands disperse linearly near the Fermi level at these locations. This feature is illustrated in the next figure where a three-dimensional plot of the highest occupied (blue) and the lowest unoccupied (red) bands are shown for the honeycomb lattice of graphene. The inset gives an expanded view of the so-called Dirac cone at the K point.

Special treatment of two dimensional structures

Although graphene is a two-dimensional material, exciting always employs three-dimensional periodicity along the directions defined by the lattice vectors. Therefore, we need to prepare an artificial three-dimensional crystal where the distance between a monolayer and its nearest periodic replica is large enough that each graphene monolayer is effectively isolated.

The procedure is illustrated in the next picture with a view parallel to the graphene monolayer.

A three-dimensional crystal is created by replicating the monolayer along the out-of-plane axis at the given interlayer distance, which corresponds to the amount of "vacuum" which should be added between graphene sheets. The interlayer distance must be large enough so that interlayer interactions are minimal.

4. First calculation

i) Creating your home directory

The very first step after you login is to create your personal directory. We advise you to do this in the following way:

  • First, define a global variable called myid which we will use throughout this tutorial. You can do this with the command below, by replacing the string "SubstituteThisWithSomethingReasonable" with something you choose, e.g. your name or username. It should be a single word without spaces and special characters.
export myid=SubstituteThisWithSomethingReasonable
  • Second, create the directory and move into it:
mkdir -p /home/$myid
cd /home/$myid
ii) Preparation of the input file

Inside the folder /home/$myid, create the directory 01-FirstRun for the first calculation and move inside it.

mkdir -p 01-FirstRun
cd 01-FirstRun

Our first example of an exciting input file (input.xml) for graphene is given below. To create the input file, just type

nano input.xml

Then, copy the content below and paste it in your screen (it is advisable to do this by clicking with the right button of your mouse, rather than using Ctrl+C and Ctrl+V).

   <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>
      <species speciesfile="C.xml" rmt="1.20">
         <atom coord="0.00000000  0.00000000  0.0"/>
         <atom coord="0.33333333  0.33333333  0.0"/>
      ngridk="8 8 1"

when finished, just type Ctrl+O to save (you may be asked to confirm the name of the file by pressing Enter), and Ctrl+X to go back to the command line.

Important: All input parameters are given in atomic units! Therefore, we specify the lattice constant as 4.65 Bohr rather than 2.46 Angstrom!

Important: Before running the calculation it is necessary to replace the string "$EXCITINGROOT" within the speciespath in input.xml by the actual value of the environment variable $EXCITINGROOT. You can do this by using the following command:
ii) Executing the calculation

To run exciting you can invoke the executable as exciting_smp :

time exciting_smp

If everything runs fine, you should get output like

   Elapsed time = 0m52s

This initial calculation self-consistently solves the Kohn-Sham equations to determine the ground-state electron density and total energy of graphene for the atomic configuration specified in the input file. More information can be found inside the file INFO.OUT. If you see a file WARNING.OUT, this means there was a problem during the calculation. Check the contents of this file for more information.

5. Lattice parameter

In this section, you will learn how to analyze the precision of your calculations.

Computational implementations of density functional theory must make certain numerical approximations. These often relate to truncating infinite summations to a finite number of terms, or discretizing integrations on a finite grid. This introduces corresponding numerical parameters that can be considered convergence parameters. Convergence parameters are not free or empirical parameters. Rather, these convergence parameters determine the quality of our calculations: when we include more terms of a series, or sample an integral more densely, we will obtain a more converged result. Here, we will consider two such convergence parameters in particular:

  • The mesh of k-points (groundstate attribute: ngridk) for performing integrals over the Brillouin zone.
  • The size of the basis set for expanding the wave function (groundstate attribute: rgkmax).

A discrete mesh of k-points replaces a continuous integral over k-space by a sum over the selected k-points. The parameter ngridk sets this k-point grid.
The parameter rgkmax controls the size of the basis in terms of the largest wavevector of the augmented plane waves.
In principle, with a denser set of k-points and a larger value of rgkmax the calculation will be more precise, but also more costly. Therefore, there is a trade-off between the computational cost and the desired precision.

i) Total Energy

We first investigate how ngridk and rgkmax impact the total energy. Create a new directory for this test and move inside it.

cd /home/$myid/
mkdir -p 02-TotalEnergy
cd 02-TotalEnergy

Inside this new directory, create two folders: one to store the tests for ngridk and the other for rgkmax.

mkdir -p rgkmax ngridk

We start by setting a fixed number of k-points and varying rgkmax, and then we inspect how the total energy changes. To do so, move inside the folder rgkmax and copy the necessary script:

cd rgkmax

The file automates the task of varying rgkmax and running the corresponding calculations. This file contains the values of rgkmax that will be tested. You can open the file


and look for the line
for rgkmax in 4.0 4.5 5.0 5.5 ; do

In this line, you can replace the default sequence by any other of your choice (be aware that a calculation with rgkmax of 7.0 might take half an hour). Note that the point '.' is the decimal separator. After editing the file, save it and exit with Ctrl+O, Ctrl+X.

When you are ready to run the calculations, just execute


For the default sequence of rgkmax values, the output will look like this:

Now running rgkmax =  4.0
Calculation for rgkmax =  4.0  finished!
Now running rgkmax =  4.5
Calculation for rgkmax =  4.5  finished!
Now running rgkmax =  5.0
Calculation for rgkmax =  5.0  finished!
Now running rgkmax =  5.5
Calculation for rgkmax =  5.5  finished!

Now we analyze the results. Create a file with the values of rgkmax and the corresponding total energy, which is printed as the last value in the file TOTENERGY.OUT. For example, for rgkmax equal to 4.0, you can see the corresponding total energy with the command:

tail -2 4.0/TOTENERGY.OUT

You will get something like

The total energy in this case is -76.1130433503 Ha (remember that exciting uses the Hartree as the unit of energy). However, a question here would be, how many digits are meaningful? To answer it, we need to investigate how the total energy changes when we change a certain parameter (e.g. rgkmax).

Make a plot of the total energy versus rgkmax (you may wish to plot the energy in eV rather than Ha, look up the conversion).

Questions: What can you conclude about the convergence of the total energy with respect to rgkmax? How many significant digits can you report for the total energy?

We can now make a similar analysis for the impact of ngridk on the total energy. To do so, choose a reasonable value for rgkmax that will be kept constant during the calculations. If you are not sure, ask for help.

Switch to the other directory you previously made for testing convergence of the total energy with respect to the k-point grid and copy the script for automating the calculations:

cd /home/$myid/02-TotalEnergy/ngridk

Edit the script to specify the value of rgkmax


Look for the lines
for ngridk in 3 4 5 6 7; do

The default value for rgkmax, 1.0, is far from good convergence and you must replace it with a more reasonable value. Be careful not to use blank spaces in this line or the script will crash (which means it should be like "rgkmax(no blank space)=(no blank space)ValueWithPointAsDecimalSeparator"). Be aware that a calculation with rgkmax of 7.0 might take half an hour.
Substitute the default sequence of ngridk with one of your choice. After editing the file, save it and exit with Ctrl+O, Ctrl+X. You may then execute the script with

Remember, you can check the total energy for a specific ngridk by inspecting the last value of TOTENERGY.OUT - for example, for ngridk equal to 6:


Now, as you did before, make a file and plot of your results for the total energy versus ngridk.

Questions: What do you conclude about the convergence of the total energy with respect to ngridk? How many significant digits can you report for the total energy?

ii) Lattice parameter

Now we will predict the lattice parameter of graphene. To prepare the calculation, execute

mkdir -p /home/$myid/03-Lattice
cd /home/$myid/03-Lattice

This should create a new directory, move inside it, and copy the script

Now, open the file with nano


In the 2nd and 3rd lines, you will find


These values make no sense. Change them to reasonable ones based on the convergence tests you just performed, then save your changes and go back to the command line.
To run the script, simply execute:


Your output will look like

Running exciting for file input-01.xml ----------------------------------
Run completed for file input-01.xml -------------------------------------

Running exciting for file input-02.xml ----------------------------------
Run completed for file input-02.xml -------------------------------------

Running exciting for file input-03.xml ----------------------------------
Run completed for file input-03.xml -------------------------------------

Running exciting for file input-04.xml ----------------------------------
Run completed for file input-04.xml -------------------------------------

Running exciting for file input-05.xml ----------------------------------
Run completed for file input-05.xml -------------------------------------

Running exciting for file input-06.xml ----------------------------------
Run completed for file input-06.xml -------------------------------------

Running exciting for file input-07.xml ----------------------------------
Run completed for file input-07.xml -------------------------------------

Running exciting for file input-08.xml ----------------------------------
Run completed for file input-08.xml -------------------------------------

Running exciting for file input-09.xml ----------------------------------
Run completed for file input-09.xml -------------------------------------

Running exciting for file input-10.xml ----------------------------------
Run completed for file input-10.xml -------------------------------------

Running exciting for file input-11.xml ----------------------------------
Run completed for file input-11.xml -------------------------------------

alat0 [Bohr]    log(chi)
    4.6856      -3.94
This script tests 11 values of the lattice parameter, and for each one, runs a calculation and stores the total energy. In the end, it makes a fit of the points to the Birch-Murnaghan equation, obtaining the lattice parameter which minimizes the total energy. This result is printed on the last line, and represents the equilibrium lattice parameter found for the selected values of ngridk and rgkmax. A plot of the energy vs. the unit cell volume together with the fit is given in the following Figure.

Your task now is to

  1. Fix ngridk and vary rgkmax, checking how this impacts the lattice parameter.
  2. Fix rgkmax and check the dependence of the lattice parameter on ngridk.

For each value of rgkmax and ngridk, you need to edit the file, execute it, and record the resulting lattice parameter.

Questions: What is your theoretical prediction of the lattice parameter? What comment can you make on the precision which you have been able to achieve? Comparing to the experimental lattice parameter, what comment can you make on the accuracy of your calculation?

6. Loss function

In this section, we will learn 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 for these calculations and move inside it.

mkdir -p /home/$myid/04-LossFunction
cd /home/$myid/04-LossFunction

Create the file input.xml with nano.

nano input.xml

Copy and paste the example input file below into your input.xml, then edit the contents to:
  • Modify the lattice parameter to one that you obtained previously.
  • Replace the question marks within rgkmax and ngridk with your choices.

See Input Reference and Excited states from TDDFT for more information about the meaning of the attributes within the xs element.

   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="YOUR_LATTICE_PARAMETER">
         <basevect> -0.5  0.8660254040  0.0 </basevect>
         <basevect>  0.5  0.8660254040  0.0 </basevect>
         <basevect>  0.0  0.0000000000  6.0 </basevect>
      <species speciesfile="C.xml" rmt="1.20">
         <atom coord="0.00000000  0.00000000  0.0"/>
         <atom coord="0.33333333  0.33333333  0.0"/>
      ngridk="? ? 1"
      ngridk="? ? 1"
      vkloff="0.097 0.273 0.493"
         intv="0.0 1.0"
         points="500" />
          <qpoint> 0.0 0.0 0.0 </qpoint>

Save this file and exit. We also must set the species path by running:

Before running the calculation, it is advisable to create a subdirectory that identifies the calculation according to the choice of convergence parameters (then enter this subdirectory). In this examples, we call it nk1_rgk1.0_gq1.0_nempty80. Please adjust the numbers (i.e. in "nk1" and "rgk1.0") to reflect your choices of ngridk and rgkmax. The remaining two numbers correspond to values of the parameters gqmax ("gq1.0") and nempty ("nempty80") that appear within the xs element.

mkdir -p nk1_rgk1.0_gq1.0_nempty80
cp input.xml nk1_rgk1.0_gq1.0_nempty80/
cd nk1_rgk1.0_gq1.0_nempty80/

To run the calculation, execute exciting the command:
time exciting_smp

Depending on your choices, the calculation can take some minutes. When it is completed, you can inspect the files LOSS_FXCRPA_OC11_QMT001.OUT and LOSS_FXCRPA_OC33_QMT001.OUT, which have the loss function for a polarization parallel and perpendicular to the graphene plane, respectively. Copy these files to your own computer and plot the loss function for each polarization. One example of a convergence test with respect to rgkmax is given in the following picture:


Now investigate how sensitive the results are to the following parameters:

  • rgkmax
  • ngridk
  • nempty
  • gqmax

To better understand the meaning of the last two parameters, please refer to this page.

For each set of parameters, you need to repeat the following procedures

cd /home/$myid/04-LossFunction
nano input.xml
mkdir -p nk1_rgk1.0_gq1.0_nempty80/
cp input.xml nk1_rgk1.0_gq1.0_nempty80/
cd nk1_rgk1.0_gq1.0_nempty80/
time /home/exciting_smp

When editing the input file, update rgkmax, ngridk, nempty and gqmax. You should vary one of them at a time and while keeping the others fixed. Adapt the name nk1_rgk1.0_gq1.0_nempty80/ to reflect your current choice.

Questions: What comment can you make about the precision of the energy positions of the main peaks in your calculation of the loss function? What about their accuracy? You can compare your calculated loss function with the experimental one presented in this work by J.C. Idrobo and W. Zhou (see Fig. 1 on page 14). What would be the effect of changing the amount of vacuum in the unit cell on your calculations of the loss function?

7. Report & Grading

You will be assigned a grade of up to 20 points for this lab. Up to 10 points will be given based on how well you execute the tasks of this lab today. For this, you are required to submit by email within 1 week the plots that you generate in sections

  • 5(i), convergence of the total energy with respect to rgkmax and with respect to ngridk. State the values of rgkmax and ngridk that you selected for use in section 5(ii); provide only the significant digits.
  • 5(ii), dependence of the lattice parameter on rgkmax and ngridk. Present plots of the dependence of the lattice parameter on these two values. State your theoretical prediction for the lattice constant, giving only the significant digits.
  • 6, dependence of the Loss function on rgkmax, ngridk, nempty and gqmax. Show plots of the Loss function as it depends on each of these four parameters.

Up to another 10 points will be given based on your lab report, which is due within two weeks of the day of your experiment. For more detailed instructions about preparing your report see the guidelines. Templates for the report are available in Word format or TeX format.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License