Skip to content

35: Silicon — Projectability-disentangled Wannier functions with custom projectors

Outline

Obtain MLWFs for silicon using projectability disentanglement, with additional \(3d\) projectors to describe high-energy conduction bands. For more details on the methodology, see Ref.1.

Input files

  • Directory: tutorial/tutorial35/

  • silicon.scf The pw.x input file for ground state calculation

    silicon.scf
    Input file
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
    &CONTROL
  calculation = 'scf'
  etot_conv_thr =   2.0000000000d-05
  forc_conv_thr =   1.0000000000d-04
  outdir = './out/'
  prefix = 'silicon'
  pseudo_dir = '../../pseudo/'
  verbosity = 'high'
/
&SYSTEM
  degauss =   1.0000000000d-02
  ecutwfc =   25
  ibrav = 0
  nat = 2
  nosym = .false.
  ntyp = 1
  occupations = 'smearing'
  smearing = 'cold'
/
&ELECTRONS
  conv_thr =   4.0000000000d-10
  electron_maxstep = 80
  mixing_beta =   4.0000000000d-01
/
ATOMIC_SPECIES
Si     28.085 Si.pbe-n-van.UPF
ATOMIC_POSITIONS angstrom
Si           4.0482056439       1.3494018813       4.0482056439
Si           0.0000000000       0.0000000000       0.0000000000
K_POINTS automatic
10 10 10 0 0 0
CELL_PARAMETERS angstrom
      0.0000000000       2.6988037626       2.6988037626
      2.6988037626       0.0000000000       2.6988037626
      2.6988037626       2.6988037626       0.0000000000

  • silicon.bands The pw.x input file for band structure calculation

    silicon.bands
    Input file
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
    &CONTROL
  calculation = 'bands'
  etot_conv_thr =   2.0000000000d-05
  forc_conv_thr =   1.0000000000d-04
  outdir = './out/'
  prefix = 'silicon'
  pseudo_dir = '../../pseudo/'
  verbosity = 'high'
/
&SYSTEM
  degauss =   1.0000000000d-02
  ecutwfc =   25
  ibrav = 0
  nat = 2
  nbnd = 40
  nosym = .false.
  ntyp = 1
  occupations = 'smearing'
  smearing = 'cold'
/
&ELECTRONS
  conv_thr =   4.0000000000d-10
  diago_full_acc = .true.
  diagonalization = 'cg'
  electron_maxstep = 80
  mixing_beta =   4.0000000000d-01
/
ATOMIC_SPECIES
Si     28.085 Si.pbe-n-van.UPF
ATOMIC_POSITIONS angstrom
Si           4.0482056439       1.3494018813       4.0482056439
Si           0.0000000000       0.0000000000       0.0000000000
CELL_PARAMETERS angstrom
      0.0000000000       2.6988037626       2.6988037626
      2.6988037626       0.0000000000       2.6988037626
      2.6988037626       2.6988037626       0.0000000000
K_POINTS crystal
         451
    0.000000    0.000000    0.000000   1.0
    0.005000    0.000000    0.005000   1.0
    0.010000    0.000000    0.010000   1.0
    0.015000    0.000000    0.015000   1.0
    0.020000    0.000000    0.020000   1.0
    0.025000    0.000000    0.025000   1.0
    0.030000    0.000000    0.030000   1.0
...
...

  • silicon.bandsx The bands.x input file for extracting band structure eigenvalues

    silicon.bandsx
    Input file
    1
2
3
4
5
6
    &bands
  prefix = 'silicon',
  outdir = './out/',
  lsym = .false.,
  filband = 'silicon.bands.dat',
/

  • silicon.nscf The pw.x input file to obtain Bloch states on a uniform grid

    silicon.nscf
    Input file
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
    &CONTROL
  calculation = 'nscf'
  etot_conv_thr =   2.0000000000d-05
  forc_conv_thr =   1.0000000000d-04
  max_seconds =   4.1040000000d+04
  outdir = './out/'
  prefix = 'silicon'
  pseudo_dir = '../../pseudo/'
  verbosity = 'high'
