Tutorials — VASPKIT 1.5 documentation (2024)

Ways to Run vaspkit¶

There are at lease five ways to run vaspkit under either the interactive user interface or command line mode. We take the generation of KPOINTS (task 102) as an example to illustrate the useage of vaspkit.

1. Just type vaspkit in the terminal to launch the interactive user interface mode;

2. vaspkit -task 102 -kpr 0.04 to generate KPOINTS file with a reciprocal space resolution of \(2\pi \times 0.04\) \(Å^{-1}\). More details can be get by run vaspkit -help. Note that part of the functions have not been implemented yet!;

3. echo -e "102\n2\n0.04\n"| vaspkit;

4. (echo 102; echo 2; echo 0.04)|vaspkit;

5. vi cmd.in (any file name if you want) including the following content:

10220.04

and then run vaspkit <cmd.in.

Run vaspkit in Batch Mode¶

if one wants to run vaspkit in batch mode, for example, to generate KPOINTS file in several sub-olders:

for i in `ls`do echo ${i} cd ${i} vaspkit -task 102 -file POSCAR -kpr 0.04 cd ..done

Or use echo -e command to input for VASPKIT:

echo -e "102\n2\n0.04" | vaspkit

means input 102, 2, 0.04 in turn in VASPKIT.

if one wants to extract the specified physical quantity, for example, to read band gap in several sub-olders:

for i in `ls`do echo ${i} cd ${i} echo -e "911\n" | vaspkit | sed -n '40p' | awk '{print $4}' cd ..done

Generate INCAR¶

Run VASPKIT in the directory containing POSCAR. Enter 1 to selectthe function VASP Input Files Generator, and then enter 101toselect customize INCAR File, you will get the following displayinformation:

101 +-------------------------- Warm Tips --------------------------+ You MUST Know What You Are Doing Some Parameters in INCAR File Neet To Be Set/Adjusted Manually +---------------------------------------------------------------+ ======================== INCAR Options ========================== ST) Static-Calculation SR) Standard Relaxation MG) Magnetic Properties SO) Spin-Orbit Coupling D3) DFT-D3 no-damping Correction H6) HSE06 Calculation PU) DFT+U Calculation MD) Molecular Dynamics GW) GW0 Calculation BS) BSE Calculation DC) Elastic Constant EL) ELF Calculation BD) Bader Charge Analysis OP) Optical Properties EC) Static Dielectric Constant PC) Decomposed Charge Density FD) Phonon-Finite-Displacement DT) Phonon-DFPT NE) Nudged Elastic Band (NEB) DM) The Dimer Method FQ) Frequence Calculations LR) Lattice Relaxation 0) Quit 9) Back ------------>> Input Key-Parameters (STH6D3 means HSE06-D3 Static-Calcualtion)

Enter the words for specific task. The generated INCAR file will containcorresponding keywords that are required for this task.

For example, to do a single-point calculation (ST) with hybridfunctional HSE06 (H6) and DFT-D3 (D3) vdW correction, enterSTH6D3.

If enter LR, one will get a INCAR for lattice relaxation task withdetail comments:

Global Parameters ISTART = 1 (Read existing wavefunction; if there) # ISPIN = 2 (Spin polarised DFT) # ICHARG = 11 (Non-self-consistent: GGA/LDA band structures) LREAL = Auto (Projection operators: automatic) # ENCUT = 400 (Cut-off energy for plane wave basis set, in eV) PREC = Normal (Precision level) LWAVE = .TRUE. (Write WAVECAR or not) LCHARG = .TRUE. (Write CHGCAR or not) ADDGRID= .TRUE. (Increase grid; helps GGA convergence) # LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not) # LVHAR = .TRUE. (Write ionic+Hartree electrostatic potential into LOCPOT or not) # NELECT = (No. of electrons: charged cells; be careful) # LPLANE = .TRUE. (Real space distribution; supercells) # NPAR = 4 (Max is no. nodes; don't set for hybrids) # NWRITE = 2 (Medium-level output) # KPAR = 2 (Divides k-grid into separate groups) # NGX = 500 (FFT grid mesh density for nice charge/potential plots) # NGY = 500 (FFT grid mesh density for nice charge/potential plots) # NGZ = 500 (FFT grid mesh density for nice charge/potential plots)Lattice Relaxation NSW = 300 (number of ionic steps) ISMEAR = 0 (gaussian smearing method ) SIGMA = 0.05 (please check the width of the smearing) IBRION = 2 (Algorithm: 0-MD; 1-Quasi-New; 2-CG) ISIF = 3 (optimize atomic coordinates and lattice parameters) EDIFFG = -1.5E-02 (Ionic convergence; eV/AA) PREC = Accurate (Precision level)

If one want a cleaner INCAR without comments, set the ~/.vaspkitfile with:

MINI_INCAR .TRUE.

The automatic generated INCAR for lattice relaxation task will be:

Global Parameters ISTART = 1 LREAL = Auto PREC = Normal LWAVE = .TRUE. LCHARG = .TRUE. ADDGRID= .TRUE.Lattice Relaxation NSW = 300 ISMEAR = 0 SIGMA = 0.05 IBRION = 2 ISIF = 3 EDIFFG = -1.5E-02 PREC = Accurate

Users can also modify some parameters as they want from these INCARfile.

If there is already an INCAR file, the original INCAR will beoverwritten as default. Edit the ~/.vaspkit to change the INCARoutput settings. Simply change the SET_INCAR_WRITE_MODE line.

SET_INCAR_WRITE_MODE OVERRIDE # OVERRIDE, APPEND,BACK-UP-OLD,BACK-UP-NEW;

Generate KPOINTS¶

For self-consistent calculation, users need to prepare a KPOINTS file tospecify density of K-points and the method for automatic k-meshgeneration.

VASPKIT can generate KPOINTS file automatically with pre-exist POSCARfile. Run VASPKIT option 1) VASP Input Files Generator, then enteroption 102 to generate KPOINTS file for SCF calculation. Then inputthe K-mesh scheme according to the following display information:

102 ======================= K-Mesh Scheme ========================== 1) Monkhorst-Pack Scheme 2) Gamma Scheme 0) Quit 9) Back ------------->>

Enter 1 to select the original Monkhorst-Pack scheme,

Enter 2 to select the Gamma centered Monkhorst-Pack scheme.

Then, VASPKIT will ask us to input KPT-Resolved Value between K-pointsin reciprocal cell in units of \(2\pi \times 0.04 Å^{-1}\).

2 -->> (01) Reading Structural Parameters from POSCAR File... +-------------------------- Warm Tips --------------------------+ * Accuracy Levels: Gamma-Only: 0; Low: 0.06~0.04; Medium: 0.04~0.03; Fine: 0.02-0.01. * 0.03-0.04 is Generally Precise Enough! +---------------------------------------------------------------+ Input KPT-Resolved Value (e.g., 0.04, in unit of 2*PI/Angstrom): ------------>>

Number of K-points increases when the KPT-resolved value (kpr) decreases.Number of K-points decreases when the kpr value increase. Foreach direction, the number is determined by

\[N=\max \left(1,\left|\vec{b}_{i}\right| / \mathrm{kpr}\right)\]

where \(\vec{b}_{i}\) are the reciprocal lattice vectors. Thesevalues are rounded to the next integer greater than or equal to N.The recommend value ~0.04 (\(2\pi \times 0.04 Å^{-1}\)) is enough for mostsystem. This parameter is similar as the parameter KSPACING in INCAR.But the unit is different. Unit of KSPACING is \(Å^{-1}\), and the unit ofVASPKIT is \(2\pi \times 0.04 Å^{-1}\).

The first line of output KPOINTS file show the users’ definedKPT-Resolved Value.

K-Mesh Generated with KP-Resolved Value (...): 0.0200Gamma 14 14 140.0 0.0 0.0

Generate POTCAR¶

When generating KPOINTS, POTCAR will also be generate automatically. Orrun VASPKIT 103 to generate POTCAR.

103 -->> (01) Reading Structural Parameters from POSCAR File... -->> (02) Written POTCAR File with the Standard Potential!

It reads the elements information from POSCAR and combines thecorresponding POTCAR from the pseudo potential folders that you set in~/.vaspkit.

GGA_PATH '~/POTCAR/GGA' # Path of GGA potential.PBE_PATH '~/POTCAR/PBE' # Path of PBE potential.LDA_PATH '~/POTCAR/LDA' # Path of LDA potential.POTCAR_TYPE PBE # PBE, GGA or LDA;RECOMMENDED_POTCAR .TRUE. # .TRUE. or .FALSE.;

Set POTCAR_TYPE to PBE, GGA, orLDA as you want. TheRECOMMENDED_POTCAR tag control whether to use recommended potentialsfrom VASP manual (Page. 195, 2018.10.29,http://cms.mpi.univie.ac.at/vasp/vasp/Recommended_PAW_potentials_DFT_calculations_using_vasp_5_2.html).If RECOMMENDED_POTCAR is .FALSE., POTCAR with no extensions willbe used. If RECOMMENDED_POTCAR is .TRUE., official recommendedPOTCAR will be used.

POTCAR types:

1. No extensions “_”

2. _d. An extension d, treat the d semi core states as valencestate.

3. _pv or _sv. The extensions _pv and _sv imply that thep and s semi-core states are treated as valence states.

4. _h and _s. An extension _h or_s implies that thepotential is harder or softer than the standard potential and hencerequires a higher or lower energy cutoff.

5. Pseudo hydrogen. ex: H.5

6. _GW. Used for GW calculation.

If one want to generate POTCAR for GW calculation, Set the GW_POTCARto .TRUE. in ~/.vaspkit.

GW_POTCAR .TRUE. # .TRUE. or .FALSE.;

VASPKIT also provide 104 option to generate POTCAR manually byselecting the type of potential for each element.

104 -->> (1) Reading Structural Parameters from POSCAR File... Auto detected POTCAR_TYPE is O, please type the one you want!O_h Auto detected POTCAR_TYPE is Ti, please type the one you want!Ti_sv -->> (2) Written POTCAR File with user specified Potential!

In this example, enter the customized type of potentials O and Ti.Ti_sv was chosen for Ti and O_h was chosen for O. If there is nouser defined potential type in the list of potentials, VASPKIT will askuser re-enter.

Generate POSCAR¶

VASPKIT can transform .cif and .xsd (Materials Studio format) files toPOSCAR format by option 105 and 106.

105 will call /vaspkit.1.00/utilities/cif2pos.py script.

105 Please type in the filename of cif->al2o3.cifPleas input the order of element, `ENTER` for default!Example: 'Al O' in this CIF -->> (01) POSCAR has been generated...

106 will call /vaspkit.1.00/utilities/xsd2pos.py scriptautomatically. Note that the atom fix information in .xsd file is keptwhen transform.

106 Build->Symmetry->Make P1,then select atoms to be fixed, Modify->Constraints->fi x fractional position-> Please type in the filename of xsd->CONTCAR-n2-3.xsd -->> (01) POSCAR has been generated...

107 can reorder the elements in POSCAR file.

107 -->> (01) Reading Structural Parameters from POSCAR File... Please Type the New Order of Elements to Sort. (Tip: The Initial Order of Elements in POSCAR File is: Al O) ------------->>O Al -->> (02) Written POSCAR_REV File!

Input File Self Check¶

VASPKIT can do format correction and pseudo potential check by option109. VASPKIT will automatically correct the INCAR and POSCAR formatsand check if the POTCAR and POSCAR are consistent.

Pre-process Band Structure (pure functional)¶

To do band structure calculation, one need to prepare a primitive celland corresponding K points path (K-path) alone Irreducible BrillouinZone. Irreducible Brillouin Zone is the first Brillouin zone reduced byall of the symmetries in the points group of the lattice (point group ofthe crystal). Recognize and select high symmetry points, and link themalong edges of Irreducible Brillouin Zone.

For example, conventional FCC metal cell, Irreducible Brillouin Zone,and high symmetry points:

Conventional BCC metal cell, Irreducible Brillouin Zone, and highsymmetry points:

K-path is not unique. Usually, It is unnecessary to select all linesbetween all the high symmetry points. Representative and important lineare selected and written as line-mode in KPOINTS file. For large-scaleand high-throughout calculation, there should be a rule to define thepath from structural information. pymatgen and seeK-path provide somesolutions but only can be used for 3D system. VASPKIT provide a tool togenerate K-path for 1D (task 301), 2D (task 302), and 3D (task 303) materials based on a systematic rule:

Here are Brillouin zones for 2D materials (V. Wang, Y.-Y. Liang, Y. Kawazeo, W.-T. Geng, High-Throughput Computational Screening of Two-Dimensional Semiconductors, arXiv:1806.04285.)

Other ways to get K-path automatically by using pymatgen（https://pymatgen.org/）Computational Materials Science 49(2010) 299–312.seek-path（https://www.materialscloud.org/work/tools/seekpath）ComputationalMaterials Science 128 (2017) 140–184.

For the suggested k-paths of bulk materials, VASPKIT uses the same algorithm as seek-path website (Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba, I. Tanaka, Band structurediagram paths based on crystallography, Comp. Mat. Sci. 128, 140 (2017).).

MoS21.0 3.1659998894 0.0000000000 0.0000000000 -1.5829999447 2.7418363326 0.0000000000 0.0000000000 0.0000000000 18.4099998474 S Mo 2 1Direct 0.000000000 0.000000000 0.413899988 0.000000000 0.000000000 0.586099982 0.666666687 0.333333343 0.500000000

92 +-------------------------- Warm Tips --------------------------+ Please Use These Features with CAUTION! +---------------------------------------------------------------+ ===================== 2D Materials Toolkit ====================== 921) Center Aomic-Layer along z direction 922) Resize Vacuum Thickness 923) Standardize 2D Crystal Cell 926) Elastic Constants for 2D Materials 927) Valence and Conduction Band Edges Referenced to Vacuum Level 929) Summary for Relaxed 2D Structure 0) Quit 9) Back ------------>>923 -->> (1) Reading Structural Parameters from POSCAR File... -->> (2) Written POSCAR_NEW File!

