Electricpotential analysis (#379)

Analysis for sampling the electric potential at arbitrary points. Related to collaboration with Lue and Curtis.

* Add externalpotential analysis

* Add non-overlapping random walk policy

* Manual update

Author: mlund

docs/_docs/analysis.md

@@ -395,6 +395,45 @@ atomic species can be saved.
 `pqrfile`            | Output PQR file (optional)
 `verbose=True`       | If `True`, add results to general output
 
+### Electric Potential
+
+`electricpotential`   | Description
+--------------------- | ------------------------------------------------------
+`nstep`               | Interval between samples
+`nskip=0`             | Number of initial steps excluded from the analysis
+`epsr`                | Dielectric constant
+`type`                | Coulomb type, `plain` etc. -- see energies
+`structure`           | Either a _filename_ (pqr, aam, gro etc) or a _list_ of positions
+`policy=fixed`        | Policy used to augment positions before each sample event, see below
+`ncalc`               | Number of potential calculations per sample event
+`stride`              | Separation between target points when using `random_walk` or `no_overlap`
+
+This calculates the mean electric potential, $\langle \phi\_i \rangle$ and correlations, $\langle \phi\_1\phi\_2 ...\rangle$
+at an arbitrary number of target positions in the simulation cell.
+The positions, given via `structure`, can be augmented using a `policy`:
+
+`policy`        | Description
+--------------- | ------------------------------------------------------------------
+`fixed`         | Expects a list of fixed positions where the potential is measured
+`random_walk`   | Assign random position to first target; while following targets are randomly placed `stride` distance from previous.
+`no_overlap`    | As `random_walk` but with no particle overlap (size defined by `sigma`, see Topology)
+
+Histograms of the correlation and the potentials at the target points are saved to disk.
+
+Example:
+
+~~~ yaml
+- electricpotential:
+    nstep: 20
+    ncalc: 10
+    epsr: 80
+    type: plain
+    policy: random_walk
+    stride: 10   # in angstrom
+    structure:
+      - [0,0,0]  # defines two target points...
+      - [0,0,0]  # ...positions are randomly set
+~~~
 
 ## Reaction Coordinate
 

docs/schema.yml

@@ -905,6 +905,36 @@ properties:
                         dr: {type: number, description: "Distance resolution along R (Å)"}
                     required: [molecules, nstep, dr]
 
+                electricpotential:
+                    type: object
+                    properties:
+                        epsr: {type: number, description: "Relative dielectric constant"}
+                        type: {type: string, description: "Coulomb potential type"}
+                        stride: {type: number, description: "Stride length for random walk"}
+                        policy:
+                            description: "Policy used to augment positions before each sample event"
+                            type: string
+                            enum: ["fixed", "random_walk", "no_overlap"]
+                            default: "fixed"
+                        ncalc: {type: integer, description: "Number of samples per event", default: 1, minimum: 1}
+                        nstep: {type: integer, description: "Sample interval"}
+                        nskip: {type: integer, default: 0, description: "Initial steps to skip"}
+                        structure:
+                            description: Target position definitions
+                            anyOf:
+                                - type: string
+                                  pattern: "(.*?)\\.(pqr|aam|xyz|gro)$"
+                                - type: array
+                                  items:
+                                      type: array
+                                      additionalProperties:
+                                          type: array
+                                          items: {type: number}
+                                          minItems: 3
+                                          maxItems: 3
+                    required: [nstep, epsr, type, structure]
+                    additionalProperties: false
+
                 chargefluctuations:
                     type: object
                     properties:

src/CMakeLists.txt

@@ -68,7 +68,7 @@ set(hdrs analysis.h average.h atomdata.h auxiliary.h bonds.h celllist.h celllist
 	aux/eigen_cerealisation.h aux/eigensupport.h aux/iteratorsupport.h aux/multimatrix.h
 	aux/eigen_cerealisation.h aux/eigensupport.h aux/equidistant_table.h aux/error_function.h
 	aux/exp_function.h aux/invsqrt_function.h aux/iteratorsupport.h aux/legendre.h aux/multimatrix.h
-	aux/pairmatrix.h aux/pow_function.h aux/table_1d.h aux/table_2d.h aux/timers.h aux/typeerasure.h)
+	aux/pairmatrix.h aux/pow_function.h aux/table_1d.h aux/table_2d.h aux/timers.h aux/typeerasure.h aux/sparsehistogram.h)
 
 set_source_files_properties(${objs} PROPERTIES LANGUAGE CXX)
 

