Loading...
Searching...
No Matches
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 397 of file GeometryAnalysis.cpp.

398 {
399 AtomMatrix r1 = m1.getPositions();
400 AtomMatrix r2 = m2.getPositions();
401 if (r1.rows() != r2.rows()) {
402 return false;
403 }
404 const int nAtoms = static_cast<int>(r1.rows());
405 const double tolerance = distanceDifference;
406
407 bool indexAligned = true;
408 for (int i = 0; i < nAtoms; i++) {
409 if (!pairInside(m1, m2, r1, r2, i, i, tolerance)) {
410 indexAligned = false;
411 break;
412 }
413 }
414 if (indexAligned) {
415 return true;
416 }
417
418 std::vector<std::vector<int>> adj(static_cast<size_t>(nAtoms));
419 for (int left = 0; left < nAtoms; left++) {
420 for (int right = 0; right < nAtoms; right++) {
421 if (pairInside(m1, m2, r1, r2, left, right, tolerance)) {
422 adj[static_cast<size_t>(left)].push_back(right);
423 }
424 }
425 }
426
427 std::vector<int> matchRight(static_cast<size_t>(nAtoms), -1);
428 int matched = 0;
429 for (int left = 0; left < nAtoms; left++) {
430 std::vector<char> seen(static_cast<size_t>(nAtoms), 0);
431 if (identicalAugment(left, adj, matchRight, seen)) {
432 matched++;
433 }
434 }
435 return matched == nAtoms;
436}
AtomMatrix
Eigen::Matrix< double, Eigen::Dynamic, 3, eOnStorageOrder > AtomMatrix
Definition Eigen.h:37
eonc::Matter::getPositions
const AtomMatrix & getPositions() const
Definition Matter.cpp:308

◆ maxAtomMotion()

double eonc::geometry::maxAtomMotion ( const AtomMatrix v1)

Definition at line 269 of file GeometryAnalysis.cpp.

269 {
270 return v1.rowwise().norm().maxCoeff();
271}

◆ maxAtomMotionApplied()

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

Definition at line 308 of file GeometryAnalysis.cpp.

309 {
310 AtomMatrix v2(v1);
311
312 double max = maxAtomMotion(v1);
313 if (max > maxMotion) {
314 v2 *= maxMotion / max;
315 }
316 return v2;
317}
eonc::geometry::maxAtomMotion
double maxAtomMotion(const AtomMatrix v1)
Definition GeometryAnalysis.cpp:269

◆ maxAtomMotionAppliedV()

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

Definition at line 319 of file GeometryAnalysis.cpp.

320 {
321 VectorXd v2(v1);
322
323 double max = maxAtomMotionV(v1);
324 if (max > maxMotion) {
325 v2 *= maxMotion / max;
326 }
327 return v2;
328}
eonc::geometry::maxAtomMotionV
double maxAtomMotionV(const VectorXd v1)
Definition GeometryAnalysis.cpp:273

◆ maxAtomMotionV()

double eonc::geometry::maxAtomMotionV ( const VectorXd v1)

Definition at line 273 of file GeometryAnalysis.cpp.

273 {
274 double max = 0.0;
275 long n = v1.rows();
276 if (n < 3) {
277 // Vector too short for 3D atom grouping; treat as single displacement
278 return v1.norm();
279 }
280 for (long i = 0; i + 3 <= n; i += 3) {
281 double norm = v1.segment<3>(i).norm();
282 if (max < norm) {
283 max = norm;
284 }
285 }
286 // Handle trailing elements (vector size not multiple of 3)
287 long rem = n % 3;
288 if (rem > 0) {
289 double norm = v1.tail(rem).norm();
290 if (max < norm) {
291 max = norm;
292 }
293 }
294 return max;
295}

◆ maxMotionApplied()

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

Definition at line 330 of file GeometryAnalysis.cpp.

331 {
332 AtomMatrix v2(v1);
333
334 double max = v1.norm();
335 if (max > maxMotion) {
336 v2 *= maxMotion / max;
337 }
338 return v2;
339}

