add a version of enforce detailed balance to the phonon scattering matrix

Author: jcoulter12

src/algebra/PMatrix.cpp

@@ -174,12 +174,10 @@ ParallelMatrix<double>::elpaDiagonalize() {
   if( eigenvalues[0] < 0 ) { // negative modes were found
     if(mpi->mpiHead()) {
       Warning("Relaxons diagonalization found negative eigenvalues."
-                "\n\tThis can happen when there's a bit of numerical noise on the scattering matrix,"
-                "\n\tand finding them may indicate the calculation is unconverged."
+                "\n\tThis can happen when there's a numerical issue with the scattering matrix."
                 "\n\tWhile we simply do not include these when computing transport, "
-                "\n\tand likely if they are small the calculation will be unaffected, "
-                "\n\tyou may want to consider using with more wavevectors or an improved DFT calculation."
-                "\n\tAdditionally, setting symmetrizeMatrix = true in your input file will help.");
+                "\n\tand likely if they are very small the calculation will be unaffected, "
+                "\n\tplease view the relaxons tutorial for advice about resolving this.");
       std::cout << "These eigenvalues are (in atomic units):" << std::endl;
       for (int i = 0; eigenvalues[i] <= 0; i++) {
         std::cout << i << " " << eigenvalues[i] << std::endl;
@@ -292,11 +290,10 @@ ParallelMatrix<double>::scalapackDiagonalize() {
   if( eigenvalues[0] < 0 ) { // negative modes were found
     if(mpi->mpiHead()) {
       Warning("Relaxons diagonalization found negative eigenvalues."
-                "\n\tThis can happen when there's a bit of numerical noise on the scattering matrix,"
-                "\n\tand finding them may indicate the calculation is unconverged."
+                "\n\tThis can happen when there's a numerical issue with the scattering matrix."
                 "\n\tWhile we simply do not include these when computing transport, "
-                "\n\tyou may want to consider using with more wavevectors or an improved DFT calculation."
-                "\n\tAdditionally, setting symmetrizeMatrix = true in your input file will help.");
+                "\n\tand likely if they are very small the calculation will be unaffected, "
+                "\n\tplease view the relaxons tutorial for advice about resolving this.");
       std::cout << "These eigenvalues are (in atomic units):" << std::endl;
       for (int i = 0; eigenvalues[i] <= 0; i++) {
         std::cout << i << " " << eigenvalues[i] << std::endl;
@@ -486,11 +483,10 @@ std::tuple<std::vector<double>, ParallelMatrix<double>>
       if(mpi->mpiHead()) {
         Warning("Relaxons diagonalization found " + std::to_string(m) +
                " in the range -1 < eigenvalues <= 0."
-                  "\n\tThis can happen when there's a bit of numerical noise on the scattering matrix,"
-                  "\n\tand finding them may indicate the calculation is unconverged."
-                  "\n\tWhile we simply do not include these when computing transport, "
-                  "\n\tyou may want to consider using with more wavevectors or an improved DFT calculation."
-                  "\n\tAdditionally, setting symmetrizeMatrix = true in your input file will help.");
+                "\n\tThis can happen when there's a numerical issue with the scattering matrix."
+                "\n\tWhile we simply do not include these when computing transport, "
+                "\n\tand likely if they are very small the calculation will be unaffected, "
+                "\n\tplease view the relaxons tutorial for advice about resolving this.");
         std::cout << "These eigenvalues are (in atomic units):" << std::endl;
 
         for (int i = 0; i < m; i++) {

src/apps/coupled_bte_app.cpp

@@ -169,7 +169,9 @@ void CoupledTransportApp::checkRequirements(Context &context) {
   throwErrorIfUnset(context.getTemperatures(), "temperatures");
   throwErrorIfUnset(context.getSmearingMethod(), "smearingMethod");
   if (context.getSmearingMethod() == DeltaFunction::gaussian) {
-    throwErrorIfUnset(context.getSmearingWidth(), "smearingWidth");
+    if (std::isnan(context.getElSmearingWidth()) || !std::isnan(context.getPhSmearingWidth())) {
+      throwErrorIfUnset(context.getSmearingWidth(), "smearingWidth");
+    }
   }
   if (!context.getElphFileName().empty()) {
     throwErrorIfUnset(context.getElectronH0Name(), "electronH0Name");

src/bte/coupled_scattering_matrix.cpp

@@ -11,7 +11,7 @@
 #include <map>
 /* 
 enum DragType {
-  elph,
+  elph, 
   phel
 };
  */
@@ -191,12 +191,12 @@ void CoupledScatteringMatrix::builder(std::shared_ptr<VectorBTE> linewidth,
   if (!std::isnan(context.getBoundaryLength())) {
     if (context.getBoundaryLength() > 0.) {
       // phonon boundary scattering
-      //addBoundaryScattering(*this, context, inPopulations, outPopulations,
-      //                      innerBandStructure, linewidth);
-      //std::cout << std::endl;
+      addBoundaryScattering(*this, context, inPopulations, outPopulations,
+                            innerBandStructure, linewidth);
+      std::cout << std::endl;
       // electron boundary scattering
-      //addBoundaryScattering(*this, context, inPopulations, outPopulations,
-      //                      outerBandStructure, linewidth);
+      addBoundaryScattering(*this, context, inPopulations, outPopulations,
+                            outerBandStructure, linewidth);
     }
   }
 
@@ -620,7 +620,7 @@ std::vector<std::vector<std::tuple<std::vector<int>, int>>>
 }
 
 // reweight the matrix quadrants
-void CoupledScatteringMatrix::reweightQuadrants() {
+/* void CoupledScatteringMatrix::reweightQuadrants() {
 
   // TODO if we use linewidths also apply 2 to them
 
@@ -668,7 +668,7 @@ void CoupledScatteringMatrix::reweightQuadrants() {
     }
   }
 }
-
+ */
 BaseBandStructure* CoupledScatteringMatrix::getPhBandStructure() { return &innerBandStructure; }
 BaseBandStructure* CoupledScatteringMatrix::getElBandStructure() { return &outerBandStructure; }
 

src/bte/coupled_scattering_matrix.h

@@ -104,7 +104,7 @@ class CoupledScatteringMatrix : virtual public BaseElScatteringMatrix,
 
   /** Function to reweight the different quadrants of the matrix
    * coupled matrix to account for spin dengeneracy */
-  void reweightQuadrants();
+  //void reweightQuadrants();
 
   // friend functions for adding scattering rates
   // see respective header files for more details

src/bte/el_scattering_matrix.cpp

@@ -49,9 +49,8 @@ void ElScatteringMatrix::builder(std::shared_ptr<VectorBTE> linewidth,
                                  getIteratorWavevectorPairs(rowMajor);
 
   // add scattering contributions ---------------------------------------
+  
   // add elph scattering
-  // TODO are we sure this should get two Fermi's and not have one of them be a Bose?
-
   { // let the interaction elph go out of scope after this, it takes a lot of memory
 
   // load the elph coupling
@@ -61,6 +60,7 @@ void ElScatteringMatrix::builder(std::shared_ptr<VectorBTE> linewidth,
   InteractionElPhWan couplingElPh =
       InteractionElPhWan::parse(context, crystal, phononH0);
 
+  if(mpi->mpiHead()) std::cout << std::endl;
   addElPhScattering(*this, context, inPopulations, outPopulations, 
                                   kPairIterator, innerFermi, //outerFermi,
                                   innerBandStructure, outerBandStructure, phononH0,

src/bte/ph_scattering_matrix.cpp

@@ -51,6 +51,8 @@ void PhScatteringMatrix::builder(std::shared_ptr<VectorBTE> linewidth,
 
   Crystal crystal = innerBandStructure.getPoints().getCrystal();
 
+  bool usePhElScattering = !context.getElphFileName().empty(); 
+
   // here we call the function to add ph-ph scattering
   if(!context.getPhFC3FileName().empty()) {
     // read this in and let it go out of scope afterwards
@@ -158,7 +160,14 @@ void PhScatteringMatrix::builder(std::shared_ptr<VectorBTE> linewidth,
 
   // recalculate the phonon linewidths from the off diagonals
   // we should do this if phel is not involved, otherwise it wipes out phel
-  //enforceDetailedBalance();
+  if(context.getEnforceDetailedBalance()) {
+    if(usePhElScattering) {
+      Warning("Ignoring request to enforce detailed balance on the scattering matrix,"
+        "as this will erase the ph-el contribution along the diagonal!"); 
+    } else {
+      enforceDetailedBalance();
+    }
+  }
 
   // some phonons like acoustic modes at the gamma, with omega = 0,
   // might have zero frequencies, and infinite populations. We set those
@@ -246,5 +255,102 @@ void PhScatteringMatrix::builder(std::shared_ptr<VectorBTE> linewidth,
   }
 }
 
+void PhScatteringMatrix::enforceDetailedBalance() {
+
+  // kill the function if it's used inappropriately
+  if(isMatrixOmega) { // If this matrix has not been symmetrized, this function won't work
+    DeveloperError("enforceDetailedBalance should not be called on a matrix without symmetrization factors in the ph only case.");
+  }
+  if(!highMemory) return;  // must be high mem, we explicitly use iCalc=1 here
+  if(context.getUseSymmetries()) return; // this is not designed for BTE syms, would need to
+                                         // change the way we are indexing this
+  if(context.getUseUpperTriangle()) {    // TODO this can be implemented without too much difficulty
+    Warning("Cannot run enforceDetailedBalance with only upper triangle for now.");
+    return;
+  };
+
+  if(mpi->mpiHead()) std::cout << "\nEnforcing detailed balance, summing off-diagonals to compute diagonal." << std::endl;
+
+  // this only happens when we have one iCalc, and no symmetries -- VectorBTE = an array,
+  // which we use here because of some trouble with coupledVectorBTE object
+  Eigen::MatrixXd newLinewidths(1,numStates); // copy matrix which is the same as internal diag
+  newLinewidths.setZero();
+
+  // rebuild the diagonal ---------------------------------
+
+  LoopPrint loopPrint("Recalculating the diagonal","matrix elements",getAllLocalStates().size());
+
+  // NOTE: if later we want to use symmetries here,
+  // these would actually be iBTE instead of iState, and we would convert
+  // sum over the v' states owned by this process
+  // TODO may want to add OMP here as well as MPI 
+  for (auto [ibte1, ibte2] : getAllLocalStates()) {
+
+    loopPrint.update();
+
+    // throw out lower triangle states
+    // FIXME DOESNT WORK!
+    if(context.getUseUpperTriangle() && ibte1 > ibte2) continue;
+
+    // we are computing the diagonal using the off diagonal,
+    // so these elements won't matter
+    if(ibte1 == ibte2) continue;
+
+    // NOTE state and bte indices are here the same, as we are blocking the use of symmetries
+    //BteIndex bteIdx1(ibte1);
+    //BteIndex bteIdx2(ibte2);
+    StateIndex is1(ibte1);
+    StateIndex is2(ibte2);
+
+    double initialEn = innerBandStructure.getEnergy(is1);
+    double finalEn = outerBandStructure.getEnergy(is2);
+
+    // remove accoustic phonons
+    if(initialEn < 1e-9 || finalEn < 1e-9) continue;
+
+    // calculate the new linewidths
+    newLinewidths(0,ibte1) -= theMatrix(ibte1,ibte2) * finalEn / initialEn;
+  }
+  loopPrint.close();
+
+  // sum the element contributions from all processes
+  mpi->allReduceSum(&newLinewidths);
+
+  if(mpi->mpiHead()) {
+    std::cout << "\nChecking the quality of recalculated ph diagonal elements." << std::endl;
+
+    for (int i = 0; i<numStates; i++) {
+
+      // don't print zeros
+      if(newLinewidths(0,i) < 1e-15 && internalDiagonal->data(0,i) < 1e-15) continue;
+
+      if(newLinewidths(0,i) < 0 || std::isnan(newLinewidths(0,i))) {
+        StateIndex sIdx(i);
+        std::cout << std::setprecision(4) << "Found a negative ph linewidth for state: " << i << " " << innerBandStructure.getEnergy(sIdx) << " " << innerBandStructure.getPoints().cartesianToCrystal(innerBandStructure.getWavevector(sIdx)).transpose() << " " << internalDiagonal->data(0,i) << " " << newLinewidths(0,i) << std::endl;
+        // replace with the standard one to avoid definite issues
+        newLinewidths(0,i) = internalDiagonal->data(0, i);
+      }
+      // flag bad linewidth ratios
+      else if(newLinewidths(0,i)/internalDiagonal->data(0,i) < 0.25 || newLinewidths(0,i)/internalDiagonal->data(0,i) > 1.75) {
+        StateIndex sIdx(i);
+        std::cout << std::setprecision(4) << "Found a bad ph linewidth for state: " << i << " " << innerBandStructure.getEnergy(sIdx) << " " << innerBandStructure.getPoints().cartesianToCrystal(innerBandStructure.getWavevector(sIdx)).transpose() << " " << newLinewidths(0,i) << " " << internalDiagonal->data(0,i)  << " " << newLinewidths(0,i)/internalDiagonal->data(0,i) << std::endl;
+        newLinewidths(0,i) = std::max(newLinewidths(0,i),internalDiagonal->data(0,i));
+      }
+    }
+  }
+
+  // reinsert the linewidths in the scattering matrix
+  internalDiagonal->data = newLinewidths;
+
+  // TODO figure out why this function doesn't work right now
+  //replaceMatrixLinewidths();
+
+  // replace the linewidths on the scattering matrix diagonal
+  int iCalc = 0;
+  for (int iMat = 0; iMat < numStates; iMat++) {
+    // zero the diagonal of the matrix
+    if(theMatrix.indicesAreLocal(iMat,iMat)) theMatrix(iMat, iMat) = newLinewidths(iCalc, iMat);
+  }
+}
 
 

src/bte/ph_scattering_matrix.h

@@ -46,6 +46,8 @@ class PhScatteringMatrix : virtual public BasePhScatteringMatrix {
                std::vector<VectorBTE> &inPopulations,
                std::vector<VectorBTE> &outPopulations) override;
 
+  void enforceDetailedBalance(); 
+
  // friend functions for adding scattering rates,
  // these live in ph_scattering.cpp
  // TODO write docstrings for these

src/bte/scattering_matrix.cpp

@@ -1470,6 +1470,7 @@ void ScatteringMatrix::enforceDetailedBalance() {
   // NOTE: if later we want to use symmetries here,
   // these would actually be iBTE instead of iState, and we would convert
   // sum over the v' states owned by this process
+  // TODO may want to add OMP here as well as MPI 
   for (auto [ibte1, ibte2] : getAllLocalStates()) {
 
     loopPrint.update();
@@ -1501,10 +1502,6 @@ void ScatteringMatrix::enforceDetailedBalance() {
     Particle initialParticle = initialBandStructure->getParticle();
     Particle finalParticle = finalBandStructure->getParticle();
 
-    // this removes the drag term contributions, for test reasons
-    //if(initialParticle.isPhonon() && finalParticle.isElectron()) continue;
-    //if(initialParticle.isElectron() && finalParticle.isPhonon()) continue;
-
     // shift the indices back to the ones used in bandstructures
     // these indices are for the full matrix, if the matrix is coupled,
     // we need to fold them back into the relevants quadrants in order
@@ -1531,7 +1528,7 @@ void ScatteringMatrix::enforceDetailedBalance() {
     double initialFFm1 = initialParticle.getPopPopPm1(initialEn, kBT, initialChemicalPotential);
     double finalFFm1 = finalParticle.getPopPopPm1(finalEn, kBT, finalChemicalPotential);
 
-    // spin degeneracy info -- TODO may need to put spin factors here
+    // spin degeneracy info -- TODO may need to put spin factors here?
     double initialD = 1;
     double finalD = 1;
 
@@ -1607,7 +1604,7 @@ void ScatteringMatrix::enforceDetailedBalance() {
       }
     }
   }
-
+/* 
   if(mpi->mpiHead()) {
 
     std::cout << "compare first 50 el states, new vs. old " << std::setw(2) << std::scientific << std::setprecision(2) << std::endl;
@@ -1637,7 +1634,7 @@ void ScatteringMatrix::enforceDetailedBalance() {
 
     }
     std::cout << std::scientific << std::setprecision(4) << std::endl;
-  }
+  } */
   // reinsert the linewidths in the scattering matrix
   internalDiagonal->data = newLinewidths;
 

src/context.cpp

@@ -1022,7 +1022,6 @@ void Context::printInputSummary(const std::string &fileName) {
     // specific to coupled scattering matrix app
     if (appName.find("oupled") != std::string::npos) {
       std::cout << "useDragTerms = " << useDragTerms << std::endl;
-      std::cout << "enforceDetailedBalance = " << enforceDetailedBalance << std::endl;
     }
 
     if (!std::isnan(constantRelaxationTime))
@@ -1079,7 +1078,8 @@ void Context::printInputSummary(const std::string &fileName) {
         if(numRelaxonsEigenvalues != 0) std::cout << "checkNegativeRelaxons = " << checkNegativeRelaxons << std::endl;
       }
       if(scatteringMatrixInMemory) {
-        std::cout << "enforcePositiveSemiDefinite = " << enforcePositiveSemiDefinite << std::endl;
+        std::cout << "enforceDetailedBalance = " << enforceDetailedBalance << std::endl;
+      //  std::cout << "enforcePositiveSemiDefinite = " << enforcePositiveSemiDefinite << std::endl;
       }
     }
 

src/observable/viscosity_io.cpp

@@ -36,6 +36,7 @@ int relaxonEigenvectorOverlap(ParallelMatrix<double>& eigenvectors,
     std::cout << "\nMaximum scalar product " << eigenvectorName << ".theta_alpha = " << maxOverlap << " at alpha = " << idxMaxOverlap << "." << std::endl;
     std::cout << "First ten products with " << eigenvectorName << ":";
     for(int gamma = 0; gamma < maxPrint; gamma++) { std::cout << " " << overlaps(gamma); }
+    std::cout << std::endl;
   }
 
   // If the best overlap isn't very good, we return -1 so nothing is skipped