Electric field response from density-functional-perturbation theory

Density-functional-perturbation theory (DFPT) provides a way to compute the second-order linear response to an external perturbation (e.g. ionic displacement, strain, electric fields). For the case of an external electric field VASP can employ DFPT for the computation of the static dielectric tensor, the Born effective charges, and the piezoelectric tensor.

Introduction

Under the presence of a small, finite external electric field, [math]\displaystyle{ \mathcal{E_\alpha} }[/math], the changes in the total energy of the system due to [math]\displaystyle{ \mathcal{E_\alpha} }[/math] can be computed from the polarization at first order, and at second order in [math]\displaystyle{ \mathcal{E_\alpha} }[/math] from the dielectric susceptibility, Born-effective charges, and piezoelectric tensor.

All quantities can be computed via finite-differences methods (see for instance perturbation expansion after discretization or the modern theory of polarisation), but they can also be computed employing DFPT. The essential quantities are the wave functions at finite electric field, which can be approximated as

[math]\displaystyle{ |\psi^{\mathcal{E}_\alpha}_\lambda\rangle = |\psi\rangle + \lambda |\partial_{\mathcal{E}_\alpha}\psi\rangle, }[/math]

with [math]\displaystyle{ \lambda }[/math] being the vanishing strength of the field, and the static dielectric tensor

[math]\displaystyle{ \varepsilon^{\infty}(\hat{\mathbf q},0)= 1 - \frac{8 \pi e^2}{\Omega} \lim _{\mathbf{q} \rightarrow 0}\frac{1}{\mathbf q^2} \sum_{c, v, \mathbf{k}} 2 w_{\mathbf{k}} \left|\langle u_{c \mathbf{k+q}}|u_{v \mathbf{k}}\rangle\right|^2, }[/math] which requires the periodic part of the wave function to be computed at different momenta. However, in the limit of very small q, this can be expanded as

[math]\displaystyle{ u_{n \mathbf{k}+\mathbf{q}}=u_{n \mathbf{k}}+\mathbf{q} \cdot \nabla_{\mathbf{k}} u_{n \mathbf{k}}+O\left(\mathbf{q}^2\right). }[/math]

So both the derivative with respect to the crystal momentum k and the finite electric field [math]\displaystyle{ \mathcal{E_\alpha} }[/math] are required. Fortunately both can be obtained from a system of coupled equations called the Sternheimer equations.

Sternheimer equations

The starting point is density-functional theory (DFT), where one has to solve the Kohn-Sham (KS) equations

[math]\displaystyle{ H(\mathbf{k}) | \psi_{n\mathbf{k}} \rangle= e_{n\mathbf{k}}S(\mathbf{k}) | \psi_{n\mathbf{k}} \rangle, }[/math]

where [math]\displaystyle{ H(\mathbf{k}) }[/math] is the DFT Hamiltonian, [math]\displaystyle{ S(\mathbf{k}) }[/math] is the overlap operator and, [math]\displaystyle{ | \psi_{n\mathbf{k}} \rangle }[/math] and [math]\displaystyle{ e_{n\mathbf{k}} }[/math] are the KS eigenpairs.

As it was mentioned before, to compute the response with respect to an external electric field [math]\displaystyle{ \mathcal{E_\alpha} }[/math] one requires the derivative of the orbitals with respect to [math]\displaystyle{ \mathbf{k} }[/math] and [math]\displaystyle{ \mathcal{E_\alpha} }[/math]. For the former we have that

[math]\displaystyle{ \left[ H(\mathbf{k}) - e_{n\mathbf{k}}S(\mathbf{k}) \right] |\partial_\mathbf{k}\tilde{u}_{n\mathbf{k}}\rangle = -\partial_\mathbf{k} \left[ H(\mathbf{k}) - e_{n\mathbf{k}}S(\mathbf{k})\right] |\tilde{u}_{n\mathbf{k}}\rangle }[/math]

while the latter is given by

[math]\displaystyle{ \left[ H(\mathbf{k}) - e_{n\mathbf{k}}S(\mathbf{k}) \right] |\partial_\mathcal{E_\alpha}\tilde{u}_{n\mathbf{k}} \rangle = -\Delta H_{\text{SCF}}(\mathbf{k}) |\tilde{u}_{n\mathbf{k}}\rangle -\mathbf{\hat{q}_\alpha}\cdot |\vec{\beta}_{n\mathbf{k}}\rangle. }[/math]

Here [math]\displaystyle{ \Delta H_{\text{SCF}}(\mathbf{k}) |\tilde{u}_{n\mathbf{k}}\rangle }[/math] is the change in the microscopic cell periodic part of the Hamiltonian that comes from the change in wave functions (self-consistent field effects). These are computed at first order only, thus ensuring that local field effects are included directly in the evaluation of the dielectric response.

The derivative with respect to [math]\displaystyle{ \mathbf{k} }[/math] in the second equation enters via the polarization vector, [math]\displaystyle{ \vec\beta_{n\mathbf k} }[/math], which in the PAW formalism is given by

[math]\displaystyle{ |\vec{\beta}_{n\mathbf{k}}\rangle= \left( 1+\sum_{ij} |\tilde{p}_{i\mathbf{k}}\rangle Q_{ij} \langle\tilde{p}_{j\mathbf{k}}| \right) |\partial_\mathbf{k} \tilde{u}_{n\mathbf{k}}\rangle+ i\left( \sum_{ij} |\tilde{p}_{i\mathbf{k}}\rangle Q_{ij} (\mathbf{r}-\mathbf{R}_i)-\vec{\tau}_{ij} \langle\tilde{p}_{j\mathbf{k}}| \right)|\tilde{u}_{n\mathbf{k}}\rangle, }[/math]

