Merge branch 'develop' into mpi_plot_wannier

.github/workflows/main.yml

@@ -3,10 +3,6 @@ name: CI
 on:
   pull_request:
   push:
-    branches:
-    - develop
-    - master
-    - gh-actions
 
 jobs:
   pre-commit:

doc/tutorial/tutorial.tex

@@ -3672,6 +3672,130 @@ \subsection*{Notes}
 \label{fig:W_fit}
 \end{figure}
 
+\sectiontitle{33: Monolayer BC$_2$N -- $k\cdot p$ expansion coefficients}
+
+\begin{itemize}
+
+\item Outline: \textit{Calculate $k\cdot p$ expansion coefficients monolayer BC$_2$N using
+quasi-degenerate (L\"owdin) 
+perturbation theory. In preparation for this example it may be useful to read Ref.
+\cite{ibanez-azpiroz-ArXiv2019} }
+
+
+
+\item Directory: \verb|examples/example33/|
+
+\item Input files:
+
+\begin{itemize}
+
+\item[--] \verb|bc2n.scf| \textit{The {\tt PWSCF} input file for ground state calculation}
+\item[--] \verb|bc2n.nscf| \textit{The {\tt PWSCF} input file to obtain Bloch states on a uniform grid}
+\item[--] \verb|bc2n.pw2wan| \textit{The input file for} \verb|pw2wannier90|
+\item[--] \verb|bc2n.win| \textit{The} \verb|wannier90| \textit{and} \verb|postw90| \textit{input file}
+
+
+\end{itemize}
+
+
+\begin{enumerate}
+
+\item Run {\tt PWSCF} to obtain the ground state of Gallium Arsenide
+
+\verb|pw.x < bc2n.scf > scf.out|
+
+
+\item Run {\tt PWSCF} to obtain the ground state of Gallium Arsenide
+
+\verb|pw.x < bc2n.nscf > nscf.out|
+
+\item Run {\tt Wannier90} to generate a list of the required overlaps (written into the \verb|GaAs.nnkp| file)
+
+\verb|wannier90.x -pp bc2n|
+
+
+\item Run {\tt pw2wannier90} to compute:
+
+\begin{itemize}
+\item[--] The overlaps $\langle u_{n\bf{k}}|u_{n\bf{k+b}}\rangle$ between spinor 
+Bloch states (written in the \verb|bc2n.mmn| file) 
+\item[--] The projections for the starting guess (written in the \verb|bc2n.amn| file)
+
+\end{itemize}
+
+
+\verb|pw2wannier90.x < bc2n.pw2wan > pw2wan.out|
+
+\item Run {\tt wannier90} to  compute MLWFs
+
+\verb|wannier90.x bc2n|
+
+\item Run {\tt postw90} to compute expansion coefficients
+
+\verb|postw90.x bc2n| 
+
+
+\end{enumerate}
+
+
+\subsection*{Expansion coefficients}
+
+For computing $k\cdot p$ expansion coefficients as given by quasi-degenerate (L\"owdin) 
+perturbation theory, set
+\begin{verbatim}
+berry = true
+berry_task = kdotp
+\end{verbatim}
+Select the k-point around which the expansion coefficients will be computed, \textit{e.g.}, the S point
+\begin{verbatim}
+kdotp_kpoint  =  0.5000 0.0000 0.5000
+\end{verbatim}
+Set number of bands that should be taken into account for the $k\cdot p$ expansion, as well
+as their band indexes within the Wannier basis
+\begin{verbatim}
+kdotp_num_bands = 2
+kdotp_bands =  2,3
+\end{verbatim}
+Since no k-space integral is needed, set
+\begin{verbatim}
+berry_kmesh = 1 1 1
+\end{verbatim}
+Although not used, we also need to input the value of the Fermi level in eV
+\begin{verbatim}
+fermi_energy = [insert your value here]
+\end{verbatim}
+
+
+
+On output, the program generates three files, namely \verb|SEED-kdotp_0.dat|, 
+\verb|SEED-kdotp_1.dat| and  \verb|SEED-kdotp_2.dat|, which correspond to the zeroth, first and second 
+order expansion coefficients, respectively. The dimension of the matrix contained in each file 
+is $3^{l}\times N^{2}$, where $N$ is the number of bands set by \verb|kdotp_num_bands|, and $l$
+is the order of the expansion term (currently $l=0,1$ or $2$). 
+
+These coefficients can be used, among other things, to compute the energy dispersion of
+the bands of interest around the chosen k-point. 
+The $k\cdot p$ band dispersion can be computed and plotted along $k_x$ (from S to X)
+using python and the file \verb|kdotp_plot.py| provided in the example folder
+\begin{verbatim}
+python kdotp_plot.py
+\end{verbatim}
+For comparison, the exact band structure calculated using Wannier90 (file \verb|bc2n_band.dat|, generated automatically) 
+is also plotted along (see Fig. \ref{fig:bc2n-bnd}).
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=12cm]{kdotp_bands_SX.pdf}
+\caption{Band dispersion of monolayer BC$_{2}$N around $S$ point. Exact results (solid dots)
+are compared to first-order (blue) and second-order (red) $k\cdot p$ model results for valence and conduction bands.} 
+\label{fig:bc2n-bnd}
+\end{center}
+\end{figure}
+
+\end{itemize}
+
+
 %\cleardoublepage
 
 \bibliographystyle{apsrev4-1}

doc/user_guide/berry.tex

@@ -250,9 +250,10 @@ \section{{\tt berry\_task=morb}: orbital magnetization}
 factor of $-1/2$ from the one in Eq.~(97) and Fig.~2 of that work.
 
 \section{{\tt berry\_task=shc}: spin Hall conductivity}
+\label{sec:shc}
 
 The Kubo-Greenwood formula for the intrinsic spin Hall conductivity (SHC) of a crystal
-in the independent-particle approximation reads \cite{qiao-prb2018,guo-prl2008}
+in the independent-particle approximation reads \cite{qiao-prb2018,ryoo-prb2019,guo-prl2008}
 %
 \begin{equation}
 \label{eq:kubo_shc}
@@ -275,7 +276,11 @@ \section{{\tt berry\_task=shc}: spin Hall conductivity}
 
 The velocity matrix element in the numerator is the same as   Eq.~(\ref{eq:velocity_mat}), 
 so the only unknown quantity is the spin current matrix $\langle n\bm{k}| \hat{j}_{\alpha}^{\gamma}|m\bm{k}\rangle$. 
-We can use Wannier interpolation technique to efficiently calculate this matrix, for a full derivation please refer to Ref.~\cite{qiao-prb2018}. 
+We can use Wannier interpolation technique to efficiently calculate this matrix, and there are two derivation according to the degree of approximation. A noteworthy difference is the way in which two {\it ab-initio} matrix elements are evaluated,
+\begin{equation}
+  \langle u_{n{\bf k}}\vert\sigma_\gamma H_{\bf k}\vert u_{m{\bf k}+{\bf b}}\rangle, \langle u_{n{\bf k}}\vert\sigma_\gamma \vert u_{m{\bf k}+{\bf b}}\rangle, \gamma = x, y, z
+\end{equation}
+These are evaluated by {\tt pw2wannier90} using Ryoo's method. In contrast, Qiao's method does not require {\tt pw2wannier90}, but it assumes an approximation $1\approx\sum_{ l\in ab-initio{\rm \,bands}}|u_{l\bm{k}}\rangle \langle u_{l\bm{k}}|$. You can choose which method to evaluate this value with {\tt shc\_method} in the input file. For a full derivation please refer to Ref.~\cite{qiao-prb2018} or Ref.~\cite{ryoo-prb2019}.
 
 The Eq.~(\ref{eq:kubo_shc}) can be further separated into 
 band-projected Berry curvature-like term
@@ -322,12 +327,166 @@ \section{{\tt berry\_task=shc}: spin Hall conductivity}
 be employed by the keyword {\tt [kubo\_]adpt\_smr} 
 (see Examples 29 and 30 in the Tutorial). 
 
-Please cite the following paper~\cite{qiao-prb2018} when publishing SHC results obtained using this method:
+Please cite the following paper~\cite{qiao-prb2018} or ~\cite{ryoo-prb2019} when publishing SHC results obtained using this method:
 \begin{quote}
 	Junfeng Qiao, Jiaqi Zhou, Zhe Yuan, and Weisheng Zhao, \\
 	\emph{Calculation of intrinsic spin Hall conductivity by Wannier interpolation},\\
 	Phys. Rev. B. 98, 214402 (2018), DOI:10.1103/PhysRevB.98.214402.
 \end{quote}
+or
+\begin{quote}
+	Ji Hoon Ryoo, Cheol-hwan Park, and Ivo Souza, \\
+	\emph{Computation of intrinsic spin Hall conductivities from first principles using maximally localized Wannier functions},\\
+	Phys. Rev. B. 99, 235113 (2019), DOI:10.1103/PhysRevB.99.235113.
+\end{quote}
+
+
+\section{{\tt berry\_task=sc}: shift current}
+
+
+The shift-current contribution to the second-order response
+is characterized by a frequency-dependent third-rank tensor~\cite{sipe-prb00}
+\begin{equation}\label{eq:shiftcurrent}
+\begin{split}
+\sigma^{abc}(0;\omega,-\omega)=&-\frac{i\pi e^3}{4\hbar^2 \Omega_c N_k}
+\sum_{\bm{k}} \sum_{n,m}(f_{n\bm{k}}-f_{m\bm{k}})
+\times
+\left(r^b_{ mn}(\bm{k})r^{c;a}_{nm}(\bm{k}) + r^c_{mn}(\bm{k})r^{b;a}_{ nm}(\bm{k})\right)\\
+&\times \left[\delta(\omega_{mn\bm{k}}-\omega)+\delta(\omega_{nm\bm{k}}-\omega)\right],
+\end{split}
+\end{equation}
+where $a,b,c$ are spatial indexes
+and $\omega_{mn\bm{k}}=(\epsilon_{n\bm{k}}-\epsilon_{m\bm{k}})/\hbar$.
+The expression in Eq.~\ref{eq:shiftcurrent} involves 
+the dipole matrix element 
+\begin{equation}
+\label{eq:r}
+r^a_{ nm}(\bm{k})=(1-\delta_{nm})A^a_{ nm}(\bm{k}),
+\end{equation}
+and its \emph{generalized
+derivative}
+%
+\begin{equation}
+\label{eq:gen-der}
+r^{a;b}_{nm}(\bm{k})=\partial_{k_{b}} r^a_{nm}(\bm{k})
+-i\left(A^b_{nn}(\bm{k})-A^b_{ mm}(\bm{k})\right)r^a_{ nm}(\bm{k}).
+\end{equation}
+The first-principles evaluation of the
+above expression is technically challenging
+due to the presence of an extra $k$-space derivative.
+The implementation in {\tt wannier90} follows the scheme proposed in Ref.~\cite{ibanez-azpiroz_ab_2018}, 
+following the spirit of the Wannier-interpolation 
+method for calculating the AHC~\cite{wang-prb06} by reformulating $k\cdot p$
+perturbation theory within the subspace of wannierized bands.
+This strategy inherits the practical advantages of the sum-over-states
+approach, but without introducing the truncation errors usually associated with this procedure~\cite{sipe-prb00}.
+
+As in the case of the optical conductivity, a broadened delta function can be 
+applied in Eq.~\ref{eq:shiftcurrent} by means of the parameter
+$\eta$ (see Eq.~\ref{eq:lorentzian}) using the keyword 
+{\tt[kubo\_]smr\_fixed\_en\_width}, and adaptive smearing can 
+be employed using the keyword {\tt [kubo\_]adpt\_smr}. 
+
+Please cite Ref.~\cite{ibanez-azpiroz_ab_2018} when publishing shift-current results using this method.
+
+\section{{\tt berry\_task=kdotp}: $k\cdot p$ coefficients}
+\label{sec:kdotp}
+
+Consider a Hamiltonian
+\beq
+\label{eq:H}
+H=H^{0}+H^{\prime}
+\eeq
+where the eigenvalues $E_{n}$ and eigenfunctions $\ket{n}$ of $H^{0}$
+are known, and $H^{\prime}$ is a perturbation. In a nutshell, quasi-degenerate
+perturbation theory assumes that the set of eigenfunctions of $H^0$
+can be divided into subsets $A$ and $B$ that are weakly coupled by $H^{\prime}$,
+and that we are only interested in subset $A$.
+This theory asserts that a transformed Hamiltonian $\tilde{H}$ exists 
+within subspace $A$ such that
+\beq
+\label{eq:pert-exp}
+\tilde{H}=\tilde{H}^{0}+\tilde{H}^{1}+\tilde{H}^{2} + \cdots
+\eeq
+where $\tilde{H}^{j}$ contain matrix elements of $H^{\prime}$ to the $j$th power.
+According to Appendix B of Ref~\cite{winkler_spin-orbit_2003}, the first three terms are
+\bea
+\label{eq:pert-matelem0}
+& \tilde{H}^{0}_{mm'} = H^{0}_{mm'},\\
+\label{eq:pert-matelem1}
+& \tilde{H}^{1}_{mm'} = H^{'}_{mm'},\\
+\label{eq:pert-matelem2}
+& \tilde{H}^{2}_{mm'} = \dfrac{1}{2}\sum_{l}H^{'}_{ml}H^{'}_{lm'}
+\left( 
+\dfrac{1}{E_{m}-E_{l}}+\dfrac{1}{E_{m'}-E_{l}}
+\right),
+\eea
+where  $m,m'\in A$ and $l\in B$.
+The approximation $\tilde{H}\sim \tilde{H}^{0}+\tilde{H}^{1}$ amounts to truncating 
+$H$ to the $A$ subspace. By further including $\tilde{H}^{2}$, the coupling to the $B$
+subspace is incorporated approximately, ``renormalizing'' the elements of the truncated matrix.
+
+We adopt the notation described in Sec. III.B of Ref.~\cite{wang-prb06}.
+We shift the origin of $k$ space to the point where the band edge (or some other
+band extremum of interest) is located, and Taylor expand
+around that point the Wannier-gauge Hamiltonian,
+\beq\label{eq:HW-exp}
+H^{(W)}(\bm{k})=H^{(W)}(0)
++\sum_{a}H_{a}^{(W)}(0)k_{a}
++\dfrac{1}{2}\sum_{ab}H_{ab}^{(W)}(0)k_{a}k_{b}
++ \mathcal{O}(k^{3})
+\eeq
+where $a,b=x,y,z$, and 
+\bea
+&H_{a}^{(W)}(0)=\left. \dfrac{\partial H^{(W)}(\bm{k})}{\partial k_{a}}\right\rvert_{\bm{k}=0}\\
+&H_{ab}^{(W)}(0)=\left. \dfrac{\partial^{2} H^{(W)}(\bm{k})}{\partial k_{a}\partial k_{b}}\right\rvert_{\bm{k}=0}
+\eea
+
+We now apply to $H^{(W)}(\bm{k})$ a similarity transformation
+$U(0)$ that diagonalizes $H^{(W)}(0)$, and call the transformed Hamiltonian $H(\bm{k})$,
+\beq\label{eq:Hbar}
+H(\bm{k})=\overbrace{\overline{H}}^{H^{0}} + \overbrace{\sum_{a}\overline{H}_{a}k_{a}
++\dfrac{1}{2}\sum_{ab}\overline{H}_{ab}k_{a}k_{b}}^{H^{\prime}} + \mathcal{O}(k^{3}),
+\eeq
+where we introduced the notation
+\beq
+\overline{\mathcal{O}}=U^{\dagger}(0)\mathcal{O}^{(W)}(0)U(0),
+\eeq
+and applied it to $\mathcal{O}=H,{H}_{a},{H}_{ab}$.
+We can now apply quasi-degenerate perturbation theory by choosing the diagonal
+matrix $\overline{H}$ as our $H^{0}$, and the remaining (nondiagonal)
+terms in Eq. \ref{eq:Hbar} as $H^{\prime}$.
+Collecting terms in \eq{pert-exp} up to second order in $k$
+we get
+\beq\label{eq:Htilde}
+ \tilde{H}_{mm'}(\bm{k}) = 
+\overline{H}_{mm'} + \sum_{a} \left(\overline{H}_{a}\right)_{mm'}k_{a}
+ + \dfrac{1}{2}\sum_{a,b}\left[
+\left(\overline{H}_{ab}\right)_{mm'} + \left({T}_{ab}\right)_{mm'}
+\right]k_{a}k_{b}+ \mathcal{O}(k^{3}),
+\eeq
+where  $m,m'\in A$ and we have defined the virtual-transition matrix
+\beq\label{eq:Tab}
+\left({T}_{ab}\right)_{mm'}=\sum_{l\in B}
+\left(\overline{H}_{a}\right)_{ml}\left(\overline{H}_{b}\right)_{lm'} 
+ \times
+\left( 
+\dfrac{1}{E_{m}-E_{l}}+\dfrac{1}{E_{m'}-E_{l}}
+\right)
+= 
+\left({T}_{ab}\right)_{m'm}^{*}.
+\eeq 
+(The $T_{ab}$ term in Eq. \ref{eq:Htilde} gives an Hermitean contribution to 
+$\tilde{H}(\bm{k})$ only after summing over $a$ and $b$, whereas the other terms
+are Hermitean already before summing.)
+
+
+The implementation in {\tt wannier90} follows the scheme proposed in Ref.~\cite{ibanez-azpiroz-ArXiv2019}, 
+and outputs $\overline{H}_{mm'}$ in {\tt   seedname-kdotp\_0.dat}, 
+$\left(\overline{H}_{a}\right)_{mm'}$ in {\tt   seedname-kdotp\_1.dat}, and
+$\left[\left(\overline{H}_{ab}\right)_{mm'} + \left({T}_{ab}\right)_{mm'}\right]/2$ in {\tt   seedname-kdotp\_2.dat}.
+
+Please cite Ref.~\cite{ibanez-azpiroz-ArXiv2019} when publishing $k\cdot p$ results using this method.
 
 \section{Needed matrix elements}
 
@@ -386,8 +545,26 @@ \section{Needed matrix elements}
 	%\end{quote}
 }
 $$
-
-
+If one uses Ryoo's method to calculate spin Hall conductivity, the further matrix elements are needed:
+%
+$$\langle u_{n{\bf k}}\vert
+\sigma_\gamma H_{\bf k}\vert u_{m{\bf k}+{\bf b}}\rangle, \langle u_{n{\bf k}}\vert
+\sigma_\gamma \vert u_{m{\bf k}+{\bf b}}\rangle,
+\gamma = x, y, z
+$$
+%
+and these are evaluated by adding to the input file {\tt
+seedname.pw2wan} the lines
+%
+$${\tt
+	write\_sHu = .true.
+}
+$$
+$$
+{\tt
+	write\_sIu = .true.
+}
+$$
 
 %A note on the implementation: the optical conductivity and orbital
 %magnetization are implemented in the {\tt berry} module in the manner