Program Listing for File fdmSpectrum.hpp#
↰ Return to documentation for file (dynamics/fdmSpectrum.hpp)
#pragma once
#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
// We don't want to inherit the nrgcore_type.
// Just want to use few of its objects as a pointer
// We don't want to reference to any specific model here.
//
// Warning: this class DOES NOT take care of the fermion sign.
// You need to provide the operator with sign
// TODO(sp): : kBt != 0 and a = b^\dag
class fdmSpectrum {
nrgcore_type *nrg_object;
double kBT{0}; // Temperature of nrg system i.e., in FDM formalism
public:
// fdmSpectrum() {}
explicit fdmSpectrum(nrgcore_type *t_nrg_object // nrgcore_type
) {
setup(t_nrg_object);
}
void setup(nrgcore_type *t_nrg_object) {
lastiteration = true;
nrg_object = t_nrg_object;
globalGroundStateEnergy = 0; //= nrg_object->all_eigenvalue[0];
}
void calcSpectrum(double energyScale) {
energyRescale = energyScale;
// Clear the operator
setCurrentIndex();
setRhoZero();
rhoDotOperators();
setReduceDensityMatrix();
lastiteration = false;
}
void setReduceDensityMatrix() {
// timer t1("setReduceDensityMatrix");
// Set reducedRho
reducedRho.clear();
for (size_t i = 0; i < nrg_object->pre_sysmQ.size(); i++) {
reducedRho.push_back( //
qmatrix<>(nrg_object->eigenvaluesQ_kept_indices[i].size(),
nrg_object->eigenvaluesQ_kept_indices[i].size(), 0));
}
// Rotate the eigen basis
// #pragma omp parallel for
for (size_t i = 0; i < nrg_object->current_sysmQ.size(); i++) {
// U. rhoZero . U.T : TODO: check
rhoZero[i] = nrg_object->current_hamiltonQ[i].dot(
rhoZero[i].dot(nrg_object->current_hamiltonQ[i].cTranspose()));
}
// Set reducedRho
for (size_t i = 0; i < nrg_object->current_sysmQ.size(); i++) {
size_t kidx = 0;
for (auto kindex : nrg_object->coupled_nQ_index[i]) {
auto ii = kindex / nrg_object->nq_bath.size(); // impurity nqi index
auto bb = kindex % nrg_object->nq_bath.size(); // bath nqi index
// create previous bath id matrix
// TODO(sp): Add parallelization here
for (size_t it : nrg_object->eigenvaluesQ_kept_indices[ii]) {
for (size_t it_p : nrg_object->eigenvaluesQ_kept_indices[ii]) {
double aa{0};
for (size_t il = 0; il < nrg_object->bath_eigenvaluesQ[bb].size();
il++) {
aa += rhoZero[i].at(
kidx + it +
nrg_object->eigenvaluesQ_kept_indices[ii].size() * il,
kidx + it_p +
nrg_object->eigenvaluesQ_kept_indices[ii].size() * il);
}
reducedRho[ii].at(it, it_p) += aa;
}
}
kidx += nrg_object->eigenvaluesQ_kept_indices[ii].size() *
nrg_object->bath_eigenvaluesQ[bb].size();
}
// End of matrix generation.
}
}
void setLocalPartitionFunction() {
localGroundStateEnergy = 0;
localPartitionFunction = 0; // Ground state degenarecy
for (auto &aa : nrg_object->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 < nrg_object->current_sysmQ.size(); i++) {
for (size_t ie = 0; ie < nrg_object->eigenvaluesQ[i].size(); ie++) {
double energy =
std::fabs(nrg_object->eigenvaluesQ[i][ie] - localGroundStateEnergy);
if (energy < energyErrorBar) {
localPartitionFunction += 1.;
}
}
} //
// std::cout << "localGroundStateEnergy" << localGroundStateEnergy
// << " localPartitionFunction: " << localPartitionFunction
// << std::endl;
BoltzmannFactor = nrg_object->eigenvaluesQ;
for (size_t i = 0; i < nrg_object->current_sysmQ.size(); i++) {
for (size_t ie = 0; ie < nrg_object->eigenvaluesQ[i].size(); ie++) {
double energy =
std::fabs(nrg_object->eigenvaluesQ[i][ie] - localGroundStateEnergy);
if (energy < energyErrorBar) {
BoltzmannFactor[i][ie] = 1. / localPartitionFunction;
} else {
BoltzmannFactor[i][ie] = 0;
}
}
} //
}
std::vector<std::vector<double>> BoltzmannFactor;
void setRhoZero() {
if (lastiteration) {
setLocalPartitionFunction();
}
double rhoTrace{0};
// Calc partition function of the shell ::
vecPartitions.push_back(localPartitionFunction);
rhoZero.clear();
for (size_t i = 0; i < nrg_object->current_sysmQ.size(); i++) {
size_t kpdim = nrg_object->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 < nrg_object->eigenvaluesQ[i].size(); ie++) {
tmat(ie, ie) = BoltzmannFactor[i][ie];
}
}
// Kept states
if (!reducedRho.empty()) { // Condition for the last Wilson site
// 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 << "rhoTrace: " << rhoTrace << std::endl;
// move the operator
} // End of update_system_operatorQ
void rhoDotOperators() { // NOLINT
double specSum = 0.0;
//
for (size_t i = 0; i < nrg_object->current_sysmQ.size(); i++) {
for (size_t j = 0; j < nrg_object->current_sysmQ.size(); j++) {
size_t kpdim = nrg_object->eigenvaluesQ[i].size();
size_t kpdim_p = nrg_object->eigenvaluesQ[j].size();
// std::cout << "--------------------------";
// std::cout << "idx" << idx << " idx_p" << idx_p << std::endl;
// TODO(sp): This
for (size_t ip = 0; ip < bOperator->size(); ip++) {
auto sys_opr_opt = (*bOperator)[ip].get(i, j);
if (sys_opr_opt) {
auto *sys_opr = sys_opr_opt.value();
qmatrix<> aMatrix;
if (aOperator == nullptr) {
aMatrix = sys_opr_opt.value()->cTranspose();
} // Else this should be just boperator
// set kept-kept part operator
for (auto iv : currentKeptIndex[i]) {
for (auto iv_p : currentKeptIndex[j]) {
aMatrix.at(iv_p, iv) = 0;
}
}
auto rhoA = sys_opr->dot(rhoZero[j]);
auto ARho = rhoZero[i].dot(*sys_opr);
for (size_t iv = 0; iv < kpdim; iv++) {
for (size_t iv_p = 0; iv_p < kpdim_p; iv_p++) {
double aa{0};
double bbv{0};
// for (size_t iv_pp = 0; iv_pp < kpdim_p; iv_pp++) {
// aa += sys_opr->at(iv, iv_pp) * rhoZero[j](iv_pp, iv_p);
//}
// for (size_t iv_pp = 0; iv_pp < kpdim; iv_pp++) {
// bbv += rhoZero[i](iv, iv_pp) * sys_opr->at(iv_pp, iv_p);
//}
aa = rhoA(iv, iv_p) * aMatrix.at(iv_p, iv);
bbv = ARho(iv, iv_p) * aMatrix.at(iv_p, iv);
int tmindex = int(
std::log(1 +
(std::fabs(energyRescale *
(nrg_object->eigenvaluesQ[i][iv] -
nrg_object->eigenvaluesQ[j][iv_p])) /
minEnergy)) *
delE);
if (tmindex < energyPts && tmindex >= 0) {
if (std::fabs(aa) > spWeightErrorBar) {
positiveWeight[ip][tmindex] += std::fabs(aa);
specSum += std::fabs(aa);
// std::cout << "Raw: " << tmindex << " " << std::fabs(aa)
// << " " << std::fabs(bbv) << std::endl;
}
if (std::fabs(bbv) > spWeightErrorBar) {
negativeWeight[ip][tmindex] += std::fabs(bbv);
specSum += std::fabs(bbv);
}
}
}
}
}
}
// all the matrices are set
}
}
std::cout << nrg_object->nrg_iterations_cnt << "specSum: " << specSum
<< " Scale: " << energyRescale << std::endl;
// set the matrix elements
// end of lm loop
// End of matrix generation.
// Rotate the c operator in the eigen basis
}
//
void setTemperature(double at) { kBT = at; }
void setCurrentIndex() {
if (currentKeptIndex.empty()) {
for (size_t i = 0; i < nrg_object->current_sysmQ.size(); i++) {
currentKeptIndex.emplace_back();
}
} else {
currentKeptIndex = previoudKeptIndex;
}
previoudKeptIndex = nrg_object->eigenvaluesQ_kept_indices;
}
void setOperator(std::vector<qOperator> *bopr,
std::vector<qOperator> *aopr = nullptr) {
aOperator = aopr;
bOperator = bopr;
for (size_t i = 0; i < bOperator->size(); i++) {
positiveWeight.emplace_back(energyPts, 0);
negativeWeight.emplace_back(energyPts, 0);
}
}
template <typename filetype> void saveFinalData(filetype *pfile) {
std::vector<double> energyPoints(energyPts, 0);
for (int i = 0; i < energyPts; i++) {
energyPoints[i] = (minEnergy * std::exp(i * 1. / delE));
}
std::string hstr = "GreenFn";
pfile->write(energyPoints, hstr + "EnergyPoints");
pfile->write(positiveWeight, hstr + "PositiveWeight");
pfile->write(negativeWeight, hstr + "NegativeWeight");
for (size_t i = 0; i < bOperator->size(); i++) {
positiveWeight[i].clear();
negativeWeight[i].clear();
}
// pfile->write(vecPartitions, "vecPartitions");
}
private:
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
double globalGroundStateEnergy{0};
std::vector<qmatrix<>> reducedRho;
std::vector<double> vecPartitions;
// rho.B and B.rho
std::vector<std::vector<double>> positiveWeight;
std::vector<std::vector<double>> negativeWeight;
// value
int energyPts = 100000;
double maxEnergy = 10; // One decade more
double minEnergy = 1e-10;
double delE = (energyPts - 1.0) / (std::log(maxEnergy / minEnergy));
public: // Give access for openchain class
std::vector<std::vector<size_t>> currentKeptIndex;
std::vector<qmatrix<>> rhoZero;
std::vector<qOperator> *aOperator{};
std::vector<qOperator> *bOperator{};
bool lastiteration{true};
double energyRescale{1};
double energyErrorBar{1e-5};
};