◆ maxMotionAppliedV()

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

Definition at line 341 of file GeometryAnalysis.cpp.

342 {
343 VectorXd v2(v1);
344
345 double max = v1.norm();
346 if (max > maxMotion) {
347 v2 *= maxMotion / max;
348 }
349 return v2;
350}

◆ numAtomsMoved()

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

Definition at line 297 of file GeometryAnalysis.cpp.

297 {
298 long num = 0;
299 for (int i = 0; i < v1.rows(); i++) {
300 double norm = v1.row(i).norm();
301 if (norm >= cutoff) {
302 num += 1;
303 }
304 }
305 return num;
306}

◆ projectOutRotTrans()

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

Definition at line 132 of file GeometryAnalysis.cpp.

133 {
134 long nAtoms = positions.rows();
135 long dof = nAtoms * 3;
136
137 // 1. Compute center of mass (unweighted geometric center)
138 Eigen::Vector3d com = Eigen::Vector3d::Zero();
139 for (long i = 0; i < nAtoms; ++i) {
140 com(0) += positions(i, 0);
141 com(1) += positions(i, 1);
142 com(2) += positions(i, 2);
143 }
144 com /= static_cast<double>(nAtoms);
145
146 // 2. Construct 6 basis vectors for rigid-body translation and rotation
147 std::vector<Eigen::VectorXd> basis;
148 basis.reserve(6);
149
150 // Translational basis vectors
151 for (int d = 0; d < 3; ++d) {
152 Eigen::VectorXd t = Eigen::VectorXd::Zero(dof);
153 for (long j = 0; j < nAtoms; ++j) {
154 t(3 * j + d) = 1.0;
155 }
156 basis.push_back(t);
157 }
158
159 // Rotational basis vectors (infinitesimal rotations around COM)
160 Eigen::VectorXd rx = Eigen::VectorXd::Zero(dof);
161 Eigen::VectorXd ry = Eigen::VectorXd::Zero(dof);
162 Eigen::VectorXd rz = Eigen::VectorXd::Zero(dof);
163
164 for (long i = 0; i < nAtoms; ++i) {
165 double x = positions(i, 0) - com(0);
166 double y = positions(i, 1) - com(1);
167 double z = positions(i, 2) - com(2);
168 // Rotation around x-axis: cross(xhat, r) = (0, -z, y)
169 rx(3 * i + 1) = -z;
170 rx(3 * i + 2) = y;
171 // Rotation around y-axis: cross(yhat, r) = (z, 0, -x)
172 ry(3 * i + 0) = z;
173 ry(3 * i + 2) = -x;
174 // Rotation around z-axis: cross(zhat, r) = (-y, x, 0)
175 rz(3 * i + 0) = -y;
176 rz(3 * i + 1) = x;
177 }
178 basis.push_back(rx);
179 basis.push_back(ry);
180 basis.push_back(rz);
181
182 // 3. Modified Gram-Schmidt orthonormalization
183 std::vector<Eigen::VectorXd> ortho;
184 ortho.reserve(6);
185
186 for (auto &v : basis) {
187 Eigen::VectorXd u = v;
188 for (const auto &e : ortho) {
189 u -= u.dot(e) * e;
190 }
191 // Handles linear molecules where one rotational mode is degenerate
192 if (u.norm() > 1e-9) {
193 u.normalize();
194 ortho.push_back(u);
195 }
196 }
197
198 // 4. Project out rigid-body components from step
199 for (const auto &e : ortho) {
200 step -= step.dot(e) * e;
201 }
202}

◆ pushApart()

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

Definition at line 507 of file GeometryAnalysis.cpp.

