ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
pressure_inline.hpp
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#pragma once
23
24#include "config/config.hpp"
25
29
30#include "BoxGeometry.hpp"
31#include "Observable_stat.hpp"
32#include "Particle.hpp"
33#include "errorhandling.hpp"
34#include "exclusions.hpp"
35#include "forces_inline.hpp"
36
37#include <utils/Vector.hpp>
39
40#include <boost/variant.hpp>
41
42#include <cstdio>
43#include <optional>
44#include <span>
45#include <string>
46#include <tuple>
47
48/** Calculate non-bonded energies between a pair of particles.
49 * @param p1 pointer to particle 1.
50 * @param p2 pointer to particle 2.
51 * @param d vector between p1 and p2.
52 * @param dist distance between p1 and p2.
53 * @param ia_params non-bonded interaction kernels.
54 * @param bonded_ias bonded interaction kernels.
55 * @param kernel_forces Coulomb force kernel.
56 * @param kernel_pressure Coulomb pressure kernel.
57 * @param[in,out] obs_pressure pressure observable.
58 */
60 Particle const &p1, Particle const &p2, Utils::Vector3d const &d,
61 double dist, IA_parameters const &ia_params,
62 [[maybe_unused]] BondedInteractionsMap const &bonded_ias,
65 Observable_stat &obs_pressure) {
66#ifdef EXCLUSIONS
67 if (do_nonbonded(p1, p2))
68#endif
69 {
70 auto const force = calc_central_radial_force(ia_params, d, dist).f +
71#ifdef THOLE
72 thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias,
73 kernel_forces) +
74#endif
75 calc_non_central_force(p1, p2, ia_params, d, dist).f;
76 auto const stress = Utils::tensor_product(d, force);
77 obs_pressure.add_non_bonded_contribution(p1.type(), p2.type(), p1.mol_id(),
78 p2.mol_id(), flatten(stress));
79 }
80
81#ifdef ELECTROSTATICS
82 if (!obs_pressure.coulomb.empty() and kernel_pressure != nullptr) {
83 /* real space Coulomb */
84 auto const p_coulomb = (*kernel_pressure)(p1.q() * p2.q(), d, dist);
85
86 for (std::size_t i = 0u; i < 3u; i++) {
87 for (std::size_t j = 0u; j < 3u; j++) {
88 obs_pressure.coulomb[i * 3u + j] += p_coulomb(i, j);
89 }
90 }
91 }
92#endif // ELECTROSTATICS
93
94#ifdef DIPOLES
95 /* real space magnetic dipole-dipole */
96 if (Dipoles::get_dipoles().impl->solver) {
97 fprintf(stderr, "calculating pressure for magnetostatics which doesn't "
98 "have it implemented\n");
99 }
100#endif // DIPOLES
101}
102
103inline std::optional<Utils::Matrix<double, 3, 3>>
105 Bonded_IA_Parameters const &iaparams, Particle const &p1,
106 Particle const &p2, BoxGeometry const &box_geo,
108 auto const dx = box_geo.get_mi_vector(p1.pos(), p2.pos());
109 auto const result = calc_bond_pair_force(p1, p2, iaparams, dx, kernel);
110 if (result) {
111 auto const &force = result.value();
112
113 return Utils::tensor_product(force, dx);
114 }
115
116 return {};
117}
118
119inline std::optional<Utils::Matrix<double, 3, 3>>
121 Particle const &p1, Particle const &p2,
122 Particle const &p3,
123 BoxGeometry const &box_geo) {
124 if ((boost::get<AngleHarmonicBond>(&iaparams) != nullptr) ||
125 (boost::get<AngleCosineBond>(&iaparams) != nullptr) ||
126#ifdef TABULATED
127 (boost::get<TabulatedAngleBond>(&iaparams) != nullptr) ||
128#endif
129 (boost::get<AngleCossquareBond>(&iaparams) != nullptr)) {
130 auto const dx21 = -box_geo.get_mi_vector(p1.pos(), p2.pos());
131 auto const dx31 = box_geo.get_mi_vector(p3.pos(), p1.pos());
132
133 auto const result =
134 calc_bonded_three_body_force(iaparams, box_geo, p1, p2, p3);
135 if (result) {
136 Utils::Vector3d force2, force3;
137 std::tie(std::ignore, force2, force3) = result.value();
138
139 return Utils::tensor_product(force2, dx21) +
140 Utils::tensor_product(force3, dx31);
141 }
142 } else {
143 runtimeWarningMsg() << "Unsupported bond type " +
144 std::to_string(iaparams.which()) +
145 " in pressure calculation.";
147 }
148
149 return {};
150}
151
152inline std::optional<Utils::Matrix<double, 3, 3>> calc_bonded_pressure_tensor(
153 Bonded_IA_Parameters const &iaparams, Particle const &p1,
154 std::span<Particle *> partners, BoxGeometry const &box_geo,
156 switch (number_of_partners(iaparams)) {
157 case 1:
158 return calc_bonded_virial_pressure_tensor(iaparams, p1, *partners[0],
159 box_geo, kernel);
160 case 2:
161 return calc_bonded_three_body_pressure_tensor(iaparams, p1, *partners[0],
162 *partners[1], box_geo);
163 default:
164 runtimeWarningMsg() << "Unsupported bond type " +
165 std::to_string(iaparams.which()) +
166 " in pressure calculation.";
168 }
169}
170
171/** Calculate kinetic pressure (aka energy) for one particle.
172 * @param[in] p1 particle for which to calculate pressure
173 * @param[out] obs_pressure pressure observable
174 */
175inline void add_kinetic_virials(Particle const &p1,
176 Observable_stat &obs_pressure) {
177 if (p1.is_virtual())
178 return;
179
180 /* kinetic pressure */
181 for (std::size_t k = 0u; k < 3u; k++)
182 for (std::size_t l = 0u; l < 3u; l++)
183 obs_pressure.kinetic[k * 3u + l] += p1.v()[k] * p1.v()[l] * p1.mass();
184}
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
int number_of_partners(Bonded_IA_Parameters const &iaparams)
Get the number of bonded partners for the specified bond.
boost::variant< NoneBond, FeneBond, HarmonicBond, QuarticBond, BondedCoulomb, BondedCoulombSR, AngleHarmonicBond, AngleCosineBond, AngleCossquareBond, DihedralBond, TabulatedDistanceBond, TabulatedAngleBond, TabulatedDihedralBond, ThermalizedBond, RigidBond, IBMTriel, IBMVolCons, IBMTribend, OifGlobalForcesBond, OifLocalForcesBond, VirtualBond > Bonded_IA_Parameters
Variant in which to store the parameters of an individual bonded interaction.
container for bonded interactions.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
Observable for the pressure and energy.
std::span< double > coulomb
Contribution(s) from Coulomb interactions.
void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2, std::span< const double > data)
std::span< double > kinetic
Contribution from linear and angular kinetic energy (accumulated).
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeWarningMsg()
bool do_nonbonded(Particle const &p1, Particle const &p2)
Determine if the non-bonded interactions between p1 and p2 should be calculated.
Force calculation.
ParticleForce calc_central_radial_force(IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist)
std::optional< Utils::Vector3d > calc_bond_pair_force(Particle const &p1, Particle const &p2, Bonded_IA_Parameters const &iaparams, Utils::Vector3d const &dx, Coulomb::ShortRangeForceKernel::kernel_type const *kernel)
Compute the bonded interaction force between particle pairs.
ParticleForce calc_non_central_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist)
std::optional< std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > > calc_bonded_three_body_force(Bonded_IA_Parameters const &iaparams, BoxGeometry const &box_geo, Particle const &p1, Particle const &p2, Particle const &p3)
Solver const & get_dipoles()
Definition dipoles.cpp:49
Matrix< T, N, M > tensor_product(const Vector< T, N > &x, const Vector< T, M > &y)
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.
std::optional< Utils::Matrix< double, 3, 3 > > calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams, Particle const &p1, Particle const &p2, Particle const &p3, BoxGeometry const &box_geo)
std::optional< Utils::Matrix< double, 3, 3 > > calc_bonded_virial_pressure_tensor(Bonded_IA_Parameters const &iaparams, Particle const &p1, Particle const &p2, BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel)
Solver::ShortRangeForceKernel kernel_type
Solver::ShortRangePressureKernel kernel_type
Parameters for non-bonded interactions.
Utils::Vector3d f
force.
Definition Particle.hpp:316
Struct holding all information for one particle.
Definition Particle.hpp:395
auto is_virtual() const
Definition Particle.hpp:518
auto const & mass() const
Definition Particle.hpp:452
auto const & q() const
Definition Particle.hpp:508
auto const & v() const
Definition Particle.hpp:433
auto const & type() const
Definition Particle.hpp:418
auto const & mol_id() const
Definition Particle.hpp:416
auto const & pos() const
Definition Particle.hpp:431
Matrix representation with static size.
Definition matrix.hpp:65
Utils::Vector3d thole_pair_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeForceKernel::kernel_type const *kernel)
Calculate Thole force.
Definition thole.hpp:44