Berry phases and finite electric fields

Modern Theory of Polarization

Berry phase expression for the macroscopic polarization

Calculating the change in dipole moment per unit cell under PBC's, is a nontrivial task. In general one cannot define it as the first moment of the induced change in charge density δ(r), through

$\Delta \mathbf {P} ={\frac {1}{\Omega _{0}}}\int _{\Omega _{0}}\mathbf {r} \delta \left(\mathbf {r} \right)d^{3}r$

without introducing a dependency on the shape of Ω₀, the chosen unit cell.

Recently King-Smith and Vanderbilt, building on the work of Resta, showed that the electronic contribution to the difference in polarization ΔP_e, due to a finite adiabatic change in the Hamiltonian of a system, can be identified as a geometric quantum phase or Berry phase of the valence wave functions. We will briefly summarize the essential results (for a review of geometric quantum phases in polarization theory see the papers of Resta).

Central to the modern theory of polarization is the proposition of Resta to write the electronic contribution to the change in polarization due to a finite adiabatic change in the Kohn-Sham Hamiltonian of the crystalline solid, as

$\Delta \mathbf {P} _{e}=\int _{\lambda _{1}}^{\lambda _{2}}{\partial \mathbf {P} _{e} \over \partial \lambda }d\lambda$

with

${\partial \mathbf {P} _{e} \over \partial \lambda }={i|e|\hbar \over N\Omega _{0}m_{e}}\sum _{\mathbf {k} }\sum _{n=1}^{M}\sum _{m=M+1}^{\infty }{\langle \psi _{n\mathbf {k} }^{\left(\lambda \right)}|\mathbf {\hat {p}} |\psi _{m\mathbf {k} }^{\left(\lambda \right)}\rangle \langle \psi _{m\mathbf {k} }^{\left(\lambda \right)}|\partial V^{\left(\lambda \right)}/\partial \lambda |\psi _{n\mathbf {k} }^{\left(\lambda \right)}\rangle \over \left(\epsilon _{n\mathbf {k} }^{\left(\lambda \right)}-\epsilon _{m\mathbf {k} }^{\left(\lambda \right)}\right)^{2}}+\mathrm {c.c.}$

where m_e and e are the electronic mass and charge, N is the number of unit cells in the crystal, Ω₀ is the unit cell volume, M is the number of occupied bands, p is the momentum operator, and the functions ψ^(λ)_nk are the usual Bloch solutions to the crystalline Hamiltonian. Within Kohn-Sham density-functional theory, the potential V^(λ) is to be interpreted as the Kohn-Sham potential V^(λ)_KS, where λ parameterizes some change in this potential, for instance due to the displacement of an atom in the unit cell.

King-Smith and Vanderbilt have cast this expression in a form in which the conduction band states ψ^(λ)_mk no longer explicitly appear, and they show that the change in polarization along an arbitrary path, can be found from only a knowledge of the system at the end points

$\Delta \mathbf {P} _{e}=\mathbf {P} _{e}^{\left(\lambda _{2}\right)}-\mathbf {P} _{e}^{\left(\lambda _{1}\right)}$

with

$\mathbf {P} _{e}^{\left(\lambda \right)}=-{if|e| \over 8\pi ^{3}}\sum _{n=1}^{M}\int _{BZ}d^{3}k\langle u_{n\mathbf {k} }^{\left(\lambda \right)}|\nabla _{\mathbf {k} }|u_{n\mathbf {k} }^{\left(\lambda \right)}\rangle$

where f is the occupation number of the states in the valence bands, u^(λ)_nk is the cell-periodic part of the Bloch function ψ^(λ)_nk, and the sum n runs over all M occupied bands.

The physics behind the equation above becomes more transparent when this expression is written in terms of the Wannier functions of the occupied bands,

$\mathbf {P} _{e}^{\left(\lambda \right)}=-{f|e| \over \Omega _{0}}\sum _{n=1}^{M}\langle W_{n}^{\left(\lambda \right)}|\mathbf {r} |W_{n}^{\left(\lambda \right)}\rangle$

where W_n is the Wannier function corresponding to valence band n.