507 {
508 if (minDistance <= 0)
509 return;
510
511 AtomMatrix r1 = m1->getPositions();
512 AtomMatrix Force(r1.rows(), 3);
513 double f = 0.025;
514 double cut = minDistance;
515 double pushAparts = 500;
516 Force.setZero();
517 for (int count = 0; count < pushAparts; count++) {
518 int moved = 0;
519 for (int i = 0; i < r1.rows(); i++) {
520 for (int j = i + 1; j < r1.rows(); j++) {
521 // distance() is the minimum-image length. The step has to use that
522 // same vector: the raw Cartesian separation is the long box diagonal
523 // when the close pair straddles a periodic face.
524 AtomMatrix delta(1, 3);
525 delta.row(0) = r1.row(i) - r1.row(j);
526 delta = m1->pbc(delta);
527 const double d = delta.norm();
528 if (!(d < cut)) {
529 continue;
530 }
531 moved++;
532 // A zero minimum-image vector has no direction. Take the finite
533 // step a unit separation would, along +x, instead of dividing.
534 if (!(d > eonc::safemath::eps)) {
535 Force(i, 0) += f;
536 Force(j, 0) -= f;
537 continue;
538 }
539 for (int axis = 0; axis <= 2; axis++) {
540 const double component = f * delta(0, axis) / d;
541 Force(i, axis) += component;
542 Force(j, axis) -= component;
543 }
544 }
545 }
546 if (moved == 0)
547 break;
548 r1 += Force;
549 Force.setZero();
550 m1->setPositions(r1);
551 // setPositions wraps when periodic boundaries are on. The next pass
552 // has to see those wrapped coordinates, not the pre-wrap copy.
553 r1 = m1->getPositions();
554 }
555}
eonc::safemath::eps
constexpr double eps
Definition SafeMath.h:19

◆ rotationExtract()

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

Definition at line 21 of file GeometryAnalysis.cpp.

22 {
23 RotationMatrix R;
24
25 // Determine optimal rotation
26 // Horn, J. Opt. Soc. Am. A, 1987
27 Matrix3d m = r1.transpose() * r2;
28
29 double sxx = m(0, 0);
30 double sxy = m(0, 1);
31 double sxz = m(0, 2);
32 double syx = m(1, 0);
33 double syy = m(1, 1);
34 double syz = m(1, 2);
35 double szx = m(2, 0);
36 double szy = m(2, 1);
37 double szz = m(2, 2);
38
39 Matrix4d n;
40 n.setZero();
41 n(0, 1) = syz - szy;
42 n(0, 2) = szx - sxz;
43 n(0, 3) = sxy - syx;
44
45 n(1, 2) = sxy + syx;
46 n(1, 3) = szx + sxz;
47
48 n(2, 3) = syz + szy;
49
50 n += n.transpose().eval();
51
52 n(0, 0) = sxx + syy + szz;
53 n(1, 1) = sxx - syy - szz;
54 n(2, 2) = -sxx + syy - szz;
55 n(3, 3) = -sxx - syy + szz;
56
57 Eigen::SelfAdjointEigenSolver<Matrix4d> es(n);
58 Eigen::Vector4d maxv = es.eigenvectors().col(3);
59
60 double aa = maxv[0] * maxv[0];
61 double bb = maxv[1] * maxv[1];
62 double cc = maxv[2] * maxv[2];
63 double dd = maxv[3] * maxv[3];
64 double ab = maxv[0] * maxv[1];
65 double ac = maxv[0] * maxv[2];
66 double ad = maxv[0] * maxv[3];
67 double bc = maxv[1] * maxv[2];
68 double bd = maxv[1] * maxv[3];
69 double cd = maxv[2] * maxv[3];
70
71 R(0, 0) = aa + bb - cc - dd;
72 R(0, 1) = 2 * (bc - ad);
73 R(0, 2) = 2 * (bd + ac);
74 R(1, 0) = 2 * (bc + ad);
75 R(1, 1) = aa - bb + cc - dd;
76 R(1, 2) = 2 * (cd - ab);
77 R(2, 0) = 2 * (bd - ac);
78 R(2, 1) = 2 * (cd + ab);
79 R(2, 2) = aa - bb - cc + dd;
80
81 return R;
82}
RotationMatrix
Eigen::Matrix< double, 3, 3, eOnStorageOrder > RotationMatrix
Definition Eigen.h:38
Matrix3d
Eigen::Matrix< double, 3, 3, eOnStorageOrder > Matrix3d
Definition Eigen.h:35
Matrix4d
Eigen::Matrix< double, 4, 4, eOnStorageOrder > Matrix4d
Definition Eigen.h:36

