Crivella · simonebrozzesi · Nov 14, 2021 · Nov 18, 2021 · Nov 18, 2021 · Nov 18, 2021
diff --git a/.DS_Store b/.DS_Store
diff --git a/0_Si_bulk/.DS_Store b/0_Si_bulk/.DS_Store
diff --git a/0_Si_bulk/0_cutoff/README.md b/0_Si_bulk/0_cutoff/README.md
@@ -1,38 +1,103 @@
-# Convergency on the kinetic cutoff
-QE expands the wavefunctions on a plane wave(PW) basis set.
-As with any set, an infinite number of function is required to perfectly represent the original function.
-We need to limit the number of PW to the smallest number that represent our system accurately
-Look at the file Docs/DFT1.pdf page 16-18
-  1. Run the calculation using the provided input file:
-      ```
-      pw.x < si.scf-manual.in > si.scf-manual.out
-      or to use multiple processors (parallel run) (require parallel version to be compiled)
-      mpirun -np 2 pw.x < si.scf-manual.in > si.scf-manual.out
-      ```
-  2. Look at the information presented on the output file ('Docs/OUTPUT_PW' run with verbosity = 'high')
-  		- An header containing information of the version of espresso used
+# Convergence on the kinetic energy cutoff
+QE expands the wavefunctions on a plane wave (PW) basis set.
+As with any basis set, an infinite number of functions is required to perfectly represent the original function.
+We need to limit the number of PW to the smallest number that represents our system accurately.
+
+  1. Look at the provided input file for bulk silicon and visualize it with `Xcrysden`. A brief explanation of variables is given in the file. Complete details are given on the [QE website](https://www.quantum-espresso.org/Doc/INPUT_PW.html).
+      ```
+      % cat si.scf.in
+      % xcrysden --pwi si.scf.in
+      ```
+      Now run the calculation, redirecting the output to a file:
+      ```
+      % pw.x < si.scf.in > si.scf.out
+      ```
+      or to use multiple processors, if quantum-espresso has been compiled in parallel
+      ```
+      % mpirun -np 2 pw.x < si.scf.in > si.scf.out
+      ```
+  2. Look at the information presented on the output file 
+      - A header containing information of the version of espresso used
       - A recap on the information of the system
       - A list of matrix representation of the symmetries of the system + the character table of the symmetries
-      - The list of the point in k-space used
+      - The list of the points in k-space used
       - The information on each self-consistent iteration
       - The eigenvalues and occupations for the requested Kohn-Sham states at every k-point
       - The total energy marked by a ! ( as a sum of different contributions)
       - The forces acting on the atoms (should be 0 for a system at the equilibrium)
-      - Information on the time for each subroutine (Usefule for big calculation)
-  3. Repeat step 1 and change each time the value of ecutwfc from 5 up to 30 and the name of the output as to not override them
+      - Information on the total time and time for each subroutine. Note that QE often uses Rydberg and bohr atomic units: 1 Ry = 13.6057eV, 1 bohr = 0.529177 Angstrom.
+  3. Repeat step 1 and change each time the value of `ecutwfc` from 5 up to 30 Ry and the name of the output (so as to not overwrite them)
+      ```
+      % pw.x < si.scf.in > si.scf.out_5Ry
+      ```
+      Modify (edit) the `si.scf.in` file and change `ecutwfc = 5` to `ecutwfc = 10`
+      ```
+      % pw.x < si.scf.in > si.scf.out_10Ry
+      ```
+      You can also make this change directly to the original input file using the `sed` command:
+      ```
+      % sed -e 's/ecutwfc   = 5/ecutwfc   = 10/' si.scf.in > si.scf.in_10Ry
+      % grep 'ecutwfc' si.scf.in_10Ry 
+      ecutwfc   = 10,
+      % pw.x < si.scf.in_10Ry > si.scf.out_10Ry
+      ```
+      Repeat for 15,20,25 and 30 Ry.
+      ```
+      [...]
+      % pw.x < si.scf.in > si.scf.out_30Ry
+      ```
+  4. Use grep to print out the total energies from all files in one command using the wildcard *
+      ```
+      % grep -e '!' *out*Ry
+      si.scf.out_10Ry:!    total energy              =     -15.77444885 Ry
+      [...]
       ```
-      pw.x < si.scf-manual.in > si.scf-manual.out_5
-      change ecutwfc 5 -> 10
-      pw.x < si.scf-manual.in > si.scf-manual.out_10
-      ...
-      change ecutwfc 25 -> 30
-      pw.x < si.scf-manual.in > si.scf-manual.out_30
-  4. Use grep to collect the total energies from all files in one command
+     Copy and paste the cutoff energies and total energies into a 2 column file (Ecut,Etot) called 'Etot_vs_Ecut.dat' and plot it to see if you have reached convergence. 
+     ```
+     gnuplot> plot "Etot_vs_Ecut.dat" w l
+     ```
+     To the eye, the total energy looks well converged at 30 Ry. This is misleading however, as it simply depends on the energy range of your plot - see the logscale plot on the right. The correct threshold to use depends on what quantity (observable) you want to compute: a bond length, a lattice parameter, an energy gap, a binding energy, an STM image, a vibrational frequency...and to what precision you want!
+
+     Here are some rough guidelines on convergence (not for production/publication!):
+     * Total energy: 1 mRy/atom = 13.6 meV/atom (0.1 mRy/atom is better)
+     * Bond length: 0.002 Å
+     * Cell parameters: 0.01 Å
+     * Energy differences: 0.37 mRy/atom = 5 meV/atom (0.1 mRy/atom is better)
+     * Forces: 10 meV/Å
+
+     _NB: you must ALWAYS perform convergence tests yourself for your system!!_
+
+
+     ![Total energy vs kinetic energy cutoff](Ref/Etot_vs_Ecut.png?raw=true "Total energy vs kinetic energy cutoff")
+
+     If the plot doesn't look right, make sure you have used the right cutoff in the appropriate input and output files:
+     ```
+     % grep 'kinetic-energy cutoff' si.scf.out_25Ry 
+     kinetic-energy cutoff     =      25.0000  Ry
+     % grep 'ecutwfc' si.scf.in_25Ry 
+     ecutwfc   = 25,
+     ```
+
+  6. Use grep on each file to extract the eigenvalues of the highest occupied and lowest unoccupied bands, and compute the band gap using the 'bc -l' program
       ```
-      grep -e "!" *
+      % grep -e 'highest' *out*Ry 
+      si.scf.out_5Ry:     highest occupied, lowest unoccupied level (ev):     6.0229    7.5083
+      % bc -l
+      bc 1.06
+      Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
+      This is free software with ABSOLUTELY NO WARRANTY.
+      For details type `warranty'. 
+      7.5083-6.0229 
+      1.4854
+      quit
       ```
-  5. Save the cutoff energies and total energies in a 2 column file and plot it to see if you have reached convergence
-  6. The script 'run_ecut' will do automatically everything from step 3 to 5 (explained using comments inside the script)
+     Save the energies in a 4 column file (Ecut, VBM, CBM, gap) called 'Gap_vs_Ecut.dat' and plot them versus the cutoff. How does the convergence compare with the value expected from the total energy run?
+     ![Eigenvalues and gap vs kinetic energy cutoff](Ref/Gap_vs_Ecut.png?raw=true "Gap vs kinetic energy cutoff")
+  7. ADVANCED USERS: The shell scripts 'run_ecut' and 'run_plots' in the 'Script' directory will do everything automatically from step 3 to 5 (explained using comments inside the script). You may first have to modify the ENVIRONMENT_VARIABLES file in the root directory. The scripts must be run from the main 0_cutoff directory (or copied there). Hit 'q' to cycle through plot windows as they appear. Inspect the PDFs or PNGs that are created.
       ```
-      ./run_ecut
+      ./Script/run_ecut
+      ./Script/run_plots
       ```
+      NB: Do not use the scripts for your own projects unless you understand well how they work!
+
+### When you have completed this tutorial, you can move on to [1_kpoints: Convergence with k-points](../1_kpoints)
diff --git a/0_Si_bulk/0_cutoff/Ref/E_vs_ecut.pdf b/0_Si_bulk/0_cutoff/Ref/E_vs_ecut.pdf
diff --git a/0_Si_bulk/0_cutoff/Ref/E_vs_ecut.dat → 0_Si_bulk/0_cutoff/Ref/Etot_vs_Ecut.dat b/0_Si_bulk/0_cutoff/Ref/E_vs_ecut.dat → 0_Si_bulk/0_cutoff/Ref/Etot_vs_Ecut.dat
@@ -1,7 +1,7 @@
 # ecut(Ry)	Energy(Ry)
 5		-15.60596061
-10		-15.77444902
-15		-15.83462803
+10		-15.77444885
+15		-15.83462812
 20		-15.84756981
 25		-15.85122085
-30		-15.85219549
+30		-15.85219551
diff --git a/0_Si_bulk/0_cutoff/Ref/Etot_vs_Ecut.pdf b/0_Si_bulk/0_cutoff/Ref/Etot_vs_Ecut.pdf
diff --git a/0_Si_bulk/0_cutoff/Ref/Etot_vs_Ecut.png b/0_Si_bulk/0_cutoff/Ref/Etot_vs_Ecut.png
diff --git a/0_Si_bulk/0_cutoff/Ref/Gap_vs_Ecut.dat b/0_Si_bulk/0_cutoff/Ref/Gap_vs_Ecut.dat
@@ -0,0 +1,7 @@
+#Ecut		VBM	CBM	Gap
+5		6.0229 7.5083 1.4854
+10		5.9615 7.1803 1.2188
+15		5.9315 7.0729 1.1414
+20		5.9210 7.0430 1.122
+25		5.9101 7.0316 1.1215
+30		5.9063 7.0268 1.1205
diff --git a/0_Si_bulk/0_cutoff/Ref/Gap_vs_Ecut.png b/0_Si_bulk/0_cutoff/Ref/Gap_vs_Ecut.png
diff --git a/0_Si_bulk/0_cutoff/Ref/Gap_vs_ecut.pdf b/0_Si_bulk/0_cutoff/Ref/Gap_vs_ecut.pdf
diff --git a/0_Si_bulk/0_cutoff/run_ecut → 0_Si_bulk/0_cutoff/Script/run_ecut b/0_Si_bulk/0_cutoff/run_ecut → 0_Si_bulk/0_cutoff/Script/run_ecut
@@ -9,20 +9,24 @@ echo "TMP_DIR:" $TMP_DIR
 echo "Parallel command:" $RUN_COMMAND
 echo "Started at: " `date`
 
-SAVE=E_vs_ecut.dat				#NAME=something assign a variable of name NAME and value something (can be a string ora  number)
+# Prepare the files to be plotted
+SAVE=Etot_vs_Ecut-script.dat				#NAME=something assign a variable of name NAME and value something (can be a string ora  number)
 						#!!!!!! There must be no spaces when assigning variable
 						#Use "... ..." for a string with spaces
-GAP=Gap_vs_ecut.dat
-echo -e "# ecut(Ry)\tEnergy(Ry)" > $SAVE	# COMMAND > NAME  write the stdout of the command COMMAND into the file NAME
+echo -e "# Ecut(Ry)\tEnergy(Ry)" > $SAVE	# COMMAND > NAME  write the stdout of the command COMMAND into the file NAME
 						#Create the file if it does not exist
 						#Overwrite the file it exist
 						#Use >> to append at the end of the file
 
+GAP=Gap_vs_Ecut-script.dat
+echo -e "# Ecut(Ry)\tVBM\tCBM\tGap (eV)" > $GAP	
+
+# Loop over cutoff energies
 for ecut in 5 10 15 20 25 30; do		#for NAME in LIST; do ... done
 						#cicle the variable NAME in a list of elements
 
 IN=si_script.scf.in
-OUT=si.scf.out_ecut$ecut
+OUT=si_script.scf.out_Ecut$ecut
 						#<<DELIM  ... DELIM   defines a block of text
 						#Will still resolve variables values using $NAME
 						#cat NAME/TEXT   read a file or string/block of text and print it into the stdout
@@ -73,36 +77,8 @@ echo -e "$ecut\t\t$ENERGY" >> $SAVE
 BANDS=`cat $OUT | grep highest | awk '{print $7,$8,$8-$7}'`
 echo -e "$ecut\t\t$BANDS" >> $GAP
 
-done
+done    # end of main loop over cutoff energy
 
 echo "Run completed at: " `date`
-
-#Print a gnuplot script to plot the 2 column files
-cat > plot.gnu <<EOF
-#!/usr/bin/gnuplot
-
-set term wxt enhanced
-
-set xlabel "E_{cut} (Ry)"
-set ylabel "E_{tot} (Ry)
-
-plot \
-	'$SAVE'  u 1:2 w lp pt 7 lw 2  lc rgb "black" title "E vs E_{cut}" 
-
-pause -1
-set term pdfcairo enhanced
-set output "${SAVE}.pdf"
-replot
-EOF
-
-#Run the gnuplot script to make the plot
-gnuplot plot.gnu
-
-
-
-exit
-
-
-
-
+echo "Now run the run_plots script" 
 
diff --git a/0_Si_bulk/0_cutoff/Script/run_plots b/0_Si_bulk/0_cutoff/Script/run_plots
@@ -0,0 +1,120 @@
+#!/usr/bin/env bash
+
+SAVE=Etot_vs_Ecut-script.dat			
+GAP=Gap_vs_Ecut-script.dat
+
+#Print a gnuplot script to plot the 2 column files
+cat > plot-energy-script.gnu <<EOF
+#!/usr/bin/gnuplot
+
+# x11 should work on any system
+# try wxt (ubuntu) or aqua (mac) for nicer plots
+reset
+set term x11 enhanced
+set out
+
+set xlabel "E_{cut} (Ry)"
+set ylabel "E_{tot} (Ry)
+thres_1mRy_per_atom=0.001*2
+thres_10mRy_per_atom=0.01*2
+
+plot \
+	'$SAVE'  u 1:2 w lp pt 7 lw 2  lc rgb "black" title "Total energy" 
+replot GPVAL_DATA_Y_MIN t "Min value" w l lt 4 dt 4 lc "red",\
+        GPVAL_DATA_Y_MIN+thres_1mRy_per_atom w l lt 2 lc "blue" t "Min + 1mRy/atom",\
+        GPVAL_DATA_Y_MIN+thres_10mRy_per_atom w l lt 3 lc "magenta" t "Min + 10mRy/atom"
+
+
+
+pause -1
+set term pdfcairo enhanced
+set output "${SAVE}.pdf"
+replot
+
+reset
+set term pngcairo enhanced size 800,400
+set output "${SAVE}.png"
+set xlabel "E_{cut} (Ry)"
+set ylabel "E_{tot} (Ry)
+thres_1mRy_per_atom=0.001*2
+thres_10mRy_per_atom=0.01*2
+stats '$SAVE' u 1:2 name 'EN' nooutput
+set multiplot 
+
+       set size 0.6,0.9
+       set origin 0.03,0.05
+
+plot \
+	'$SAVE'  u 1:2 w lp pt 7 lw 2  lc rgb "black" title "Total energy",\
+        EN_min_y t "Min value" w l lt 4 dt 4 lc "red",\
+        EN_min_y+thres_1mRy_per_atom w l lt 2 lc "blue" t "Min + 1mRy/atom",\
+        EN_min_y+thres_10mRy_per_atom w l lt 3 lc "magenta" t "Min + 10mRy/atom"
+
+       set size 0.38,0.9
+       set origin 0.63,0.05
+set ylabel "log_{10} (E_{tot}-E_{min}) (Ry)
+set xlabel "E_{cut} (Ry)"
+set grid y
+set autoscale y
+set logscale y
+set yrange [0.0001:1]
+plot \
+	'$SAVE'  u 1:(\$2-EN_min_y) w lp pt 7 lw 2  lc rgb "black" title "Total energy" 
+set nomulti
+pause -1
+set term x11
+set out
+replot
+pause -1
+
+EOF
+
+
+#Run the gnuplot script to make the plot of the total energy
+gnuplot plot-energy-script.gnu
+
+#Print a gnuplot script to plot the 4 column files
+cat > plot-gap-script.gnu <<EOF
+#!/usr/bin/gnuplot
+
+set xlabel "E_{cut} (Ry)"
+set ylabel "E (eV)
+
+# Generate a PDF image
+set term pdfcairo enhanced
+set output "${GAP}.pdf"
+set multiplot layout 1,2
+plot \
+	'$GAP'  u 1:3 w lp pt 7 lw 2  lc rgb "red" title "CBM",\
+	'$GAP'  u 1:2 w lp pt 6 lw 2  lc rgb "blue" title "VBM"
+plot \
+	'$GAP'  u 1:4 w lp pt 5 lw 2  lc rgb "black" title "Gap"
+set nomulti
+
+# Generate a PNG image (useful for GitHub)
+set term pngcairo enhanced size 800,400
+set output "${GAP}.png"
+set multiplot layout 1,2
+plot \
+	'$GAP'  u 1:3 w lp pt 7 lw 2  lc rgb "red" title "CBM",\
+	'$GAP'  u 1:2 w lp pt 6 lw 2  lc rgb "blue" title "VBM"
+plot \
+	'$GAP'  u 1:4 w lp pt 5 lw 2  lc rgb "black" title "Gap"
+set nomulti
+replot
+
+# Plot onscreen just the gap.
+# x11 should work on most system. Try wxt (ubuntu) or aqua (mac) for nicer plots.
+set term x11 enhanced
+set out
+plot \
+	'$GAP'  u 1:4 w lp pt 5 lw 2  lc rgb "black" title "Gap"
+replot
+pause -1
+EOF
+
+#Run the gnuplot script to make the plot
+gnuplot plot-gap-script.gnu
+
+exit
+
diff --git a/0_Si_bulk/0_cutoff/si.scf-manual.in → 0_Si_bulk/0_cutoff/si.scf.in b/0_Si_bulk/0_cutoff/si.scf-manual.in → 0_Si_bulk/0_cutoff/si.scf.in
@@ -1,33 +1,32 @@
 &CONTROL
   calculation  = "scf",			!Type of calculation to be run
   prefix       = "Si",			!Prefix that will mark the temporary datafiles. 
-					!KEEP THE SAME FOR CALCULATION ON THE SAME MATERIAL
-  pseudo_dir   = "../../Pseudo",	!Directory where psuedopotential files are located
+					!KEEP THE SAME FOR ALL CALCULATIONS ON THE SAME MATERIAL
+  pseudo_dir   = "../../Pseudo",	!Directory where pseudopotential files are located
   outdir       = "./tmp",		!Directory where temporary files are stored for further post-processing
   restart_mode = "from_scratch"		
   verbosity    = "high",		!Amount of details to be printed in the output
   tprnfor = .true.			!If .true. calculate forces on atoms
   wf_collect=.true.			!Collect all temporary files (strongly advise to use this expecally for parallel runs)
-					!Only downside, it will consume some time to collect and write the data (not relevant for this tutorial)
 /
 &SYSTEM
   ibrav     = 2,			!Number that specifies the Bravais lattice  2=FCC	
-  celldm(1) = 10.21,			!Specify the side lenght of the FCC lattice in [bohr]
+  celldm(1) = 10.21,			!Specify the side length of the FCC lattice in [bohr]
 					!Look at the documentation for the meaning of celldim(1/2/3) for all lattices
   nat       = 2,			!Total number of atoms in the unit cell
 					!The number of lines under the NAMELIST 'ATOMIC_POSITIONS' must be the same as 'nat'
   ntyp      = 1,			!The number of types of atoms (pseudopotentials loaded)
 					!The number of lines under the NAMELIST 'ATOMIC_SPECIES' must be the same as 'ntyp'
-  ecutwfc   = 5, 			!The kinetict cutoff on the wave functions in [Ry]
-  nbnd      = 5				!The number of Kohn-Sham states to be calculated
+  ecutwfc   = 5, 			!The kinetic energy cutoff on the wave functions in [Ry]
+  nbnd      = 6				!The number of Kohn-Sham states to be calculated (here 4 occupied, 1 empty)
 /
 &ELECTRONS
   conv_thr    = 1.D-8,			!Convergency threshold for self-consistency
   mixing_beta = 0.7D0	,		!Mix current charge density with previous. (Can influence convergence speed)
 /
 ATOMIC_SPECIES
  Si  28.086  Si.pz-vbc.UPF
-ATOMIC_POSITIONS
+ATOMIC_POSITIONS {alat}
  Si 0.00 0.00 0.00
  Si 0.25 0.25 0.25
 K_POINTS {automatic}

diff --git a/0_Si_bulk/1_alat/README.md b/0_Si_bulk/1_alat/README.md