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