Segmentation fault due to division by zero in helper funciton `coords_to_index` in neighbor module

[ipcamit]

Josh stumbled upon the bug where for a system where maximum coordinate == minimum coordinate != 0, then there can be cases where following snippet will result in division by zero:

File kliff/neighbor/helper.hpp:96

inline void coords_to_index(double const * x,
                            int const * size,
                            double const * max,
                            double const * min,
                            int * const index)
{
  index[0] = static_cast<int>(((x[0] - min[0]) / (max[0] - min[0])) * size[0]);
  index[0] = std::min(index[0], size[0] - 1);  // edge case when x[0] = max[0]
  index[1] = static_cast<int>(((x[1] - min[1]) / (max[1] - min[1])) * size[1]);
  index[1] = std::min(index[1], size[1] - 1);  // edge case when x[1] = max[1]
  index[2] = static_cast<int>(((x[2] - min[2]) / (max[2] - min[2])) * size[2]);
  index[2] = std::min(index[2], size[2] - 1);  // edge case when x[2] = max[2]
}

As then (max + eps) - min == 0.

Reproducible minimal example:

import numpy as np

cutoff = 5.0
r0 = 2.0
vac = 100
# Build the atoms object
# Check KLIFF results
from kliff.neighbor import NeighborList as KLIFFNeighborList
from kliff.dataset import Configuration as KLIFFConfig

kliff_atoms = KLIFFConfig(
    cell=np.diag(np.diag(np.array([202.0, 202.0, 202.0]))),
    species=["H", "H"],
    coords=np.array([[100.0, 2.0, 0.0], [101.0, 2.0, 0.0]]),
    PBC=np.array([False, False, False]),
)
nl_kliff = KLIFFNeighborList(kliff_atoms, infl_dist=cutoff)
nl_kliff.get_neigh(0)

It can be mitigated simply by checking for zero

inline void coords_to_index(double const * x,
                            int const * size,
                            double const * max,
                            double const * min,
                            int * const index)
{
  index[0] = static_cast<int>(((x[0] - min[0]) / std::max(SMALL,max[0] - min[0])) * size[0]);
  index[0] = std::min(index[0], size[0] - 1);  // edge case when x[0] = max[0]
  index[1] = static_cast<int>(((x[1] - min[1]) / std::max(SMALL,max[1] - min[1])) * size[1]);
  index[1] = std::min(index[1], size[1] - 1);  // edge case when x[1] = max[1]
  index[2] = static_cast<int>(((x[2] - min[2]) / std::max(SMALL,max[2] - min[2])) * size[2]);
  index[2] = std::min(index[2], size[2] - 1);  // edge case when x[2] = max[2]
}

Or more accurately

#include <limits>
...
inline void coords_to_index(double const * x,
                            int const * size,
                            double const * max,
                            double const * min,
                            int * const index)
{
  index[0] = static_cast<int>(((x[0] - min[0]) / std::max(std::numeric_limits<double>::epsilon(),max[0] - min[0])) * size[0]);
  index[0] = std::min(index[0], size[0] - 1);  // edge case when x[0] = max[0]
  index[1] = static_cast<int>(((x[1] - min[1]) / std::max(std::numeric_limits<double>::epsilon(),max[1] - min[1])) * size[1]);
  index[1] = std::min(index[1], size[1] - 1);  // edge case when x[1] = max[1]
  index[2] = static_cast<int>(((x[2] - min[2]) / std::max(std::numeric_limits<double>::epsilon(),max[2] - min[2])) * size[2]);
  index[2] = std::min(index[2], size[2] - 1);  // edge case when x[2] = max[2]
}

Following above changes, above example works correctly.

[mjwen]

Yeah, this has been fixed in the latest release v0.3.3.

The problem in the previous version is that std::numeric_limits<double>::epsilon() is used as eps to be added to max. Then when we do (max+eps) - min, it got evaluated to 0, instead of std::numeric_limits<double>::epsilon() due to numerical error.

We fix it is by using a larger eps (1e-10 now). But I think something like index[0] = static_cast<int>(((x[0] - min[0]) / std::max(std::numeric_limits<double>::epsilon(),max[0] - min[0])) * size[0]); should be useful to avoid similar problems. And we should include it.

Do you want to make a PR for this?

[ipcamit]

Sure. Will submit the pr tomorrow morning.

On Thu, Mar 31, 2022, 00:07 Mingjian Wen @.***> wrote:

Yeah, this has been fixed in the latest release v0.3.3.

The problem in the previous version is that std::numeric_limits::epsilon() is used as eps to be added to max. Then when we do (max+eps) - min, it got evaluated to 0, instead of std::numeric_limits::epsilon() due to numerical error.

The way we fixed it is by using a larger eps (1e-10 now). But I think some thing like index[0] = static_cast(((x[0] - min[0]) / std::max(std::numeric_limits::epsilon(),max[0] - min[0])) * size[0]); should be useful to avoid similar problems. And we should include it.

Do you want to make a PR for this?

— Reply to this email directly, view it on GitHub https://github.com/openkim/kliff/issues/43#issuecomment-1084080731, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABTZ2D4JUTMP5RQB4CRFLVTVCUXJ5ANCNFSM5SEDGZOA . You are receiving this because you authored the thread.Message ID: @.***>