diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index a4214956f67..9ec1d9daf46 100644 --- a/doc/sphinx/analysis.rst +++ b/doc/sphinx/analysis.rst @@ -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 diff --git a/doc/sphinx/figures/analysis_cylindrical_average.png b/doc/sphinx/figures/analysis_cylindrical_average.png deleted file mode 100644 index 9b792a0d220..00000000000 Binary files a/doc/sphinx/figures/analysis_cylindrical_average.png and /dev/null differ diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index cd44c59ab0f..1388de3a72d 100644 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -506,117 +506,6 @@ std::vector> modify_stucturefactor(int order, return structure_factor; } -int calc_cylindrical_average( - PartCfg &partCfg, std::vector const ¢er_, - std::vector const &direction_, double length, double radius, - int bins_axial, int bins_radial, std::vector types, - std::map>>> - &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>>>( - "density", - std::vector>>(types.size()))); - distribution.insert( - std::pair>>>( - "v_r", std::vector>>(types.size()))); - distribution.insert( - std::pair>>>( - "v_t", std::vector>>(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(floor(height / binwd_radial)); - index_axial = - static_cast(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 ****************************************************************************************/ diff --git a/src/core/statistics.hpp b/src/core/statistics.hpp index bacca55913d..71e918fa307 100644 --- a/src/core/statistics.hpp +++ b/src/core/statistics.hpp @@ -185,13 +185,6 @@ std::vector calc_structurefactor(PartCfg &, int const *p_types, std::vector> modify_stucturefactor(int order, double const *sf); -int calc_cylindrical_average( - PartCfg &, std::vector const ¢er, - std::vector const &direction, double length, double radius, - int bins_axial, int bins_radial, std::vector types, - std::map>>> - &distribution); - template double min_distance2(Utils::Vector const &pos1, Utils::Vector const &pos2) { diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index 7760d3b8cb8..d96dd697761 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -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) diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 4774c48cb0b..4d32021b4db 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -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