VaspBandUnfolding

Introduction

VaspBandUnfolding consists of a collection of python scripts that deal with VASP output files.

• vaspwfc.py can be used to read the plane-wave coefficients in the WAVECAR file and then generate real-space representation of the pseudo-wavefunction by Fourier transform. Other wavefunction related quantites, e.g. transition dipole moment (pseudo-wavefunction part), inverse participation ratio and electron localization function et al. can also be conveniently obtained.

  Moreover, a command line tool named wfcplot in the bin directory can be used to output real-space pseudo-wavefunctions.

• paw.py can be used to parse the PAW POTCAR file. A command line tool named potplot in the bin directory can be used to visualize the projector function and partial waves contained in the POTCAR.

• aewfc.py can be used to generate all-electron (AE) wavefunction.

• unfold.py can be used to perform band unfolding from supercell calculations.

Publications

A list of publications utilizing VaspBandUnfolding can be found here!

Installation

• Requirements

  • Numpy
  • Scipy
  • Matplotlib
  • ASE
  • pySBT for spherical Bessel transform

• Manual Installation

  git clone https://github.com/QijingZheng/VaspBandUnfolding
  # set the 
  cd VaspBandUnfolding
  python setup.py install --prefix=/the/path/of/your/dir/
  export PYTHONPATH=/the/path/of/your/dir:${PYTHONPATH}
  

• Using Pip

  pip install git+https://github.com/QijingZheng/VaspBandUnfolding
  

Examples

Reading VASP WAVECAR

• Pseudo-wavefunction.

  As is well known, VASP WAVECAR is a binary file and contains the plane-wave coefficients for the pseudo-wavefunctions. The pseudo-wavefunction in real space can be obtained by 3D Fourier transform on the plane-wave coefficients and represented on a 3D uniform grid which can be subsequently visualized with software such as VESTA.

  • For a normal WAVECAR, i.e. not gamma-only or non-collinear WAVECAR, one can write a small python script and convert the desired Kohn-Sham states to real space.

    #/usr/bin/env python
    from vaspwfc import vaspwfc
    
    pswfc = vaspwfc('./examples/wfc_r/WAVECAR')
    # KS orbital in real space, double the size of the FT grid
    phi = pswfc.get_ps_wfc(ikpt=2, iband=27, ngrid=pswfc._ngrid * 2)
    # Save the orbital into files. Since the wavefunction consist of complex
    # numbers, the real and imaginary part are saved separately.
    pswfc.save2vesta(phi, poscar='./examples/wfc_r/POSCAR')
    

    • In the above script, pswfc._ngrid is the default 3D grid size and phi is a numpy 3D array of size 2*pswfc._ngrid, with the first dimensiton being x and the last "z".
    • The spin, k-point and band index for the KS state are designated by the argumnt ispin, ikpt and iband, respectively, all of which start from 1.
    • Generally, the pseudo-wavefunction is complex, pswfc.save2vesta will export both the real and imaginary part of the wavefunction, with the file name "wfc_r.vasp" and "wfc_i.vasp", respectively.

    Below are the real (left) and imaginary (right) part of the selected KS orbital:

    

  • For gamma-only WAVECAR, one must pass the argument lgamma=True when reading WAVECAR in the vaspwfc method. Moreover, as VASP only stores half of the full plane-wave coefficients for gamma-only WAVECAR and VASP changes the idea about which half to save from version 5.2 to 5.4. An addition argument must be passed.

    #/usr/bin/env python
    from vaspwfc import vaspwfc
    
    # For VASP <= 5.2.x, check
    # which FFT VASP uses by the following command:
    #
    #     $ grep 'use.* FFT for wave' OUTCAR
    #
    # Then
    #
    #     # for parallel FFT, VASP <= 5.2.x
    #     pswfc = vaspwfc('WAVECAR', lgamma=True, gamma_half='z')
    #
    #     # for serial FFT, VASP <= 5.2.x
    #     pswfc = vaspwfc('WAVECAR', lgamma=True, gamma_half='x')
    #
    # For VASP >= 5.4, WAVECAR is written with x-direction half grid regardless of
    # parallel or serial FFT.
    #
    #     # "gamma_half" default to "x" for VASP >= 5.4
    #     pswfc = vaspwfc('WAVECAR', lgamma=True, gamma_half='x')
    
    pswfc = vaspwfc('WAVECAR', lgamma=True, gamma_half='x')
    

  • For non-collinear WAVECAR, however, one must pass the argument lsorbit=True when reading WAVECAR. Note that in the non-collinear case, the wavefunction now is a two-component spinor.

    #/usr/bin/env python
    from vaspwfc import vaspwfc
    
    # for WAVECAR from a noncollinear run, the wavefunction at each k-piont/band is
    # a two component spinor. Turn on the lsorbit flag when reading WAVECAr.
    pswfc = vaspwfc('examples/wfc_r/wavecar_mose2-wse2', lsorbit=True)
    phi_spinor = pswfc.get_ps_wfc(1, 1, 36, ngrid=pswfc._ngrid*2)
    for ii in range(2):
        phi = phi_spinor[ii]
        prefix = 'spinor_{:02d}'.format(ii)
        pswfc.save2vesta(phi, prefix=prefix,
                poscar='examples/wfc_r/poscar_mose2-wse2')
    

  • If only real-space representation of the pseudo-wavefunction is needed, a helping script wfcplot in the bin directory comes to rescue.

    $ wfcplot -w WAVECAR -p POSCAR -s spin_index -k kpoint_index -n band_index             # for normal WAVECAR
    $ wfcplot -w WAVECAR -p POSCAR -s spin_index -k kpoint_index -n band_index  -lgamma    # for gamma-only WAVECAR
    $ wfcplot -w WAVECAR -p POSCAR -s spin_index -k kpoint_index -n band_index  -lsorbit   # for noncollinear WAVECAR
    

    Please refer to wfcplot -h for more information of the usage.

