Part 1:¶
By the end of this tutorial, you will be able to:
- name and compute static response to a small electric field, strain and ionic displacement
- set INCAR tags related to the finite differences approach for ionic displacement and dielectric response
- use the perturbation expression after discretization (PEAD) method
- distinguish ion-clamped and relaxed-ion polarizability, elastic modulus and piezoelectric tensor
1.1 Task¶
Compute the response to an applied static electric field, strain and ionic displacement for SiC using the finite differences approach.
The total energy $E$ can be expanded around the unperturbed system with total energy $E^0$ up to second order of small applied electric field $\mathcal{E}$, strain $\eta$ and ionic displacement $u$ in the following way: $$ E(\mathbf{u}, \mathbf{\mathcal{E}}, \mathbf{\eta}) = E^{0} + \frac{\partial E}{\partial u_i} u_i + \frac{\partial E}{\partial \mathcal{E}_i} \mathcal{E}_i + \frac{\partial E}{\partial \eta_{ij}} \eta_{ij} $$ $$ + \frac{1}{2} \frac{\partial^2 E}{\partial u_i\partial u_j} u_i u_j + \frac{1}{2} \frac{\partial^2 E}{\partial \mathcal{E}_i\partial \mathcal{E}_j} \mathcal{E}_i \mathcal{E}_j + \frac{1}{2} \frac{\partial^2 E}{\partial \eta_{ij}\partial \eta_{km}} \eta_{ij} \eta_{km} $$ $$\tag{1.1} + \frac{\partial^2 E}{\partial u_i\partial \mathcal{E}_j} u_i \mathcal{E}_j + \frac{\partial^2 E}{\partial u_i\partial \eta_{jk}} u_i \eta_{jk} + \frac{\partial^2 E}{\partial \mathcal{E}_i\partial \eta_{jk}} \mathcal{E}_i \eta_{jk}, $$ with an implied sum rule over repeated indices. Note that linear response to, e.g., an applied electric field also appears in second order in combination with other perturbations.
The coefficients have names and sign conventions according to their physical interpretation. In particular, the force $F$, polarization $P$, and stress $\sigma$ are defined as: $$\tag{1.2} F_i = - \Omega_0 \left. \frac{\partial E}{\partial u_i}\right|_{\mathcal{E},\eta}, \quad P_i = - \left. \frac{\partial E}{\partial \mathcal{E}_i} \right|_{u,\eta} = -\frac{fe}{(2\pi)^3} \sum_n \int_{BZ} \mathrm{d}\mathbf{k} \left\langle u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right| \mathrm{i} \nabla_{\mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle , \quad \sigma_{ij} = \left. \frac{\partial E}{\partial \eta_{ij}} \right|_{u,\mathcal{E}}, $$ where $\Omega_0$ is the volume of the unperturbed unit cell, $f=2$ is the spin multiplicity in a spin-unpolarized calculation, $\left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle$ is the cell-periodic part of a field-polarized Kohn-Sham (KS) orbital. Forces and stresses are computed using the Hellmann-Feynman theorem. Further, the Hessian or force-constant matrix $K$, frozen ion dielectric susceptibility or ion-clamped polarizability $\chi$, and the zero-field, ion-clamped elastic modulus $C$ read: $$\tag{1.3} K_{ij} = \Omega_0 \left. \frac{\partial^2 E}{\partial u_i\partial u_j}\right|_{\mathcal{E},\eta}, \quad \chi_{ij} = - \left. \frac{\partial^2 E}{\partial \mathcal{E}_i\partial \mathcal{E}_j} \right|_{u,\eta} = \left. \frac{\partial P_i}{\partial \mathcal{E}_j} \right|_{u,\eta} , \quad C_{ijkm} = \left. \frac{\partial^2 E}{\partial \eta_{ij}\partial \eta_{km}}\right|_{u,\mathcal{E}}. $$ Finally, the Born dynamical effective charge $Z^*$, the (force-response) internal-strain tensor $\Lambda$, and the ion-clamped piezoelectric tensor $e^{(0)}$ can be written as: $$\tag{1.4} Z^*_{ij} = - \frac{\Omega_0}{e} \left. \frac{\partial^2 E}{\partial u_i\partial \mathcal{E}_j} \right|_\eta = \frac{\Omega_0}{e} \left. \frac{\partial F_i}{\partial \mathcal{E}_j}\right|_\eta , \quad \Lambda_{ijk} = - \Omega_0 \left. \frac{\partial^2 E}{\partial u_i\partial \eta_{jk}} \right|_{\mathcal{E}} = - \Omega_0 \left. \frac{\partial \sigma_{jk}}{\partial u_i} \right|_{\mathcal{E}}, \quad e^{(0)}_{ijk} = - \left. \frac{\partial^2 E}{\partial \mathcal{E}_i\partial \eta_{jk}} \right|_{u} = - \left. \frac{\partial \sigma_{jk}}{\partial \mathcal{E}_i}\right|_{u} . $$
There are two independent levels on which we can take a finite differences approach:
- the ionic displacement and strain with IBRION=5 and 6,
- the applied electric field with LCALCEPS that uses the perturbation expression after discretization (PEAD) method with LPEAD
Mind: This can only be used when the system has a band gap, i.e., not for metals. The derivative of the cell-periodic part of the KS orbital with respect to $\mathbf{k}$, $\nabla_{\mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle$, is computed using a discretized version of the Sternheimer equation. There is no sum over empty bands.
1.2 Input¶
The input files to run this example are prepared at $TUTORIALS/response/e01_SiC-static-finite-differences. Check them out!
fcc SiC
4.35
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
Si C
1 1
Cartesian
0.00 0.00 0.00 Si
0.25 0.25 0.25 C
SYSTEM = fcc SiC
ISMEAR = 0 ! Gaussian smearing
SIGMA = 0.01 ! smearing in eV
EDIFF = 1E-6 ! convergence criterion
ENCUT = 400 ! energy cutoff
electric/INCAR
SYSTEM = fcc SiC
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1E-6
ENCUT = 400
! finite differences for polarization
LCALCPOL = T
LCALCEPS = T
LPEAD = T
EFIELD_PEAD = 0.1 0.1 0.1
electric-ionic/INCAR
SYSTEM = fcc SiC
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1E-6
ENCUT = 400
! finite differences for polarization
LCALCPOL = T
LCALCEPS = T
LPEAD = T
EFIELD_PEAD = 0.01 0.01 0.01
IPEAD = 4
! finite differences for ionic displacements
IBRION = 6
POTIM = 0.015
NFREE = 2
ISIF = 3
K-Points
0
Gamma
4 4 4
0 0 0
k points for band structure
10 ! intersections
line mode
reciprocal
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.0000000000 0.5000000000 X
0.5000000000 0.0000000000 0.5000000000 X
0.6250000000 0.2500000000 0.6250000000 U
0.3750000000 0.3750000000 0.7500000000 K
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.5000000000 0.5000000000 L
0.5000000000 0.5000000000 0.5000000000 L
0.5000000000 0.2500000000 0.7500000000 W
Pseudopotentials of Si and C.
The POSCAR file defines the face-centered-cubic structure of SiC.
The first INCAR file specifies a density-functional-theory (DFT) run using the default generalized-gradient-approximation functional of Perdew, Burke, and Ernzerhof (GGA-PBE) with Gaussian smearing set by ISMEAR and SIGMA. The energy cutoff is set relatively high to avoid Pulay stress. The electronic relaxation will stop when the change in band energies and total energy is each below the value set by the EDIFF tag.
The electric/INCAR file, defines calculations at small, static electric fields $\mathcal{E}_i$ on top of DFT with LCALCPOL=T to compute the electronic polarization. The electronic polarization of a periodic system is $-fe/(2\pi)^3$ times the sum over the Berry phases $\int_{BZ} \mathrm{d}\mathbf{k} \left\langle u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right| \mathrm{i} \nabla_{\mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle$ of all occupied states, see Equation (1.2).
With LCALCEPS=T, we use the perturbation expression after discretization (PEAD) method, which yields the following finite differences expression for the derivative of the cell-periodic part of the field-polarized KS orbitals with respect to $\mathbf{k}$:
$$\tag{1.5} \frac{\partial \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}_j} \right\rangle}{\partial k} = \frac{1}{2\Delta k} \sum_{m=1}^N \left[ \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}_{j+1}} \right\rangle S_{mn}^{-1}(\mathbf{k}_j, \mathbf{k}_{j+1}) - \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}_{j-1}} \right\rangle S_{mn}^{-1}(\mathbf{k}_j, \mathbf{k}_{j-1}) \right], $$
where $S_{mn}(\mathbf{k}_j, \mathbf{k}_{j+1}) = \left\langle u^{(\mathcal{E}_i)}_{n\mathbf{k}_{j}} \right. \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}_{j+1}} \right\rangle$ is the overlap operator. What are the default values for the tags related to LPEAD? What do you need to consider, to set these tags appropriately for your system?
Click to see the answer!
If LPEAD=T, then
EFIELD_PEAD = 0.01 0.01 0.01
IPEAD = 4
Here, EFIELD_PEAD is the step size and IPEAD is the number of steps for the finite difference operator. You must consider the size of the band gap and the number of $\mathbf{k}$ points when setting EFIELD_PEAD. Read the VASP Wiki for more details!
The electric-ionic/INCAR file, additionally switches on ionic displacements using finite differences with IBRION. The step size and number of steps is set by POTIM and NFREE, respectively. ISIF=3 allows for changes in volume and cell shape. This is necessary to compute the elastic modulus.
The KPOINTS file defines an equally spaced Γ-centered grid of $\mathbf{k}$ points, while the KPOINTS_OPT file contains high-symmetry k points for the band structure.
The POTCAR file contains the pseudopotential for Si and C.
1.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/response/e01_*
vasp_std
This calculation computes the DFT ground state. Plot the band structure with the $\mathbf{k}$ points specified in the KPOINTS_OPT file using py4vasp! What is the size of the band gap?
import py4vasp
import numpy as np
my_calc = py4vasp.Calculation.from_path("./e01_SiC-static-finite-differences")
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
print("Band gap of SiC: ", gap(my_calc), "eV")
my_calc.band.plot("kpoints_opt")
# Optionally, you can adjust the height of the band-structure plot as follows
# my_calc.band.plot("kpoints_opt").to_plotly().update_layout(height=1000)
Next, let us perturb the system with a small, static electric field. Based on the previous calculation, run VASP in the subdirectory electric!
cd $TUTORIALS/response/e01_*/electric
cp ../WAVECAR .
mpirun -np 2 vasp_std
What warning does this calculation produce? Why and what can you do to avoid it?
Click to see the answer!
The electric field is too strong such that Zener tunneling can occur. You can reduce the number of $\mathbf{k}$ points and/or reduce the electric field set by EFIELD_PEAD. In the next calculation EFIELD_PEAD is set one order of magnitude smaller.
In the next calculation, we additionally consider ionic displacements and changes to the cell volume and cell shape. This gives access to all the static response properties defined in Equation (2) - (4). Start from the WAVECAR of the first DFT calcultion and run VASP in the subdirectory electric-ionic!
cd $TUTORIALS/response/e01_*/electric-ionic
cp ../WAVECAR .
mpirun -np 2 vasp_std
You can print the response properties using py4vasp in the following way:
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e01_SiC-static-finite-differences/electric-ionic")
my_calc.force.print()
print(" ")
my_calc.polarization.print()
print(" ")
my_calc.stress.print()
# If you just executed the code box above, you can omit the following two lines
#import py4vasp
#my_calc = py4vasp.Calculation.from_path("./e01_SiC-static-finite-differences/electric-ionic")
my_calc.force_constant.print()
print(" ")
my_calc.elastic_modulus.print()
For the comparison with experimental data, the ion-clamped polarizability, elastic modulus and piezoelectric tensor defined in Equation(1.3) is of limited interest. As you can see in the output of above code box, an relaxed-ion version of these quantities is computed. Those are defined as
$$\tag{1.6} \left. \chi_{ij} \right|_{\eta} = \left. \chi_{ij} \right|_{u,\eta} + \Omega_0^{-1} Z_{mi} (K^{-1})_{mn} Z_{nj}, $$
$$\tag{1.7} \left. C_{ijkm} \right|_{\mathcal{E}} = \left. C_{ijkm} \right|_{u,\mathcal{E}} - \Omega_0^{-1} \Lambda_{nij} (K^{-1})_{no} \Lambda_{okm} $$
$$\tag{1.8} e^{(0)}_{ijk} = \left. e^{(0)}_{ijk} \right|_{u} - \Omega_0^{-1} Z_{mi} (K^{-1})_{mn} \Lambda_{njk} $$
with an implied sum over repeated indices. Note that we did not explicitly print the polarizability here because VASP directly uses it to construct the dielectric tensor. This is typically more convenient in practice in will be discussed in more detail in the next section.
# If you just executed the code box above, you can omit the following two lines
#import py4vasp
#my_calc = py4vasp.Calculation.from_path("./e01_SiC-static-finite-differences/electric-ionic")
my_calc.born_effective_charge.print()
print(" ")
my_calc.internal_strain.print()
print(" ")
my_calc.piezoelectric_tensor.print()
You can find the same information towards the end of the OUTCAR file. Search for the corresponding lines!
By the end of this tutorial, you will be able to:
- relate the dielectric constant to the polarizability
- use density-functional-perturbation theory (DFPT) to compute static dielectric properties
- explain the long-wave-length approximation and the difference between macroscopic and microscopic response
- state the advantages and disadvantages of LEPSILON vs LCALCEPS
2.1 Task¶
Compute the static dielectric constant of AlP within the long-wave-length approximation neglecting local-field effects using density-functional-perturbation theory.
The dielectric constant couples the electric field in a material $\mathcal{E}_i(\mathbf{r}, \omega)$ to an applied, external electric field $\mathcal{E}_{\mathrm{ext}, i}(\mathbf{r}', \omega)$. This yields
$$\tag{2.1} \mathcal{E}_i(\mathbf{r}, \omega) = \int \mathrm{d}\mathbf{r}'\epsilon^{-1}_{ij}(\mathbf{r}, \mathbf{r}', \omega) \mathcal{E}_{\mathrm{ext}, j}(\mathbf{r}',\omega) $$
and in momentum space we obtain
$$\tag{2.2} \mathcal{E}_i(\mathbf{q}+\mathbf{G}, \omega) = \sum_{\mathbf{G}'} \epsilon^{-1}_{\mathbf{G}\mathbf{G'}ij}(\mathbf{q}, \omega) \mathcal{E}_{\mathrm{ext}, j}(\mathbf{q}+\mathbf{G}', \omega). $$
We can neglect the details of the local environment at $\mathbf{r}$ and $\mathbf{r}'$, if the material is homogeneous and/or the electric field varies on much larger length scales than atomic distances. This is the so-called long-wave-length approximation, where $\epsilon_{ij}(\mathbf{r}-\mathbf{r}', \omega)$ depends only on the distance $\mathbf{r}-\mathbf{r}'$, and Equation (2.2) takes on the following simple form:
$$\tag{2.3} \mathcal{E}_i(\mathbf{q}, \omega) = \epsilon^{-1}_{00\,ij}(\mathbf{q}, \omega) \mathcal{E}_{\mathrm{ext}, j}(\mathbf{q}, \omega). $$
In other words, the off-diagonal elements for $\mathbf{G}\neq \mathbf{G'}$ vanish, if local-field effects can be neglected.
The ion-clamped, electronic dielectric constant is closely related to the ion-clamped polarizability through the following relation:
$$\tag{2.4} \epsilon^{\infty}_{ij} = \delta_{ij} + \frac{4\pi}{\epsilon_0} \chi_{ij}, $$
where $\delta_{ij}$ is the Kronecker delta, here $\chi_{ij}$ is defined in Equation (1.3), and $\epsilon_0$ is the vacuum permittivity. $\epsilon^{\infty}$ refers to the response to alternating current (AC), or optical fields at frequencies well above the phonon-frequency range. This corresponds to the limit $\mathbf{q}\to 0$ and is indicated by the $\infty$ label, which is a common notation in the polariton community.
2.2 Input¶
The input files to run this example are prepared at $TUTORIALS/response/e02_AlP-static-DFPT. Check them out!
Al P F-43m
1.0
3.372404 0.000000 1.947058
1.124135 3.179533 1.947058
0.000000 0.000000 3.894116
Al P
1 1
direct
0.000000 0.000000 0.000000 Al
0.250000 0.250000 0.250000 P
SYSTEM = AlP
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1E-6
ENCUT = 400
electric-ionic/INCAR
SYSTEM = AlP
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1E-6
ENCUT = 400
LASPH = T
! finite differences for electric field
LPEAD = T
! dfpt for polarization
LEPSILON = T
! dfpt for ionic displacements
IBRION = 8
K-Points
0
Gamma
4 4 4
0 0 0
Pseudopotential of Al and P.
The POSCAR file defines the structure of AlP.
The first INCAR file specifies a DFT run using the default GGA-PBE. Make comments in the INCAR file to recall the meaning of each tag!
The electric-ionic/INCAR file, defines a calculation using DFPT for both, the dielectric response with LEPSILON=T and the ionic displacement and strain with IBRION=8.
Within density-functional-perturbation theory (DFPT), the linear response of the cell-periodic part of the KS orbitals to an electric field, $\left| \xi_{n\mathbf{k}} \right\rangle$, can be found in terms of the unperturbed system by solving the following linear Sternheimer equation:
$$\tag{2.5} \left[ H(\mathbf{k}) - \epsilon_{n\mathbf{k}} S(\mathbf{k}) \right] \left| \xi_{n\mathbf{k}} \right\rangle = - \Delta H_{\mathrm{SCF}}(\mathbf{k}) \left| u_{n\mathbf{k}} \right\rangle - \hat{\mathbf{q}} \left| \nabla_{\mathbf{k}} u_{n\mathbf{k}} \right\rangle, $$
where $\Delta H_{\mathrm{SCF}}(\mathbf{k})$ is the local-field effect, i.e., the microscopic cell-periodic change in the KS Hamiltonian due to changes in the KS orbitals. When LRPA=F, then $\Delta H_{\mathrm{SCF}}(\mathbf{k})=0$.
Next, still within DFPT, the derivative of the cell-periodic part of the KS orbitals with respect to $\mathbf{k}$ is determined using a second linear Sternheimer equation, that reads
$$\tag{2.6} \left[ H(\mathbf{k}) - \epsilon_{n\mathbf{k}} S(\mathbf{k}) \right] \left| \nabla_{\mathbf{k}} u_{n\mathbf{k}} \right\rangle = \frac{\partial \left[ H(\mathbf{k}) - \epsilon_{n\mathbf{k}} S(\mathbf{k}) \right]}{\partial \mathbf{k}} \left| u_{n\mathbf{k}} \right\rangle. $$
Finally, the static, macroscopic dielectric matrix is given by
$$\tag{2.7} \hat{\mathbf{q}} \cdot \epsilon^{\infty} \cdot \hat{\mathbf{q}} = 1 - \frac{f 8\pi e^2}{\Omega_0} \sum_{n\mathbf{k}} w_\mathbf{k} \left\langle \hat{\mathbf{q}} \nabla_{\mathbf{k}} u_{n\mathbf{k}} \right. \left| \xi_{n\mathbf{k}} \right\rangle. $$
Note that, here we performed the same kind of expansion as in Equation (1.1), but instead of a finite differences approach with LCALCEPS, we compute the response within DFPT with LEPSILON. In contrast to LCALCEPS, LEPSILON is applicable to metallic systems. In both cases, there is no sum over unoccupied states. Currently LEPSILON cannot be used for exchange-correlation functionals that depend explicitly on the orbitals, e.g., HF-like and hybrid functionals.
The POTCAR file contains the pseudopotential for Al and P.
2.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/response/e02_*
vasp_std
Then, compute the dielectric response in the subdirectory electric-ionic:
cd $TUTORIALS/response/e02_*/electric-ionic
cp ../WAVECAR .
mpirun -np 4 vasp_std
In the long-wave-length approximation, where local-field effects can be neglected, the macroscopic, AC dielectric response $\epsilon^{\infty}$ is given by
$$\tag{2.5} \epsilon^{\infty}_{ij}(\hat{\mathbf{q}},\omega) \approx \lim_{\mathbf{q}\to 0} \epsilon_{00\,ij}(\mathbf{q}, \omega). $$
On the other hand, if local-field effects cannot be neglected, the inverse macroscopic, AC dielectric response is related to the $\mathbf{q}=0$ limit of the microscopic, electronic dielectric constant as follows:
$$\tag{2.6} \frac{1}{\epsilon^{\infty}_{ij}(\hat{\mathbf{q}},\omega)} \equiv \lim_{\mathbf{q}\to 0} \epsilon^{-1}_{00\,ij}(\mathbf{q}, \omega). $$
Print the static dielectric response using py4vasp!
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e02_AlP-static-DFPT/electric-ionic")
my_calc.dielectric_tensor.print()
print(" ")
my_calc.born_effective_charge.print()
print(" ")
my_calc.piezoelectric_tensor.print()
2.4 Questions¶
- What is the AC dielectric constant? How is it connected to the polarizability?
- What pros and cons do the two methods, density-functional-perturbation theory and the finite differences approach, have when computing the dielectric constant? Which INCAR tags control these methods?
- When can local-field effects be neglected?
- What is static response?
By the end of this tutorial, you will be able to:
- compute the frequency-dependent dielectric response
- relate the imaginary and real part of the dielectric function using Kramers-Kronig relations
- relate the dielectric function to optical properties
- consider ionic contributions to the frequency-dependent dielectric response
- compute phonon frequencies
- convert units of frequencies
3.1 Task¶
Compute the frequency-dependent dielectric function of NaCl in the independent particle approximation.
The full frequency-dependent dielectric function comprises electronic and ionic contributions. These have to be computed separately on top of a DFT run, and then added together:
$$\tag{3.1} \varepsilon (\omega )=\varepsilon _{{{\rm {elec}}}}(\omega )+\varepsilon _{{{\rm {ion}}}}(\omega ). $$
The imaginary part and the real part of the dielectric function are related by the Kramers-Kronig relation. Starting from the frequency-dependent dielectric function, other optical properties such as the reflectivity, absorption and the optical conductivity can be computed.
For more information read the documentation of the LOPTICS tag on the VASP Wiki!
3.2 Input¶
The input files to run this example are prepared at $TUTORIALS/response/e03_NaCl-optics. Check them out!
NaCl fcc
5.55596202
0.500 0.500 0.000
0.000 0.500 0.500
0.500 0.000 0.500
Na Cl
1 1
Direct
0.000 0.000 0.000 Na
0.500 0.500 0.500 Cl
SYSTEM = NaCl
ALGO = Normal
ISMEAR = -5
EDIFF = 1.E-6
LORBIT = T ! compute DOS
electron/INCAR
SYSTEM = NaCl
ALGO = Exact
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1.E-6
LASPH = T
! electronic dielectric function with Kramers-Kronig
LOPTICS = T ! frequency-dependent dielectric matrix
CSHIFT = 0.1 ! complex shift (default)
NBANDS = 32 ! sum over unoccupied bands
NEDOS = 2000 ! frequency grid
ion-dfpt/INCAR
SYSTEM = NaCl
ALGO = Normal
PREC = High
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1.E-8
! ionic dielectric function with Kramers-Kronig
! option 1: using purely DFPT
IBRION = 8
LEPSILON = T
ion-finite-differences/INCAR
SYSTEM = NaCl
ALGO = Normal
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1.E-6
! ionic dielectric function with Kramers-Kronig
! option 2: using purely a finite differences approach
IBRION = 6
NFREE = 2
POTIM = 0.015
TIME = 0.1
LPEAD = T
LCALCEPS = T
K-Points
0
Gamma centered
4 4 4
0 0 0
GW pseudopotential of Na, and Cl.
Check the meaning of all tags used in the INCAR file and make comments in the file!
The electronic polarization of a periodic system is $-fe/(2\pi)^3$ times the sum over the Berry phases $\int_{BZ} \mathrm{d}\mathbf{k} \left\langle u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right| \mathrm{i} \nabla_{\mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle$ of all occupied states, see Equation (1.2). The derivative of the cell-periodic part of the field-polarized KS orbitals with respect to $\mathbf{k}$, $\nabla_{\mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle$, can be computed within second-order perturbation theory:
$$\tag{3.2} \nabla_{\mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle = \sum_{n'\neq n} \frac{ \left| u^{(\mathcal{E}_i)}_{n'\mathbf{k}} \right\rangle \left\langle u^{(\mathcal{E}_i)}_{n'\mathbf{k}} \right| \frac{\partial H(\mathbf{k}) - \epsilon_{n\mathbf{k}} S(\mathbf{k})}{\partial \mathbf{k}} \left| u^{(\mathcal{E}_i)}_{n\mathbf{k}} \right\rangle}{\epsilon_{n\mathbf{k}}-\epsilon_{n'\mathbf{k}}} $$
Here, $H(\mathbf{k})$ is the KS Hamiltonian of the cell-periodic part of the KS orbitals with eigenenergies $\epsilon_{n\mathbf{k}}$, and $S(\mathbf{k})$ is the overlap operator. Equation (3.2) is used when LOPTICS=T. This introduces a sum over unoccupied states! Therefore, NBANDS must be set such that sufficient unoccupied bands are included. In other words, the polarization should be converged with respect to increasing NBANDS.
The POTCAR file contains the pseudopotential for Na and Cl, that are recommended for GW calculations. Why are we not using the pseudopotential recommended for DFT calculations?
Click to see the answer!
In standard DFT calculations, unoccupied bands do not contribute to any ground-state property, e.g., the total energy. Pseudopotentials that are constructed for GW calculations are also tested to describe properties above the Fermi level. So when you increase NBANDS, you should generally opt for pseudopotential, that are recommended for GW calculations.
3.3 Calculation¶
Open a terminal, navigate to this example's directory and run the DFT calculation with VASP by entering the following:
cd $TUTORIALS/response/e03_*
vasp_std
You can print the band gap and plot the density of states (DOS) using py4vasp!
import py4vasp
import numpy as np
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
my_calc = py4vasp.Calculation.from_path("./e03_NaCl-optics")
print("Band gap of NaCl:", gap(my_calc), "eV")
my_calc.dos.plot()
Next, compute the electronic contribution to the frequency-dependent dielectric function. Copy the WAVECAR file to the subdirectory electron and run VASP there!
cd $TUTORIALS/response/e03_*/electron
cp ../WAVECAR .
mpirun -np 4 vasp_std
Finally, extract the dielectric function from the vasprun.xml file using the following script:
get_dielectric_function.sh
#!/bin/bash
file="vasprun.xml"
outfile="dielectric_function.dat"
echo "Extracting dielectric function from " $file
awk '
/<dielectricfunction comment="density-density"/ { on=1 }
on==1 && /<r>/ { print $2,($3+$4+$5)/3 }
on==1 && /<real>/ { print " " }
on==1 && /<imag>/ { print " " }
/<\/dielectricfunction/ { on=0 }
' < $file > $outfile
echo "Data written to " $outfile
Run the script and plot the result!
cd $TUTORIALS/response/e03_*/electron
bash get_dielectric_function.sh
import numpy as np
from py4vasp import plot
filename = "./e03_NaCl-optics/electron/dielectric_function.dat"
nedos = 2000
# read the dielectric function from file
w, eps = np.loadtxt(filename, unpack=True)
w = w[0:nedos]
eps_imag = eps[0:nedos]
eps_real = eps[nedos:2*nedos]
# plot the dielectric function as a function of frequency
plot(
(w, eps_real, "Dielectric function (real)"),
(w, eps_imag, "Dielectric function (imaginary)"),
xlabel = "Frequency (eV)",
)
Alternatively, you can plot the dielectric function using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path("./e03_NaCl-optics/electron")
mycalc.dielectric_function.plot()
For the ionic contribution, you can choose to use DFPT or take a finite difference approach. These are prepared in the subdirectories ion-dfpt and ion-finite-differences, respectively. Run both calculations and compare the results!
cd $TUTORIALS/response/e03_*/ion-dfpt
cp ../WAVECAR .
mpirun -np 4 vasp_std
bash get_dielectric_function.sh
import numpy as np
from py4vasp import plot
filename1 = "./e03_NaCl-optics/ion-dfpt/dielectric_function.dat"
# read the dielectric function from file
w, eps = np.loadtxt(filename1, unpack=True)
nw = int(len(w)/2)
w = w[0:nw]
eps_imag = eps[0:nw]
eps_real = eps[nw:2*nw]
# plot the dielectric function as a function of frequency
# DFPT
plot(
(w, eps_real, "Dielectric function (real)"),
(w, eps_imag, "Dielectric function (imaginary)"),
xlabel = "Frequency (2π THz)",
)
The finite differences-approach will take much longer.
cd $TUTORIALS/response/e03_*/ion-finite-differences
cp ../WAVECAR .
mpirun -np 4 vasp_std
bash get_dielectric_function.sh
import numpy as np
from py4vasp import plot
filename2 = "./e03_NaCl-optics/ion-finite-differences/dielectric_function.dat"
w, eps = np.loadtxt(filename2, unpack=True)
nw = int(len(w)/2)
w = w[0:nw]
eps_imag = eps[0:nw]
eps_real = eps[nw:2*nw]
# plot the dielectric function as a function of frequency
# finite differences
plot(
(w, eps_real, "Dielectric function (real)"),
(w, eps_imag, "Dielectric function (imaginary)"),
xlabel = "Frequency (2π THz)",
)
Do the two approaches, DFPT and finite differences, yield the identical result? Why?
With the ionic displacement, we obtain vibration equations of the form: $$\tag{3.2} m_{\mathrm{ion}} \omega^2 \mathbf{u} = \left( \frac{\partial^2 E}{\partial u_i \partial u_j} \right) \mathbf{u}, $$ where $\omega$ are the phonon frequencies, $m_{\mathrm{ion}}$ is the mass of the ion, $\mathbf{u}$ is the ionic displacement, and the Hessian matrix is the same as in Equation (1.1) and (1.3). Open the OUTCAR file and find the phonon frequencies below:
Eigenvectors and eigenvalues of the dynamical matrix
----------------------------------------------------
How are the phonon frequencies related to the ionic contribution to the dielectric function?
Click to see the answer!
If a phonon mode is dipole-active, it will appear in the frequency-dependent dielectric function. A phonon mode is dipole-active, if the electric dipole moment of all ions involved in the mode changes in size during one period. This is the case here as we can see in the eigenvectors of the dynamical matrix.
Besides the three dipole-active modes, there are three imaginary (f/i) phonon modes. Here these modes do not indicate that our structure is dynamically unstable. The finite plane-wave basis set that we have used can give slightly different answers for the total energy depending on where in the cell we place the center of mass of all ions combined.
For each phonon mode, the ionic position and the mass-weighted, ionic displacement is listed. Are the displacement vectors for the two ions, Na and Cl, in, e.g., mode 4 parallel? How do their lengths compare?
Click to see the answer!
The two vectors are parallel, but do not have the same length. That is the result of mass-weighting of the eigenvectors. This is a consequence of the ionic mass occurring in Equation (3.2). You can get the masses from the POMASS tag in the POTCAR file. Then, calculate the mass ratio $$ \sqrt{ \frac{ m_{\mathrm{Na}} }{ m_{\mathrm{Cl}} } }. $$ The mass ratio is the same as the ratio between the two displacement vectors.
The phonon frequencies are given in various different units. And in fact the frequency axes of the ionic and the electronic dielectric function, that we calculated in this example, are not given in the same units here. Compare the plots! How can the units be converted? How do the energy scales of the ionic and the electronic contributions compare?
Click to see a hint!
From Equation (3.2) we see that in our calculation the unit of frequency is $$ [\omega] = \sqrt{ \frac{1}{ [ m_{\mathrm{ion}}] } \left[ \frac{\partial^2 E}{\partial u_i \partial u_j} \right] }. $$ POMASS is given in atomic mass units (a.m.u.): $$ [ m_{\mathrm{ion}}] = 1.660599\cdot 10^{-27} \mathrm{kg} $$ energy is given in eV and the ionic displacement is given in Å. So that the unit of frequency can be written in SI as $$ 1 [\omega] = \sqrt{ \frac{\mathrm{eV}/Å^2}{ \mathrm{a.m.u.} } } = 9.822517\cdot 10^{13} \mathrm{s}^{-1}. $$ How can this be converted to other typical frequency units, such as eV, THz, and cm$^{-1}$?
Click to see the answer!
In order to convert the unit of frequency in our calculation to SI units, we did the following calculation:
$$ \sqrt{ \frac{\mathrm{eV}/Å^2}{ \mathrm{a.m.u.} } } = \sqrt{ \frac{1.602176487\cdot10^{-19}\mathrm{J}/(10^{-20} \mathrm{m}^2)}{ 1.660599\cdot 10^{-27} \mathrm{kg} } } = \sqrt{ \frac{16.02176487\mathrm{J}/\mathrm{m}^2}{ 1.660599\cdot 10^{-27} \mathrm{kg} } } = 9.822517\cdot 10^{13} \mathrm{s}^{-1}. $$
Let us do a similar manipulation for the other typical frequency units, i.e., eV, THz, and cm$^{-1}$. To get eV we use $E_{\mathrm{phonon}}=\hbar \omega$ and yield
$$ x [E_{\mathrm{phonon}}] = 1 [\omega] \cdot \hbar = 9.822517\cdot10^{13} s^{−1} \times 1.054572\cdot10^{−34} \mathrm{Js} = 1.035855\cdot 10^{−20} \mathrm{J} = 0.064652976\, \mathrm{eV} $$
Then, the angular frequency converted using $\nu = \omega/2\pi$
$$ x [\nu] = 1 [\omega]/2\pi = 98.22517/ (2\pi) \cdot 10^{12} \mathrm{Hz} = 15.6330214 \,\mathrm{THz}. $$
Finally, the inverse wavelength $\lambda^{-1} = \nu/c$, where $c$ is the speed of light. That yields
$$ x [\lambda^{-1}] = 1 [\nu]/c = \frac{ 15.6330214 \cdot 10^{12} \mathrm{s}^{-1}}{ 29979245800\, \mathrm{cm/s}} = 521.46146389 \,\mathrm{cm}^{-1}. $$
In summary, the frequency units can be converted using the following conversion:
$$ 1 \, \mathrm{THz} = 4.1357 \, \mathrm{meV} = 33.356 \, \mathrm{cm}^{-1} $$
$$ 1 \, \mathrm{meV} = 0.242 \, \mathrm{THz} = 8.066 \, \mathrm{cm}^{-1} $$
$$ 1 \, \mathrm{cm}^{-1} = 0.030 \, \mathrm{THz} = 0.124 \, \mathrm{meV}. $$
Plot the ionic contribution to the total dielectric function, which was computed using DFPT, as a function of frequencies in meV!
import numpy as np
from py4vasp import plot
filename = "./e03_NaCl-optics/ion-dfpt/dielectric_function.dat"
# read the dielectric function from file
w, eps = np.loadtxt(filename, unpack=True)
nw = int(len(w)/2)
w = w[0:nw]
eps_imag = eps[0:nw]
eps_real = eps[nw:2*nw]
# unit conversion from 2piTHz to meV
w = w*4.135667 / (2*np.pi)
# plot the dielectric function as a function of frequency
plot(
(w, eps_real, "Dielectric function (real)"),
(w, eps_imag, "Dielectric function (imaginary)"),
xlabel = "Frequency (meV)",
)
Note that the maximum frequency at which the ionic, dielectric function is computed is determined automatically based on the phonon frequencies. Specifically, it is $1.2\times \omega_{max}$.
Alternatively, you can use py4vasp!
import py4vasp
dfpt_calc = py4vasp.Calculation.from_path("./e03_NaCl-optics/ion-dfpt")
dfpt_calc.dielectric_function.plot("ion")
3.4 Questions¶
- How do phonons enter the frequency-dependent dielectric function?
- What is computed when LOPTICS=T?
- How can you set the number of frequencies at which the electronic contribution to the dielectric function is computed?
- How can the unit of frequency be converted from $\mathrm{cm}^{-1}$ to $\mathrm{eV}$?