Physical Properties of Graphene

by Ronaldo Rodrigues Pela, Pasquale Pavone, & Caterina Cocchi for exciting nitrogen

Purpose: In this computational experiment, you will:

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

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

All the experiments will be done remotely in 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 problem, please write an email to this ed.nilreb-uh.kisyhp|odlanor#sserdda.

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

1. Preliminary notes

Read the following paragraphs before starting with the rest of this tutorial!

Text editor

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

Programs for making plots

There are many possibilities. We recommend matplotlib, although it requires some knowledge of Python Programming. If you want something easier and more straightforward, you may fell more comfortable with Scidavis, or Qtiplot (recommended by the PC-Pool - check it here). In any case, you will need to install the programm of your choice in your own computer.

2. Basic background on graphene

Graphene is a two-dimensional crystal, which is built out of carbon atoms arranged in 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

\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 at the corners of the graphene Brillouin zone . These are named Dirac points for reasons which are related to a particular featrure of the electronic band-structure of graphene. This feature is illustrated in the next figure where a three.dimensional plot of the uppermost occupied and the lowermost empty bands are shown for the honeycomb lattice of graphene (taken from Castro Neto et al, Rev. Mod. Phys. 81, 109 (2009) ).

Special treatment of two dimensional structures

In order to perform the calculation, we have to let exciting know that we want to consider a two-dimensional system. However, exciting always assumes a three-dimensional periodicity along the directions defined by the lattice vectors. Therefore, we need to prepare a special three-dimensional crystal where the distances between a monolayer and each of its replicas are large enough so that each monolayer behaves as isolated.

The procedure 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, which corresponds to the amount of "vacuum" which should be added to the system. This vacuum has to be chosen in such a way to prevent interlayer interactions.

3. Remote access to the PC-Pool

4. Creating your directory

The very first step after you login is to create a directory which will employed to carry out our calculations. We advise to perform it the following way:

  • First, we define a global variable named as myid which will be used along this tutorial. You can do this with the command below, by replacing the string "SubstituteThisWithSomethingReasonable" with something you choose, e.g. your name or a nickname. It should be a single word without spaces and special characters.
export myid=SubstituteThisWithSomethingReasonable
  • Second, we create a directory and enter it:
mkdir -p /home/$myid
cd /home/$myid

5. First calculation

i) Preparation of the input file

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

mkdir -p /home/$myid/01-test
cd /home/$myid/01-test

The first example of an exciting input file (input.xml) for graphene is given in the following. 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="./">
      <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.

ii) Execute calculations

We still need a species file to run the calculation. You can copy from here, or better, just type in the current directory:


Before running exciting, it is useful to better specify the time format:

export TIMEFORMAT="   Elapsed time = %0lR"

And since the provided executable has been compiled with an Intel compiler, you need to load the specific libraries:

source /usr/global/intel-2019/bin/ intel64> /dev/null

Now you can run exciting. The code has already been compiled and you could invoke the executable as /home/excitingser as below:

time /home/excitingser

If everything runs fine, you should get an output like

   Elapsed time = 0m28s

This (preliminary) calculation solve self-consistently the Kohn-Sham equations and allows to determine the ground-state total energy of graphene in the configuration specified in the input file. More information can be found inside the file INFO.OUT.

6. Lattice parameter

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

To have a numerical feasible implementation, some approximations, specially truncations, are considered. Their quality can be tuned by the choice of the computational parameters, and here we select:

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

A discrete mesh of k-points replaces an integral over the k-space by a sum over the selected k-points. And rgkmax controls the size of the basis in terms of the largest wavevector of the augmented plane waves. In principle, the more k-points and the larger rgkmax, the more precise and more costly a calculation. There is therefore a trade-off between the computational cost and the desired precision.

i) Total Energy

We check initially how ngridk and rgkmax impact the total energy. We create a new directory for this test and move inside it.

mkdir -p /home/$myid/02-test-Energy
cd /home/$myid/02-test-Energy

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

mkdir -p rgkmax ngridk

We start setting a fixed number of k-points, and varying rgkmax, and then we inspect how the total energy changes. To do so, we move inside the folder rgkmax, and download the script; and then we make it executable:

