20 Eigen::MatrixXd pointSet(n, 3);
25 for (
int k = 0; k < 3; k++) {
37 for (
int i = 0; i < n; i++) {
38 pointSet(i, dims[0]) = cos((2.0 *
gen::pi * i) / n);
39 pointSet(i, dims[1]) = sin((2.0 *
gen::pi * i) / n);
41 pointSet(i, axialDim) = 0.0;
57 int nop = ringSize * 2;
58 Eigen::MatrixXd pointSet(nop, 3);
83 for (
int i = 0; i < ringSize; i++) {
85 iBasalTwo = i + ringSize;
87 for (
int k = 0; k < 3; k++) {
90 pointSet(iBasalOne, k) = 0.0;
91 pointSet(iBasalTwo, k) = avgHeight;
94 pointSet(iBasalOne, k) = avgRadius * refPoints(i, k);
95 pointSet(iBasalTwo, k) = avgRadius * refPoints(i, k);
111 double avgRadius = 0.0;
114 int ringSize = basal1.
size();
123 for (
int i = 0; i < ringSize; i++) {
125 centroid1[0] += yCloud->pts[basal1[i]].x;
126 centroid1[1] += yCloud->pts[basal1[i]].y;
127 centroid1[2] += yCloud->pts[basal1[i]].z;
130 centroid2[0] += yCloud->pts[basal2[i]].x;
131 centroid2[1] += yCloud->pts[basal2[i]].y;
132 centroid2[2] += yCloud->pts[basal2[i]].z;
136 for (
int k = 0; k < 3; k++) {
137 centroid1[k] /= ringSize;
138 centroid2[k] /= ringSize;
143 for (
int i = 0; i < ringSize; i++) {
171 double avgHeight = 0.0;
174 int ringSize = basal1.
size();
179 for (
int i = 0; i < ringSize; i++) {
209 int ringSize = basal1.
size();
213 bool isClock, isAntiClock;
216 int currentIatom, currentJatom;
222 for (
int l = 0; l < ringSize; l++) {
227 for (
int m = 0; m < ringSize; m++) {
231 auto it =
std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
234 if (it != nList[l_k].end()) {
249 std::cerr <<
"Something is wrong with your deformed prism.\n";
259 int tempJfor, tempJback;
261 tempJfor = jatom + 1;
262 tempJback = jatom - 1;
264 if (jatom == ringSize - 1) {
266 tempJback = ringSize - 2;
270 tempJback = ringSize - 1;
273 int forwardJ = basal2[tempJfor];
274 int backwardJ = basal2[tempJback];
275 int currentI = basal1[iatom];
282 if (distClock < distAntiClock) {
286 if (distAntiClock < distClock) {
290 if (isClock ==
false && isAntiClock ==
false) {
297 for (
int i = 0; i < ringSize; i++) {
298 currentIatom = iatom + i;
299 if (currentIatom >= ringSize) {
300 currentIatom -= ringSize;
305 currentJatom = jatom + i;
306 if (currentJatom >= ringSize) {
307 currentJatom -= ringSize;
311 currentJatom = jatom - i;
312 if (currentJatom < 0) {
313 currentJatom += ringSize;
318 (*outBasal1).push_back(basal1[currentIatom]);
319 (*outBasal2).push_back(basal2[currentJatom]);
338 int ringSize = basal1.
size();
342 bool isClock, isAntiClock;
345 int currentIatom, currentJatom;
349 double infHugeNumber = 100000;
350 double leastDist = infHugeNumber;
358 for (
int m = 0; m < ringSize; m++) {
365 if (leastDist > dist) {
373 if (leastDist < infHugeNumber) {
379 std::cerr <<
"Something is wrong with your deformed prism.\n";
389 int tempJfor, tempJback;
391 tempJfor = jatom + 1;
392 tempJback = jatom - 1;
394 if (jatom == ringSize - 1) {
396 tempJback = ringSize - 2;
400 tempJback = ringSize - 1;
403 int forwardJ = basal2[tempJfor];
404 int backwardJ = basal2[tempJback];
405 int currentI = basal1[iatom];
412 if (distClock < distAntiClock) {
416 if (distAntiClock < distClock) {
420 if (isClock ==
false && isAntiClock ==
false) {
427 for (
int i = 0; i < ringSize; i++) {
428 currentIatom = iatom + i;
429 if (currentIatom >= ringSize) {
430 currentIatom -= ringSize;
435 currentJatom = jatom + i;
436 if (currentJatom >= ringSize) {
437 currentJatom -= ringSize;
441 currentJatom = jatom - i;
442 if (currentJatom < 0) {
443 currentJatom += ringSize;
448 (*outBasal1).push_back(basal1[currentIatom]);
449 (*outBasal2).push_back(basal2[currentJatom]);
466 int ringSize = basalRing.
size();
467 int nop = basalRing.
size();
468 Eigen::MatrixXd pointSet(nop, dim);
474 if (startingIndex >= ringSize || startingIndex < 0) {
479 for (
int i = 0; i < nop; i++) {
482 currentPosition = startingIndex + i;
484 if (currentPosition >= ringSize) {
485 currentPosition -= ringSize;
490 index = basalRing[currentPosition];
491 pointSet(i, 0) = yCloud->pts[index].x;
492 pointSet(i, 1) = yCloud->pts[index].y;
493 pointSet(i, 2) = yCloud->pts[index].z;
510 int ringSize = basal1.
size();
511 int nop = ringSize * 2;
512 Eigen::MatrixXd pointSet(nop, dim);
516 int iatomIndex, jatomIndex;
522 if (startingIndex >= ringSize || startingIndex < 0) {
527 for (
int i = 0; i < ringSize; i++) {
530 currentPosition = startingIndex + i;
532 if (currentPosition >= ringSize) {
533 currentPosition -= ringSize;
539 basal1[currentPosition];
541 pointSet(iatom, 0) = yCloud->pts[iatomIndex].x;
542 pointSet(iatom, 1) = yCloud->pts[iatomIndex].y;
543 pointSet(iatom, 2) = yCloud->pts[iatomIndex].z;
547 basal2[currentPosition];
548 jatom = i + ringSize;
549 pointSet(jatom, 0) = yCloud->pts[jatomIndex].x;
550 pointSet(jatom, 1) = yCloud->pts[jatomIndex].y;
551 pointSet(jatom, 2) = yCloud->pts[jatomIndex].z;
566 Eigen::MatrixXd pointSet(12, 3);
575 int n = setCloud.
nop;
576 Eigen::MatrixXd pointSet(n, 3);
579 for (
int i = 0; i < n; i++) {
580 pointSet(i, 0) = setCloud.
pts[i].x;
581 pointSet(i, 1) = setCloud.
pts[i].y;
582 pointSet(i, 2) = setCloud.
pts[i].z;
613 int ringSize = basal1.
size();
614 bool neighOne, neighTwo;
619 int nextBasal1Element;
620 int nextBasal2Element;
621 bool isReversedOrder;
624 (*matchedBasal1).resize(ringSize);
625 (*matchedBasal2).resize(ringSize);
630 for (
int i = 0; i < ringSize; i++) {
633 auto it =
std::find(nList[l1].begin() + 1, nList[l1].end(), iatom);
635 if (it != nList[l1].end()) {
636 neighbourFound =
true;
640 nextBasal1Element = basal1[2];
645 auto it1 =
std::find(nList[l2].begin() + 1, nList[l2].end(), iatom);
647 if (it1 != nList[l2].end()) {
648 neighbourFound =
true;
651 nextBasal1Element = basal1[3];
661 if (!neighbourFound) {
672 index = m_kIndex + 2;
674 if (index >= ringSize) {
677 nextBasal2Element = basal2[index];
679 auto it =
std::find(nList[nextBasal1Element].begin() + 1,
680 nList[nextBasal1Element].end(), nextBasal2Element);
682 if (it != nList[nextBasal1Element].end()) {
684 for (
int i = 0; i < ringSize; i++) {
685 index = m_kIndex + i;
687 if (index >= ringSize) {
690 (*matchedBasal2)[i] = basal2[index];
695 index = m_kIndex - 2;
700 nextBasal2Element = basal2[index];
702 auto it =
std::find(nList[nextBasal1Element].begin() + 1,
703 nList[nextBasal1Element].end(), nextBasal2Element);
705 if (it != nList[nextBasal1Element].end()) {
707 for (
int i = 0; i < ringSize; i++) {
708 index = m_kIndex - i;
713 (*matchedBasal2)[i] = basal2[index];
726 (*matchedBasal1) = basal1;
739 Eigen::MatrixXd pointSet(12, 3);
740 int iatomIndex, jatomIndex;
750 if (startingIndex > 5 || startingIndex < 0) {
755 if (startingIndex > 0) {
756 for (
int k = 0; k < 6; k++) {
757 iPnt = k + startingIndex;
773 iatomOne = newBasal1[0];
774 pointSet(0, 0) = yCloud->pts[iatomOne].x;
775 pointSet(0, 1) = yCloud->pts[iatomOne].y;
776 pointSet(0, 2) = yCloud->pts[iatomOne].z;
778 jatomIndex = newBasal2[0];
783 pointSet(6, 0) = yCloud->pts[iatomOne].x + dr[0];
784 pointSet(6, 1) = yCloud->pts[iatomOne].y + dr[1];
785 pointSet(6, 2) = yCloud->pts[iatomOne].z + dr[2];
788 for (
int i = 1; i < 6; i++) {
790 iatomIndex = newBasal1[i];
791 jatomIndex = newBasal2[i];
796 pointSet(i, 0) = yCloud->pts[iatomOne].x + dr[0];
797 pointSet(i, 1) = yCloud->pts[iatomOne].y + dr[1];
798 pointSet(i, 2) = yCloud->pts[iatomOne].z + dr[2];
803 pointSet(i + 6, 0) = yCloud->pts[iatomOne].x + dr[0];
804 pointSet(i + 6, 1) = yCloud->pts[iatomOne].y + dr[1];
805 pointSet(i + 6, 2) = yCloud->pts[iatomOne].z + dr[2];
854 int nextI, prevI, nextJ, prevJ;
855 int jatomIndex, atomIndex;
859 iring = cageList[index]
863 for (
int i = 0; i < ringSize; i++) {
871 for (
int i = 0; i < ringSize; i++) {
876 if (i == ringSize - 1) {
878 prevI = ddcOrder[i - 1];
881 nextI = ddcOrder[i + 1];
884 nextI = ddcOrder[i + 1];
885 prevI = ddcOrder[i - 1];
891 for (
int j = 1; j <= ringSize; j++) {
893 jring = cageList[index].rings[j];
896 auto it =
std::find(rings[jring].begin(), rings[jring].end(), iatom);
898 if (it != rings[jring].end()) {
906 atomIndex = jatomIndex + 1;
908 if (atomIndex >= ringSize) {
909 atomIndex -= ringSize;
911 nextJ = rings[jring][atomIndex];
913 atomIndex = jatomIndex - 1;
915 atomIndex += ringSize;
917 prevJ = rings[jring][atomIndex];
922 if (prevJ != nextI && prevJ != prevI) {
930 atomIndex = jatomIndex - 2;
933 atomIndex += ringSize;
935 apex1 = rings[jring][atomIndex];
944 atomIndex = jatomIndex - 2;
947 atomIndex += ringSize;
949 apex2 = rings[jring][atomIndex];
958 if (nextJ != nextI && nextJ != prevI) {
966 atomIndex = jatomIndex + 2;
968 if (atomIndex >= ringSize) {
969 atomIndex -= ringSize;
971 apex1 = rings[jring][atomIndex];
980 atomIndex = jatomIndex + 2;
982 if (atomIndex >= ringSize) {
983 atomIndex -= ringSize;
985 apex2 = rings[jring][atomIndex];
1000 for (
int i = 0; i < peripheral1.
size(); i++) {
1008 for (
int i = 0; i < peripheral2.
size(); i++) {
1043 Eigen::MatrixXd pointSet(nop, 3);
1044 int peripheralStartingIndex;
1050 int iatomIndex, jatomIndex;
1053 if (startingIndex == 0) {
1054 wrappedDDC = ddcOrder;
1061 for (
int i = 0; i < ringSize; i++) {
1062 currentIndex = startingIndex + i;
1064 if (currentIndex >= ringSize) {
1065 currentIndex -= ringSize;
1068 wrappedDDC[i] = ddcOrder[currentIndex];
1074 if (startingIndex <= 1) {
1075 peripheralStartingIndex = 0;
1076 }
else if (startingIndex > 1 && startingIndex <= 3) {
1077 peripheralStartingIndex = 1;
1079 peripheralStartingIndex = 2;
1083 for (
int i = 6; i < 9; i++) {
1084 currentIndex = i + peripheralStartingIndex;
1086 if (currentIndex <= 9) {
1089 wrappedDDC[i] = ddcOrder[currentIndex];
1093 for (
int i = 10; i < 13; i++) {
1094 currentIndex = i + peripheralStartingIndex;
1096 if (currentIndex <= 13) {
1099 wrappedDDC[i] = ddcOrder[currentIndex];
1103 wrappedDDC[9] = ddcOrder[9];
1104 wrappedDDC[13] = ddcOrder[13];
1109 iatomIndex = wrappedDDC[0];
1110 pointSet(0, 0) = yCloud->pts[iatomIndex].x;
1111 pointSet(0, 1) = yCloud->pts[iatomIndex].y;
1112 pointSet(0, 2) = yCloud->pts[iatomIndex].z;
1114 for (
int i = 1; i < 14; i++) {
1116 jatomIndex = wrappedDDC[i];
1121 pointSet(i, 0) = yCloud->pts[iatomIndex].x + dr[0];
1122 pointSet(i, 1) = yCloud->pts[iatomIndex].y + dr[1];
1123 pointSet(i, 2) = yCloud->pts[iatomIndex].z + dr[2];
std::array< double, 3 > relDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
double periodicDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Inline generic function for obtaining the unwrapped periodic distance between two particles,...
double unWrappedDistFromPoint(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, std::vector< double > singlePoint)
double getAverageWithoutOutliers(std::vector< double > inpVec)
Get the average, after excluding the outliers, using quartiles.
int nop
Current frame number.
@ HCbasal
The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an H...
int readXYZ(std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Function for reading in atom coordinates from an XYZ file.
Eigen::MatrixXd changeDiaCageOrder(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > ddcOrder, int startingIndex=0)
double getRadiusFromRings(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2)
Eigen::MatrixXd getPointSetCage(ring::strucType type)
Eigen::MatrixXd fillPointSetPrismRing(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basalRing, int startingIndex)
Eigen::MatrixXd getPointSetRefRing(int n, int axialDim)
Eigen::MatrixXd createPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, int ringSize, std::vector< int > basal1, std::vector< int > basal2)
std::vector< int > relOrderDDC(int index, std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList)
Matches the order of the basal rings of an DDC or a potential HC.
int relOrderHC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, std::vector< std::vector< int >> nList, std::vector< int > *matchedBasal1, std::vector< int > *matchedBasal2)
Matches the order of the basal rings of an HC or a potential HC.
Eigen::MatrixXd changeHexCageOrder(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex=0)
int relOrderPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, std::vector< std::vector< int >> nList, std::vector< int > *outBasal1, std::vector< int > *outBasal2)
double getAvgHeightPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2)
Eigen::MatrixXd fillPointSetPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex)
This contains a collection of points; contains information for a particular frame.
This contains per-particle information.