Voronoi analysis using Voronota-LT (#444)

* Add voronotalt headers

* Add voronoi analysis

* Add test

* Make voronotalt header a system header

* Update documentation

* DOI to voronota-lt

* Detect and log PBC information

* Support PBC

Author: mlund

docs/_docs/analysis.md

@@ -364,6 +364,22 @@ is filed to disk, `sasa_histogram.dat`.
 `file`         | Optionally stream area for each `nstep` to file (`.dat|.dat.gz`)
 `radius=1.4`   | Probe radius (Å)
 
+### Voronoi Tessellation (experimental)
+
+Performs Voronoi tessellation to detect contacts and surface areas using
+the [Voronota-LT library](https://doi.org/10/mq8k).
+Only the total surface area is currently reported, but more options are planned for future
+releases.
+The algorithm is significantly faster than the above `sasa` analysis.
+
+`voronoi`      | Description
+-------------- | ---------------------------
+`nstep`        | Interval between samples
+`nskip=0`      | Number of initial steps excluded from the analysis
+`file`         | Optionally stream surface area for each `nstep` to file (`.dat|.dat.gz`)
+`radius=1.4`   | Probe radius (Å)
+
+
 ## Charge Properties
 
 ### Molecular Multipoles

docs/schema.yml

@@ -1263,6 +1263,23 @@ properties:
                         nskip: {type: integer, default: 0, description: Number of steps to initially skip}
                     required: [dL, molecule, nstep]
                     additionalProperties: false
+
+                voronoi:
+                    description: "Voronoi tessellation using Voronota"
+                    type: object
+                    properties:
+                        nstep: {type: integer, description: "Interval between samples"}
+                        nskip: {type: integer, default: 0, description: "Initial steps to skip"}
+                        radius:
+                            type: number
+                            default: 1.4
+                            description: "Probe radius (Å)"
+                        file:
+                            type: string
+                            pattern: "(.*?)\\.(dat|dat.gz)$"
+                            description: "Streaming filename (.dat|.dat.gz)"
+                    required: [nstep]
+                    additionalProperties: false
 
                 widom:
                     description: "Widom ghost particle insertion"

examples/sasa/sasa.out.json

@@ -6,12 +6,21 @@
         "nstep": 1,
         "radius": 2.0,
         "reference": "doi:10/dbjh",
-        "relative time": 0.592,
         "samples": 1,
         "slices_per_atom": 20,
         "⟨SASA²⟩-⟨SASA⟩²": 0.0,
         "⟨SASA⟩": 2177.145143881812
       }
+    },
+    {
+      "voronota": {
+        "nstep": 1,
+        "radius": 2.0,
+        "reference": "doi:10/mq8k",
+        "samples": 1,
+        "⟨SASA²⟩-⟨SASA⟩²": 0.0,
+        "⟨SASA⟩": 2186.5484868984963
+      }
     }
   ],
   "compiler": "Apple LLVM 15.0.0 (clang-1500.0.40.1)",

examples/sasa/sasa.yml

@@ -26,4 +26,8 @@ analysis:
         molecule: mymolecule
         radius: 2.0
         nstep: 1
+    - voronoi:
+        file: voronoi_sasa.dat
+        radius: 2.0
+        nstep: 1
 

src/CMakeLists.txt

@@ -8,6 +8,8 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 target_link_libraries(project_options INTERFACE ZLIB::ZLIB)
 target_link_libraries(project_options INTERFACE ${CMAKE_THREAD_LIBS_INIT})
 
+include_directories(SYSTEM ${CMAKE_SOURCE_DIR}/voronotalt)
+
 # ========== pre-compiler system headers ==========
 
 option(ENABLE_PCH "Enable Precompiled Headers" OFF)
