Getting the band structure of graphene

To get the band structure for graphene (or h-BN), only a few changes are required compared to the previous example for calculating the PDOS:

graphene_kp_dos.inp

&GLOBAL
  PROJECT graphene_kp_dos
  RUN_TYPE ENERGY 
  PRINT_LEVEL MEDIUM
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  POTENTIAL

    &POISSON
      PERIODIC XYZ
    &END POISSON
    &QS
      EXTRAPOLATION USE_GUESS ! required for K-Point sampling
    &END QS
    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-6
      MAX_SCF 300

      ADDED_MOS 2
      &SMEAR ON
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE [K] 300
      &END SMEAR
      &DIAGONALIZATION
        ALGORITHM STANDARD
        EPS_ADAPT 0.01
      &END DIAGONALIZATION
      &MIXING
        METHOD BROYDEN_MIXING
        ALPHA 0.2
        BETA 1.5
        NBROYDEN 8
      &END MIXING
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
    &KPOINTS
      SCHEME MONKHORST-PACK 3 3 1
      SYMMETRY OFF 
      WAVEFUNCTIONS REAL
      FULL_GRID .TRUE.
      PARALLEL_GROUP_SIZE  0
      &BAND_STRUCTURE
        ADDED_MOS 2
        FILE_NAME graphene.bs
        &KPOINT_SET
          UNITS CART_BOHR ! work around a bug in CP2K, should be B_VECTOR
          SPECIAL_POINT 0.0   0.0   0.0
          SPECIAL_POINT 1./2. 0.0   0.0
          NPOINTS 5
        &END
      &END BAND_STRUCTURE
    &END KPOINTS
  &END DFT

  &SUBSYS
    &CELL
      ABC [angstrom] 2.4612 2.4612 8
      ALPHA_BETA_GAMMA 90. 90. 60.
      SYMMETRY HEXAGONAL
      PERIODIC XYZ
      MULTIPLE_UNIT_CELL 1 1 1
    &END CELL
    &TOPOLOGY
      MULTIPLE_UNIT_CELL 1 1 1
    &END TOPOLOGY
    &COORD
      SCALED
      C  1./3.  1./3.  0.
      C  2./3.  2./3.  0.
    &END
    &KIND C
      ELEMENT C
      BASIS_SET TZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE
    &END KIND
  &END SUBSYS

&END FORCE_EVAL

At present (CP2K 4.1) it is not possible to get the projected density of states when doing a K-Point calculation. Plus there is currently an issue with the UNITS specification for the special point coordinates: even though the unit is set to Cartesian coordinates (in Bohr), the special points are multiplied by the reciprocal vectors and must therefore be given in terms of the b-vectors.

Some notes on the input file:

By specifying the KPOINT section you are enabling the K-Point calculation.
While you could specify the K-Points directly, we are using the Monkhorst-Pack scheme ¹⁾ to generate them. The numbers following the keyword MONKHORST-PACK specify the tiling of the brillouin zone.
After the basic calculation, CP2K calculates the energies along certain lines, denoted as KPOINT_SET (when you check the documentation you will note that this section can be repeated).
The keyword NPOINTS specifies how many points (in the addition to the starting point) should be sampled between two special points.
The SPECIAL_POINT keyword is used to specify the start-, mid- and endpoints of the line. Those points usually denote special points in the reciprocal lattice/unit cell, like the $\Gamma$, $M$ or $K$ point. You can find the definition for these in the appendix section of this paper. This keyword can also be specified multiple times, making it possible to directly get the band structure for a complete path.

