42#include <boost/format.hpp>
43#include <boost/mpi/collectives/all_reduce.hpp>
44#include <boost/mpi/collectives/broadcast.hpp>
45#include <boost/mpi/communicator.hpp>
66 throw std::domain_error(
"Invalid particle id: " + std::to_string(p_id));
69 if (std::isnan(
pos[0]) or std::isnan(
pos[1]) or std::isnan(
pos[2]) or
70 std::isinf(
pos[0]) or std::isinf(
pos[1]) or std::isinf(
pos[2])) {
71 throw std::domain_error(
"Particle position must be finite");
77 auto bitfield =
static_cast<uint8_t
>(0
u);
79 bitfield |=
static_cast<uint8_t
>(1u);
81 bitfield |=
static_cast<uint8_t
>(2u);
83 bitfield |=
static_cast<uint8_t
>(4u);
87static auto error_msg(std::string
const &name, std::string
const &reason) {
88 std::stringstream msg;
89 msg <<
"attribute '" << name <<
"' of 'ParticleHandle' " << reason;
98 auto const q = get_value<Utils::Vector4d>(value);
99 if (q.norm2() == 0.) {
100 throw std::domain_error(
error_msg(name,
"must be non-zero"));
105#ifdef THERMOSTAT_PER_PARTICLE
107#ifdef PARTICLE_ANISOTROPY
111 return get_value<Utils::Vector3d>(value);
114 return get_value<double>(value);
121 throw std::domain_error(
"Invalid particle id: " + std::to_string(p_id));
124 auto &cell_structure = *system.cell_structure;
125 auto ptr = cell_structure.get_local_particle(p_id);
126 if (ptr !=
nullptr and ptr->is_ghost()) {
129 auto const n_found = boost::mpi::all_reduce(
130 comm,
static_cast<int>(ptr !=
nullptr), std::plus<>());
132 throw std::runtime_error(
"Particle with id " + std::to_string(p_id) +
138template <
typename T,
class F>
139T ParticleHandle::get_particle_property(F
const &fun)
const {
142 std::optional<T> ret;
143 if (ptr ==
nullptr) {
152T ParticleHandle::get_particle_property(T
const &(
Particle::*getter)()
154 return get_particle_property<T>(
155 [getter](
Particle const &p) {
return (p.*getter)(); });
159void ParticleHandle::set_particle_property(F
const &fun)
const {
160 auto const &comm = context()->get_comm();
162 if (ptr !=
nullptr) {
169void ParticleHandle::set_particle_property(T &(
Particle::*setter)(),
171 set_particle_property(
172 [&value, setter](
Particle &p) { (p.*setter)() = get_value<T>(value); });
175ParticleHandle::ParticleHandle() {
177 {
"id", AutoParameter::read_only, [
this]() {
return m_pid; }},
181 auto const new_type = get_value<int>(value);
183 throw std::domain_error(
184 error_msg(
"type",
"must be an integer >= 0"));
193 auto const pos = get_value<Utils::Vector3d>(value);
199 auto const pos = p.pos();
200 auto const image_box = p.image_box();
216 if (get_value<double>(value) <= 0.) {
217 throw std::domain_error(
error_msg(
"mass",
"must be a float > 0"));
224 if (std::abs(get_value<double>(value) - default_mass) > 1e-10) {
225 throw std::runtime_error(
"Feature MASS not compiled in");
237 if (get_value<double>(value) != 0.) {
238 throw std::runtime_error(
"Feature ELECTROSTATICS not compiled in");
246 set_particle_property([&value](
Particle &p) {
247 auto const director = get_value<Utils::Vector3d>(value).normalized();
258 set_particle_property([&quat](
Particle &p) { p.
quat() = quat; });
268 set_particle_property([&value](
Particle &p) {
269 auto const rotation_flag =
277 ::detail::get_nth_bit(rotation_bits, 1),
278 ::detail::get_nth_bit(rotation_bits, 2)}};
282 set_particle_property([&value](
Particle &p) {
283 auto const omega = get_value<Utils::Vector3d>(value);
293 set_particle_property([&value](
Particle &p) {
294 auto const torque = get_value<Utils::Vector3d>(value);
306 set_particle_property([&value](
Particle &p) {
307 auto const dip = get_value<Utils::Vector3d>(value);
318#ifdef DIPOLE_FIELD_TRACKING
325#ifdef ROTATIONAL_INERTIA
332#ifdef LB_ELECTROHYDRODYNAMICS
339#ifdef EXTERNAL_FORCES
342 set_particle_property([&value](
Particle &p) {
343 auto const fix_flag =
351 ::detail::get_nth_bit(fixed, 1),
352 ::detail::get_nth_bit(fixed, 2)}};
369 auto const propagation = get_value<int>(value);
372 "propagation",
"propagation combination not accepted: " +
378#ifdef THERMOSTAT_PER_PARTICLE
394 {
"pos_folded", AutoParameter::read_only,
399 {
"lees_edwards_offset",
404 {
"lees_edwards_flag", AutoParameter::read_only,
406 {
"image_box", AutoParameter::read_only,
408 {
"node", AutoParameter::read_only,
414 auto const mol_id = get_value<int>(value);
416 throw std::domain_error(
417 error_msg(
"mol_id",
"must be an integer >= 0"));
422#ifdef VIRTUAL_SITES_RELATIVE
426 set_particle_property(
436 auto const array = get_value<std::vector<Variant>>(value);
437 if (array.size() != 3) {
440 vs_relative.
distance = get_value<double>(array[1]);
441 vs_relative.to_particle_id = get_value<int>(array[0]);
442 vs_relative.rel_orientation =
446 "vs_relative",
"must take the form [id, distance, quaternion]"));
448 set_particle_property([&vs_relative](
Particle &p) {
450 if (vs_relative.to_particle_id != -1) {
451 p.propagation() = PropagationMode::TRANS_VS_RELATIVE |
452 PropagationMode::ROT_VS_RELATIVE;
459 return std::vector<Variant>{{vs_rel.to_particle_id, vs_rel.distance,
466 set_particle_property([&value](
Particle &p) {
469 auto const dict = get_value<VariantMap>(value);
470 if (dict.count(
"f_swim") != 0) {
471 swim.f_swim = get_value<double>(dict.at(
"f_swim"));
473 if (dict.count(
"is_engine_force_on_fluid") != 0) {
474 auto const is_engine_force_on_fluid =
475 get_value<bool>(dict.at(
"is_engine_force_on_fluid"));
476 swim.is_engine_force_on_fluid = is_engine_force_on_fluid;
484 {
"f_swim", swim.f_swim},
485 {
"is_engine_force_on_fluid", swim.is_engine_force_on_fluid},
500 auto &cell_structure = *system.cell_structure;
501 if (
auto p1 = cell_structure.get_local_particle(pid1)) {
504 if (
auto p2 = cell_structure.get_local_particle(pid2)) {
516 auto &cell_structure = *system.cell_structure;
517 if (
auto p1 = cell_structure.get_local_particle(pid1)) {
520 if (
auto p2 = cell_structure.get_local_particle(pid2)) {
525void ParticleHandle::particle_exclusion_sanity_checks(
int pid1,
528 throw std::runtime_error(
"Particles cannot exclude themselves (id " +
529 std::to_string(pid1) +
")");
531 std::ignore = get_real_particle(context()->get_comm(), pid1);
532 std::ignore = get_real_particle(context()->get_comm(), pid2);
536Variant ParticleHandle::do_call_method(std::string
const &name,
538 if (name ==
"set_param_parallel") {
539 auto const param_name = get_value<std::string>(
params,
"name");
540 if (
params.count(
"value") == 0) {
541 throw Exception(
"Parameter '" + param_name +
"' is missing.");
543 auto const &value =
params.at(
"value");
544 context()->parallel_try_catch(
545 [&]() { do_set_parameter(param_name, value); });
548 if (name ==
"get_bonds_view") {
549 if (not context()->is_head_node()) {
553 std::vector<std::vector<int>> bonds_flat;
554 for (
auto const &&bond_view : bond_list) {
555 std::vector<int> bond_flat;
556 bond_flat.emplace_back(bond_view.bond_id());
557 for (
auto const pid : bond_view.partner_ids()) {
558 bond_flat.emplace_back(pid);
560 bonds_flat.emplace_back(std::move(bond_flat));
564 if (name ==
"add_bond") {
566 auto const bond_id = get_value<int>(
params,
"bond_id");
567 auto const part_id = get_value<std::vector<int>>(
params,
"part_id");
568 auto const bond_view =
569 BondView(bond_id, {part_id.data(), part_id.size()});
570 p.
bonds().insert(bond_view);
572 }
else if (name ==
"del_bond") {
574 auto const bond_id = get_value<int>(
params,
"bond_id");
575 auto const part_id = get_value<std::vector<int>>(
params,
"part_id");
576 auto const bond_view =
577 BondView(bond_id, {part_id.data(), part_id.size()});
578 auto &bond_list = p.
bonds();
579 auto it = std::find(bond_list.begin(), bond_list.end(), bond_view);
580 if (it != bond_list.end()) {
584 }
else if (name ==
"delete_all_bonds") {
585 set_particle_property([&](
Particle &p) { p.
bonds().clear(); });
586 }
else if (name ==
"is_valid_bond_id") {
587 auto const bond_id = get_value<int>(
params,
"bond_id");
588 return ::bonded_ia_params.get_zero_based_type(bond_id) != 0;
590 if (name ==
"remove_particle") {
591 context()->parallel_try_catch([&]() {
595 }
else if (name ==
"is_virtual") {
596 if (not context()->is_head_node()) {
600#ifdef VIRTUAL_SITES_RELATIVE
601 }
else if (name ==
"vs_relate_to") {
602 if (not context()->is_head_node()) {
605 auto const other_pid = get_value<int>(
params,
"pid");
606 auto const override_cutoff_check =
607 get_value<bool>(
params,
"override_cutoff_check");
608 if (m_pid == other_pid) {
609 throw std::invalid_argument(
"A virtual site cannot relate to itself");
612 throw std::domain_error(
"Invalid particle id: " +
613 std::to_string(other_pid));
628 p_current, p_relate_to, *system.box_geo, system.get_min_global_cut(),
629 override_cutoff_check);
630 set_parameter(
"vs_relative",
Variant{std::vector<Variant>{
634 }
else if (name ==
"has_exclusion") {
635 auto const other_pid = get_value<int>(
params,
"pid");
638 return p->has_exclusion(other_pid);
641 if (name ==
"add_exclusion") {
642 auto const other_pid = get_value<int>(
params,
"pid");
643 context()->parallel_try_catch(
644 [&]() { particle_exclusion_sanity_checks(m_pid, other_pid); });
647 }
else if (name ==
"del_exclusion") {
648 auto const other_pid = get_value<int>(
params,
"pid");
649 context()->parallel_try_catch(
650 [&]() { particle_exclusion_sanity_checks(m_pid, other_pid); });
653 }
else if (name ==
"set_exclusions") {
654 std::vector<int> exclusion_list;
656 auto const pid = get_value<int>(
params,
"p_ids");
657 exclusion_list.push_back(pid);
659 exclusion_list = get_value<std::vector<int>>(
params,
"p_ids");
661 context()->parallel_try_catch([&]() {
662 for (
auto const pid : exclusion_list) {
663 particle_exclusion_sanity_checks(m_pid, pid);
666 set_particle_property([
this, &exclusion_list](
Particle &p) {
670 for (
auto const pid : exclusion_list) {
676 }
else if (name ==
"get_exclusions") {
677 if (not context()->is_head_node()) {
681 return Variant{std::vector<int>{excl_list.begin(), excl_list.end()}};
685 if (name ==
"rotate_particle") {
687 auto const axis = get_value<Utils::Vector3d>(
params,
"axis");
688 auto const angle = get_value<double>(
params,
"angle");
692 if (name ==
"convert_vector_body_to_space") {
693 return get_particle_property<std::vector<double>>(
695 auto const vec = get_value<Utils::Vector3d>(
params,
"vec");
699 if (name ==
"convert_vector_space_to_body") {
700 return get_particle_property<std::vector<double>>(
702 auto const vec = get_value<Utils::Vector3d>(
params,
"vec");
712 std::array<std::string, 3>>{{
714 "Setting 'dip' is sufficient as it defines the scalar dipole moment."}},
715 {{
"quat",
"director",
716 "Setting 'quat' is sufficient as it defines the director."}},
718 "Setting 'dip' would overwrite 'quat'. Set 'quat' and 'dipm' instead."}},
720 "Setting 'dip' would overwrite 'director'. Set 'director' and "
726 auto const n_extra_args =
params.size() -
params.count(
"id");
727 auto const has_param = [&
params](std::string
const key) {
728 return params.count(key) == 1;
730 m_pid = (has_param(
"id")) ? get_value<int>(
params,
"id")
734 if (!has_param(
"id")) {
735 auto head_node_reference = m_pid;
736 boost::mpi::broadcast(context()->get_comm(), head_node_reference, 0);
737 assert(m_pid == head_node_reference &&
"global max_seen_pid has diverged");
742 if (n_extra_args == 0) {
746 auto const pos = get_value<Utils::Vector3d>(
params,
"pos");
747 context()->parallel_try_catch([&]() {
750 auto const &cell_structure = *system.cell_structure;
751 auto ptr = cell_structure.get_local_particle(m_pid);
752 if (ptr !=
nullptr) {
753 throw std::invalid_argument(
"Particle " + std::to_string(m_pid) +
759 context()->parallel_try_catch([&]() {
762 if (not has_param(
"__cpt_sentinel")) {
764 boost::format(
"Contradicting particle attributes: '%s' and '%s'. %s");
766 if (has_param(prop1) and has_param(prop2)) {
767 auto const err_msg = boost::str(formatter % prop1 % prop2 % reason);
768 throw std::invalid_argument(err_msg);
779 context()->parallel_try_catch([&]() {
782 std::set<std::string>
const skip = {
783 "pos_folded",
"pos",
"quat",
"director",
"id",
784 "exclusions",
"dip",
"node",
"image_box",
"bonds",
785 "lees_edwards_flag",
"__cpt_sentinel",
793 for (std::string name : {
"quat",
"director",
"dip"}) {
794 if (has_param(name)) {
795 do_set_parameter(name,
params.at(name));
800 for (
auto const &name : get_parameter_insertion_order()) {
801 if (
params.count(name) != 0ul and skip.count(name) == 0ul) {
802 do_set_parameter(name,
params.at(name));
805 if (not has_param(
"type")) {
806 do_set_parameter(
"type", 0);
Vector implementation and trait types for boost qvm interoperability.
__global__ float float * torque
__shared__ int pos[MAXDEPTH *THREADS5/WARPSIZE]
Data structures for bonded interactions.
Immutable view on a bond.
virtual boost::mpi::communicator const & get_comm() const =0
Context * context() const
Responsible context.
void on_particle_change()
Called every time a particle property changes.
std::shared_ptr< BoxGeometry > box_geo
std::shared_ptr< InteractionsNonBonded > nonbonded_ias
static Vector< T, N > broadcast(T const &s)
Create a vector that has all entries set to one value.
std::vector< T > as_vector() const
This file contains the defaults for ESPResSo.
void delete_exclusion(Particle &p, int p_id)
Remove exclusion from particle if possible.
void add_exclusion(Particle &p, int p_id)
Insert an exclusion if not already set.
static uint8_t bitfield_from_flag(Utils::Vector3i const &flag)
static void particle_checks(int p_id, Utils::Vector3d const &pos)
static auto get_quaternion_safe(std::string const &name, Variant const &value)
static auto const contradicting_arguments_quat
static void local_add_exclusion(int pid1, int pid2)
Locally add an exclusion to a particle.
static auto get_gamma_safe(Variant const &value)
static auto quat2vector(Utils::Quaternion< double > const &q)
static void local_remove_exclusion(int pid1, int pid2)
Locally remove an exclusion to a particle.
static auto get_real_particle(boost::mpi::communicator const &comm, int p_id)
static auto error_msg(std::string const &name, std::string const &reason)
std::unordered_map< std::string, Variant > VariantMap
auto make_vector_of_variants(std::vector< T > const &v)
boost::make_recursive_variant< None, bool, int, std::size_t, double, std::string, ObjectRef, Utils::Vector3b, Utils::Vector3i, Utils::Vector2d, Utils::Vector3d, Utils::Vector4d, std::vector< int >, std::vector< double >, std::vector< boost::recursive_variant_ >, std::unordered_map< int, boost::recursive_variant_ >, std::unordered_map< std::string, boost::recursive_variant_ > >::type Variant
Possible types for parameters.
T reduce_optional(boost::mpi::communicator const &comm, std::optional< T > const &result)
Reduce an optional on the head node.
Quaternion< T > convert_director_to_quaternion(Vector< T, 3 > const &d)
Convert director to quaternion.
Vector< T, 3 > convert_quaternion_to_director(Quaternion< T > const &quat)
Convert quaternion to director.
Various procedures concerning interactions between particles.
void make_new_particle(int p_id, Utils::Vector3d const &pos)
Create a new particle and attach it to a cell.
const Particle & get_particle_data(int p_id)
Get particle data.
int get_particle_node(int p_id)
Get the MPI rank which owns the a specific particle.
void set_particle_pos(int p_id, Utils::Vector3d const &pos)
Move particle to a new position.
void remove_particle(int p_id)
Remove particle with a given identity.
int get_maximal_particle_id()
Get maximal particle id.
void on_particle_type_change(int p_id, int old_type, int new_type)
Particles creation and deletion.
std::string propagation_bitmask_to_string(int propagation)
Convert a propagation modes bitmask to a string.
bool is_valid_propagation_combination(int propagation)
Note for developers: when enabling new propagation mode combinations, make sure every single line of ...
This file contains all subroutines required to process rotational motion.
Utils::Vector3d convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec)
std::pair< Utils::Quaternion< double >, double > convert_dip_to_quat(const Utils::Vector3d &dip)
convert a dipole moment to quaternions and dipolar strength
Utils::Vector3d convert_vector_space_to_body(const Particle &p, const Utils::Vector3d &v)
void local_rotate_particle(Particle &p, const Utils::Vector3d &axis_space_frame, const double phi)
Rotate the particle p around the NORMALIZED axis aSpaceFrame by amount phi.
static SteepestDescentParameters params
Currently active steepest descent instance.
Properties of a self-propelled particle.
bool swimming
Is the particle a swimmer.
The following properties define, with respect to which real particle a virtual site is placed and at ...
Struct holding all information for one particle.
auto const & dip_fld() const
bool has_exclusion(int pid) const
auto const & swimming() const
auto const & lees_edwards_offset() const
auto const & propagation() const
auto const & rinertia() const
auto const & mass() const
Utils::compact_vector< int > & exclusions()
auto const & quat() const
auto const & rotation() const
auto const & vs_relative() const
auto const & gamma() const
auto const & gamma_rot() const
auto const & lees_edwards_flag() const
auto const & torque() const
auto const & fixed() const
auto const & ext_force() const
auto const & omega() const
auto const & image_box() const
auto const & type() const
auto const & ext_torque() const
auto const & bonds() const
auto const & dipm() const
auto const & mol_id() const
auto const & mu_E() const
auto const & force() const
Quaternion representation.
std::tuple< Utils::Quaternion< double >, double > calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, BoxGeometry const &box_geo, double min_global_cut, bool override_cutoff_check)
Calculate the rotation quaternion and distance between two particles.