ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
virtual_sites.cpp
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#include "config/config.hpp"
23
24#ifdef VIRTUAL_SITES_RELATIVE
25
26#include "virtual_sites.hpp"
27
28#include "BoxGeometry.hpp"
29#include "Particle.hpp"
30#include "communication.hpp"
31#include "errorhandling.hpp"
33
34#include <utils/Vector.hpp>
36#include <utils/quaternion.hpp>
37
38#include <cmath>
39#include <cstdio>
40#include <tuple>
41
42std::tuple<Utils::Quaternion<double>, double>
43calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to,
44 BoxGeometry const &box_geo, double min_global_cut,
45 bool override_cutoff_check) {
46 // get the distance between the particles
47 auto d = box_geo.get_mi_vector(p_vs.pos(), p_relate_to.pos());
48
49 // Check if the distance between virtual and non-virtual particles is larger
50 // than minimum global cutoff. If so, warn user.
51 auto const dist = d.norm();
52 if (dist > min_global_cut and ::communicator.size > 1 and
53 not override_cutoff_check) {
55 << "Warning: The distance between virtual and non-virtual particle ("
56 << dist << ") is larger than the minimum global cutoff ("
57 << min_global_cut << "). This may lead to incorrect simulations under "
58 << "certain conditions. Adjust the property system.min_global_cut to "
59 << "increase the minimum cutoff.";
60 }
61
62 // If the distance between real & virtual particle is 0 we just set the
63 // relative orientation to {1 0 0 0}, as it is irrelevant but needs to be
64 // a valid quaternion
65 if (dist == 0.) {
66 return std::make_tuple(Utils::Quaternion<double>::identity(), dist);
67 }
68
69 // Now, calculate the quaternions which specify the angle between
70 // the director of the particle we relate to and the vector
71 // (particle_we_relate_to - this_particle)
72 // The vs_relative implementation later obtains the director by multiplying
73 // the quaternions representing the orientation of the real particle
74 // with those in the virtual particle. The resulting quaternion is then
75 // converted to a director.
76 // We have quat_(real particle) * quat_(virtual particle)
77 // = quat_(obtained from desired director)
78 // Resolving this for the quat_(virtual particle)
79
80 d.normalize();
81
82 // Obtain quaternion from desired director
83 Utils::Quaternion<double> quat_director =
85
86 // Define quaternion as described above
87 auto relate_to_quat = p_relate_to.quat();
88 auto quat =
89 Utils::Quaternion<double>{{{{Utils::dot(relate_to_quat, quat_director),
90 -quat_director[0] * relate_to_quat[1] +
91 quat_director[1] * relate_to_quat[0] +
92 quat_director[2] * relate_to_quat[3] -
93 quat_director[3] * relate_to_quat[2],
94 relate_to_quat[1] * quat_director[3] +
95 relate_to_quat[0] * quat_director[2] -
96 relate_to_quat[3] * quat_director[1] -
97 relate_to_quat[2] * quat_director[0],
98 quat_director[3] * relate_to_quat[0] -
99 relate_to_quat[3] * quat_director[0] +
100 relate_to_quat[2] * quat_director[1] -
101 relate_to_quat[1] * quat_director[2]}}}};
102 quat /= relate_to_quat.norm2();
103
104 // Verify result
105 Utils::Quaternion<double> qtemp = relate_to_quat * quat;
106 for (unsigned int i = 0; i < 4; i++) {
107 if (fabs(qtemp[i] - quat_director[i]) > 1E-9) {
108 fprintf(stderr, "vs_relate_to: component %u: %f instead of %f\n", i,
109 qtemp[i], quat_director[i]);
110 }
111 }
112 return std::make_tuple(quat, dist);
113}
114
115#endif // VIRTUAL_SITES_RELATIVE
Vector implementation and trait types for boost qvm interoperability.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
Quaternion algebra.
Quaternion< T > convert_director_to_quaternion(Vector< T, 3 > const &d)
Convert director to quaternion.
Various procedures concerning interactions between particles.
Quaternion implementation and trait types for boost qvm interoperability.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & quat() const
Definition Particle.hpp:477
auto const & pos() const
Definition Particle.hpp:431
Quaternion representation.
value_type norm2() const
Retrieve the square of the norm of the quaternion.
std::tuple< Utils::Quaternion< double >, double > calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, BoxGeometry const &box_geo, double min_global_cut, bool override_cutoff_check)
Calculate the rotation quaternion and distance between two particles.