This shows the change in polarization of a solid, induced by an adiabatic change in the Hamiltonian, to be proportional to the displacement of the charge centers r_n=⟨ W^(λ)_n|r| W^(λ)_n⟩, of the Wannier functions corresponding to the valence bands.

It is important to realize that the polarization in terms of Bloch or Wannier functions, and consequently the change in polarization, is only well-defined modulo feR/Ω₀, where R is a lattice vector. This indeterminacy stems from the fact that the charge center of a Wannier function is only invariant modulo R, with respect to the choice of phase of the Bloch functions.

In practice one is usually interested in polarization changes |ΔP_e| << |feR₁/Ω₀|, where R₁ is the shortest nonzero lattice vector. An arbitrary term feR/Ω₀ can therefore often be removed by simple inspection of the results. In cases where |ΔP_e| is of the same order of magnitude as feR₁/Ω₀ any uncertainty can always be removed by dividing the total change in the Hamiltonian λ₁→λ₂ into a number of intervals.

Computational aspects

In general, the direct evaluation of P_e^(λ) is useless, because there is no specific relationship between the phases of the eigenvectors u^(λ)_kn generated by a numerical diagonalization routine. This problem is circumvented by dividing the Brillouin zone integration in two parts, a two-dimensional integral and a line integral, and by transforming the above into three equations, which separately provide the components of P_e^(λ) along the directions of three reciprocal lattice vectors G₁, G₂, and G₃, which together span a unit cell of the reciprocal lattice. The component of P_e^(λ) along for instance G₁ can be found from

$\mathbf {G} _{1}\cdot \mathbf {P} _{e}^{\left(\lambda \right)}=-{if|e| \over 8\pi ^{3}}\int _{A}dk_{2}dk_{3}\sum _{n=1}^{M}\int _{0}^{|\mathbf {G} _{1}|}dk_{1}\langle u_{n\mathbf {k} }^{\left(\lambda \right)}|\partial /\partial k_{1}|u_{n\mathbf {k} }^{\left(\lambda \right)}\rangle$

where the two-dimensional integral is taken over the area A, spanned by G₂ and G₃, and the line integral runs over a line segment parallel to G₁. Interchanging the indices 1, 2 and 3 in the equation above yields the expressions for two other components of P_e^(λ).

Thus the electronic part of the polarization P_e^(λ) is given (modulo feR/Ω₀) by the sum

$\sum _{i=1}^{3}(\mathbf {P} _{e}^{\left(\lambda \right)})_{i}=\sum _{i=1}^{3}\left(\mathbf {G} _{i}\cdot \mathbf {P} _{e}^{\left(\lambda \right)}\right){\mathbf {R} _{i} \over 2\pi }$

where the lattice vectors R_i obey the relationship R_i·G_j=2πδ_ij.

The integration over A in the above, is straightforward and can be performed by sampling a 2D Monkhorst-Pack mesh of k-points, termed the perpendicular mesh or k_⊥-mesh by King-Smith and Vanderbilt. However, to remove the influence of the random phase of the functions u^(λ)_nk, introduced by the diagonalization routine, King-Smith and Vanderbilt propose to replace the line integral alias integration in the parallel or G_|| direction by,

$\phi _{J}^{\left(\lambda \right)}\left(\mathbf {k} _{\perp }\right)=\mathrm {Im} \left\{\ln \prod _{j=0}^{J-1}\mathrm {det} \left(\langle u_{m\mathbf {k} _{j}}^{\left(\lambda \right)}|u_{n\mathbf {k} _{j+1}}^{\left(\lambda \right)}\rangle \right)\right\}$

which is evaluated by calculating the cell-periodic parts of the wave functions at a string of J k-points, k_j= k_⊥+jG_||/J (with j=0,..,J-1), and where for sufficiently large J one has that

$\phi _{J}^{\left(\lambda \right)}\left(\mathbf {k} _{\perp }\right)=-i\sum _{n=1}^{M}\int _{0}^{|\mathbf {G} _{\parallel }|}dk_{\parallel }\langle u_{n\mathbf {k} }^{\left(\lambda \right)}|\partial /\partial k_{\parallel }|u_{n\mathbf {k} }^{\left(\lambda \right)}\rangle$

