Skip to content
Permalink
Browse files

Merge branch 'python' into local_cells_global

  • Loading branch information
kodiakhq committed Jan 8, 2020
2 parents 6886c49 + baa014f commit a6b46573738070cf12fcced1c146b44eadca9eb9
Showing with 315 additions and 375 deletions.
  1. +0 −1 src/core/communication.cpp
  2. +8 −13 src/core/diamond.cpp
  3. +0 −2 src/core/diamond.hpp
  4. +2 −1 src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp
  5. +1 −2 src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp
  6. +6 −5 src/core/observables/BondAngles.hpp
  7. +7 −7 src/core/observables/BondDihedrals.hpp
  8. +8 −11 src/core/observables/ComForce.hpp
  9. +9 −8 src/core/observables/ComPosition.hpp
  10. +9 −8 src/core/observables/ComVelocity.hpp
  11. +4 −3 src/core/observables/CosPersistenceAngles.hpp
  12. +8 −7 src/core/observables/Current.hpp
  13. +7 −12 src/core/observables/CylindricalDensityProfile.hpp
  14. +9 −21 src/core/observables/CylindricalFluxDensityProfile.hpp
  15. +11 −22 src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp
  16. +3 −1 src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp
  17. +11 −22 src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp
  18. +4 −1 src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.hpp
  19. +10 −22 src/core/observables/CylindricalVelocityProfile.hpp
  20. +7 −3 src/core/observables/DensityProfile.hpp
  21. +8 −6 src/core/observables/DipoleMoment.hpp
  22. +7 −5 src/core/observables/FluxDensityProfile.hpp
  23. +5 −5 src/core/observables/ForceDensityProfile.hpp
  24. +8 −5 src/core/observables/MagneticDipoleMoment.hpp
  25. +10 −5 src/core/observables/ParticleAngularVelocities.hpp
  26. +8 −5 src/core/observables/ParticleBodyAngularVelocities.hpp
  27. +5 −6 src/core/observables/ParticleBodyVelocities.hpp
  28. +6 −3 src/core/observables/ParticleDistances.hpp
  29. +7 −5 src/core/observables/ParticleForces.hpp
  30. +7 −5 src/core/observables/ParticlePositions.hpp
  31. +7 −5 src/core/observables/ParticleVelocities.hpp
  32. +8 −1 src/core/observables/PidObservable.cpp
  33. +3 −5 src/core/observables/PidObservable.hpp
  34. +36 −48 src/core/polymer.cpp
  35. +2 −23 src/core/polymer.hpp
  36. +1 −1 src/core/statistics.hpp
  37. +7 −27 src/core/statistics_chain.cpp
  38. +35 −28 src/core/statistics_chain.hpp
  39. +3 −6 src/python/espressomd/analyze.pxd
  40. +15 −6 src/python/espressomd/analyze.pyx
  41. +1 −1 testsuite/python/dawaanr-and-bh-gpu.py
  42. +1 −1 testsuite/python/dawaanr-and-dds-gpu.py
  43. +1 −1 testsuite/python/dds-and-bh-gpu.py
@@ -57,7 +57,6 @@
#include "pressure.hpp"
#include "rotation.hpp"
#include "statistics.hpp"
#include "statistics_chain.hpp"
#include "virtual_sites.hpp"

#include "electrostatics_magnetostatics/coulomb.hpp"
@@ -27,19 +27,14 @@
* The corresponding header file is diamond.hpp.
*/

#include <cmath>
#include <cstddef>

#include "PartCfg.hpp"
#include "bonded_interactions/bonded_interaction_data.hpp"
#include "communication.hpp"
#include "constraints.hpp"
#include "constraints/ShapeBasedConstraint.hpp"
#include "diamond.hpp"
#include "global.hpp"
#include "integrate.hpp"
#include "polymer.hpp"
#include "particle_data.hpp"
#include "random.hpp"
#include "statistics.hpp"

#include <utils/Vector.hpp>

#include <cmath>
#include <stdexcept>

