36#include "constraints/Constraints.hpp"
37#include "constraints/ShapeBasedConstraint.hpp"
39#include "system/System.hpp"
45#include <boost/mpi/collectives/all_reduce.hpp>
46#include <boost/optional.hpp>
58 for (
int i = 0; i < 3; ++i)
59 v[i] = box_geo.
length()[i] * rng();
65 double const phi = acos(1. - 2. * rng());
66 double const theta = 2. *
Utils::pi() * rng();
67 v[0] = sin(phi) * cos(theta);
68 v[1] = sin(phi) * sin(theta);
88 std::vector<std::vector<Utils::Vector3d>>
const &positions,
90 BoxGeometry const &box_geo,
double const min_distance,
91 int const respect_constraints) {
94 auto operator()(
double const a,
double const b)
const {
95 return std::min(a, b);
100 if (respect_constraints) {
105 std::dynamic_pointer_cast<const Constraints::ShapeBasedConstraint>(c);
110 cs->calc_dist(folded_pos, d, v);
119 if (min_distance > 0.) {
121 auto local_mindist_sq = std::numeric_limits<double>::infinity();
124 local_mindist_sq = std::min(local_mindist_sq, d.norm2());
126 auto const global_mindist_sq =
127 boost::mpi::all_reduce(
::comm_cart, local_mindist_sq, reduce_min{});
128 if (std::sqrt(global_mindist_sq) < min_distance) {
132 auto const min_distance_sq = min_distance * min_distance;
133 for (
auto const &p : positions) {
134 for (
auto const &m : p) {
144std::vector<std::vector<Utils::Vector3d>>
146 int const beads_per_chain,
double const bond_length,
147 std::vector<Utils::Vector3d>
const &start_positions,
148 double const min_distance,
int const max_tries,
149 int const use_bond_angle,
double const bond_angle,
150 int const respect_constraints,
int const seed) {
152 auto const &box_geo = *system.
box_geo;
154 dist = std::uniform_real_distribution<double>(
155 0.0, 1.0)]()
mutable {
return dist(mt); };
157 std::vector<std::vector<Utils::Vector3d>> positions(n_polymers);
158 for (
auto &p : positions) {
159 p.reserve(beads_per_chain);
164 min_distance, respect_constraints);
167 for (std::size_t p = 0; p < start_positions.size(); p++) {
168 if (is_valid_pos(start_positions[p])) {
169 positions[p].push_back(start_positions[p]);
171 throw std::runtime_error(
"Invalid start positions.");
177 auto draw_monomer_position = [&](
int p,
int m) {
179 return (p < start_positions.size()) ? start_positions[p]
183 if (not use_bond_angle or m < 2) {
187 auto const last_vec = positions[p][m - 1] - positions[p][m - 2];
188 return positions[p][m - 1] +
190 bond_angle, -last_vec);
194 auto draw_valid_monomer_position =
195 [&](
int p,
int m) -> boost::optional<Utils::Vector3d> {
196 for (
auto i = 0; i < max_tries; i++) {
197 auto const trial_pos = draw_monomer_position(p, m);
198 if (is_valid_pos(trial_pos)) {
207 for (
int p = 0; p < n_polymers; ++p) {
208 for (
int attempts_poly = 0; attempts_poly < max_tries; attempts_poly++) {
210 while (positions[p].size() < beads_per_chain) {
211 auto pos = draw_valid_monomer_position(
212 p,
static_cast<int>(positions[p].size()));
216 positions[p].push_back(*
pos);
217 }
else if (not positions[p].empty()) {
219 positions[p].pop_back();
221 if (rejections > max_tries) {
233 if (positions[p].size() == beads_per_chain) {
240 if (positions[p].size() < beads_per_chain)
241 throw std::runtime_error(
"Failed to create polymer positions.");
Vector implementation and trait types for boost qvm interoperability.
__shared__ int pos[MAXDEPTH *THREADS5/WARPSIZE]
Utils::Vector3d const & length() const
Box length.
auto folded_position(Utils::Vector3d const &p) const
Calculate coordinates folded to primary simulation box.
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< BoxGeometry > box_geo
boost::mpi::communicator comm_cart
The communicator.
__device__ void vector_product(float const *a, float const *b, float *out)
Constraints< ParticleRange, Constraint > constraints
std::mt19937 mt19937(T &&seed)
Mersenne Twister with warmup.
DEVICE_QUALIFIER constexpr T pi()
Ratio of diameter and circumference of a circle.
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(Utils::Vector3d const &pos, std::vector< std::vector< Utils::Vector3d > > const &positions, CellStructure const &cell_structure, 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.
Describes a cell structure / cell system.
ParticleRange local_particles() const