Note: the determinant appearing in the equation above is the determinant of the M×M matrix formed by letting n and m run over all valence bands.

The crucial step, instrumental in removing the random phase, is that the functions u^(λ)_nkJ are not obtained from an independent diagonalization, but found through their relationship with the functions u^(λ)_nk0,

$u_{n\mathbf {k} _{J}}^{\left(\lambda \right)}(\mathbf {r} )=e^{-i\mathbf {G} _{\parallel }\cdot \mathbf {r} }u_{n\mathbf {k} _{0}}^{\left(\lambda \right)}(\mathbf {r} )$

This way the product φ_J^(λ) becomes cyclic, and contains both u^(λ)_nkj as well as its complex conjugate for every k-point in the string, thus removing the random phase.

In practice G_||·P_e^(λ) is evaluated by way of the following summation over the k_⊥-mesh,

$(\mathbf {P} _{e}^{\left(\lambda \right)})_{i}={\frac {f|e|\mathbf {R} _{i}}{2\pi \Omega _{0}}}\left({\frac {1}{N_{k_{\perp }}}}\sum _{\mathbf {k} _{\perp }}\mathrm {Im} \ln {\frac {D_{J}^{\left(\lambda \right)}\left(\mathbf {k} _{\perp }\right)}{\langle D\rangle }}+\mathrm {Im} \ln \langle D\rangle \right)$

where

$D_{J}^{\left(\lambda \right)}\left(\mathbf {k} _{\perp }\right)=\prod _{j=0}^{J-1}\mathrm {det} \left(\langle u_{m\mathbf {k} _{j}}^{\left(\lambda \right)}|u_{n\mathbf {k} _{j+1}}^{\left(\lambda \right)}\rangle \right)$

with k_j= k_⊥+jG_||/J (with j=0,..,J-1), and

$\langle D\rangle ={\frac {1}{N_{k_{\perp }}}}\sum _{\mathbf {k} _{\perp }}D_{J}^{\left(\lambda \right)}\left(\mathbf {k} _{\perp }\right),$

and where we used R_i·G_j=2πδ_ij.

Assuming the D_J^(λ)(k_⊥) are reasonably well-clustered around ⟨D⟩, all terms D_J^(λ)(k_⊥)/⟨D⟩ will lie on the same branch of the logarithm. This makes it less likely that (P_e^(λ))_i will pick up a spurious contribution (of nR_i/N_k⊥).

Self-consistent response to finite electric fields

As of version 5.2, VASP can calculate the ground state of an insulating system under the application of a finite homogeneous electric field. The VASP implementation closely follows the PEAD (Perturbation Expression After Discretization) approach of Nunes and Gonze and the work of Souza et al..

In short: to determine the ground state of an insulating system under the application of a finite homogeneous electric field ε, VASP solves for the field-polarized Bloch functions {ψ^(ε)} by minimizing the electric enthalpy functional:

$E[\{\psi ^{({\mathcal {E}})}\},{\mathcal {E}}]=E_{0}[\{\psi ^{({\mathcal {E}})}\}]-\Omega {\mathcal {E}}\cdot \mathbf {P} [\{\psi ^{({\mathcal {E}})}\}],$

where P[{ψ^(ε)}] is the macroscopic polarization as defined in the modern theory of polarization:

$\mathbf {P} [\{\psi ^{({\mathcal {E}})}\}]=-{\frac {2ie}{(2\pi )^{3}}}\sum _{n}\int _{\mathrm {BZ} }d\mathbf {k} \langle u_{n\mathbf {k} }^{({\mathcal {E}})}|\nabla _{\mathbf {k} }|u_{n\mathbf {k} }^{({\mathcal {E}})}\rangle$

and u^(ε)_nk is the cell-periodic part of ψ^(ε)_nk.

The second term on the right-hand side of the electric enthalpy functional introduces a corresponding additional term to the Hamiltonian