302 +-------------------------- Warm Tips --------------------------+ See An Example in vaspkit/examples/seek_kpath/graphene_2D. This feature is still experimental & check the PRIMCELL.vasp file. +---------------------------------------------------------------+ -->> (1) Reading Structural Parameters from POSCAR File... +-------------------------- Summary ----------------------------+ The vacuum slab is supposed to be along c axis Prototype: AB2 Total Atoms in Input Cell: 3 Lattice Constants in Input Cell: 3.166 3.166 18.410 Lattice Angles in Input Cell: 90.000 90.000 120.000 Total Atoms in Primitive Cell: 3 Lattice Constants in Primitive Cell: 3.166 3.166 18.410 Lattice Angles in Primitive Cell: 90.000 90.000 120.000 2D Bravais Lattice: Hexagonal Space Group: 187 Point Group: 26 [ D3h ] International: P-6m2 Symmetry Operations: 12 Suggested K-Path: (shown in the next line) [ GAMMA-M-K-GAMMA ] +---------------------------------------------------------------+ -->> (2) Written PRIMCELL.vasp file. -->> (3) Written KPATH.in File for Band-Structure Calculation. -->> (4) Written HIGH_SYMMETRY_POINTS File for Reference.

K-Path Generated by VASPKIT 20Line-ModeReciprocal 0.0000000000 0.0000000000 0.0000000000 GAMMA 0.5000000000 0.0000000000 0.0000000000 M 0.5000000000 0.0000000000 0.0000000000 M 0.3333333333 0.3333333333 0.0000000000 K 0.3333333333 0.3333333333 0.0000000000 K 0.0000000000 0.0000000000 0.0000000000 GAMMA

High-symmetry points (in fractional coordinates). You can check them with the seekpath database [https://www.materialscloud.org/work/tools/seekpath]. 0.0000000000 0.0000000000 0.0000000000 GAMMA 0.3333333333 0.3333333333 0.0000000000 K 0.5000000000 0.0000000000 0.0000000000 MIf you use this module, please cite the following work:[1] V. Wang, N. Xu, J.-C. Liu, G. Tang, W.-T. Geng, VASPKIT: A User-friendly Interface Facilitating High-throughput Computing and Analysis Using VASP Code, Computer Physics Communications 267, 108033 (2021).[2] V. Wang, Y.-Y. Liang, Y. Kawazeo, W.-T. Geng, High-Throughput Computational Screening of Two-Dimensional Semiconductors, arXiv:1806.04285.

##### initial I/O #####SYSTEM = MoS2ICHARG = 11LWAVE = .TRUE.LCHARG = .TRUE.LVTOT = .FALSE.LVHAR = .FALSE.LELF = .FALSE.LORBIT = 11NEDOS = 1000##### SCF #####ENCUT = 500ISMEAR = 0SIGMA = 0.05EDIFF = 1E-6NELMIN = 5NELM = 300GGA = PELREAL = .FALSE.PREC = Accurate

Post-process Band Structure (pure functional)¶

After calculation, do band structure post-process by VASPKIT option21

Use 211 get the basic band structure. If there are pythonenvironment with matplotlib module. VASPKIT can output a band.pngfigure automatically. By default, the Fermi energy will be shifted tothe 0 eV.

211 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Energy-Levels From EIGENVAL File... -->> (04) Reading Structural Parameters from POSCAR File... -->> (05) Reading K-Paths From KPOINTS File... -->> (06) Written BAND.dat File! -->> (07) Written BAND_REFORMATTED.dat File! -->> (08) Written KLINES.dat File! -->> (09) Written KLABELS File! -->> (10) Written BAND_GAP File! If you want use the default setting, type 0, if modity type 10/public1/home/pg2059/.local/lib/python2.7/site-packages/matplotlib/font_manager.py:1328: UserWarning: findfont: Font family [u'arial'] not found. Falling back to DejaVu Sans (prop.get_family(), self.defaultFamily[fontext])) -->> (11) Graph has been generated!

OutputBAND.dat，BAND_REFORMATTED.dat，KLINES.dat，KLABELS，BAND_GAPfiles，

BAND.dat，BAND_REFORMATTED.dat files save the bandinformation, which can be open by ORIGIN directly.

BAND_REFORMATTED.dat: First column is the length of K-path in unitof Å-1, and following columns are the energies of energy of each bands.

#K-Path Energy-Level 0.000 -14.278 -13.063 -5.798 -2.813 -2.813 -1.962 -1.693 ... 0.060 -14.268 -13.058 -5.787 -2.836 -2.802 -2.055 -1.710 ... 0.120 -14.239 -13.044 -5.752 -2.906 -2.769 -2.226 -1.761 ... 0.180 -14.190 -13.021 -5.695 -3.015 -2.716 -2.416 -1.840 ... 0.240 -14.123 -12.989 -5.619 -3.157 -2.642 -2.614 -1.943 ... 0.300 -14.039 -12.948 -5.528 -3.321 -2.817 -2.551 -2.064 ... 0.360 -13.938 -12.900 -5.428 -3.496 -3.023 -2.443 -2.195 ... 0.420 -13.823 -12.846 -5.329 -3.671 -3.232 -2.333 -2.322 ... 0.480 -13.696 -12.786 -5.244 -3.830 -3.441 -2.471 -2.190 ... ...

KLABELS file saves the positions of high symmetry points on bandstructure figures:

K-Label K-Coordinate in band-structure plotsGAMMA 0.000M 1.140K 1.798GAMMA 3.114* Give the label for each high symmetry point in KPOINTS (KPATH.in) file. Otherwise, they will be identified as 'Undefined or XX' in KLABELS file

The BAND_GAP file save the information of band gap, VBM, CBM, and itslocations at reciprocal lattice,

+-------------------------- Summary ----------------------------+ Band Character: Direct Band Gap (eV): 1.6743 Eigenvalue of VBM (eV): -0.2257 Eigenvalue of CBM (eV): 1.4485 hom*o & LUMO Bands: 9 10 Location of VBM: 0.333333 0.333333 0.000000 Location of CBM: 0.333333 0.333333 0.000000+---------------------------------------------------------------+NOTE: The VBM and CBM are subtracted by the Fermi Energy.

If users have python environment with matplotlib module, VASPKIT canoutput a band.png figure automatically：

212，213，214 get projected band structures：

Ensure the LORBIT = 10 or LORBIT = 11 parameter in INCAR tooutput projection information.

212) Projected Band-Structure for Selected Atoms. select projectedatoms:

Input number of selected atoms：1-4 7 8 24 orelements：C Fe H

For example if you want to draw projected band structure on atom 1 and2, input 1-2：

212 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Structural Parameters from POSCAR File... -->> (04) Reading Energy-Levels From EIGENVAL File... -->> (05) Reading Band-Weights From PROCAR File... -->> (06) Reading K-Paths From KPOINTS File... |---------------------------------------------------------------| Input the element-symbol and/or atom-index to SUM [Total-atom: 3] (Free Format is OK, e.g., C Fe H 1-4 7 8 24) ------------>>1-2 -->> (07) Written SELECTED_ATOM_LIST File! -->> (08) Written PBAND_A1.dat File! -->> (09) Written PBAND_A2.dat File! -->> (10) Written KLINES.dat File! -->> (11) Written KLABELS File!

VASPKIT will output two files: PBAND_A1.dat and PBAND_A2.dat.

Each file include the projection information on select atoms, and thereweight on each angular momentum,s py pz px dxy dyz dz2 dxz x2-y2 tot.The first column is the length of K-path in unit of Å-1. Second columnis the energy of band. Following column are the projection oflm-orbitals on this band. Last column is the total projection ofselected atom on this band.

#K-Path Energy s py pz px dxy dyz dz2 dxz x2-y2 tot#Band-index 1 0.000 -14.278 0.275 0.000 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.281 0.060 -14.268 0.276 0.000 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.281 0.120 -14.239 0.276 0.000 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.281 0.180 -14.190 0.277 0.000 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.282 0.240 -14.123 0.278 0.000 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.284 0.300 -14.039 0.280 0.000 0.005 0.001 0.000 0.000 0.000 0.000 0.000 0.285 0.360 -13.938 0.281 0.000 0.005 0.001 0.000 0.000 0.000 0.000 0.000 0.288 0.420 -13.823 0.283 0.000 0.005 0.001 0.000 0.000 0.000 0.000 0.000 0.290 ...#Band-index 2 3.114 -13.063 0.315 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.316 3.045 -13.056 0.315 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.316 ...#Band-index 3 0.000 -5.798 0.018 0.000 0.125 0.000 0.000 0.000 0.000 0.000 0.000 0.143 0.060 -5.787 0.018 0.000 0.125 0.000 0.000 0.000 0.000 0.000 0.000 0.143 ...

213) Projected Band-Structure for Each Element，projection on eachelements:

213 -->> (1) Reading Input Parameters From INCAR File... -->> (2) Reading Fermi-Level From DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (3) Reading Structural Parameters from POSCAR File... -->> (4) Reading Energy-Levels From EIGENVAL File... -->> (5) Reading Band-Weights From PROCAR File... -->> (6) Reading K-Paths From KPOINTS File... -->> (7) Written PBAND_S.dat File! -->> (8) Written PBAND_Mo.dat File! -->> (9) Written KLINES.dat File! -->> (*) Written KLABELS File!

Format of PBAND_S.dat and PBAND_Mo.dat is same asPBAND_A1.dat. These files can be open by origin:

Select the first and second line, draw the original band structurefirst:

Then find the position of high symmetry points from KLABELS:

GAMMA 0.000M 1.139K 1.797GAMMA 3.113

Label those points:

Then, adjust the energy range. The original band structure withoutprojection information shows as follow.

From ORIGIN plot setup tag, add the atoms, elements, or orbitalsprojections on band.

Select Bubble tag, and add the projection. click OK.

214) The Sum of Projected Band-Structure for Selected Atoms：

214 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Structural Parameters from POSCAR File... -->> (04) Reading Energy-Levels From EIGENVAL File... -->> (05) Reading Band-Weights From PROCAR File... -->> (06) Reading K-Paths From KPOINTS File... |---------------------------------------------------------------| Input the element-symbol and/or atom-index to SUM [Total-atom: 3] (Free Format is OK, e.g., C Fe H 1-4 7 8 24) ------------>>1-2 -->> (07) Written SELECTED_ATOM_LIST File! -->> (08) Written PBAND_SUM.dat File! -->> (09) Written KLINES.dat File! -->> (10) Written KLABELS File!

The output PBAND_SUM.dat file includes the sum projection ofselected atoms or elements. This is useful to investigate the layeredband and compare surface band and inner band.

Single-shot Band Structue¶

The VASP.6.3.0 or latter version can calculate band structure through a single-shot.Note the PROCAR_OPT file does not output unless we compiled VASP after adding -DVASP_HDF5 to the CPP_OPTIONS in the makefile.include.Alternatively, the PROCAR_OPT and EIGENVAL_OPT (new added) can automatically be generated after we modified the vasp.6.3.0/src/linear_response.F and recomplied VASP again.Specifically, first delete the lines from 1804-1827 in vasp.6.3.0/src/linear_response.F file,

#ifdef VASP_HDF5 CALL VH5_WRITE_DOS(IH5OUTFILEID, WDES_INTER, KPOINTS_INTER, DOS, DOSI, DOSPAR, EFERMI, & T_INFO%NIONP, -1,SUBGROUP="electron_dos_kpoints_opt") CALL VH5_WRITE_EIGENVAL(IH5OUTFILEID, WDES_INTER, W_INTER, KPOINTS_INTER, & SUBGROUP="electron_eigenvalues_kpoints_opt") ! Write projections IF (JOBPAR/=0 .OR. IO%LORBIT>=10) THEN IF (IO%LORBIT==11.OR.IO%LORBIT==12) THEN CALL VH5_WRITE_PROJECTORS(IH5OUTFILEID, T_INFO, WDES_INTER, W_INTER, IO%LORBIT, LPAR, PAR_INTER, LMCHAR, PHAS_INTER, & SUBGROUP="projectors_kpoints_opt") ELSE CALL VH5_WRITE_PROJECTORS(IH5OUTFILEID, T_INFO, WDES_INTER, W_INTER, IO%LORBIT, LPAR, PAR_INTER, LCHAR, PHAS_INTER, & SUBGROUP="projectors_kpoints_opt") ENDIF IF (IO%IU6>=0) THEN CALL XML_PROCAR(PAR, W_INTER%CELTOT, W_INTER%FERTOT, WDES_INTER%NB_TOT, WDES_INTER%NKPTS, LPAR,& T_INFO%NIONP, WDES_INTER%NCDIJ, TAG="projected_kpoints_opt") OPEN(UNIT=99,FILE='PROCAR_OPT',STATUS='UNKNOWN') CALL WRITE_PROCAR(99,W_INTER,WDES_INTER,IO%LORBIT,LPAR,LMDIMP,PAR_INTER,PHAS_INTER,T_INFO,KPOINTS%EMIN,KPOINTS%EMAX) CLOSE(99) ENDIF ENDIF#endif

and then add the following code

