From e21c3132b07b969f468bbfd0a06d0bd367ddbfc7 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 20 Jan 2026 14:06:03 +0100 Subject: [PATCH 1/4] Refactored implementation of self-energy Commit implements a cleaner version of the self-energy contribution to the independent particle bands, as well as enabling outputing (to a file with extension .selfenergy) of the k-resolved self energy contribution to each band. Implementation is currently only implemented when the interactions are computed in real-space, although the generalization to reciprocal-space computation should be trivial. --- include/xatu/Exciton.hpp | 3 + include/xatu/Exciton.tpp | 11 +++ include/xatu/ExcitonConfiguration.hpp | 4 + include/xatu/ExcitonTB.hpp | 4 + include/xatu/Result.hpp | 1 + main/xatu.cpp | 13 +++ src/ExcitonConfiguration.cpp | 21 +++++ src/ExcitonTB.cpp | 126 ++++++++++++++++++++++++-- 8 files changed, 175 insertions(+), 8 deletions(-) diff --git a/include/xatu/Exciton.hpp b/include/xatu/Exciton.hpp index eb4a5c2..b73706c 100755 --- a/include/xatu/Exciton.hpp +++ b/include/xatu/Exciton.hpp @@ -41,6 +41,7 @@ class Exciton { // Flags bool exchange = false; + bool selfenergy = false; // Internals for BSE arma::mat eigvalKStack_, eigvalKQStack_; @@ -101,6 +102,8 @@ class Exciton { void setCutoff(double); void setScissor(double); void setExchange(bool); + void setSelfEnergy(bool); + protected: void initializeBasis(); diff --git a/include/xatu/Exciton.tpp b/include/xatu/Exciton.tpp index c385aaf..3a26e3d 100755 --- a/include/xatu/Exciton.tpp +++ b/include/xatu/Exciton.tpp @@ -116,6 +116,17 @@ void Exciton::setExchange(bool exchange){ this->exchange = exchange; } + +/** + * To toggle on or off the self-energy term in the interaction matrix elements. + * @param selfenergy Either true of false + * @return void + */ +template +void Exciton::setSelfEnergy(bool selfenergy){ + this->selfenergy = selfenergy; +} + /*------------------------------------ Electron-hole pair basis ------------------------------------*/ /** diff --git a/include/xatu/ExcitonConfiguration.hpp b/include/xatu/ExcitonConfiguration.hpp index 0d51e27..d957e94 100755 --- a/include/xatu/ExcitonConfiguration.hpp +++ b/include/xatu/ExcitonConfiguration.hpp @@ -52,6 +52,10 @@ class ExcitonConfiguration : public ConfigurationBase{ std::string exchangePotential = "keldysh"; // Regularization distance double regularization = 0.0; + // Flag to compute exciton energy while taking into account the Self-Energy contribution + bool selfenergy = false; + // Potential to use in selfenergy term + std::string selfenergyPotential = "keldysh"; }; public: diff --git a/include/xatu/ExcitonTB.hpp b/include/xatu/ExcitonTB.hpp index b38b3ef..d7d4f06 100644 --- a/include/xatu/ExcitonTB.hpp +++ b/include/xatu/ExcitonTB.hpp @@ -36,9 +36,11 @@ class ExcitonTB : public Exciton { std::string mode_ = "realspace"; std::string potential_ = "keldysh"; std::string exchangePotential_ = "keldysh"; + std::string selfenergyPotential_ = "keldysh"; // Internals for BSE arma::cx_mat ftMotifQ; + arma::cx_cube ftMotifQ3; arma::cx_cube ftMotifStack; std::complex ftX; arma::mat potentialMat; @@ -142,6 +144,8 @@ class ExcitonTB : public Exciton { public: // BSE initialization and energies void initializeHamiltonian(); + std::complex selfenergyTerm(bool, uint32_t, uint32_t, const arma::cx_vec&, const arma::cx_vec&); + void writeBandSelfEnergy(FILE*); void BShamiltonian(); void BShamiltonian(const arma::imat& basis); std::unique_ptr diagonalize(std::string method = "diag", int nstates = 8); diff --git a/include/xatu/Result.hpp b/include/xatu/Result.hpp index 9161cb7..339519b 100755 --- a/include/xatu/Result.hpp +++ b/include/xatu/Result.hpp @@ -59,6 +59,7 @@ class Result { void writeRealspaceAmplitude(int, int, const arma::rowvec&, FILE*, int); virtual void writeRealspaceAmplitude(const arma::cx_vec&, int, const arma::rowvec&, FILE*, int) = 0; virtual void writeAbsorptionSpectrum() = 0; + void writeSelfEnergy(FILE*); protected: int findExcitonPeak(int); diff --git a/main/xatu.cpp b/main/xatu.cpp index f4e43c2..2b7f145 100755 --- a/main/xatu.cpp +++ b/main/xatu.cpp @@ -28,6 +28,7 @@ int main(int argc, char* argv[]){ TCLAP::ValueArg w90Arg("w", "w90", "Indicates that the system file is a tb.dat Wannier90 file.", false, -1, "No. electrons ", cmd); TCLAP::SwitchArg absorptionArg("a", "absorption", "Computes the absorption spectrum.", cmd, false); TCLAP::ValueArg formatArg("f", "format", "Format of the input system file.", false, "model", "model or hdf5", cmd); + TCLAP::SwitchArg selfenergyprint("i", "printSelfEnergy", "Prints the contribution of self energy to the single-particle bands.", cmd, false); TCLAP::AnyOf outputOptions; TCLAP::SwitchArg energyArg("e", "energy", "Write energies.", false); @@ -173,6 +174,18 @@ int main(int argc, char* argv[]){ fclose(textfile_en); } + bool writeSelfEnergy = selfenergyprint.isSet(); + if (writeSelfEnergy){ + std::string filename_selfen = output + ".selfenergy"; + FILE* textfile_selfen = fopen(filename_selfen.c_str(), "w"); + + std::cout << "Writing self energy to file: " << filename_selfen << std::endl; + + bulkExciton.writeBandSelfEnergy(textfile_selfen); + + fclose(textfile_selfen); + } + bool writeStates = eigenstatesArg.isSet(); if(writeStates){ std::string filename_st = output + ".states"; diff --git a/src/ExcitonConfiguration.cpp b/src/ExcitonConfiguration.cpp index 7a37951..63b5311 100755 --- a/src/ExcitonConfiguration.cpp +++ b/src/ExcitonConfiguration.cpp @@ -97,6 +97,20 @@ void ExcitonConfiguration::parseContent(){ std::cout << arg << " " << str << std::endl; excitonInfo.exchangePotential = str; } + else if(arg == "selfenergy"){ + std::string str = parseWord(content[0]); + if ((str != "true") && (str != "false")){ + throw std::invalid_argument("Self-Energy option must be set to 'true' or 'false'."); + } + if (str == "true"){ + excitonInfo.selfenergy = true; + } + } + else if(arg == "selfenergy.potential"){ + std::string str = parseWord(content[0]); + std::cout << arg << " " << str << std::endl; + excitonInfo.selfenergyPotential = str; + } else if(arg == "potential"){ std::string str = parseWord(content[0]); excitonInfo.potential = str; @@ -136,6 +150,7 @@ void ExcitonConfiguration::checkContentCoherence(){ bool potentialFound = false; bool exchangePotentialFound = false; + bool selfenergyPotentialFound = false; for (auto potential : supportedPotentials){ if(excitonInfo.potential == potential){ potentialFound = true; @@ -143,6 +158,9 @@ void ExcitonConfiguration::checkContentCoherence(){ if(excitonInfo.exchange && excitonInfo.exchangePotential == potential){ exchangePotentialFound = true; } + if(excitonInfo.selfenergy && excitonInfo.selfenergyPotential == potential){ + selfenergyPotentialFound = true; + } } if (!potentialFound){ throw std::invalid_argument("Specified 'potential' not supported. Use 'keldysh' or 'coulomb'"); @@ -150,6 +168,9 @@ void ExcitonConfiguration::checkContentCoherence(){ if (excitonInfo.exchange && !exchangePotentialFound){ throw std::invalid_argument("Specified 'exchange.potential' not supported. Use 'keldysh' or 'coulomb'"); } + if (excitonInfo.selfenergy && !selfenergyPotentialFound){ + throw std::invalid_argument("Specified 'selfenergy.potential' not supported. Use 'keldysh' or 'coulomb'"); + } if (excitonInfo.mode != "realspace" && excitonInfo.mode != "reciprocalspace"){ throw std::invalid_argument("Invalid mode. Use 'realspace' or 'reciprocalspace'"); } diff --git a/src/ExcitonTB.cpp b/src/ExcitonTB.cpp index db64500..3360a25 100644 --- a/src/ExcitonTB.cpp +++ b/src/ExcitonTB.cpp @@ -83,6 +83,7 @@ void ExcitonTB::initializeExcitonAttributes(const ExcitonConfiguration& cfg){ // Set flags this->exchange = cfg.excitonInfo.exchange; + this->selfenergy = cfg.excitonInfo.selfenergy; this->scissor_ = cfg.excitonInfo.scissor; this->mode_ = cfg.excitonInfo.mode; this->nReciprocalVectors_ = cfg.excitonInfo.nReciprocalVectors; @@ -92,6 +93,8 @@ void ExcitonTB::initializeExcitonAttributes(const ExcitonConfiguration& cfg){ } this->potential_ = cfg.excitonInfo.potential; this->exchangePotential_ = cfg.excitonInfo.exchangePotential; + this->selfenergyPotential_ = cfg.excitonInfo.selfenergyPotential; + } /** @@ -769,6 +772,9 @@ void ExcitonTB::initializeResultsH0(){ this->ftMotifStack = arma::cx_cube(natoms, natoms, system->meshBZ.n_rows); this->ftMotifQ = arma::cx_mat(natoms, natoms); + if(this->selfenergy){ + this->ftMotifQ3 = arma::cx_cube(natoms, natoms, system->meshBZ.n_rows); + } vec auxEigVal(basisdim); arma::cx_mat auxEigvec(basisdim, basisdim); arma::cx_mat h; @@ -813,7 +819,18 @@ void ExcitonTB::initializeResultsH0(){ #pragma omp parallel for for (unsigned int i = 0; i < system->meshBZ.n_rows; i++){ // BIGGEST BOTTLENECK OF THE CODE - initializeMotifFT(i, cells, directPotential); + initializeMotifFT(i, cells, directPotential); + if(this->selfenergy){ + if(arma::norm(Q) != 0){ + potptr selfenergyPotential = selectPotential(this->selfenergyPotential_); + this->ftMotifQ3.slice(i) = motifFTMatrix(system->kpoints.row(i) - this->Q, cells, selfenergyPotential); + } + else { + potptr selfenergyPotential = selectPotential(this->selfenergyPotential_); + this->ftMotifQ3.slice(i) = this->ftMotifStack.slice(i); + }; + } + /* AJU 24-11-23: Progress bar does not work properly with parallel for Fix it or remove it direcly? */ @@ -872,6 +889,55 @@ void ExcitonTB::BShamiltonian(){ BShamiltonian(basis); } +/** + * Routine to compute the self-energy contribution to the bands used in the exciton computation. + * @param Qtoggle toggle on whether to use Q in calculations. it is only set to zero when computing the correction to the band structure for exporting. + * @param k/kp_index index of k/k' point for the computation of the shifted mot + * @param coefsK/Kp k/k' points in the definition of self-energy contribution + * @return self energy contribution + */ +std::complex ExcitonTB::selfenergyTerm(bool Qtoggle, uint32_t k_index, uint32_t kp_index, const arma::cx_vec& coefsK, const arma::cx_vec& coefsKp){ + + std::complex self = 0.0; + if (this->selfenergy){ + uint32_t k_index_0 = system_->findEquivalentPointBZ(arma::rowvec{0,0,0}, ncell); + arma::cx_mat motifFT0 = ftMotifStack.slice(k_index_0); + arma::cx_vec coefsK3; + arma::cx_mat motifFT3; + if (mode == "realspace"){ + for (int v3 = 0; v3 < (int)valenceBands.n_elem; v3++){ + for (uint32_t k3_index = 0; k3_index < system->nk; k3_index++){ + if (gauge == "atomic"){ + coefsK3 = system_->latticeToAtomicGauge(eigvecKStack.slice(k3_index).col(v3), system->kpoints.row(k3_index)); + } + else{ + coefsK3 = eigvecKStack.slice(k3_index).col(v3); + } + uint32_t effective_k3_index = system_->findEquivalentPointBZ(system->kpoints.row(k3_index) - system->kpoints.row(k_index), ncell); + arma::cx_mat motifFT3 = (Qtoggle) ? this->ftMotifQ3.slice(effective_k3_index) : ftMotifStack.slice(effective_k3_index); + // if (Qtoggle){ + // arma::cx_mat motifFT3 = this->ftMotifQ3.slice(effective_k3_index); + // } + // else{ + // arma::cx_mat motifFT3 = ftMotifStack.slice(effective_k3_index); + // + // } + self = self + realSpaceInteractionTerm(coefsK, coefsK3 , coefsKp, coefsK3, motifFT0) - realSpaceInteractionTerm(coefsK, coefsK3, coefsK3, coefsKp, motifFT3); + + }; + }; + return self; + } + else if (mode == "reciprocalspace"){ + std::cout << self << std::endl; + return self; + } + } + else{ + return self; + } +}; + /** * Initialize BSE hamiltonian matrix and kinetic matrix. * @details Instead of calculating the energies and coeficients dinamically, which @@ -891,7 +957,6 @@ void ExcitonTB::BShamiltonian(const arma::imat& basis){ if (!basis.is_empty()){ basisStates = basis; }; - uint64_t basisDimBSE = basisStates.n_rows; std::cout << "BSE dimension: " << basisDimBSE << std::endl; std::cout << "Initializing Bethe-Salpeter matrix... " << std::flush; @@ -921,7 +986,6 @@ void ExcitonTB::BShamiltonian(const arma::imat& basis){ int v2 = bandToIndex[basisStates(j, 0)]; int c2 = bandToIndex[basisStates(j, 1)]; uint32_t k2Q_index = k2_index; - // Using the atomic gauge if(gauge == "atomic"){ coefsK = system_->latticeToAtomicGauge(eigvecKStack.slice(k_index).col(v), system->kpoints.row(k_index)); @@ -936,14 +1000,27 @@ void ExcitonTB::BShamiltonian(const arma::imat& basis){ coefsK2Q = eigvecKQStack.slice(k2Q_index).col(c2); } - std::complex D, X = 0.0; + std::complex D, X, selfcond, selfval = 0.0; if (mode == "realspace"){ uint32_t effective_k_index = system_->findEquivalentPointBZ(system->kpoints.row(k2_index) - system->kpoints.row(k_index), ncell); arma::cx_mat motifFT = ftMotifStack.slice(effective_k_index); D = realSpaceInteractionTerm(coefsKQ, coefsK2, coefsK2Q, coefsK, motifFT); if(this->exchange){ X = realSpaceInteractionTerm(coefsKQ, coefsK2, coefsK, coefsK2Q, this->ftMotifQ); - } + } + //self-energy terms + if(this->selfenergy){ + bool testval = (kQ_index==k2Q_index and v==v2); + bool testcond = (k_index==k2_index and c==c2); + if (testval or testcond){ + if (testval){ + selfcond = selfenergyTerm(true, kQ_index, k2_index, coefsKQ, coefsK2Q); + } + if (testcond){ + selfval = selfenergyTerm(false, k2_index, k_index, coefsK2, coefsK); + } + } + } } else if (mode == "reciprocalspace"){ arma::rowvec k = system->kpoints.row(k_index); @@ -955,9 +1032,8 @@ void ExcitonTB::BShamiltonian(const arma::imat& basis){ } if (i == j){ - HBS_(i, j) = (this->scissor + - eigvalKQStack.col(kQ_index)(c) - eigvalKStack.col(k_index)(v))/2. - - (D - X)/2.; + HBS_(i, j) = (this->scissor + (eigvalKQStack.col(kQ_index)(c) + selfcond) - (eigvalKStack.col(k_index)(v) + selfval))/2. + - (D - X)/2.; } else{ HBS_(i, j) = - (D - X); @@ -968,6 +1044,36 @@ void ExcitonTB::BShamiltonian(const arma::imat& basis){ std::cout << "Done" << std::endl; }; +/** + * Routine to write the self energy contribution to each band to a file. Creates a matrix named selfen with nk rows and bands columns. + * Each row is organized as selfenergy(band0) selfenergy(band1) .... + * @param file Pointer to file. + * + * + * @return void + */ +void ExcitonTB::writeBandSelfEnergy(FILE* file){ + arma::cx_vec coefsK; + arma::cx_mat selfen = arma::zeros(system->nk, (int)bands.n_elem); + for (uint32_t k_index = 0; k_index < system->nk; k_index++){ + arma::rowvec kpoint = system->kpoints.row(k_index); + fprintf(file, "%11.7lf\t%11.7lf\t%11.7lf\t", kpoint(0), kpoint(1), kpoint(2)); + + for (int bandindex = 0; bandindex < (int)bands.n_elem; bandindex++){ + if (gauge == "atomic"){ + coefsK = system_->latticeToAtomicGauge(eigvecKStack.slice(k_index).col(bandindex), system->kpoints.row(k_index)); + } + else{ + coefsK = eigvecKStack.slice(k_index).col(bandindex); + } + selfen(k_index, bandindex) = selfenergyTerm(false, k_index, k_index, coefsK, coefsK); + fprintf(file, "%11.7lf\t%11.7lf\t", real(selfen(k_index, bandindex)), imag(selfen(k_index, bandindex))); + } + fprintf(file, "\n"); + } +}; + + /** * Routine to diagonalize the BSE and return a Result object. * @param method Method to diagonalize the BSE, either 'diag' (standard diagonalization) @@ -1428,6 +1534,10 @@ void ExcitonTB::printInformation(){ cout << std::left << std::setw(30) << "Exchange: " << (exchange ? "True" : "False") << endl; cout << std::left << std::setw(30) << "Exchange potential: " << exchangePotential_ << endl; } + if(selfenergy){ + cout << std::left << std::setw(30) << "Self-Energy: " << (selfenergy ? "True" : "False") << endl; + cout << std::left << std::setw(30) << "Self-Energy potential: " << selfenergyPotential_ << endl; + } if(arma::norm(Q) > 1E-7){ cout << std::left << std::setw(30) << "Q: "; for (auto qi : Q){ From 98ccedd94eac3c90891d05ebd11b3b476952ee33 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 20 Jan 2026 14:31:14 +0100 Subject: [PATCH 2/4] Add test draft for self-energy in hBN Simple test --- test/tests.cpp | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/test/tests.cpp b/test/tests.cpp index 934c905..7120398 100644 --- a/test/tests.cpp +++ b/test/tests.cpp @@ -801,6 +801,39 @@ TEST(hBNTest, AnisotropicConductivity){ } +TEST(hBNTest, SelfEnergyContribution){ + + std::cout.clear(); + std::cout << std::setw(50) << std::left << "Testing TB hBN self-energy contribution (fulldiag)... " << std::endl; + std::cout.setstate(std::ios_base::failbit); + + int ncell = 30; + int nstates = 6; + + std::string modelfile = "../examples/material_models/hBN.model"; + xatu::SystemConfiguration config(modelfile); + exciton.setExchange(true); + exciton.setSelfEnergy(true); + + xatu::ExcitonTB exciton = xatu::ExcitonTB(config, ncell, 1, 0, {1, 1, 10}); + + exciton.brillouinZoneMesh(ncell); + exciton.initializeHamiltonian(); + exciton.BShamiltonian(); + auto results = exciton.diagonalize("diag", nstates); + + auto energies = xatu::detectDegeneracies(results->eigval, nstates, 6); + + std::vector> expectedEnergies = {{4.978280, 2}, + {5.790652, 1}, + {5.807516, 1}}; + + for(uint i = 0; i < energies.size(); i++){ + EXPECT_NEAR(energies[i][0], expectedEnergies[i][0], 1E-4); + EXPECT_EQ(energies[i][1], expectedEnergies[i][1]); + } +} + TEST(MoS2Test, FullDiagonalization){ std::cout.clear(); From ecd7c641b3fc2270d273c90a3d04a8a22a7661ee Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 20 Jan 2026 14:57:55 +0100 Subject: [PATCH 3/4] fixed me being dumb --- test/tests.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/tests.cpp b/test/tests.cpp index 7120398..5339d86 100644 --- a/test/tests.cpp +++ b/test/tests.cpp @@ -812,10 +812,10 @@ TEST(hBNTest, SelfEnergyContribution){ std::string modelfile = "../examples/material_models/hBN.model"; xatu::SystemConfiguration config(modelfile); - exciton.setExchange(true); - exciton.setSelfEnergy(true); xatu::ExcitonTB exciton = xatu::ExcitonTB(config, ncell, 1, 0, {1, 1, 10}); + exciton.setExchange(true); + exciton.setSelfEnergy(true); exciton.brillouinZoneMesh(ncell); exciton.initializeHamiltonian(); From 7f190f9933444b0974cedd6aff35a5126085e47f Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 20 Jan 2026 15:11:24 +0100 Subject: [PATCH 4/4] fixed me being dumb v2 --- test/tests.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tests.cpp b/test/tests.cpp index 5339d86..1cc6111 100644 --- a/test/tests.cpp +++ b/test/tests.cpp @@ -808,7 +808,7 @@ TEST(hBNTest, SelfEnergyContribution){ std::cout.setstate(std::ios_base::failbit); int ncell = 30; - int nstates = 6; + int nstates = 4; std::string modelfile = "../examples/material_models/hBN.model"; xatu::SystemConfiguration config(modelfile);