$H|\psi _{n\mathbf {k} }^{({\mathcal {E}})}\rangle =H_{0}|\psi _{n\mathbf {k} }^{({\mathcal {E}})}\rangle -\Omega {\mathcal {E}}\cdot {\frac {\delta \mathbf {P} \left[\{\psi ^{({\mathcal {E}})}\}\right]}{\delta \langle \psi _{n\mathbf {k} }^{({\mathcal {E}})}|}}.$

Following the work of Nunes and Gonze we write,

${\frac {\delta \mathbf {P} \left[\{\psi ^{({\mathcal {E}})}\}\right]}{\delta \langle \psi _{n\mathbf {k} }^{({\mathcal {E}})}|}}=-{\frac {ie}{2\Delta k}}\sum _{m=1}^{N}\left[|u_{m\mathbf {k} _{j+1}}^{({\mathcal {E}})}\rangle S_{mn}^{-1}(\mathbf {k} _{j},\mathbf {k} _{j+1})-|u_{m\mathbf {k} _{j-1}}^{({\mathcal {E}})}\rangle S_{mn}^{-1}(\mathbf {k} _{j},\mathbf {k} _{j-1})\right]$

where m runs over the N occupied bands of the system, Δk=|k_j+1-k_j|, and

$S_{nm}(\mathbf {k} _{j},\mathbf {k} _{j+1})=\langle u_{n\mathbf {k} _{j}}^{({\mathcal {E}})}|u_{m\mathbf {k} _{j+1}}^{({\mathcal {E}})}\rangle .$

This Hamiltonian allows one to solve for {ψ^(ε)} by means of a direct optimization method.

Note: By analogy, it can be shown that

${\frac {\partial |u_{n\mathbf {k} _{j}}\rangle }{\partial k}}={\frac {ie}{2\Delta k}}\sum _{m=1}^{N}\left[|u_{m\mathbf {k} _{j+1}}^{({\mathcal {E}})}\rangle S_{mn}^{-1}(\mathbf {k} _{j},\mathbf {k} _{j+1})-|u_{m\mathbf {k} _{j-1}}^{({\mathcal {E}})}\rangle S_{mn}^{-1}(\mathbf {k} _{j},\mathbf {k} _{j-1})\right]$

in the sense of a first-order finite difference scheme (higher-order stencils may be similarly defined).

Note: One should be aware that when the electric field is chosen to be too large, the electric enthalpy functional will lose its minima, and VASP will not be able to find a stationary solution for the field-polarized orbitals. This is discussed in some detail by Souza et al.. VASP will produce a warning if:

$e|{\mathcal {E}}\cdot \mathbf {a} _{i}|>{\frac {1}{10}}E_{\mathrm {gap} }/N_{i},$

where E_gap is the bandgap, a_i are the lattice vectors, and N_i is the number of k-points along the reciprocal lattice vector i, in the regular (N₁×N₂×N₃) k-mesh. The factor 1/10 is chosen to be on the safe side. If one does not include unoccupied bands, VASP is obviously not able to determine the bandgap and can not check whether the electric field might be too large. This will also produce a warning message.

Response properties

The change in the macroscopic polarization due to the electric field ε defines the ion-clamped static dielectric tensor

$\epsilon _{ij}^{\infty }=\delta _{ij}+{\frac {4\pi }{\epsilon _{0}}}{\frac {\partial P_{i}}{\partial {\mathcal {E}}_{j}}},\qquad {i,j=x,y,z},$

the change in the Hellmann-Feynman forces due to ε, the Born effective charge tensors

$Z_{ij}^{*}={\frac {\Omega }{e}}{\frac {\partial P_{i}}{\partial u_{j}}}={\frac {1}{e}}{\frac {\partial F_{i}}{\partial {\mathcal {E}}_{j}}},\qquad {i,j=x,y,z},$

and the ion-clamped piezoelectric tensor of the system

$e_{ij}^{(0)}=-{\frac {\partial \sigma _{i}}{\partial {\mathcal {E}}_{j}}},\qquad {i=xx,yy,zz,xy,yz,zx}\quad {j=x,y,z},$

is found as the change in the stress tensor.

Related Tags and Sections

LCALCPOL, LCALCEPS, EFIELD_PEAD, LPEAD, IPEAD, LBERRY, IGPAR, NPPSTR, DIPOL

---