Skip to content

Commit

Permalink
Fix bugs in espressomd.analyze module (espressomd#4534)
Browse files Browse the repository at this point in the history 
Description of changes:
- fix particle double-counting bug in the structure factor analysis code
- fix distance calculation bug in the particle distribution analysis code when `r_max > np.sum(system.box_l)`
- remove unused v_kappa analysis code
  • Loading branch information
@kodiakhq @jngrad
kodiakhq[bot] authored and jngrad committed Nov 25, 2022
1 parent 4870af6 commit e5d7313
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 103 deletions.
73 changes: 33 additions & 40 deletions src/core/statistics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,8 +184,9 @@ void calc_part_distribution(PartCfg &partCfg, std::vector<int> const &p1_types,
double inv_bin_width = 0.0;
double min_dist, min_dist2 = 0.0, start_dist2;

start_dist2 = Utils::sqr(box_geo.length()[0] + box_geo.length()[1] +
box_geo.length()[2]);
auto const r_max2 = Utils::sqr(r_max);
auto const r_min2 = Utils::sqr(r_min);
start_dist2 = Utils::sqr(r_max + 1.);
/* bin preparation */
*low = 0.0;
for (int i = 0; i < r_bins; i++)
Expand All @@ -197,40 +198,36 @@ void calc_part_distribution(PartCfg &partCfg, std::vector<int> const &p1_types,

/* particle loop: p1_types */
for (auto const &p1 : partCfg) {
for (int t1 : p1_types) {
if (p1.type() == t1) {
min_dist2 = start_dist2;
/* particle loop: p2_types */
for (auto const &p2 : partCfg) {
if (p1 != p2) {
for (int t2 : p2_types) {
if (p2.type() == t2) {
auto const act_dist2 =
box_geo.get_mi_vector(p1.pos(), p2.pos()).norm2();
if (act_dist2 < min_dist2) {
min_dist2 = act_dist2;
}
}
if (Utils::contains(p1_types, p1.type())) {
min_dist2 = start_dist2;
/* particle loop: p2_types */
for (auto const &p2 : partCfg) {
if (p1 != p2) {
if (Utils::contains(p2_types, p2.type())) {
auto const act_dist2 =
box_geo.get_mi_vector(p1.pos(), p2.pos()).norm2();
if (act_dist2 < min_dist2) {
min_dist2 = act_dist2;
}
}
}
min_dist = sqrt(min_dist2);
if (min_dist <= r_max) {
if (min_dist >= r_min) {
/* calculate bin index */
if (log_flag)
ind = (int)((log(min_dist / r_min)) * inv_bin_width);
else
ind = (int)((min_dist - r_min) * inv_bin_width);
if (ind >= 0 && ind < r_bins) {
dist[ind] += 1.0;
}
} else {
*low += 1.0;
}
if (min_dist2 <= r_max2) {
if (min_dist2 >= r_min2) {
min_dist = sqrt(min_dist2);
/* calculate bin index */
if (log_flag)
ind = (int)((log(min_dist / r_min)) * inv_bin_width);
else
ind = (int)((min_dist - r_min) * inv_bin_width);
if (ind >= 0 && ind < r_bins) {
dist[ind] += 1.0;
}
} else {
*low += 1.0;
}
cnt++;
}
cnt++;
}
}
if (cnt == 0)
Expand Down Expand Up @@ -260,13 +257,10 @@ void calc_structurefactor(PartCfg &partCfg, std::vector<int> const &p_types,
if ((static_cast<std::size_t>(n) <= order_sq) && (n >= 1)) {
double C_sum = 0.0, S_sum = 0.0;
for (auto const &p : partCfg) {
for (int t : p_types) {
if (p.type() == t) {
auto const qr =
twoPI_L * (Utils::Vector3i{{i, j, k}} * p.pos());
C_sum += cos(qr);
S_sum += sin(qr);
}
if (Utils::contains(p_types, p.type())) {
auto const qr = twoPI_L * (Utils::Vector3i{{i, j, k}} * p.pos());
C_sum += cos(qr);
S_sum += sin(qr);
}
}
ff[2 * n - 2] += C_sum * C_sum + S_sum * S_sum;
Expand All @@ -278,9 +272,8 @@ void calc_structurefactor(PartCfg &partCfg, std::vector<int> const &p_types,

long n_particles = 0l;
for (auto const &p : partCfg) {
for (int t : p_types) {
if (p.type() == t)
n_particles++;
if (Utils::contains(p_types, p.type())) {
n_particles++;
}
}

Expand Down
68 changes: 5 additions & 63 deletions src/python/espressomd/analyze.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -541,8 +541,8 @@ class Analysis:
type B, disregarding the fact that a spherical shell of a larger radius
covers a larger volume). The distance is defined as the minimal distance
between a particle of group ``type_list_a`` to any of the group
``type_list_b``. Returns two arrays, the bins and the (normalized)
distribution.
``type_list_b`` (excluding the self-contribution). Returns two arrays,
the bins and the (normalized) distribution.
Parameters
----------
Expand All @@ -555,7 +555,9 @@ class Analysis:
r_min : :obj:`float`
Minimum distance.
r_max : :obj:`float`
Maximum distance.
Maximum distance. By default, it is half the box size.
A value larger than half the box size is allowed for systems
with :ref:`open boundary conditions <Open boundaries>`.
r_bins : :obj:`int`
Number of bins.
log_flag : :obj:`bool`
Expand Down Expand Up @@ -743,63 +745,3 @@ class Analysis:
MofImatrix_np[i] = MofImatrix[i]

return MofImatrix_np.reshape((3, 3))

#
# Vkappa
#

_Vkappa = {
"Vk1": 0.0,
"Vk2": 0.0,
"avk": 0.0
}

def v_kappa(self, mode=None, Vk1=None, Vk2=None, avk=None):
"""
.. todo:: Looks to be incomplete
Calculate the compressibility through volume fluctuations.
Parameters
----------
mode : :obj:`str`, \{'read', 'set' or 'reset'\}
Mode.
Vk1 : :obj:`float`
Volume.
Vk2 : :obj:`float`
Volume squared.
avk : :obj:`float`
Number of averages.
"""

utils.check_type_or_throw_except(
mode, 1, str, "mode has to be a string")

if mode == "reset":
self._Vkappa["Vk1"] = 0.0
self._Vkappa["Vk2"] = 0.0
self._Vkappa["avk"] = 0.0
elif mode == "read":
return self._Vkappa
elif mode == "set":
utils.check_type_or_throw_except(
Vk1, 1, float, "Vk1 has to be a float")
self._Vkappa["Vk1"] = Vk1
utils.check_type_or_throw_except(
Vk2, 1, float, "Vk2 has to be a float")
self._Vkappa["Vk2"] = Vk2
utils.check_type_or_throw_except(
avk, 1, float, "avk has to be a float")
self._Vkappa["avk"] = avk
if self._Vkappa["avk"] <= 0.0:
result = self._Vkappa["Vk1"] = self._Vkappa[
"Vk2"] = self._Vkappa["avk"] = 0.0
raise ValueError(
"# of averages <avk> has to be positive! Resetting values.")
else:
result = self._Vkappa["Vk2"] / self._Vkappa["avk"] - \
(self._Vkappa["Vk1"] / self._Vkappa["avk"])**2
return result
else:
raise ValueError(f"Unknown mode {mode!r}")

0 comments on commit e5d7313

Please sign in to comment.