Program Listing for File fdmBackwardIteration.hpp

↰ Return to documentation for file (dynamics/fdmBackwardIteration.hpp)

#pragma once
#include "nrgcore/nrgData.hpp"
#include "nrgcore/qOperator.hpp"
#include "utils/qmatrix.hpp"
#include "utils/timer.hpp"
#include <algorithm> // std::min_element
#include <cmath>
#include <iostream>
#include <iterator>
#include <map>
#include <numeric>
#include <optional>
#include <tuple>
#include <vector>
template <typename nrgcore_type> // nrgcore_type is a type of
class fdmBackwardIteration {
public:
  nrgcore_type *nrgObject;
  double        kBT{0}; // Temperature of nrg system i.e., in FDM formalism
  // fdmBackwardIteration() {}
  explicit fdmBackwardIteration(nrgcore_type *t_nrgObject // nrgcore_type
  ) {
    setup(t_nrgObject);
    this->clearKeptIndex();
  }
  void setup(nrgcore_type *t_nrgObject) {
    lastiteration = true;
    nrgObject     = t_nrgObject;
  }
  void calcSpectrum(double energyScale) {
    // Clear the operator
    setCurrentIndex();
    // Order of these functions are important
    setRhoZero(energyScale);
    // rhoDotOperators();
    setReduceDensityMatrix();
  }
  void setReduceDensityMatrix() {
    // Set reducedRho
    reducedRho.clear();
    for (size_t i = 0; i < nrgObject->pre_sysmQ.size(); i++) {
      reducedRho.push_back( //
          qmatrix<>(nrgObject->eigenvaluesQ_kept_indices[i].size(),
                    nrgObject->eigenvaluesQ_kept_indices[i].size(), 0));
    }
    // Rotate the eigen basis
    for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
      // U. rhoZero . U.T : TODO: check
      rhoZero[i] = nrgObject->current_hamiltonQ[i].dot(
          rhoZero[i].dot(nrgObject->current_hamiltonQ[i].cTranspose()));
    }
    // Set reducedRho
    for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
      size_t kidx = 0;
      for (auto kindex : nrgObject->coupled_nQ_index[i]) {
        auto ii = kindex / nrgObject->nq_bath.size(); // impurity nqi index
        auto bb = kindex % nrgObject->nq_bath.size(); // bath nqi index
        // create previous bath id matrix
        for (size_t it : nrgObject->eigenvaluesQ_kept_indices[ii]) {
          for (size_t it_p : nrgObject->eigenvaluesQ_kept_indices[ii]) {
            double aa{0};
            for (size_t il = 0; il < nrgObject->bath_eigenvaluesQ[bb].size();
                 il++) {
              aa += rhoZero[i].at(
                  kidx + it +
                      nrgObject->eigenvaluesQ_kept_indices[ii].size() * il,
                  kidx + it_p +
                      nrgObject->eigenvaluesQ_kept_indices[ii].size() * il);
            }
            reducedRho[ii].at(it, it_p) += aa;
          }
        }
        kidx += nrgObject->eigenvaluesQ_kept_indices[ii].size() *
                nrgObject->bath_eigenvaluesQ[bb].size();
      }
      // End of matrix generation.
    }
    // double trace = std::accumulate(
    //                    reducedRho.begin(), reducedRho.end(), 0.0,
    // [](double a, const qmatrix<> &b) {
    //     return a + b.trace();
    // });
    // std::cout << "NrgItr: " << nrgObject->nrg_iterations_cnt
    //           << "rhoTrace: " << trace << std::endl;
    // once the reduce density matrix is defined We
    // set the lastiteration to be false for the next iteration
    lastiteration = false;
  }
  void setLocalPartitionFunction() {
    localGroundStateEnergy = 0;
    localPartitionFunction = 0; // Ground state degenarecy
    for (const auto &aa : nrgObject->eigenvaluesQ) {
      if (aa.size() != 0) {
        double result          = *std::min_element(aa.begin(), aa.end());
        localGroundStateEnergy = std::min(localGroundStateEnergy, result);
      }
    }
    for (size_t i = 0; i < nrgObject->eigenvaluesQ.size(); i++) {
      for (size_t ie = 0; ie < nrgObject->eigenvaluesQ[i].size(); ie++) {
        double energy =
            std::fabs(nrgObject->eigenvaluesQ[i][ie] - localGroundStateEnergy);
        if (energy < energyErrorBar) {
          localPartitionFunction += 1.;
        }
      }
    } //
    std::cout << "localGroundStateEnergy" << localGroundStateEnergy
              << " localPartitionFunction: " << localPartitionFunction
              << std::endl;
    BoltzmannFactor = nrgObject->eigenvaluesQ;
    for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
      for (size_t ie = 0; ie < nrgObject->eigenvaluesQ[i].size(); ie++) {
        double energy =
            std::fabs(nrgObject->eigenvaluesQ[i][ie] - localGroundStateEnergy);
        if (energy < energyErrorBar) {
          BoltzmannFactor[i][ie] = 1. / localPartitionFunction;
        } else {
          BoltzmannFactor[i][ie] = 0;
        }
      }
    } //
  }
  auto rhoDotStaticOperators(std::vector<qOperator> *bOperator) {
    // timer               t1("rhoDotStaticOperators");
    std::vector<double> specSum(bOperator->size(), 0.0);
    //
    for (size_t ip = 0; ip < bOperator->size(); ip++) {
      for (size_t i = 0; i < nrgObject->eigenvaluesQ.size(); i++) {
        size_t kpdim = nrgObject->eigenvaluesQ[i].size();
        // std::cout << "--------------------------";
        // std::cout << "idx" << idx << " idx_p" << idx_p << std::endl;
        // TODO(sp): This
        auto sys_opr_opt = (*bOperator)[ip].get(i, i);
        if (sys_opr_opt) {
          auto *sys_opr = sys_opr_opt.value();
          // set kept-kept part operator
          for (auto iv : currentKeptIndex[i]) {
            for (auto iv_p : currentKeptIndex[i]) {
              sys_opr->at(iv, iv_p) = 0;
            }
          }
          for (size_t iv = 0; iv < kpdim; iv++) {
            for (size_t iv_p = 0; iv_p < kpdim; iv_p++) {
              specSum[ip] += sys_opr->at(iv, iv_p) * rhoZero[i].at(iv_p, iv);
            }
          }
        }
        // all the matrices are set
      }
    }
    // set the matrix elements
    // end of lm loop
    // End of matrix generation.
    // Rotate the c operator in the eigen basis
    return specSum;
  }
  void setRhoZero(const std::vector<std::vector<double>> &tBoltzmannFactor) {
    // Clear the operator
    // Order of these functions are important
    rhoZero.clear();
    double rhoTrace = 0;
    for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
      size_t kpdim = nrgObject->eigenvaluesQ[i].size();
      // std::cout << "--------------------------";
      // std::cout << "idx" << idx << " idx_p" << idx_p << std::endl;
      qmatrix<> tmat(kpdim, kpdim, 0);
      // Discarded states
      for (size_t ie = currentKeptIndex[i].size();
           ie < nrgObject->eigenvaluesQ[i].size(); ie++) {
        tmat(ie, ie) = tBoltzmannFactor[i][ie];
      }
      // Kept states
      if (!lastiteration) { // Condition for the last Wilson site
        // std::cout << "Not lastiteration" << std::endl;
        //  Just Override the kept states
        for (auto ik : currentKeptIndex[i]) {
          for (auto ikp : currentKeptIndex[i]) {
            tmat(ik, ikp) = reducedRho[i](ik, ikp);
          }
        }
      }
      rhoTrace += tmat.trace();
      // Save the matrix
      rhoZero.push_back(tmat);
    }
    std::cout << "NRG Itr: " << nrgObject->nrg_iterations_cnt
              << "rhoTrace: " << rhoTrace << std::endl;
    // move the operator
  } // End of  update_system_operatorQ
  void setRhoZero() {
    // Clear the operator
    // Order of these functions are important
    // This function is called
    // for setting up current rho
    //
    //
    //
    if (lastiteration) {
      setLocalPartitionFunction();
    }
    double rhoTrace{0};
    // Calc partition function of the shell ::
    vecPartitions.push_back(localPartitionFunction);
    rhoZero.clear();
    std::cout << "Size : " << reducedRho.size() << " "
              << currentKeptIndex.size() << std::endl;
    for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
      size_t kpdim = nrgObject->eigenvaluesQ[i].size();
      // std::cout << "--------------------------";
      // std::cout << "idx" << idx << " idx_p" << idx_p << std::endl;
      qmatrix<> tmat(kpdim, kpdim, 0);
      // Discarded states
      if (lastiteration) { // Only for the laast iterationT= 0
        for (size_t ie = 0; ie < nrgObject->eigenvaluesQ[i].size(); ie++) {
          tmat(ie, ie) = BoltzmannFactor[i][ie];
        }
      }
      // Kept states
      if (!lastiteration) { // Condition for the last Wilson site
        // Just Override the kept states
        // std::cout << reducedRho[i].size() << " " <<
        // currentKeptIndex[i].size()
        //           << std::endl;
        for (auto ik : currentKeptIndex[i]) {
          for (auto ikp : currentKeptIndex[i]) {
            tmat(ik, ikp) = reducedRho[i](ik, ikp);
          }
        }
      }
      rhoTrace += tmat.trace();
      // Save the matrix
      rhoZero.push_back(tmat);
    }
    std::cout << "NrgItr: " << nrgObject->nrg_iterations_cnt
              << "rhoTrace: " << rhoTrace << std::endl;
    // move the operator
  } // End of  update_system_operatorQ
  //
  void setTemperature(double mkBT) { kBT = mkBT; }
  void setCurrentIndex() {
    // This only work for the backward iteration
    if (lastiteration) {
      // every state is Discarded
      for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
        currentKeptIndex.emplace_back();
      }
    } else {
      currentKeptIndex = previoudKeptIndex;
    }
    previoudKeptIndex = nrgObject->eigenvaluesQ_kept_indices;
  }
  void clearKeptIndex() {
    lastiteration = true;
    BoltzmannFactor.clear();
    vecPartitions.clear();
    rhoZero.clear();
    reducedRho.clear();
    currentKeptIndex.clear();
    previoudKeptIndex.clear();
  }

private:
  std::vector<std::vector<double>> BoltzmannFactor;
  double                           localGroundStateEnergy{0};
  double localPartitionFunction{0}; // Ground state degenarecy
  std::vector<std::vector<size_t>> previoudKeptIndex;
  double spWeightErrorBar{1e-20}; // We dont care for the lower value
  std::vector<qmatrix<>> reducedRho;
  std::vector<double>    vecPartitions;

public: // Give access for openchain class
  std::vector<std::vector<size_t>> currentKeptIndex;
  std::vector<qmatrix<>>           rhoZero;
  bool                             lastiteration{true};
  double                           energyErrorBar{1e-5};
};