/
&SYSTEM
  degauss =   1.0000000000d-02
  ecutwfc =   25
  ibrav = 0
  nat = 2
  nbnd = 40
  noinv = .true.
  nosym = .true.
  ntyp = 1
  occupations = 'smearing'
  smearing = 'cold'
/
&ELECTRONS
  conv_thr =   4.0000000000d-10
  diago_full_acc = .true.
  electron_maxstep = 80
  mixing_beta =   4.0000000000d-01
  startingpot = 'file'
/
ATOMIC_SPECIES
Si     28.085 Si.pbe-n-van.UPF
ATOMIC_POSITIONS angstrom
Si           4.0482056439       1.3494018813       4.0482056439
Si           0.0000000000       0.0000000000       0.0000000000
CELL_PARAMETERS angstrom
      0.0000000000       2.6988037626       2.6988037626
      2.6988037626       0.0000000000       2.6988037626
      2.6988037626       2.6988037626       0.0000000000
K_POINTS crystal
1000
  0.00000000  0.00000000  0.00000000  1.000000e-03
  0.00000000  0.00000000  0.10000000  1.000000e-03
  0.00000000  0.00000000  0.20000000  1.000000e-03
  0.00000000  0.00000000  0.30000000  1.000000e-03
  0.00000000  0.00000000  0.40000000  1.000000e-03
...
...

  • silicon.pw2wan Input file for pw2wannier90.x

    silicon.pw2wan
    Input file
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
    &INPUTPP
  atom_proj = .true.
  ! use external projectors
  atom_proj_ext = .true.
  ! dir for external projectors
  atom_proj_dir = './ext_proj'
  outdir = './out/'
  prefix = 'silicon'
  seedname = 'silicon'
  ! for excluding specific projectors
  ! this excludes 3d projectors, then the results are similar
  ! to that of using UPF file, i.e., project onto Si s+p orbitals
  ! for the indices of orbitals, see the pw2wan stdout
  ! atom_proj_exclude = 5 6 7 8 9 14 15 16 17 18
/

  • silicon.win The wannier90.x input file

    silicon.win
    Input file
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
    num_wann = 18
num_bands = 40

auto_projections = .true.

dis_froz_proj = .true.
dis_proj_max =   0.95
dis_proj_min =   0.01

! you can also add an energy window to be safe
dis_froz_max = 8.7

fermi_energy =   6.7332

num_iter = 4000
conv_tol =   2.0000000000d-10
conv_window = 3
dis_conv_tol =   2.0000000000d-10
dis_num_iter = 4000
num_cg_steps = 200

mp_grid = 10 10 10

bands_plot = .true.
begin atoms_cart
ang
Si         4.0482056439       1.3494018813       4.0482056439
Si         0.0000000000       0.0000000000       0.0000000000
end atoms_cart

begin unit_cell_cart
ang
      0.0000000000       2.6988037626       2.6988037626
      2.6988037626       0.0000000000       2.6988037626
      2.6988037626       2.6988037626       0.0000000000
end unit_cell_cart

begin kpoint_path
G  0.0000000000  0.0000000000  0.0000000000  X  0.5000000000  0.0000000000  0.5000000000
X  0.5000000000  0.0000000000  0.5000000000  U  0.6250000000  0.2500000000  0.6250000000
K  0.3750000000  0.3750000000  0.7500000000  G  0.0000000000  0.0000000000  0.0000000000
G  0.0000000000  0.0000000000  0.0000000000  L  0.5000000000  0.5000000000  0.5000000000
L  0.5000000000  0.5000000000  0.5000000000  W  0.5000000000  0.2500000000  0.7500000000
W  0.5000000000  0.2500000000  0.7500000000  X  0.5000000000  0.0000000000  0.5000000000
end kpoint_path

begin kpoints
  0.00000000  0.00000000  0.00000000
  0.00000000  0.00000000  0.10000000
  0.00000000  0.00000000  0.20000000
...
...

  • silicon_bandsdiff.gnu The gnuplot script to compare DFT and Wannier bands

    silicon_bandsdiff.gnu
    Gnuplot script
     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
    #!/usr/bin/env gnuplot