! Added support to VASPKIT OPEN(UNIT=99,FILE='EIGENVAL_OPT',STATUS='UNKNOWN') WRITE(99,'(4I5)') T_INFO%NIONS,T_INFO%NIONS,0.0,WDES_INTER%ISPIN WRITE(99,'(5E15.7)') 0.d0,0.d0,0.d0,0.d0,0.d0 WRITE(99,'(1E15.7)') 0.d0 WRITE(99,*) ' CAR ' WRITE(99,*) INFO%SZNAM1 WRITE(99,'(3I7)') NINT(INFO%NELECT),WDES_INTER%NKPTS,WDES_INTER%NB_TOT DO IK=1,WDES_INTER%NKPTS WRITE(99,*) WRITE(99,'((4E15.7))') WDES_INTER%VKPT(1,IK),WDES_INTER%VKPT(2,IK),WDES_INTER%VKPT(3,IK),WDES_INTER%WTKPT(IK) DO I=1,WDES_INTER%NB_TOT IF (WDES_INTER%ISPIN==1) WRITE(99,852) I,REAL( W_INTER%CELTOT(I,IK,1) ,KIND=q),W%FERTOT(I,IK,1) IF (WDES_INTER%ISPIN==2) WRITE(99,853) I,(REAL( W_INTER%CELTOT(I,IK,ISP) ,KIND=q),ISP=1,WDES_INTER%ISPIN),(W%FERTOT(I,IK,ISP),ISP=1,WDES_INTER%ISPIN) ENDDO ENDDO CLOSE(99) 852 FORMAT(1X,I6,4X,F14.6,2X,F9.6) 853 FORMAT(1X,I6,4X,F14.6,2X,F14.6,2X,F9.6,2X,F9.6) IF (JOBPAR/=0 .OR. IO%LORBIT>=10) THEN IF (IO%IU6>=0) THEN CALL XML_PROCAR(PAR, W_INTER%CELTOT, W_INTER%FERTOT, WDES_INTER%NB_TOT, WDES_INTER%NKPTS, LPAR,& T_INFO%NIONP, WDES_INTER%NCDIJ, TAG="projected_kpoints_opt") OPEN(UNIT=99,FILE='PROCAR_OPT',STATUS='UNKNOWN') CALL WRITE_PROCAR(99,W_INTER,WDES_INTER,IO%LORBIT,LPAR,LMDIMP,PAR_INTER,PHAS_INTER,T_INFO,KPOINTS%EMIN,KPOINTS%EMAX) CLOSE(99) ENDIF ENDIF#ifdef VASP_HDF5 CALL VH5_WRITE_DOS(IH5OUTFILEID, WDES_INTER, KPOINTS_INTER, DOS, DOSI, DOSPAR, EFERMI, & T_INFO%NIONP, -1,SUBGROUP="electron_dos_kpoints_opt") CALL VH5_WRITE_EIGENVAL(IH5OUTFILEID, WDES_INTER, W_INTER, KPOINTS_INTER, & SUBGROUP="electron_eigenvalues_kpoints_opt") ! Write projections IF (JOBPAR/=0 .OR. IO%LORBIT>=10) THEN IF (IO%LORBIT==11.OR.IO%LORBIT==12) THEN CALL VH5_WRITE_PROJECTORS(IH5OUTFILEID, T_INFO, WDES_INTER, W_INTER, IO%LORBIT, LPAR, PAR_INTER, LMCHAR, PHAS_INTER, & SUBGROUP="projectors_kpoints_opt") ELSE CALL VH5_WRITE_PROJECTORS(IH5OUTFILEID, T_INFO, WDES_INTER, W_INTER, IO%LORBIT, LPAR, PAR_INTER, LCHAR, PHAS_INTER, & SUBGROUP="projectors_kpoints_opt") ENDIF! IF (IO%IU6>=0) THEN! CALL XML_PROCAR(PAR, W_INTER%CELTOT, W_INTER%FERTOT, WDES_INTER%NB_TOT, WDES_INTER%NKPTS, LPAR,&! T_INFO%NIONP, WDES_INTER%NCDIJ, TAG="projected_kpoints_opt")! OPEN(UNIT=99,FILE='PROCAR_OPT',STATUS='UNKNOWN')! CALL WRITE_PROCAR(99,W_INTER,WDES_INTER,IO%LORBIT,LPAR,LMDIMP,PAR_INTER,PHAS_INTER,T_INFO,KPOINTS%EMIN,KPOINTS%EMAX)! CLOSE(99)! ENDIF ENDIF#endif

Finally recomplie VASP again. The VASPKIT.1.3.2 or latter version can read the raw data of plain and projected band structures from the EIGENVAL_OPT and PROCAR_OPT files, respectively.

Pre-process Band Structure (hybrid functional)¶

Besides the line-mode K-path, VASPKIT can also generate K-path withuniformed spacing between K points in units of 2\(\pi\)*Å-1,which can be used for hybrid functional band structure calculations.Such KPOINTS file contains two parts. First part is same as theself-consistent calculation with symmetry weighted K-points inIrreducilbe Brillouin Zone. And Second part is the 0-weighted K-pointsalone the k-path. To generate this KPOINTS file:

1. Same as pure functional calculation. Do geometry optimization. Run302 (for 2D materials) or 303 (for 3D materials) to get:

1. standardized primitive cell (PRIMCELL.vasp)

2. line-mode k-path for pure functional band structure calculation(KPATH.in)

3. High-symmetry points in fractional coordinates.(HIGH_SYMMETRY_POINTS) You can check them with the seekpath database[https://www.materialscloud.org/work/tools/seekpath].

2. Run 251 to generate KPOINTS file for hybrid functionalband-structure calculations. Input the KPT resolution values todetermine density of k-mesh for SCF calculation and k-path for bandstructure calculation. Then VASPKIT will read KPATH.in file andgenerate the KPOINTS file for hybrid functional band-structurecalculation.

3. Optional. Do a PBE SCF calculations based on the new generatedKPOINTS file and get the wavefunctions, which can be read for nextstep hybrid functional calculation. Sometimes, this step reduce theSCF time of the next step hybrid functional calculation.

4. Do hybrid functional band structure calculation.

302 +-------------------------- Warm Tips --------------------------+ See An Example in vaspkit/examples/seek_kpath/graphene_2D. This feature is still experimental & check the PRIMCELL.vasp file. +---------------------------------------------------------------+ -->> (1) Reading Structural Parameters from POSCAR File... +-------------------------- Summary ----------------------------+ The vacuum slab is supposed to be along c axis Prototype: AB2 Total Atoms in Input Cell: 3 Lattice Constants in Input Cell: 3.184 3.184 18.410 Lattice Angles in Input Cell: 90.000 90.000 120.000 Total Atoms in Primitive Cell: 3 Lattice Constants in Primitive Cell: 3.184 3.184 18.410 Lattice Angles in Primitive Cell: 90.000 90.000 120.000 2D Bravais Lattice: Hexagonal Space Group: 187 Point Group: 26 [ D3h ] International: P-6m2 Symmetry Operations: 12 Suggested K-Path: (shown in the next line) [ GAMMA-M-K-GAMMA ] +---------------------------------------------------------------+ -->> (2) Written PRIMCELL.vasp file. -->> (3) Written KPATH.in File for Band-Structure Calculation. -->> (4) Written HIGH_SYMMETRY_POINTS File for Reference.

Primitive Cell 1.000000 3.18401832481292 0.00000000000000 0.00000000000000 -1.59200916240646 2.75744075540316 0.00000000000000 0.00000000000000 0.00000000000000 18.40999984740000 S Mo 2 1DIRECT 0.0000000000000000 0.0000000000000000 0.4151505217091287 S1 0.0000000000000000 0.0000000000000000 0.5848494782908713 S2 0.6666666666666666 0.3333333333333333 0.5000000000000000 Mo1

K-Path Generated by VASPKIT. 20Line-ModeReciprocal 0.0000000000 0.0000000000 0.0000000000 GAMMA 0.5000000000 0.0000000000 0.0000000000 M 0.5000000000 0.0000000000 0.0000000000 M 0.3333333333 0.3333333333 0.0000000000 K 0.3333333333 0.3333333333 0.0000000000 K 0.0000000000 0.0000000000 0.0000000000 GAMMA

251 -->> (1) Reading Structural Parameters from POSCAR File... ======================= K-Mesh Scheme ========================== 1) Monkhorst-Pack Scheme 2) Gamma Scheme 0) Quit 9) Back ------------->>2 +-------------------------- Warm Tips --------------------------+ Input Resolution Value to Determine K-Mesh for SCF Calculation: (Typical Value: 0.03-0.04 is Generally Precise Enough) ------------>>0.05 Input Resolution Value along K-Path for Band Calculation: (Typical Value: 0.02-0.04 for DFT and 0.04-0.06 for hybrid DFT) ------------>>0.05 +---------------------------------------------------------------+ -->> (2) Reading K-Paths From KPATH.in File... +-------------------------- Summary ----------------------------+ K-Mesh for SCF Calculation: 7 7 1 The Number of K-Point along K-Path 1: 22 The Number of K-Point along K-Path 2: 13 The Number of K-Point along K-Path 3: 26 +---------------------------------------------------------------+ -->> (3) Written KPOINTS File!

Parameters to Generate KPOINTS (Don't Edit This Line): 0.050 0.050 8 61 3 22 13 26 69Reciprocal lattice 0.00000000000000 0.00000000000000 0.00000000000000 1 0.14285714285714 0.00000000000000 0.00000000000000 6 0.28571428571429 0.00000000000000 0.00000000000000 6 0.42857142857143 0.00000000000000 0.00000000000000 6 0.14285714285714 0.14285714285714 0.00000000000000 6 0.28571428571429 0.14285714285714 0.00000000000000 12 0.42857142857143 0.14285714285714 0.00000000000000 6 0.28571428571429 0.28571428571429 0.00000000000000 6 0.00000000000000 0.00000000000000 0.00000000000000 0 0.02380952380952 0.00000000000000 0.00000000000000 0 ... ... ... 0.02666666666400 0.02666666666400 0.00000000000000 0 0.01333333333200 0.01333333333200 0.00000000000000 0 0.00000000000000 0.00000000000000 0.00000000000000 0

##### initial I/O #####SYSTEM = Hybird FunctionalISTART = 1ICHARG = 1LWAVE = .TRUE.LCHARG = .TRUE.LVTOT = .FALSE.LVHAR = .FALSE.LELF = .FALSE.LORBIT = 11NEDOS = 1000##### SCF #####ENCUT = 500ISMEAR = 0SIGMA = 0.05EDIFF = 1E-6NELMIN = 5NELM = 300GGA = PELREAL = .FALSE.# PREC = Accurate# ISYM = 3##### Geo opt #####EDIFFG = -0.01IBRION = 2POTIM = 0.2NSW = 0 # no Geo_optISIF = 3 # 2 not optimize cell#### HSE ####LHFCALC = .TRUE.AEXX = 0.25HFSCREEN = 0.2ALGO = ALL # or DampedTIME = 0.4

Post-process Band Structure (hybrid functional)¶

5. Extract band structure information by 252. band data is saved inBAND.dat and BAND-REFORMATTED.dat files. High-symmetrypoints positions on band structure figures are saved in KLABELS.

(Note: bug in 0.73 version, please use new version)

252 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Energy-Levels From EIGENVAL File... -->> (04) Reading Structural Parameters from POSCAR File... -->> (05) Reading K-Paths From KPATH.in File... -->> (06) Written BAND.dat File! -->> (07) Written BAND_REFORMATTED.dat File! -->> (08) Written KLINES.dat File! -->> (09) Written KLABELS File! -->> (10) Written BAND_GAP File! If you want use the default setting, type 0, if modity type 10 -->> (11) Graph has been generated!

if the python and matplotlib enviornment is set correctly. VASPKIT willdraw a band figure band.png automatically. Following parametersshould be set in ~/.vaspkit file.

PYTHON_BIN ~/anaconda3/bin/python3PLOT_MATPLOTLIB .TURE.

You can also draw the figure from BAND.dat orBAND_REFORMATTED.dat by ORIGIN, which is described in Section 3.2.

By comparing with the line-mode of KPOINTS (option 302 and 303).The biggest advantages of 25 is that the k spacing along K-path isaveraged, so that saves computational cost.

Versions prior to 0.72 took the same number of K-points for energy line,resulting in uneven K-distribution on different paths, as shown in theleft figure below. The latest version of VASPKIT supports automaticdetermination of K-points number on different energy band paths based ona given k-point interval, thus ensuring uniform spattering throughoutthe band calculation, as shown in right Figure below.

KPOINTS generated from vaspkit 251 keeps the spacing between every kpoints. So, one can use less 0-weighted k points to get a similarqualified band structure, and thus, take less time when doing bandstructure calculations.

To get projected hybrid functional band structure, please use253，254，255：

Ensure the LORBIT = 10 or LORBIT = 11 parameter in INCAR tooutput projection information.

253) Get Projected Band-Structure for Selected Atoms254) Get Projected Band-Structure for Each ELement255) Get the Sum of Projected Band-Structure for Selected Atoms

For example: by using 253, a free format input is available. You caninput any atoms by its index and any elements symbol. Projected bandstructure will saved for each elements one by one.

 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Structural Parameters from POSCAR File... -->> (04) Reading Energy-Levels From EIGENVAL File... -->> (05) Reading Band-Weights From PROCAR File... -->> (06) Reading K-Paths From KPATH.in File... |---------------------------------------------------------------| Input the element-symbol and/or atom-index to SUM [Total-atom: 3] (Free Format is OK, e.g., C Fe H 1-4 7 8 24) ------------>>1-3 -->> (07) Written SELECTED_ATOM_LIST File! -->> (08) Written PBAND_A1.dat File! -->> (09) Written PBAND_A2.dat File! -->> (10) Written PBAND_A3.dat File! -->> (11) Written KLINES.dat File! -->> (12) Written KLABELS File!

Elements projected band structure obtained by using 254,

254 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Structural Parameters from POSCAR File... -->> (04) Reading Energy-Levels From EIGENVAL File... -->> (05) Reading Band-Weights From PROCAR File... -->> (06) Reading K-Paths From KPATH.in File... -->> (07) Written PBAND_S.dat File! -->> (08) Written PBAND_Mo.dat File! -->> (09) Written KLINES.dat File! -->> (10) Written KLABELS File!

