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> class fdmBackwardIteration {
public:
  nrgcore_type *nrgObject;
  double        kBT{0};

  explicit fdmBackwardIteration(nrgcore_type *t_nrgObject) {
    setup(t_nrgObject);
    this->clearKeptIndex();
  }

  void setup(nrgcore_type *t_nrgObject) {
    lastiteration = true;
    nrgObject     = t_nrgObject;
  }

  void calcSpectrum([[maybe_unused]] double energyScale) {
    // Clear the operator
    setCurrentIndex();
    // Order of these functions are important
    setRhoZero();
    // 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.
    }
    // once the reduced density matrix is defined we
    // set the last iteration flag for the next shell.
    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();
        auto sys_opr_opt = (*bOperator)[ip].get(i, i);
        if (sys_opr_opt) {
          auto *sys_opr = sys_opr_opt.value();
          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);
            }
          }
        }
      }
    }
    return specSum;
  }

  void setRhoZero(const std::vector<std::vector<double>> &tBoltzmannFactor) {
    rhoZero.clear();
    double rhoTrace = 0;
    for (size_t i = 0; i < nrgObject->current_sysmQ.size(); i++) {
      size_t kpdim = nrgObject->eigenvaluesQ[i].size();
      qmatrix<> tmat(kpdim, kpdim, 0);
      for (size_t ie = currentKeptIndex[i].size();
           ie < nrgObject->eigenvaluesQ[i].size(); ie++) {
        tmat(ie, ie) = tBoltzmannFactor[i][ie];
      }
      if (!lastiteration) {
        for (auto ik : currentKeptIndex[i]) {
          for (auto ikp : currentKeptIndex[i]) {
            tmat(ik, ikp) = reducedRho[i](ik, ikp);
          }
        }
      }
      rhoTrace += tmat.trace();
      rhoZero.push_back(tmat);
    }
    std::cout << "NRG Itr: " << nrgObject->nrg_iterations_cnt
              << "rhoTrace: " << rhoTrace << std::endl;
  }

  void setRhoZero() {
    if (lastiteration) {
      setLocalPartitionFunction();
    }
    double rhoTrace{0};
    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();
      qmatrix<> tmat(kpdim, kpdim, 0);
      if (lastiteration) {
        for (size_t ie = 0; ie < nrgObject->eigenvaluesQ[i].size(); ie++) {
          tmat(ie, ie) = BoltzmannFactor[i][ie];
        }
      }
      if (!lastiteration) {
        for (auto ik : currentKeptIndex[i]) {
          for (auto ikp : currentKeptIndex[i]) {
            tmat(ik, ikp) = reducedRho[i](ik, ikp);
          }
        }
      }
      rhoTrace += tmat.trace();
      rhoZero.push_back(tmat);
    }
    std::cout << "NrgItr: " << nrgObject->nrg_iterations_cnt
              << "rhoTrace: " << rhoTrace << std::endl;
  }

  void setTemperature(double mkBT) { kBT = mkBT; }

  void setCurrentIndex() {
    if (lastiteration) {
      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};
};