ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
magnetostatics/scafacos_impl.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 SCAFACOS_DIPOLES
25
28
29#include "BoxGeometry.hpp"
31#include "communication.hpp"
32#include "system/System.hpp"
33
34#include <utils/Vector.hpp>
35#include <utils/matrix.hpp>
36
37#include <cassert>
38#include <iterator>
39#include <memory>
40#include <span>
41#include <string>
42
43std::shared_ptr<DipolarScafacos>
44make_dipolar_scafacos(std::string const &method,
45 std::string const &parameters) {
46 return std::make_shared<DipolarScafacosImpl>(comm_cart, method, parameters);
47}
48
50 auto const &system = get_system();
51 auto const &box_geo = *system.box_geo;
52 auto const &cell_structure = *system.cell_structure;
53
54 positions.clear();
55 dipoles.clear();
56
57 for (auto const &p : cell_structure.local_particles()) {
58 auto const pos = box_geo.folded_position(p.pos());
59 positions.push_back(pos[0]);
60 positions.push_back(pos[1]);
61 positions.push_back(pos[2]);
62 auto const dip = p.calc_dip();
63 dipoles.push_back(dip[0]);
64 dipoles.push_back(dip[1]);
65 dipoles.push_back(dip[2]);
66 }
67}
68
70 if (positions.empty())
71 return;
72
73 auto const &cell_structure = *get_system().cell_structure;
74
75 auto it_potentials = potentials.begin();
76 auto index = std::size_t{0ul};
77 for (auto &p : cell_structure.local_particles()) {
78 // The scafacos term "potential" here in fact refers to the magnetic
79 // field. So, the torques are given by m \times B
80 auto const dip = p.calc_dip();
81 auto const t = vector_product(
82 dip, Utils::Vector3d(std::span<const double>(&*it_potentials, 3ul)));
83 // The force is given by G m, where G is a matrix
84 // which comes from the "fields" output of scafacos like this
85 // 0 1 2
86 // 1 3 4
87 // 2 4 5
88 // where the numbers refer to indices in the "field" output from scafacos
89 auto const G = Utils::Matrix<double, 3, 3>{
90 {fields[index + 0ul], fields[index + 1ul], fields[index + 2ul]},
91 {fields[index + 1ul], fields[index + 3ul], fields[index + 4ul]},
92 {fields[index + 2ul], fields[index + 4ul], fields[index + 5ul]}};
93 auto const f = G * dip;
94
95 // Add to particles
96 p.force() += prefactor * f;
97 p.torque() += prefactor * t;
98 index += 6ul;
99 std::advance(it_potentials, 3);
100 }
101
102 /* Check that the particle number did not change */
103 assert(it_potentials == potentials.end());
104}
105
106#endif // SCAFACOS_DIPOLES
Vector implementation and trait types for boost qvm interoperability.
double prefactor
Magnetostatics prefactor.
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
__device__ void vector_product(float const *a, float const *b, float *out)
std::shared_ptr< DipolarScafacos > make_dipolar_scafacos(std::string const &method, std::string const &parameters)
Matrix implementation and trait types for boost qvm interoperability.
void update_particle_forces() const override
Matrix representation with static size.
Definition matrix.hpp:65