ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
polymer.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21/** \file
22 * This file contains everything needed to create a start-up configuration
23 * of (possibly charged) polymer chains with counterions and salt molecules,
24 * assigning velocities to the particles and cross-linking the polymers if
25 * necessary.
26 *
27 * The corresponding header file is polymer.hpp.
28 */
29
30#include "polymer.hpp"
31
32#include "BoxGeometry.hpp"
34#include "communication.hpp"
35#include "constraints/Constraints.hpp"
36#include "constraints/ShapeBasedConstraint.hpp"
37#include "random.hpp"
38#include "system/System.hpp"
39
40#include <utils/Vector.hpp>
42
43#include <boost/mpi/collectives/all_reduce.hpp>
44
45#include <cmath>
46#include <cstddef>
47#include <functional>
48#include <memory>
49#include <numbers>
50#include <optional>
51#include <stdexcept>
52#include <vector>
53
54template <class RNG>
55static Utils::Vector3d random_position(BoxGeometry const &box_geo, RNG &rng) {
57 for (auto i = 0u; i < 3u; ++i)
58 v[i] = box_geo.length()[i] * rng();
59 return v;
60}
61
62template <class RNG> static Utils::Vector3d random_unit_vector(RNG &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);
68 v[2] = cos(phi);
69 v /= v.norm();
70 return v;
71}
72
73/** Determines whether a given position @p pos is valid, i.e., it doesn't
74 * collide with existing or buffered particles, nor with existing constraints
75 * (if @c respect_constraints).
76 * @param system the system containing the particles
77 * @param pos the trial position in question
78 * @param positions buffered positions to respect
79 * @param box_geo Box geometry
80 * @param min_distance threshold for the minimum distance between
81 * trial position and buffered/existing particles
82 * @param respect_constraints whether to respect constraints
83 * @return true if valid position, false if not.
84 */
85static bool
87 std::vector<std::vector<Utils::Vector3d>> const &positions,
88 BoxGeometry const &box_geo, double const min_distance,
89 int const respect_constraints) {
90
91 struct reduce_min {
92 auto operator()(double const a, double const b) const {
93 return std::min(a, b);
94 }
95 };
96
97 // check if constraint is violated
98 if (respect_constraints) {
99 Utils::Vector3d const folded_pos = box_geo.folded_position(pos);
100
101 for (auto &c : *system.constraints) {
102 auto cs =
103 std::dynamic_pointer_cast<const Constraints::ShapeBasedConstraint>(c);
104 if (cs) {
105 double d;
107
108 cs->calc_dist(folded_pos, d, v);
109
110 if (d <= 0) {
111 return false;
112 }
113 }
114 }
115 }
116
117 if (min_distance > 0.) {
118 // check for collision with existing particles
119 auto local_mindist_sq = std::numeric_limits<double>::infinity();
120 for (auto const &p : system.cell_structure->local_particles()) {
121 auto const d = box_geo.get_mi_vector(pos, p.pos());
122 local_mindist_sq = std::min(local_mindist_sq, d.norm2());
123 }
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) {
127 return false;
128 }
129
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) {
134 return false;
135 }
136 }
137 }
138 }
139 return true;
140}
141
142std::vector<std::vector<Utils::Vector3d>>
143draw_polymer_positions(System::System const &system, int const n_polymers,
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) {
149
150 auto const &box_geo = *system.box_geo;
151 auto rng = [mt = Random::mt19937(static_cast<unsigned>(seed)),
152 dist = std::uniform_real_distribution<double>(
153 0.0, 1.0)]() mutable { return dist(mt); };
154
155 std::vector<std::vector<Utils::Vector3d>> positions(n_polymers);
156 for (auto &p : positions) {
157 p.reserve(beads_per_chain);
158 }
159
160 auto is_valid_pos = [&](Utils::Vector3d const &pos) {
161 return is_valid_position(system, pos, positions, box_geo, min_distance,
162 respect_constraints);
163 };
164
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]);
168 } else {
169 throw std::runtime_error("Invalid start positions.");
170 }
171 }
172
173 /* Draw a monomer position, obeying angle and starting position
174 * constraints where appropriate. */
175 auto draw_monomer_position = [&](int p, int m) {
176 if (m == 0) {
177 return (p < start_positions.size()) ? start_positions[p]
178 : random_position(box_geo, rng);
179 }
180
181 if (not use_bond_angle or m < 2) {
182 return positions[p][m - 1] + bond_length * random_unit_vector(rng);
183 }
184
185 auto const last_vec = positions[p][m - 1] - positions[p][m - 2];
186 return positions[p][m - 1] +
188 bond_angle, -last_vec);
189 };
190
191 /* Try up to max_tries times to draw a valid position */
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)) {
197 return trial_pos;
198 }
199 }
200
201 return std::nullopt;
202 };
203
204 // create remaining monomers' positions by backtracking.
205 for (int p = 0; p < n_polymers; ++p) {
206 for (int attempts_poly = 0; attempts_poly < max_tries; attempts_poly++) {
207 int rejections = 0;
208 while (positions[p].size() < beads_per_chain) {
209 auto pos = draw_valid_monomer_position(
210 p, static_cast<int>(positions[p].size()));
211
212 if (pos) {
213 /* Move on one position */
214 positions[p].push_back(*pos);
215 } else if (not positions[p].empty()) {
216 /* Go back one position and try again */
217 positions[p].pop_back();
218 rejections++;
219 if (rejections > max_tries) {
220 /* Give up for this try. */
221 break;
222 }
223 } else {
224 /* Give up for this try. */
225 break;
226 }
227 }
228
229 /* If the polymer has not full length, we try again. Otherwise we can
230 * move on to the next polymer. */
231 if (positions[p].size() == beads_per_chain) {
232 break;
233 }
234 }
235
236 /* We did not get a complete polymer, but have exceeded the maximal
237 * number of tries, which means failure. */
238 if (positions[p].size() < beads_per_chain)
239 throw std::runtime_error("Failed to create polymer positions.");
240 }
241 return positions;
242}
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.
Main system class.
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< Constraints::Constraints > constraints
std::shared_ptr< BoxGeometry > box_geo
T norm() const
Definition Vector.hpp:138
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.
Definition random.hpp:194
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)
Definition polymer.cpp:55
static Utils::Vector3d random_unit_vector(RNG &rng)
Definition polymer.cpp:62
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.
Definition polymer.cpp:143
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 ...
Definition polymer.cpp:86
This file contains everything needed to create a start-up configuration of polymer chains which may r...
Random number generation using Philox.