eonc::geometry Namespace Reference

Classes
struct	atom
struct	by_atom

Functions
RotationMatrix	rotationExtract (const AtomMatrix r1, const AtomMatrix r2)
bool	rotationMatch (const Matter &m1, const Matter &m2, const double max_diff)
void	projectOutRotTrans (Eigen::VectorXd &step, const AtomMatrix &positions)
void	rotationRemove (const AtomMatrix r1, std::shared_ptr< Matter > m2)
void	rotationRemove (const std::shared_ptr< Matter > m1, std::shared_ptr< Matter > m2)
void	translationRemove (Matter &m1, const AtomMatrix r1)
void	translationRemove (Matter &m1, const Matter &m2)
double	maxAtomMotion (const AtomMatrix v1)
double	maxAtomMotionV (const VectorXd v1)
long	numAtomsMoved (const AtomMatrix v1, double cutoff)
AtomMatrix	maxAtomMotionApplied (const AtomMatrix v1, double maxMotion)
VectorXd	maxAtomMotionAppliedV (const VectorXd v1, double maxMotion)
AtomMatrix	maxMotionApplied (const AtomMatrix v1, double maxMotion)
VectorXd	maxMotionAppliedV (const VectorXd v1, double maxMotion)
bool	identical (const Matter &m1, const Matter &m2, const double distanceDifference)
bool	sortedR (const Matter &m1, const Matter &m2, const double distanceDifference)
void	pushApart (std::shared_ptr< Matter > m1, double minDistance)

Function Documentation

identical()

bool eonc::geometry::identical	(	const Matter &	m1,
		const Matter &	m2,
		const double	distanceDifference )

Definition at line 369 of file GeometryAnalysis.cpp.

                                                                {
 
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = m2.getPositions();
 
  std::set<int> matched;
  double tolerance = distanceDifference;
 
  if (r1.rows() != r2.rows()) {
    return false;
  }
  int N = r1.rows();
 
  for (int i = 0; i < N; i++) {
    if (std::fabs((m1.pbc(r1.row(i) - r2.row(i))).norm()) < tolerance &&
        m1.getAtomicNr(i) == m2.getAtomicNr(i)) {
      matched.insert(i);
    }
  }
 
  for (int j = 0; j < N; j++) {
 
    if (matched.count(j) == 1)
      continue;
 
    for (int k = 0; k < N; k++) {
      if (matched.count(j) == 1)
        break;
 
      if (std::fabs((m1.pbc(r1.row(j) - r2.row(k))).norm()) < tolerance &&
          m1.getAtomicNr(j) == m2.getAtomicNr(k)) {
        matched.insert(j);
      }
    }
  }
 
  if (matched.size() == static_cast<unsigned>(N)) {
    return true;
  } else {
    return false;
  }
}

maxAtomMotion()

double eonc::geometry::maxAtomMotion

(

const AtomMatrix

v1

)

Definition at line 269 of file GeometryAnalysis.cpp.

                                                      {
  return v1.rowwise().norm().maxCoeff();
}

maxAtomMotionApplied()

AtomMatrix eonc::geometry::maxAtomMotionApplied	(	const AtomMatrix	v1,
		double	maxMotion )