• All-electron wavefunction in real space

  Refer to this post for detail formulation.

  PAW All-Electron Wavefunction in VASP

  #/usr/bin/env python
  
  from vaspwfc import vaspwfc
  from aewfc import vasp_ae_wfc
  
  # the pseudo-wavefunction
  ps_wfc = vaspwfc('WAVECAR', lgamma=True)
  # the all-electron wavefunction
  # here 25x Encut, or 5x grid size is used
  ae_wfc = vasp_ae_wfc(ps_wfc, aecut=-25)
  
  phi_ae = ae_wfc.get_ae_wfc(iband=8)
  

  The comparison of All-electron and pseudo wavefunction of CO₂ HOMO can be found in examples/aewfc/co2.

  

• Dipole transition matrix

  Refer to this post for detail formulation.

  Light-Matter Interaction and Dipole Transition Matrix

  Under the electric-dipole approximation (EDA), The dipole transition matrix elements in the length gauge is given by:

        <psi_nk | e r | psi_mk>
  

  where | psi_nk > is the pseudo-wavefunction. In periodic systems, the position operator "r" is not well-defined. Therefore, we first evaluate the momentum operator matrix in the velocity gauge, i.e.

        <psi_nk | p | psi_mk>
  

  And then use simple "p-r" relation to apprimate the dipole transition matrix element

                                    -i⋅h
      <psi_nk | r | psi_mk> =  -------------- ⋅ <psi_nk | p | psi_mk>
                                 m⋅(En - Em)
  

  Apparently, the above equaiton is not valid for the case Em == En. In this case, we just set the dipole matrix element to be 0.

  NOTE that, the simple "p-r" relation only applies to molecular or finite system, and there might be problem in directly using it for periodic system. Please refer to this paper for more details.

  Relation between the interband dipole and momentum matrix elements in semiconductors

  The momentum operator matrix in the velocity gauge

          <psi_nk | p | psi_mk> = hbar <u_nk | k - i nabla | u_mk>
  

  In PAW, the matrix element can be divided into plane-wave parts and one-center parts, i.e.

      <u_nk | k - i nabla | u_mk> = <tilde_u_nk | k - i nabla | tilde_u_mk>
                                   - \sum_ij <tilde_u_nk | p_i><p_j | tilde_u_mk>
                                     \times i [
                                       <phi_i | nabla | phi_j>
                                       -
                                       <tilde_phi_i | nabla | tilde_phi_j>
                                     ]
  

  where | u_nk > and | tilde_u_nk > are cell-periodic part of the AE/PS wavefunctions, | p_j > is the PAW projector function and | phi_j > and | tilde_phi_j > are PAW AE/PS partial waves.

  The nabla operator matrix elements between the pseudo-wavefuncitons

      <tilde_u_nk | k - i nabla | tilde_u_mk>
  
     = \sum_G C_nk(G).conj() * C_mk(G) * [k + G]
  

  where C_nk(G) is the plane-wave coefficients for | u_nk >.

  import numpy as np
  
  from vaspwfc import vaspwfc
  from aewfc import vasp_ae_wfc
  
  # the pseudo-wavefunction
  ps_wfc = vaspwfc('WAVECAR', lgamma=True)
  # the all-electron wavefunction
  ae_wfc = vasp_ae_wfc(ps_wfc)
  
  # (ispin, ikpt, iband) for initial and final states
  ps_dp_mat = ps_wfc.get_dipole_mat((1,1,1), (1, 1, 9))
  ae_dp_mat = ae_wfc.get_dipole_mat((1,1,1), (1, 1, 9))
  

• Inverse Participation Ratio

  IPR is a measure of the localization of Kohn-Sham states. For a particular KS state \phi_j, it is defined as

                  \sum_n |\phi_j(n)|^4 
  IPR(\phi_j) = -------------------------
                |\sum_n |\phi_j(n)|^2||^2
  

  where n iters over the number of grid points.

• Electron Loca