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.

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.

### 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 EXCITINGTOOLS=$EXCITINGROOT/tools export EXCITINGSCRIPTS=$EXCITINGROOT/scripts
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).

<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>  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: SETUP-excitingroot.sh  ##### 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
cp $EXCITINGSCRIPTS/script-Energy-rgkmax.sh .  The file script-Energy-rgkmax.sh 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 nano script-Energy-rgkmax.sh  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 ./script-Energy-rgkmax.sh  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  -76.1130433503  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
cp $EXCITINGSCRIPTS/script-Energy-ngridk.sh .  Edit the script to specify the value of rgkmax nano script-Energy-ngridk.sh  Look for the lines rgkmax=1.0 ... 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 ./script-Energy-ngridk.sh  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: tail -2 6/TOTENERGY.OUT  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 cp$EXCITINGSCRIPTS/script-Lattice.sh .
cp $EXCITINGSCRIPTS/pre-processing.py . cp$EXCITINGSCRIPTS/post-processing.py .


This should create a new directory, move inside it, and copy the script script-Lattice.sh.

Now, open the file script-Lattice.sh with nano

nano script-Lattice.sh


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

ngridk=1
rgkmax=1.0


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:

./script-Lattice.sh


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.

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 script-Lattice.sh, 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.

<input>

<title>Graphene</title>

<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> </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="? ? 1" rgkmax="?" swidth="0.01"> </groundstate> <xs xstype="TDDFT" rgkmax="?" ngridk="? ? 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> </input>  Save this file and exit. We also must set the species path by running: SETUP-excitingroot.sh  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?