src/analysis.cpp

@@ -3,6 +3,7 @@
 #include "energy.h"
 #include "reactioncoordinate.h"
 #include "multipole.h"
+#include "potentials.h"
 #include "aux/iteratorsupport.h"
 #include "aux/eigensupport.h"
 #include <range/v3/numeric/accumulate.hpp>
@@ -11,6 +12,7 @@
 #include <zstr.hpp>
 #include <cereal/types/memory.hpp>
 #include <cereal/archives/binary.hpp>
+#include <range/v3/numeric/accumulate.hpp>
 
 #include <iomanip>
 #include <iostream>
@@ -114,6 +116,8 @@ std::shared_ptr<Analysisbase> createAnalysis(const std::string& name, const json
             return std::make_shared<AtomDipDipCorr>(j, spc);
         } else if (name == "density") {
             return std::make_shared<Density>(j, spc);
+        } else if (name == "electricpotential") {
+            return std::make_shared<ElectricPotential>(j, spc);
         } else if (name == "chargefluctuations") {
             return std::make_shared<ChargeFluctuations>(j, spc);
         } else if (name == "molrdf") {
@@ -1776,4 +1780,145 @@ void SpaceTrajectory::_sample() {
 void SpaceTrajectory::_to_json(json& j) const { j = {{"file", filename}}; }
 
 void SpaceTrajectory::_to_disk() { stream->flush(); }
+
+// -----------------------------
+
+ElectricPotential::ElectricPotential(const json& j, Space& spc)
+    : Analysisbase(spc, "electricpotential"), potential_correlation_histogram(histogram_resolution) {
+    from_json(j);
+    coulomb = std::make_shared<Potential::NewCoulombGalore>();
+    coulomb->from_json(j);
+    getTargets(j);
+    setPolicy(j);
+    calculations_per_sample_event = j.value("ncalc", 1);
+}
+
+void ElectricPotential::setPolicy(const json& j) {
+    output_information.clear();
+    policy = j.value("policy", FIXED);
+    auto stride = 0.0;
+    switch (policy) {
+    case FIXED:
+        applyPolicy = []() {};
+        break;
+    case RANDOM_WALK:
+        stride = j.at("stride").get<double>();
+        output_information["stride"] = stride;
+        applyPolicy = [&, stride] {
+            auto origin = targets.begin();
+            spc.geo.randompos(origin->position, random);
+            std::for_each(std::next(origin), targets.end(), [&](Target& target) {
+                target.position = origin->position + stride * ranunit(random);
+                std::advance(origin, 1);
+            });
+        };
+        break;
+    case RANDOM_WALK_NO_OVERLAP:
+        stride = j.at("stride").get<double>();
+        output_information["stride"] = stride;
+        applyPolicy = [&, stride] {
+            auto origin = targets.begin();
+            do {
+                spc.geo.randompos(origin->position, random);
+            } while (overlapWithParticles(origin->position));
+            std::for_each(std::next(origin), targets.end(), [&](Target& target) {
+                do {
+                    target.position = origin->position + stride * ranunit(random);
+                } while (overlapWithParticles(target.position));
+                std::advance(origin, 1);
+            });
+        };
+        break;
+    default:
+        throw ConfigurationError("unknown policy");
+    }
+}
+
+void ElectricPotential::getTargets(const json& j) {
+    if (const auto& structure = j.find("structure"); structure == j.end()) {
+        throw ConfigurationError("missing structure");
+    } else {
+        PointVector positions;
+        if (structure->is_string()) { // load positions from chemical structure file
+            auto particles = loadStructure(structure->get<std::string>(), false);
+            positions = particles | ranges::cpp20::views::transform(&Particle::pos) | ranges::to_vector;
+        } else if (structure->is_array()) { // list of positions
+            positions = structure->get<PointVector>();
+        }
+        std::transform(positions.begin(), positions.end(), std::back_inserter(targets), [&](auto& position) {
+            Target target;
+            target.position = position;
+            target.potential_histogram = std::make_shared<SparseHistogram<double>>(histogram_resolution);
+            return target;
+        });
+        if (targets.empty()) {
+            throw ConfigurationError("no targets defined");
+        }
+    }
+}
+
+void ElectricPotential::_sample() {
+    for (unsigned int i = 0; i < calculations_per_sample_event; i++) {
+        applyPolicy();
+        auto potential_correlation = 1.0; // phi1 * phi2 * ...
+        for (auto& target : targets) {    // loop over each target point
+            auto potential = calcPotentialOnTarget(target);
+            target.potential_histogram->add(potential);
+            target.mean_potential += potential;
+            potential_correlation *= potential;
+        }
+        mean_potential_correlation += potential_correlation;        // <phi1 * phi2 * ...>
+        potential_correlation_histogram.add(potential_correlation); // P(<phi1 * phi2 * ...>)
+    }
+}
+
+double ElectricPotential::calcPotentialOnTarget(const ElectricPotential::Target& target) {
+    auto potential_from_particle = [&](const Particle& particle) {
+        const auto distance_to_target = std::sqrt(spc.geo.sqdist(particle.pos, target.position));
+        return coulomb->getCoulombGalore().ion_potential(particle.charge, distance_to_target);
+    };
+    auto potentials = spc.activeParticles() | ranges::cpp20::views::transform(potential_from_particle);
+    return std::accumulate(potentials.begin(), potentials.end(), 0.0);
+}
+
+void ElectricPotential::_to_json(json& j) const {
+    j = output_information;
+    coulomb->to_json(j["coulomb"]);
+    j["policy"] = policy;
+    j["number of targets"] = targets.size();
+    j["calculations per sample"] = calculations_per_sample_event;
+    if (number_of_samples > 0) {
+        j["correlation (βe)ⁱ⟨ϕ₁ϕ₂...ϕᵢ⟩"] = mean_potential_correlation.avg();
+        auto& mean_potential_j = j["mean potentials βe⟨ϕᵢ⟩"] = json::array();
+        std::transform(targets.begin(), targets.end(), std::back_inserter(mean_potential_j),
+                       [](const auto& target) { return target.mean_potential.avg(); });
+    }
+}
+void ElectricPotential::_to_disk() {
+    if (auto stream = std::ofstream(MPI::prefix + "potential_correlation_histogram.dat")) {
+        stream << potential_correlation_histogram;
+    }
+    int filenumber = 1;
+    for (const auto& target : targets) {
+        if (auto stream = std::ofstream(fmt::format("{}potential_histogram{}.dat", MPI::prefix, filenumber++))) {
+            stream << *(target.potential_histogram);
+        }
+    }
+}
+
+/**
+ * Checks if position lies within the spheres of diameters `sigma` defined
+ * by each active particle in the system. Complexity: N
+ *
+ * @return True if overlap with any particle
+ */
+bool ElectricPotential::overlapWithParticles(const Point& position) const {
+    auto overlap = [&position, &geometry = spc.geo](const Particle& particle) {
+        const auto radius = 0.5 * particle.traits().sigma;
+        return geometry.sqdist(particle.pos, position) < radius * radius;
+    };
+    auto particles = spc.activeParticles();
+    return std::any_of(particles.begin(), particles.end(), overlap);
+}
+
 } // namespace Faunus::Analysis

src/analysis.h

@@ -7,19 +7,23 @@
 #include "aux/timers.h"
 #include "aux/table_2d.h"
 #include "aux/equidistant_table.h"
+#include <aux/sparsehistogram.h>
 #include <set>
 
 namespace cereal {
 class BinaryOutputArchive;
 }
 
-namespace Faunus {
+namespace Faunus::Potential {
+class NewCoulombGalore;
+}
 
-namespace Energy {
+namespace Faunus::Energy {
 class Hamiltonian;
 class Energybase;
-} // namespace Energy
+} // namespace Faunus::Energy
 
+namespace Faunus {
 /**
  * Adding a new analysis requires the following steps:
  *
@@ -162,6 +166,48 @@ class WidomInsertion : public PerturbationAnalysisBase {
     WidomInsertion(const json& j, Space& spc, Energy::Hamiltonian& pot);
 };
 
+/**
+ * @brief Samples the electric potential at arbitrary positions in the simulation box
+ * @todo Add policies for random mass-center position and orientation
+ */
+class ElectricPotential : public Analysisbase {
+  public:
+    enum Policies { FIXED, RANDOM_WALK, RANDOM_WALK_NO_OVERLAP, INVALID };
+
+  private:
+    double histogram_resolution = 0.05; //!< Angstrom
+    unsigned int calculations_per_sample_event = 1;
+    struct Target {
+        Point position;                                               //!< Target position
+        Average<double> mean_potential;                               //!< mean potential at position
+        std::shared_ptr<SparseHistogram<double>> potential_histogram; //!< Histogram of observed potentials
+    };
+    std::vector<Target> targets;                             //!< List of target points where to sample the potential
+    Policies policy;                                         //!< Policy to apply to targets before each sample event
+    std::shared_ptr<Potential::NewCoulombGalore> coulomb;    //!< Class for calculating the potential
+    Average<double> mean_potential_correlation;              //!< Correlation between targets, <phi1 x phi2 x ... >
+    SparseHistogram<double> potential_correlation_histogram; //!< Distribution of correlations, P(<phi1 x phi2 x ... >)
+    void getTargets(const json& j);                          //!< Get user defined target positions
+    void setPolicy(const json& j);                           //!< Set user defined position setting policy
+    double calcPotentialOnTarget(const Target& target);      //!< Evaluate net potential of target position
+    bool overlapWithParticles(const Point& position) const;  //!< Check if position is within the radius of any particle
+    std::function<void()> applyPolicy;                       //!< Lambda for position setting policy
+    json output_information;                                 //!< json output generated during construction
+    void _to_json(json& j) const override;
+    void _sample() override;
+    void _to_disk() override;
+
+  public:
+    ElectricPotential(const json& j, Space& spc);
+};
+
+NLOHMANN_JSON_SERIALIZE_ENUM(ElectricPotential::Policies, {
+                                                              {ElectricPotential::INVALID, nullptr},
+                                                              {ElectricPotential::FIXED, "fixed"},
+                                                              {ElectricPotential::RANDOM_WALK_NO_OVERLAP, "no_overlap"},
+                                                              {ElectricPotential::RANDOM_WALK, "random_walk"},
+                                                          })
+
 /**
  * @brief Density of atom along axis
  *

src/aux/sparsehistogram.h

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <map>
+#include <fstream>
+
+namespace Faunus {
+
+/**
+ * @brief Histogram for an arbitrary set of values using a sparse memory layout (map)
+ *
+ * Builds a histogram by binning given values to a specified resolution. Values are stored
+ * in a memory efficient map-structure with log(N) lookup complexity.
+ */
+template <typename T = double> class SparseHistogram {
+    T resolution;
+    std::map<int, unsigned int> data;
+
+  public:
+    explicit SparseHistogram(T resolution) : resolution(resolution) {}
+    void add(const T value) {
+        if (std::isfinite(value)) {
+            data[static_cast<int>(std::round(value / resolution))]++;
+        } else {
+            faunus_logger->warn("histogram: skipping inf/nan number");
+        }
+    }
+    friend auto& operator<<(std::ostream& stream, const SparseHistogram& histogram) {
+        std::for_each(histogram.data.begin(), histogram.data.end(), [&](const auto& sample) {
+            stream << fmt::format("{:.6E} {}\n", T(sample.first) * histogram.resolution, sample.second);
+        });
+        return stream;
+    }
+};
+} // namespace Faunus