Fixed Bug in GF calculation: Call to band diagonalization through Lap…

…ackpp was not returning the right eigenvectors, but the eigenvectors to an intermediat tridiagonal matrix.

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

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,13 +326,15 @@ 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);
           std::vector<std::complex<double>> ws(nws,
                                                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>(