eonc::IRACompare Class Reference

C++ wrapper for the IRA (Iterative Rotations and Assignments) library. More...

#include <IRACompare.h>

Classes
struct	MatchResult
struct	SymmetryResult

Static Public Member Functions
static MatchResult	match (const Matter &m1, const Matter &m2, double distThreshold)
	Match two structures using CShDA + SVD (optimal rotation + assignment).
static MatchResult	matchArrays (int nat1, const int *typ1, const double *pos1, int nat2, const int *typ2, const double *pos2, double distThreshold)
	Same as match(), from packed (n,3) row-major coordinates and Z arrays.
static MatchResult	matchArrays (int nat1, const int *typ1, const double *pos1, int nat2, const int *typ2, const double *pos2, double distThreshold, IIRAResource &res)
	Test seam: injected resource, no process-default libira load.
static MatchResult	matchPBC (const Matter &m1, const Matter &m2, double distThreshold)
	Atom assignment under periodic boundary conditions (CShDA only, no rotation/SVD).
static MatchResult	matchPBC (const Matter &m1, const Matter &m2, double distThreshold, IIRAResource &res)
static SymmetryResult	findSymmetry (const Matter &m, double threshold, bool prescreenIh=true)
	Find all symmetry operations of a structure (SOFI algorithm).
static SymmetryResult	findSymmetry (const Matter &m, double threshold, bool prescreenIh, IIRAResource &res)
static MatchResult	alignReactantToProduct (Matter &reactant, const Matter &product, double distThreshold)
	Rigid-align + permute reactant onto product.

Detailed Description

C++ wrapper for the IRA (Iterative Rotations and Assignments) library.

Provides structure comparison via CShDA and symmetry detection via SOFI.

Definition at line 25 of file IRACompare.h.

Member Function Documentation

alignReactantToProduct()

IRACompare::MatchResult eonc::IRACompare::alignReactantToProduct ( Matter & reactant, const Matter & product, double distThreshold )

static

Rigid-align + permute reactant onto product.

Returns the match (error != 0 means reactant was left unchanged).

Definition at line 339 of file IRACompare.cpp.

                                                         {
  MatchResult result = match(reactant, product, distThreshold);
  if (result.error != 0) {
    return result;
  }
  const long n = reactant.numberOfAtoms();
  if (n != product.numberOfAtoms() ||
      static_cast<long>(result.permutation.size()) != n) {
    result.error = -2;
    return result;
  }
  AtomMatrix pos = reactant.getPositionsCopy();
  for (long i = 0; i < n; ++i) {
    Eigen::Vector3d p = pos.row(i);
    pos.row(i) = (result.rotation * p + result.translation).transpose();
  }
  bool oneBased = false;
  for (int idx : result.permutation) {
    if (idx == n) {
      oneBased = true;
      break;
    }
  }
  AtomMatrix reordered = pos;
  for (long i = 0; i < n; ++i) {
    int src = result.permutation[static_cast<size_t>(i)];
    if (oneBased) {
      --src;
    }
    if (src < 0 || src >= n) {
      result.error = -3;
      return result;
    }
    reordered.row(i) = pos.row(src);
  }
  reactant.setPositions(reordered);
  return result;
}

findSymmetry() [1/2]

IRACompare::SymmetryResult eonc::IRACompare::findSymmetry ( const Matter & m, double threshold, bool prescreenIh, IIRAResource & res )

static

