Fix implementation of Wigner thermal conductivity

doc/sphinx/source/theory/electronBTE.rst

@@ -47,7 +47,7 @@ with
    \times
    \bigg[
    (1-\bar{f}_{\boldsymbol{k}'b'} + \bar{n}_{\boldsymbol{q}\nu})
-   \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} - \hbar \omega_{\boldsymbol{q}\nu}) \\ 
+   \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} - \hbar \omega_{\boldsymbol{q}\nu}) \\
    &+
    (\bar{f}_{\boldsymbol{k}'b'} + \bar{n}_{\boldsymbol{q}\nu})
    \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} + \hbar \omega_{\boldsymbol{q}\nu})
@@ -75,7 +75,7 @@ which results in the matrix with diagonal matrix elements:
    \Omega_{\lambda \lambda} = \frac{1}{\tau_{\boldsymbol{k}b}}
 
 and for the off-diagonal terms:
-   
+
 .. math::
    \tilde{\Omega}_{\boldsymbol{k}b,\boldsymbol{k}'b'} =&
    -
@@ -93,11 +93,11 @@ This matrix is symmetric and has a number of interesting physical properties (e.
 Computationally, the symmetric matrix can be used in a conjugate gradient method that maximises the electrical and thermal conductivity, and guarantees the existence of eigenvalues.
 
 In Phoebe, instead of solving the original BTE problem in the form :math:`\sum_{\lambda'} \Omega_{\lambda,\lambda'} \delta f_{\lambda'} = b_{\lambda}`, we solve the symmetrized problem:
-   
+
 .. math::
    \sum_{\lambda'} \tilde{\Omega}_{\lambda,\lambda'} \delta \tilde{f}_{\lambda'} = \tilde{b}_{\lambda}
 
-with 
+with
 
 .. math::
    \delta f_{\lambda} = ( \bar{f}_{\lambda} (1-\bar{f}_{\lambda}) )^{-\frac{1}{2}} \delta f_{\lambda}
@@ -108,12 +108,12 @@ and
    \tilde{b}_{\lambda} = ( \bar{f}_{\lambda} (1-\bar{f}_{\lambda}) )^{-\frac{1}{2}} b_{\lambda}
 
 
-   
+
 
 Onsager coefficients
 --------------------
 
-In the electronic case, there are a handful of transport coefficients we'd like to solve the BTE to find, such as the electrical conductivity, :math:`\sigma`, the electronic part of the thermal conductivity, :math:`\kappa_e`, and the Seebeck coeffieint, :math:`S`. These quantities are defined in terms of the Onsager coefficients. 
+In the electronic case, there are a handful of transport coefficients we'd like to solve the BTE to find, such as the electrical conductivity, :math:`\sigma`, the electronic part of the thermal conductivity, :math:`\kappa_e`, and the Seebeck coeffieint, :math:`S`. These quantities are defined in terms of the Onsager coefficients.
 
 We assume that the response to the applied electric field and thermal gradient is linear in these external fields:
 
@@ -133,7 +133,7 @@ and
 
 where :math:`g_s` is the spin degeneracy.
 
-We can decompose these to write, 
+We can decompose these to write,
 
 .. math::
    \boldsymbol{J} = L_{EE} \boldsymbol{E} + L_{ET} \boldsymbol{\nabla} T
@@ -162,7 +162,7 @@ where :math:`n` is the carrier concentration.
 Solutions of the electron BTE
 --------------------------------------
 
-Largely, these solvers follow the equivalent section in the phonon BTE section, where they may be described in more detail. For further details and references on any specific solver, we suggest you visit the equivalent phonon sections, as well. Here, we again establish methods of finding the solution vector to the BTE, :math:`f`, but in this case, we have two: :math:`f^T` and :math:`f^E`, for each field. 
+Largely, these solvers follow the equivalent section in the phonon BTE section, where they may be described in more detail. For further details and references on any specific solver, we suggest you visit the equivalent phonon sections, as well. Here, we again establish methods of finding the solution vector to the BTE, :math:`f`, but in this case, we have two: :math:`f^T` and :math:`f^E`, for each field.
 
 RTA Solution
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -216,6 +216,8 @@ Iterative solution
 .. note::
    Generally, we recommend the variational method over this for numerical stability. However, the iterative solution allows us to use symmetries as in `Chaput, PRL (2013). <https://doi.org/10.1103/PhysRevLett.110.265506>`_, which can decrease cost.  
 
+This is an adaptation of the Omini-Sparavigna method to electrons. To better understand this method, please have a look first at the counterpart phonon section.
+
 To better understand this method, please have a look first at the counterpart iterative solution in the phonon documentation. 
 In short, the electron BTE consists in two linear algebra problems:
 
@@ -352,23 +354,23 @@ The related transport coefficients are defined as:
 
 .. math::
    L_{EE}^{ij} =
-   \frac{e g_s}{V N_k} \sum_{\boldsymbol{k}b} \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{E_j}(\boldsymbol{k}) \Big\}_{bb}
+   \frac{e g_s}{V N_k} \sum_{\boldsymbol{k}b} \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{E_j}(\boldsymbol{k}) \Big\}_{bb'}
 
 .. math::
    L_{ET}^{ij} =
-   \frac{e g_s}{V N_k} \sum_{\boldsymbol{k}b} \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{T_j}(\boldsymbol{k}) \Big\}_{bb}
+   \frac{e g_s}{V N_k} \sum_{\boldsymbol{k}b} \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{T_j}(\boldsymbol{k}) \Big\}_{bb'}
 
 .. math::
    L_{TE}^{ij} =
    \frac{g_s}{V N_k}
    \sum_{\boldsymbol{k}b}
    \big( \epsilon_{b}(\boldsymbol{k})-\mu \big)
-   \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{E_j}(\boldsymbol{k}) \Big\}_{bb}
+   \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{E_j}(\boldsymbol{k}) \Big\}_{bb'}
 
 .. math::
    L_{TT}^{ij} =
    \frac{g_s}{V N_k}
    \sum_{\boldsymbol{k}b}
    \big( \epsilon_{b}(\boldsymbol{k})-\mu \big)
-   \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{T_j}(\boldsymbol{k}) \Big\} _{bb}
+   \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{T_j}(\boldsymbol{k}) \Big\} _{bb'}
 

src/harmonic/phonon_h0.cpp

@@ -639,6 +639,8 @@ PhononH0::diagonalizeVelocity(Point &point) {
   return diagonalizeVelocityFromCoordinates(coordinates);
 }
 
+// TODO I think it would be easier to 
+// apply these phases to the dynmat during the transform
 Eigen::Tensor<std::complex<double>, 3>
 PhononH0::diagonalizeVelocityFromCoordinates(Eigen::Vector3d &coordinates) {
 
@@ -654,12 +656,29 @@ PhononH0::diagonalizeVelocityFromCoordinates(Eigen::Vector3d &coordinates) {
   auto energies = std::get<0>(tup);
   auto eigenvectors = std::get<1>(tup);
 
+  Eigen::MatrixXcd eigenvectors_bar(numBands, numBands);
+  Eigen::MatrixXcd eigenvectors_bar_minus(numBands, numBands);
+  Eigen::MatrixXcd eigenvectors_bar_plus(numBands, numBands);
+  eigenvectors_bar.setZero();
+  eigenvectors_bar_minus.setZero();
+  eigenvectors_bar_plus.setZero();
+
+  Eigen::VectorXcd matrix_U(numBands);
+  Eigen::VectorXcd matrix_U_plus(numBands);
+  Eigen::VectorXcd matrix_U_minus(numBands);
+  matrix_U.setZero();
+  matrix_U_plus.setZero();
+  matrix_U_minus.setZero();
+
   // now we compute the velocity operator, diagonalizing the expectation
   // value of the derivative of the dynamical matrix.
   // This works better than doing finite differences on the frequencies.
   double deltaQ = 1.0e-8;
+  // kx, ky, kz directions 
   for (int i : {0, 1, 2}) {
+
     // define q+ and q- from finite differences.
+    Eigen::Vector3d qcenter = coordinates;
     Eigen::Vector3d qPlus = coordinates;
     Eigen::Vector3d qMinus = coordinates;
     qPlus(i) += deltaQ;
@@ -673,6 +692,27 @@ PhononH0::diagonalizeVelocityFromCoordinates(Eigen::Vector3d &coordinates) {
     auto enMinus = std::get<0>(tup1);
     auto eigMinus = std::get<1>(tup1);
 
+    // apply the Wallace phase 
+    for (int iat = 0; iat < numAtoms; iat++){
+      for (int i : {0, 1, 2}) {
+        auto ip = i + iat * 3;
+        double arg = qcenter.dot(atomicPositions.row(iat));
+        matrix_U(ip) = exp(-complexI*arg);
+        arg = qPlus.dot(atomicPositions.row(iat));
+        matrix_U_plus(ip) = exp(-complexI*arg);
+        arg = qMinus.dot(atomicPositions.row(iat));
+        matrix_U_minus(ip) = exp(-complexI*arg);
+      }
+    }
+
+    for(int i = 0; i < numBands; i++){
+      for(int j = 0; j < numBands; j++){
+        eigenvectors_bar(i,j) = matrix_U(i) * eigenvectors(i,j);
+        eigenvectors_bar_plus(i,j) = matrix_U_plus(i) * eigPlus(i,j);
+        eigenvectors_bar_minus(i,j) = matrix_U_minus(i) * eigMinus(i,j);
+      }
+    }
+
     // build diagonal matrices with frequencies
     Eigen::MatrixXd enPlusMat(numBands, numBands);
     Eigen::MatrixXd enMinusMat(numBands, numBands);
@@ -684,19 +724,18 @@ PhononH0::diagonalizeVelocityFromCoordinates(Eigen::Vector3d &coordinates) {
     // build the dynamical matrix at the two wavevectors
     // since we diagonalized it before, A = M.U.M*
     Eigen::MatrixXcd sqrtDPlus(numBands, numBands);
-    sqrtDPlus = eigPlus * enPlusMat * eigPlus.adjoint();
+    sqrtDPlus = eigenvectors_bar_plus * enPlusMat * eigenvectors_bar_plus.adjoint();
     Eigen::MatrixXcd sqrtDMinus(numBands, numBands);
-    sqrtDMinus = eigMinus * enMinusMat * eigMinus.adjoint();
+    sqrtDMinus = eigenvectors_bar_minus * enMinusMat * eigenvectors_bar_minus.adjoint();
 
     // now we can build the velocity operator
     Eigen::MatrixXcd der(numBands, numBands);
     der = (sqrtDPlus - sqrtDMinus) / (2. * deltaQ);
 
-    // and to be safe, we reimpose hermiticity
-    der = 0.5 * (der + der.adjoint());
-
-    // now we rotate in the basis of the eigenvectors at q.
-    der = eigenvectors.adjoint() * der * eigenvectors;
+    // option below can probably be decommented, this will enforce hermiticity numerically, 
+    // results should be equivalent within numerical noise
+    //der = 0.5 * (der + der.adjoint()); 
+    der = eigenvectors_bar.adjoint() * der * eigenvectors_bar;
 
     // copy to output
     for (int ib2 = 0; ib2 < numBands; ib2++) {
@@ -705,65 +744,6 @@ PhononH0::diagonalizeVelocityFromCoordinates(Eigen::Vector3d &coordinates) {
       }
     }
   }
-
-  // turns out that the above algorithm has problems with degenerate bands
-  // so, we diagonalize the velocity operator in the degenerate subspace,
-
-  for (int ib = 0; ib < numBands; ib++) {
-    // first, we check if the band is degenerate, and the size of the
-    // degenerate subspace
-    int sizeSubspace = 1;
-    for (int ib2 = ib + 1; ib2 < numBands; ib2++) {
-      // I consider bands degenerate if their frequencies are the same
-      // within 0.0001 cm^-1
-      if (abs(energies(ib) - energies(ib2)) > 0.0001 / ryToCmm1) {
-        break;
-      }
-      sizeSubspace += 1;
-    }
-
-    if (sizeSubspace > 1) {
-      Eigen::MatrixXcd subMat(sizeSubspace, sizeSubspace);
-      // we have to repeat for every direction
-      for (int iCart : {0, 1, 2}) {
-        // take the velocity matrix of the degenerate subspace
-        for (int j = 0; j < sizeSubspace; j++) {
-          for (int i = 0; i < sizeSubspace; i++) {
-            subMat(i, j) = velocity(ib + i, ib + j, iCart);
-          }
-        }
-
-        // reinforce hermiticity
-        subMat = 0.5 * (subMat + subMat.adjoint());
-
-        // diagonalize the subMatrix
-        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigenSolver(subMat);
-        Eigen::MatrixXcd newEigenVectors = eigenSolver.eigenvectors();
-        //newEigenVectors = 3*subMat; // TODO: undo
-
-        // rotate the original matrix in the new basis
-        // that diagonalizes the subspace.
-        subMat = newEigenVectors.adjoint() * subMat * newEigenVectors;
-
-        // reinforce hermiticity
-        subMat = 0.5 * (subMat + subMat.adjoint());
-
-        // substitute back
-        for (int j = 0; j < sizeSubspace; j++) {
-          for (int i = 0; i < sizeSubspace; i++) {
-            velocity(ib + i, ib + j, iCart) = subMat(i, j);
-          }
-        }
-      }
-    }
-
-    // we skip the bands in the subspace, since we corrected them already
-    ib += sizeSubspace - 1;
-  }
-  // if we are working at gamma, we set all velocities to zero.
-  if (coordinates.norm() < 1.0e-6) {
-    velocity.setZero();
-  }
   Kokkos::Profiling::popRegion(); // diagonalizeVelocityFromCoordinates
   return velocity;
 }