Analysis to disk (#227)

* Move file writing from destructor to _to_disk() throughout Analysis
namespace. Save to disk every macro-step. Minor check for atom array
when reading molecule types

* Add templated inv_sqrt() function and test; clarify manual

* Add ideal cluster test and update cluster documentation which did not
show properly on readthedocs

* Allow control of expensive debug check in virtual volume move and minor
code cleanup

* Fix json compilation error introduced in previous commit

Related to #196

docs/_docs/energy.md

@@ -557,21 +557,21 @@ customexternal:
 ## Bonded Interactions
 
 Bonds and angular potentials are added via the keyword `bondlist` either directly
-in a molecule definition (topology) or in energy/bonded where the latter can be
-used to add inter-molecular bonds:
+in a molecule definition (topology) for intra-molecular bonds, or in `energy->bonded`
+where the latter can be used to add inter-molecular bonds:
 
 ~~~ yaml
-energy:
-    - bonded:
-        bondlist: # absolute index
-           - harmonic: { index: [56,921], k: 10, req: 15 }
 moleculelist:
     - water: # TIP3P
         structure: "water.xyz"
         bondlist: # index relative to molecule
             - harmonic: { index: [0,1], k: 5024, req: 0.9572 }
             - harmonic: { index: [0,2], k: 5024, req: 0.9572 }
             - harmonic_torsion: { index: [1,0,2], k: 628, aeq: 104.52 }
+energy:
+    - bonded:
+        bondlist: # absolute index; can be between molecules
+           - harmonic: { index: [56,921], k: 10, req: 15 }
 ~~~
 
 Bonded potential types:
@@ -588,7 +588,7 @@ $\mu V T$ ensembles and Widom insertion are currently unsupported for molecules
 `index`        | Array with _exactly two_ indices (relative to molecule)
 
 $$
-u(r) = \frac{1}{2}k(r-r_{\mathrm{eq}})^2
+u(r) = \frac{1}{2}k(r-r\_{\mathrm{eq}})^2
 $$
 
 ### Finite Extensible Nonlinear Elastic
@@ -604,8 +604,8 @@ Finite extensible nonlinear elastic potential long range repulsive potential.
 $$
      u(r) =
   \begin{cases} 
-   -\frac{1}{2} k r_{\mathrm{max}}^2 \ln \left [ 1-(r/r_{\mathrm{max}})^2 \right ],       & \text{if } r < r_{\mathrm{max}} \\
-   \infty, & \text{if } r \geq r_{\mathrm{max}}
+   -\frac{1}{2} k r\_{\mathrm{max}}^2 \ln \left [ 1-(r/r\_{\mathrm{max}})^2 \right ],       & \text{if } r < r\_{\mathrm{max}} \\
+   \infty, & \text{if } r \geq r\_{\mathrm{max}}
   \end{cases}
 $$
 

docs/_docs/moves.md

@@ -86,7 +86,7 @@ Atomic _rotation_ affects only anisotropic particles such as dipoles, spherocyli
 ### Cluster Move
 
 `cluster`       | Description
---------------– | -----------------------
+--------------- | --------------------------------------------
 `molecules`     | Array of molecule names; `[*]` selects all
 `threshold`     | Mass-center threshold for forming a cluster
 `dir=[1,1,1]`   | Directions to translate

examples/cluster/cluster-ideal.out.json

@@ -0,0 +1,28 @@
+{
+    "moves": [
+        {
+            "moltransrot": {
+                "acceptance": 1.0,
+                "√⟨r²⟩": 11.5
+            }
+        },
+        {
+            "moltransrot": {
+                "acceptance": 1.0,
+                "√⟨r²⟩": 11.5
+            }
+        },
+        {
+            "cluster": {
+                "acceptance": 0.85,
+                "bias rejection rate": 0.15,
+                "satellites": [
+                    "other"
+                ],
+                "√⟨r²⟩": 10.4,
+                "√⟨θ²⟩/°": 45.7,
+                "⟨N⟩": 1.22
+            }
+        }
+    ]
+}

examples/cluster/cluster-ideal.yml

@@ -0,0 +1,36 @@
+#!/usr/bin/env yason.py
+#
+# In an ideal system the average cluster size equals
+#
+#     N_c = 1 + 4 * pi * Rc**3 * N / (3 * L**3) = 1.22 in this case
+#
+temperature: 300
+random: {seed: fixed}
+geometry: {type: cuboid, length: 40}
+mcloop: {macro: 10, micro: 20000}
+
+atomlist:
+    - OW: {q: -0.8476, sigma: 3.166, eps: 0.650, mw: 15.999}
+    - HW: {q: 0.4238,  sigma: 2, eps: 0, mw: 1.007}
+
+moleculelist:
+    - water: {structure: water.xyz, rigid: true}
+    - other: {structure: water.xyz, rigid: true}
+
+insertmolecules:
+    - water: {N: 5}
+    - other: {N: 5}
+
+energy:
+    - bonded: {} # just a dummy energy
+
+moves:
+    - moltransrot: {molecule: water, dp: 20, dprot: 1, repeat: 1}
+    - moltransrot: {molecule: other, dp: 20, dprot: 1, repeat: 1}
+    - cluster: {molecules: [water, other], dp: 20, dprot: 3, threshold: 7, satellites: [other]}
+
+analysis:
+    - savestate: {file: confout.pqr}
+    - savestate: {file: state.json}
+    - sanity: {nstep: 100}
+

src/CMakeLists.txt

@@ -266,6 +266,13 @@ add_test(
     && ${PYTHON_EXECUTABLE} ${JSON_COMPARE} cluster.out.json out.json --tol 0.1"
     WORKING_DIRECTORY ${EXAMPLES_DIR}/cluster)
 
+add_test(
+    NAME cluster-ideal
+    COMMAND sh -c "${PYTHON_EXECUTABLE} ${YASON} cluster-ideal.yml\
+    | $<TARGET_FILE:faunus> --quiet\
+    && ${PYTHON_EXECUTABLE} ${JSON_COMPARE} cluster-ideal.out.json out.json --tol 0.01"
+    WORKING_DIRECTORY ${EXAMPLES_DIR}/cluster)
+
 add_test(
     NAME bulk
     COMMAND sh -c "${PYTHON_EXECUTABLE} ${YASON} bulk.yml\

src/analysis.cpp

@@ -135,12 +135,16 @@ SystemEnergy::SystemEnergy(const json &j, Energy::Hamiltonian &pot) {
     auto u = energyFunc();
     uinit = std::accumulate(u.begin(), u.end(), 0.0); // initial energy
 }
+void SystemEnergy::_to_disk() {
+    if (f)
+        f.flush(); // empty buffer
+}
 
 void SaveState::_to_json(json &j) const { j["file"] = file; }
 
 void SaveState::_sample() { writeFunc(file); }
 
-SaveState::~SaveState() {
+void SaveState::_to_disk() {
     if (steps == -1)
         _sample();
 }
@@ -204,26 +208,6 @@ SaveState::SaveState(const json &j, Space &spc) {
 
 PairFunctionBase::PairFunctionBase(const json &j) { from_json(j); }
 
-PairFunctionBase::~PairFunctionBase() {
-    std::ofstream f(MPI::prefix + file);
-    if (f) {
-        double Vr = 1, sum = hist.sumy();
-        hist.stream_decorator = [&](std::ostream &o, double r, double N) {
-            if (dim == 3)
-                Vr = 4 * pc::pi * std::pow(r, 2) * dr;
-            else if (dim == 2) {
-                Vr = 2 * pc::pi * r * dr;
-                if (Rhypersphere > 0)
-                    Vr = 2.0 * pc::pi * Rhypersphere * std::sin(r / Rhypersphere) * dr;
-            } else if (dim == 1)
-                Vr = dr;
-            if (Vr > 0)
-                o << r << " " << N * V / (Vr * sum) << "\n";
-        };
-        f << hist;
-    }
-}
-
 void PairFunctionBase::_to_json(json &j) const {
     j = {{"dr", dr / 1.0_angstrom}, {"name1", name1},        {"name2", name2}, {"file", file}, {"dim", dim},
          {"slicedir", slicedir},    {"thickness", thickness}};
@@ -243,10 +227,29 @@ void PairFunctionBase::_from_json(const json &j) {
     hist.setResolution(dr, 0);
     Rhypersphere = j.value("Rhyper", -1.0);
 }
+void PairFunctionBase::_to_disk() {
+    std::ofstream f(MPI::prefix + file);
+    if (f) {
+        double Vr = 1, sum = hist.sumy();
+        hist.stream_decorator = [&](std::ostream &o, double r, double N) {
+            if (dim == 3)
+                Vr = 4 * pc::pi * std::pow(r, 2) * dr;
+            else if (dim == 2) {
+                Vr = 2 * pc::pi * r * dr;
+                if (Rhypersphere > 0)
+                    Vr = 2.0 * pc::pi * Rhypersphere * std::sin(r / Rhypersphere) * dr;
+            } else if (dim == 1)
+                Vr = dr;
+            if (Vr > 0)
+                o << r << " " << N * V / (Vr * sum) << "\n";
+        };
+        f << hist;
+    }
+}
 
 PairAngleFunctionBase::PairAngleFunctionBase(const json &j) : PairFunctionBase(j) { from_json(j); }
 
-PairAngleFunctionBase::~PairAngleFunctionBase() {
+void PairAngleFunctionBase::_to_disk() {
     std::ofstream f(MPI::prefix + file);
     if (f) {
         hist2.stream_decorator = [&](std::ostream &o, double r, double N) {
@@ -268,18 +271,22 @@ void VirtualVolume::_sample() {
         scaleVolume(Vold);                               // restore saved system
 
         double du = Unew - Uold;          // system energy change
-        if (-du < pc::max_exp_argument) { // does minus energy change fits exp() function?
+        if (-du < pc::max_exp_argument) { // does minus energy change fit exp() function?
             double exp_du = std::exp(-du);
             assert(not std::isnan(exp_du));
             duexp += exp_du; // collect average, <exp(-du)>
-            if (f)           // write sample event to output file
-                f << cnt << " " << dV << " " << du << " " << exp_du << " " << std::log(duexp.avg()) / dV << "\n";
-#ifndef NDEBUG
-            if (Uold != 0) { // check if volume and particle positions are properly restored
-                double should_be_small = std::fabs((Uold - pot.energy(c)) / Uold);
-                assert(should_be_small < 1e-4); // expensive
+            if (output_file) // write sample event to output file
+                output_file << cnt << " " << dV << " " << du << " " << exp_du << " " << std::log(duexp.avg()) / dV
+                            << "\n";
+
+            // Check if volume and particle positions are properly restored.
+            // Expensive and one would normally not perform this test and we trigger it
+            // only when using log-level "debug" or lower
+            if (faunus_logger->level() <= spdlog::level::debug and Uold != 0) {
+                double should_be_small = std::fabs((Uold - pot.energy(c)) / Uold); // expensive!
+                if (should_be_small > 1e-6)
+                    faunus_logger->error("{} failed to restore system", name);
             }
-#endif
         } else {   // energy change too large (negative) to fit exp() function
             cnt--; // cnt is incremented by sample() so we need to decrease
             faunus_logger->warn("{0}: skipping sample event due to excessive energy, dU/kT={1}", name, du);
@@ -291,14 +298,15 @@ void VirtualVolume::_from_json(const json &j) {
     dV = j.at("dV");
     file = MPI::prefix + j.value("file", std::string());
     if (not file.empty()) { // if filename is given, create output file
-        if (f)
-            f.close();
-        f.open(file);
-        if (!f)
+        if (output_file)
+            output_file.close();
+        output_file.open(file);
+        if (!output_file)
             throw std::runtime_error(name + ": cannot open output file " + file);
-        f << "# steps dV/" + u8::angstrom + u8::cubed + " du/kT exp(-du/kT) <Pex>/kT/" + u8::angstrom + u8::cubed
-          << "\n"; // file header
-        f.precision(14);
+        output_file << "# steps dV/" + u8::angstrom + u8::cubed + " du/kT exp(-du/kT) <Pex>/kT/" + u8::angstrom +
+                           u8::cubed
+                    << "\n"; // file header
+        output_file.precision(14);
     }
 }
 
@@ -316,6 +324,10 @@ VirtualVolume::VirtualVolume(const json &j, Space &spc, Energy::Energybase &pot)
     getVolume = [&spc]() { return spc.geo.getVolume(); };
     scaleVolume = [&spc](double Vnew) { spc.scaleVolume(Vnew); };
 }
+void VirtualVolume::_to_disk() {
+    if (output_file)
+        output_file.flush(); // empty buffer
+}
 
 void QRtraj::_sample() { write_to_file(); }
 
@@ -340,14 +352,17 @@ QRtraj::QRtraj(const json &j, Space &spc) {
         f << "\n";               // newline for every frame
     };
 }
+void QRtraj::_to_disk() {
+    if (f)
+        f.flush(); // empty buffer
+}
 
 void CombinedAnalysis::sample() {
     for (auto &ptr : this->vec)
         ptr->sample();
 }
 
-CombinedAnalysis::~CombinedAnalysis() {
-    // this is really a hack; the destructor should not be in charge of this
+void CombinedAnalysis::to_disk() {
     for (auto &ptr : this->vec)
         ptr->to_disk();
 }
@@ -446,6 +461,10 @@ FileReactionCoordinate::FileReactionCoordinate(const json &j, Space &spc) {
     type = j.at("type").get<std::string>();
     rc = ReactionCoordinate::createReactionCoordinate({{type, j}}, spc);
 }
+void FileReactionCoordinate::_to_disk() {
+    if (file)
+        file.flush(); // empty buffer
+}
 
 void WidomInsertion::_sample() {
     if (!change.empty()) {
@@ -603,7 +622,7 @@ Density::Density(const json &j, Space &spc) : spc(spc) {
         }
     }
 }
-Density::~Density() {
+void Density::_to_disk() {
     for (auto &m : atmdhist) { // atomic molecules
         std::string file = "rho-"s + molecules.at(m.first).name + ".dat";
         std::ofstream f(MPI::prefix + file);
@@ -657,6 +676,7 @@ Density::~Density() {
         }
     }
 }
+
 void SanityCheck::_sample() {
     // loop over all groups
     for (auto &g : spc.groups) {
@@ -905,7 +925,7 @@ MultipoleDistribution::MultipoleDistribution(const json &j, Space &spc) : spc(sp
         throw std::runtime_error("specify exactly two molecules");
 }
 
-MultipoleDistribution::~MultipoleDistribution() { save(); }
+void MultipoleDistribution::_to_disk() { save(); }
 
 // =============== AtomInertia ===============
 
@@ -930,6 +950,10 @@ AtomInertia::AtomInertia(const json &j, Space &spc) : spc(spc) {
     file.open(filename); // output file
     index = j.at("index").get<size_t>(); // atom id
 }
+void AtomInertia::_to_disk() {
+    if (file)
+        file.flush(); // empty buffer
+}
 
 // =============== InertiaTensor ===============
 
@@ -961,6 +985,10 @@ InertiaTensor::InertiaTensor(const json &j, Space &spc) : spc(spc) {
     index = j.at("index").get<size_t>(); // group index
     indexes = j.value("indexes", std::vector<size_t>({0, spc.groups[index].size()})); // whole molecule by default
 }
+void InertiaTensor::_to_disk() {
+    if (file)
+        file.flush(); // empty buffer
+}
 
 // =============== MultipoleMoments ===============
 
@@ -1011,6 +1039,10 @@ MultipoleMoments::MultipoleMoments(const json &j, Space &spc) : spc(spc) {
     indexes = j.value("indexes", std::vector<size_t>({0, spc.groups[index].size()-1})); // whole molecule by default
     mol_cm = j.value("mol_cm", true); // use the mass center of the whole molecule
 }
+void MultipoleMoments::_to_disk() {
+    if (file)
+        file.flush(); // empty buffer
+}
 
 // =============== PolymerShape ===============
 
@@ -1168,7 +1200,7 @@ SlicedDensity::SlicedDensity(const json &j, Space &spc) : spc(spc) {
     name = "sliceddensity";
     from_json(j);
 }
-SlicedDensity::~SlicedDensity() {
+void SlicedDensity::_to_disk() {
     std::ofstream f(MPI::prefix + file);
     if (f and cnt > 0) {
         f << "# z rho/M\n";