Refactored Self-Energy implementation

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();

include/xatu/Exciton.tpp

@@ -116,6 +116,17 @@ void Exciton<T>::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 <typename T>
+void Exciton<T>::setSelfEnergy(bool selfenergy){
+    this->selfenergy = selfenergy;
+}
+
 /*------------------------------------ Electron-hole pair basis ------------------------------------*/
 
 /**

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:

include/xatu/ExcitonTB.hpp

@@ -36,9 +36,11 @@ class ExcitonTB : public Exciton<SystemTB> {
         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<double> ftX;
         arma::mat potentialMat;
@@ -142,6 +144,8 @@ class ExcitonTB : public Exciton<SystemTB> {
     public:
         // BSE initialization and energies
         void initializeHamiltonian();
+        std::complex<double> 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<ResultTB> diagonalize(std::string method = "diag", int nstates = 8);

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);

main/xatu.cpp

@@ -28,6 +28,7 @@ int main(int argc, char* argv[]){
     TCLAP::ValueArg<int>    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<std::string> 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";

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,20 +150,27 @@ void ExcitonConfiguration::checkContentCoherence(){
 
     bool potentialFound = false;
     bool exchangePotentialFound = false;
+    bool selfenergyPotentialFound = false;
     for (auto potential : supportedPotentials){
         if(excitonInfo.potential == potential){
             potentialFound = true;
         }
         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'");
     }
     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'");
     }

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<double> ExcitonTB::selfenergyTerm(bool Qtoggle, uint32_t k_index, uint32_t kp_index, const arma::cx_vec& coefsK, const arma::cx_vec& coefsKp){
+
+    std::complex<double> 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<double> D, X = 0.0;
+        std::complex<double> 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<cx_mat>(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){

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 = 4;
+
+    std::string modelfile = "../examples/material_models/hBN.model";    
+    xatu::SystemConfiguration config(modelfile);
+
+    xatu::ExcitonTB exciton = xatu::ExcitonTB(config, ncell, 1, 0, {1, 1, 10});
+    exciton.setExchange(true);
+    exciton.setSelfEnergy(true);
+
+    exciton.brillouinZoneMesh(ncell);
+    exciton.initializeHamiltonian();
+    exciton.BShamiltonian();
+    auto results = exciton.diagonalize("diag", nstates);
+
+    auto energies = xatu::detectDegeneracies(results->eigval, nstates, 6);
+
+    std::vector<std::vector<double>> 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();