All your input atoms will be sumed up and are projected to the bandstructure by using 255:

255 -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading Structural Parameters from POSCAR File... -->> (04) Reading Energy-Levels From EIGENVAL File... -->> (05) Reading Band-Weights From PROCAR File... -->> (06) Reading K-Paths From KPATH.in File... |---------------------------------------------------------------| Input the element-symbol and/or atom-index to SUM [Total-atom: 3] (Free Format is OK, e.g., C Fe H 1-4 7 8 24) ------------>>1-2 -->> (07) Written SELECTED_ATOM_LIST File! -->> (08) Written PBAND_SUM.dat File! -->> (09) Written KLINES.dat File! -->> (10) Written KLABELS File!

supplementary: band structure can also be down by pymatgen package:

Effective Mass¶

For electrons or holes in a solid, the effective mass (m*) is aquantity that is used to simplify band structures by modeling thebehavior of a free particle with that mass. At the highest energies ofthe valence band in semiconductors (Ge, Si, GaAs, …), and the lowestenergies of the conduction band in semiconductors (GaAs, …), the bandstructure E(k) can be locally approximated as:

\[E(\mathbf{k})=E_{0}+\frac{\hbar^{2} \mathbf{k}^{2}}{2 m^{*}}\]

where E(k) is the energy of an electron at wavevector k inthat band, E0 is a constant giving the edge of energy of that band.So m* can be calculated by following equation:

\[\frac{1}{m^{*}}=\frac{1}{\hbar^{2}} \frac{\partial^{2} E}{\partial k_{i} \partial k_{j}}\]

Usually, we are interested in the m* from one high symmetry pointto another high symmetry point. So the direction of \(k_{i}\) and\(k_{j}\) are same, such as K -> M, K -> \(\Gamma\). To get thesecond-order partial derivative, KPOINTS file should include K pointsvery near to that high symmetry point. VASPKIT can generate KPOINTS filefor VASP calculation and get the effective mass automatically. (Note:Now, VASPKIT only support effective mass calculation for Non-Charged &Non-Magnetic Semiconductor!)

1 # "1" for pre-process (generate KPOINTS), "2" for post-process(calculate m*)6 # number of points for quadratic function fitting.0.015 # k-cutoff, unit Å-1.2 # number of tasks for effective mass calculation0.333333 0.3333333 0.000 0.000 0.000 0.000 K->Γ # Specified two K-points and direction0.333333 0.3333333 0.000 0.500 0.000 0.000 K->M # Specified two K-points and direction

 +-------------------------- Warm Tips --------------------------+ Test ONLY for Non-Charged & Non-Magnetic Semiconductor! +---------------------------------------------------------------+ -->> (01) Reading VPKIT.in File... -->> (02) Reading Structural Parameters from POSCAR File... ======================= K-Mesh Scheme ========================== 1) Monkhorst-Pack Scheme 2) Gamma Scheme 0) Quit 9) Back ------------->>2 +-------------------------- Warm Tips --------------------------+- Accuracy Levels: Gamma-Only: 0; Low: 0.06~0.04; Medium: 0.04~0.03; Fine: 0.02-0.01.- 0.03-0.04 is Generally Precise Enough! +---------------------------------------------------------------+ Input KPT-Resolved Value (e.g., 0.04, in unit of 2*PI/Angstrom): ------------>> 0.04 1 -->> (03) Written KPOINTS File! -->> (04) Written INCAR file! -->> (05) Written POTCAR File with the Standard Potential!

Parameters to Generate KPOINTS (Don't Edit This Line): 0.0400 12 2 6 24Reciprocal lattice 0.00000000000000 0.00000000000000 0.00000000000000 1 0.11111111111111 0.00000000000000 0.00000000000000 6 0.22222222222222 0.00000000000000 0.00000000000000 6 0.33333333333333 0.00000000000000 0.00000000000000 6 0.44444444444444 0.00000000000000 0.00000000000000 6 0.11111111111111 0.11111111111111 0.00000000000000 6 0.22222222222222 0.11111111111111 0.00000000000000 12 0.33333333333333 0.11111111111111 0.00000000000000 12 0.44444444444444 0.11111111111111 0.00000000000000 6 0.22222222222222 0.22222222222222 0.00000000000000 6 0.33333333333333 0.22222222222222 0.00000000000000 12 0.33333333333333 0.33333333333333 0.00000000000000 2 0.33333300000000 0.33333330000000 0.00000000000000 0 # K->Γ 0.32795808307727 0.33136593967371 0.00000000000000 0 0.32258316615454 0.32939857934742 0.00000000000000 0 0.31720824923180 0.32743121902113 0.00000000000000 0 0.31183333230907 0.32546385869484 0.00000000000000 0 0.30645841538634 0.32349649836855 0.00000000000000 0 0.33333300000000 0.33333330000000 0.00000000000000 0 # K->M 0.33673240318272 0.32574567174541 0.00000000000000 0 0.34013180636544 0.31815804349082 0.00000000000000 0 0.34353120954815 0.31057041523622 0.00000000000000 0 0.34693061273087 0.30298278698163 0.00000000000000 0 0.35033001591359 0.29539515872704 0.00000000000000 0

913 +-------------------------- Warm Tips --------------------------+ Test ONLY for Non-Charged & Non-Magnetic Semiconductor! +---------------------------------------------------------------+ -->> (01) Reading VPKIT.in File... -->> (02) Reading Structural Parameters from POSCAR File... -->> (03) Reading Input Parameters From INCAR File... -->> (04) Reading Energy-Levels From EIGENVAL File... +-------------------------- Summary ----------------------------+ Effective-Mass is obtained by fitting a third order polynomial, which yields masses are less sensitive to the employed k-points. Band Index: LUMO = 10 hom*o = 9 Effective-Mass (in m0) Electron (Prec.) Hole (Prec.) K-Path Index 1: K->Γ 0.480 (0.2E-07) -0.577 (0.2E-07) K-Path Index 2: K->M 0.527 (0.8E-07) -0.705 (0.2E-06) +---------------------------------------------------------------+

Band-unfolding¶

The electronic structures of real materials are perturbed by variousstructural defects, impurities, fl?uctuations of the chemical compositionin compound alloys and so on. In DFT calcualtions, these defects andincommensurate structures are usually investigated by using supercell(SC) models. However, interpretation of \(E(\vec{k})\) versus \(\vec{k}\)band structures is most effective within the primitive cell (PC).Popescu and Zunger [see PRL 104, 236403 (2010) and PRB 85, 085201(2012)] proposed the effective band structures (EBS) method which canunfold the band structures of supercells into the Brillouin zone (BZ) ofthe correspoind primitive cell. Such band unfolding techniques greatlysimplify the analysis of the results and enabling direct comparisonswith experimental measurements, for example, angle-resolvedphotoemission spectroscopy (ARPES), often represented along thehigh-symmetry directions of the primtive cell BZ.

Next we will build a rectangle-like MoS\(_2\) supercell withsingle sulfur vacancy and perfom band unfolding calcualtion using VASPtogether with VASPKIT codes.

First step: Prepare the following files;

1. POSCAR for monolayer MoS\(_2\)

   primitive cell

MoS2-H.POSCAR 1.00000000000000000 3.1904063769892548 0.0000000000000000 0.0000000000000000 -1.5952031884946274 2.7629729709066906 0.0000000000000000 0.0000000000000000 0.0000000000000000 17.6585355706663840 Mo S 1 2Direct 0.3333333429999996 0.6666666870000029 0.5000000000000000 0.6666666870000029 0.3333333429999996 0.4117002781430051 0.6666666870000029 0.3333333429999996 0.5882996918569924

2. KPATH.in file for MoS\(_2\)

   primitive cell (not for supercell), you can genenate it by runvaspkit with task 302 for 2D or 303 for bulk, or edit it manually. Wechoose the K-Path of M-K-\(\Gamma\) for monolayer MoS2.

KPATH for MoS220Line modeRec 0.00000000000 0.50000000000 0.00000000000 # M 0.33333333333 0.33333333333 0.00000000000 # K 0.33333333333 0.33333333333 0.00000000000 # K 0.00000000000 0.00000000000 0.00000000000 # GAMMA

3. INCAR file for static calculation，set LWAVE =.TRUE.. To aviod theerror ’ERROR! Computed NPLANE= **** != input NPLANE= ****’, pleaseset PREC=N or a medium/lower value of energy cuttoff.

SYSTEM = MoS2ISTART = 1LREAL = FPREC = NLWAVE = .TRUE.LCHARG = FISMEAR = 0SIGMA = 0.05NELM = 60NELMIN = 6EDIFF = 1E-08

4. TRANSMAT.in file (optional). The content of this file is

Read transformation matrix from the TRANSMAT.in file if it exists. 2 0 0 # must be three integers 0 3 0 # must be three integers 0 0 3 # must be three integers

where the first line is the comment line, and the integers from secondline to four lines are the elements of 3$:raw-latex:times`$3 thetransmission matrix :math:`M, the relationship of lattice vectorsbetween supercell \(E(\vec{A})\) satisfies \(E(\vec{a})\)

\[\begin{split}\left(\begin{array}{l}{\vec{A}_{1}} \\ {\vec{A}_{2}} \\ {\vec{A}_{3}}\end{array}\right)=M \cdot\left(\begin{array}{l}{\vec{a}_{1}} \\ {\vec{a}_{2}} \\ {\vec{a}_{3}}\end{array}\right)=\left(\begin{array}{lll}{m_{11}} & {m_{12}} & {m_{13}} \\ {m_{21}} & {m_{22}} & {m_{23}} \\ {m_{31}} & {m_{32}} & {m_{33}}\end{array}\right) \cdot\left(\begin{array}{l}{\vec{a}_{1}} \\ {\vec{a}_{2}} \\ {\vec{a}_{3}}\end{array}\right)\end{split}\]

If the off-diagonal elements in \(M\) are zero, the lattice vectorsof supercell will be parallel to those of primtivel cell.

Second step: Run vapskit with task 400 to generate SUPERCELL.vasp andTRANSMAT (if TRANSMAT.in not exist). We modify the SUPERCELL.vasp fileby removing a sulfur atom to creat signle S vacancy(V:math:`_text{S}`) point defect;

The primitive cell of MoS\(_2\) monolayer belongs to hexagonalcrystal system. We can build a rectangle-supercell by using vaspkit withtask 400. If the TRANSMAT.in file does not exist, vaspkit will read the9 elements of \(M\) by dialogue menu (see the following commands).Otherwise, the elements of \(M\) will be read from the TRANSMAT.infile.

vaspkit -task 400 +---------------------------------------------------------------+ | VASPKIT Version: 1.00.RC (17 Jun. 2019) | | A Pre- and Post-Processing Program for VASP Code | | Running VASPKIT Under Command-Line Mode | +---------------------------------------------------------------+ -->> (01) Reading Structural Parameters from POSCAR File... Enter the new lattice verctor a in terms of old: (MUST be three integers, e.g., 1 2 3) 4 0 0 Enter the new lattice verctor b in terms of old: 2 4 0 Enter the new lattice verctor c in terms of old: 0 0 1 +-------------------------- Summary ----------------------------+ The Transformation Matrix P is: 4 0 0 2 4 0 0 0 1 Lattice Constants in New Cell: 12.762 11.052 17.659 Lattice Angles in New Cell: 90.00 90.00 90.00 Total Atoms in New Cell: 48 Volume of New Cell is 16 times of the Old Cell +---------------------------------------------------------------+ -->> (02) Written SUPERCELL.vasp File

Third step: cp SUPERCELL.vasp POSCAR and cp TRANSMAT TRANSMAT.in (ifTRANSMAT.in not exist) and run vaspkit with task 281 to generate KPOINTSfile;

 ------------>>28 +-------------------------- Warm Tips --------------------------+ See some examples in vaspkit/examples/band_unfolding. Only Support KPOINTS file Generated by VASPKIT. Please Set LWAVE= .TRUE. in the INCAR file. +---------------------------------------------------------------+ ================== Band-Unfolding ============================== 281) Generate KPOINTS File for Band-Unfolding Calculation 282) Calculate Effective Band Structure 0) Quit 9) Back ------------>>281 -->> (01) Reading Structural Parameters from POSCAR File... +---------------------------------------------------------------+ | Selective Dynamics is Activated! | +---------------------------------------------------------------+ ======================= K-Mesh Scheme ========================== 1) Monkhorst-Pack Scheme 2) Gamma Scheme 0) Quit 9) Back ------------->>2 +-------------------------- Warm Tips --------------------------+ Input Resolution Value to Determine K-Mesh for SCF Calculation: (Typical Value: 0.02-0.03 is Generally Precise Enough) ------------>>0.03 Input Resolution Value along K-Path for Band Calculation: (Typical Value: 0.02-0.04 for DFT and 0.04-0.06 for hybrid DFT) ------------>>0.03 +---------------------------------------------------------------+ -->> (02) Readin Transformation Matrix from TRANSMAT.in file... -->> (03) Reading K-Paths From KPATH.in File... +-------------------------- Summary ----------------------------+ K-Mesh for SCF Calculation: 3 3 1 # Based on SC reciprocal lattice length The Number of K-Point along K-Path 1: 21 # Based on PC reciprocal lattice length The Number of K-Point along K-Path 2: 43 # Based on PC reciprocal lattice length +---------------------------------------------------------------+

Forth step: Submit vasp job；