◆ rotationMatch()

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

Definition at line 84 of file GeometryAnalysis.cpp.

85 {
86 AtomMatrix r1 = m1.getPositions();
87 AtomMatrix r2 = m2.getPositions();
88
89 // Align centroids
90 Eigen::VectorXd c1(3);
91 Eigen::VectorXd c2(3);
92
93 c1[0] = r1.col(0).sum();
94 c1[1] = r1.col(1).sum();
95 c1[2] = r1.col(2).sum();
96 c2[0] = r2.col(0).sum();
97 c2[1] = r2.col(1).sum();
98 c2[2] = r2.col(2).sum();
99 c1 /= r1.rows();
100 c2 /= r2.rows();
101
102 for (int i = 0; i < r1.rows(); i++) {
103 r1(i, 0) -= c1[0];
104 r1(i, 1) -= c1[1];
105 r1(i, 2) -= c1[2];
106
107 r2(i, 0) -= c2[0];
108 r2(i, 1) -= c2[1];
109 r2(i, 2) -= c2[2];
110 }
111
112 RotationMatrix R = rotationExtract(r1, r2);
113
114 // Eigen is transposed relative to numpy
115 r2 = r2 * R;
116
117 for (int i = 0; i < r1.rows(); i++) {
118 double diff = (r2.row(i) - r1.row(i)).norm();
119 if (diff > max_diff) {
120 return false;
121 }
122 }
123 return true;
124}
eonc::geometry::rotationExtract
RotationMatrix rotationExtract(const AtomMatrix r1, const AtomMatrix r2)
Definition GeometryAnalysis.cpp:21

◆ rotationRemove() [1/2]

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

Definition at line 204 of file GeometryAnalysis.cpp.

205 {
206 // Skip for extended systems (slabs/surfaces with frozen atoms).
207 // Rigid-body rotation and translation are not well-defined when the
208 // system is anchored by frozen atoms.
209 if (m2->numberOfFixedAtoms() > 0) {
210 return;
211 }
212
213 AtomMatrix r2 = m2->getPositions();
214 long n = r1_passed.size();
215
216 // Compute displacement as flat 3N vector
217 Eigen::VectorXd step(n);
218 Eigen::Map<const Eigen::VectorXd> r1_flat(r1_passed.data(), n);
219 Eigen::Map<const Eigen::VectorXd> r2_flat(r2.data(), n);
220 step = r2_flat - r1_flat;
221
222 // Project out rigid-body translation and rotation
223 projectOutRotTrans(step, r1_passed);
224
225 // Reconstruct positions: r1 + projected step
226 Eigen::VectorXd result = r1_flat + step;
227 AtomMatrix resultMat(r1_passed.rows(), 3);
228 Eigen::Map<Eigen::VectorXd>(resultMat.data(), n) = result;
229
230 m2->setPositions(resultMat);
231}
eonc::geometry::projectOutRotTrans
void projectOutRotTrans(Eigen::VectorXd &step, const AtomMatrix &positions)
Definition GeometryAnalysis.cpp:132

◆ 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.

234 {
235 AtomMatrix r1 = m1->getPositions();
236 rotationRemove(r1, m2);
237 return;
238}
eonc::geometry::rotationRemove
void rotationRemove(const AtomMatrix r1, std::shared_ptr< Matter > m2)
Definition GeometryAnalysis.cpp:204

◆ sortedR()

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

Definition at line 438 of file GeometryAnalysis.cpp.

