diff --git a/doc/sphinx/source/theory/electronBTE.rst b/doc/sphinx/source/theory/electronBTE.rst index 50b0d186..ed384c44 100644 --- a/doc/sphinx/source/theory/electronBTE.rst +++ b/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). `_, 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'} diff --git a/src/harmonic/phonon_h0.cpp b/src/harmonic/phonon_h0.cpp index fe75b66a..dd7d0086 100644 --- a/src/harmonic/phonon_h0.cpp +++ b/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, 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 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; } diff --git a/src/harmonic/phonon_h0_kokkos.cpp b/src/harmonic/phonon_h0_kokkos.cpp index d948d37c..cba9ea17 100644 --- a/src/harmonic/phonon_h0_kokkos.cpp +++ b/src/harmonic/phonon_h0_kokkos.cpp @@ -378,15 +378,19 @@ PhononH0::kokkosBatchedDiagonalizeWithVelocities( int numK = cartesianCoordinates.extent(0); double delta = 1.0e-8; + Kokkos::complex complexI(0.0, 1.0); // prepare all the wavevectors at which we need the hamiltonian, // (kx, ky, kz), (kx+dk, ky, kz), (kx - dk, ky, kz), (kx, ky + dk, kz), etc. + // 7 here is for the seven possible k+/- values DoubleView2D allVectors("wavevectors_k", numK * 7, 3); Kokkos::parallel_for( - "elHamiltonianShifted_d", numK, KOKKOS_LAMBDA(int iK) { + "phHamiltonianShifted_d", numK, KOKKOS_LAMBDA(int iK) { + // central points, i = x,y,z for (int i = 0; i < 3; ++i) { allVectors(iK * 7, i) = cartesianCoordinates(iK, i); } + // iDir for each kx, ky, kz +/-, i = x,y,z for (int iDir = 0; iDir < 3; ++iDir) { for (int i = 0; i < 3; ++i) { allVectors(iK * 7 + iDir * 2 + 1, i) = cartesianCoordinates(iK, i); @@ -411,66 +415,119 @@ PhononH0::kokkosBatchedDiagonalizeWithVelocities( numBands, 1, numBands, numBands); - // save energies and eigenvectors to results - DoubleView2D resultEnergies("energies", numK, numBands); StridedComplexView3D resultEigenvectors("eigenvectors", eigenvectorLayout); ComplexView4D resultVelocities("velocities", numK, numBands, numBands, 3); + ComplexView2D phases("wallacePhases", numK, numBands); + // copy slice of eigenvectors, energies into views + Kokkos::parallel_for( + "phase", Range3D({0, 0, 0}, {numK, numAtoms, 3}), KOKKOS_LAMBDA(int iK, int iat, int i) { + // prepare the Wallace phase + int m = i + iat * 3; // generate the band index + auto R = atomicPositions.row(iat); // TODO need atomic positions as a view... + + // i for each kx, ky, kz +/-, j = x,y,z + double arg = 0.0; + for(int j = 0; j < 3; ++j) { + arg += R(j) * allVectors(iK * 7 + i * 2 + 1, j); + } + phases(iK,m) = exp( -complexI * arg ); + }); + // apply phases -- here we use the results eigenvectors container + // to hold phase * U(k), to be replaced with standard ones before return Kokkos::parallel_for( "der", Range2D({0, 0}, {numK, numBands}), KOKKOS_LAMBDA(int iK, int m) { - auto E = Kokkos::subview(allEnergies, iK * 7, Kokkos::ALL); + auto X = Kokkos::subview(allEigenvectors, iK * 7, Kokkos::ALL, Kokkos::ALL); - resultEnergies(iK, m) = E(m); + + // copy resultEigenvectors to return at the end for (int n=0; n x(0.,0.); for (int l=0; l tmp(0.,0.); for (int l = 0; l < numBands; ++l) { - tmp += Kokkos::conj(L(l,m)) * 0.5 * (R(l, n) + Kokkos::conj(R(n, l))); + // average (der^* + der / 2) to enforce Hermiticity + //tmp += Kokkos::conj(L(l,m)) * 0.5 * (R(l, n) + Kokkos::conj(R(n, l))); + tmp += Kokkos::conj(U(l,m)) * dER(l, n); } - A(m, n) = tmp; + tempV(m, n) = tmp; }); Kokkos::fence(); + // then applying v_final = v1 * U(k) Kokkos::parallel_for( "vel", Range3D({0, 0, 0}, {numK, numBands, numBands}), KOKKOS_LAMBDA(int iK, int m, int n) { double norm = 0.; + // TODO why skip gamma here but not above? for (int i=0; i<3; ++i) { norm += cartesianCoordinates(iK,i) * cartesianCoordinates(iK,i); } @@ -488,14 +545,24 @@ PhononH0::kokkosBatchedDiagonalizeWithVelocities( Kokkos::resize(der, 0, 0, 0); Kokkos::resize(tmpV, 0, 0, 0); - //print4DComplex("new = ", resultVelocities); + // copy resultEigenvectors and energies to return + DoubleView2D resultEnergies("energies", numK, numBands); + Kokkos::parallel_for( + "copyEigenvectors", Range2D({0, 0}, {numK, numBands}), KOKKOS_LAMBDA(int iK, int m) { + auto E = Kokkos::subview(allEnergies, iK * 7, Kokkos::ALL); + auto X = Kokkos::subview(allEigenvectors, iK * 7, Kokkos::ALL, Kokkos::ALL); + + resultEnergies(iK, m) = E(m); + + for (int n=0; n