QuantumEspresso(Open)

[started by Pourya 2014.12.3]

Introduction†

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:

QE-Tutorial.pdf

↑

Prerequisites†

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): http://qe-forge.org/gf/download/frsrelease/116/403/espresso-5.0.2.tar.gz
If you use quantum espresso(v.5.0.2), please download the patch (e.g. as "espresso-5.0.2-5.0.3.diff"): http://qe-forge.org/gf/project/q-e/frs/?action=FrsReleaseView&release_id=128

↑

How to 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†

The structure of the input file of quantum espresso is

&CONTROL 
/
&SYSTEM
/
&ELECTRONS
/
[ &IONS
 ...
/ ]
[ &CELL
 ...
/ ]
ATOMIC_SPECIES
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)
  nks
  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) ]
[ OCCUPATIONS
  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) ] ]
[ CONSTRAINTS
  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:†

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 ?:†

http://www.quantum-espresso.org/pseudopotentials/

↑

Where can we find the Brillouin zone and High symmetry points?†

Bilbao Crystallographic Server: http://www.cryst.ehu.es/cryst/get_kvec.html

↑

Some Conversions:†

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†

A single layer graphene is shown as following figure:

The electron structure of graphene is : &tex(): Error! cURL: Could not resolve host: chart.apis.google.com; &tex(): Error! cURL: Could not resolve host: chart.apis.google.com; &tex(): Error! cURL: Could not resolve host: chart.apis.google.com; and its radial wave function for each orbital is :

We use the following input file as an example in this calculation:

&CONTROL
      calculation = 'scf',
     restart_mode = 'from_scratch',
       pseudo_dir = 'pseudo/',
           outdir = 'tmp/',
           prefix = 'graphene',
/
&SYSTEM
            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,
/
&ELECTRONS
         conv_thr = 1.0d-10,
      mixing_mode = 'plain',
      mixing_beta = 0.7,
  diagonalization = 'cg',
/
ATOMIC_SPECIES
C 12.0107 C.pbe-rrkjus.UPF
ATOMIC_POSITIONS {crystal}
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†

&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.
```

&SYSTEM:
```
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†

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
1. an executable file "pw.x" which is the main program of quantum espresso
2. a file "input.in"
3. a directory "pseudo" which consists of pseudopotential files
4. 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†

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:

/home/students/pourya/espresso/silicon

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.

&control
   prefix='silicon',
   pseudo_dir='/home/pourya/Desktop/espresso/pseudo/'
   outdir = '/home/pourya/tmp/',
/
&system    
   ibrav=  2, celldm(1) =10.2, nat=  2, ntyp= 1,
   ecutwfc = 12.0, 
/
&electrons
/

ATOMIC_SPECIES

Si  28.086  Si.pz-vbc.UPF

ATOMIC_POSITIONS

Si 0.00 0.00 0.00 
Si 0.25 0.25 0.25

K_POINTS

  2
  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 
/plotrho.x
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†

#ref(): File not found: "si.charge001.png" at page "QuantumEspresso(Open)"

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†

One LA phonon mode and one LO phonon mode is shown at &tex(): Error! cURL: Could not resolve host: chart.apis.google.com; point. The vectors indicate of movement direction.