[started by 2014.12.3]
[started by Pourya 2014.12.3]
**Introduction [#intro]
This is a Quantum Espresso tutorial in Saito Lab.

We can also learn from other tutorial provided by Saito(Susumu)-sensei's group at Titech, Tokyo:
http://www.stat.phys.titech.ac.jp/SATL_qe_tutorial/ (in Japanese!)

(Very rough) Translation is available here: 

**Prerequisites [#pre]
- We will install and run the program in our own home directory.
- Therefore it's better to prepare a directory specific for quantum espresso calculation, e.g. "espresso"
 mkdir ~/espresso
- Download quantum espresso package from (for version 5.0.2):
- If you use quantum espresso(v.5.0.2), please download the patch (e.g. as "espresso-5.0.2-5.0.3.diff"):

**How to Install [#install]
- We assume that espresso-5.0.2.tar.gz is already in the "espresso" directory
- Login to FLEX and then TUBE system, e.g.
 ssh tube50
- Run the following command under the "espresso" directory (to make sure, don't forget execute "cd ~/espresso" firstly):
 tar -xzvf espresso-5.0.2.tar.gz
 mv espresso-5.0.2 00main
 cd 00main
 patch -p1 < espresso-5.0.2-5.0.3.diff
You will be asked for input file to be patched, neglect it ans press enter some times, then
 ./configure --with-scalapack=NO
 make all
(note that we don't use scalapack in TUBE, therefore we set the option as "NO", otherwise the compilation will be failed)

**Structure of the input data  [#input]
- The structure of the input file of quantum espresso is 
 [ &IONS
 / ]
 [ &CELL
 / ]
 X  Mass_X  PseudoPot_X
 Y  Mass_Y  PseudoPot_Y
 Z  Mass_Z  PseudoPot_Z
 ATOMIC_POSITIONS { alat | bohr | crystal | angstrom }
  X 0.0  0.0  0.0  {if_pos(1) if_pos(2) if_pos(3)}
  Y 0.5  0.0  0.0
  Z O.0  0.2  0.2
 K_POINTS { tpiba | automatic | crystal | gamma | tpiba_b | crystal_b | tpiba_c | crystal_c }
 if (gamma)
   nothing to read
 if (automatic)
   nk1, nk2, nk3, k1, k2, k3
 if (not automatic)
   xk_x, xk_y, xk_z,  wk
 [ CELL_PARAMETERS { alat | bohr | angstrom }
   v1(1) v1(2) v1(3)
   v2(1) v2(2) v2(3)
   v3(1) v3(2) v3(3) ]
   f_inp1(1)  f_inp1(2)  f_inp1(3) ... f_inp1(10)
   f_inp1(11) f_inp1(12) ... f_inp1(nbnd)
 [ f_inp2(1)  f_inp2(2)  f_inp2(3) ... f_inp2(10)
   f_inp2(11) f_inp2(12) ... f_inp2(nbnd) ] ]
   nconstr  { constr_tol }
   constr_type(.)   constr(1,.)   constr(2,.) [ constr(3,.)   constr(4,.) ] {   constr_target(.) } ]

- These notations mean: Input data format: { } = optional, [ ] = it depends, | = or.

** Compile file in quantum espresso: [#compile]
- pw.x:
 calculates electronic structure, structural optimization, molecular dynamics, barriers with NEB.

- ph.x:
 calculates phonon frequencies and displacement patterns, dielectric tensors, effective charges (uses data produced by pw.x).

- pp.x:
 extracts the specified data from files produced by pw.x, prepare
data for plotting by writing them into formats that can be read by
several plotting programs.

- bands.x:
 extracts and reorders eigenvalues from files produced by pw.x for band structure plotting.

- projwfc.x:
 calculates projections of wavefunction over atomic orbitals,performs L?owdin population analysis and calculates projected density of states. These can be summed using auxiliary code sumpdos.x

- plotrho.x:
 produces PostScript 2-d contour plots.

- plotband.x:
 reads the output of bands.x, produces band structure PostScript plots.

- dos.x:
 calculates electronic Density of States (DOS).

** Where is the pseudopotentials ?: [#pseudo]
- http://www.quantum-espresso.org/pseudopotentials/

** Where can we find the Brillouin zone and High symmetry points? [#p52c2f19]
- Bilbao Crystallographic Server:

** Some Conversions: [#q42f9c19]
- These are some helpful conversions for calculation in quantum espresso
 1 bohr = 0.529177249 angstrom
 1 Rydberg (R ) = 13.6056981 eV
 1 eV =1.60217733 x 10-19 Joules

**Graphene [#Graphene]

- A single layer graphene is shown as following figure:


- The electron structure of graphene is : $1s^2$ $2s^2$  $2p^4$ and its radial wave function for each orbital is :


- We use the following input file as an example in this calculation:
       calculation = 'scf',
      restart_mode = 'from_scratch',
        pseudo_dir = 'pseudo/',
            outdir = 'tmp/',
            prefix = 'graphene',
             ibrav = 4,
                 a = 2.4623,
                 c = 10
               nat = 2,
              ntyp = 1,
       occupations = 'smearing',
          smearing = 'methfessel-paxton',
           degauss = 0.02,
           ecutwfc = 60,
           ecutrho = 720,
              nbnd = 8,
          conv_thr = 1.0d-10,
       mixing_mode = 'plain',
       mixing_beta = 0.7,
   diagonalization = 'cg',
 C 12.0107 C.pbe-rrkjus.UPF
 C  0.333333333  0.666666666  0.500000000
 C  0.666666666  0.333333333  0.500000000
 K_POINTS {automatic}
  42 42 1   0 0 0

** The meaning of each command [#command]
- &control:
 declares the control block.

- calculation = 'scf':
 tells PWSCF that this will be a self-consistent field calculation.

- restart_mode = 'from_scratch':
 declares that we will be generation a new structure.

- pseudo_dir = 'pseudo/':
 defines the location of the directory where you store the pseudo-potentials

- outdir='tmp/':
 defines the location of the temporary files. This should always be a local scratch disk so that large I/O operations do not occur across the network.

- prefix='graphene':
 declares the file name prefix to be used for temporary files.

- / :
 denotes the end of a block.

 declares the system block.

- ibrav: 
 gives the crystal system. ibrav=4 : refers to the hexagonal lattice.

- nat: 
 number of atoms (each individual unique atom). Note that graphene for graphene nat=2 because of quantum berry phase. 

- ntyp: 
 number of types of atoms
- ecutwfc:
 Energy cutoff for pseudo-potentials. This one is important; we change this to find stable structure.

- Occupations, Smearing, deguass: 
 these keywords are particular details for the Brillioun zone integration for metals. Since there is a discontinuity of the occupation number for the bands around the Fermi energy, total energy with respect to the number of k-points converges very slowly. Adding electronic temperature (degauss) smooth out the abrupt change of the occupation number and as a result total energy converges with fewer number of k-points.

- nbnd: 
 number of electronic states (bands) to be calculated.Note that in spin-polarized calculations the number of k-point, not the number of bands per k-point, is doubled

- Atomic species declaration.

- After the keyword ATOMIC_SPECIES, for each ntyp enter:
 atomic symbol      atomic weight        pseudo-potential

- Atomic positions:
 after the keyword ATOMIC_POSITIONS, for each nat enter 
 atomic symbol   x      y      z 
 where x,y,z are given as fractional coordinates of the conventional cell.

- k-point selection:
 after the keywork K_POINTS, 'automatic' tells PWSCF to automatically generate a k-point grid.
 The format of the next line is 
 nkx nky nkz offx offy offz 
 where nk* is the number of intervals in a direction and off* is the offset of the origin of the grid.   

- The more details can be found on this website: http://www.quantum-espresso.org/wp-content/uploads/Doc/INPUT_PW.html
**How to Run [#run]
- Create a directory for running the program inside the espresso directory, for example:
 cd ~/espresso
 mkdir graphene
- Under "graphene" directory, put the necessary input files and pseudopotential files needed for calculation. Therefore, at least in the graphene directory we should have
++ an executable file "pw.x" which is the main program of quantum espresso
++ a file "input.in"
++ a directory "pseudo" which consists of pseudopotential files
++ a temporary directory "tmp" for calculation.
- The calculation can be performed by the following command inside graphene directory:
 cd ~/espresso/graphene
 nohup ./pw.x < input.in > output.out &
(we use "nohup" command for running the calculation as a background process)
- To check the calculation is running, we can see the iteration process by following command:
 grep iteration output.out


- Some important output about this calculation:

     bravais-lattice index     =            4
     lattice parameter (alat)  =       4.6090  a.u.
     unit-cell volume          =     347.6569 (a.u.)^3
     number of atoms/cell      =            2
     number of atomic types    =            1
     number of electrons       =         8.00
     number of Kohn-Sham states=            8
     kinetic-energy cutoff     =      40.0000  Ry
     charge density cutoff     =     480.0000  Ry
     Exchange-correlation      =  SLA  PZ   NOGX NOGC ( 1 1 0 0 0)
     EXX-fraction              =        0.00

     Largest allocated arrays     est. size (Mb)     dimensions
        Kohn-Sham Wavefunctions         0.03 Mb     (    264,    8)
        NL pseudopotentials             0.06 Mb     (    264,   16)
        Each V/rho on FFT grid          0.45 Mb     (  29808)
        Each G-vector array             0.08 Mb     (  10283)
        G-vector shells                 0.02 Mb     (   2599)
     Largest temporary arrays     est. size (Mb)     dimensions
        Each subspace H/S matrix        0.00 Mb     (   8,   8)
        Each <psi_i|beta_j> matrix      0.00 Mb     (     16,    8)

** Example silicon  [#si]
- The crystal shape of Silicon is plotted as followig:


- let's make a simple example to reproduce silicon electronic structure. All calculations and results can be found in:


- The structure of si is shown in si.eps which was ploted by xcrysden software.I will explain how to plot at this program later.

    outdir = '/home/pourya/tmp/',
    ibrav=  2, celldm(1) =10.2, nat=  2, ntyp= 1,
    ecutwfc = 12.0, 
 Si  28.086  Si.pz-vbc.UPF
 Si 0.00 0.00 0.00 
 Si 0.25 0.25 0.25 
   0.25 0.25 0.75 3.0
   0.25 0.25 0.25 1.0

- In order to compile this code, you need type following command:
 /pw.x < si.scf.in > si.scf.out
- Here the cutoff energy is consider 12.0 R; but, the total energy versus different cutoff energy is plotted in figure si_cut.pdf in order to show how total energy convergence by different cuttoff value.


- In order to find band at high-symmetry point in silicon please compile the si.nscf.in. The symmetry point is add manually.
 / pw.x < si.nscf.in> si.nscf.out

- Now, the band structure of silicon will be calculated by following command. Please note that the k-point is add manually. It can be calculated in the BZ along high-symmetry points direction as defined by in k-point-path.
 / pw.x < si.bands.in> si.bands.out

- Collect the band results for plotting:
 / bands.x < bands.in > bands.out

- Therefore, there will be 2 files bands.out and bands.dat. for plotting:
 / plotband.x
 input file > bands.dat
 Range: -5.6680 16.4950eV Emin, Emax > -6.0 10.0
 output file (xmgr) > si.bands.xmgr
 output file (ps) > si.bands.ps
 Efermi > 6.337
 deltaE, reference E (for tics) 1.0, 6.337

- hence, you can see the band structure of silicon
 /evince si.bands.ps


- Charge density plot for silicon plane:
 /pp.x < si.pp_rho.in 
 input file > si.rho.dat
 output file > si.rho.ps
 Logarithmic scale (y/n)? > n
 /evince si.rho.ps

- Finally, the general information about electron structure of silicon is produced.

- Some important outputs about this calculation:

     bravais-lattice index     =            2
     lattice parameter (alat)  =      10.2000  a.u.
     unit-cell volume          =     265.3020 (a.u.)^3
     number of atoms/cell      =            2
     number of atomic types    =            1
     number of electrons       =         8.00
     number of Kohn-Sham states=            8
     kinetic-energy cutoff     =      12.0000  Ry
     charge density cutoff     =      48.0000  Ry
     Exchange-correlation      =  SLA  PZ   NOGX NOGC ( 1 1 0 0 0)
     EXX-fraction              =        0.00 

     Largest allocated arrays     est. size (Mb)     dimensions
        Kohn-Sham Wavefunctions         0.02 Mb     (    200,    8)
        NL pseudopotentials             0.02 Mb     (    200,    8)
        Each V/rho on FFT grid          0.05 Mb     (   3375)
        Each G-vector array             0.01 Mb     (   1459)
        G-vector shells                 0.00 Mb     (     43)
     Largest temporary arrays     est. size (Mb)     dimensions
        Auxiliary wavefunctions         0.10 Mb     (    200,   32)
        Each subspace H/S matrix        0.02 Mb     (  32,  32)
        Each <psi_i|beta_j> matrix      0.00 Mb     (      8,    8)

** Charge Density of Silicon  [#si_ph]


- The charge density of silicon can be calculated by QE as following commands the source file can be found in :

- pw.x < si.scf.in > si.scf.out
- pp.x < si.pp.in > si.pp.out

** Phonon  [#si_ph]

- One LA phonon mode and one LO phonon mode is shown at $\Gamma$ point. The vectors indicate of movement direction.


- LO Mode :


Front page   New List of pages Search Recent changes   Help   RSS of recent changes