CP2K Open Source Molecular Dynamics

Here we give some notes on how to use QUANTUM ESPRESSO to perform one of the standard tasks e.g. band structure calculation. In this exercise, we calculate the band structure of graphene along the high-symmetry lines (Γ-M-K-Γ) of the Brillouin zone (BZ). First, we do a self-consistent field (scf) calculation and then 'bands' calculation with a set of k-points (along high-symmetry lines). Please have a look at the end of the input files, scf.in, and bands.in. In general, to study the electronic band structure one must 'relax' the ions in the system (not required for this exercise!) and then do the post-processing.

Steps:

1. Before starting the calculations please load following modules. Please type the commands at the terminal.

 $ module load intel/14.0.1 open_mpi/1.6.5
 $ module load quantum_espresso

2. Download all the commented files from the media manager: exercise-12.tar.gz or

Self-Consistent Field (SCF) calculation:

The SCF calculation means that the program iteratively updates the orbitals until certain criterion is reached. Inspect the input file “scf.in” for graphene. The most important input parameters in scf.in are: definition of a crystal system, lattice constant, k-grid, and vacuum region to avoid interactions between two graphene sheets. In our case, a vacuum along the perpendicular direction is ~12 Å. The symmetry has been included (ibrav) to reduce the number of calculations.

 ibrav = 4                                  ! definition of the crystal system (hexagonal for graphene
 celldm(1) = 4.65303296132007217931         ! lattice constant in a.u. (2.462/0.529117249) where 2.462 is the lattice const. in Angstroms
 celldm(3) = 4.87408610885458976441         ! c/a (12/2.462), where c is the vacuum along z-direction

The types of species (only “C” for graphene) is defined by the section,

ATOMIC_SPECIES
C   6    C.pbe-rrkjus.UPF

Quantum espresso will recognize the type of pseudopotential and exchange correlation from the file. Here, for graphene, we use ultrasoft pseudopotential (rrkjus) with PBE exchange-correlation functional (pbe). For more details on the naming conventions please visit QEPPS

The atomic positions are defined in crystal format as,

 ATOMIC_POSITIONS crystal
 C        0.666666666   0.333333333   0.000000000
 C        0.333333333   0.666666666   0.000000000

This is a ground-state calculation for graphene with two atoms per unit cell. For many properties, e.g. density of states (DOS), response functions, charge density etc., the integrals over the Brillouin Zone (BZ) are necessary. To evaluate the integrals computationally one must replace them with weighted sum over special k-points. The number of points in the k-space (i.e. k-grid) at which BZ is sampled determines the quality of the calculation. The choice of the k-point mesh is system dependent. More grid points (smaller grid spacing), gives better convergence of the total energy. The line “K_POINTS automatic” defines the way in which the BZ is to be sampled. Here we sample the BZ using Monkhorst-Pack scheme to get equally spaced mesh.

 K_POINTS automatic
 KX KY KZ 0 0 0

3. Insert the Γ-centred k-grid KX KY KZ (e.g. 10 10 1). As we are dealing with two-dimensional (2D) system (and so 2D BZ), KZ should be 1. There are two choices for the centre of the mesh: 1) centred on Γ (Γ belongs to the mesh) and 2) centred around Γ (Γ does not belong to the mesh). The second choice can break the symmetry!!. Here we use Γ centred k-grid.

4. Run the ground-state calculation using pw.x code of the QUANTUM ESPRESSO suit,

 bsub -n 4 " mpirun pw.x < scf.in > SCF.out "

5. When the calculation finishes, check for errors and/or warnings. If everything is correct then note down the Fermi energy from the output file (SCF.out).

 $ grep “Fermi energy”  SCF.out | tail -1

6. Copy the mol.save directory to scf-mol.save. The scf-mol.save is required for the 'bands' calculation

 $ cp -rf mol.save/ scf-mol.save/

Non-Self-Consistent Field (NSCF) calculation:

In non-self-consistent calculation, the potential is constructed from some “input” charge density and remains fixed. Calculation of DOS can be done in two ways:

a) For finite geometries (e.g. molecules), where a single k-point (centre of BZ) is sufficient, the simple way is to perform self-consistent calculation and then DOS calculation.