Now, when you run this input file you will get in addition the the output file, a file named graphene.bs which will look similar to the following:

 SET:       1                 TOTAL POINTS:      6
   POINT   1                     0.000000    0.000000    0.000000
   POINT   2                     0.500000    0.000000    0.000000
       Nr.    1    Spin 1        K-Point  0.00000000  0.00000000  0.00000000
                8
           -15.30752034     -3.31285773      0.93143545      1.03651421
             8.71874068     12.74920179     12.83785311     15.50144316
       Nr.    2    Spin 1        K-Point  0.02500000  0.00000000  0.00000000
                8
           -15.29453364     -3.29547462      0.87472486      1.00321991
             8.31998068     12.81500348     12.93001933     15.45108207
       Nr.    3    Spin 1        K-Point  0.05000000  0.00000000  0.00000000
[...]

For each set there is a block named SET with the special points listed as POINT, followed by sub-blocks for each K-Point containing the energies for each MO.

Your tasks:

Lookup the special points for the $\Gamma$, $M$, $K$ points in the mentioned paper (make sure you choose the right lattice). Calculate and plot the band structure for graphene from $\Gamma$ over $M$, $K$ back to $\Gamma$ (you are free to decide whether to use multiple K-Point sets are multiple special points in a single set). Mark the special points. Choose an appropriate number of interpolation points to get a smooth plot.
Compare your plot with plots from literature. What is different?
Why do you get 8 orbital energies? Try to change the input to get more unoccupied orbitals.

To convert the band structure file to a file which can be loaded directly into MATLAB for example, you can use the script cp2k_bs2csv.py from below, which when passed a band structure file graphene.bs as an argument will write files graphene.bs-setN.csv for each set containing the K-Point coordinates and the energies in one line.

cp2k_bs2csv.py

#!/usr/bin/env python
"""
Convert the CP2K band structure output to CSV files
"""
 
import re
import argparse
 
SET_MATCH = re.compile(r'''
[ ]*
  SET: [ ]* (?P<setnr>\d+) [ ]*
  TOTAL [ ] POINTS: [ ]* (?P<totalpoints>\d+) [ ]*
  \n
(?P<content>
  [\s\S]*?(?=\n.*?[ ] SET|$)  # match everything until next 'SET' or EOL
)
''', re.VERBOSE)
 
SPOINTS_MATCH = re.compile(r'''
[ ]*
  POINT [ ]+ (?P<pointnr>\d+) [ ]+ (?P<a>\S+) [ ]+ (?P<b>\S+) [ ]+ (?P<c>\S+)
''', re.VERBOSE)
 
POINTS_MATCH = re.compile(r'''
[ ]*
  Nr\. [ ]+ (?P<nr>\d+) [ ]+
  Spin [ ]+ (?P<spin>\d+) [ ]+
  K-Point [ ]+ (?P<a>\S+) [ ]+ (?P<b>\S+) [ ]+ (?P<c>\S+) [ ]*
  \n
[ ]* (?P<npoints>\d+) [ ]* \n
(?P<values>
  [\s\S]*?(?=\n.*?[ ] Nr|$)  # match everything until next 'Nr.' or EOL
)
''', re.VERBOSE)
 
if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('bsfilename', metavar='bandstructure-file', type=str,
                        help="the band structure file generated by CP2K")
 
    args = parser.parse_args()
 
    with open(args.bsfilename, 'r') as fhandle:
        for kpoint_set in SET_MATCH.finditer(fhandle.read()):
            filename = "{}.set-{}.csv".format(args.bsfilename,
                                              kpoint_set.group('setnr'))
            set_content = kpoint_set.group('content')
 
            with open(filename, 'w') as csvout:
                print(("writing point set {}"
                       " (total number of k-points: {totalpoints})"
                       .format(filename, **kpoint_set.groupdict())))
 
                print("  with the following special points:")
                for point in SPOINTS_MATCH.finditer(set_content):
                    print("  {pointnr}: {a}/{b}/{c}".format(
                        **point.groupdict()))
 
                for point in POINTS_MATCH.finditer(set_content):
                    results = point.groupdict()
                    results['values'] = " ".join(results['values'].split())
                    csvout.write("{a} {b} {c} {values}\n".format(**results))

¹⁾ http://journals.aps.org/prb/abstract/10.1103/PhysRevB.13.5188