cd rgkmax
chmod +x

The file automatizes the task of varying rgkmax and running a corresponding calculation. You will find in this file some values of rgkmax that can 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, as already pointed out before.

When you are ready to run the calculations, just execute


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

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, you should analyze the results. You can make a table 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 a rgkmax equal to 4.0, you can see the corresponding total energy with the command:

nano -v 4.0/TOTENERGY.OUT

the argument "-v" passed to nano ensures that we open the file in a read-only mode. You will get as output something like

To escape the file back to the command line, make use of Ctrl+X. From the last line, we conclude that the total energy in this case is -76.1130433503 Ha. 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).

Now, you can make a table and/or a plot of rgkmax and the total energy.

Questions: What do you conclude about the convergence of this property (total energy) with respect to rgkmax? Is the convergence fast or slow? How reliable are the digits of 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 along the calculations. If you are not sure, ask for help. Here are some

cd /home/$myid/02-test-Energy/ngridk
chmod +x

Now, it is time to edit the script and adapt it to your taste.


Look for the lines
for ngridk in 3 4 5 6 7 8 9 10 11 12; do

The default value for rgkmax, 1.0, makes no sense, and you must replace it by a reasonable one. Be careful not to let blank spaces in this line, otherwise the script will crash (which means, it should be like "rgkmax(no blank space)=(no blank space)ValueWithPointAsDecimalSeparator").
Substitute the default sequence of ngridk with another one of your choice (be aware that a calculation with rgkmax of 7.0 might take half an hour). After editing the file, save it and exit with Ctrl+O, Ctrl+X, as already pointed out before. 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 table and/or a plot of your results of ngridk and total energy.

Questions: What do you conclude about the convergence of this property (total energy) with respect to ngridk? Is the convergence fast or slow?

ii) Lattice parameter

In reality, it does not make much sense to select the total energy as the quantity to analyze the convergence behavior of ngridk and rgkmax. The major reason is that the reference (the zero) for the energy is usually not well defined when we deal with periodic boundary conditions, and, therefore, energy differences are the relevant quantities. This is what we will do now, theoretically predicting the lattice parameter of graphene.

To prepare the calculation, execute

mkdir -p /home/$myid/03-test-Lattice
cd /home/$myid/03-test-Lattice
chmod +x

This should create a new directory, move inside it, download the file, and make it executable.

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, 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 different lattice parameters, and for each of them, run 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 is printed in the last line, and represents the equilibrium lattice parameter found for ngridk and rgkmax selected before. A plot of the energy vs. the lattice parameter 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, and then execute it. In the end, for the selected rgkmax and ngridk, you will obtain the corresponding lattice parameter.

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

7. 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 -p /home/$myid/04-test-LossFunction
cd /home/$myid/04-test-LossFunction
cp /home/$myid/01-test/C.xml .

Create the file input.xml with nano.

nano input.xml

Your input file should be like the one given below. Do not forget to
  • Modify the lattice parameter to one that you obtained previously.
  • Replace question marks listed for rgkmax and ngridk with your choices.
  • Include the element xs just after groundstate (see Input Reference and Excited states from TDDFT for more information about the meaning of the different attributes).
   <structure speciespath="./">
      <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"
      xstype ="TDDFT"
      ngridk="? ? 1"
      vkloff="0.097 0.273 0.493"
         intv="0.0 1.0"
         points="500" />
          <qpoint> 0.0 0.0 0.0 </qpoint>

Then save this file, exit, going back to the command line.
Before running a calculation, it is advisable to create a directory that identifies and enter it. In this examples, we call it nk1_rgk1.0_gq1.0_nempty80. Please adjust the numbers to reflect your choices.

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

This procedure avoids that you discard the output files, after varying the parameters when checking for convergence. To run the calculation, execute exciting the command:
time /home/excitingser

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 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, it is time to verify how sensitive your results are with respect to the following parameters:

  • rgkmax
  • ngridk
  • nempty
  • gqmax

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

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

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

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

Questions: What can you comment about the precision of your loss function? What about the accuracy? You can compare your calculated loss function with the experimental one, presented in the paper of T. Eberlein et al. How can you improve the accuracy of your calculations?

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