From 7f1a0f5b45fbbb37d94e49059af07f6d033ae316 Mon Sep 17 00:00:00 2001
From: Mikael Lund <Mlund@mac.com>
Date: Tue, 12 Dec 2023 20:38:24 +0000
Subject: [PATCH] Save minimum energy configuration from `systemenergy`
 analysis (#440)

* Add `save_min_conf` to `systemenergy` analysis

This enables saving of the minimum observed energy configuration.

- [x] docs
- [x] json schema

* fix decimal

* Fix clang17 compilation (constexpr name)
---
 docs/_docs/analysis.md |  2 ++
 docs/schema.yml        |  4 ++++
 src/analysis.cpp       | 35 ++++++++++++++++++++++++++++++++---
 src/analysis.h         |  3 +++
 4 files changed, 41 insertions(+), 3 deletions(-)

diff --git a/docs/_docs/analysis.md b/docs/_docs/analysis.md
index 5752e17e2..34fd0d8bb 100644
--- a/docs/_docs/analysis.md
+++ b/docs/_docs/analysis.md
@@ -644,6 +644,8 @@ All units in $k\_BT$.
 -------------- | -------------------------------------------
 `file`         | Output filename (`.dat`, `.csv`, `.dat.gz`)
 `nstep`        | Interval between samples
+`nskip=0`      | Number of initial steps excluded from the analysis
+`save_min_conf=false` | Dump minimum energy configuration to `PQR` file
 
 
 ### Penalty function
diff --git a/docs/schema.yml b/docs/schema.yml
index bf922079c..74d2c7173 100644
--- a/docs/schema.yml
+++ b/docs/schema.yml
@@ -1206,6 +1206,10 @@ properties:
                         file: {type: string}
                         nstep: {type: integer}
                         nskip: {type: integer, default: 0, description: Initial steps to skip}
+                        save_min_conf:
+                            type: boolean
+                            default: false
+                            description: "Save minimum energy configuration to PQR file"
                     required: [file, nstep]
                     additionalProperties: false
                     type: object
diff --git a/src/analysis.cpp b/src/analysis.cpp
index 6afcde9eb..7807a6013 100644
--- a/src/analysis.cpp
+++ b/src/analysis.cpp
@@ -12,6 +12,7 @@
 #include "aux/eigensupport.h"
 #include "aux/arange.h"
 #include "aux/matrixmarket.h"
+#include <cmath>
 #include <exception>
 #include <iterator>
 #include <range/v3/numeric/accumulate.hpp>
@@ -241,12 +242,33 @@ void SystemEnergy::normalize() {
     }
 }
 
+/**
+ * @brief Checks if the current energy is the lowest so far and saves the configuration if so.
+ *
+ * The output is hardcoded to PQR format, tagged with the step number and energy.
+ *
+ * @return True if a new minimum energy was encountered
+ */
+bool SystemEnergy::updateMinimumEnergy(const double current_energy) {
+    if (!dump_minimum_energy_configuration || current_energy >= minimum_energy) {
+        return false;
+    }
+    minimum_energy = current_energy;
+    const auto filename = MPI::prefix + "mininum_energy.pqr";
+    faunus_logger->debug("{}: saving {} ({:.2f} kT) at step {}", name, filename, minimum_energy, getNumberOfSteps());
+    PQRWriter(PQRWriter::Style::PQR).save(filename, spc.groups, spc.geometry.getLength());
+    return true;
+}
+
 void SystemEnergy::_sample() {
     const auto energies = calculateEnergies(); // current energy from all terms in Hamiltonian
     const auto total_energy = ranges::accumulate(energies, 0.0);
+    updateMinimumEnergy(total_energy);
     if (std::isfinite(total_energy)) {
         mean_energy += total_energy;
         mean_squared_energy += total_energy * total_energy;
+    } else {
+        faunus_logger->warn("{}: non-finite energy excluded from averages at step {}", name, getNumberOfSteps());
     }
     *output_stream << fmt::format("{:10d}{}{:.6E}", getNumberOfSteps(), separator, total_energy);
     for (auto energy : energies) {
@@ -256,7 +278,10 @@ void SystemEnergy::_sample() {
 }
 
 void SystemEnergy::_to_json(json& j) const {
-    j = {{"file", file_name}, {"init", initial_energy}, {"final", calculateEnergies()}};
+    j = {{"file", file_name},
+         {"init", initial_energy},
+         {"final", calculateEnergies()},
+         {"save_min_conf", dump_minimum_energy_configuration}};
     if (!mean_energy.empty()) {
         j["mean"] = mean_energy.avg();
         j["Cv/kB"] = mean_squared_energy.avg() - std::pow(mean_energy.avg(), 2);
@@ -264,11 +289,15 @@ void SystemEnergy::_to_json(json& j) const {
     roundJSON(j, 5);
 }
 
-void SystemEnergy::_from_json(const json& j) { file_name = MPI::prefix + j.at("file").get<std::string>(); }
+void SystemEnergy::_from_json(const json& j) {
+    file_name = MPI::prefix + j.at("file").get<std::string>();
+    dump_minimum_energy_configuration = j.value("save_min_conf", false);
+}
 
 void SystemEnergy::createOutputStream() {
     output_stream = IO::openCompressedOutputStream(file_name, true);
-    if (const auto suffix = file_name.substr(file_name.find_last_of('.') + 1); suffix == "csv") {
+    const bool is_csv_file = file_name.substr(file_name.find_last_of('.') + 1) == "csv";
+    if (is_csv_file) {
         separator = ",";
     } else {
         separator = " ";
diff --git a/src/analysis.h b/src/analysis.h
index 8582efd00..39c036b26 100644
--- a/src/analysis.h
+++ b/src/analysis.h
@@ -510,7 +510,10 @@ class SystemEnergy : public Analysis {
     Average<double> mean_squared_energy;
     Table2D<double, double> energy_histogram; // Density histograms
     double initial_energy = 0.0;
+    double minimum_energy = std::numeric_limits<double>::infinity(); //!< Tracks minimum energy
+    bool dump_minimum_energy_configuration = false;                  //!< Dump minimum energy config to disk
 
+    bool updateMinimumEnergy(double current_energy);
     void createOutputStream();
     void normalize();
     void _sample() override;