Definition at line 229 of file IRACompare.cpp.

                                                                       {
  SymmetryResult result;
  std::lock_guard<std::mutex> lock(res.library_mutex);
 
  try {
    res.require_loaded();
  } catch (const std::exception &e) {
    result.error = -1;
    return result;
  }
 
  const int nat = m.numberOfAtoms();
  std::vector<int> typ(nat);
  auto nrs = m.getAtomicNrs();
  for (int i = 0; i < nat; i++)
    typ[i] = nrs[i];
 
  // Prepare coordinates using direct Eigen data access
  const AtomMatrix &pos = m.getPositions();
  Eigen::Map<const AtomMatrixF> coords_map(pos.data(), 3, nat);
 
  // libira_compute_all requires pre-allocated output arrays of size nmax
  int nmax = res.get_get_nmax_fn()();
 
  int n_mat = 0;
  std::vector<double> mat_buf(9 * nmax);
  std::vector<int> perm_buf(nat * nmax);
  std::vector<char> op_buf(nmax + 1);
  std::vector<int> n_buf(nmax);
  std::vector<int> p_buf(nmax);
  std::vector<double> ax_buf(3 * nmax);
  std::vector<double> angle_buf(nmax);
  std::vector<double> dH_buf(nmax);
  std::vector<char> pg_buf(11);
  int n_prin_ax = 0;
  std::vector<double> prin_ax_buf(3 * nmax);
  int cerr = 0;
 
  double *mat_data = mat_buf.data();
  int *perm_data = perm_buf.data();
  char *op_data = op_buf.data();
  int *n_data = n_buf.data();
  int *p_data = p_buf.data();
  double *ax_data = ax_buf.data();
  double *angle_data = angle_buf.data();
  double *dH_data = dH_buf.data();
  char *pg = pg_buf.data();
  double *prin_ax = prin_ax_buf.data();
 
  res.get_compute_all_fn()(nat, typ.data(), coords_map.data(), threshold,
                           prescreenIh ? 1 : 0, &n_mat, &mat_data, &perm_data,
                           &op_data, &n_data, &p_data, &ax_data, &angle_data,
                           &dH_data, &pg, &n_prin_ax, &prin_ax, &cerr);
 
  result.error = cerr;
  if (cerr == 0) {
    result.nOperations = n_mat;
    if (pg) {
      result.pointGroup = std::string(pg);
    }
 
    result.operations.resize(n_mat);
    for (int k = 0; k < n_mat; k++) {
      Eigen::Map<const Matrix3d> op_map(&mat_data[k * 9]);
      result.operations[k] =
          op_map.transpose(); // Transpose to get proper row-major layout
    }
 
    if (angle_data) {
      result.angles.assign(angle_data, angle_data + n_mat);
    }
    if (ax_data) {
      result.axes.resize(n_mat);
      for (int k = 0; k < n_mat; k++) {
        result.axes[k] = Eigen::Map<const Vector3d>(&ax_data[k * 3]);
      }
    }
  }
 
  // Only free if Fortran reallocated (pointer changed from our pre-allocated
  // buffer). If pointer is unchanged, the vector destructor handles cleanup.
  // Reallocated buffers come from libira's c_malloc side, so std::free is the
  // matching deallocator; never free the vector-backed originals.
  if (mat_data != mat_buf.data())
    std::free(mat_data);
  if (perm_data != perm_buf.data())
    std::free(perm_data);
  if (op_data != op_buf.data())
    std::free(op_data);
  if (n_data != n_buf.data())
    std::free(n_data);
  if (p_data != p_buf.data())
    std::free(p_data);
  if (ax_data != ax_buf.data())
    std::free(ax_data);
  if (angle_data != angle_buf.data())
    std::free(angle_data);
  if (dH_data != dH_buf.data())
    std::free(dH_data);
  if (pg != pg_buf.data())
    std::free(pg);
  if (prin_ax != prin_ax_buf.data())
    std::free(prin_ax);
  return result;
}

findSymmetry() [2/2]

IRACompare::SymmetryResult eonc::IRACompare::findSymmetry ( const Matter & m, double threshold, bool prescreenIh = true )

static

Find all symmetry operations of a structure (SOFI algorithm).

Definition at line 219 of file IRACompare.cpp.

                                                                            {
#ifdef WITH_IRA
  return findSymmetry(m, threshold, prescreenIh, get_ira_resource());
#else
  SymmetryResult result;
  result.error = -1;
  return result;
#endif
}

match()

IRACompare::MatchResult eonc::IRACompare::match ( const Matter & m1, const Matter & m2, double distThreshold )

static

Match two structures using CShDA + SVD (optimal rotation + assignment).

Both structures must have the same atom types.

