ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
thermalized_bond_kernel.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#include "thermalized_bond.hpp"
25
26#include "Particle.hpp"
27#include "random.hpp"
28#include "thermostat.hpp"
29
30#include <utils/Vector.hpp>
31
32#include <cmath>
33#include <optional>
34#include <tuple>
35
36/** Separately thermalizes the com and distance of a particle pair.
37 * @param[in] p1 First particle.
38 * @param[in] p2 Second particle.
39 * @param[in] dx Distance between the particles.
40 * @return the forces on @p p1 and @p p2
41 */
42inline std::optional<std::tuple<Utils::Vector3d, Utils::Vector3d>>
44 Utils::Vector3d const &dx) const {
45 // Bond broke?
46 if (r_cut > 0.0 && dx.norm() > r_cut) {
47 return {};
48 }
49
50 auto const mass_tot = p1.mass() + p2.mass();
51 auto const mass_tot_inv = 1.0 / mass_tot;
52 auto const sqrt_mass_tot = sqrt(mass_tot);
53 auto const sqrt_mass_red = sqrt(p1.mass() * p2.mass() / mass_tot);
54 auto const com_vel = mass_tot_inv * (p1.mass() * p1.v() + p2.mass() * p2.v());
55 auto const dist_vel = p2.v() - p1.v();
56 auto const thermostat_view = m_thermostat.lock();
57 assert(thermostat_view);
58 auto const &thermalized_bond = *thermostat_view->thermalized_bond;
59
60 Utils::Vector3d force1{};
61 Utils::Vector3d force2{};
62 auto const noise = Random::noise_uniform<RNGSalt::THERMALIZED_BOND>(
63 thermalized_bond.rng_counter(), thermalized_bond.rng_seed(), p1.id(),
64 p2.id());
65
66 for (unsigned int i = 0u; i < 3u; ++i) {
67 double force_lv_com, force_lv_dist;
68
69 // Langevin thermostat for center of mass
70 if (pref2_com > 0.0) {
71 force_lv_com =
72 -pref1_com * com_vel[i] + sqrt_mass_tot * pref2_com * noise[i];
73 } else {
74 force_lv_com = -pref1_com * com_vel[i];
75 }
76
77 // Langevin thermostat for distance p1->p2
78 if (pref2_dist > 0.0) {
79 force_lv_dist =
80 -pref1_dist * dist_vel[i] + sqrt_mass_red * pref2_dist * noise[i];
81 } else {
82 force_lv_dist = -pref1_dist * dist_vel[i];
83 }
84 // Add forces
85 force1[i] = p1.mass() * mass_tot_inv * force_lv_com - force_lv_dist;
86 force2[i] = p2.mass() * mass_tot_inv * force_lv_com + force_lv_dist;
87 }
88
89 return std::make_tuple(force1, force2);
90}
Vector implementation and trait types for boost qvm interoperability.
T norm() const
Definition Vector.hpp:138
Random number generation using Philox.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & mass() const
Definition Particle.hpp:452
auto const & v() const
Definition Particle.hpp:433
auto const & id() const
Definition Particle.hpp:414
std::optional< std::tuple< Utils::Vector3d, Utils::Vector3d > > forces(Particle const &p1, Particle const &p2, Utils::Vector3d const &dx) const
Separately thermalizes the com and distance of a particle pair.
Routines to thermalize the center of mass and distance of a particle pair.