Exercise 3: Electronic band structure of monolayer MoS$_\text{2}$

The electronic band structure gives the relation between crystal momentum $\mathbf{k}$, band index $n$ and the energy $\varepsilon_{n \mathbf{k}}$ of an electron in a crystal, according to Bloch's theorem.

The electronic band structure can be computed using DFT in an approximate way via the Kohn-Sham equations

$$ \left( -\frac{\nabla^2}{2m} + v_\text{ext}(\mathbf{r}) + v_\text{Hartree}(\mathbf{r}) + v_\text{xc}(\mathbf{r}) \right) \psi_{n\mathbf{k}}(\mathbf{r}) = \varepsilon_{n\mathbf{k}} \psi_{n\mathbf{k}}(\mathbf{r}) $$

When using the standard exchange-correlation (xc) functionals like PBE, the band gap between the occupied valence bands and the empty conduction bands is usually underestimated with respect to experiment.

Nevertheless, PBE often gives the correct band ordering, dispersions (i.e., curvature as function of $\mathbf{k}$), and orbital character of the bands.

In this exercise, we compute the band structure of monolayer MoS$_2$, a two-dimensional crystal which has been discovered in 2010 (doi:10.1103/PhysRevLett.105.136805).

The following input file can be used as a starting point:

DFT_bandstructure.inp
&GLOBAL
 PROJECT  MoS2
 RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
 &DFT
   BASIS_SET_FILE_NAME BASIS_MOLOPT
   POTENTIAL_FILE_NAME GTH_POTENTIALS
   &MGRID
     CUTOFF  500
     REL_CUTOFF  100
   &END MGRID
   &QS
     METHOD GPW
     EPS_DEFAULT 1.0E-12
     EPS_PGF_ORB 1.0E-12
   &END QS
   &SCF
     SCF_GUESS ATOMIC
     EPS_SCF 1.0E-9
     MAX_SCF 500
     &MIXING
         METHOD BROYDEN_MIXING
         ALPHA 0.1
         BETA 1.5
         NBROYDEN 8
     &END
   &END SCF
   &XC
     &XC_FUNCTIONAL PBE
     &END XC_FUNCTIONAL
   &END XC
   &KPOINTS
     SCHEME MONKHORST-PACK 8 8 1
   &END KPOINTS
   &PRINT
     &BAND_STRUCTURE
       ADDED_MOS 10
       FILE_NAME bandstructure.bs
       &KPOINT_SET
         NPOINTS 49
         SPECIAL_POINT GAMMA 0.0 0.0 0.0
         SPECIAL_POINT K     0.333333 0.333333 0.0
         SPECIAL_POINT M     0.0 0.5 0.0
         SPECIAL_POINT GAMMA 0.0 0.0 0.0
       &END KPOINT_SET
     &END BAND_STRUCTURE
   &END PRINT
 &END DFT
 &SUBSYS
   &CELL
     ABC                A B C            ! Comment: FILL HERE "LENGTHS" A B C from C2DB (unit: Angström)
     ALPHA_BETA_GAMMA   ALPHA BETA GAMMA ! Comment: FILL HERE "ANGLES" ALPHA BETA GAMMA from C2DB
     PERIODIC XY
   &END CELL
       &KIND S
     BASIS_SET  DZVP-MOLOPT-GTH
     POTENTIAL  GTH-PBE
   &END KIND
 
   &KIND Mo
     BASIS_SET  DZVP-MOLOPT-SR-GTH
     POTENTIAL  GTH-PBE
   &END KIND
 
   &COORD
Mo       X Y Z  ! Comment: FILL HERE POSITION OF Mo FROM XYZ FILE FROM C2DB
S        X Y Z  ! Comment: FILL HERE POSITION OF S  FROM XYZ FILE FROM C2DB
S        X Y Z  ! Comment: FILL HERE POSITION OF S  FROM XYZ FILE FROM C2DB
   &END COORD
 
   &TOPOLOGY
     &CENTER_COORDINATES
     &END
   &END
 
 &END SUBSYS
&END FORCE_EVAL

To complete it, the atomic positions and cell size need to be filled for MoS$_\text{2}$.

Task: Navigate to the Computational 2D Materials Database (C2DB) and search for MoS$_\text{2}$. Download the “XYZ” file which contains the atomic positions. Also fill in the lengths of the cell and the angles, as given on the CSDB website. Then run CP2K: 
cp2k.psmp DFT_Bandstructure.inp | tee cp2k.out

You can validate the output by comparing to the solution available here.

The PBE band structure is contained in the file bandstructure.bs. You can plot the PBE band structure using a plotting script available via github; obtain the script via 

git clone https://github.com/stefabat/cp2k-scripts
Task: Run the plotting script via 
python3 cp2k-scripts/bin/cp2k_plot_bands.py bandstructure.bs --energy_range -2 3

and compare your plot to the png file provided in the solution directory.

You can also compare the $k$-path we have chosen in the input to the $k$-path of the hexagonal crystal structure available in the appendix of this paper.

Task: Compare your DFT/PBE band structure to reference data, for example from C2DB. Evaluate the convergence parameters, specifically:
  • Basis set size: You can test larger basis sets like TZV2P-MOLOPT-PBE-GTH-q6 (for S) and TZV2P-MOLOPT-PBE-GTH-q14 (for Mo), contained in the basis file BASIS_MOLOPT_UZH (also found here).
  • $k$-point mesh: You can test larger $k$-grids, for example SCHEME MONKHORST-PACK 16 16 1 instead of SCHEME MONKHORST-PACK 8 8 1
Task: Compare your DFT/PBE band structure to the $GW$ method, a higher level method deriving from Green's function theory 10.3389/fchem.2019.00377, for example from C2DB. What is the main difference of the $GW$ bandstructure compared to DFT/PBE?

In case you would like to execute the $GW$ band structure calculations of 2D materials with CP2K, you can check out 10.48550/arXiv.2507.18411 (Fig. 1,4,10 and Table I) and the corresponding input and output files on github.

Can you identify what is the effect of spin-orbit coupling on the band structure described in this paper?