int create_counterions(PartCfg &partCfg, int const N_CI, int part_id,
int const mode, double const shield, int const max_try,
@@ -53,7 +48,7 @@ int create_counterions(PartCfg &partCfg, int const N_CI, int part_id,
pos[0] = box_geo.length()[0] * d_random();
pos[1] = box_geo.length()[1] * d_random();
pos[2] = box_geo.length()[2] * d_random();
if ((mode != 0) or (mindist(partCfg, pos) > shield))
if ((mode != 0) or (distto(partCfg, pos) > shield))
break;
}
if (cnt1 >= max_try)
@@ -30,8 +30,6 @@
*/

#include "PartCfg.hpp"
#include "Particle.hpp"
#include "utils/Vector.hpp"

/** C implementation of 'counterions \<N_CI\> [options]'.
* @param N_CI number of counterions to create
@@ -20,7 +20,8 @@
*/
#include "p3m_send_mesh.hpp"

#ifdef P3M
#if defined(P3M) || defined(DP3M)

#include "fft.hpp"

#include <utils/mpi/cart_comm.hpp>
@@ -23,7 +23,7 @@

#include "p3m-common.hpp"

#ifdef P3M
#if defined(P3M) || defined(DP3M)

#include <utils/Span.hpp>

@@ -78,6 +78,5 @@ class p3m_send_mesh {
spread_grid(Utils::make_span(&mesh, 1), comm, dim);
}
};

