diff --git a/src/analysis.cpp b/src/analysis.cpp index 4c0a75774..fa80b4627 100644 --- a/src/analysis.cpp +++ b/src/analysis.cpp @@ -2263,13 +2263,26 @@ void ElectricPotential::setPolicy(const json& j) { stride = j.at("stride").get(); output_information["stride"] = stride; applyPolicy = [&, stride] { + unsigned int cnt = 0; auto origin = targets.begin(); do { spc.geometry.randompos(origin->position, random); + cnt++; + if (cnt > max_overlap_trials) { + faunus_logger->warn("{}: Max number of overlaps reached. Using overlapping point", name); + break; + } } while (overlapWithParticles(origin->position)); + + cnt = 0; std::for_each(std::next(origin), targets.end(), [&](Target& target) { do { target.position = origin->position + stride * randomUnitVector(random); + cnt++; + if (cnt > max_overlap_trials) { + faunus_logger->warn("{}: Max number of overlaps reached. Using overlapping point", name); + break; + } } while (overlapWithParticles(target.position)); std::advance(origin, 1); }); diff --git a/src/analysis.h b/src/analysis.h index 3e8032eb6..09f5ad1ef 100644 --- a/src/analysis.h +++ b/src/analysis.h @@ -228,7 +228,8 @@ class ElectricPotential : public Analysis { enum class Policies { FIXED, RANDOM_WALK, RANDOM_WALK_NO_OVERLAP, INVALID }; private: - double histogram_resolution = 0.05; //!< Angstrom + unsigned int max_overlap_trials = 100; //!< Maximum number of overlap checks before bailing out + double histogram_resolution = 0.01; //!< Potential resolution unsigned int calculations_per_sample_event = 1; struct Target { Point position; //!< Target position diff --git a/src/move.cpp b/src/move.cpp index 0c632d32d..6daa740b4 100644 --- a/src/move.cpp +++ b/src/move.cpp @@ -725,11 +725,11 @@ void ParallelTempering::exchangeState(Change& change) { void ParallelTempering::_move(Change& change) { mpi.world.barrier(); // wait until all ranks reach here - if (!MPI::checkRandomEngineState(mpi.world, Move::slump)) { + if (!MPI::checkRandomEngineState(mpi.world, slump)) { faunus_logger->error("Random numbers out of sync across MPI nodes. Do not use 'hardware' seed."); mpi.world.abort(1); // neighbor search *requires* that random engines are in sync } - partner->generate(mpi.world, Move::slump); + partner->generate(mpi.world, slump); if (partner->rank.has_value()) { exchangeState(change); } diff --git a/src/move.h b/src/move.h index f03ca06b2..7129fa737 100644 --- a/src/move.h +++ b/src/move.h @@ -532,7 +532,7 @@ class ParallelTempering : public Move { std::unique_ptr partner; //!< Policy for finding MPI partners Geometry::VolumeMethod volume_scaling_method = Geometry::VolumeMethod::ISOTROPIC; //!< How to scale volumes std::map> acceptance_map; //!< Exchange statistics - + Random slump; // static instance of Random (shared for all in ParallelTempering) void _to_json(json& j) const override; void _from_json(const json& j) override; void _move(Change& change) override; diff --git a/src/tabulate.h b/src/tabulate.h index 3dfa4709a..ea2cd8fc3 100644 --- a/src/tabulate.h +++ b/src/tabulate.h @@ -121,7 +121,7 @@ template class Andrea : public TabulatorBase * - `[0]==true`: tolerance is approved, * - `[1]==true` Repulsive part is found. */ - std::vector CheckUBuffer(std::vector &ubuft, T rlow, T rupp, std::function f) const { + std::vector CheckUBuffer(std::vector& ubuft, T rlow, T rupp, std::function f) const { // Number of points to control int ncheck = 11; @@ -161,7 +161,7 @@ template class Andrea : public TabulatorBase * @param r2 value * @note Auto-vectorization in Clang: https://llvm.org/docs/Vectorizers.html */ - inline T eval(const typename base::data &d, T r2) const { + inline T eval(const typename base::data& d, T r2) const { size_t pos = std::lower_bound(d.r2.begin(), d.r2.end(), r2) - d.r2.begin() - 1; size_t pos6 = 6 * pos; assert((pos6 + 5) < d.c.size()); @@ -183,13 +183,13 @@ template class Andrea : public TabulatorBase * @param d Table data * @param r2 value */ - T evalDer(const typename base::data &d, T r2) const { + T evalDer(const typename base::data& d, T r2) const { size_t pos = std::lower_bound(d.r2.begin(), d.r2.end(), r2) - d.r2.begin() - 1; size_t pos6 = 6 * pos; T dz = r2 - d.r2[pos]; return (d.c[pos6 + 1] + - dz * (2.0 * d.c[pos6 + 2] + - dz * (3.0 * d.c[pos6 + 3] + dz * (4.0 * d.c[pos6 + 4] + dz * (5.0 * d.c[pos6 + 5]))))); + dz * (2.0 * d.c[pos6 + 2] + + dz * (3.0 * d.c[pos6 + 3] + dz * (4.0 * d.c[pos6 + 4] + dz * (5.0 * d.c[pos6 + 5]))))); } /** @@ -268,12 +268,12 @@ template class Andrea : public TabulatorBase throw std::runtime_error("Andrea spline: try to increase utol/ftol"); // create final reversed c and r2 - assert( td.c.size() % 6 == 0 ); - assert( td.c.size() / (td.r2.size()-1) == 6 ); - assert( std::is_sorted(td.r2.rbegin(), td.r2.rend())); - std::reverse(td.r2.begin(), td.r2.end()); // reverse all elements - for (size_t i = 0; i < td.c.size() / 2; i += 6) // reverse knot order in packets of six - std::swap_ranges(td.c.begin()+i, td.c.begin()+i+6, td.c.end()-i-6); // c++17 only + assert(td.c.size() % 6 == 0); + assert(td.c.size() / (td.r2.size() - 1) == 6); + assert(std::is_sorted(td.r2.rbegin(), td.r2.rend())); + std::reverse(td.r2.begin(), td.r2.end()); // reverse all elements + for (size_t i = 0; i < td.c.size() / 2; i += 6) // reverse knot order in packets of six + std::swap_ranges(td.c.begin() + i, td.c.begin() + i + 6, td.c.end() - i - 6); // c++17 only return td; } }; @@ -290,6 +290,21 @@ TEST_CASE("[Faunus] Andrea") { spline.setTolerance(2e-6, 1e-4); // ftol carries no meaning auto d = spline.generate(f, 0, 10); + CHECK(d.r2.size() == 19); + CHECK(d.c.size() == 108); + CHECK(d.numKnots() == 19); + CHECK(d.rmin2 == Approx(0.0)); + CHECK(d.rmax2 == Approx(10.0)); + CHECK(d.r2.at(0) == Approx(0.0)); + CHECK(d.r2.at(1) == Approx(0.212991)); + CHECK(d.r2.at(2) == Approx(0.782554)); + CHECK(d.r2.back() == Approx(10.0)); + + CHECK(d.c.at(0) == Approx(2.0)); + CHECK(d.c.at(1) == Approx(0.0)); + CHECK(d.c.at(2) == Approx(0.5)); + CHECK(d.c.back() == Approx(-0.0441931)); + CHECK(spline.eval(d, 1e-9) == Approx(f(1e-9))); CHECK(spline.eval(d, 5) == Approx(f(5))); CHECK(spline.eval(d, 10) == Approx(f(10))); @@ -317,4 +332,3 @@ TEST_CASE("[Faunus] Andrea") { CHECK(spline.evalDer(d, x) == Approx(f_prime_exact(x))); } #endif -