Fifth step: Run vaspkit with task 282 to get EBS；

 ------------>>282 +-------------------------- Warm Tips --------------------------+ Current Version Only Support the Stardard Version of VASP code. +---------------------------------------------------------------+ -->> (01) Reading the header information in WAVECAR file... +--------------------- WAVECAR Header --------------------------+ SPIN = 1 NKPTS = 68 NBANDS = 252 ENCUT = 280.00 Coefficients Precision: Complex*8 Maximum number of G values: GX = 18, GY = 16, GZ = 25 Estimated maximum number of plane-waves: 30159 +---------------------------------------------------------------+ -->> (02) Start to read WAVECAR file, take your time ^.^ Percentage complete: 25.0% Percentage complete: 50.0% Percentage complete: 75.0% Percentage complete: 100.0% -->> (03) Readin Transformation Matrix from TRANSMAT.in file... -->> (04) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (05) Reading KPOINTS_MAPPING_TABLE.in file... -->> (06) Reading K-Paths From KPATH.in File... -->> (07) Start to Calculate Effective Band Structure... Percentage complete: 25.0% Percentage complete: 50.0% Percentage complete: 75.0% Percentage complete: 100.0% -->> (08) Written EBS.dat File! -->> (09) Written KLABELS File!

Last step: You can use vaspkit/examples/band_unfolding/eps_plot.py orother scripts to plot.

3D Band Structure¶

In solid-state physics, the most usual electronic band structure (orsimply band structure) is 2D band structure, which shows the energychange along the K-path line.

VASPKIT can also do 3D band structure, which select K-path on a surfaceof Irreducilbe Brillouin zone. And calculate the K-dependent bandenergies on those K-points. This method is often used to 2D materials,such as graphene.

231 +-------------------------- Warm Tips --------------------------+ * Accuracy Levels: (1) Low: 0.04~0.03; (2) Medium: 0.03~0.02; (3) Fine: 0.02~0.01. * 0.015 is Generally Precise Enough! +---------------------------------------------------------------+ Input KP-Resolved Value (unit: 2*PI/Ang): ------------>>0.015 -->> (01) Reading Structural Parameters from POSCAR File... -->> (02) Written KPOINTS File!

------------>>233 +-------------------------- Warm Tips --------------------------+ ONLY Reliable for Band-Structure Calculations! +---------------------------------------------------------------+ -->> (1) Reading Input Parameters From INCAR File... -->> (2) Reading Structural Parameters from POSCAR File... -->> (3) Reading Fermi-Level From DOSCAR File...ooooooooo The Fermi Energy will be set to zero eV oooooooooooooooo -->> (4) Reading Energy-Levels From EIGENVAL File... -->> (5) Reading Kmesh From KPOINTS File... -->> (6) Written BAND.hom*o.grd File. -->> (7) Written BAND.LUMO.grd File. -->> (8) Written KX.grd and KY.grd Files.

Extract and output DOS and PDOS¶

Option 11 is used to DOS post-process.

 Input the Symbol and/or Number of Atoms to Sum [Total-atom: 7] (Free Format is OK, e.g., C Fe H 1-4 7 8 24),Press "Enter" if you want to end e ntry! ------------>>O Input the Orbitals to Sum Which orbital? s py pz px dxy dyz dz2 dxz dx2 f-3 ~ f3, "all" for summing ALL.s......

 #Energy O_s O_p C_s C_p Ni_d-27.10266 0.00000 0.00000 0.00000 0.00000 0.00000-26.92966 0.00000 0.00000 0.00000 0.00000 0.00000-26.75566 0.00000 0.00000 0.00000 0.00000 0.00000-26.58266 0.00000 0.00000 0.00000 0.00000 0.00000......

d-band Center¶

The DOS projected onto the d-states that interact with the adsorbatestate can be characterized by the center of d-projected DOS.

\[\varepsilon_{\mathrm{d}}=\frac{\int_{-\infty}^{\infty} n_{\mathrm{d}}(\varepsilon) \varepsilon d \varepsilon}{\int_{-\infty}^{\infty} n_{\mathrm{d}}(\varepsilon) d \varepsilon}\]

VASPKIT option 503 can calculate d-band center for every atomsdirectly from VASP DOS outputs. And users can set the the energy rangefor calculation.

Example: Pt metal cell, the energy window is set up to the Fermi energy.

503 +-------------------------- Warm Tips --------------------------+ See An Example in vaspkit/examples/d_band_center. d-Band Center is Sensitive to the Number of Unoccupied Band. Anyway, the Trends are More Important than the Absolute Energies. +---------------------------------------------------------------+ -->> (01) Reading Input Parameters From INCAR File... -->> (02) Reading Fermi-Energy from DOSCAR File... ooooooooo The Fermi Energy will be set to zero eV ooooooooooooooo -->> (03) Reading DOS Data From DOSCAR File... -->> (04) Reading Structural Parameters from POSCAR File... +------------------------ Your Option --------------------------+ The default energy window of integration is [-11.82 24.30] Do you want to change the change (Y/N)?y Please enter the new energy window with blank (e.g., -10 10)-11.8 0 +---------------------------------------------------------------+ -->> (05) Written D_BAND_CENTER File!

the d-band center for every atoms and the total d-band center is writenin D_BAND_CENTER file as following:

# Atom ID d-Band-Center (eV) Pt1: -3.1184 Total: -3.1184Notice:The energy window of integration is [-11.80 0.00].The calculated d-band-center presented here are referenced to Fermi level, i.e., E_F=0 eV.

Gas Molecule Free Energy Correction¶

