ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
pressure.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#include "config/config.hpp"
23
24#include "BoxGeometry.hpp"
25#include "Observable_stat.hpp"
26#include "Particle.hpp"
27#include "ParticleRange.hpp"
32#include "pressure_inline.hpp"
33#include "short_range_loop.hpp"
34#include "system/System.hpp"
36
37#include <utils/Vector.hpp>
38#include <utils/flatten.hpp>
39
40#include <algorithm>
41#include <cstddef>
42#include <memory>
43#include <span>
44
45namespace System {
46std::shared_ptr<Observable_stat> System::calculate_pressure() {
47
48 auto obs_pressure_ptr = std::make_shared<Observable_stat>(
49 9ul, static_cast<std::size_t>(bonded_ias->get_next_key()),
50 nonbonded_ias->get_max_seen_particle_type());
51
52 if (long_range_interactions_sanity_checks()) {
53 return obs_pressure_ptr;
54 }
55
56 auto &obs_pressure = *obs_pressure_ptr;
57
58 on_observable_calc();
59
60 auto const volume = box_geo->volume();
61 auto const local_parts = cell_structure->local_particles();
62
63 for (auto const &p : local_parts) {
64 add_kinetic_virials(p, obs_pressure);
65 }
66
67 auto const coulomb_force_kernel = coulomb.pair_force_kernel();
68 auto const coulomb_pressure_kernel = coulomb.pair_pressure_kernel();
69
71 [this, coulomb_force_kernel_ptr = get_ptr(coulomb_force_kernel),
72 &obs_pressure](Particle const &p1, int bond_id,
73 std::span<Particle *> partners) {
74 auto const &iaparams = *bonded_ias->at(bond_id);
75 auto const result = calc_bonded_pressure_tensor(
76 iaparams, p1, partners, *box_geo, coulomb_force_kernel_ptr);
77 if (result) {
78 auto const &tensor = result.value();
79 /* pressure tensor part */
80 for (std::size_t k = 0u; k < 3u; k++)
81 for (std::size_t l = 0u; l < 3u; l++)
82 obs_pressure.bonded_contribution(bond_id)[k * 3u + l] +=
83 tensor(k, l);
84
85 return false;
86 }
87 return true;
88 },
89 [coulomb_force_kernel_ptr = get_ptr(coulomb_force_kernel),
90 coulomb_pressure_kernel_ptr = get_ptr(coulomb_pressure_kernel), this,
91 &obs_pressure](Particle const &p1, Particle const &p2,
92 Distance const &d) {
93 auto const &ia_params =
94 nonbonded_ias->get_ia_param(p1.type(), p2.type());
95 add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), ia_params,
96 *bonded_ias, coulomb_force_kernel_ptr,
97 coulomb_pressure_kernel_ptr, obs_pressure);
98 },
99 *cell_structure, maximal_cutoff(), bonded_ias->maximal_cutoff());
100
101#ifdef ELECTROSTATICS
102 /* calculate k-space part of electrostatic interaction. */
103 auto const coulomb_pressure = coulomb.calc_pressure_long_range(local_parts);
104 std::ranges::copy(coulomb_pressure, obs_pressure.coulomb.begin() + 9u);
105#endif
106#ifdef DIPOLES
107 /* calculate k-space part of magnetostatic interaction. */
109#endif
110
111#ifdef VIRTUAL_SITES_RELATIVE
112 if (!obs_pressure.virtual_sites.empty()) {
113 auto const vs_pressure = vs_relative_pressure_tensor(*cell_structure);
114 std::ranges::copy(Utils::flatten(vs_pressure),
115 obs_pressure.virtual_sites.begin());
116 }
117#endif
118
119 obs_pressure.rescale(volume);
120
121 obs_pressure.mpi_reduce();
122 return obs_pressure_ptr;
123}
124} // namespace System
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
std::shared_ptr< Observable_stat > calculate_pressure()
Calculate the pressure from a virial expansion.
Definition pressure.cpp:46
This file contains the defaults for ESPResSo.
const T * get_ptr(std::optional< T > const &opt)
Solver const & get_dipoles()
Definition dipoles.cpp:49
void flatten(Range const &v, OutputIterator out)
Flatten a range of ranges.
Definition flatten.hpp:64
Various procedures concerning interactions between particles.
void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, Utils::Vector3d const &d, double dist, IA_parameters const &ia_params, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeForceKernel::kernel_type const *kernel_forces, Coulomb::ShortRangePressureKernel::kernel_type const *kernel_pressure, Observable_stat &obs_pressure)
Calculate non-bonded energies between a pair of particles.
std::optional< Utils::Matrix< double, 3, 3 > > calc_bonded_pressure_tensor(Bonded_IA_Parameters const &iaparams, Particle const &p1, std::span< Particle * > partners, BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel)
void add_kinetic_virials(Particle const &p1, Observable_stat &obs_pressure)
Calculate kinetic pressure (aka energy) for one particle.
Utils::Matrix< double, 3, 3 > vs_relative_pressure_tensor(CellStructure const &cell_structure)
Definition relative.cpp:194
void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, CellStructure &cell_structure, double pair_cutoff, double bond_cutoff, VerletCriterion const &verlet_criterion={})
void calc_pressure_long_range() const
Definition dipoles.cpp:193
Distance vector and length handed to pair kernels.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & type() const
Definition Particle.hpp:418