ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
harmonic.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 * Routines to calculate the harmonic bond potential between particle pairs.
26 */
27
28#include "config/config.hpp"
29#include "errorhandling.hpp"
30
31#include <utils/Vector.hpp>
32#include <utils/math/sqr.hpp>
33
34#include <optional>
35
36/** Parameters for harmonic bond Potential */
38 /** spring constant */
39 double k;
40 /** equilibrium bond length */
41 double r;
42 /** cutoff length */
43 double r_cut;
44
45 double cutoff() const { return r_cut; }
46
47 static constexpr int num = 1;
48
49 HarmonicBond(double k, double r, double r_cut) {
50 this->k = k;
51 this->r = r;
52 this->r_cut = r_cut;
53 }
54
55 std::optional<Utils::Vector3d> force(Utils::Vector3d const &dx) const;
56 std::optional<double> energy(Utils::Vector3d const &dx) const;
57};
58
59/** Compute the harmonic bond force.
60 * @param[in] dx Distance between the particles.
61 */
62inline std::optional<Utils::Vector3d>
64 auto const dist = dx.norm();
65
66 if ((r_cut > 0.0) && (dist > r_cut)) {
67 return {};
68 }
69
70 auto const dr = dist - r;
71 auto fac = -k * dr;
72 if (dist > ROUND_ERROR_PREC) { /* Regular case */
73 fac /= dist;
74 } else {
75 if (r > 0.) {
76 runtimeErrorMsg() << "Harmonic bond: Particles have zero distance. "
77 "This is most likely an error in the system setup.";
78 }
79 }
80 return fac * dx;
81}
82
83/** Compute the harmonic bond energy.
84 * @param[in] dx Distance between the particles.
85 */
86inline std::optional<double>
88 auto const dist = dx.norm();
89
90 if ((r_cut > 0.0) && (dist > r_cut)) {
91 return {};
92 }
93
94 return 0.5 * k * Utils::sqr(dist - r);
95}
Vector implementation and trait types for boost qvm interoperability.
T norm() const
Definition Vector.hpp:138
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
Definition config.hpp:66
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Parameters for harmonic bond Potential.
Definition harmonic.hpp:37
std::optional< Utils::Vector3d > force(Utils::Vector3d const &dx) const
Compute the harmonic bond force.
Definition harmonic.hpp:63
double cutoff() const
Definition harmonic.hpp:45
static constexpr int num
Definition harmonic.hpp:47
std::optional< double > energy(Utils::Vector3d const &dx) const
Compute the harmonic bond energy.
Definition harmonic.hpp:87
double r
equilibrium bond length
Definition harmonic.hpp:41
double r_cut
cutoff length
Definition harmonic.hpp:43
double k
spring constant
Definition harmonic.hpp:39
HarmonicBond(double k, double r, double r_cut)
Definition harmonic.hpp:49