ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
debye_hueckel.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/**
23 * @file
24 * Calculate the Debye-Hückel energy and force for a particle pair.
25 */
26
27#pragma once
28
29#include "config/config.hpp"
30
31#ifdef ELECTROSTATICS
32
34
35#include "Particle.hpp"
36
37#include <utils/Vector.hpp>
39
40#include <cmath>
41#include <stdexcept>
42
43/** @brief Debye-Hückel parameters. */
44struct DebyeHueckel : public Coulomb::Actor<DebyeHueckel> {
45 /** @brief Ionic strength. */
46 double kappa;
47 /** @brief Interaction cutoff. */
48 double r_cut;
49
50 DebyeHueckel(double prefactor, double kappa, double r_cut) {
51 if (kappa < 0.0) {
52 throw std::domain_error("Parameter 'kappa' must be >= 0");
53 }
54 if (r_cut < 0.0) {
55 throw std::domain_error("Parameter 'r_cut' must be >= 0");
56 }
57
59 this->kappa = kappa;
60 this->r_cut = r_cut;
61 }
62
63 void on_activation() const { sanity_checks(); }
64 void on_boxl_change() const {}
65 void on_node_grid_change() const {}
66 void on_periodicity_change() const {}
68 void init() const {}
69
71
72 /** @brief Compute the pair force.
73 * @param[in] q1q2 Product of the charges on p1 and p2.
74 * @param[in] d Vector pointing from p1 to p2.
75 * @param[in] dist Distance between p1 and p2.
76 */
77 Utils::Vector3d pair_force(double const q1q2, Utils::Vector3d const &d,
78 double const dist) const {
79 if (dist >= r_cut) {
80 return {};
81 }
82 // pure Coulomb case
83 auto fac = prefactor * q1q2 / Utils::int_pow<3>(dist);
84 if (kappa > 0.) {
85 // Debye-Hueckel case
86 fac *= std::exp(-kappa * dist) * (1. + kappa * dist);
87 }
88 return fac * d;
89 }
90
91 /** @brief Compute the pair energy.
92 * @param q1q2 Product of the charges on p1 and p2.
93 * @param dist Distance between p1 and p2.
94 */
95 double pair_energy(double const q1q2, double const dist) const {
96 if (dist >= r_cut) {
97 return 0.;
98 }
99 // pure Coulomb case
100 auto energy = prefactor * q1q2 / dist;
101 if (kappa > 0.) {
102 // Debye-Hueckel case
103 energy *= std::exp(-kappa * dist);
104 }
105 return energy;
106 }
107};
108
109#endif // ELECTROSTATICS
Vector implementation and trait types for boost qvm interoperability.
void set_prefactor(double new_prefactor)
double prefactor
Electrostatics prefactor.
This file contains the defaults for ESPResSo.
Debye-Hückel parameters.
double kappa
Ionic strength.
void on_node_grid_change() const
void on_activation() const
void on_periodicity_change() const
double r_cut
Interaction cutoff.
double pair_energy(double const q1q2, double const dist) const
Compute the pair energy.
void on_boxl_change() const
DebyeHueckel(double prefactor, double kappa, double r_cut)
void on_cell_structure_change() const
void init() const
Utils::Vector3d pair_force(double const q1q2, Utils::Vector3d const &d, double const dist) const
Compute the pair force.
void sanity_checks() const