set terminal pdf enhanced color dashed lw 1 size 6in,6in
set output "silicon_bandsdiff.pdf"
set size ratio 1
fermi = 6.7332
# set style data dots
set key
set xrange [0: 5.22358]
set yrange [ (-6.69831 - fermi) : (34.57817 - fermi)]
set arrow from  1.16407,  (-6.69831 - fermi) to  1.16407,  (34.57817 - fermi) nohead
set arrow from  1.57563,  (-6.69831 - fermi) to  1.57563,  (34.57817 - fermi) nohead
set arrow from  2.81031,  (-6.69831 - fermi) to  2.81031,  (34.57817 - fermi) nohead
set arrow from  3.81842,  (-6.69831 - fermi) to  3.81842,  (34.57817 - fermi) nohead
set arrow from  4.64154,  (-6.69831 - fermi) to  4.64154,  (34.57817 - fermi) nohead
set xtics ("G"  0.00000,  "X"  1.16407,  "U|K"  1.57563,  "G"  2.81031, \
           "L"  3.81842,  "W"  4.64154,  "X"  5.22358)
# scale QE x-axis to be consistent with w90
plot "silicon.bands.dat.gnu" u ($1*1.6459):($2 - fermi) w l title "QE", \
     "silicon_band.dat" u 1:($2 - fermi) w l title "W90"

Steps

  1. Run pw.x to obtain the ground state of silicon

    Terminal
    pw.x < silicon.scf > scf.out

  2. Run pw.x to obtain the band structure of silicon

    Terminal
    pw.x < silicon.bands > bands.out

  3. Run bands.x to obtain a silicon.bands.dat file containing the band structure of silicon

    Terminal
    bands.x < silicon.bandsx > bandsx.out

  4. Run pw.x to obtain the Bloch states on a uniform k-point grid

    Terminal
    pw.x < silicon.nscf > nscf.out

  5. Run pw.x to generate a list of the required overlaps (written into the silicon.nnkp file).

    Note

    See win input file, no need to specify initial projections, they are chosen from the pseudo-atomic orbitals inside the ext_proj/Si.dat file.

    Terminal
    wannier90.x -pp silicon

  6. Run pw2wannier90.x to compute the overlap between Bloch states and the projections for the starting guess (written in the silicon.mmn and silicon.amn files).

    Terminal
    pw2wannier90.x < silicon.pw2wan > pw2wan.out

  7. Run pw.x to compute the MLWFs.

    Terminal
    wannier90.x silicon

  8. Run gnuplot to compare DFT and Wannier-interpolated bands, this will generate a PDF file silicon_bandsdiff.pdf, see Fig.Bands comparison.

    Terminal
    ./silicon_bandsdiff.gnu

    Bands diff
    Comparison of DFT and Wannier bands for silicon.

  9. (Optional) Clean up all output files

    Terminal
    make clean

Further ideas

  1. Try changing the atom_proj_exclude in silicon.pw2wan file, i.e., these commented lines

    Input file
    10
11
12
13
14
      ! for excluding specific projectors
  ! this excludes 3d projectors, then the results are similar
  ! to that of using UPF file, i.e., project onto Si s+p orbitals
  ! for the indices of orbitals, see the pw2wan stdout
  ! atom_proj_exclude = 5 6 7 8 9 14 15 16 17 18

    Hint

    You need to set num_wann = 8 in the silicon.win file as well, to be consistent with the reduced number of projectors.

  2. Now that \(3d\) projectors provide us a larger space for optimization, you can try increasing the dis_froz_max to freeze higher energy bands, if you are targeting at reproducing those eigenvalues.

    Note

    The dis_proj_min/max and dis_froz_min/max can be enabled simultaneously: the union of inner energy window and high-projectability states will be freezed, and the union of states outside outer energy window and having low projectability will be discarded. Thus, you can still use energy window to make sure near-Fermi energy states are well reproduced, and use "projectability window" to selectively freeze atomic-like states in the conduction region.

  3. The default dis_proj_max = 0.95 might not freeze all the states you want, try changing this value and see the band interpolation results. For other materials, it might worth trying decreasing this value to freeze more states.

  1. Junfeng Qiao, Giovanni Pizzi, and Nicola Marzari. Projectability disentanglement for accurate and automated electronic-structure Hamiltonians. npj Comput. Mater., 9(1):208, Nov 2023. doi:10.1038/s41524-023-01146-w. ↩