ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
energy_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/** \file
25 * Energy calculation.
26 */
27
28#include "config/config.hpp"
29
50
51#include "BoxGeometry.hpp"
52#include "Observable_stat.hpp"
53#include "Particle.hpp"
54#include "bond_error.hpp"
55#include "errorhandling.hpp"
56#include "exclusions.hpp"
57
58#include <utils/Vector.hpp>
59
60#include <boost/variant.hpp>
61
62#include <optional>
63#include <span>
64#include <string>
65
66/** Calculate non-bonded energies between a pair of particles.
67 * @param p1 particle 1.
68 * @param p2 particle 2.
69 * @param ia_params the interaction parameters between the two particles
70 * @param d vector between p1 and p2.
71 * @param dist distance between p1 and p2.
72 * @param bonded_ias bonded interaction kernels.
73 * @param coulomb_kernel Coulomb energy kernel.
74 * @return the short-range interaction energy between the two particles
75 */
77 Particle const &p1, Particle const &p2, IA_parameters const &ia_params,
78 Utils::Vector3d const &d, double const dist,
79 [[maybe_unused]] BondedInteractionsMap const &bonded_ias,
81
82 double ret = 0;
83
84#ifdef LENNARD_JONES
85 /* Lennard-Jones */
86 ret += lj_pair_energy(ia_params, dist);
87#endif
88
89#ifdef WCA
90 /* WCA */
91 ret += wca_pair_energy(ia_params, dist);
92#endif
93
94#ifdef LENNARD_JONES_GENERIC
95 /* Generic Lennard-Jones */
96 ret += ljgen_pair_energy(ia_params, dist);
97#endif
98
99#ifdef SMOOTH_STEP
100 /* smooth step */
101 ret += SmSt_pair_energy(ia_params, dist);
102#endif
103
104#ifdef HERTZIAN
105 /* Hertzian potential */
106 ret += hertzian_pair_energy(ia_params, dist);
107#endif
108
109#ifdef GAUSSIAN
110 /* Gaussian potential */
111 ret += gaussian_pair_energy(ia_params, dist);
112#endif
113
114#ifdef BMHTF_NACL
115 /* BMHTF NaCl */
116 ret += BMHTF_pair_energy(ia_params, dist);
117#endif
118
119#ifdef MORSE
120 /* Morse */
121 ret += morse_pair_energy(ia_params, dist);
122#endif
123
124#ifdef BUCKINGHAM
125 /* Buckingham */
126 ret += buck_pair_energy(ia_params, dist);
127#endif
128
129#ifdef SOFT_SPHERE
130 /* soft-sphere */
131 ret += soft_pair_energy(ia_params, dist);
132#endif
133
134#ifdef HAT
135 /* hat */
136 ret += hat_pair_energy(ia_params, dist);
137#endif
138
139#ifdef LJCOS2
140 /* Lennard-Jones */
141 ret += ljcos2_pair_energy(ia_params, dist);
142#endif
143
144#ifdef THOLE
145 /* Thole damping */
146 ret +=
147 thole_pair_energy(p1, p2, ia_params, d, dist, bonded_ias, coulomb_kernel);
148#endif
149
150#ifdef TABULATED
151 /* tabulated */
152 ret += tabulated_pair_energy(ia_params, dist);
153#endif
154
155#ifdef LJCOS
156 /* Lennard-Jones cosine */
157 ret += ljcos_pair_energy(ia_params, dist);
158#endif
159
160#ifdef GAY_BERNE
161 /* Gay-Berne */
162 ret += gb_pair_energy(p1.quat(), p2.quat(), ia_params, d, dist);
163#endif
164
165 return ret;
166}
167
168/** Add non-bonded and short-range Coulomb energies between a pair of particles
169 * to the energy observable.
170 * @param[in] p1 particle 1.
171 * @param[in] p2 particle 2.
172 * @param[in] d vector between p1 and p2.
173 * @param[in] dist distance between p1 and p2.
174 * @param[in] dist2 distance squared between p1 and p2.
175 * @param[in] ia_params non-bonded interaction kernels.
176 * @param[in] bonded_ias bonded interaction kernels.
177 * @param[in] coulomb_kernel Coulomb energy kernel.
178 * @param[in] dipoles_kernel Dipolar energy kernel.
179 * @param[in,out] obs_energy energy observable.
180 */
182 Particle const &p1, Particle const &p2, Utils::Vector3d const &d,
183 double const dist, double const dist2, IA_parameters const &ia_params,
184 [[maybe_unused]] BondedInteractionsMap const &bonded_ias,
187 Observable_stat &obs_energy) {
188
189#ifdef EXCLUSIONS
190 if (do_nonbonded(p1, p2))
191#endif
193 p1.type(), p2.type(), p1.mol_id(), p2.mol_id(),
194 calc_non_bonded_pair_energy(p1, p2, ia_params, d, dist, bonded_ias,
195 coulomb_kernel));
196
197#ifdef ELECTROSTATICS
198 if (!obs_energy.coulomb.empty() and coulomb_kernel != nullptr) {
199 auto const q1q2 = p1.q() * p2.q();
200 obs_energy.coulomb[0] += (*coulomb_kernel)(p1, p2, q1q2, d, dist);
201 }
202#endif
203
204#ifdef DIPOLES
205 if (!obs_energy.dipolar.empty() and dipoles_kernel != nullptr)
206 obs_energy.dipolar[0] += (*dipoles_kernel)(p1, p2, d, dist, dist2);
207#endif
208}
209
210inline std::optional<double>
212 std::span<Particle *> partners, BoxGeometry const &box_geo,
214 auto const n_partners = static_cast<int>(partners.size());
215
216 auto p2 = (n_partners > 0) ? partners[0] : nullptr;
217 auto p3 = (n_partners > 1) ? partners[1] : nullptr;
218 auto p4 = (n_partners > 2) ? partners[2] : nullptr;
219
220 if (n_partners == 1) {
221 auto const dx = box_geo.get_mi_vector(p1.pos(), p2->pos());
222 if (auto const *iap = boost::get<FeneBond>(&iaparams)) {
223 return iap->energy(dx);
224 }
225 if (auto const *iap = boost::get<HarmonicBond>(&iaparams)) {
226 return iap->energy(dx);
227 }
228 if (auto const *iap = boost::get<QuarticBond>(&iaparams)) {
229 return iap->energy(dx);
230 }
231#ifdef ELECTROSTATICS
232 if (auto const *iap = boost::get<BondedCoulomb>(&iaparams)) {
233 return iap->energy(p1.q() * p2->q(), dx);
234 }
235 if (auto const *iap = boost::get<BondedCoulombSR>(&iaparams)) {
236 return iap->energy(p1, *p2, dx, *kernel);
237 }
238#endif
239#ifdef BOND_CONSTRAINT
240 if (boost::get<RigidBond>(&iaparams)) {
241 return {0.};
242 }
243#endif
244#ifdef TABULATED
245 if (auto const *iap = boost::get<TabulatedDistanceBond>(&iaparams)) {
246 return iap->energy(dx);
247 }
248#endif
249 if (boost::get<VirtualBond>(&iaparams)) {
250 return {0.};
251 }
252 throw BondUnknownTypeError();
253 } // 1 partner
254 if (n_partners == 2) {
255 auto const vec1 = box_geo.get_mi_vector(p2->pos(), p1.pos());
256 auto const vec2 = box_geo.get_mi_vector(p3->pos(), p1.pos());
257 if (auto const *iap = boost::get<AngleHarmonicBond>(&iaparams)) {
258 return iap->energy(vec1, vec2);
259 }
260 if (auto const *iap = boost::get<AngleCosineBond>(&iaparams)) {
261 return iap->energy(vec1, vec2);
262 }
263 if (auto const *iap = boost::get<AngleCossquareBond>(&iaparams)) {
264 return iap->energy(vec1, vec2);
265 }
266 if (auto const *iap = boost::get<TabulatedAngleBond>(&iaparams)) {
267 return iap->energy(vec1, vec2);
268 }
269 if (boost::get<IBMTriel>(&iaparams)) {
270 runtimeWarningMsg() << "Unsupported bond type " +
271 std::to_string(iaparams.which()) +
272 " in energy calculation.";
273 return 0.;
274 }
275 throw BondUnknownTypeError();
276 } // 2 partners
277 if (n_partners == 3) {
278 // note: particles in a dihedral bond are ordered as p2-p1-p3-p4
279 auto const v12 = box_geo.get_mi_vector(p1.pos(), p2->pos());
280 auto const v23 = box_geo.get_mi_vector(p3->pos(), p1.pos());
281 auto const v34 = box_geo.get_mi_vector(p4->pos(), p3->pos());
282 if (auto const *iap = boost::get<DihedralBond>(&iaparams)) {
283 return iap->energy(v12, v23, v34);
284 }
285 if (auto const *iap = boost::get<TabulatedDihedralBond>(&iaparams)) {
286 return iap->energy(v12, v23, v34);
287 }
288 if (boost::get<IBMTribend>(&iaparams)) {
289 runtimeWarningMsg() << "Unsupported bond type " +
290 std::to_string(iaparams.which()) +
291 " in energy calculation.";
292 return 0.;
293 }
294 throw BondUnknownTypeError();
295 } // 3 partners
296 if (n_partners == 0) {
297 return 0.;
298 }
299
300 throw BondInvalidSizeError(n_partners);
301}
302
303/** Calculate kinetic energies from translation for one particle.
304 * @param p particle for which to calculate energies
305 */
306inline double translational_kinetic_energy(Particle const &p) {
307 return p.is_virtual() ? 0. : 0.5 * p.mass() * p.v().norm2();
308}
309
310/** Calculate kinetic energies from rotation for one particle.
311 * @param p particle for which to calculate energies
312 */
313inline double rotational_kinetic_energy(Particle const &p) {
314#ifdef ROTATION
315 return (p.can_rotate() and not p.is_virtual())
316 ? 0.5 * (hadamard_product(p.omega(), p.omega()) * p.rinertia())
317 : 0.0;
318#else
319 return 0.0;
320#endif
321}
322
323/** Calculate kinetic energies for one particle.
324 * @param p particle for which to calculate energies
325 */
Vector implementation and trait types for boost qvm interoperability.
Routines to calculate the Born-Meyer-Huggins-Tosi-Fumi potential between particle pairs.
double BMHTF_pair_energy(IA_parameters const &ia_params, double dist)
Calculate BMHTF potential energy.
Data structures for bonded interactions.
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.
Routines to calculate the Buckingham potential between particle pairs.
double buck_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Buckingham energy.
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 > dipolar
Contribution(s) from dipolar interactions.
This file contains the defaults for ESPResSo.
Routines to calculate the Gaussian potential between particle pairs.
double gaussian_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Gaussian energy.
std::optional< double > calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, std::span< Particle * > partners, BoxGeometry const &box_geo, Coulomb::ShortRangeEnergyKernel::kernel_type const *kernel)
double translational_kinetic_energy(Particle const &p)
Calculate kinetic energies from translation for one particle.
double calc_kinetic_energy(Particle const &p)
Calculate kinetic energies for one particle.
double calc_non_bonded_pair_energy(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_kernel)
Calculate non-bonded energies between a pair of particles.
double rotational_kinetic_energy(Particle const &p)
Calculate kinetic energies from rotation for one particle.
void add_non_bonded_pair_energy(Particle const &p1, Particle const &p2, Utils::Vector3d const &d, double const dist, double const dist2, IA_parameters const &ia_params, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_kernel, Dipoles::ShortRangeEnergyKernel::kernel_type const *dipoles_kernel, Observable_stat &obs_energy)
Add non-bonded and short-range Coulomb energies between a pair of particles to the energy observable.
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.
Routines to calculate the Gay-Berne potential between particle pairs.
double gb_pair_energy(Utils::Quaternion< double > const &qi, Utils::Quaternion< double > const &qj, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist)
Calculate Gay-Berne energy.
Routines to calculate the hat potential between particle pairs.
double hat_pair_energy(IA_parameters const &ia_params, double dist)
Calculate hat energy.
Definition hat.hpp:46
Routines to calculate the Hertzian potential between particle pairs.
double hertzian_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Hertzian energy.
Definition hertzian.hpp:49
Routines to calculate the Lennard-Jones potential between particle pairs.
double lj_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Lennard-Jones energy.
Definition lj.hpp:52
Routines to calculate the Lennard-Jones with cosine tail potential between particle pairs.
double ljcos2_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Lennard-Jones cosine squared energy.
Definition ljcos2.hpp:66
Routines to calculate the Lennard-Jones+cosine potential between particle pairs.
double ljcos_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Lennard-Jones cosine energy.
Definition ljcos.hpp:63
Routines to calculate the generalized Lennard-Jones potential between particle pairs.
double ljgen_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Lennard-Jones energy.
Definition ljgen.hpp:81
Routines to calculate the Morse potential between particle pairs.
double morse_pair_energy(IA_parameters const &ia_params, double dist)
Calculate Morse energy.
Definition morse.hpp:53
Various procedures concerning interactions between particles.
Routines to calculate the energy and/or force for particle pairs via interpolation of lookup tables.
double tabulated_pair_energy(IA_parameters const &ia_params, double dist)
Calculate a non-bonded pair energy by linear interpolation from a table.
Routines to calculate the smooth step potential between particle pairs.
double SmSt_pair_energy(IA_parameters const &ia_params, double dist)
Calculate smooth step energy.
Routines to calculate the soft-sphere potential between particle pairs.
double soft_pair_energy(IA_parameters const &ia_params, double dist)
Calculate soft-sphere energy.
Exception indicating that a bond with an unexpected number of partners was encountered.
Exception indicating that a bond type was unknown.
Solver::ShortRangeEnergyKernel kernel_type
Solver::ShortRangeEnergyKernel kernel_type
Parameters for non-bonded interactions.
Struct holding all information for one particle.
Definition Particle.hpp:395
bool can_rotate() const
Definition Particle.hpp:460
auto const & rinertia() const
Definition Particle.hpp:502
auto is_virtual() const
Definition Particle.hpp:518
auto const & mass() const
Definition Particle.hpp:452
auto const & quat() const
Definition Particle.hpp:477
auto const & q() const
Definition Particle.hpp:508
auto const & v() const
Definition Particle.hpp:433
auto const & omega() const
Definition Particle.hpp:481
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
Routines to calculate the Thole damping potential between particle pairs.
double thole_pair_energy(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeEnergyKernel::kernel_type const *kernel)
Calculate Thole energy.
Definition thole.hpp:67
Routines to calculate the Weeks-Chandler-Andersen potential between particle pairs.
double wca_pair_energy(IA_parameters const &ia_params, double dist)
Calculate WCA energy.
Definition wca.hpp:49