Compatibility for new and old grid files

src/green/sc/common_defs.h

@@ -22,6 +22,7 @@
 #define GREEN_SC_COMMON_DEFS_H
 
 #include <green/ndarray/ndarray.h>
+#include <green/params/params.h>
 
 namespace green::sc {
   /**
@@ -82,5 +83,35 @@ namespace green::sc {
     p.define<std::string>("input_file,", "File with input data", "input.h5");
     p.define<bool>("const_density", "Maintain constant number of electrons through iterations", true);
   }
+
+  inline int compare_version_strings(const std::string& v1, const std::string& v2) {
+    int major_V1 = 0, minor_V1 = 0, patch_V1 = 0;
+    int major_V2 = 0, minor_V2 = 0, patch_V2 = 0;
+
+    char suffix_V1[32] = "";
+    char suffix_V2[32] = "";
+
+    int parsed_1 = std::sscanf(v1.c_str(), "%d.%d.%d%30s", &major_V1, &minor_V1, &patch_V1, suffix_V1);
+    int parsed_2 = std::sscanf(v2.c_str(), "%d.%d.%d%30s", &major_V2, &minor_V2, &patch_V2, suffix_V2);
+
+    if (parsed_1 < 3) {
+      throw std::runtime_error("First version string (v1) failed to parse: '" + v1 + "'. Expected format: major.minor.patch[suffix]");
+    }
+    if (parsed_2 < 3) {
+      throw std::runtime_error("Second version string (v2) failed to parse: '" + v1 + "'. Expected format: major.minor.patch[suffix]");
+    }
+
+    if (major_V1 != major_V2) {
+      return major_V1 > major_V2 ? 1 : -1;
+    }
+    if (minor_V1 != minor_V2) {
+      return minor_V1 > minor_V2 ? 1 : -1;
+    }
+    if (patch_V1 != patch_V2) {
+      return patch_V1 > patch_V2 ? 1 : -1;
+    }
+
+    return 0;
+  }
 }  // namespace green::sc
 #endif  // GREEN_SC_COMMON_DEFS_H

src/green/sc/except.h

@@ -35,6 +35,11 @@ namespace green::sc {
     explicit sc_diis_vsp_error(const std::string& what) : std::runtime_error(what) {}
   };
 
+  class outdated_results_file_error : public std::runtime_error {
+  public:
+    explicit outdated_results_file_error(const std::string& what) : std::runtime_error(what) {}
+  };
+
 }  // namespace green::sc
 
 #endif  // GREEN_SC_EXCEPT_H

src/green/sc/sc_loop.h

@@ -27,6 +27,7 @@
 #include <green/params/params.h>
 #include <green/utils/mpi_utils.h>
 #include <green/utils/timing.h>
+#include <green/grids/common_defs.h>
 #include <mpi.h>
 
 #include <cstdio>
@@ -37,6 +38,7 @@
 #include "common_defs.h"
 #include "mixing.h"
 #include "solver.h"
+#include "except.h"
 
 namespace green::sc {
   /**
@@ -74,7 +76,11 @@ namespace green::sc {
   public:
     sc_loop(MPI_Comm comm, params::params& p) :
         _itermax(p["itermax"]), _iter(0), _e_thr(p["threshold"]), _e_thr_sp(p["E_thr_sp"]), _input_path(p["input_file"]),
-        _results_file(p["results_file"]), _restart(p["restart"]), _dyson_solver(p), _mix(p), _context(comm) {}
+        _results_file(p["results_file"]), _restart(p["restart"]), _dyson_solver(p), _mix(p), _context(comm) {
+      if (_restart) {
+        check_grids_version_consistency();
+      }
+    }
 
     virtual ~sc_loop() = default;
 
@@ -164,13 +170,64 @@ namespace green::sc {
      */
     void dump_iteration(size_t iter, G& g_tau, S1& sigma_1, St& sigma_tau) const {
       h5pp::archive ar(_results_file, "a");
+      if (!ar.has_attribute("__grids_version__")) {
+        // Grids attribute -- will only be performed once
+        const std::string& grid_version = _dyson_solver.get_grids_version();
+        ar.set_attribute<std::string>("__grids_version__", grid_version);
+      }
       ar["iter"] << iter;
       internal::write(sigma_1, "iter" + std::to_string(iter) + "/Sigma1", ar);
       internal::write(sigma_tau, "iter" + std::to_string(iter) + "/Selfenergy/data", ar);
       internal::write(g_tau, "iter" + std::to_string(iter) + "/G_tau/data", ar);
       ar.close();
     }
 
+    /**
+     * @brief Checks consistency between grid-file version used in current run vs. the version
+     * used to generate the results file that is being restarted from.
+     * This is to prevent users from accidentally restarting from a results file that was generated with
+     * an older version of green-grids, which can lead to silent errors in the results.
+     * 
+     * 1. If the results file does not have a green-grids version attribute, treat it as having been
+     *    generated with the baseline grids version (GRIDS_MIN_VERSION, currently 0.2.4) and check
+     *    that the current grid-file version is not newer; otherwise an error is raised.
+     * 2. If the results file has a green-grids version attribute, compare that version with the current
+     *    grid-file version and distinguish older, equal, and newer cases, allowing only compatible
+     *    combinations and throwing specific errors when there is a mismatch.
+     */
+    void check_grids_version_consistency() {
+      // If results file does not exist, nothing to check
+      if (!std::filesystem::exists(_results_file)) return;
+
+      // Read grid-file version
+      const std::string& grid_file_version = _dyson_solver.get_grids_version();
+
+      h5pp::archive ar(_results_file, "r");
+      if (ar.has_attribute("__grids_version__")) {
+        std::string grids_version_in_results;
+        grids_version_in_results = ar.get_attribute<std::string>("__grids_version__");
+        ar.close(); // safely close before throwing
+        if (compare_version_strings(grid_file_version, grids_version_in_results) < 0) {
+          throw green::grids::outdated_grids_file_error("The current green-grids version (" + grid_file_version +
+                                            ") is older than the green-grids version used to create the original results file ("
+                                            + grids_version_in_results +
+                                            "). Please update green-grids to version " + grids_version_in_results);
+        } else if (compare_version_strings(grid_file_version, grids_version_in_results) > 0) {
+          throw outdated_results_file_error("The green-grids version used to create the results file (" + grids_version_in_results +
+                                            ") is older than the current specified grid file (" + grid_file_version +
+                                            "). Please download the appropriate version from: " +
+                                            "https://github.com/Green-Phys/green-grids/releases/ or https://github.com/Green-Phys/green-grids/tags");
+        }
+      } else if (compare_version_strings(grid_file_version, green::grids::GRIDS_MIN_VERSION) > 0) {
+        ar.close(); // safely close before throwing
+        throw outdated_results_file_error("The results file was created using un-versioned grid file (equiv. to " + green::grids::GRIDS_MIN_VERSION +
+                                          ") and the current green-grids version (" + grid_file_version + ") is newer.\n" + 
+                                          "Please use old grid files from: https://github.com/Green-Phys/green-grids/releases/tag/v0.2.4.");
+      } else {
+        ar.close();
+      }
+    }
+
     DysonSolver& dyson_solver() { return _dyson_solver; }
   };
 }  // namespace green::sc