@@ -60,7 +62,7 @@ set(objs actions.cpp analysis.cpp average.cpp atomdata.cpp auxiliary.cpp bonds.c
 	chainmove.cpp clustermove.cpp core.cpp forcemove.cpp units.cpp energy.cpp externalpotential.cpp
 	geometry.cpp group.cpp io.cpp molecule.cpp montecarlo.cpp move.cpp mpicontroller.cpp
 	particle.cpp penalty.cpp potentials.cpp random.cpp reactioncoordinate.cpp regions.cpp rotate.cpp sasa.cpp
-        scatter.cpp smart_montecarlo.cpp space.cpp speciation.cpp spherocylinder.cpp tensor.cpp)
+        scatter.cpp smart_montecarlo.cpp space.cpp speciation.cpp spherocylinder.cpp tensor.cpp voronota.cpp)
 
 set(hdrs actions.h analysis.h average.h atomdata.h auxiliary.h bonds.h celllist.h celllistimpl.h
 	chainmove.h clustermove.h core.h forcemove.h energy.h externalpotential.h geometry.h group.h io.h

src/analysis.cpp

@@ -25,7 +25,6 @@
 #include <cereal/archives/binary.hpp>
 #include <range/v3/numeric/accumulate.hpp>
 #include <range/v3/view/zip.hpp>
-
 #include <iomanip>
 #include <iostream>
 #include <memory>
@@ -194,6 +193,8 @@ std::unique_ptr<Analysis> createAnalysis(const std::string& name, const json& j,
             return std::make_unique<VirtualVolumeMove>(j, spc, pot);
         } else if (name == "virtualtranslate") {
             return std::make_unique<VirtualTranslate>(j, spc, pot);
+        } else if (name == "voronoi") {
+            return std::make_unique<Voronota>(j, spc);
         } else if (name == "widom") {
             return std::make_unique<WidomInsertion>(j, spc, pot);
         } else if (name == "xtcfile") {
@@ -1879,14 +1880,14 @@ SASAAnalysis::SASAAnalysis(const double probe_radius, const int slices_per_atom,
         break;
     }
     setPolicy(selected_policy);
+    cite = "doi:10/dbjh";
 }
 
 SASAAnalysis::SASAAnalysis(const json& j, const Space& spc)
     : SASAAnalysis(j.value("radius", 1.4_angstrom), j.value("slices", 20), j.value("resolution", 50.0),
                    j.value("policy", Policies::INVALID), spc) {
     from_json(j);
     policy->from_json(j);
-    cite = "doi:10/dbjh";
 }
 
 void SASAAnalysis::_to_json(json& json_ouput) const {

src/analysis.h

@@ -262,7 +262,7 @@ class ElectricPotential : public Analysis {
     unsigned int calculations_per_sample_event = 1;
     std::string file_prefix; //!< Output filename prefix for potential histogram and correlation
     struct Target {
-        Point position;                                               //!< Target position
+        Point position;                                                //!< Target position
         Average<double> mean_potential;                               //!< mean potential at position
         std::unique_ptr<SparseHistogram<double>> potential_histogram; //!< Histogram of observed potentials
     };
@@ -1086,6 +1086,36 @@ class SavePenaltyEnergy : public Analysis {
     SavePenaltyEnergy(const json& j, const Space& spc, const Energy::Hamiltonian& pot);
 };
 
+/**
+ * @brief Analysis of Vorononoi tessellation using the Voronota-LT library
+ *
+ * @todo Currently only SASA information is reported in the output. Add more information!
+ *
+ * https://doi.org/10/mq8k
+ */
+class Voronota : public Analysis {
+  private:
+    struct Averages {
+        Average<double> area;
+        Average<double> area_squared;
+    }; //!< Placeholder class for average properties
+    Averages average_data; //!< Stores all averages for the selected molecule
+
+    std::unique_ptr<std::ostream> output_stream; //!< output stream
+    double probe_radius;                         //!< radius of the probe sphere
+    std::string filename;                        //!< output file name
+    bool use_pbc = false;                        //!< Is the cell periodic?
+
+    void _to_json(json& j) const override;
+    void _from_json(const json& input) override;
+    void _to_disk() override;
+    void _sample() override;
+
+  public:
+    Voronota(const json& j, const Space& spc);
+    Voronota(double probe_radius, const Space& spc);
+};
+
 } // namespace analysis
 
 } // namespace Faunus

src/voronota.cpp

@@ -0,0 +1,85 @@
+#include "analysis.h"
+#include "mpicontroller.h"
+#include <voronotalt/voronotalt.h>
+#include <numeric>
+
+namespace Faunus::analysis {
+
+Voronota::Voronota(double probe_radius, const Faunus::Space& spc)
+    : Analysis(spc, "voronota")
+    , probe_radius(probe_radius) {
+    cite = "doi:10/mq8k";
+    auto n_pbc = spc.geometry.asSimpleGeometry()->boundary_conditions.isPeriodic().count();
+    switch (n_pbc) {
+    case 0:
+        faunus_logger->debug("{}: No PBC detected", name);
+        use_pbc = false;
+        break;
+    case 3:
+        faunus_logger->debug("{}: 3D PBC detected", name);
+        use_pbc = true;
+        break;
+    default:
+        faunus_logger->warn("{}: Non-uniform PBC is currently ignored - be careful!", name);
+        use_pbc = false;
+    }
+}
+
+Voronota::Voronota(const Faunus::json& input, const Faunus::Space& spc)
+    : Voronota(input.value("radius", 1.4_angstrom), spc) {
+    from_json(input);
+}
+
+void Voronota::_from_json(const Faunus::json& input) {
+    if (filename = input.value("file", ""s); !filename.empty()) {
+        output_stream = IO::openCompressedOutputStream(MPI::prefix + filename);
+        *output_stream << "# step SASA\n";
+    }
+}
+
+void Voronota::_to_json(json& json_ouput) const {
+    if (!average_data.area.empty()) {
+        json_ouput = {{"⟨SASA⟩", average_data.area.avg()},
+                      {"⟨SASA²⟩-⟨SASA⟩²", average_data.area_squared.avg() - std::pow(average_data.area.avg(), 2)}};
+    }
+    json_ouput["radius"] = probe_radius;
+}
+
+void Voronota::_to_disk() {
+    if (output_stream) {
+        output_stream->flush();
+    }
+}
+
+void Voronota::_sample() {
+    using voronotalt::SimplePoint;
+    using namespace ranges::cpp20::views;
+
+    // Convert single `Particle` to Voronota's `SimpleSphere`
+    auto to_sphere = [&](const Particle& p) -> voronotalt::SimpleSphere {
+        return {{p.pos.x(), p.pos.y(), p.pos.z()}, 0.5 * p.traits().sigma + probe_radius};
+    };
+    const auto spheres = spc.activeParticles() | transform(to_sphere) | ranges::to_vector;
+
+    voronotalt::RadicalTessellation::Result result;
+
+    if (use_pbc) {
+        auto to_point = [](const Point& p) -> SimplePoint { return {p.x(), p.y(), p.z()}; };
+        const auto corner = 0.5 * spc.geometry.getLength();
+        const std::vector<SimplePoint> box_corners = {to_point(-corner), to_point(corner)};
+        voronotalt::RadicalTessellation::construct_full_tessellation(spheres, box_corners, result);
+    } else {
+        voronotalt::RadicalTessellation::construct_full_tessellation(spheres, result);
+    }
+
+    auto areas = result.cells_summaries | transform([](auto& s) { return s.sas_area; });
+    const auto total_area = std::accumulate(areas.begin(), areas.end(), 0.0);
+    average_data.area.add(total_area);
+    average_data.area_squared.add(total_area * total_area);
+
+    if (output_stream) {
+        *output_stream << this->getNumberOfSteps() << " " << total_area << "\n";
+    }
+}
+
+} // namespace Faunus::analysis

src/voronotalt/LICENSE.txt

@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2023 Kliment Olechnovic
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.

src/voronotalt/VERSION.txt

@@ -0,0 +1 @@
+Voronota-LT version 0.9.2