Electric field response from density-functional-perturbation theory

Density-functional-perturbation theory provides a way to compute the second-order linear response to ionic displacement, strain, and electric fields. The equations are derived as follows.

In density-functional theory, we 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 eigenstates.

To compute the response with respect to an external electric field one has to solve two Sternheimer equations. The first one to obtain the derivative of the orbitals with respect to [math]\displaystyle{ \mathbf{k} }[/math]

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

and a second one with respect to an external electric field

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

with

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

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

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

from this last step 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]

Knowing the derivative of the Kohn-Sham orbitals with respect to an electric field one can compute the orbital with a finite electric field

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

and then compute born effective charges and piezoelectric tensor in the same way as what is done in when applying finite electric fields.