ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
steepest_descent.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
23
24#include "Particle.hpp"
25#include "ParticleRange.hpp"
27#include "communication.hpp"
28#include "config/config.hpp"
29#include "rotation.hpp"
30#include "system/System.hpp"
31
32#include <utils/Vector.hpp>
33#include <utils/math/sqr.hpp>
34
35#include <boost/mpi/collectives/all_reduce.hpp>
36#include <boost/mpi/operations.hpp>
37
38#include <algorithm>
39#include <cmath>
40#include <limits>
41#include <stdexcept>
42
43/** Currently active steepest descent instance */
45
46bool steepest_descent_step(const ParticleRange &particles) {
47 // Maximal force encountered on node
48 auto f_max = -std::numeric_limits<double>::max();
49
50 // Iteration over all local particles
51 for (auto &p : particles) {
52 auto f = 0.0;
53
54 // For all Cartesian coordinates
55 for (auto j = 0u; j < 3u; ++j) {
56 // Skip, if coordinate is fixed
57 if (!p.is_fixed_along(j)) {
58 // Square of force on particle
59 f += Utils::sqr(p.force()[j]);
60
61 // Positional increment, crop to maximum allowed by user
62 auto const dp =
63 std::clamp(params.gamma * p.force()[j], -params.max_displacement,
65
66 // Move particle
67 p.pos()[j] += dp;
68 }
69 }
70#ifdef ROTATION
71 {
72 // Rotational increment
73 auto const dq = params.gamma * p.torque(); // Vector parallel to torque
74 auto const t = p.torque().norm2();
75
76 // Normalize rotation axis and compute amount of rotation
77 auto const l = dq.norm();
78 if (l > 0.0) {
79 auto const axis = dq / l;
80 auto const angle =
82
83 // Rotate the particle around axis dq by amount l
84 local_rotate_particle(p, axis, angle);
85 }
86
87 f_max = std::max(f_max, t);
88 }
89#endif
90 // Note maximum force/torque encountered
91 f_max = std::max(f_max, f);
92 }
93
94 auto &cell_structure = *System::get_system().cell_structure;
95 cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
96
97 // Synchronize maximum force/torque encountered
98 auto const f_max_global =
99 boost::mpi::all_reduce(comm_cart, f_max, boost::mpi::maximum<double>());
100
101 return sqrt(f_max_global) < params.f_max;
102}
103
105 ::params = obj;
106}
107
109 const double f_max, const double gamma, const double max_displacement)
110 : f_max{f_max}, gamma{gamma}, max_displacement{max_displacement} {
111 if (f_max < 0.0) {
112 throw std::runtime_error("The maximal force must be positive.");
113 }
114 if (gamma < 0.0) {
115 throw std::runtime_error("The dampening constant must be positive.");
116 }
117 if (max_displacement < 0.0) {
118 throw std::runtime_error("The maximal displacement must be positive.");
119 }
120}
Vector implementation and trait types for boost qvm interoperability.
A range of particles.
std::shared_ptr< CellStructure > cell_structure
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
System & get_system()
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
This file contains all subroutines required to process rotational motion.
void local_rotate_particle(Particle &p, const Utils::Vector3d &axis_space_frame, const double phi)
Rotate the particle p around the NORMALIZED axis aSpaceFrame by amount phi.
Definition rotation.hpp:133
bool steepest_descent_step(const ParticleRange &particles)
Steepest descent integrator.
void register_integrator(SteepestDescentParameters const &obj)
static SteepestDescentParameters params
Currently active steepest descent instance.
Parameters for the steepest descent algorithm.
double gamma
Dampening constant.
double max_displacement
Maximal particle displacement.
double f_max
Maximal particle force.
SteepestDescentParameters(double f_max, double gamma, double max_displacement)