The equations used for computing thermochemical data for Gas in VASPKITis equivalent to those in Gaussian.(https://gaussian.com/thermo/) The partition function from any componentcan be used to determine the entropy contribution S from thatcomponent:

\[S=N k_{B}+N k_{B} \ln \left(\frac{q(V, T)}{N}\right)+N k_{B} T\left(\frac{\partial \ln q}{\partial T}\right)_{V}\]

When molar values are given \(n=N / N_{A}\), and based on Ideal gasapproximation, \(N_{A} k_{B}=R\), We write the S as:

\[\begin{split}\begin{aligned} S &=R+R \ln (q(V, T))+R T\left(\frac{\partial \ln q}{\partial T}\right)_{V} \\ &=R \ln (q(V, T) e)+R T\left(\frac{\partial \ln q}{\partial T}\right)_{V} \\ &=R\left(\ln \left(q_{\mathrm{t}} q_{\mathrm{e}} q_{\mathrm{r}} q_{\mathrm{v}} e\right)+T\left(\frac{\partial \ln q}{\partial T}\right)_{V}\right) \end{aligned}\end{split}\]

The internal thermal energy U can also be obtained from the partitionfunction:

\[U=N k_{B} T^{2}\left(\frac{\partial \ln q}{\partial T}\right)_{V}=R T^{2}\left(\frac{\partial \ln q}{\partial T}\right)_{V}\]

The partition functions of translational, electronic, rotational,vibrational contribution are calculate as the equations list in(https://gaussian.com/thermo/). Then the entropy correction and internalthermal energy correction are calculated at specific temperature andpresure.

\[U=E_{t}+E_{e}+E_{v}+E_{r}\]

\[S=S_{t}+S_{e}+S_{v}+S_{r}\]

For linear molecules, degree of vibrational freedom is 3N - 5, VASPKITwill neglect smallest 5 frequencies. For non-linear molecules, degree ofvibrational freedom is 3N - 6, VASPKIT will neglect smallest 6frequencies.

VASPKIT also include zero-point energy (ZPE) in thermo energy correctionfrom the frequency calculation OUTCAR.

\[\varepsilon_{\mathrm{ZPE}}=\frac{h \nu}{2}\]

The Gibbs free energy (G) and enthalpy (H) include ∆PV = ∆NRT basedon ideal gas approximation.

\[H=U+k_{B} T\]

The Gibbs free energy G can be derived with the help of the totalentropy S,

\[G=H-T S\]

At last, VASPKIT will output:

1. Zero-point energy \(\varepsilon_{\mathrm{ZPE}}\);

2. Thermal correction to U(T), = \(\varepsilon_{\mathrm{ZPE}}\)+ \(\Delta U_{0 \rightarrow T}\)

3. Thermal correction to H(T), = \(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta U_{0 \rightarrow T}\) + PV =\(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta H_{0 \rightarrow T}\)

4. Thermal correction to G(T), = \(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta U_{0 \rightarrow T}\) + PV + TS =\(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta G_{0 \rightarrow T}\)

where \(\Delta U_{0 \rightarrow T}\),\(\Delta H_{0 \rightarrow T}\) , and\(\Delta G_{0 \rightarrow T}\) are interal energy, enthalpy, andGibbs free energy difference between 0 K and T K.

1 f = 46.979964 THz 295.183821 2PiTHz 1567.082879 cm-1 194.293584 meV2 f = 1.952595 THz 12.268518 2PiTHz 65.131565 cm-1 8.075288 meV3 f = 1.120777 THz 7.042052 2PiTHz 37.385108 cm-1 4.635164 meV4 f/i= 0.006984 THz 0.043884 2PiTHz 0.232971 cm-1 0.028885 meV5 f/i= 1.046106 THz 6.572876 2PiTHz 34.894327 cm-1 4.326347 meV6 f/i= 1.294104 THz 8.131095 2PiTHz 43.166663 cm-1 5.351986 meV

502 +-------------------------- Warm Tips --------------------------+ See An Example in vaspkit/examples/thermo_correction/ethanol. This Feature Was Contributed by Nan XU, Jincheng Liu and Sobereva. Included Vibrations, Translation, Rotation & Electron contributions. GAS molecules should not be with any fix. -->> (01) Reading Structural Parameters from CONTCAR File... -->> (02) Analyzing Molecular Symmetry Information... Molecular Symmetry is: D(inf)h Linear molecule found! +---------------------------------------------------------------+ Please input Temperature(K)!298.15 Please input Pressure(Atm)!1 Please input Spin multiplicity!--(Number of Unpaired electron + 1)3 ------------>> -->> (03) Extracting frequencies from OUTCAR... -->> (04) Reading OUTCAR File... -->> (05) Calculating Thermal Corrections... U(T) = ZPE + Delta_U(0->T) H(T) = U(T) + PV = ZPE + Delta_U(0->T) + PV G(T) = H(T) - TS = ZPE + Delta_U(0->T) + PV - TS Zero-point energy E_ZPE : 2.240 kcal/mol 0.097144 eV Thermal correction to U(T): 3.724 kcal/mol 0.161475 eV Thermal correction to H(T): 4.316 kcal/mol 0.187167 eV Thermal correction to G(T): -10.302 kcal/mol -0.446723 eV Entropy S : 205.141 J/(mol*K) 0.002126 eV

Adsorbed Molecular Free Energy Correction¶

Unlike gas molecules, adsorbed molecules form chemical bonds withsubstrate, which limits the translational and rotational freedom of theadsorbed molecules. So the contribution of translation and rotation toentropy and enthalpy is significantly reduced (so called hinderedtranslator / hindered rotor model). This does not mean no translationalor rotational contribution.

One common method is to attribute the translational or rotational partof the contribution to vibration, that is, the 3N vibrations of thesurface-adsorbing molecules (except the virtual frequency) are all usedto calculate the correction of the thermo energy.

The small the vibration frequencies have large contribution to entropy.It is very likely that a small vibration frequency will lead to abnormalentropy and free energy correction. So, it suggests that when the freeenergy of the surface adsorption molecule is corrected, the contributionof the frequency below 50 cm-1 is calculated as 50 cm-1, and VASPKITalso does this. And VASPKIT neglects the electron motion because of itssmall contribution.

And the PV term is also delete. At last VASPKIT will output:

1. Zero-point energy \(\varepsilon_{\mathrm{ZPE}}\);

2. Thermal correction to U(T) = H(T) =\(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta U_{0 \rightarrow T}\)

3. Thermal correction to G(T) = \(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta U_{0 \rightarrow T}\) + TS =\(\varepsilon_{\mathrm{ZPE}}\) +\(\Delta G_{0 \rightarrow T}\)

$ grep cm OUTCAR 1 f = 10.904836 THz 68.517103 2PiTHz 363.746154 cm-1 45.098792 meV 2 f = 10.827330 THz 68.030119 2PiTHz 361.160834 cm-1 44.778253 meV 3 f = 10.751668 THz 67.554721 2PiTHz 358.637024 cm-1 44.465340 meV

501 +-------------------------- Warm Tips --------------------------+ This Feature Was Contributed by Nan XU, Qiang LI and Jincheng LIU. See An Example in vaspkit/examples/thermo_correction/ORR. Only vibrations! No Translation & Rotation & Electron contributions. +---------------------------------------------------------------+ Please Enter The Temperature (K): ------------>>298.15 -->> (01) Reading OUTCAR File... +-------------------------- Summary ----------------------------+ Neglect PV contribution to translation for adsorbed molecules. To avoid abnormal entropy contribution, frequencies less than 50 cm-1 are set to 50 cm-1. U(T) = H(T) = ZPE + Delat_U(0->T) G(T) = H(T) - TS = ZPE + Delat_U(0->T) - TS Temperature (K): 298.1 Zero-point energy E_ZPE : 1.549 kcal/mol 0.067171 eV Thermal correction to U(T): 2.206 kcal/mol 0.095670 eV Thermal correction to H(T): 2.206 kcal/mol 0.095670 eV Thermal correction to G(T): 1.207 kcal/mol 0.052342 eV Entropy S : 14.021 J/(mol*K) 0.000145 eV

Determine elastic constants based on energy-strain method¶

Within the linear elastic region, the stress \(\boldsymbol{\sigma}\)response of solids to external loading strain\(\boldsymbol{\varepsilon}\) satisfies the generalized Hooke’s lawand can be simplified in the Voigt notation

\[\sigma_{i}=\sum_{j=1}^{6} c_{i j} \varepsilon_{j},\]

where the strain or stress are represented as a vector with 6independent components respectively, i.e., \(1 \leq i, j \leq 6\).\(C_{i j}\) is the second order elastic stiffness tensor expressedby a 6 × 6 symmetric matrix in units of GPa. The elastic stiffnesstensor \(C_{i j}\) can be determined based on the first-orderderivative of the stress-strain curves.

An alternative approach to calculate theoretical elastic constants isbased on the energy variation by applying small strains to theequilibrium lattice configuration. The elastic energy\(\Delta E\left(V,\left\{\varepsilon_{i}\right\}\right)\) of a solidunder strain in the harmonic approximation is given by

\[\Delta E\left(V,\left\{\varepsilon_{i}\right\}\right)=E\left(V,\left\{\varepsilon_{i}\right\}\right)-E\left(V_{0}, 0\right)=\frac{V_{0}}{2} \sum_{i, j=1}^{6} C_{i j} \varepsilon_{j} \varepsilon_{i},\]

where \(E\left(V,\left\{\varepsilon_{i}\right\}\right)\) and \(E\left(V_{0}, 0\right)\)are the total energies of the distorted and undistorted lattice cells,with the volume of \(V\) and \(V_0\) respectively. Theenergy-strain method corresponds that the elastic stiffness tensor isderived from the second-order derivative of the total energies versusstrain. In general, the stress-strain method require highercomputational precise to achieve the same accuracy as energy-strainmethod. Nevertheless, the former requires much smaller set ofdistortions than the latter. Considering that the energy-strain has lessstress sensitivity than the stress-strain one, the former method hasbeen implemented into the VASPKIT code.

The number of independent elastic constants depends on the symmetry of crystal. The lower the symmetry means the more the independent elastic constants. For example, the cubic crystals have three but the triclinic ones have 21 independent elastic constants. Next we take diamond as an example and calculate its elastic contants(vaspkit/examples/elastic/diamond_3D). Thera are three independentelastic constants for FCC crystal, \(\mathrm{C}_{11}\), \(\mathrm{C}_{12}\) and \(\mathrm{C}_{44}\). Thus, the elastic energy is given by

\[\begin{split}\Delta E\left(V,\left\{\varepsilon_{i}\right\}\right)=\frac{V_{0}}{2}\left[\begin{array}{cccccc}\varepsilon_{1} & \varepsilon_{2} & \varepsilon_{3} & \varepsilon_{4} & \varepsilon_{5} & \varepsilon_{6}\end{array}\right]\left[ \begin{array}{cccccc}{C_{11}} & {C_{12}} & {C_{12}} & {0} & {0} & {0} \\ {C_{12}} & {C_{11}} & {C_{12}} & {0} & {0} & {0} \\ {C_{12}} & {C_{12}} & {C_{11}} & {0} & {0} & {0} \\ {0} & {0} & {0} & {C_{44}} & {0} & {0} \\ {0} & {0} & {0} & {0} & {C_{44}} & {0} \\ {0} & {0} & {0} & {0} & {0} & {C_{44}}\end{array}\right]\left[ \begin{array}{l}{\varepsilon_{1}} \\ {\varepsilon_{2}} \\ {\varepsilon_{3}} \\ {\varepsilon_{4}} \\ {\varepsilon_{5}} \\ {\varepsilon_{6}}\end{array}\right]\end{split}\]

To determine the elastic constants for the cubic phase, we apply thetri-axial shear strain\(\varepsilon=(0,0,0, \delta, \delta, \delta)\), to the crystal. Then \(\Delta E=\frac{V_{0}}{2}\left(C_{44} \varepsilon_{4} \varepsilon_{4}+C_{44} \varepsilon_{5} \varepsilon_{5}+C_{44} \varepsilon_{6} \varepsilon_{6}\right)\)The \(C_{44}\) can be calculated from\(\frac{\Delta E}{V_0}=\frac{3}{2} C_{44} \delta^{2}\).Similarly, we appply \(\varepsilon=(\delta, \delta, 0,0,0,0)\) andget \(\Delta E=\frac{V}{2}\left(C_{11} \varepsilon_{1} \varepsilon_{1}+C_{11} \varepsilon_{2} \varepsilon_{2}+C_{12} \varepsilon_{1} \varepsilon_{2}+C_{12} \varepsilon_{2} \varepsilon_{1}\right)\)

Apply \(\varepsilon=(\delta, \delta, \delta, 0,0,0)\) to get\(\frac{\Delta E}{V_0}=\frac{3}{2}\left(C_{11}+2 C_{12}\right) \delta^{2}\)。

Finally, math::` mathrm{C}_{11}`, math::` mathrm{C}_{12}` and math::` mathrm{C}_{44}` can be calcauclated by solving these equations.

The relationship of lattice vectors between distorted and fully relaxedlattice cells is:

\[\begin{split}\left( \begin{array}{l}{\boldsymbol{a}^{\prime}} \\ {\boldsymbol{b}^{\prime}} \\ {\boldsymbol{c}^{\prime}}\end{array}\right)=\left( \begin{array}{l}{\boldsymbol{a}} \\ {\boldsymbol{b}} \\ {\boldsymbol{c}}\end{array}\right) \cdot(\boldsymbol{I}+\boldsymbol{\varepsilon})\end{split}\]

where math::` boldsymbol{I}` is the identity matrix, The strain tensor

\(\boldsymbol{\varepsilon}\) is defined by

\[\begin{split}\boldsymbol{\varepsilon}=\left( \begin{array}{ccc}{\varepsilon_{1}} & {\frac{\varepsilon_{6}}{2}} & {\frac{\varepsilon_{5}}{2}} \\ {\frac{\varepsilon_{6}}{2}} & {\varepsilon_{2}} & {\frac{\varepsilon_{4}}{2}} \\ {\frac{\varepsilon_{5}}{2}} & {\frac{\varepsilon_{4}}{2}} & {\varepsilon_{3}}\end{array}\right)\end{split}\]

Fisrt Step：Prepare the following files：

1. Prepare the fully-relaxed POSCAR file including the standardconventional cell of diamond. You can get the standard conventionalcell by running vaspkit with task of 603/604;

2. Run vaspkit with task of 102 to generate KPOINTS with fine k-mesh.

3. Prepare INCAR as follows.

Global Parameters ISTART = 0 LREAL = F PREC = High LWAVE = F LCHARG = F ADDGRID= .TRUE.Electronic Relaxation ISMEAR = 0 SIGMA = 0.05 NELM = 40 NELMIN = 4 EDIFF = 1E-08Ionic Relaxation NELMIN = 6 NSW = 100 IBRION = 2 ISIF = 2 （Must be 2） EDIFFG = -1E-02

4，Prepare VPKIT.in file and set the value of first line to 1 (1 means pre-processing)；

1 ! 1 for pre-processing; 2 for post-processing3D ! 2D for two-dimentional, 3D for bulk7 ! number of strain case-0.015 -0.010 -0.005 0.000 0.005 0.010 0.015 ! Strain range

and run vaspkit-201

 ------------>>201 -->> (01) Reading VPKIT.in File... +-------------------------- Warm Tips --------------------------+ See an example in vaspkit/examples/elastic/diamond_3D, Require the fully-relaxed and standardized Convertional cell. +---------------------------------------------------------------+ -->> (02) Reading Structural Parameters from POSCAR File... -> C44 folder created successfully! -> strain_-0.015 folder created successfully! -> strain_-0.010 folder created successfully! -> strain_-0.005 folder created successfully! -> strain_0.000 folder created successfully! -> strain_+0.005 folder created successfully! -> strain_+0.010 folder created successfully! -> strain_+0.015 folder created successfully! -> C11_C12_I folder created successfully! -> strain_-0.015 folder created successfully! -> strain_-0.010 folder created successfully! -> strain_-0.005 folder created successfully! -> strain_0.000 folder created successfully! -> strain_+0.005 folder created successfully! -> strain_+0.010 folder created successfully! -> strain_+0.015 folder created successfully! -> C11_C12_II folder created successfully! -> strain_-0.015 folder created successfully! -> strain_-0.010 folder created successfully! -> strain_-0.005 folder created successfully! -> strain_0.000 folder created successfully! -> strain_+0.005 folder created successfully! -> strain_+0.010 folder created successfully! -> strain_+0.015 folder created successfully!

Second step: Batch performing DFT calucaltions using VASP code；

Third step: Modify the value of the first line in VPKIT.in file to 2(2 means post-processing); and run vaspkit-201 again.You will getthe following information,

 ------------>>201 -->> (01) Reading VPKIT.in File.. \\ =o) (o> +------------------------------------------------------+/\\_(()_|VPKIT.in will be renamed to INPUT.in in future version| \_V // +------------------------------------------------------+ \\ -->> (01) Reading VPKIT.in File...+-------------------------- Warm Tips --------------------------+ See some examples in vaspkit/examples/elastic, Require the fully-relaxed and standardized Conventional cell.+---------------------------------------------------------------+ -->> (02) Reading Structural Parameters from POSCAR File... -->> (03) Calculating Fitting Coefficients of Energy vs. Strain. -->> Current directory: Fitting Precision C44 Folder: 0.817E-09 C11_C12_I Folder: 0.814E-08 C11_C12_II Folder: 0.135E-07+-------------------------- Summary ----------------------------+Crystal Class: m-3mSpace Group: Fd-3mCrystal System: Cubic systemIncluding Point group classes: 23, 2/m-3, 432, -43m, 4/m-32/mThere are 3 independent elastic constants C11 C12 C12 0 0 0 C12 C11 C12 0 0 0 C12 C12 C11 0 0 0 0 0 0 C44 0 0 0 0 0 0 C44 0 0 0 0 0 0 C44Stiffness Tensor C_ij (in GPa): 1050.640 126.640 126.640 0.000 0.000 0.000 126.640 1050.640 126.640 0.000 0.000 0.000 126.640 126.640 1050.640 0.000 0.000 0.000 0.000 0.000 0.000 559.861 0.000 0.000 0.000 0.000 0.000 0.000 559.861 0.000 0.000 0.000 0.000 0.000 0.000 559.861Compliance Tensor S_ij (in GPa^{-1}): 0.000977 -0.000105 -0.000105 0.000000 0.000000 0.000000 -0.000105 0.000977 -0.000105 0.000000 0.000000 0.000000 -0.000105 -0.000105 0.000977 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001786 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001786 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001786Anisotropic mechanical properties of bulk singlecrystal:+---------------------------------------------------------------+| Mechanical Properties | Min | Max | Anisotropy|+---------------------------------------------------------------+| Bulk Modulus B (GPa) | 434.640 | 434.640 | 1.000 || Young's Modulus E (GPa) | 1023.395 | 1175.046 | 1.148 || Shear Modulus G (GPa) | 462.000 | 559.861 | 1.212 || Poisson's Ratio v | 0.012 | 0.119 | 9.987 || Linear compressibility | 0.767 | 0.767 | 1.000 |+---------------------------------------------------------------+Linear compressibility beta is in unit of TPa^-1References:[1] Marmier A, Comput. Phys. Commun. 181, 2102–2115 (2010)[2] Gaillac R, J. Phys. Condens. Matter 28, 275201 (2016)Average mechanical properties of bulk polycrystal:+---------------------------------------------------------------+| Mechanical Properties | Voigt | Reuss | Hill |+---------------------------------------------------------------+| Bulk Modulus B (GPa) | 434.64 | 434.640 | 434.640 || Young's Modulus E (GPa) | 1116.34 | 1109.298 | 1112.824 || Shear Modulus G (GPa) | 520.72 | 516.130 | 518.423 || Poisson's Ratio v | 0.07 | 0.075 | 0.073 || P-wave Modulus (GPa) | 1128.93 | 1122.814 | 1125.871 || Pugh's Ratio (B/G) | 0.83 | 0.842 | 0.838 || Vickers Hardness (GPa)[6] | 91.69 | 90.237 | 90.962 || Vickers Hardness (GPa)[7] | 94.70 | 93.170 | 93.935 |+---------------------------------------------------------------+Pugh's Ratio (B/G): 0.84 --> Brittle region (< 1.75)Cauchy Pressure Pc (GPa): -433.2 --> Metallic-like bonding (< 0)Kleinman’s parameter: 0.29 --> Bonding Bending dominated (< 0.5)Universal Elastic Anisotropy: 0.04Chung-Buessem Anisotropy: 0.0Isotropic Poisson's Ratio: 0.07Longitudinal wave velocity (m/s): 17945.173Transverse wave velocity (m/s): 12177.146Average wave velocity (m/s): 13280.911Debye temperature (K): 2212.9References:[1] Voigt W, Lehrbuch der Kristallphysik (1928)[2] Reuss A, Z. Angew. Math. Mech. 9 49-58 (1929)[3] Hill R, Proc. Phys. Soc. A 65 349-54 (1952)[4] Debye temperature J. Phys. Chem. Solids 24, 909-917 (1963)[5] Elastic wave velocities calculated using Navier's equation[6] Chen X-Q, Intermetallics 19, 1275 (2011)[7] Tian Y-J, Int. J. Refract. Hard Met. 33, 93–106 (2012)Eigenvalues of the stiffness matrix (in GPa):Eigenvalue lamda_1 = 559.861 > 0 meeted.Eigenvalue lamda_2 = 559.861 > 0 meeted.Eigenvalue lamda_3 = 559.861 > 0 meeted.Eigenvalue lamda_4 = 924.000 > 0 meeted.Eigenvalue lamda_5 = 924.000 > 0 meeted.Eigenvalue lamda_6 = 1303.921 > 0 meeted.Elastic stability criteria as seen in PRB 90, 224104 (2014).Criteria (i) C11 - C12 > 0 meeted.Criteria (ii) C11 + 2C12 > 0 meeted.Criteria (iii) C44 > 0 meeted.This Structure is Mechanically Stable.+---------------------------------------------------------------+ -->> (04) Written ELASTIC_TENSOR File!

The calcualted elastic contants are well agreement with the availableexperimental values, \(C_{11}=1079\) GPa, \(C_{12}=124\) GPa and\(C_{44}=578\) GPa. VASPKIT can also read elastic constants fromthe OUTCAR (task 203) or ELASTIC_TENSOR.in (task 202) file and determine various mechanicalproperites and elastic stability criterion. For more details, seevaspkit/examples/elastic.

Linear Optical Properties¶

The linear optical properties can be obtained from thefrequency-dependent complex dielectric function\(\varepsilon(\omega)\)

\[\varepsilon(\omega)=\varepsilon_{1}(\omega)+i \varepsilon_{2}(\omega),\]

where \(\varepsilon_{1}(\omega)\) and\(\varepsilon_{2}(\omega)\) are the real and imaginary parts of thedielectric function, and \(\omega\) is the photon frequency.

The frequency-dependent linear optical spectra, e.g., refractive index \(n(\omega)\),extinction coefficient \(\kappa(\omega)\),absorption coefficient \(\alpha(\omega)\),energy-loss function \(L(\omega)\),reflectivity \(R(\omega)\) can be calculated from the real \(\varepsilon_{1}(\omega)\) and \(\varepsilon_{2}(\omega)\) parts [See Ref. A. M. Fox, Optical Properties of Solids]:

\[n(\omega)=\left[\frac{\sqrt{\varepsilon_{1}^{2}+\varepsilon_{2}^{2}}+\varepsilon_{1}}{2}\right]^{\frac{1}{2}}\]

\[k(\omega)=\left[\frac{\sqrt{\varepsilon_{1}^{2}+\varepsilon_{2}^{2}}-\varepsilon_{1}}{2}\right]^{\frac{1}{2}}\]

\[\alpha(\omega)=\frac{\sqrt{2} \omega}{c}\left[\sqrt{\varepsilon_{1}^{2}+\varepsilon_{2}^{2}}-\varepsilon_{1}\right]^{\frac{1}{2}}\]

\[L(\omega)=\operatorname{Im}\left(\frac{-1}{\varepsilon(\omega)}\right)=\frac{\varepsilon_{2}}{\varepsilon_{1}^{2}+\varepsilon_{2}^{2}}\]

\[R(\omega)=\frac{(n-1)^{2}+k^{2}}{(n+1)^{2}+k^{2}}\]

In VASPKIT Ver. 1.00 or later, it is not necessary to run optical.sh toextract image and real parts of dielectric function from vasprun.xml atfirst step (optional). You can run vaspkit -task 711 to get the theselinear optical spectrums at one shot. The program will read dielectricfunction from vasprun.xml directely if both REAL.in and IMAG.in filesnot exist. We take Si as an example and calcuculate its opticalproperties at GW0+BSE leve, as shown in figure below.

 ------------>>71 ============================ Optical Options ==================== 711) Linear Optical Spectrums 712) Transition Dipole Moment at Single kpoint 713) Transition Dipole Moment for DFT Band-Structure 0) Quit 9) Back ------------>>711 +-------------------------- Warm Tips --------------------------+ See an example in vaspkit/examples/Si_bse_optical. This module is NOT suitable for low-dimensional materials. +---------------------------------------------------------------+ ===================== Energy Unit =============================== Which Energy Unit Do You Want to Adopt? 1) eV 2) nm 3) THz ------------>>1 -->> (01) Reading Dielectric Function From vasprun.xml File... -->> (02) Reading IMAG.in and REAL.in Files... -->> (03) Written Optical Files Succesfully!

