ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
CosPersistenceAngles.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2019-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19#ifndef OBSERVABLES_PERSISTENCEANGLES_HPP
20#define OBSERVABLES_PERSISTENCEANGLES_HPP
21
22#include "BoxGeometry.hpp"
23#include "PidObservable.hpp"
24#include "system/System.hpp"
25
26#include <utils/Vector.hpp>
27
28#include <cassert>
29#include <cmath>
30#include <cstddef>
31#include <stdexcept>
32#include <utility>
33#include <vector>
34
35namespace Observables {
36
37/** Calculate bond angles in a polymer.
38 * The \a ith entry in the result vector corresponds to the
39 * averaged cosine of the angle between bonds that are \a i bonds apart.
40 */
42public:
44 explicit CosPersistenceAngles(std::vector<int> ids)
45 : PidObservable(std::move(ids)) {
46 if (this->ids().size() < 3)
47 throw std::runtime_error("At least 3 particles are required");
48 }
49
50 std::vector<double>
51 evaluate(boost::mpi::communicator const &comm,
52 ParticleReferenceRange const &local_particles,
53 const ParticleObservables::traits<Particle> &traits) const override {
54
55 auto const positions_sorted = detail::get_all_particle_positions(
56 comm, local_particles, ids(), traits, false);
57
58 if (comm.rank() != 0) {
59 return {};
60 }
61
62 auto const &box_geo = *System::get_system().box_geo;
63 auto const no_of_angles = n_values();
64 auto const no_of_bonds = no_of_angles + 1;
65 std::vector<double> angles(no_of_angles);
66 std::vector<Utils::Vector3d> bond_vectors(no_of_bonds);
67 auto get_bond_vector = [&](auto index) {
68 return box_geo.get_mi_vector(positions_sorted[index + 1],
69 positions_sorted[index]);
70 };
71 for (std::size_t i = 0; i < no_of_bonds; ++i) {
72 auto const tmp = get_bond_vector(i);
73 bond_vectors[i] = tmp / tmp.norm();
74 }
75 // calculate angles between neighbouring bonds, next neighbours, etc...
76 for (std::size_t i = 0; i < no_of_angles; ++i) {
77 auto average = 0.0;
78 for (std::size_t j = 0; j < no_of_angles - i; ++j) {
79 average += bond_vectors[j] * bond_vectors[j + i + 1];
80 }
81 angles[i] = average / static_cast<double>(no_of_angles - i);
82 }
83
84 return angles;
85 }
86 std::vector<std::size_t> shape() const override {
87 assert(ids().size() >= 2);
88 return {ids().size() - 2};
89 }
90};
91
92} // Namespace Observables
93
94#endif
Vector implementation and trait types for boost qvm interoperability.
Calculate bond angles in a polymer.
std::vector< double > evaluate(boost::mpi::communicator const &comm, ParticleReferenceRange const &local_particles, const ParticleObservables::traits< Particle > &traits) const override
std::vector< std::size_t > shape() const override
CosPersistenceAngles(std::vector< int > ids)
std::size_t n_values() const
Size of the flat array returned by the observable.
std::vector< int > const & ids() const
std::shared_ptr< BoxGeometry > box_geo
std::vector< std::reference_wrapper< Particle const > > ParticleReferenceRange
System & get_system()