scattering_matrix.h

#ifndef SCATTERING_H
#define SCATTERING_H
 
#include "Matrix.h"
#include "context.h"
#include "delta_function.h"
#include "vector_bte.h"
 
// 3 cases:
// we compute and store in memory the scattering matrix and the diagonal
// we compute the action of the scattering matrix on the in vector, returning outVec = sMatrix*vector
// we compute only the linewidths
enum MatrixCase {
  fullMatrix,
  matrixVectorProduct,
  linewidthOnly
};
 
class ScatteringMatrix {
public:
  ScatteringMatrix(Context &context_, StatisticsSweep &statisticsSweep_,
                   BaseBandStructure &innerBandStructure_,
                   BaseBandStructure &outerBandStructure_);
 
  ScatteringMatrix();
 
  void setup();
 
  VectorBTE diagonal();
 
  double& operator()(const int &row, const int &col);
 
  VectorBTE offDiagonalDot(VectorBTE &inPopulation);
  std::vector<VectorBTE> offDiagonalDot(std::vector<VectorBTE> &inPopulations);
 
  VectorBTE dot(VectorBTE &inPopulation);
  std::vector<VectorBTE> dot(std::vector<VectorBTE> &inPopulations);
 
//  /** Computes the product A*B, where A is the scattering matrix, and
//   * B is an Eigen::MatrixXd. This can be used to compute products of the
//   * scattering matrix with other vectors.
//   */
//  ParallelMatrix<double> dot(const ParallelMatrix<double> &otherMatrix);
 
  VectorBTE getSingleModeTimes(); 
  
  VectorBTE getLinewidths();
 
  void setLinewidths(VectorBTE &linewidths);
  
   VectorBTE getSingleModeTimes(const VectorBTE& anyInternalDiagonal);
 
   VectorBTE getLinewidths(const VectorBTE& anyInternalDiagonal);
 
  void a2Omega();
 
  void omega2A();
 
  std::tuple<Eigen::VectorXd, ParallelMatrix<double>>
                                diagonalize(int numEigenvalues = 0);
 
  int getSMatrixIndex(BteIndex &bteIndex, CartIndex &cartIndex);
 
  std::tuple<BteIndex, CartIndex> getSMatrixIndex(const int &iMat);
 
  void symmetrize();
 
  void relaxonsToJSON(const std::string& fileName, const Eigen::VectorXd& eigenvalues);
 
  static void symmetrizeCoupling(Eigen::Tensor<double,3>& coupling,
                        const Eigen::VectorXd& energies1,
                        const Eigen::VectorXd& energies2,
                        const Eigen::VectorXd& energies3);
 
  std::vector<std::tuple<int, int>> getAllLocalStates();
 
  // IO operations ---------
 
  void outputToHDF5(const std::string &outFileName);
 
 protected:
 
  Context &context;
  StatisticsSweep &statisticsSweep;
 
  // Smearing is a pointer created in constructor with a smearing factory
  // Used in the construction of the scattering matrix to approximate the
  // Dirac-delta function in transition rates.
  //DeltaFunction *smearing; // moved to specific ph and el base scatteringMatrix objs
 
  BaseBandStructure &innerBandStructure;
  BaseBandStructure &outerBandStructure;
  
  enum MatrixCase matrixCase;
  
  // constant relaxation time approximation -> the matrix is just a scalar
  // and there are simplified evaluations taking place
  bool constantRTA = false;
  bool highMemory = true;     // whether the matrix is kept in memory
  bool outputUNTimes = false;    // whether to output U and N processes in RTA
 
  double boundaryLength;
  bool doBoundary;
 
  // NOTE: about the definition of A
  // A acts on the canonical population, omega on the population
  // A acts on f, Omega on n, with n = bose(bose+1)f e.g. for phonons
  //
  // Generally: A_out = diagonal of the scattering matrix.
  // In the case of phonons, A_out = linewidth * n(n+1)
  // In the case of electrons, A_out = linewidht
  bool isMatrixOmega = false; // whether the matrix is Omega or A
 
  // if we're using the coupled matrix, we need to use this to
  // save the scattering rates differently
  bool isCoupled = false;
 
  /* function which shifts the indices to the correct quadrant in the CBTE case. 
  * for a pure ph or el matrix, does nothing (default defined here)
  * @param iBte1, iBte2: the bte indices used to index this quadrant.
  * @param p1, p2: the particle types indicating the desired quadrant
  * @return: a tuple containing the shifted indices
  */
  std::function<std::tuple<long,long>(long, long, const Particle&, const Particle &)> shiftToCoupledIndices 
      = [](long iBte1, long iBte2, [[maybe_unused]] const Particle &p1, [[maybe_unused]] const Particle &p2){ 
          return std::make_tuple(iBte1, iBte2); }; 
 
  // we save the diagonal matrix element in a dedicated vector
  std::shared_ptr<VectorBTE> internalDiagonal, internalDiagonalUmklapp, internalDiagonalNormal;
  // the scattering matrix, initialized if highMemory==true
  ParallelMatrix<double> theMatrix;
 
