ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
nonbonded_interaction_data.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 * Various procedures concerning interactions between particles.
26 */
27
29#include "config/config.hpp"
30#include "system/Leaf.hpp"
31
32#include <utils/index.hpp>
34
35#include <algorithm>
36#include <cassert>
37#include <cmath>
38#include <memory>
39#include <vector>
40
41/** Cutoff for deactivated interactions. Must be negative, so that even
42 * particles on top of each other don't interact by chance.
43 */
44constexpr double INACTIVE_CUTOFF = -1.;
45
46/** Lennard-Jones with shift */
48 double eps = 0.0;
49 double sig = 0.0;
51 double shift = 0.0;
52 double offset = 0.0;
53 double min = 0.0;
54 LJ_Parameters() = default;
55 LJ_Parameters(double epsilon, double sigma, double cutoff, double offset,
56 double min, double shift);
57 double get_auto_shift() const {
58 auto auto_shift = 0.;
59 if (cut != 0.) {
60 auto_shift = Utils::int_pow<6>(sig / cut) - Utils::int_pow<12>(sig / cut);
61 }
62 return auto_shift;
63 }
64 double max_cutoff() const { return cut + offset; }
65 double min_cutoff() const { return min + offset; }
66};
67
68/** WCA potential */
70 double eps = 0.0;
71 double sig = 0.0;
73 WCA_Parameters() = default;
74 WCA_Parameters(double epsilon, double sigma);
75 double max_cutoff() const { return cut; }
76};
77
78/** Generic Lennard-Jones with shift */
80 double eps = 0.0;
81 double sig = 0.0;
83 double shift = 0.0;
84 double offset = 0.0;
85 double lambda = 1.0;
86 double softrad = 0.0;
87 double a1 = 0.0;
88 double a2 = 0.0;
89 double b1 = 0.0;
90 double b2 = 0.0;
91 LJGen_Parameters() = default;
92 LJGen_Parameters(double epsilon, double sigma, double cutoff, double shift,
93 double offset,
94#ifdef LJGEN_SOFTCORE
95 double lam, double delta,
96#endif
97 double e1, double e2, double b1, double b2);
98 double get_auto_shift() const {
99 auto auto_shift = 0.;
100 if (cut != 0.) {
101 auto_shift = b2 * std::pow(sig / cut, a2) - b1 * std::pow(sig / cut, a1);
102 }
103 return auto_shift;
104 }
105 double max_cutoff() const { return cut + offset; }
106};
107
108/** smooth step potential */
110 double eps = 0.0;
111 double sig = 0.0;
113 double d = 0.0;
114 int n = 0;
115 double k0 = 0.0;
117 SmoothStep_Parameters(double eps, double sig, double cutoff, double d, int n,
118 double k0);
119 double max_cutoff() const { return cut; }
120};
121
122/** Hertzian potential */
124 double eps = 0.0;
127 Hertzian_Parameters(double eps, double sig);
128 double max_cutoff() const { return sig; }
129};
130
131/** Gaussian potential */
133 double eps = 0.0;
134 double sig = 1.0;
137 Gaussian_Parameters(double eps, double sig, double cutoff);
138 double max_cutoff() const { return cut; }
139};
140
141/** BMHTF NaCl potential */
143 double A = 0.0;
144 double B = 0.0;
145 double C = 0.0;
146 double D = 0.0;
147 double sig = 0.0;
149 double computed_shift = 0.0;
150 BMHTF_Parameters() = default;
151 BMHTF_Parameters(double A, double B, double C, double D, double sig,
152 double cut);
153 double max_cutoff() const { return cut; }
154};
155
156/** Morse potential */
158 double eps = 0.;
163 Morse_Parameters() = default;
164 Morse_Parameters(double eps, double alpha, double rmin, double cutoff);
165 double max_cutoff() const { return cut; }
166};
167
168/** Buckingham potential */
170 double A = 0.0;
171 double B = 0.0;
172 double C = 0.0;
173 double D = 0.0;
175 double discont = 0.0;
176 double shift = 0.0;
177 double F1 = 0.0;
178 double F2 = 0.0;
180 Buckingham_Parameters(double a, double b, double c, double d, double cutoff,
181 double discont, double shift);
182 double max_cutoff() const { return cut; }
183};
184
185/** soft-sphere potential */
187 double a = 0.0;
188 double n = 0.0;
190 double offset = 0.0;
192 SoftSphere_Parameters(double a, double n, double cutoff, double offset);
193 double max_cutoff() const { return cut + offset; }
194};
195
196/** hat potential */
198 double Fmax = 0.0;
200 Hat_Parameters() = default;
201 Hat_Parameters(double F_max, double cutoff);
202 double max_cutoff() const { return r; }
203};
204
205/** Lennard-Jones+Cos potential */
207 double eps = 0.0;
208 double sig = 0.0;
210 double offset = 0.0;
211 double alfa = 0.0;
212 double beta = 0.0;
213 double rmin = 0.0;
214 LJcos_Parameters() = default;
215 LJcos_Parameters(double epsilon, double sigma, double cutoff, double offset);
216 double max_cutoff() const { return cut + offset; }
217};
218
219/** Lennard-Jones with a different Cos potential */
221 double eps = 0.0;
222 double sig = 0.0;
224 double offset = 0.0;
225 double w = 0.0;
226 double rchange = 0.0;
227 LJcos2_Parameters() = default;
228 LJcos2_Parameters(double epsilon, double sigma, double offset, double width);
229 double max_cutoff() const { return cut + offset; }
230};
231
232/** Gay-Berne potential */
234 double eps = 0.0;
235 double sig = 0.0;
237 double k1 = 0.0;
238 double k2 = 0.0;
239 double mu = 0.0;
240 double nu = 0.0;
241 double chi1 = 0.0;
242 double chi2 = 0.0;
244 GayBerne_Parameters(double eps, double sig, double cut, double k1, double k2,
245 double mu, double nu);
246 double max_cutoff() const { return cut; }
247};
248
249/** Thole potential */
251 double scaling_coeff = 0.; // inactive cutoff is 0
252 double q1q2 = 0.;
253 Thole_Parameters() = default;
256};
257
258/** DPD potential */
260 double gamma = 0.;
261 double k = 1.;
263 int wf = 0;
264 double pref = 0.0;
265};
266
270 DPD_Parameters() = default;
271 DPD_Parameters(double gamma, double k, double r_c, int wf, double tgamma,
272 double tr_c, int twf) {
273 radial = DPDParameters{gamma, k, r_c, wf, -1.};
274 trans = DPDParameters{tgamma, k, tr_c, twf, -1.};
275 }
276 double max_cutoff() const { return std::max(radial.cutoff, trans.cutoff); }
277};
278
279/** @brief Parameters for non-bonded interactions. */
281 /** maximal cutoff for this pair of particle types. This contains
282 * contributions from the short-ranged interactions, plus any
283 * cutoffs from global interactions like electrostatics.
284 */
286
287#ifdef LENNARD_JONES
289#endif
290
291#ifdef WCA
293#endif
294
295#ifdef LENNARD_JONES_GENERIC
297#endif
298
299#ifdef SMOOTH_STEP
301#endif
302
303#ifdef HERTZIAN
305#endif
306
307#ifdef GAUSSIAN
309#endif
310
311#ifdef BMHTF_NACL
313#endif
314
315#ifdef MORSE
317#endif
318
319#ifdef BUCKINGHAM
321#endif
322
323#ifdef SOFT_SPHERE
325#endif
326
327#ifdef HAT
329#endif
330
331#ifdef LJCOS
333#endif
334
335#ifdef LJCOS2
337#endif
338
339#ifdef GAY_BERNE
341#endif
342
343#ifdef TABULATED
345#endif
346
347#ifdef DPD
349#endif
350
351#ifdef THOLE
353#endif
354};
355
356class InteractionsNonBonded : public System::Leaf<InteractionsNonBonded> {
357 /** @brief List of pairwise interactions. */
358 std::vector<std::shared_ptr<IA_parameters>> m_nonbonded_ia_params{};
359 /** @brief Maximal particle type seen so far. */
360 int max_seen_particle_type = -1;
361
362 void realloc_ia_params(int type) {
363 assert(type >= 0);
364 auto const old_size = m_nonbonded_ia_params.size();
365 m_nonbonded_ia_params.resize(Utils::lower_triangular(type, type) + 1);
366 auto const new_size = m_nonbonded_ia_params.size();
367 if (new_size > old_size) {
368 for (auto &data : m_nonbonded_ia_params) {
369 if (data == nullptr) {
370 data = std::make_shared<IA_parameters>();
371 }
372 }
373 }
374 }
375
376public:
378 /* make sure interaction 0<->0 always exists */
380 }
381
382 /**
383 * @brief Make sure the interaction parameter list is large enough to cover
384 * interactions for this particle type.
385 * New interactions are initialized with values such that no physical
386 * interaction occurs.
387 */
389 assert(type >= 0);
390 if (type > max_seen_particle_type) {
391 realloc_ia_params(type);
392 max_seen_particle_type = type;
393 }
394 }
395
396 auto get_ia_param_key(int i, int j) const {
397 assert(i >= 0 and i <= max_seen_particle_type);
398 assert(j >= 0 and j <= max_seen_particle_type);
399 auto const key = static_cast<unsigned int>(
400 Utils::lower_triangular(std::max(i, j), std::min(i, j)));
401 assert(key < m_nonbonded_ia_params.size());
402 return key;
403 }
404
405 /**
406 * @brief Get interaction parameters between particle types i and j
407 *
408 * This is symmetric, e.g. it holds that `get_ia_param(i, j)` and
409 * `get_ia_param(j, i)` point to the same data.
410 *
411 * @param i First type, must exist
412 * @param j Second type, must exist
413 *
414 * @return Reference to interaction parameters for the type pair.
415 */
416 auto &get_ia_param(int i, int j) {
417 return *m_nonbonded_ia_params[get_ia_param_key(i, j)];
418 }
419
420 auto const &get_ia_param(int i, int j) const {
421 return *m_nonbonded_ia_params[get_ia_param_key(i, j)];
422 }
423
424 auto get_ia_param_ref_counted(int i, int j) const {
425 return m_nonbonded_ia_params[get_ia_param_key(i, j)];
426 }
427
428 void set_ia_param(int i, int j, std::shared_ptr<IA_parameters> const &ia) {
429 m_nonbonded_ia_params[get_ia_param_key(i, j)] = ia;
430 }
431
432 auto get_max_seen_particle_type() const { return max_seen_particle_type; }
433
434 /** @brief Recalculate cutoff of each interaction struct. */
436
437 /** @brief Get maximal cutoff. */
438 double maximal_cutoff() const;
439
440 /** @brief Notify system that non-bonded interactions changed. */
441 void on_non_bonded_ia_change() const;
442};
auto get_ia_param_ref_counted(int i, int j) const
auto get_ia_param_key(int i, int j) const
void set_ia_param(int i, int j, std::shared_ptr< IA_parameters > const &ia)
auto const & get_ia_param(int i, int j) const
void recalc_maximal_cutoffs()
Recalculate cutoff of each interaction struct.
void on_non_bonded_ia_change() const
Notify system that non-bonded interactions changed.
double maximal_cutoff() const
Get maximal cutoff.
auto & get_ia_param(int i, int j)
Get interaction parameters between particle types i and j.
void make_particle_type_exist(int type)
Make sure the interaction parameter list is large enough to cover interactions for this particle type...
Abstract class that represents a component of the system.
This file contains the defaults for ESPResSo.
T lower_triangular(T i, T j)
Linear index into a lower triangular matrix.
Definition index.hpp:71
constexpr double INACTIVE_CUTOFF
Cutoff for deactivated interactions.
BMHTF_Parameters()=default
Buckingham_Parameters()=default
DPD_Parameters()=default
DPD_Parameters(double gamma, double k, double r_c, int wf, double tgamma, double tr_c, int twf)
Gaussian_Parameters()=default
GayBerne_Parameters()=default
Hat_Parameters()=default
Hertzian_Parameters()=default
Parameters for non-bonded interactions.
Gaussian_Parameters gaussian
double max_cut
maximal cutoff for this pair of particle types.
GayBerne_Parameters gay_berne
SoftSphere_Parameters soft_sphere
SmoothStep_Parameters smooth_step
Hertzian_Parameters hertzian
Buckingham_Parameters buckingham
Generic Lennard-Jones with shift.
LJGen_Parameters()=default
Lennard-Jones with shift.
LJ_Parameters()=default
Lennard-Jones with a different Cos potential.
LJcos2_Parameters()=default
Lennard-Jones+Cos potential.
LJcos_Parameters()=default
Morse_Parameters()=default
SmoothStep_Parameters()=default
SoftSphere_Parameters()=default
Evaluate forces and energies using a custom potential profile.
Thole_Parameters()=default
Thole_Parameters(double scaling_coeff, double q1q2)
WCA_Parameters()=default