ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
fene.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 FENE 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 <cmath>
35#include <optional>
36
37/** Parameters for FENE bond Potential. */
38struct FeneBond {
39 /** spring constant */
40 double k;
41 /** maximal bond stretching */
42 double drmax;
43 /** equilibrium bond length */
44 double r0;
45 /** square of @p drmax (internal parameter) */
46 double drmax2;
47 /** inverse square of @p drmax (internal parameter) */
48 double drmax2i;
49
50 double cutoff() const { return r0 + drmax; }
51
52 static constexpr int num = 1;
53
54 FeneBond(double k, double drmax, double r0) {
55 this->k = k;
56 this->drmax = drmax;
57 this->r0 = r0;
58
59 this->drmax2 = Utils::sqr(this->drmax);
60 this->drmax2i = 1. / this->drmax2;
61 }
62
63 std::optional<Utils::Vector3d> force(Utils::Vector3d const &dx) const;
64 std::optional<double> energy(Utils::Vector3d const &dx) const;
65};
66
67/** Compute the FENE bond force.
68 * @param[in] dx Distance between the particles.
69 */
70inline std::optional<Utils::Vector3d>
72 auto const len = dx.norm();
73 auto const dr = len - r0;
74
75 if (dr >= drmax) {
76 return {};
77 }
78
79 auto fac = -k * dr / (1.0 - dr * dr * drmax2i);
80 if (len > ROUND_ERROR_PREC) {
81 fac /= len;
82 } else {
83 if (r0 > 0.) {
84 runtimeErrorMsg() << "FENE bond: Particles have zero distance. "
85 "This is most likely an error in the system setup.";
86 }
87 }
88
89 return fac * dx;
90}
91
92/** Compute the FENE bond energy.
93 * @param[in] dx Distance between the particles.
94 */
95inline std::optional<double> FeneBond::energy(Utils::Vector3d const &dx) const {
96 /* compute bond stretching (r-r0) */
97 double const dr = dx.norm() - r0;
98
99 /* check bond stretching */
100 if (dr >= drmax) {
101 return {};
102 }
103
104 return -0.5 * k * drmax2 * log(1.0 - dr * dr * drmax2i);
105}
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 FENE bond Potential.
Definition fene.hpp:38
double drmax2i
inverse square of drmax (internal parameter)
Definition fene.hpp:48
double drmax
maximal bond stretching
Definition fene.hpp:42
FeneBond(double k, double drmax, double r0)
Definition fene.hpp:54
double drmax2
square of drmax (internal parameter)
Definition fene.hpp:46
double cutoff() const
Definition fene.hpp:50
static constexpr int num
Definition fene.hpp:52
double r0
equilibrium bond length
Definition fene.hpp:44
std::optional< double > energy(Utils::Vector3d const &dx) const
Compute the FENE bond energy.
Definition fene.hpp:95
std::optional< Utils::Vector3d > force(Utils::Vector3d const &dx) const
Compute the FENE bond force.
Definition fene.hpp:71
double k
spring constant
Definition fene.hpp:40