Definition at line 308 of file GeometryAnalysis.cpp.

                                                                  {
  AtomMatrix v2(v1);
 
  double max = maxAtomMotion(v1);
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

maxAtomMotionAppliedV()

VectorXd eonc::geometry::maxAtomMotionAppliedV	(	const VectorXd	v1,
		double	maxMotion )

Definition at line 319 of file GeometryAnalysis.cpp.

                                                                 {
  VectorXd v2(v1);
 
  double max = maxAtomMotionV(v1);
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

maxAtomMotionV()

double eonc::geometry::maxAtomMotionV

(

const VectorXd

v1

)

Definition at line 273 of file GeometryAnalysis.cpp.

                                                     {
  double max = 0.0;
  long n = v1.rows();
  if (n < 3) {
    // Vector too short for 3D atom grouping; treat as single displacement
    return v1.norm();
  }
  for (long i = 0; i + 3 <= n; i += 3) {
    double norm = v1.segment<3>(i).norm();
    if (max < norm) {
      max = norm;
    }
  }
  // Handle trailing elements (vector size not multiple of 3)
  long rem = n % 3;
  if (rem > 0) {
    double norm = v1.tail(rem).norm();
    if (max < norm) {
      max = norm;
    }
  }
  return max;
}

maxMotionApplied()

AtomMatrix eonc::geometry::maxMotionApplied	(	const AtomMatrix	v1,
		double	maxMotion )

Definition at line 330 of file GeometryAnalysis.cpp.

                                                              {
  AtomMatrix v2(v1);
 
  double max = v1.norm();
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

maxMotionAppliedV()

VectorXd eonc::geometry::maxMotionAppliedV	(	const VectorXd	v1,
		double	maxMotion )

Definition at line 341 of file GeometryAnalysis.cpp.

                                                             {
  VectorXd v2(v1);
 
  double max = v1.norm();
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

numAtomsMoved()

long eonc::geometry::numAtomsMoved	(	const AtomMatrix	v1,
		double	cutoff )

Definition at line 297 of file GeometryAnalysis.cpp.

                                                                   {
  long num = 0;
  for (int i = 0; i < v1.rows(); i++) {
    double norm = v1.row(i).norm();
    if (norm >= cutoff) {
      num += 1;
    }
  }
  return num;
}

projectOutRotTrans()

void eonc::geometry::projectOutRotTrans	(	Eigen::VectorXd &	step,
		const AtomMatrix &	positions )

Definition at line 132 of file GeometryAnalysis.cpp.

                                                                     {
  long nAtoms = positions.rows();
  long dof = nAtoms * 3;
 
  // 1. Compute center of mass (unweighted geometric center)
  Eigen::Vector3d com = Eigen::Vector3d::Zero();
  for (long i = 0; i < nAtoms; ++i) {
    com(0) += positions(i, 0);
    com(1) += positions(i, 1);
    com(2) += positions(i, 2);
  }
  com /= static_cast<double>(nAtoms);
 
  // 2. Construct 6 basis vectors for rigid-body translation and rotation
  std::vector<Eigen::VectorXd> basis;
  basis.reserve(6);
 
  // Translational basis vectors
  for (int d = 0; d < 3; ++d) {
    Eigen::VectorXd t = Eigen::VectorXd::Zero(dof);
    for (long j = 0; j < nAtoms; ++j) {
      t(3 * j + d) = 1.0;
    }
    basis.push_back(t);
  }
 
  // Rotational basis vectors (infinitesimal rotations around COM)
  Eigen::VectorXd rx = Eigen::VectorXd::Zero(dof);
  Eigen::VectorXd ry = Eigen::VectorXd::Zero(dof);
  Eigen::VectorXd rz = Eigen::VectorXd::Zero(dof);
 
  for (long i = 0; i < nAtoms; ++i) {
    double x = positions(i, 0) - com(0);
    double y = positions(i, 1) - com(1);
    double z = positions(i, 2) - com(2);
    // Rotation around x-axis: cross(xhat, r) = (0, -z, y)
    rx(3 * i + 1) = -z;
    rx(3 * i + 2) = y;
    // Rotation around y-axis: cross(yhat, r) = (z, 0, -x)
    ry(3 * i + 0) = z;
    ry(3 * i + 2) = -x;
    // Rotation around z-axis: cross(zhat, r) = (-y, x, 0)
    rz(3 * i + 0) = -y;
    rz(3 * i + 1) = x;
  }
  basis.push_back(rx);
  basis.push_back(ry);
  basis.push_back(rz);
 
  // 3. Modified Gram-Schmidt orthonormalization
  std::vector<Eigen::VectorXd> ortho;
  ortho.reserve(6);
 
  for (auto &v : basis) {
    Eigen::VectorXd u = v;
    for (const auto &e : ortho) {
      u -= u.dot(e) * e;
    }
    // Handles linear molecules where one rotational mode is degenerate
    if (u.norm() > 1e-9) {
      u.normalize();
      ortho.push_back(u);
    }
  }
 
  // 4. Project out rigid-body components from step
  for (const auto &e : ortho) {
    step -= step.dot(e) * e;
  }
}

pushApart()

void eonc::geometry::pushApart	(	std::shared_ptr< Matter >	m1,
		double	minDistance )

Definition at line 485 of file GeometryAnalysis.cpp.

                                                                         {
  if (minDistance <= 0)
    return;
 
  AtomMatrix r1 = m1->getPositions();
  AtomMatrix Force(r1.rows(), 3);
  double f = 0.025;
  double cut = minDistance;
  double pushAparts = 500;
  for (int p = 0; p < r1.rows(); p++) {
    for (int axis = 0; axis <= 2; axis++) {
      Force(p, axis) = 0;
    }
  }
  for (int count = 0; count < pushAparts; count++) {
    int moved = 0;
    for (int i = 0; i < r1.rows(); i++) {
      for (int j = i + 1; j < r1.rows(); j++) {
        double d = m1->distance(i, j);
        if (d < cut) {
          moved++;
          for (int axis = 0; axis <= 2; axis++) {
            double componant = f * (r1(i, axis) - r1(j, axis)) / d;
            Force(i, axis) += componant;
            Force(j, axis) -= componant;
          }
        }
      }
    }
    if (moved == 0)
      break;
    for (int k = 0; k < r1.rows(); k++) {
      for (int axis = 0; axis <= 2; axis++) {
        r1(k, axis) += Force(k, axis);
        Force(k, axis) = 0;
      }
    }
    m1->setPositions(r1);
  }
}

rotationExtract()

RotationMatrix eonc::geometry::rotationExtract	(	const AtomMatrix	r1,
		const AtomMatrix	r2 )

Definition at line 21 of file GeometryAnalysis.cpp.

                                                                    {
  RotationMatrix R;
 
  // Determine optimal rotation
  // Horn, J. Opt. Soc. Am. A, 1987
  Matrix3d m = r1.transpose() * r2;
 
  double sxx = m(0, 0);
  double sxy = m(0, 1);
  double sxz = m(0, 2);
  double syx = m(1, 0);
  double syy = m(1, 1);
  double syz = m(1, 2);
  double szx = m(2, 0);
  double szy = m(2, 1);
  double szz = m(2, 2);
 
  Matrix4d n;
  n.setZero();
  n(0, 1) = syz - szy;
  n(0, 2) = szx - sxz;
  n(0, 3) = sxy - syx;
 
  n(1, 2) = sxy + syx;
  n(1, 3) = szx + sxz;
 
  n(2, 3) = syz + szy;
 
  n += n.transpose().eval();
 
  n(0, 0) = sxx + syy + szz;
  n(1, 1) = sxx - syy - szz;
  n(2, 2) = -sxx + syy - szz;
  n(3, 3) = -sxx - syy + szz;
 
  Eigen::SelfAdjointEigenSolver<Matrix4d> es(n);
  Eigen::Vector4d maxv = es.eigenvectors().col(3);
 
  double aa = maxv[0] * maxv[0];
  double bb = maxv[1] * maxv[1];
  double cc = maxv[2] * maxv[2];
  double dd = maxv[3] * maxv[3];
  double ab = maxv[0] * maxv[1];
  double ac = maxv[0] * maxv[2];
  double ad = maxv[0] * maxv[3];
  double bc = maxv[1] * maxv[2];
  double bd = maxv[1] * maxv[3];
  double cd = maxv[2] * maxv[3];
 
  R(0, 0) = aa + bb - cc - dd;
  R(0, 1) = 2 * (bc - ad);
  R(0, 2) = 2 * (bd + ac);
  R(1, 0) = 2 * (bc + ad);
  R(1, 1) = aa - bb + cc - dd;
  R(1, 2) = 2 * (cd - ab);
  R(2, 0) = 2 * (bd - ac);
  R(2, 1) = 2 * (cd + ab);
  R(2, 2) = aa - bb - cc + dd;
 
  return R;
}

rotationMatch()

bool eonc::geometry::rotationMatch	(	const Matter &	m1,
		const Matter &	m2,
		const double	max_diff )

Definition at line 84 of file GeometryAnalysis.cpp.

                                                          {
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = m2.getPositions();
 
  // Align centroids
  Eigen::VectorXd c1(3);
  Eigen::VectorXd c2(3);
 
  c1[0] = r1.col(0).sum();
  c1[1] = r1.col(1).sum();
  c1[2] = r1.col(2).sum();
  c2[0] = r2.col(0).sum();
  c2[1] = r2.col(1).sum();
  c2[2] = r2.col(2).sum();
  c1 /= r1.rows();
  c2 /= r2.rows();
 
  for (int i = 0; i < r1.rows(); i++) {
    r1(i, 0) -= c1[0];
    r1(i, 1) -= c1[1];
    r1(i, 2) -= c1[2];
 
    r2(i, 0) -= c2[0];
    r2(i, 1) -= c2[1];
    r2(i, 2) -= c2[2];
  }
 
  RotationMatrix R = rotationExtract(r1, r2);
 
  // Eigen is transposed relative to numpy
  r2 = r2 * R;
 
  for (int i = 0; i < r1.rows(); i++) {
    double diff = (r2.row(i) - r1.row(i)).norm();
    if (diff > max_diff) {
      return false;
    }
  }
  return true;
}

rotationRemove() [1/2]

void eonc::geometry::rotationRemove	(	const AtomMatrix	r1,
		std::shared_ptr< Matter >	m2 )

Definition at line 204 of file GeometryAnalysis.cpp.

                                                              {
  // Skip for extended systems (slabs/surfaces with frozen atoms).
  // Rigid-body rotation and translation are not well-defined when the
  // system is anchored by frozen atoms.
  if (m2->numberOfFixedAtoms() > 0) {
    return;
  }
 
  AtomMatrix r2 = m2->getPositions();
  long n = r1_passed.size();
 
  // Compute displacement as flat 3N vector
  Eigen::VectorXd step(n);
  Eigen::Map<const Eigen::VectorXd> r1_flat(r1_passed.data(), n);
  Eigen::Map<const Eigen::VectorXd> r2_flat(r2.data(), n);
  step = r2_flat - r1_flat;
 
  // Project out rigid-body translation and rotation
  projectOutRotTrans(step, r1_passed);
 
  // Reconstruct positions: r1 + projected step
  Eigen::VectorXd result = r1_flat + step;
  AtomMatrix resultMat(r1_passed.rows(), 3);
  Eigen::Map<Eigen::VectorXd>(resultMat.data(), n) = result;
 
  m2->setPositions(resultMat);
}

rotationRemove() [2/2]

void eonc::geometry::rotationRemove	(	const std::shared_ptr< Matter >	m1,
		std::shared_ptr< Matter >	m2 )

Definition at line 233 of file GeometryAnalysis.cpp.

                                                              {
  AtomMatrix r1 = m1->getPositions();
  rotationRemove(r1, m2);
  return;
}

sortedR()

bool eonc::geometry::sortedR	(	const Matter &	m1,
		const Matter &	m2,
		const double	distanceDifference )

Definition at line 413 of file GeometryAnalysis.cpp.

                                                              {
  EONC_LOG_INFO("In sortedR");
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = m2.getPositions();
  double tolerance = distanceDifference;
  int matches = 0;
 
  if (r1.rows() != r2.rows()) {
    return false;
  }
 
  // Allocate memory for rdf1 and rdf2
  std::vector<std::set<atom, by_atom>> rdf1(r1.rows());
  std::vector<std::set<atom, by_atom>> rdf2(r2.rows());
 
  for (int i2 = 0; i2 < r2.rows(); i2++) {
    rdf2[i2].clear();
    for (int j2 = 0; j2 < r2.rows(); j2++) {
      if (j2 == i2)
        continue;
      atom a2;
      a2.r = m2.distance(i2, j2);
      a2.z = m2.getAtomicNr(j2);
      rdf2[i2].insert(a2);
      rdf2[j2].insert(a2);
    }
  }
 
  for (int i1 = 0; i1 < r1.rows(); i1++) {
    if (matches == i1 - 2) {
      return false;
    }
    for (int j1 = 0; j1 < r1.rows(); j1++) {
      if (j1 == i1)
        continue;
      atom a;
      a.r = m1.distance(i1, j1);
      a.z = m1.getAtomicNr(j1);
      rdf1[i1].insert(a);
      rdf1[j1].insert(a);
    }
    for (int x = 0; x < r2.rows(); x++) {
      auto it2 = rdf2[x].begin();
      auto it = rdf1[i1].begin();
      int c = 0;
      int counter = 0;
      for (; c < r1.rows(); c++) {
        if (it == rdf1[i1].end() || it2 == rdf2[x].end())
          break;
        atom k1 = *it;
        atom k2 = *it2;
        if (std::fabs(k1.r - k2.r) < tolerance && k1.z == k2.z) {
          counter++;
        } else {
          EONC_LOG_INFO("No match");
          break;
        }
        ++it;
        ++it2;
      }
      if (counter == r1.rows()) {
        matches++;
      } else {
        EONC_LOG_INFO("No match");
      }
    }
  }
 
  return matches >= r1.rows();
}

translationRemove() [1/2]

void eonc::geometry::translationRemove	(	Matter &	m1,
		const AtomMatrix	r1 )

Definition at line 240 of file GeometryAnalysis.cpp.

                                                                           {
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = r2_passed;
 
  // net displacement
  Eigen::VectorXd disp(3);
  AtomMatrix r12 = m1.pbc(r2 - r1);
 
  disp[0] = r12.col(0).sum();
  disp[1] = r12.col(1).sum();
  disp[2] = r12.col(2).sum();
  disp /= r1.rows();
 
  for (int i = 0; i < r1.rows(); i++) {
    r1(i, 0) += disp[0];
    r1(i, 1) += disp[1];
    r1(i, 2) += disp[2];
  }
 
  m1.setPositions(r1);
  return;
}

translationRemove() [2/2]

void eonc::geometry::translationRemove	(	Matter &	m1,
		const Matter &	m2 )

Definition at line 263 of file GeometryAnalysis.cpp.

                                                                 {
  AtomMatrix r2 = m2.getPositions();
  translationRemove(m1, r2);
  return;
}