ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
oif_global_forces.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#include "oif_global_forces.hpp"
21
22#include "BoxGeometry.hpp"
23#include "Particle.hpp"
26#include "communication.hpp"
28#include "system/System.hpp"
29
30#include <utils/Vector.hpp>
32
33#include <boost/mpi/collectives/all_reduce.hpp>
34#include <boost/serialization/utility.hpp>
35
36#include <cassert>
37#include <cmath>
38#include <functional>
39#include <span>
40
41/** Calculate the mesh volume and area. */
42static auto calc_oif_mesh(int molType, BoxGeometry const &box_geo,
43 CellStructure &cs,
44 BondedInteractionsMap const &bonded_ias) {
45
46 double area = 0.;
47 double volume = 0.;
48
49 cs.bond_loop([&area, &volume, &box_geo, &bonded_ias, molType](
50 Particle &p1, int bond_id, std::span<Particle *> partners) {
51 if (p1.mol_id() != molType)
52 return false;
53
54 if (boost::get<OifGlobalForcesBond>(bonded_ias.at(bond_id).get())) {
55 auto const p11 = box_geo.unfolded_position(p1.pos(), p1.image_box());
56 auto const p22 = p11 + box_geo.get_mi_vector(partners[0]->pos(), p11);
57 auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11);
58
59 auto const VOL_A = Utils::area_triangle(p11, p22, p33);
60 area += VOL_A;
61
62 auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33);
63 auto const VOL_dn = VOL_norm.norm();
64 auto const VOL_hz = (1.0 / 3.0) * (p11[2] + p22[2] + p33[2]);
65 volume -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz;
66 }
67
68 return false;
69 });
70
71 return Utils::Vector2d{{area, volume}};
72}
73
74/** Distribute the OIF global forces to all particles in the mesh. */
75static void add_oif_global_forces(double area, double volume, int molType,
76 BoxGeometry const &box_geo, CellStructure &cs,
77 BondedInteractionsMap const &bonded_ias) {
78
79 cs.bond_loop([&box_geo, &bonded_ias, area, volume, molType](
80 Particle &p1, int bond_id, std::span<Particle *> partners) {
81 if (p1.mol_id() != molType)
82 return false;
83
84 auto const *bond_ptr = bonded_ias.at(bond_id).get();
85 if (auto const *bond = boost::get<OifGlobalForcesBond>(bond_ptr)) {
86 auto const p11 = box_geo.unfolded_position(p1.pos(), p1.image_box());
87 auto const p22 = p11 + box_geo.get_mi_vector(partners[0]->pos(), p11);
88 auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11);
89
90 // unfolded positions correct
91 // starting code from volume force
92 auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33).normalize();
93 auto const VOL_A = Utils::area_triangle(p11, p22, p33);
94 auto const VOL_vv = (volume - bond->V0) / bond->V0;
95
96 auto const VOL_force = (1. / 3.) * bond->kv * VOL_vv * VOL_A * VOL_norm;
97
98 auto const h = (1. / 3.) * (p11 + p22 + p33);
99
100 auto const deltaA = (area - bond->A0_g) / bond->A0_g;
101
102 auto const m1 = h - p11;
103 auto const m2 = h - p22;
104 auto const m3 = h - p33;
105
106 auto const m1_length = m1.norm();
107 auto const m2_length = m2.norm();
108 auto const m3_length = m3.norm();
109
110 auto const fac = bond->ka_g * VOL_A * deltaA /
111 (m1_length * m1_length + m2_length * m2_length +
112 m3_length * m3_length);
113
114 p1.force() += fac * m1 + VOL_force;
115 partners[0]->force() += fac * m2 + VOL_force;
116 partners[1]->force() += fac * m3 + VOL_force;
117 }
118
119 return false;
120 });
121}
122
123void OifGlobal::run_force_kernel() {
124 auto &system = get_system();
125 auto &box_geo = *system.box_geo;
126 auto &bonded_ias = *system.bonded_ias;
127 auto &cell_structure = *system.cell_structure;
128 for (int i = 0; i < max_oif_objects; ++i) {
129 // There are two global quantities that need to be evaluated:
130 // object's surface and object's volume.
131 auto const local = calc_oif_mesh(i, box_geo, cell_structure, bonded_ias);
132 auto const global = boost::mpi::all_reduce(comm_cart, local, std::plus());
133 auto const area = std::abs(global[0]);
134 auto const volume = std::abs(global[1]);
135 if (area < 1e-100 and volume < 1e-100) {
136 break;
137 }
138 add_oif_global_forces(area, volume, i, box_geo, cell_structure, bonded_ias);
139 }
140}
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
container for bonded interactions.
mapped_type at(key_type const &key) const
auto unfolded_position(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Unfold particle coordinates to image box.
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.
Vector & normalize()
Definition Vector.hpp:147
boost::mpi::communicator comm_cart
The communicator.
double area_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the area of triangle between vectors P1,P2,P3, by computing the cross product P1P2 x P1P3 an...
Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the normal vector of a triangle.
static void add_oif_global_forces(double area, double volume, int molType, BoxGeometry const &box_geo, CellStructure &cs, BondedInteractionsMap const &bonded_ias)
Distribute the OIF global forces to all particles in the mesh.
static auto calc_oif_mesh(int molType, BoxGeometry const &box_geo, CellStructure &cs, BondedInteractionsMap const &bonded_ias)
Calculate the mesh volume and area.
Routines to calculate the OIF global forces for a particle triple (triangle from mesh).
Describes a cell structure / cell system.
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & image_box() const
Definition Particle.hpp:444
auto const & mol_id() const
Definition Particle.hpp:416
auto const & pos() const
Definition Particle.hpp:431
auto const & force() const
Definition Particle.hpp:435