ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
energy.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 "BoxGeometry.hpp"
23#include "Observable_stat.hpp"
25#include "constraints/Constraints.hpp"
26#include "energy_inline.hpp"
28#include "short_range_loop.hpp"
29#include "system/System.hpp"
30
33
34#include <memory>
35#include <span>
36
37namespace System {
38
39std::shared_ptr<Observable_stat> System::calculate_energy() {
40
41 auto obs_energy_ptr = std::make_shared<Observable_stat>(
42 1ul, static_cast<std::size_t>(bonded_ias->get_next_key()),
43 nonbonded_ias->get_max_seen_particle_type());
44
45 if (long_range_interactions_sanity_checks()) {
46 return obs_energy_ptr;
47 }
48
49 auto &obs_energy = *obs_energy_ptr;
50#if defined(CUDA) and (defined(ELECTROSTATICS) or defined(DIPOLES))
51 gpu.clear_energy_on_device();
52 gpu.update();
53#endif
54 on_observable_calc();
55
56 auto const local_parts = cell_structure->local_particles();
57
58 for (auto const &p : local_parts) {
59 obs_energy.kinetic[0] += calc_kinetic_energy(p);
60 }
61
62 auto const coulomb_kernel = coulomb.pair_energy_kernel();
63 auto const dipoles_kernel = dipoles.pair_energy_kernel();
64
66 [this, coulomb_kernel_ptr = get_ptr(coulomb_kernel), &obs_energy](
67 Particle const &p1, int bond_id, std::span<Particle *> partners) {
68 auto const &iaparams = *bonded_ias->at(bond_id);
69 auto const result = calc_bonded_energy(iaparams, p1, partners, *box_geo,
70 coulomb_kernel_ptr);
71 if (result) {
72 obs_energy.bonded_contribution(bond_id)[0] += result.value();
73 return false;
74 }
75 return true;
76 },
77 [coulomb_kernel_ptr = get_ptr(coulomb_kernel),
78 dipoles_kernel_ptr = get_ptr(dipoles_kernel), this,
79 &obs_energy](Particle const &p1, Particle const &p2, Distance const &d) {
80 auto const &ia_params =
81 nonbonded_ias->get_ia_param(p1.type(), p2.type());
82 add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2,
83 ia_params, *bonded_ias, coulomb_kernel_ptr,
84 dipoles_kernel_ptr, obs_energy);
85 },
86 *cell_structure, maximal_cutoff(), bonded_ias->maximal_cutoff());
87
88#ifdef ELECTROSTATICS
89 /* calculate k-space part of electrostatic interaction. */
90 obs_energy.coulomb[1] = coulomb.calc_energy_long_range(local_parts);
91#endif
92
93#ifdef DIPOLES
94 /* calculate k-space part of magnetostatic interaction. */
95 obs_energy.dipolar[1] = dipoles.calc_energy_long_range(local_parts);
96#endif
97
98 constraints->add_energy(local_parts, get_sim_time(), obs_energy);
99
100#if defined(CUDA) and (defined(ELECTROSTATICS) or defined(DIPOLES))
101 auto const energy_host = gpu.copy_energy_to_host();
102 if (!obs_energy.coulomb.empty())
103 obs_energy.coulomb[1] += static_cast<double>(energy_host.coulomb);
104 if (!obs_energy.dipolar.empty())
105 obs_energy.dipolar[1] += static_cast<double>(energy_host.dipolar);
106#endif
107
108 obs_energy.mpi_reduce();
109 return obs_energy_ptr;
110 // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
111}
112
114 if (cell_structure->get_resort_particles()) {
115 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
116 }
117
118 auto ret = 0.0;
119 if (auto const p = cell_structure->get_local_particle(pid)) {
120 auto const coulomb_kernel = coulomb.pair_energy_kernel();
121 auto kernel = [coulomb_kernel_ptr = get_ptr(coulomb_kernel), &ret,
122 this](Particle const &p, Particle const &p1,
123 Utils::Vector3d const &vec) {
124#ifdef EXCLUSIONS
125 if (not do_nonbonded(p, p1))
126 return;
127#endif
128 auto const &ia_params = nonbonded_ias->get_ia_param(p.type(), p1.type());
129 // Add energy for current particle pair to result
130 ret += calc_non_bonded_pair_energy(p, p1, ia_params, vec, vec.norm(),
131 *bonded_ias, coulomb_kernel_ptr);
132 };
133 cell_structure->run_on_particle_short_range_neighbors(*p, kernel);
134 }
135 return ret;
136}
137
138#ifdef DIPOLE_FIELD_TRACKING
140 dipoles.calc_long_range_field(cell_structure->local_particles());
141}
142#endif
143
144} // namespace System
double particle_short_range_energy_contribution(int pid)
Compute the short-range energy of a particle.
Definition energy.cpp:113
void calculate_long_range_fields()
Calculate dipole fields.
Definition energy.cpp:139
std::shared_ptr< Observable_stat > calculate_energy()
Calculate the total energy.
Definition energy.cpp:39
const T * get_ptr(std::optional< T > const &opt)
Energy calculation.
std::optional< double > calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, std::span< Particle * > partners, BoxGeometry const &box_geo, Coulomb::ShortRangeEnergyKernel::kernel_type const *kernel)
double calc_kinetic_energy(Particle const &p)
Calculate kinetic energies for one particle.
double calc_non_bonded_pair_energy(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_kernel)
Calculate non-bonded energies between a pair of particles.
void add_non_bonded_pair_energy(Particle const &p1, Particle const &p2, Utils::Vector3d const &d, double const dist, double const dist2, IA_parameters const &ia_params, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_kernel, Dipoles::ShortRangeEnergyKernel::kernel_type const *dipoles_kernel, Observable_stat &obs_energy)
Add non-bonded and short-range Coulomb energies between a pair of particles to the energy observable.
bool do_nonbonded(Particle const &p1, Particle const &p2)
Determine if the non-bonded interactions between p1 and p2 should be calculated.
Various procedures concerning interactions between particles.
void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, CellStructure &cell_structure, double pair_cutoff, double bond_cutoff, VerletCriterion const &verlet_criterion={})
Distance vector and length handed to pair kernels.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & type() const
Definition Particle.hpp:418