35#include "constraints/Constraints.hpp"
36#include "constraints/ShapeBasedConstraint.hpp"
38#include "system/System.hpp"
43#include <boost/mpi/collectives/all_reduce.hpp>
57 for (
auto i = 0u; i < 3u; ++i)
58 v[i] = box_geo.
length()[i] * rng();
64 double const phi = acos(1. - 2. * rng());
65 double const theta = 2. * std::numbers::pi * rng();
66 v[0] = sin(phi) * cos(theta);
67 v[1] = sin(phi) * sin(theta);
87 std::vector<std::vector<Utils::Vector3d>>
const &positions,
88 BoxGeometry const &box_geo,
double const min_distance,
89 int const respect_constraints) {
92 auto operator()(
double const a,
double const b)
const {
93 return std::min(a, b);
98 if (respect_constraints) {
103 std::dynamic_pointer_cast<const Constraints::ShapeBasedConstraint>(c);
108 cs->calc_dist(folded_pos, d, v);
117 if (min_distance > 0.) {
119 auto local_mindist_sq = std::numeric_limits<double>::infinity();
122 local_mindist_sq = std::min(local_mindist_sq, d.norm2());
124 auto const global_mindist_sq =
125 boost::mpi::all_reduce(
::comm_cart, local_mindist_sq, reduce_min{});
126 if (std::sqrt(global_mindist_sq) < min_distance) {
130 auto const min_distance_sq = min_distance * min_distance;
131 for (
auto const &p : positions) {
132 for (
auto const &m : p) {
133 if (box_geo.
get_mi_vector(pos, m).norm2() < min_distance_sq) {
142std::vector<std::vector<Utils::Vector3d>>
144 int const beads_per_chain,
double const bond_length,
145 std::vector<Utils::Vector3d>
const &start_positions,
146 double const min_distance,
int const max_tries,
147 int const use_bond_angle,
double const bond_angle,
148 int const respect_constraints,
int const seed) {
150 auto const &box_geo = *system.
box_geo;
152 dist = std::uniform_real_distribution<double>(
153 0.0, 1.0)]()
mutable {
return dist(mt); };
155 std::vector<std::vector<Utils::Vector3d>> positions(n_polymers);
156 for (
auto &p : positions) {
157 p.reserve(beads_per_chain);
162 respect_constraints);
165 for (std::size_t p = 0; p < start_positions.size(); p++) {
166 if (is_valid_pos(start_positions[p])) {
167 positions[p].push_back(start_positions[p]);
169 throw std::runtime_error(
"Invalid start positions.");
175 auto draw_monomer_position = [&](
int p,
int m) {
177 return (p < start_positions.size()) ? start_positions[p]
181 if (not use_bond_angle or m < 2) {
185 auto const last_vec = positions[p][m - 1] - positions[p][m - 2];
186 return positions[p][m - 1] +
188 bond_angle, -last_vec);
192 auto draw_valid_monomer_position =
193 [&](
int p,
int m) -> std::optional<Utils::Vector3d> {
194 for (
auto i = 0; i < max_tries; i++) {
195 auto const trial_pos = draw_monomer_position(p, m);
196 if (is_valid_pos(trial_pos)) {
205 for (
int p = 0; p < n_polymers; ++p) {
206 for (
int attempts_poly = 0; attempts_poly < max_tries; attempts_poly++) {
208 while (positions[p].size() < beads_per_chain) {
209 auto pos = draw_valid_monomer_position(
210 p,
static_cast<int>(positions[p].size()));
214 positions[p].push_back(*pos);
215 }
else if (not positions[p].empty()) {
217 positions[p].pop_back();
219 if (rejections > max_tries) {
231 if (positions[p].size() == beads_per_chain) {
238 if (positions[p].size() < beads_per_chain)
239 throw std::runtime_error(
"Failed to create polymer positions.");
Vector implementation and trait types for boost qvm interoperability.
auto folded_position(Utils::Vector3d const &pos) const
Calculate coordinates folded to primary simulation box.
Utils::Vector3d const & length() const
Box length.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< Constraints::Constraints > constraints
std::shared_ptr< BoxGeometry > box_geo
boost::mpi::communicator comm_cart
The communicator.
__device__ void vector_product(float const *a, float const *b, float *out)
std::mt19937 mt19937(T &&seed)
Mersenne Twister with warmup.
Vector3d vec_rotate(const Vector3d &axis, double angle, const Vector3d &vector)
Rotate a vector around an axis.
static Utils::Vector3d random_position(BoxGeometry const &box_geo, RNG &rng)
static Utils::Vector3d random_unit_vector(RNG &rng)
std::vector< std::vector< Utils::Vector3d > > draw_polymer_positions(System::System const &system, int const n_polymers, int const beads_per_chain, double const bond_length, std::vector< Utils::Vector3d > const &start_positions, double const min_distance, int const max_tries, int const use_bond_angle, double const bond_angle, int const respect_constraints, int const seed)
Determines valid polymer positions and returns them.
static bool is_valid_position(System::System const &system, Utils::Vector3d const &pos, std::vector< std::vector< Utils::Vector3d > > const &positions, BoxGeometry const &box_geo, double const min_distance, int const respect_constraints)
Determines whether a given position pos is valid, i.e., it doesn't collide with existing or buffered ...
This file contains everything needed to create a start-up configuration of polymer chains which may r...
Random number generation using Philox.