b) For periodic geometries, a high quality DOS might require very fine meshes and for large cells one might need many k-points (depending on the system). Therefore, to save the computational time it is a good idea to calculate the self-consistent charge density with few k-points and then non-self-consistent calculation using fixed self-consistent charge density.

7. For density of states calculation, do the non-self-consistent calculation using the input file “nscf.in” (change the k-grid KX KY KZ) and prefix = './mol'.

 $ bsub -n 4 " mpirun pw.x < nscf.in > NSCF.out "

8. and then use dos.in input file to get the density of states from -20 to 10 eV.

 $ bsub -n 4 " mpirun dos.x < dos.in > DOS.out "

First column of “graphene.dos” file is the energy and second column is the total DOS. Open plot_dos.plt script to add Fermi energy and save the plot using gnuplot. The DOS should look like this,

9. For band structure calculation, use the input file 'bands.in' and the ground-state density obtained from the 'scf' run (prefix = './scf-mol'). The BZ for graphene is shown in figure below. The points Γ, K and M are called the zone centre, the corner and the centre of the edge respectively. The green lines show the borders of the irreducible BZ along which the band extrema occurs. Therefore, we move along these lines to get the energy that an electron can have within the solid. Now the k-grid is replaced with a list of high-symmetry points along Γ-M-K-Γ directions. You can run the kpoints program and get a list of kpoints along the high-symmetry lines. For this you need to compile the kpoints.c program using gcc compiler

 $ gcc -Wall kpoints.c -o kpoints

and then

 $ ./kpoints

and enter starting and ending points. For example if you want to generate a string of points from Γ to M, the coordinates of Γ- and M- points should be (0 0 0) and (0.5 0 0) and number of points = 100. These coordinates are in units of reciprocal lattice vectors. In the present case the reciprocal lattice vectors have been determined using the known formulae [see QE-exercise.pdf] and the coordinates of Γ-, K- and M-points are (0,0,0), (1/3, 1/3,0) and (1/2,0,0) respectively [see the figure of BZ for graphene]. For the present band structure calculation there is no need to enter a string of kpoints.

 $ bsub -n 4 " mpirun pw.x < bands.in > BANDS.out "

10. Ones the job is over, arrange the band energies according to the number of k-points. This means you need to write the data (from BANDS.out)

        k = 0.0000 0.0000 0.0000     band energies (ev):
        
        -20.6085  -8.7205  -4.1559  -4.1559   2.1035   3.2561   3.9941   7.1819
         7.2173   7.2173   8.4394  10.3628
         
         
        k = 0.0051 0.0029 0.0000     band energies (ev):
        
        -20.6079  -8.7198  -4.1583  -4.1572   2.1044   3.2569   3.9950   7.1827
         7.2177   7.2200   8.4402  10.3608
         
         ......
         ......
         ......

in two columns. The first column is the k-point number (N = 230 in bands.in) and second column is the energy eigenvalue (E) e.g. in present case,

   Sr. No    Energy [eV]
     1        	E1 (K1)
     2        	E1 (K2)
     3        	E1 (K3)
     …        	...
     …        	...
     N        	E1 (KN)
    
     1        	E2 (K1)
     2        	E2 (K2)
     3        	E2 (K3)
     …        	...
     …        	...
     N        	E2 (KN)
     .           .
     .           .
     .           .
     .           .
     .           .
    
     1        	EN (K1)
     2        	EN (K2)
     3        	EN (K3)
     …        	...
     …        	...
     N        	EN (KN)

11. Use Qe_Bands.sh script to do this job

 $ bash Qe_Bands.sh BANDS.out <FermiEnergy>

12. Plot (column 1 Vs column 2 of EnergyValues.txt) and save the band structure using gnu plot. Compare your band structure with the image provided here. The Fermi energy is set to zero. Do you see the linear dispersion of bands at K-point (+- 1 eV around the Fermi energy)?

 $ gnuplot plot_band.plt

13. Before going for the different k-grid please delete all mol.save directories.

  $ rm -rf mol.save/
  $ rm -rf scf-mol.save/

Don't forget to visit QUANTUM ESPRESSO for more details about the input parameters!