#endif
#endif // ESPRESSO_P3M_SEND_MESH_HPP
@@ -34,14 +34,15 @@ namespace Observables {
class BondAngles : public PidObservable {
public:
using PidObservable::PidObservable;
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::vector<double> res(n_values());
auto v1 =
get_mi_vector(partCfg[ids()[1]].r.p, partCfg[ids()[0]].r.p, box_geo);
auto v1 = get_mi_vector(particles[1]->r.p, particles[0]->r.p, box_geo);
auto n1 = v1.norm();
for (int i = 0, end = n_values(); i < end; i++) {
auto v2 = get_mi_vector(partCfg[ids()[i + 2]].r.p,
partCfg[ids()[i + 1]].r.p, box_geo);
auto v2 =
get_mi_vector(particles[i + 2]->r.p, particles[i + 1]->r.p, box_geo);
auto n2 = v2.norm();
auto cosine = (v1 * v2) / (n1 * n2);
// sanitize cosine value
@@ -39,16 +39,16 @@ namespace Observables {
class BondDihedrals : public PidObservable {
public:
using PidObservable::PidObservable;
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::vector<double> res(n_values());
auto v1 =
get_mi_vector(partCfg[ids()[1]].r.p, partCfg[ids()[0]].r.p, box_geo);
auto v2 =
get_mi_vector(partCfg[ids()[2]].r.p, partCfg[ids()[1]].r.p, box_geo);
auto v1 = get_mi_vector(particles[1]->r.p, particles[0]->r.p, box_geo);
auto v2 = get_mi_vector(particles[2]->r.p, particles[1]->r.p, box_geo);
auto c1 = vector_product(v1, v2);
for (int i = 0, end = n_values(); i < end; i++) {
auto v3 = get_mi_vector(partCfg[ids()[i + 3]].r.p,
partCfg[ids()[i + 2]].r.p, box_geo);
auto v3 =
get_mi_vector(particles[i + 3]->r.p, particles[i + 2]->r.p, box_geo);
auto c2 = vector_product(v2, v3);
/* the 2-argument arctangent returns an angle in the range [-pi, pi] that
* allows for an unambiguous determination of the 4th particle position */
@@ -21,26 +21,23 @@

#include "PidObservable.hpp"

#include "integrate.hpp"

namespace Observables {

class ComForce : public PidObservable {
public:
using PidObservable::PidObservable;
int n_values() const override { return 3; }
std::vector<double> evaluate(PartCfg &partCfg) const override {
std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::vector<double> res(n_values());
for (int i : ids()) {
if (partCfg[i].p.is_virtual)
for (auto const &p : particles) {
if (p->p.is_virtual)
continue;
res[0] += partCfg[i].f.f[0];
res[1] += partCfg[i].f.f[1];
res[2] += partCfg[i].f.f[2];
res[0] += p->f.f[0];
res[1] += p->f.f[1];
res[2] += p->f.f[2];
}
return res;
};
}
};

} // Namespace Observables
#endif
@@ -24,21 +24,22 @@
#include <vector>

namespace Observables {

class ComPosition : public PidObservable {
public:
using PidObservable::PidObservable;
int n_values() const override { return 3; }
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::vector<double> res(n_values());
double total_mass = 0;
for (int i : ids()) {
if (partCfg[i].p.is_virtual)
for (auto p : particles) {
if (p->p.is_virtual)
continue;
double mass = partCfg[i].p.mass;
res[0] += mass * partCfg[i].r.p[0];
res[1] += mass * partCfg[i].r.p[1];
res[2] += mass * partCfg[i].r.p[2];
double mass = p->p.mass;
res[0] += mass * p->r.p[0];
res[1] += mass * p->r.p[1];
res[2] += mass * p->r.p[2];
total_mass += mass;
}
res[0] /= total_mass;
@@ -24,21 +24,22 @@
#include <vector>

namespace Observables {

class ComVelocity : public PidObservable {
public:
using PidObservable::PidObservable;
int n_values() const override { return 3; }
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::vector<double> res(n_values());
double total_mass = 0;
for (int i : ids()) {
if (partCfg[i].p.is_virtual)
for (auto p : particles) {
if (p->p.is_virtual)
continue;
double mass = partCfg[i].p.mass;
res[0] += mass * partCfg[i].m.v[0];
res[1] += mass * partCfg[i].m.v[1];
res[2] += mass * partCfg[i].m.v[2];
double mass = p->p.mass;
res[0] += mass * p->m.v[0];
res[1] += mass * p->m.v[1];
res[2] += mass * p->m.v[2];
total_mass += mass;
}
res[0] /= total_mass;
@@ -34,14 +34,15 @@ namespace Observables {
class CosPersistenceAngles : public PidObservable {
public:
using PidObservable::PidObservable;
std::vector<double> evaluate(PartCfg &partCfg) const override {
std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
auto const no_of_angles = n_values();
std::vector<double> angles(no_of_angles);
auto const no_of_bonds = n_values() + 1;
std::vector<Utils::Vector3d> bond_vectors(no_of_bonds);
auto get_bond_vector = [&](auto index) {
return get_mi_vector(partCfg[ids()[index + 1]].r.p,
partCfg[ids()[index]].r.p, box_geo);
return get_mi_vector(particles[index + 1]->r.p, particles[index]->r.p,
box_geo);
};
for (int i = 0; i < no_of_bonds; ++i) {
auto const tmp = get_bond_vector(i);
@@ -24,19 +24,20 @@
#include <vector>

namespace Observables {

class Current : public PidObservable {
public:
using PidObservable::PidObservable;
int n_values() const override { return 3; };
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::vector<double> res(n_values());
#ifdef ELECTROSTATICS
for (int i : ids()) {
double charge = partCfg[i].p.q;
res[0] += charge * partCfg[i].m.v[0];
res[1] += charge * partCfg[i].m.v[1];
res[2] += charge * partCfg[i].m.v[2];
for (auto p : particles) {
double charge = p->p.q;
res[0] += charge * p->m.v[0];
res[1] += charge * p->m.v[1];
res[2] += charge * p->m.v[2];
};
#endif
return res;
@@ -27,26 +27,21 @@ namespace Observables {
class CylindricalDensityProfile : public CylindricalPidProfileObservable {
public:
using CylindricalPidProfileObservable::CylindricalPidProfileObservable;
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::array<size_t, 3> n_bins{{static_cast<size_t>(n_r_bins),
static_cast<size_t>(n_phi_bins),
static_cast<size_t>(n_z_bins)}};
std::array<std::pair<double, double>, 3> limits{
{std::make_pair(min_r, max_r), std::make_pair(min_phi, max_phi),
std::make_pair(min_z, max_z)}};
Utils::CylindricalHistogram<double, 3> histogram(n_bins, 1, limits);
std::vector<::Utils::Vector3d> folded_positions;
std::transform(ids().begin(), ids().end(),
std::back_inserter(folded_positions), [&partCfg](int id) {
return ::Utils::Vector3d(
folded_position(partCfg[id].r.p, box_geo));
});
for (auto &p : folded_positions) {
p -= center;
histogram.update(
Utils::transform_coordinate_cartesian_to_cylinder(p, axis));

for (auto p : particles) {
histogram.update(Utils::transform_coordinate_cartesian_to_cylinder(
folded_position(p->r.p, box_geo) - center, axis));
}

histogram.normalize();
return histogram.get_histogram();
}
@@ -27,35 +27,23 @@ namespace Observables {
class CylindricalFluxDensityProfile : public CylindricalPidProfileObservable {
public:
using CylindricalPidProfileObservable::CylindricalPidProfileObservable;
std::vector<double> evaluate(PartCfg &partCfg) const override {

std::vector<double>
evaluate(Utils::Span<const Particle *const> particles) const override {
std::array<size_t, 3> n_bins{{static_cast<size_t>(n_r_bins),
static_cast<size_t>(n_phi_bins),
static_cast<size_t>(n_z_bins)}};
std::array<std::pair<double, double>, 3> limits{
{std::make_pair(min_r, max_r), std::make_pair(min_phi, max_phi),
std::make_pair(min_z, max_z)}};
Utils::CylindricalHistogram<double, 3> histogram(n_bins, 3, limits);
std::vector<::Utils::Vector3d> folded_positions;
std::transform(ids().begin(), ids().end(),
std::back_inserter(folded_positions), [&partCfg](int id) {
return ::Utils::Vector3d(
folded_position(partCfg[id].r.p, box_geo));
});
std::vector<::Utils::Vector3d> velocities;
std::transform(ids().begin(), ids().end(), std::back_inserter(velocities),
[&partCfg](int id) {
return ::Utils::Vector3d{{partCfg[id].m.v[0],
partCfg[id].m.v[1],
partCfg[id].m.v[2]}};
});
for (auto &p : folded_positions)
p -= center;

// Write data to the histogram
for (size_t ind = 0; ind < ids().size(); ++ind) {
histogram.update(Utils::transform_coordinate_cartesian_to_cylinder(
folded_positions[ind], axis),
Utils::transform_vector_cartesian_to_cylinder(
velocities[ind], axis, folded_positions[ind]));
for (auto p : particles) {
auto const pos = folded_position(p->r.p, box_geo) - center;
histogram.update(
Utils::transform_coordinate_cartesian_to_cylinder(pos, axis),
Utils::transform_vector_cartesian_to_cylinder(p->m.v, axis, pos));
}
histogram.normalize();
return histogram.get_histogram();
@@ -21,13 +21,11 @@
#include <utils/Histogram.hpp>
#include <utils/math/coordinate_transformation.hpp>

#include <boost/range/algorithm.hpp>

namespace Observables {

std::vector<double>
CylindricalLBFluxDensityProfileAtParticlePositions::evaluate(
PartCfg &partCfg) const {
Utils::Span<const Particle *const> particles) const {

std::array<size_t, 3> n_bins{{static_cast<size_t>(n_r_bins),
static_cast<size_t>(n_phi_bins),
static_cast<size_t>(n_z_bins)}};
@@ -37,26 +35,17 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate(
Utils::CylindricalHistogram<double, 3> histogram(n_bins, 3, limits);
// First collect all positions (since we want to call the LB function to
// get the fluid velocities only once).
std::vector<Utils::Vector3d> folded_positions(ids().size());
boost::transform(ids(), folded_positions.begin(),
[&partCfg](int id) -> Utils::Vector3d {
return folded_position(partCfg[id].r.p, box_geo);
});

std::vector<Utils::Vector3d> velocities(folded_positions.size());
boost::transform(folded_positions, velocities.begin(),
[](const Utils::Vector3d &pos) {
return lb_lbfluid_get_interpolated_velocity(pos) *
lb_lbfluid_get_lattice_speed();
});
for (auto &p : folded_positions)
p -= center;
for (int ind = 0; ind < ids().size(); ++ind) {
histogram.update(Utils::transform_coordinate_cartesian_to_cylinder(
folded_positions[ind], axis),
Utils::transform_vector_cartesian_to_cylinder(
velocities[ind], axis, folded_positions[ind]));
for (auto p : particles) {
auto const pos = folded_position(p->r.p, box_geo);
auto const v = lb_lbfluid_get_interpolated_velocity(pos) *
lb_lbfluid_get_lattice_speed();

histogram.update(
Utils::transform_coordinate_cartesian_to_cylinder(pos - center, axis),
Utils::transform_vector_cartesian_to_cylinder(v, axis, pos - center));
}

histogram.normalize();
return histogram.get_histogram();
}

0 comments on commit a6b4657

Please sign in to comment.
You can’t perform that action at this time.