ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
BondAngles.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_BONDANGLES_HPP
20#define OBSERVABLES_BONDANGLES_HPP
21
22#include "BoxGeometry.hpp"
23#include "PidObservable.hpp"
24#include "system/System.hpp"
25
26#include <utils/Vector.hpp>
27
28#include <algorithm>
29#include <cassert>
30#include <cmath>
31#include <cstddef>
32#include <stdexcept>
33#include <utility>
34#include <vector>
35
36namespace Observables {
37
38/** Calculate bond angles between particles in a polymer.
39 * For @f$ n @f$ bonded particles, return the @f$ n-2 @f$ angles along the
40 * chain, in radians.
41 */
42class BondAngles : public PidObservable {
43public:
45 explicit BondAngles(std::vector<int> ids) : 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 auto const positions_sorted = detail::get_all_particle_positions(
55 comm, local_particles, ids(), traits, false);
56
57 if (comm.rank() != 0) {
58 return {};
59 }
60
61 auto const &box_geo = *System::get_system().box_geo;
62 std::vector<double> res(n_values());
63 auto v1 = box_geo.get_mi_vector(positions_sorted[1], positions_sorted[0]);
64 auto n1 = v1.norm();
65 for (std::size_t i = 0, end = n_values(); i < end; i++) {
66 auto v2 = box_geo.get_mi_vector(positions_sorted[i + 2],
67 positions_sorted[i + 1]);
68 auto const n2 = v2.norm();
69 auto const cosine =
70 std::clamp((v1 * v2) / (n1 * n2), -TINY_COS_VALUE, TINY_COS_VALUE);
71 /* to reduce computational time, after calculating an angle ijk, the
72 * vector r_ij takes the value r_jk, but to orient it correctly, it has
73 * to be multiplied -1; it's cheaper to do this operation on a double
74 * than on a vector of doubles
75 */
76 res[i] = acos(-cosine);
77 v1 = v2;
78 n1 = n2;
79 }
80 return res;
81 }
82 std::vector<std::size_t> shape() const override {
83 assert(ids().size() >= 2);
84 return {ids().size() - 2};
85 }
86};
87
88} // Namespace Observables
89
90#endif
Vector implementation and trait types for boost qvm interoperability.
Calculate bond angles between particles in a polymer.
std::vector< std::size_t > shape() const override
BondAngles(std::vector< int > ids)
std::vector< double > evaluate(boost::mpi::communicator const &comm, ParticleReferenceRange const &local_particles, const ParticleObservables::traits< Particle > &traits) const override
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
#define TINY_COS_VALUE
Tiny angle cutoff for cosine calculations.
Definition config.hpp:75
std::vector< std::reference_wrapper< Particle const > > ParticleReferenceRange
System & get_system()