ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
BondDihedrals.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_PARTICLEDIHEDRALS_HPP
20#define OBSERVABLES_PARTICLEDIHEDRALS_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 dihedral angles between particles in a polymer.
38 * For @f$ n @f$ bonded particles, return the @f$ n-3 @f$ dihedral angles
39 * along the chain, in radians.
40 *
41 * The sign of the dihedral angles follow the IUPAC nomenclature for the
42 * Newman projection, specifically section "Torsion Angle" pages 2220-2221
43 * in @cite moss96a.
44 *
45 */
47public:
49 explicit BondDihedrals(std::vector<int> ids) : PidObservable(std::move(ids)) {
50 if (this->ids().size() < 4)
51 throw std::runtime_error("At least 4 particles are required");
52 }
53
54 std::vector<double>
55 evaluate(boost::mpi::communicator const &comm,
56 ParticleReferenceRange const &local_particles,
57 const ParticleObservables::traits<Particle> &traits) const override {
58 auto const positions_sorted = detail::get_all_particle_positions(
59 comm, local_particles, ids(), traits, false);
60
61 if (comm.rank() != 0) {
62 return {};
63 }
64
65 auto const &box_geo = *System::get_system().box_geo;
66 std::vector<double> res(n_values());
67 auto v1 = box_geo.get_mi_vector(positions_sorted[1], positions_sorted[0]);
68 auto v2 = box_geo.get_mi_vector(positions_sorted[2], positions_sorted[1]);
69 auto c1 = Utils::vector_product(v1, v2);
70 for (std::size_t i = 0, end = n_values(); i < end; i++) {
71 auto v3 = box_geo.get_mi_vector(positions_sorted[i + 3],
72 positions_sorted[i + 2]);
73 auto c2 = vector_product(v2, v3);
74 /* the 2-argument arctangent returns an angle in the range [-pi, pi] that
75 * allows for an unambiguous determination of the 4th particle position */
76 res[i] = atan2((Utils::vector_product(c1, c2) * v2) / v2.norm(), c1 * c2);
77 v1 = v2;
78 v2 = v3;
79 c1 = c2;
80 }
81 return res;
82 }
83 std::vector<std::size_t> shape() const override {
84 assert(ids().size() >= 3);
85 return {ids().size() - 3};
86 }
87};
88
89} // Namespace Observables
90
91#endif
Vector implementation and trait types for boost qvm interoperability.
Calculate dihedral angles between particles in a polymer.
BondDihedrals(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::vector< std::size_t > shape() 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
__device__ void vector_product(float const *a, float const *b, float *out)
std::vector< std::reference_wrapper< Particle const > > ParticleReferenceRange
System & get_system()
Vector< T, 3 > vector_product(Vector< T, 3 > const &a, Vector< T, 3 > const &b)
Definition Vector.hpp:368