Add on-the-fly swapping between gillespie and standard dynamics

src/main.cpp

@@ -70,8 +70,11 @@ int main(int argc, char *argv[])
             "exit rates at once and flips a spin every iteration of the "
             "algorithm, but with a waiting time not necessarily equal to 1. The "
             "auto selection runs quick simulations of both types to see which "
-            "is faster, and selects that one."
-        )->check(CLI::IsMember({"standard", "gillespie", "auto"}));
+            "is faster, and selects that one. Finally, the dynamic option "
+            "switches between standard and gillespie dynamics on the fly depending "
+            "on the energy of the current tracer: when it is below the threshold "
+            "energy, we use gillespie; above, we use standard."
+        )->check(CLI::IsMember({"standard", "gillespie", "auto", "dynamic"}));
 
         app.add_option(
             "-n, --n_tracers", p.n_tracers,

src/spin.cpp

@@ -34,6 +34,10 @@ void SpinSystem::_first_time_state_initialization_()
 
     if (params.dynamics == "standard"){_init_standard();}
     else if (params.dynamics == "gillespie"){_init_gillespie();}
+    else if (params.dynamics == "dynamic"){
+        _init_standard();
+        _init_gillespie();
+    }
     else
     {
         throw std::runtime_error("Uknown dynamics during setup");
@@ -255,6 +259,29 @@ double SpinSystem::step()
     {
         waiting_time = _step_gillespie();
     }
+    else if (params.dynamics == "dynamic")
+    {
+        const double current_energy = energy();
+        const double threshold_E = params.energetic_threshold;
+        const double threshold_A = params.entropic_attractor;
+        double threshold;
+        if (params.valid_entropic_attractor)
+        {
+            threshold = std::max(threshold_E, threshold_A);
+        }
+        else
+        {
+            threshold = threshold_E;
+        }
+        if (current_energy < threshold)
+        {
+            waiting_time = _step_gillespie();
+        }
+        else
+        {
+            waiting_time = _step_standard();
+        }
+    }
     else
     {
         throw std::runtime_error("Uknown dynamics during step");
@@ -270,6 +297,10 @@ SpinSystem::~SpinSystem()
 {
     if (params.dynamics == "standard"){_teardown_standard();}
     else if (params.dynamics == "gillespie"){_teardown_gillespie();}
+    else if (params.dynamics == "dynamic"){
+        _teardown_standard();
+        _teardown_gillespie();
+    }
 
     // Else both to be safe
     else{_teardown_standard(); _teardown_gillespie();}