CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
======= 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_pdos|calculating the PDOS]]:

<code - 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
</code>

<note important>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.</note>

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 [(http://journals.aps.org/prb/abstract/10.1103/PhysRevB.13.5188)] 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 [[https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/KPOINTS/BAND_STRUCTURE/KPOINT_SET.html|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 [[http://www.sciencedirect.com/science/article/pii/S0927025610002697|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:

<code>
 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
[...]
</code>

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.

<file python 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))

</file>