Definition at line 26 of file IRACompare.cpp.

                                                                {
  const int nat1 = m1.numberOfAtoms();
  const int nat2 = m2.numberOfAtoms();
  if (nat1 <= 0 || nat2 <= 0) {
    MatchResult result;
    result.error = -1;
    return result;
  }
  std::vector<int> typ1(static_cast<size_t>(nat1)),
      typ2(static_cast<size_t>(nat2));
  auto nrs1 = m1.getAtomicNrs();
  auto nrs2 = m2.getAtomicNrs();
  for (int i = 0; i < nat1; i++)
    typ1[static_cast<size_t>(i)] = nrs1[i];
  for (int i = 0; i < nat2; i++)
    typ2[static_cast<size_t>(i)] = nrs2[i];
  return matchArrays(nat1, typ1.data(), m1.getPositions().data(), nat2,
                     typ2.data(), m2.getPositions().data(), distThreshold);
}

matchArrays() [1/2]

IRACompare::MatchResult eonc::IRACompare::matchArrays ( int nat1, const int * typ1, const double * pos1, int nat2, const int * typ2, const double * pos2, double distThreshold )

static

Same as match(), from packed (n,3) row-major coordinates and Z arrays.

Used by the Python server so rot_match does not need a dummy Potential. Production: IRAResource::instance() when built with WITH_IRA.

Definition at line 47 of file IRACompare.cpp.

                                                                      {
#ifdef WITH_IRA
  return matchArrays(nat1, typ1, pos1, nat2, typ2, pos2, distThreshold,
                     get_ira_resource());
#else
  MatchResult result;
  result.error = -1;
  return result;
#endif
}

matchArrays() [2/2]

IRACompare::MatchResult eonc::IRACompare::matchArrays ( int nat1, const int * typ1, const double * pos1, int nat2, const int * typ2, const double * pos2, double distThreshold, IIRAResource & res )

static

Test seam: injected resource, no process-default libira load.

Definition at line 63 of file IRACompare.cpp.

                                                                 {
  MatchResult result;
  std::lock_guard<std::mutex> lock(res.library_mutex);
 
  try {
    res.require_loaded();
  } catch (const std::exception &e) {
    result.error = -1;
    return result;
  }
 
  if (nat1 <= 0 || nat2 <= 0 || typ1 == nullptr || typ2 == nullptr ||
      pos1 == nullptr || pos2 == nullptr) {
    result.error = -1;
    return result;
  }
 
  // Use Eigen::Map to reinterpret the row-major data as column-major for
  // Fortran
  Eigen::Map<const AtomMatrixF> coords1_map(pos1, 3, nat1);
  Eigen::Map<const AtomMatrixF> coords2_map(pos2, 3, nat2);
 
  // Candidate arrays: -1 means "use geometric center as origin" (good
  // initial guess for translation with equal-size structures)
  std::vector<int> cand1(nat1, 0), cand2(nat2, 0);
  if (nat1 == nat2) {
    cand1[0] = -1;
    cand2[0] = -1;
  } else {
    cand1[0] = 1;
    for (int i = 0; i < nat2; i++)
      cand2[i] = i + 1;
  }
 
  // Pre-allocate output buffers (libira_match fills these)
  std::vector<double> rmat_buf(9);
  std::vector<double> tr_buf(3);
  std::vector<int> perm_buf(nat2);
  double hd = 0.0;
  int ierr = 0;
 
  // libira_match expects double** for rotation and translation
  // (it may reallocate internally)
  double *rmat_ptr = rmat_buf.data();
  double *tr_ptr = tr_buf.data();
  int *perm_ptr = perm_buf.data();
 
  res.get_match_fn()(nat1, typ1, coords1_map.data(), cand1.data(), nat2, typ2,
                     coords2_map.data(), cand2.data(), distThreshold, &rmat_ptr,
                     &tr_ptr, &perm_ptr, &hd, &ierr);
 
  result.error = ierr;
  if (ierr == 0) {
    result.permutation.assign(perm_ptr, perm_ptr + nat2);
    result.hausdorffDistance = hd;
    // rmat is 3x3 flat array (column-major from Fortran)
    // Use Eigen maps for efficient conversion
    Eigen::Map<const Matrix3d> rot_map(rmat_ptr);
    result.rotation =
        rot_map.transpose(); // Transpose to get proper row-major layout
    result.translation = Eigen::Map<const Vector3d>(tr_ptr);
  }
 
  // If libira reallocated (pointers changed), free the new buffers.
  // libira allocates replacement output buffers via c_malloc on the Fortran
  // side (see libira's lib_match wrapper), so std::free is the matching
  // deallocator. Never free the original stack/vector-backed buffers.
  if (rmat_ptr != rmat_buf.data())
    std::free(rmat_ptr);
  if (tr_ptr != tr_buf.data())
    std::free(tr_ptr);
  if (perm_ptr != perm_buf.data())
    std::free(perm_ptr);
  return result;
}

