ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
ObservableStat.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2013-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
20#include "config/config.hpp"
21
22#include "ObservableStat.hpp"
23
26
28
29#include <cstddef>
30#include <span>
31#include <string>
32#include <vector>
33
34namespace ScriptInterface {
35namespace Analysis {
36
37/**
38 * @brief Generate an observable summary.
39 * @param[in] system The system to analyze.
40 * @param[in] obs The observable handle.
41 * @param[in] calc_sp Whether to compute a scalar pressure.
42 */
44 Observable_stat const &obs, bool const calc_sp) {
45 auto const obs_dim = obs.get_chunk_size();
46
47 auto const get_obs_contribs = [obs_dim,
48 calc_sp](std::span<double> const &views) {
49 if (obs_dim == 1ul) {
50 return std::vector<Variant>(views.begin(), views.end());
51 }
52 assert(obs_dim == 9ul);
53 assert(views.size() % 9ul == 0ul);
54 std::vector<Variant> out;
55 for (std::size_t i = 0ul; i < views.size() / 9ul; ++i) {
56 auto const view = views.subspan(i * 9ul, 9ul);
57 if (calc_sp) {
58 auto const trace = view[0ul] + view[4ul] + view[8ul];
59 out.emplace_back(trace / 3.);
60 } else {
61 auto const flat_matrix = std::vector<double>(view.begin(), view.end());
62 out.emplace_back(flat_matrix);
63 }
64 }
65 return out;
66 };
67
68 auto const get_obs_contrib =
69 [&get_obs_contribs](std::span<double> const &views) -> Variant {
70 return get_obs_contribs(views)[0ul];
71 };
72
73 std::unordered_map<std::string, Variant> dict;
74 dict["kinetic"] = get_obs_contrib(obs.kinetic);
75 dict["external_fields"] = get_obs_contrib(obs.external_fields);
76
77 {
78 auto values = std::vector<double>(obs_dim);
79 for (std::size_t i = 0ul; i < obs_dim; ++i) {
80 values[i] = obs.accumulate(0., i);
81 }
82 dict["total"] = get_obs_contrib({values.data(), obs_dim});
83 }
84
85 auto const n_bonds = static_cast<int>(system.bonded_ias->get_next_key());
86 for (int bond_id = 0; bond_id < n_bonds; ++bond_id) {
87 if (system.bonded_ias->get_zero_based_type(bond_id) != 0) {
88 dict["bonded," + std::to_string(bond_id)] =
89 get_obs_contrib(obs.bonded_contribution(bond_id));
90 }
91 }
92
93 auto const n_nonbonded =
94 system.nonbonded_ias->get_max_seen_particle_type() + 1;
95 for (int i = 0; i < n_nonbonded; ++i) {
96 for (int j = i; j < n_nonbonded; ++j) {
97 auto const indices = std::to_string(i) + "," + std::to_string(j);
98 dict["non_bonded_intra," + indices] =
99 get_obs_contrib(obs.non_bonded_intra_contribution(i, j));
100 dict["non_bonded_inter," + indices] =
101 get_obs_contrib(obs.non_bonded_inter_contribution(i, j));
102 }
103 }
104
105#ifdef ELECTROSTATICS
106 {
107 auto const values = get_obs_contribs(obs.coulomb);
108 for (std::size_t i = 0ul; i < values.size(); ++i) {
109 dict["coulomb," + std::to_string(i)] = values[i];
110 }
111 }
112#endif // ELECTROSTATICS
113
114#ifdef DIPOLES
115 {
116 auto const values = get_obs_contribs(obs.dipolar);
117 for (std::size_t i = 0ul; i < values.size(); ++i) {
118 dict["dipolar," + std::to_string(i)] = values[i];
119 }
120 }
121#endif // DIPOLES
122
123#ifdef VIRTUAL_SITES
124 {
125 auto const values = get_obs_contribs(obs.virtual_sites);
126 for (std::size_t i = 0ul; i < values.size(); ++i) {
127 dict["virtual_sites," + std::to_string(i)] = values[i];
128 }
129 }
130#endif // VIRTUAL_SITES
131
132 return dict;
133}
134
136 VariantMap const &parameters) {
137 auto &system = get_system();
138 if (name == "calculate_energy") {
139 auto const obs = system.calculate_energy();
140 return get_summary(system, *obs, false);
141 }
142 if (name == "calculate_scalar_pressure") {
143 auto const obs = system.calculate_pressure();
144 return get_summary(system, *obs, true);
145 }
146 if (name == "calculate_pressure_tensor") {
147 auto const obs = system.calculate_pressure();
148 return get_summary(system, *obs, false);
149 }
150 return {};
151}
152
153} // namespace Analysis
154} // namespace ScriptInterface
Data structures for bonded interactions.
Observable for the pressure and energy.
Variant do_call_method(std::string const &name, VariantMap const &parameters) override
boost::string_ref name() const
Main system class.
This file contains the defaults for ESPResSo.
static auto get_summary(::System::System const &system, Observable_stat const &obs, bool const calc_sp)
Generate an observable summary.
T get_value(Variant const &v)
Extract value of specific type T from a Variant.
std::unordered_map< std::string, Variant > VariantMap
Definition Variant.hpp:69
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.
Definition Variant.hpp:67
Various procedures concerning interactions between particles.