Skip to content

Commit

Permalink
Save minimum energy configuration from systemenergy analysis (#440)
Browse files Browse the repository at this point in the history
* 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)
  • Loading branch information
mlund authored Dec 12, 2023
1 parent 7bb33d4 commit 7f1a0f5
Show file tree
Hide file tree
Showing 4 changed files with 41 additions and 3 deletions.
2 changes: 2 additions & 0 deletions docs/_docs/analysis.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions docs/schema.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
35 changes: 32 additions & 3 deletions src/analysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -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) {
Expand All @@ -256,19 +278,26 @@ 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);
}
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 = " ";
Expand Down
3 changes: 3 additions & 0 deletions src/analysis.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit 7f1a0f5

Please sign in to comment.