Added minor modification to make Band Lanczos computation of the GF compatible with comlpex wave functions. The normal Lanczos calculation of off-diagonal elements is not appropriate for complex wave functions, though.

Author: CarlosMejZ

include/macis/gf/gf.hpp

@@ -638,6 +638,7 @@ void RunGFCalc(std::vector<std::vector<std::complex<double>>> &GF,
       for(int iw = 0; iw < ws.size(); iw++) GF[iw][i * nvecs + i] = tGF[iw];
       for(int j = i + 1; j < nvecs; j++) {
         // OFF DIAGONAL ELEMENTS
+	// NOTE: WORKS ONLY FOR REAL-VALUED WAVE FUNCTIONS.
         std::cout << "DOING ELEMENT (" << i << ", " << j << ")" << std::endl;
         for(size_t iii = 0; iii < nterms; iii++)
           twfn(iii) = wfns[i * nterms + iii] + wfns[j * nterms + iii];

src/macis/gf/bandlan.cxx

@@ -325,14 +325,32 @@ void BandResolvent(
       ofile.close();
     }
     std::cout << "DONE! COMPUTING RESOLVENT ...";
+    if( ispart )
+    {
 #pragma omp parallel for
-    for(int iw = 0; iw < ws.size(); iw++) {
-      for(int k = 0; k < nvecs; k++) {
-        for(int l = 0; l < nvecs; l++) {
-          for(int i_lan = 0; i_lan < nLanIts; i_lan++) {
-            res[iw][k * nvecs + l] += S[i_lan * nvecs + k] * 1. /
-                                      (ws[iw] + eigvals[i_lan]) *
-                                      S[i_lan * nvecs + l];
+      for(int iw = 0; iw < ws.size(); iw++) {
+        for(int k = 0; k < nvecs; k++) {
+          for(int l = 0; l < nvecs; l++) {
+            for(int i_lan = 0; i_lan < nLanIts; i_lan++) {
+              res[iw][k * nvecs + l] += S[i_lan * nvecs + k] * 1. /
+                                        (ws[iw] + eigvals[i_lan]) *
+                                        S[i_lan * nvecs + l];
+            }
+          }
+        }
+      }
+    }
+    else 
+    {
+#pragma omp parallel for
+      for(int iw = 0; iw < ws.size(); iw++) {
+        for(int k = 0; k < nvecs; k++) {
+          for(int l = 0; l < nvecs; l++) {
+            for(int i_lan = 0; i_lan < nLanIts; i_lan++) {
+              res[iw][l * nvecs + k] += S[i_lan * nvecs + k] * 1. /
+                                        (ws[iw] + eigvals[i_lan]) *
+                                        S[i_lan * nvecs + l];
+            }
           }
         }
       }