Transition Dipole Moment¶

The transition dipole moment (TDM) or transition moment, usually denotedfor a transition between an initial state a and a final state b, isthe electric dipole moment associated with the transition between thetwo states. In general the TDM is a complex vector quantity thatincludes the phase factors associated with the two states. Its directiongives the polarization of the transition, which determines how thesystem will interact with an electromagnetic wave of a givenpolarization, while the square of the magnitude gives the strength ofthe interaction due to the distribution of charge within the system. TheSI unit of the transition dipole moment is the Coulomb-meter (Cm); amore conveniently sized unit is the Debye (D). The transitionprobabilities between the valence and the conduction band are revealedby the calculated sum of the squares of TDM (in unit of Debey\(^2\)).

We task cubic CsPbI\(_3\) perovskite phase as an example andpresent the PBE calculated band structures and corresponding TDM infigure below. This system has 44 valence electrons in all. Next wecalcualte the squares of transition dipole moment from the highestvalence band (band index 22) to the lowest conduction band (band index23). It is found the results are well agreement with the previousfinding (see Ref. J. Phys. Chem. Lett. 2017, 8, 2999−3007).

The calculating process is very similar with the band structurecalcuation except for setting LWAVE= .TRUE. in the INCAR file. Moredetails on these calculations are given in vaspkit/examples/tdm.

 ------------>>71 ============================ Optical Options ==================== 711) Linear Optical Spectrums 712) Transition Dipole Moment at Single kpoint 713) Transition Dipole Moment for DFT Band-Structure 0) Quit 9) Back ------------>>713 +-------------------------- Warm Tips --------------------------+ See an example in vaspkit/examples/tdm. ONLY Support KPOINTS created by VASPKIT for hybrid bandstructure. +---------------------------------------------------------------+ ======================= Band-Structure Mode ==================== 1) DFT Band-Structure Calculation 2) Hybrid-DFT Band-Structure Calculation ------------->>1 Please input TWO bands separated by spaces. (e.g., 4 5 means to get TDM between 4 and 5 bands at all kpts) ------------>>22 23 +-------------------------- Warm Tips --------------------------+ Current Version Only Support the Stardard Version of VASP code. +---------------------------------------------------------------+ -->> (01) Reading the header information in WAVECAR file... +--------------------- WAVECAR Header --------------------------+ SPIN = 1 NKPTS = 80 NBANDS = 36 ENCUT = 400.00 Coefficients Precision: Complex*8 Maximum number of G values: GX = 11, GY = 11, GZ = 11 Estimated maximum number of plane-waves: 5575 +---------------------------------------------------------------+ -->> (02) Start to read WAVECAR file, take your time ^.^ Percentage complete: 25.0% Percentage complete: 50.0% Percentage complete: 75.0% Percentage complete: 100.0% -->> (03) Reading K-Paths From KPOINTS File... -->> (04) Start to calculate Transition Dipole Moment... Percentage complete: 25.0% Percentage complete: 50.0% Percentage complete: 75.0% Percentage complete: 100.0% -->> (05) Written TDM.dat File!

If one want to calculate the square of TDM between two state at singlek-point. Just run vaspkit with task 712, for example,

71 ============================ Optical Options ==================== 711) Linear Optical Spectrums 712) Transition Dipole Moment at Single kpoint 713) Transition Dipole Moment for DFT Band-Structure 0) Quit 9) Back ------------>>712 +-------------------------- Warm Tips --------------------------+ See an example in vaspkit/examples/tdm. Please Set LWAVE= .TRUE. in the INCAR file. +---------------------------------------------------------------+ Input ONE kpoint and TWO bands separated by spaces. (e.g., 1 4 5 means to get TDM between 4 and 5 band at 1st kpt) ------------>>1 22 23 +-------------------------- Warm Tips --------------------------+ Current Version Only Support the Stardard Version of VASP code. +---------------------------------------------------------------+ -->> (01) Reading the header information in WAVECAR file... +--------------------- WAVECAR Header --------------------------+ SPIN = 1 NKPTS = 80 NBANDS = 36 ENCUT = 400.00 Coefficients Precision: Complex*8 Maximum number of G values: GX = 11, GY = 11, GZ = 11 Estimated maximum number of plane-waves: 5575 +---------------------------------------------------------------+ Square of TDM (Debye^2): X Y Z Total 0.000 201.617 0.000 201.617

Cleary, the calculate total-square of TDM at \(\text{X}\) point forcubic CsPbI\(_3\) is 201.62 Debey\(^2\).

Redefine Lattice¶

VASPKIT option 400 can build supercell and rotate the latticevector. It requires three coefficients for each new vector by followingequation:

\[\begin{split}\left(\begin{array}{l}{A} \\ {B} \\ {C}\end{array}\right)=\left(\begin{array}{lll}{c_1} & {c_2} & {c_3} \\ {c_4} & {c_5} & {c_6} \\ {c_7} & {c_8} & {c_9}\end{array}\right)\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)\end{split}\]

A, B, and C are new lattice vectors, a, b, and care old lattice vectors.

Example: Build \(\operatorname{Au}(111)-(\sqrt{3} \times \sqrt{3}) R 30^{\circ}\) slab from primitive slab

Build a supercell with its length of new lattice vector equal to\(\sqrt{3}\)a.

400 -->> (01) Reading Structural Parameters from POSCAR File... Enter the new lattice verctor a in terms of old: (MUST be three integers, e.g., 1 2 3)1 -1 0 Enter the new lattice verctor b in terms of old:2 1 0 Enter the new lattice verctor c in terms of old:0 0 1 +-------------------------- Summary ----------------------------+ The Transformation Matrix P is: 1 -1 0 2 1 0 0 0 1 Lattice Constants in New Cell: 4.995 4.995 17.064 Lattice Angles in New Cell: 90.00 90.00 60.00 Total Atoms in New Cell: Volume of New Cell is 3 times of the Old Cell

Build Supercell¶

VASPKIT 401 can build supercell by three number C1, C2, C3.

\[\begin{split}\left(\begin{array}{l}{A} \\ {B} \\ {C}\end{array}\right)=\left(\begin{array}{lll}{c_{1}} & {c_{2}} & {c_{3}}\end{array}\right)\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)\end{split}\]

For example, build supercell with length is 2 times in the a andb direction and keep the c vector.

401 ===================== Structural File =========================== 1) POSCAR 2) CONTCAR 3) The Specified Filename 0) Quit 9) Back ------------>>1 -->> (01) Reading Structural Parameters from POSCAR File... +-------------------------- Warm Tips --------------------------+ Enter the repeated unit along a, b and c directions with space! (MUST be three integers, e.g., 1 2 3) ------------>>2 2 1 -->> (02) Written SC221.vasp File

Fix Atoms by Layers¶

VASP default structural optimization allows all atoms to move freely inall directions (x, y, z). In some cases, structural optimizationcalculations require to fix some atoms. For example, surfacecalculations with slab model need to fix some bottom atoms.

To fix atoms, the POSCAR seventh line should be switched toselective dynamics (only the first character is relevant and must beS or s). This mode allows to provide extra flags for each atomsignaling whether the respective coordinate(s) of this atom will beallowed to change during the ionic relaxation. F means fix, andT means not fix.

VASPKIT provides options 402 to fix atoms by layers:

1. Input the file name that we want to fix. 1 for POSCAR, 2 forCONTCAR, 3 for specified name but should also be POSCAR format.

2. Input a threshold (unit: Å) to separate layers. Only atomic layers’z coordinates difference is bigger than the threshold, those layerscan be recognized. Otherwise, it will be recognized to be one layer.VASPKIT shows some suggested threshold and its corresponding layersnumber.

Then VASPKIT will print the recognized layers number.

3. Choose the layers that we want to fix. VASPKIT will automiticallyfix those bottom layers and output a POSCAR_FIX file. BeforeVASP calcualtion, cp POSCAR_FIX POSCAR.

 402 ===================== Structural File =========================== 1) POSCAR 2) CONTCAR 3) The Specified Filename 0) Quit 9) Back ------------>>1 -->> (01) Reading Structural Parameters from POSCAR File... Threshold: 0.3 layers: 4 Threshold: 0.6 layers: 4 Threshold: 0.9 layers: 4 Threshold: 1.2 layers: 4 Threshold: 1.5 layers: 4 Please choose a threshold to separate layers->1 Found 4 layers, choose how many layers to be fixed=>2Fixed atoms are: 1 4 -->> (02) Written POSCAR_FIX File!

Fix Atoms by Height¶

VASPKIT 403 can also fix atoms, but by a range of z coordinates.Similar as the option 402, first select a structural file, and theninput two numbers to set the z coordinate range. Then, select the outputformat Fractional or Cartesian. it will output POSCAR_FIXfile.

If the (1) Fractional coordinates is select, the inputz_min z_max numbers should be the z fractional coordination and theoutput file POSCAR_FIX also is written as fractional style.

If the (2) Cartesian coordinates is select, the inputz_min z_max numbers should be the z Cartesian coordination and theoutput file POSCAR_FIX also is written as Cartesian style.

403 ===================== Structural File =========================== 1) POSCAR 2) CONTCAR 3) The Specified Filename 0) Quit 9) Back ------------>>1 -->> (01) Reading Structural Parameters from POSCAR File... +---------------------------------------------------------------+ | Selective Dynamics is Activated! | +---------------------------------------------------------------+ Atoms between your chosen section in z direction will be fixed Type z_min z_max, Note: z_min < z_max0 5 Read and Write atomic coordinates as: (1) Fractional coordinates (2) Cartesian coordinates2Fixed atoms are: 1 4 -->> (02) Written POSCAR_FIX File!

Swap Axis of Lattice Vector¶

VASPKIT 407 can do swap two axis of lattice vector.

Symmetry Tools¶

VASPKIT can identify symmetry of crystal from POSCAR or CONTCAR by601 or608, respectively.

For example \(\theta\)-Al2O3:

601 -->> (01) Reading Structural Parameters from POSCAR File... +-------------------------- Summary ----------------------------+ Prototype: A2B3 Total Atoms: 10 Formula Unit: Al2O3 [ Alphabetically Listed: Al2O3 ] Full Formula Unit: Al4O6 Crystal System: Monoclinic Crystal Class: 2/m Bravais Lattice: mC Lattice Constants: 6.142 6.142 5.671 Lattice Angles: 76.387 103.613 152.294 Volume: 96.475 Density (g/cm3): 3.510 Space Group: 12 Point Group: 5 [ C2h ] International: C2/m Symmetry Operations: 4

If POSCAR is conventional cell, VASPKIT 602 can transfer it toprimitive Cell.

If POSCAR is primitive cell, VASPKIT 602 can transfer it toconventional Cell.

VASPKIT can also identify point group of molecule from POSCAR 609,

For example: point group of O2 molecule is \(D_{\infty h}\):