test/CMakeLists.txt

@@ -12,6 +12,7 @@ FetchContent_MakeAvailable(Catch2)
 list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras)
 
 add_executable(sc_test sc_test.cpp)
+target_compile_definitions(sc_test PRIVATE TEST_PATH="${CMAKE_CURRENT_SOURCE_DIR}/data")
 target_link_libraries(sc_test
         PRIVATE
         Catch2::Catch2

test/sc_test.cpp

@@ -33,6 +33,10 @@ inline std::string random_name() {
   return str.substr(0, 32) + ".h5";  // assumes 32 < number of characters in str
 }
 
+// Grid file paths for testing
+const std::string GRID_FILE_NEW = std::string(TEST_PATH) + "/1e4.h5";
+const std::string GRID_FILE_OLD = std::string(TEST_PATH) + "/old_1e4.h5";
+
 // Try to solve equation B*X^4 + C*X^2 - A*X + D = 0
 // It leads to following itertative scheme:
 // X = \alpha + \beta X^4 + \gamma X^2,
@@ -51,7 +55,7 @@ class fourth_power_equation_dyson {
   using Sigma1    = double;
   using Sigma_tau = double;
 
-  fourth_power_equation_dyson(green::params::params& p) : _alpha(p["alpha"]), _beta(p["beta"]) {}
+  fourth_power_equation_dyson(green::params::params& p) : _alpha(p["alpha"]), _beta(p["beta"]), _ft(p) {}
 
   void solve(G& g, Sigma1& sigma1, Sigma_tau& sigma_tau) {
     double g_new = _alpha + _alpha * (sigma1 + sigma_tau) * g;
@@ -62,12 +66,15 @@ class fourth_power_equation_dyson {
   void   dump_iteration(size_t, const G&, const Sigma1&, const Sigma_tau&, const std::string&) {};
   double mu() const { return 0; }
   double& mu() { return _mu; }
+  const green::grids::transformer_t& ft() const { return _ft; }
+  const std::string& get_grids_version() const { return _ft.get_version(); }
 
 private:
   double _alpha;
   double _beta;
   double _diff;
   double _mu;
+  green::grids::transformer_t _ft;
 };
 
 class fourth_power_equation_solver {
@@ -102,8 +109,9 @@ void solve_with_mixing(const std::string& mixing_type, const std::string& mixing
   auto        p        = green::params::params("DESCR");
   std::string res_file = random_name();
   std::string args     = "test --restart 0 --itermax 1000 --E_thr 1e-13 --mixing_type " + mixing_type + " --mixing_weight " + mixing_weight +
-                     " --verbose 1 --results_file=" + res_file;
+                     " --verbose 1 --results_file=" + res_file + " --grid_file ir/1e4.h5 --BETA 100";
   green::sc::define_parameters(p);
+  green::grids::define_parameters(p);
   p.define<double>("alpha", "", 0.45);
   p.define<double>("beta", "", 0.5);
   p.define<double>("gamma", "", 0.25);
@@ -150,8 +158,10 @@ TEST_CASE("Self-consistency") {
   SECTION("Solve simple") {
     auto        p        = green::params::params("DESCR");
     std::string res_file = random_name();
-    std::string args     = "test --verbose 1 --restart 0 --mixing_type NO_MIXING --itermax 100 --E_thr 1e-13 --results_file=" + res_file;
+    std::string args     = "test --verbose 1 --restart 0 --mixing_type NO_MIXING --itermax 100 "s
+        + "--grid_file ir/1e4.h5 --BETA 100 --E_thr 1e-13 --results_file=" + res_file;
     green::sc::define_parameters(p);
+    green::grids::define_parameters(p);
     p.define<double>("alpha", "", 0.45);
     p.define<double>("beta", "", 0.5);
     p.define<double>("gamma", "", 0.25);
@@ -210,10 +220,14 @@ TEST_CASE("Self-consistency") {
     auto        p2         = green::params::params("DESCR");
     std::string res_file_1 = random_name();
     std::string res_file_2 = random_name();
-    std::string args_1     = "test --restart 0 --itermax 4 --E_thr 1e-13 --results_file=" + res_file_1;
-    std::string args_2     = "test --restart 1 --itermax 2 --E_thr 1e-13 --results_file=" + res_file_2;
+    std::string args_1     = "test --restart 0 --itermax 4 --E_thr 1e-13 --grid_file " + std::string(GRID_FILE_NEW)
+                              + " --BETA 100 --results_file=" + res_file_1;
+    std::string args_2     = "test --restart 1 --itermax 2 --E_thr 1e-13 --grid_file " + std::string(GRID_FILE_NEW)
+                              + " --BETA 100 --results_file=" + res_file_2;
     green::sc::define_parameters(p);
     green::sc::define_parameters(p2);
+    green::grids::define_parameters(p);
+    green::grids::define_parameters(p2);
     p.define<double>("alpha", "", 0.45);
     p.define<double>("beta", "", 0.5);
     p.define<double>("gamma", "", 0.25);
@@ -232,30 +246,50 @@ TEST_CASE("Self-consistency") {
     double g2          = p2["x0"];
     double sigma1_2    = 0;
     double sigma_tau_2 = 0;
+
+    /**
+     * CASE 1: Compare a fresh start with 4 iterations against a 2x (restart + 2 iterations)
+     *         to validate restart correctness and stability.
+     */
     {
+      // First run: generate a results file from a fresh start with 4 iterations
       green::sc::sc_loop<fourth_power_equation_dyson> sc(MPI_COMM_WORLD, p);
       fourth_power_equation_solver                    solver(p["alpha"], p["beta"]);
       sc.solve(solver, h0, ovlp, g, sigma1, sigma_tau);
     }
     {
+      // Second run: restart from a results file that does not exist -- should work as a fresh start -- and run 2 iterations.
       green::sc::sc_loop<fourth_power_equation_dyson> sc(MPI_COMM_WORLD, p2);
       fourth_power_equation_solver                    solver(p2["alpha"], p2["beta"]);
       sc.solve(solver, h0, ovlp, g2, sigma1_2, sigma_tau_2);
     }
     {
+      // Third run: restart again to confirm restart is repeatable and stable. Run 2 more iterations
       green::sc::sc_loop<fourth_power_equation_dyson> sc(MPI_COMM_WORLD, p2);
       fourth_power_equation_solver                    solver(p2["alpha"], p2["beta"]);
       sc.solve(solver, h0, ovlp, g2, sigma1_2, sigma_tau_2);
     }
+    // Final value of 2 x (restart + 2 iterations) should be equal to (fresh start + 4 iterations)
     REQUIRE(std::abs(sigma_tau_2 - sigma_tau) < 1e-14);
+
+    /**
+     * CASE 2: Restart from an empty file
+     *
+     */
     std::filesystem::remove(res_file_2);
     {
-      // create empty file to check that
+      // Read __grids_version__ from the new grid files
+      green::h5pp::archive ar_grid(GRID_FILE_NEW, "r");
+      std::string         new_grid_file_version = ar_grid.get_attribute<std::string>("__grids_version__");
+      ar_grid.close();
+      // Create an empty results file to validate restart handling on minimal data.
       green::h5pp::archive ar(res_file_2, "w");
       ar["test"] << 1;
+      ar.set_attribute<std::string>("__grids_version__", new_grid_file_version);
       ar.close();
     }
     {
+      // Restart from the empty file and ensure the loop can be reinitialized
       p2["itermax"] = 4;
       g2            = p2["x0"];
       sigma1_2      = 0;
@@ -265,6 +299,61 @@ TEST_CASE("Self-consistency") {
       sc.solve(solver, h0, ovlp, g2, sigma1_2, sigma_tau_2);
     }
 
+    /**
+     * CASE 3: Restart with an old, incompatible grid file
+     *          Should throw an error due to outdated grid file (no version or version < 0.2.4)
+     */
+    // CASE: Restart with an old, incompatible grid file should throw an error
+    {
+      auto        p3         = green::params::params("DESCR");
+      std::string args_3     = "test --restart 1 --itermax 1 --E_thr 1e-13 --grid_file " + GRID_FILE_OLD
+                               + " --BETA 100 --results_file=" + res_file_2;
+      green::sc::define_parameters(p3);
+      green::grids::define_parameters(p3);
+      p3.define<double>("alpha", "", 0.45);
+      p3.define<double>("beta", "", 0.5);
+      p3.define<double>("gamma", "", 0.25);
+      p3.define<double>("x0", "", 0.2);
+      p3.parse(args_3);
+      REQUIRE_THROWS_AS(green::sc::sc_loop<fourth_power_equation_dyson>(MPI_COMM_WORLD, p3), green::grids::outdated_grids_file_error);
+    }
+
+    /**
+     * CASE 4: Restart with a results file created with an old grid file, while using a new grid file in the current run
+     *         Should throw an error due to outdated grid in the results file
+     */
+    std::filesystem::remove(res_file_2);
+    {
+      auto        p4         = green::params::params("DESCR");
+      std::string args_4    = "test --restart 1 --itermax 1 --E_thr 1e-13 --grid_file " + GRID_FILE_OLD
+                              + " --BETA 100 --results_file=" + res_file_2;
+      green::sc::define_parameters(p4);
+      green::grids::define_parameters(p4);
+      p4.define<double>("alpha", "", 0.45);
+      p4.define<double>("beta", "", 0.5);
+      p4.define<double>("gamma", "", 0.25);
+      p4.define<double>("x0", "", 0.2);
+      p4.parse(args_4);
+      // Fresh start (even though --restart 1) using old grid file in p4
+      green::sc::sc_loop<fourth_power_equation_dyson> sc(MPI_COMM_WORLD, p4);
+      fourth_power_equation_solver                    solver(p4["alpha"], p4["beta"]);
+      REQUIRE_NOTHROW(sc.solve(solver, h0, ovlp, g2, sigma1_2, sigma_tau_2));
+    }
+    {
+      auto        p5         = green::params::params("DESCR");
+      green::sc::define_parameters(p5);
+      green::grids::define_parameters(p5);
+      p5.define<double>("alpha", "", 0.45);
+      p5.define<double>("beta", "", 0.5);
+      p5.define<double>("gamma", "", 0.25);
+      p5.define<double>("x0", "", 0.2);
+      p5.parse(args_2);
+      // Initialize p5 from args_2 (using the new grid file) and attempt to restart from res_file_2,
+      // which was created with an older grid file; this should trigger outdated_results_file_error.
+      REQUIRE_THROWS_AS(green::sc::sc_loop<fourth_power_equation_dyson>(MPI_COMM_WORLD, p5), green::sc::outdated_results_file_error);
+    }
+
+    // Clean up
     std::filesystem::remove(res_file_1);
     std::filesystem::remove(res_file_2);
   }