439 {
440 EONC_LOG_INFO("In sortedR");
441 AtomMatrix r1 = m1.getPositions();
442 AtomMatrix r2 = m2.getPositions();
443 double tolerance = distanceDifference;
444 int matches = 0;
445
446 if (r1.rows() != r2.rows()) {
447 return false;
448 }
449
450 // Allocate memory for rdf1 and rdf2
451 std::vector<std::set<atom, by_atom>> rdf1(r1.rows());
452 std::vector<std::set<atom, by_atom>> rdf2(r2.rows());
453
454 for (int i2 = 0; i2 < r2.rows(); i2++) {
455 rdf2[i2].clear();
456 for (int j2 = 0; j2 < r2.rows(); j2++) {
457 if (j2 == i2)
458 continue;
459 atom a2;
460 a2.r = m2.distance(i2, j2);
461 a2.z = m2.getAtomicNr(j2);
462 rdf2[i2].insert(a2);
463 rdf2[j2].insert(a2);
464 }
465 }
466
467 for (int i1 = 0; i1 < r1.rows(); i1++) {
468 if (matches == i1 - 2) {
469 return false;
470 }
471 for (int j1 = 0; j1 < r1.rows(); j1++) {
472 if (j1 == i1)
473 continue;
474 atom a;
475 a.r = m1.distance(i1, j1);
476 a.z = m1.getAtomicNr(j1);
477 rdf1[i1].insert(a);
478 rdf1[j1].insert(a);
479 }
480 for (int x = 0; x < r2.rows(); x++) {
481 auto it2 = rdf2[x].begin();
482 auto it = rdf1[i1].begin();
483 int counter = 0;
484 // Shell length is the neighbor set, which is not the atom count.
485 for (; it != rdf1[i1].end() && it2 != rdf2[x].end(); ++it, ++it2) {
486 atom k1 = *it;
487 atom k2 = *it2;
488 if (std::fabs(k1.r - k2.r) < tolerance && k1.z == k2.z) {
489 counter++;
490 } else {
491 EONC_LOG_INFO("No match");
492 break;
493 }
494 }
495 if (static_cast<size_t>(counter) == rdf1[i1].size() &&
496 it2 == rdf2[x].end()) {
497 matches++;
498 } else {
499 EONC_LOG_INFO("No match");
500 }
501 }
502 }
503
504 return matches >= r1.rows();
505}
EONC_LOG_INFO
#define EONC_LOG_INFO(...)
Definition EonLogger.h:249
eonc::Matter::distance
double distance(long index1, long index2) const
Definition Matter.cpp:448
eonc::Matter::getAtomicNr
long getAtomicNr(long int atom) const
Definition Matter.cpp:495

◆ translationRemove() [1/2]

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

Definition at line 240 of file GeometryAnalysis.cpp.

240 {
241 AtomMatrix r1 = m1.getPositions();
242 AtomMatrix r2 = r2_passed;
243
244 // net displacement
245 Eigen::VectorXd disp(3);
246 AtomMatrix r12 = m1.pbc(r2 - r1);
247
248 disp[0] = r12.col(0).sum();
249 disp[1] = r12.col(1).sum();
250 disp[2] = r12.col(2).sum();
251 disp /= r1.rows();
252
253 for (int i = 0; i < r1.rows(); i++) {
254 r1(i, 0) += disp[0];
255 r1(i, 1) += disp[1];
256 r1(i, 2) += disp[2];
257 }
258
259 m1.setPositions(r1);
260 return;
261}
eonc::Matter::setPositions
void setPositions(const AtomMatrix &pos)
Definition Matter.cpp:350
eonc::Matter::pbc
AtomMatrix pbc(const AtomMatrix &diff) const
Definition Matter.cpp:810

◆ translationRemove() [2/2]

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

Definition at line 263 of file GeometryAnalysis.cpp.

263 {
264 AtomMatrix r2 = m2.getPositions();
265 translationRemove(m1, r2);
266 return;
267}
eonc::geometry::translationRemove
void translationRemove(Matter &m1, const AtomMatrix r1)
Definition GeometryAnalysis.cpp:240