6 =================== Symmetry Options ============================ 601) Symmetry Finding for Crystals 602) Primitive Cell Finding 603) Conventional Cell Finding (experimental) 604) Standardize Crystal Cell 608) Summary for Relaxed-Structure 609) Symmetry Finding for Molecules and Clusters 0) Quit 9) Back ------------>>609 -->> (01) Reading Structural Parameters from CONTCAR File... -->> (02) Analyzing Molecular Symmetry Information... Molecular Symmetry is: D(inf)h

Planar Average Potential¶

Planar average for grid files is one of the basic function of VASPKIT.It supports all grid files of VASP output, such as CHGCAR, PARCHG,LOCPOT, ELFCAR. All such grid files are written by following commands inFortran:

WRITE(IU,FORM) (((C(NX,NY,NZ),NX=1,NGXC),NY=1,NGYZ),NZ=1,NGZC)

For one selected direction, such as z direction. VASPKIT do an averagein XY plane on each NZ,

\[\sum_{i j} \Delta x_{i} \Delta y_{j} \rho_{i, j} / \sum_{i j} \Delta x_{i} \Delta y_{j}\]

The output POTPAVG.dat contains two columns. First column is zcoordinate (in Å), Second column is Planar Average-Potential (in eV) orDensitiy (in e).

Here, take LOCPOT as an example. It contains total local potential (ineV) when LVTOT = .TRUE. exists in the INCAR file, or containselectrostatic potential (in eV) when LVHAR = .TRUE. exists in theINCAR file. Electrostatic potential is desirable for the evaluation ofthe work-function, because the electrostatic potential converges morerapidly to the vacuum level than the total potential. To get thework-function, potential along z direction (i.e.slab normal direction)must be averaged in the slab planes.

Hamiltonian of Kohn-Sham is:

\[\hat{H}=-\frac{1}{2} \nabla^{2}-\sum_{A} \frac{Z_{A}}{\left|\boldsymbol{r}-\boldsymbol{R}_{A}\right|}+\int_{\infty} \frac{\rho\left(\boldsymbol{r}^{\prime}\right)}{\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|} d \boldsymbol{r}^{\prime}+v_{x c}[\rho(\boldsymbol{r})]\]

Electrostatic potential is the summation of second (ionic) and third(hartree) parts of Hamiltonian .

##### initial parameters I/O #####SYSTEM = AuNCORE = 5ISTART = 1ICHARG = 1LWAVE = .TRUE.LCHARG = .TRUE.LVTOT = .FALSE.LVHAR = .TRUE.LELF = .FALSE.#### Electronic Relaxation ####ENCUT = 400ISMEAR = 1SIGMA = 0.2EDIFF = 1E-6NELMIN = 5NELM = 300GGA = PELREAL = AutoIDIPOL = 3

426 +-------------------------- Warm Tips --------------------------+ You MUST Know What You are Doing Check Convergence of Planar Average Value on Vacuum-Thickness! +---------------------------------------------------------------+ ===================== Specified Direction ======================= Which Direction is for the Planar Average? 1) x Direction 2) y Direction 3) z Direction 0) Quit 9) Back ------------>>3 -->> (1) Reading Structural Parameters from LOCPOT File... -->> (2) Reading Local Potential From LOCPOT File... -->> (3) Written POTPAVG.dat File! +-------------------------- Warm Tips --------------------------+ Check the Convergence of Vacuum-Level Versus Vacuum Thickness! +---------------------------------------------------------------+ +-------------------------- Summary ----------------------------+ Vacuum-Level (eV): 6.695 +---------------------------------------------------------------+

# Planar Distance (in Angstrom), Planar Average-Potential (in eV) or -Densitiy (in e) 0.0000 0.27377302E+01 0.1025 0.94832051E+00 0.2049 -0.13842911E+01 0.3074 -0.40872693E+01 0.4099 -0.69356914E+01 0.5123 -0.96882699E+01 0.6148 -0.12150765E+02 0.7172 -0.14164741E+02 ...

grep E-fermi OUTCARE-fermi : 1.5199 XC(G=0): -6.7056 alpha+bet : -6.4642

Macroscopic Averaged Potential¶

Besides planar averaged potential , VASPKIT can also calculatemacroscopic-averaged potential, which is obtained by converging theplanar averaged potential with a moving average. e.g.If one want to doplanar-averaged potential in z-direction, use z = average(x, y) forevery z. If one want to do macroscopic-averaged potential in z-dirction,use z = average(x, y, z ± \(\Delta\)z’) for every z, where2\(\Delta\)z’ is the layer distance for repeat structures.

For example: MoS2/WS2 heterogeneous junction.

1. Do a single-point calculation with LVHAR = .TRUE.,

2. Run VASPKIT 427 to obtain planar averaged potential andmacroscopic averaged potential from LOCPOT.

1. Select which direction to do planar average?

2. Input the period length to calculate macroscopic average. (Toaccelerate the convergence, this value should be close to and smallerthan the layer distance)

3. Input iteration number to get smoother macroscopic curve.

427 +-------------------------- Warm Tips --------------------------+ You MUST Know What You are Doing Check Convergence of Planar Average Value on Vacuum-Thickness! +---------------------------------------------------------------+ ===================== Specified Direction ======================= Which Direction is for the Planar Average? 1) x Direction 2) y Direction 3) z Direction 0) Quit 9) Back ------------>>1 -->> (01) Reading Structural Parameters from LOCPOT File... -->> (02) Reading Local Potential From LOCPOT File... Please input the period length to calculate macroscopic average: (Typical value: the minimum length of repeated unit.) ------------>>1.5 Please input iteration number to get smoother macroscopic curve: (Must be integer. You can gradually increase from 1 to 3.) ------------>>3 -->> (03) Written PLANAR_AVERAGE.dat File! +-------------------------- Summary ----------------------------+ Maximum of macroscopic average: 0.099 Minimum of macroscopic average: -0.101 +---------------------------------------------------------------+ -->> (04) Written MACROSCOPIC_AVERAGE.dat File!

Step 3. Draw graphs from the output MACROSCOPIC_AVERAGE.dat, andPLANAR_AVERAGE.dat

Charge Density and Spin Density¶

When use spin-polarized parameter (SPIN = 2), the output CHGCAR willcontain charge density and spin density. VASPKIT can extract the chargedensity and spin density by options 311 and 312 respectively.Outputs are saved in CHARGE.vasp and SPIN.vasp.

VASPKIT can also calculate the Spin-Up & Spin-Down density by option313. The outputs are saved in SPIN_UP.vasp & SPIN_DW.vaspfiles

Charge Density Difference¶

VASPKIT is powerful for post-processing grid file from VASP, such asCHGCAR, ELFCAR, LOCPOT, PARCHG.

CHGCAR contains the lattice vectors, atomic coordinates, the totalcharge density multiplied by the volume on the fine FFT-grid(NG(X,Y,Z)F), and the PAW one-center occupancies. By subtract two ormore CHGCAR, One can get the charge density difference.

Note: From electrostatic potential (LOCPOT), electron localizationfunction (ELFCAR), and partial charge (PARCHG) density differencecalculation can also be accomplished by the same method of VASPKIT.

Charge density difference is one of the important way to investigateelectronic structure. It can intuitively get the electron flow from theinteraction of two fragments. And explore the essence of chemical bond.There are several forms of charge density difference.：

（1）The charge density of the total system subtract the density of twoor more segments that make up it：

\[\Delta \rho=\rho_{A B}-\rho_{A}-\rho_{B}\]

（2）Charge density difference before and after self-consistentcalculation. Also known as deformation charge density:

\[\Delta \rho=\rho\left(A B_{\text { self-consistent }}\right)-\rho\left(A B_{\text { atomic }}\right)\]

（3）Density in one electronic state subtracts density in another. Forexample, the charge density under an applied electric field subtractsthe charge density without an external electric field. Another exampleis the density of the excited state minus the density of the groundstate.

\[\Delta \rho=\rho\left(A B_{\text { state 1 }}\right)-\rho\left(A B_{\text { state 0 }}\right)\]

Among the three charge density differences mentioned above, no matterwhich calculation, must ensure that the cell parameters andcoordinates of atoms must be consistent.

314 ======================= File Options ============================ Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR) ------------>>./CHGCAR ./co/CHGCAR ./slab/CHGCAR -->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./co/CHGCAR File... -->> (04) Reading Charge Density From ./co/CHGCAR File... -->> (05) Reading Structural Parameters from ./slab/CHGCAR File... -->> (06) Reading Charge Density From ./slab/CHGCAR File... -->> (07) Written CHGDIFF.vasp File!

314 ======================= File Options ============================ Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR) ------------>>./CHGCAR ./deform/CHGCAR -->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./deform/CHGCAR File... -->> (04) Reading Charge Density From ./deform/CHGCAR File... -->> (05) Written CHGDIFF.vasp File!

EFIELD = 0.05IDIPOL = 3LDIPOL = .TRUE.

WaveFunction Visualization in Real-Space¶

VASPKIT can extract the planewave coefficients of Kohn-Sham (KS) orbitalfrom the WAVECAR file and output the real-space wavefunction. Users needto input which K point to plot, and which band to plot.

Note: Now, VASPKIT can only output wavefunction for specified oneK-point and one band at one time. To sum several K points or seting theenergy range, the partial charge density calculation in the VASP can doit.

 Which K-POINT do you want to plot? (1<= ikpt <=1) ------------>>1

511 +-------------------------- Warm Tips --------------------------+ Open Real-Space WaveFunction Files with VESTA/VMD Package. +---------------------------------------------------------------+ Which K-POINT do you want to plot? (1<= ikpt <=1) ------------>>1 Which BAND do you want to plot? (1<= iband <=36) ------------>>5 +-------------------------- Warm Tips --------------------------+ Current Version Only Support the Stardard Version of VASP code. +---------------------------------------------------------------+ -->> (01) Reading the header information in WAVECAR file... +--------------------- WAVECAR Header --------------------------+ SPIN = 2 NKPTS = 1 NBANDS = 36 ENCUT = 400.00 Coefficients Precision: Complex*8 Maximum number of G values: GX = 25, GY = 25, GZ = 25 Estimated maximum number of plane-waves: 65450 +---------------------------------------------------------------+ -->> (02) Start to Post-Process Wavefunctions... -->> (03) Reading Structural Parameters from POSCAR File... -->> (04) Written RWAV_B0005_K0001_UP.vasp File! -->> (05) Written IWAV_B0005_K0001_UP.vasp File! -->> (06) Written RWAV_B0005_K0001_DW.vasp File! -->> (07) Written IWAV_B0005_K0001_DW.vasp File!Note: Now, VASPKIT can only output wavefunction for specified one K points and one band at one time. To sum several K points or seting the energy range, VASP partial charge calculation can do it.

... 287 -2.799024 1.000000 288 -2.799005 1.000000 289 -0.201544 0.000000 290 -0.201541 0.000000...

Molecular Dynamics process¶

721 calculate the Mean Squared Displacement (MSD) based on VASP MDresults. First prepare a reference structure file POSCAR.ref. Thenrun VASPKIT 721 to calculate MSD. The output MSD.dat containsdisplacement information |x|, |y|, |z|, and MSD, RMSD.

# ion_step |dx| |dy| |dz| MSD sqrt(MSD) 1 0.0000016 0.0000011 0.0000019 0.0000000 0.0000001 2 0.1724941 0.2677518 0.2270937 0.0001325 0.0115120 3 0.3367914 0.5283185 0.4461284 0.0005124 0.0226364 4 0.4918222 0.7752755 0.6499894 0.0010922 0.0330486 5 0.6370632 1.0074569 0.8382370 0.0018188 0.0426475 6 0.7959320 1.2340219 1.0093023 0.0026613 0.0515876 7 0.9545851 1.4665808 1.1630884 0.0036234 0.0601945 8 1.1064692 1.6957629 1.3186163 0.0047342 0.0688057 9 1.2477962 1.9293770 1.4885940 0.0060285 0.0776432

ATOM_DISPLACEMENT.dat contains displacement and RMSD information foreach atoms.

722 calculate the Radial Distribution Function (RDF) also based onVASP MD results.

File Format Transfer¶

USER interface¶

Cautionary Notes of VASP Wavefunctions¶

1. The calculated plane-wave coefficients of pseudo periodic part (Bloch function) are read from the given WAVECAR file and the corresponding G-vectors are generated using the algorithm developed in WaveTrans code. It should be noted that the pseudopotential augmentation is not included in the WAVECAR file. As a result, some caution should be exercised when deriving value from this information. For example, the partial charge density obtained using vaspkit-515 will differ from the PARCHG from VASP, but the qualitative shape of the charge density should match.

2. The actual number of plane waves (NPWS) of each k value is determined by the maximum plane-wave energy it is used for the VASP run, i.e. all G values for which \(\frac{ ℏ^2|k+G|^2 }{2m}\) is less than the cut-off energy. Unfortunately, the imprecisely-known scaling factor \(\frac{2m}{ℏ^2}\) leads to a discrepancy between the NPWS computed by VASPKIT and that contained in the WAVECAR file occasionally. If so, you will encounter an error message saying that “Error: the calculated NPWS is not equal to the read NPWS”. Then please adjust the final decimal place (or one beyond that) of FACTOR_ENCUT2NPWS (default value: 0.26246583) in the ~/.vaspkit file (1.3 and newer version). If such an adjustment is successfully determined, please contact wangvei@icloud.com with the new value, so that it can be incorporated into a future version of VASPKIT.

Usage FAQs¶

1. The fermi energy is set to zero in the DOS.dat and BAND.dat files by default, unless the paramerter SET_FERMI_ENERGY_ZERO in the ~/.vaspkit file is set to .FALSE.;

2. The fermi energy is read from the DOSCAR file by default, unless the FERMI_ENERGY.in file exists. You can also set the fermi energy manually in the second line of the FERMI_ENERGY.in file.