Skip to content

Commit

Permalink
Remove cylindrical_average. (#3470)
Browse files Browse the repository at this point in the history
Fixes #3450
  • Loading branch information
kodiakhq[bot] authored Feb 11, 2020
2 parents e0a0e7c + 9e25bb0 commit 4d1512c
Showing 6 changed files with 0 additions and 293 deletions.
71 changes: 0 additions & 71 deletions doc/sphinx/analysis.rst
Original file line number Diff line number Diff line change
@@ -28,7 +28,6 @@ The direct analysis commands can be classified into two types:
- :ref:`Minimal distances between particles`
- :ref:`Particles in the neighborhood`
- :ref:`Particle distribution`
- :ref:`Cylindrical Average`
- :ref:`Radial distribution function` with ``rdf_type='rdf'``
- :ref:`Structure factor`
- :ref:`Center of mass`
@@ -186,76 +185,6 @@ Two arrays are returned corresponding to the normalized distribution and the bin
central axis, which one would think is natural.


.. _Cylindrical average:

Cylindrical Average
~~~~~~~~~~~~~~~~~~~

:meth:`espressomd.analyze.Analysis.cylindrical_average`

Calculates the particle distribution using cylindrical binning.

The volume considered is inside a cylinder defined by the parameters ``center``, ``axis``, ``length`` and ``radius``.

The geometrical details of the cylindrical binning is defined using ``bins_axial`` and ``bins_radial`` which are the number bins in the axial and radial directions (respectively).
See figure :ref:`cylindrical_average` for a visual representation of the binning geometry.

.. _cylindrical_average:

.. figure:: figures/analysis_cylindrical_average.png
:alt: Geometry for the cylindrical binning
:align: center
:height: 6.00000cm

Geometry for the cylindrical binning


The command returns a list of lists. The outer list contains all data
combined whereas each inner list contains one line. Each lines stores a
different combination of the radial and axial index. The output might
look something like this

.. code-block:: numpy
[ [ 0 0 0.05 -0.25 0.0314159 0 0 0 0 0 0 ]
[ 0 1 0.05 0.25 0.0314159 31.831 1.41421 1 0 0 0 ]
... ]
In this case two different particle types were present.
The columns of the respective lines are coded like this

============= ============ =========== ========== ========= ======= ======== ======== ======= ========= =======
index_radial index_axial pos_radial pos_axial binvolume density v_radial v_axial density v_radial v_axial
============= ============ =========== ========== ========= ======= ======== ======== ======= ========= =======
0 0 0.05 -0.25 0.0314159 0 0 0 0 0 0
0 1 0.05 0.25 0.0314159 31.831 1.41421 1 0 0 0
============= ============ =========== ========== ========= ======= ======== ======== ======= ========= =======

As one can see the columns **density**, **v_radial** and **v_axial** appear twice.
The order of appearance corresponds to the order of the types in the argument ``types``.
For example if was set to ``types=[0, 1]`` then the first triple is associated to type 0 and
the second triple to type 1.

..
.. _Vkappa:
Vkappa
~~~~~~
:meth:`espressomd.analyze.Analysis.v_kappa`

.. todo:: Implementation appears to be incomplete

Calculates the compressibility :math:`V \times \kappa_T` through the
Volume fluctuations
:math:`V \times \kappa_T = \beta \left(\langle V^2\rangle - \langle V \rangle^2\right)`
:cite:`kolb99a`. Given no arguments this function calculates
and returns the current value of the running average for the volume
fluctuations. The ``mode=reset`` argument clears the currently stored values. With ``mode=read`` the
cumulative mean volume, cumulative mean squared volume and how many
samples were used can be retrieved. Likewise the option ``mode=set`` enables you to
set those.


.. _Radial distribution function:

Radial distribution function
Binary file removed doc/sphinx/figures/analysis_cylindrical_average.png
Binary file not shown.
111 changes: 0 additions & 111 deletions src/core/statistics.cpp
Original file line number Diff line number Diff line change
@@ -506,117 +506,6 @@ std::vector<std::vector<double>> modify_stucturefactor(int order,
return structure_factor;
}

int calc_cylindrical_average(
PartCfg &partCfg, std::vector<double> const &center_,
std::vector<double> const &direction_, double length, double radius,
int bins_axial, int bins_radial, std::vector<int> types,
std::map<std::string, std::vector<std::vector<std::vector<double>>>>
&distribution) {
int index_axial;
int index_radial;
double binwd_axial = length / bins_axial;
double binwd_radial = radius / bins_radial;

auto center = Utils::Vector3d{center_};
auto direction = Utils::Vector3d{direction_};

// Select all particle types if the only entry in types is -1
bool all_types = false;
if (types.size() == 1 && types[0] == -1)
all_types = true;

distribution.insert(
std::pair<std::string, std::vector<std::vector<std::vector<double>>>>(
"density",
std::vector<std::vector<std::vector<double>>>(types.size())));
distribution.insert(
std::pair<std::string, std::vector<std::vector<std::vector<double>>>>(
"v_r", std::vector<std::vector<std::vector<double>>>(types.size())));
distribution.insert(
std::pair<std::string, std::vector<std::vector<std::vector<double>>>>(
"v_t", std::vector<std::vector<std::vector<double>>>(types.size())));

for (unsigned int type = 0; type < types.size(); type++) {
distribution["density"][type].resize(bins_radial);
distribution["v_r"][type].resize(bins_radial);
distribution["v_t"][type].resize(bins_radial);
for (int index_radial = 0; index_radial < bins_radial; index_radial++) {
distribution["density"][type][index_radial].assign(bins_axial, 0.0);
distribution["v_r"][type][index_radial].assign(bins_axial, 0.0);
distribution["v_t"][type][index_radial].assign(bins_axial, 0.0);
}
}

auto const norm_direction = direction.norm();

for (auto const &p : partCfg) {
for (unsigned int type_id = 0; type_id < types.size(); type_id++) {
if (types[type_id] == p.p.type || all_types) {
auto const pos = folded_position(p.r.p, box_geo);

Utils::Vector3d vel{p.m.v};

auto const diff = pos - center;

// Find the height of the particle above the axis (height) and
// the distance from the center point (dist)
auto const hat = vector_product(direction, diff);
auto const height = hat.norm();
auto const dist = direction * diff / norm_direction;

// Determine the components of the velocity parallel and
// perpendicular to the direction vector
double v_radial;
if (height == 0)
v_radial = vector_product(vel, direction).norm() / norm_direction;
else
v_radial = vel * hat / height;

auto const v_axial = vel * direction / norm_direction;

// Work out relevant indices for x and y
index_radial = static_cast<int>(floor(height / binwd_radial));
index_axial =
static_cast<int>(floor((dist + 0.5 * length) / binwd_axial));

if ((index_radial < bins_radial && index_radial >= 0) &&
(index_axial < bins_axial && index_axial >= 0)) {
distribution["density"][type_id][index_radial][index_axial] += 1;
distribution["v_r"][type_id][index_radial][index_axial] += v_radial;
distribution["v_t"][type_id][index_radial][index_axial] += v_axial;
}
}
}
}

// Now we turn the counts into densities by dividing by one radial
// bin (binvolume). We also divide the velocities by the counts.
double binvolume;
for (unsigned int type_id = 0; type_id < types.size(); type_id++) {
for (int index_radial = 0; index_radial < bins_radial; index_radial++) {
// All bins are cylindrical shells of thickness binwd_radial.
// The volume is thus: binvolume = pi*(r_outer - r_inner)^2 * length
if (index_radial == 0)
binvolume = M_PI * binwd_radial * binwd_radial * length;
else
binvolume = M_PI * (index_radial * index_radial + 2 * index_radial) *
binwd_radial * binwd_radial * length;
for (int index_axial = 0; index_axial < bins_axial; index_axial++) {
if (distribution["density"][type_id][index_radial][index_axial] != 0) {
distribution["v_r"][type_id][index_radial][index_axial] /=
distribution["density"][type_id][index_radial][index_axial];
distribution["v_t"][type_id][index_radial][index_axial] /=
distribution["density"][type_id][index_radial][index_axial];
distribution["density"][type_id][index_radial][index_axial] /=
binvolume;
}
}
}
}

return ES_OK;
}

/****************************************************************************************
* config storage functions
****************************************************************************************/
7 changes: 0 additions & 7 deletions src/core/statistics.hpp
Original file line number Diff line number Diff line change
@@ -185,13 +185,6 @@ std::vector<double> calc_structurefactor(PartCfg &, int const *p_types,
std::vector<std::vector<double>> modify_stucturefactor(int order,
double const *sf);

int calc_cylindrical_average(
PartCfg &, std::vector<double> const &center,
std::vector<double> const &direction, double length, double radius,
int bins_axial, int bins_radial, std::vector<int> types,
std::map<std::string, std::vector<std::vector<std::vector<double>>>>
&distribution);

template <typename T>
double min_distance2(Utils::Vector<T, 3> const &pos1,
Utils::Vector<T, 3> const &pos2) {
4 changes: 0 additions & 4 deletions src/python/espressomd/analyze.pxd
Original file line number Diff line number Diff line change
@@ -79,10 +79,6 @@ cdef extern from "statistics.hpp":
cdef double * obsstat_nonbonded_intra(Observable_stat_non_bonded * stat, int i, int j)
cdef vector[double] calc_linear_momentum(int include_particles, int include_lbfluid)
cdef vector[double] centerofmass(PartCfg &, int part_type)
cdef int calc_cylindrical_average(
PartCfg &, vector[double] center, vector[double] direction,
double length, double radius, int bins_axial, int bins_radial,
vector[int] types, map[string, vector[vector[vector[double]]]] & distribution)

void calc_rdf(PartCfg &, vector[int] p1_types, vector[int] p2_types,
double r_min, double r_max, int r_bins, vector[double] rdf)
100 changes: 0 additions & 100 deletions src/python/espressomd/analyze.pyx
Original file line number Diff line number Diff line change
@@ -274,106 +274,6 @@ class Analysis:

return create_nparray_from_int_list(ids)

def cylindrical_average(self, center=None, axis=None,
length=None, radius=None,
bins_axial=None, bins_radial=None,
types=[-1]):
"""
Calculate the particle distribution using cylindrical binning.
Parameters
----------
center : (3,) array_like of :obj:`float`
Coordinates of the centre of the cylinder.
axis : (3,) array_like of :obj:`float`
Axis vectory of the cylinder, does not need to be normalized.
length : :obj:`float`
Length of the cylinder.
radius : :obj:`float`
Radius of the cylinder.
bins_axial : :obj:`int`
Number of axial bins.
bins_radial : :obj:`int`
Number of radial bins.
types : lists of :obj:`int`
Particle :attr:`~espressomd.particle_data.ParticleHandle.type`
Returns
-------
list of lists
columns indicate ``index_radial``, ``index_axial``, ``pos_radial``,
``pos_axial``, ``binvolume``, ``density``, ``v_radial``,
``v_axial``, ``density``, ``v_radial`` and ``v_axial``.
Note that the columns ``density``, ``v_radial`` and ``v_axial``
appear for each type indicated in ``types``, in the same order.
"""

# Check the input types
check_type_or_throw_except(
center, 3, float, "center has to be 3 floats")
check_type_or_throw_except(
axis, 3, float, "direction has to be 3 floats")
check_type_or_throw_except(
length, 1, float, "length has to be a float")
check_type_or_throw_except(
radius, 1, float, "radius has to be a float")
check_type_or_throw_except(
bins_axial, 1, int, "bins_axial has to be an int")
check_type_or_throw_except(
bins_radial, 1, int, "bins_radial has to be an int")

# Convert Python types to C++ types
cdef vector[double] c_center = center
cdef vector[double] c_direction = axis
cdef double c_length = length
cdef double c_radius = radius
cdef int c_bins_axial = bins_axial
cdef int c_bins_radial = bins_radial
cdef vector[int] c_types = types

cdef map[string, vector[vector[vector[double]]]] distribution
analyze.calc_cylindrical_average(
analyze.partCfg(), c_center, c_direction, c_length,
c_radius, c_bins_axial, c_bins_radial, c_types, distribution)

cdef double binwd_axial = c_length / c_bins_axial
cdef double binwd_radial = c_radius / c_bins_radial
cdef double binvolume, pos_radial, pos_axial

cdef vector[string] names = [b"density", b"v_r", b"v_t"]

buffer = np.empty(
[bins_radial * bins_axial, 5 + c_types.size() * names.size()])

cdef int index_radial, index_axial
cdef unsigned int type_id, name
for index_radial in range(0, bins_radial):
for index_axial in range(0, bins_axial):
pos_radial = (index_radial + .5) * binwd_radial
pos_axial = (index_axial + .5) * binwd_axial - .5 * c_length

if index_radial == 0:
binvolume = np.pi * binwd_radial * binwd_radial * c_length
else:
binvolume = np.pi * \
(index_radial * index_radial + 2 * index_radial) * \
binwd_radial * binwd_radial * c_length

buffer[index_axial + bins_axial *
index_radial, 0] = index_radial
buffer[index_axial + bins_axial *
index_radial, 1] = index_axial
buffer[index_axial + bins_axial * index_radial, 2] = pos_radial
buffer[index_axial + bins_axial * index_radial, 3] = pos_axial
buffer[index_axial + bins_axial * index_radial, 4] = binvolume
for type_id in range(0, c_types.size()):
for name in range(0, names.size()):
buffer[index_axial + bins_axial * index_radial, 5 + name + type_id * names.size()] \
= distribution[names[name]][type_id][index_radial][index_axial]

return buffer

def pressure(self, v_comp=False):
"""Calculate the instantaneous pressure (in parallel). This is only
sensible in an isotropic system which is homogeneous (on average)! Do

0 comments on commit 4d1512c

Please sign in to comment.