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]

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

[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]

The perturbed orbitals, [math]\displaystyle{ \psi^{\mathcal{E}_\alpha}_\lambda }[/math], under a small, finite electric field of magnitude [math]\displaystyle{ \lambda }[/math] are given by

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

and can be used to compute 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]

where [math]\displaystyle{ \mathbf{F}[\psi] }[/math] are the forces on each atom, [math]\displaystyle{ a }[/math], for a given set of perturbed KS orbitals.

Finally, the piezoelectric tensor [math]\displaystyle{ \overline{e}_{ij} }[/math] can be computed 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]

with [math]\displaystyle{ \mathbf{\sigma}[\psi_{n\mathbf{k}}] }[/math] being the strain tensor.

Related tags and sections

LEPSILON, LRPA, LPEAD, Born effective charges