matchPBC() [1/2]

IRACompare::MatchResult eonc::IRACompare::matchPBC ( const Matter & m1, const Matter & m2, double distThreshold )

static

Atom assignment under periodic boundary conditions (CShDA only, no rotation/SVD).

Returns permutation and per-atom distances.

Definition at line 141 of file IRACompare.cpp.

                                                                   {
#ifdef WITH_IRA
  return matchPBC(m1, m2, distThreshold, get_ira_resource());
#else
  MatchResult result;
  result.error = -1;
  return result;
#endif
}

matchPBC() [2/2]

IRACompare::MatchResult eonc::IRACompare::matchPBC ( const Matter & m1, const Matter & m2, double distThreshold, IIRAResource & res )

static

Definition at line 152 of file IRACompare.cpp.

                                                                {
  MatchResult result;
  std::lock_guard<std::mutex> lock(res.library_mutex);
 
  try {
    res.require_loaded();
  } catch (const std::exception &e) {
    result.error = -1;
    return result;
  }
 
  const int nat1 = m1.numberOfAtoms();
  const int nat2 = m2.numberOfAtoms();
  if (nat1 <= 0 || nat2 <= 0) {
    result.error = -1;
    return result;
  }
 
  std::vector<int> typ1(nat1), typ2(nat2);
  auto nrs1 = m1.getAtomicNrs();
  auto nrs2 = m2.getAtomicNrs();
  for (int i = 0; i < nat1; i++)
    typ1[i] = nrs1[i];
  for (int i = 0; i < nat2; i++)
    typ2[i] = nrs2[i];
 
  // Prepare coordinates using direct Eigen data access
  const AtomMatrix &pos1 = m1.getPositions();
  const AtomMatrix &pos2 = m2.getPositions();
 
  Eigen::Map<const AtomMatrixF> coords1_map(pos1.data(), 3, nat1);
  Eigen::Map<const AtomMatrixF> coords2_map(pos2.data(), 3, nat2);
 
  // Lattice vectors (3x3 column-major)
  Matrix3d cell = m2.getCell();
  double lat[9];
  for (int i = 0; i < 3; i++)
    for (int j = 0; j < 3; j++)
      lat[j * 3 + i] = cell(i, j);
 
  // Caller must pre-allocate output arrays; Fortran receives pointers only
  std::vector<int> found_buf(nat1);
  std::vector<double> dists_buf(nat1);
  int *found_ptr = found_buf.data();
  double *dists_ptr = dists_buf.data();
 
  // NOTE: this calls cshda_pbc only (assignment, no rotation/SVD matching).
  // It finds the best atom assignment under PBC but does not compute the
  // optimal rotation or translation.
  res.get_cshda_pbc_fn()(nat1, typ1.data(), coords1_map.data(), nat2,
                         typ2.data(), coords2_map.data(), lat, distThreshold,
                         &found_ptr, &dists_ptr);
 
  result.permutation.assign(found_ptr, found_ptr + nat1);
  result.hausdorffDistance = 0.0;
  for (int i = 0; i < nat1; i++) {
    result.hausdorffDistance = std::max(result.hausdorffDistance, dists_ptr[i]);
  }
  result.rotation = Eigen::Matrix3d::Identity();
  result.translation = Eigen::Vector3d::Zero();
  result.error = 0;
  return result;
}

---

The documentation for this class was generated from the following files:

• /home/runner/work/eOn/eOn/include/eon/IRACompare.h
• /home/runner/work/eOn/eOn/client/IRACompare.cpp