From ff65ee7ed5316da640329a5701025f43cce4c304 Mon Sep 17 00:00:00 2001
From: Carlos Mejuto Zaera <cmejutoz@sissa.it>
Date: Thu, 5 Oct 2023 08:58:43 +0200
Subject: [PATCH] Fixed Bug in GF calculation: Call to band diagonalization
 through Lapackpp was not returning the right eigenvectors, but the
 eigenvectors to an intermediat tridiagonal matrix.

---
 src/macis/gf/bandlan.cxx    | 3 ++-
 tests/standalone_driver.cxx | 9 +++++++++
 2 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/src/macis/gf/bandlan.cxx b/src/macis/gf/bandlan.cxx
index 411f6363..99bfb658 100644
--- a/src/macis/gf/bandlan.cxx
+++ b/src/macis/gf/bandlan.cxx
@@ -138,7 +138,7 @@ bool GetEigsysBand(std::vector<double> &mat, int nSupDiag,
   // PREPARE VARIABLES FOR LAPACK
   lapack::Uplo UPLO = lapack::Uplo::Upper;
   lapack::Job VECT = lapack::Job::Vec;
-  lapack::Job COMPZ = lapack::Job::Vec;
+  lapack::Job COMPZ = lapack::Job::UpdateVec;
   int N = matsize, LDQ = matsize, LDAB = nSupDiag + 1;
   std::vector<double> AB((nSupDiag + 1) * N, 0.);
   std::vector<double> D(N, 0.), E(N - 1, 0.), Q(N * N, 0.);
@@ -225,6 +225,7 @@ void BandResolvent(
   int nbands = nvecs;
   BandLan<double>(Hop, vecs, bandH, nLanIts, nbands, len_vec, 1.E-6, print);
   std::cout << "DONE! ";
+  
   if(print) {
     std::ofstream ofile("BLH.dat", std::ios::out);
     ofile.precision(dbl::max_digits10);
diff --git a/tests/standalone_driver.cxx b/tests/standalone_driver.cxx
index e0f32f65..5a3f35c2 100644
--- a/tests/standalone_driver.cxx
+++ b/tests/standalone_driver.cxx
@@ -18,6 +18,7 @@
 #include <macis/asci/refine.hpp>
 #include <macis/gf/gf.hpp>
 #include <macis/hamiltonian_generator/double_loop.hpp>
+#include <macis/hamiltonian_generator/sd_build.hpp>
 #include <macis/util/cas.hpp>
 #include <macis/util/detail/rdm_files.hpp>
 #include <macis/util/fcidump.hpp>
@@ -325,6 +326,7 @@ int main(int argc, char** argv) {
           double wmin = -8.;
           double wmax = 8.;
           size_t nws = 1001;
+	  double beta = 157.;
           double eta = 0.1;
           std::complex<double> w0(wmin, eta);
           std::complex<double> wf(wmax, eta);
@@ -332,6 +334,7 @@ int main(int argc, char** argv) {
                                                std::complex<double>(0., 0.));
           for(int i = 0; i < nws; i++)
             ws[i] = w0 + (wf - w0) / double(nws - 1) * double(i);
+            //ws[i] = std::complex<double>(0., 1.) * (2. * i + 1.) * M_PI / beta; 
 
           // MCSCF Settings
           macis::GFSettings gf_settings;
@@ -361,6 +364,12 @@ int main(int argc, char** argv) {
 
           // Occupation numbers
           std::vector<double> occs(n_active, 1.);
+	  for( int i = 0; i < n_active; i++ )
+	    occs[i] = active_ordm[i + n_active * i] / 2.;
+	  std::cout << "Occupation Nrs.: ";
+	  for( int i = 0; i < n_active; i++ )
+	    std::cout << "  " << occs[i];
+	  std::cout << std::endl; 
 
           // GS vector
           Eigen::VectorXd psi0 = Eigen::Map<Eigen::VectorXd, Eigen::Unaligned>(