where [math]\displaystyle{ Q_{ij} }[/math] and [math]\displaystyle{ \vec\tau_{ij} }[/math] are the norm and the dipole moments of the one-center charge densities, each respectively expressed as

[math]\displaystyle{ Q_{ij} = \int_{\Omega_\text{PAW}} \left[ \phi_i(\mathbf{r}) \phi_j(\mathbf{r}) - \tilde{\phi}_i(\mathbf{r}) \tilde{\phi}_j(\mathbf{r}) \right] d^3 \mathbf{r} }[/math]

and

[math]\displaystyle{ \vec{\tau}_{ij} = \int_{\Omega_\text{PAW}} (\mathbf{r}-\mathbf{R}_i) \left[ \phi_i(\mathbf{r}) \phi_j(\mathbf{r}) - \tilde{\phi}_i(\mathbf{r}) \tilde{\phi}_j(\mathbf{r}) \right] d^3 \mathbf{r}. }[/math]

For norm-conserving pseudopotentials the expression for [math]\displaystyle{ |\vec{\beta}_{n\mathbf{k}}\rangle }[/math] is reduced to

[math]\displaystyle{ |\vec{\beta}_{n\mathbf{k}}\rangle=|\partial_\mathbf{k} \tilde{u}_{n\mathbf{k}}\rangle. }[/math]

To employ the Sternheimer formalism instead of that of the perturbation expansion after discretisation you should set LEPSILON = TRUE  in the INCAR file. Control over how the changes in the wave functions due to the finite electric field are included in the Hamiltonian can be set via the LRPA tag:

• LRPA = TRUE  - only changes in the Hartree potential are included
• LRPA = FALSE  - both changes in the Hartree potential and in the exchange correlation potential are included

Computation of physical quantities

After obtaining the derivatives of the Kohn-Sham orbitals with respect to the electric field [math]\displaystyle{ |\partial_\mathcal{E_\alpha}\tilde{u}_{n\mathbf{k}} \rangle }[/math], the dielectric tensor can then be computed using

[math]\displaystyle{ \epsilon^\infty(\hat{\mathbf{q}},0) = 1-\frac{8\pi e^2}{\Omega} \sum_{v\mathbf{k}} 2 w_\mathbf{k} \langle \mathbf{\hat{q}}\vec{\beta}_{n\mathbf{k}} | \partial_{\mathcal{E}_\alpha} \tilde{u}_{n\mathbf{k}} \rangle. }[/math]

Note that the sum runs over occupied states only! Meaning that very few empty states are required in this calculation. The macroscopic dielectric tensor is reported in the OUTCAR in the field

MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in DFT) 

By computing the the forces on each atom, [math]\displaystyle{ a }[/math], for a given set of perturbed KS orbitals it is also possible to obtain the Born effective charges using

[math]\displaystyle{ Z^a_{ij}=\frac{\Omega}{e}\frac{\partial P_i}{\partial u^a_j} =\frac{1}{e}\frac{\partial F^a_j}{\partial \mathcal{E}_i} \approx\frac{1}{e}\frac{ \left( \mathbf{F}[\{\psi^{\mathcal{E}_i}_{\lambda/2}\}]- \mathbf{F}[\{\psi^{\mathcal{E}_i}_{-\lambda/2}\}] \right)^a_j}{\lambda}. }[/math]

These are also reported in the OUTCAR, after

BORN EFFECTIVE CHARGES (in e, cummulative output) 

Finally, the ion-clamped piezoelectric tensor [math]\displaystyle{ \overline{e}_{ij} }[/math] can be computed from the strain tensor, [math]\displaystyle{ \mathbf{\sigma}[\psi_{n\mathbf{k}}] }[/math], using

[math]\displaystyle{ \overline{e}_{ij} =-\frac{\partial \sigma_j}{\partial \mathcal{E}_i} \approx -\frac{ \left( \sigma[\{\psi^{\mathcal{E}_i}_{\lambda/2}\}]- \sigma[\{\psi^{\mathcal{E}_i}_{-\lambda/2}\}] \right)_j }{\lambda}, }[/math]

which can be found in

PIEZOELECTRIC TENSOR for field in x, y, z (e Angst) 

Ionic contributions

Piezoelectric tensor

The ionic contributions to the piezoelectric tensor can be computed by setting IBRION=6, 8, and LEPSILON=.TRUE. or LCALCEPS = TRUE  in the INCAR. This will correct the piezoelectric tensor using

[math]\displaystyle{ e_{\alpha j} = \overline{e}_{\alpha j} + \Omega_0^{-1}Z_{m\alpha}(\Phi)^{-1}_{mn}\Xi_{nj}, }[/math]

where [math]\displaystyle{ \Phi_{mn} }[/math] and [math]\displaystyle{ \Xi_{nj} }[/math] are the force constants and force response to the internal strain tensors, respectively. The second term on the right-hand side will be written to the OUTCAR in

PIEZOELECTRIC TENSOR IONIC CONTR  for field in x, y, z (C/m^2)

while toe contributions from the ion-clamped tensor are written in

PIEZOELECTRIC TENSOR (including local field effects) (e Angst)

Dielectric tensor

Setting the same tags in the INCAR as for the corrected piezoelectric tensor, it is possible to compute the ionic contributions to the dielectric tensor. These are given by

[math]\displaystyle{ \chi_{\alpha\beta} = \overline{\chi_{\alpha\beta}} + \Omega_0^{-1}Z_{m\alpha}(\Phi)^{-1}_{mn}Z_{n\beta}. }[/math]

The ionic contributions (second term on the right-hand side) are written in the OUTCAR to

MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC CONTRIBUTION

while the ion-clamped tensor is written to

MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects)

Related tags and sections

LEPSILON, LRPA, LPEAD, Born effective charges