Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Reuse distto from statistics in polymer/diamond code, fix minor perfo… #3402

Merged
merged 2 commits into from
Jan 8, 2020
Merged

Reuse distto from statistics in polymer/diamond code, fix minor perfo… #3402

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c928db6
Reuse distto from statistics in polymer/diamond code, fix minor perfo…
fweik Jan 7, 2020
1b11d43
Merge branch 'python' into polymer_
kodiakhq[bot] Jan 8, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 8 additions & 13 deletions src/core/diamond.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,19 +27,14 @@
* The corresponding header file is diamond.hpp.
*/

#include <cmath>
#include <cstddef>

#include "PartCfg.hpp"
#include "bonded_interactions/bonded_interaction_data.hpp"
#include "communication.hpp"
#include "constraints.hpp"
#include "constraints/ShapeBasedConstraint.hpp"
#include "diamond.hpp"
#include "global.hpp"
#include "integrate.hpp"
#include "polymer.hpp"
#include "particle_data.hpp"
#include "random.hpp"
#include "statistics.hpp"

#include <utils/Vector.hpp>

#include <cmath>
#include <stdexcept>

int create_counterions(PartCfg &partCfg, int const N_CI, int part_id,
int const mode, double const shield, int const max_try,
Expand All @@ -53,7 +48,7 @@ int create_counterions(PartCfg &partCfg, int const N_CI, int part_id,
pos[0] = box_geo.length()[0] * d_random();
pos[1] = box_geo.length()[1] * d_random();
pos[2] = box_geo.length()[2] * d_random();
if ((mode != 0) or (mindist(partCfg, pos) > shield))
if ((mode != 0) or (distto(partCfg, pos) > shield))
break;
}
if (cnt1 >= max_try)
Expand Down
2 changes: 0 additions & 2 deletions src/core/diamond.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,6 @@
*/

#include "PartCfg.hpp"
#include "Particle.hpp"
#include "utils/Vector.hpp"

/** C implementation of 'counterions \<N_CI\> [options]'.
* @param N_CI number of counterions to create
Expand Down
84 changes: 36 additions & 48 deletions src/core/polymer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,64 +27,54 @@
* The corresponding header file is polymer.hpp.
*/

#include <cmath>
#include <cstddef>
#include "polymer.hpp"

#include "PartCfg.hpp"
#include "bonded_interactions/bonded_interaction_data.hpp"
#include "communication.hpp"
#include "constraints.hpp"
#include "constraints/ShapeBasedConstraint.hpp"
#include "global.hpp"
#include "integrate.hpp"
#include "polymer.hpp"
#include "random.hpp"

#include <utils/Vector.hpp>
#include <utils/constants.hpp>
#include <utils/math/sqr.hpp>
#include <utils/math/vec_rotate.hpp>

Utils::Vector3d random_position(std::function<double()> const &generate_rn) {
#include <cmath>
#include <random>
#include <stdexcept>

template <class RNG> static Utils::Vector3d random_position(RNG &rng) {
Utils::Vector3d v;
for (int i = 0; i < 3; ++i)
v[i] = box_geo.length()[i] * generate_rn();
v[i] = box_geo.length()[i] * rng();
return v;
}

Utils::Vector3d random_unit_vector(std::function<double()> const &generate_rn) {
template <class RNG> static Utils::Vector3d random_unit_vector(RNG &rng) {
Utils::Vector3d v;
double const phi = acos(1. - 2. * generate_rn());
double const theta = 2. * Utils::pi() * generate_rn();
double const phi = acos(1. - 2. * rng());
double const theta = 2. * Utils::pi() * rng();
v[0] = sin(phi) * cos(theta);
v[1] = sin(phi) * sin(theta);
v[2] = cos(phi);
v /= v.norm();
return v;
}

double mindist(PartCfg &partCfg, Utils::Vector3d const &pos) {
if (partCfg.empty()) {
return std::min(std::min(box_geo.length()[0], box_geo.length()[1]),
box_geo.length()[2]);
}

auto const mindist = std::accumulate(
partCfg.begin(), partCfg.end(), std::numeric_limits<double>::infinity(),
[&pos](double mindist, Particle const &p) {
return std::min(mindist, get_mi_vector(pos, p.r.p, box_geo).norm2());
});

if (mindist < std::numeric_limits<double>::infinity())
return std::sqrt(mindist);
return -1.0;
}

bool is_valid_position(
Utils::Vector3d const *pos,
std::vector<std::vector<Utils::Vector3d>> const *positions,
PartCfg &partCfg, double const min_distance,
int const respect_constraints) {
/** Determines whether a given position @p pos is valid, i.e., it doesn't
* collide with existing or buffered particles, nor with existing constraints
* (if @c respect_constraints).
* @param pos the trial position in question
* @param positions buffered positions to respect
* @param partCfg existing particles to respect
* @param min_distance threshold for the minimum distance between
* trial position and buffered/existing particles
* @param respect_constraints whether to respect constraints
* @return true if valid position, false if not.
*/
static bool
is_valid_position(Utils::Vector3d const *pos,
std::vector<std::vector<Utils::Vector3d>> const *positions,
PartCfg &partCfg, double const min_distance,
int const respect_constraints) {
Utils::Vector3d const folded_pos = folded_position(*pos, box_geo);
// check if constraint is violated
if (respect_constraints) {
Expand All @@ -106,14 +96,14 @@ bool is_valid_position(

if (min_distance > 0) {
// check for collision with existing particles
if (mindist(partCfg, *pos) < min_distance) {
if (distto(partCfg, *pos, -1) < min_distance) {
return false;
}
// check for collision with buffered positions
double buff_mindist = std::numeric_limits<double>::infinity();
double h;
for (auto const p : *positions) {
for (auto m : p) {
for (auto const &p : *positions) {
for (auto const &m : p) {
h = (folded_position(*pos, box_geo) - folded_position(m, box_geo))
.norm2();
buff_mindist = std::min(h, buff_mindist);
Expand All @@ -136,8 +126,9 @@ draw_polymer_positions(PartCfg &partCfg, int const n_polymers,
std::vector<std::vector<Utils::Vector3d>> positions(
n_polymers, std::vector<Utils::Vector3d>(beads_per_chain));

std::mt19937 mt(seed);
std::uniform_real_distribution<double> dist(0.0, 1.0);
auto rng = [mt = std::mt19937{static_cast<unsigned>(seed)},
dist = std::uniform_real_distribution<double>(
0.0, 1.0)]() mutable { return dist(mt); };

Utils::Vector3d trial_pos;
int attempts_mono, attempts_poly;
Expand All @@ -158,7 +149,7 @@ draw_polymer_positions(PartCfg &partCfg, int const n_polymers,
// first monomer for all polymers
if (start_positions.empty()) {
do {
trial_pos = random_position([&]() { return dist(mt); });
trial_pos = random_position(rng);
attempts_mono++;
} while ((not is_valid_position(&trial_pos, &positions, partCfg,
min_distance, respect_constraints)) and
Expand All @@ -180,21 +171,18 @@ draw_polymer_positions(PartCfg &partCfg, int const n_polymers,
do {
if (m == 0) {
// m == 0 is only true after a failed attempt to position a polymer
trial_pos = random_position([&]() { return dist(mt); });
trial_pos = random_position(rng);
} else if (not use_bond_angle or m < 2) {
// random step, also necessary if angle is set, placing the second
// bead
trial_pos =
positions[p][m - 1] +
bond_length * random_unit_vector([&]() { return dist(mt); });
positions[p][m - 1] + bond_length * random_unit_vector(rng);
} else {
// use prescribed angle
Utils::Vector3d last_vec = positions[p][m - 1] - positions[p][m - 2];
trial_pos = positions[p][m - 1] +
Utils::vec_rotate(
vector_product(last_vec, random_unit_vector([&]() {
return dist(mt);
})),
vector_product(last_vec, random_unit_vector(rng)),
bond_angle, -last_vec);
}
attempts_mono++;
Expand Down
25 changes: 2 additions & 23 deletions src/core/polymer.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,31 +31,10 @@

#include "PartCfg.hpp"
#include "Particle.hpp"
#include "utils/Vector.hpp"

Utils::Vector3d random_position(std::function<double()> const &generate_rn);
Utils::Vector3d random_unit_vector(std::function<double()> const &generate_rn);
#include <utils/Vector.hpp>

/** Returns the minimum distance between position @p pos and all existing
* particles.
*/
double mindist(PartCfg &partCfg, Utils::Vector3d const &pos);

/** Determines whether a given position @p pos is valid, i.e., it doesn't
* collide with existing or buffered particles, nor with existing constraints
* (if @c respect_constraints).
* @param pos the trial position in question
* @param positions buffered positions to respect
* @param partCfg existing particles to respect
* @param min_distance threshold for the minimum distance between
* trial position and buffered/existing particles
* @param respect_constraints whether to respect constraints
* @return true if valid position, false if not.
*/
bool is_valid_position(
Utils::Vector3d const *pos,
std::vector<std::vector<Utils::Vector3d>> const *positions,
PartCfg const &partCfg, double min_distance, int respect_constraints);
#include <vector>

/** Determines valid polymer positions and returns them.
* @param n_polymers how many polymers to create
Expand Down
2 changes: 1 addition & 1 deletion src/core/statistics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ IntList nbhood(PartCfg &partCfg, const Utils::Vector3d &pos, double r_catch,
* position of a particle).
* @return the minimal distance of a particle to coordinates @p pos
*/
double distto(PartCfg &partCfg, const Utils::Vector3d &pos, int pid);
double distto(PartCfg &partCfg, const Utils::Vector3d &pos, int pid = -1);

/** Append particles' positions in %p partCfg to #configs
* @param partCfg @copybrief PartCfg
Expand Down