From e906bdaf551cba37ff3d338503a914ac66849201 Mon Sep 17 00:00:00 2001
From: FotisKalk <fokalkavouras@gmail.com>
Date: Thu, 16 May 2024 13:52:40 +0200
Subject: [PATCH] Removed QSTLS bug fix

---
 include/vector3D.hpp |  2 --
 src/qstls.cpp        | 16 +++++-----------
 src/vector3D.cpp     |  5 -----
 3 files changed, 5 insertions(+), 18 deletions(-)

diff --git a/include/vector3D.hpp b/include/vector3D.hpp
index 4e6a5f3..7421270 100644
--- a/include/vector3D.hpp
+++ b/include/vector3D.hpp
@@ -59,8 +59,6 @@ class Vector3D {
 
   // Set all values in the vector to num
   void fill(const double &num);
-  // Set all values of row j and column k to num of a specific index i
-  void fill(const size_t i, const double &num);
   // Set all entries of row i and column j to num
   void fill(const size_t i, const size_t j, const double &num);
   // Set all entries of row i and column j to be a copy of vector num
diff --git a/src/qstls.cpp b/src/qstls.cpp
index 99a7c28..5420c3d 100644
--- a/src/qstls.cpp
+++ b/src/qstls.cpp
@@ -318,15 +318,15 @@ void Qstls::readAdrFixedFile(Vector3D &res,
     res.resize(nx, nl, nx);
   }
   readDataFromBinary<Vector3D>(file, res);
-  if (!file) { throwError("Error in reading from file " + fileName); }
   file.close();
+  if (!file) { throwError("Error in reading from file " + fileName); }
   if (checkAdrFixed(wvg_, Theta_, nl_) != 0) {
     throwError("Fixed component of the auxiliary density response"
                " loaded from file is incompatible with input");
   }
 }
 
-int Qstls::checkAdrFixed(const std::vector<double> &wvg_,
+int Qstls::checkAdrFixed(const vector<double> &wvg_,
                          const double Theta_,
                          const int nl_) const {
   constexpr double tol = 1e-15;
@@ -334,13 +334,8 @@ int Qstls::checkAdrFixed(const std::vector<double> &wvg_,
   const double &wvgMaxDiff =
       std::abs(*std::max_element(wvgDiff.begin(), wvgDiff.end()));
   const bool consistentMatsubara = nl_ == in.getNMatsubara();
-  const bool consistentTheta = std::abs(Theta_ - in.getDegeneracy()) <= tol;
+  const bool consistentTheta = abs(Theta_ - in.getDegeneracy()) <= tol;
   const bool consistentGrid = wvg_.size() == wvg.size() && wvgMaxDiff <= tol;
-  if (!consistentMatsubara) {
-    std::cout << "Inconsistent Matsubara" << std::endl;
-  }
-  if (!consistentTheta) { std::cout << "Inconsistent Theta" << std::endl; }
-  if (!consistentGrid) { std::cout << "Inconsistent grid values" << std::endl; }
   if (!consistentMatsubara || !consistentTheta || !consistentGrid) { return 1; }
   return 0;
 }
@@ -579,7 +574,7 @@ void Adr::get(const vector<double> &wvg, const Vector3D &fixed, Vector2D &res) {
 void AdrFixed::get(vector<double> &wvg, Vector3D &res) const {
   const int nx = wvg.size();
   const int nl = res.size(1);
-  if (x == 0.0) { res.fill(0, 0.0); };
+  if (x == 0.0) { res.fill(0.0); };
   const double x2 = x * x;
   auto it = find(wvg.begin(), wvg.end(), x);
   assert(it != wvg.end());
@@ -609,20 +604,19 @@ double AdrFixed::integrand1(const double &q, const double &l) const {
 double
 AdrFixed::integrand2(const double &t, const double &y, const double &l) const {
   const double q = itg.getX();
+  if (q == 0 || t == 0 || y == 0) { return 0; };
   const double x2 = x * x;
   const double y2 = y * y;
   const double q2 = q * q;
   const double txq = 2.0 * x * q;
   if (l == 0) {
     if (t == txq) { return 2.0 * q2 / (y2 + 2.0 * txq - x2); };
-    if (x == y && t == 0.0) { return 1.0 / (y * q); };
     const double t2 = t * t;
     double logarg = (t + txq) / (t - txq);
     logarg = (logarg < 0.0) ? -logarg : logarg;
     return 1.0 / (2.0 * t + y2 - x2) *
            ((q2 - t2 / (4.0 * x2)) * log(logarg) + q * t / x);
   }
-  if (x == y && t == 0.0) { return 0.0; };
   const double tplT = 2.0 * M_PI * l * Theta;
   const double tplT2 = tplT * tplT;
   const double txqpt = txq + t;
diff --git a/src/vector3D.cpp b/src/vector3D.cpp
index e97141a..9170f2f 100644
--- a/src/vector3D.cpp
+++ b/src/vector3D.cpp
@@ -61,11 +61,6 @@ void Vector3D::fill(const double &num) {
   std::for_each(v.begin(), v.end(), [&](double &vi) { vi = num; });
 }
 
-void Vector3D::fill(const size_t i, const double &num) {
-  const auto dest = v.begin() + i * s2 * s3;
-  std::for_each(dest, dest + s2 * s3, [&](double &vi) { vi = num; });
-}
-
 void Vector3D::fill(const size_t i, const size_t j, const double &num) {
   const auto &dest = v.begin() + j * s3 + i * s2 * s3;
   std::for_each(dest, dest + s3, [&](double &vi) { vi = num; });