  int numStates; // number of Bloch states (i.e. the size of theMatrix)
  int numCalculations;  // number of "Calculations", i.e. number of temperatures and
  // chemical potentials on which we compute scattering.
  int dimensionality_;
 
  // number of states for shifting to the coupled matrix
  int numElStates = 0; // this will never be used unless it's a coupled matrix
  // however, it still needs to be here because the coupled
  // shift incides function unfortunately also has to be here. For the case of the
  // ph matrix, we call phScattering using the BasePh object -- even if the
  // matrix supplied is a Coupled matrix, the base class's function is called.
 
  // this is to exclude problematic Bloch states, e.g. acoustic phonon modes
  // at gamma, which have zero frequencies and thus have several non-analytic
  // behaviors (e.g. Bose--Einstein population=infinity). We set to zero
  // terms related to these states.
  std::vector<int> excludeIndices;
 
  virtual void builder(std::shared_ptr<VectorBTE> linewidth,
                       std::vector<VectorBTE> &inPopulations,
                       std::vector<VectorBTE> &outPopulations) = 0;
 
  std::vector<std::tuple<std::vector<int>, int>> getIteratorWavevectorPairs(
                                                const bool &rowMajor = false);
 
  void degeneracyAveragingLinewidths(std::shared_ptr<VectorBTE> linewidth);
 
  Eigen::MatrixXd precomputeOccupations(BaseBandStructure &bandStructure);
 
  std::vector<int> getExcludeIndices(BaseBandStructure& bandStructure);
 
  /* Set the use case of the scattering matrix into a class enum. 
  * Choices are 
  * we compute and store in memory the scattering matrix and the diagonal
  * we compute the action of the scattering matrix on the in vector, returning outVec = sMatrix*vector
  * we compute only the linewidths
  * @param linewidths pointer to a vectorBTE containing the internal diagonal of the scattering matrix 
  * @param inPopulations particle population the scattering matrix acts on 
  * @param outPopulations particle population Smatrix * inPopulation 
  */
  void setMatrixCase(std::shared_ptr<VectorBTE> linewidth, 
    std::vector<VectorBTE> &inPopulations,
    std::vector<VectorBTE> &outPopulations); 
 
  void addRateToMatrix(const Context &context, double linewidthRate, double matrixRate, 
                                      int iCalc, int is1, int is2Irr, int iBte1, int iBte2, 
                                      const Particle &p1, const Particle &p2, 
                                      const Eigen::Matrix3d &rotation, 
                                      std::shared_ptr<VectorBTE> linewidth, 
                                      const std::vector<VectorBTE> &inPopulations,
                                      std::vector<VectorBTE> &outPopulations); 
 
  void enforceDetailedBalance();
 
  void replaceMatrixLinewidths();
 
  std::tuple<BaseBandStructure*,BaseBandStructure*> getStateBandStructures(const BteIndex& iBte1,
                                                                    const BteIndex& iBte2);
 
  std::tuple<int,int> coupledToBandStructureIndices(const int& iBte1, const int& iBte2,
                                                    const Particle& p1, const Particle& p2);
 
  template <size_t n>
  void addUNRates(int iCalc, int iBte, double rate, 
                                    const std::array<Point, n> &crysPoints, auto momentumConsExpr) {
 
    if(outputUNTimes) {
      // Note : for the purpose of folding bzToWs for points, 
      // the bandstructure we use doesn't matter, only the points object
 
      // get vectors in ws cell 
      std::array<Eigen::Vector3d, n> wsVectors;  
      for(size_t ipt = 0; ipt < n; ipt++) {
        Eigen::Vector3d wvCart = crysPoints[ipt].getCoordinates(Points::cartesianCoordinates); 
        Eigen::Vector3d wvWS = outerBandStructure.getPoints().bzToWs(wvCart, Points::cartesianCoordinates); 
        wsVectors[ipt] = wvWS; 
      }
      
      // check if this process is umklapp
      Eigen::Vector3d wvFinalCart = std::apply(momentumConsExpr, wsVectors); 
      Eigen::Vector3d wvFinalFold = outerBandStructure.getPoints().bzToWs(wvFinalCart, Points::cartesianCoordinates);
      if(abs((wvFinalCart-wvFinalFold).norm()) > 1e-6) {
        internalDiagonalUmklapp->operator()(iCalc, 0, iBte) += rate;
      } else {
        internalDiagonalNormal->operator()(iCalc, 0, iBte) += rate;
      }
    }
  }
  
  // friend functions for scattering
  friend void addBoundaryScattering(ScatteringMatrix &matrix, Context &context,
    std::vector<VectorBTE> &inPopulations,
    std::vector<VectorBTE> &outPopulations,
    BaseBandStructure &outerBandStructure,
    std::shared_ptr<VectorBTE> linewidth);
 
};
 
#endif