pntToPnt Namespace Reference

Functions

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)
 
double getRadiusFromRings (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2)
 
double getAvgHeightPrismBlock (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2)
 
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)
 
int relOrderPrismBlock (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, std::vector< int > *outBasal1, std::vector< int > *outBasal2)
 
Eigen::MatrixXd fillPointSetPrismRing (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basalRing, int startingIndex)
 
Eigen::MatrixXd fillPointSetPrismBlock (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex)
 
Eigen::MatrixXd getPointSetCage (ring::strucType type)
 
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. More...
 
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. More...
 
Eigen::MatrixXd changeHexCageOrder (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex=0)
 
Eigen::MatrixXd changeDiaCageOrder (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > ddcOrder, int startingIndex=0)
 

Function Documentation

◆ changeDiaCageOrder()

Eigen::MatrixXd pntToPnt::changeDiaCageOrder ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  ddcOrder,
int  startingIndex = 0 
)

Fills up an eigen matrix point set using information of the equatorial ring and peripheral rings, embedded in a vector, already filled in relOrderDDC.

Fills up an eigen matrix point set using information of the equatorial ring and peripheral rings, embedded in a vector, already filled in relOrderDDC. The order is such that the next three atoms in the vector are the first nearest neighbours in the peripheral rings of one triplet (peripheral1), followed by apex1, and three nearest neighbours of the other triplet, followed by apex2. Basically: ((6 equatorial ring atoms), (3 nearest neighbours of [l1,l3,l5] in the peripheral rings), (1 second shell neighbour of [l1,l3,l5]), (first nearest neighbours of [l2,l4,l6]), (second nearest neighbour of [l2,l4,l6]) )

Thus, when you want to change the order of the DDC, this is how the order should be wrapped:

  1. The first 6 atoms should be wrapped around (i.e. the sixth atom should become the first and so on)
  2. The next 3 atoms should be wrapped but only when moved by multiples of 2.
  3. The 9th element should not be changed.
  4. The next three elements should be wrapped around (multiples of 2), since alternate elements of the equatorial ring are bonded.

Definition at line 1038 of file pntCorrespondence.cpp.

1040  {
1041  int nop = 14; // Number of elements in the DDC
1042  int ringSize = 6; // Six nodes in the rings
1043  Eigen::MatrixXd pointSet(nop, 3); // Output point set
1044  int peripheralStartingIndex; // Index using which the elements of peripheral1
1045  // and peripheral2 will be wrapped
1046  std::vector<int> wrappedDDC; // Changed order of the DDC: should be the same
1047  // size as the original ddcOrder vector
1048  int currentIndex; // Current index
1049  // Variables for filling the point set
1050  int iatomIndex, jatomIndex;
1051  std::array<double, 3> dr; // Components of the distance
1052 
1053  if (startingIndex == 0) {
1054  wrappedDDC = ddcOrder;
1055  } // if the order does not have to be changed
1056  else {
1057  wrappedDDC.resize(nop); // Should be the same size as ddcOrder
1058  // ------------------
1059  // EQUATORIAL RING
1060  // Change the order of the equatorial ring
1061  for (int i = 0; i < ringSize; i++) {
1062  currentIndex = startingIndex + i;
1063  // Wrap-around
1064  if (currentIndex >= ringSize) {
1065  currentIndex -= ringSize;
1066  } // end of wrap-around
1067  // Update the equatorial ring
1068  wrappedDDC[i] = ddcOrder[currentIndex];
1069  } // end of wrapping the equatorial ring
1070  // ------------------
1071  // PERIPHERAL RINGS
1072  //
1073  // The index of the peripheral ring starting indices
1074  if (startingIndex <= 1) {
1075  peripheralStartingIndex = 0;
1076  } else if (startingIndex > 1 && startingIndex <= 3) {
1077  peripheralStartingIndex = 1;
1078  } else {
1079  peripheralStartingIndex = 2;
1080  }
1081  //
1082  // Update the portions of the wrappedDDC vector
1083  for (int i = 6; i < 9; i++) {
1084  currentIndex = i + peripheralStartingIndex;
1085  // wrap-around
1086  if (currentIndex <= 9) {
1087  currentIndex -= 3;
1088  } // end of wrap-around
1089  wrappedDDC[i] = ddcOrder[currentIndex];
1090  } // peripheral1
1091  //
1092  // Update the peripheral2 portions of the wrappedDDC vector
1093  for (int i = 10; i < 13; i++) {
1094  currentIndex = i + peripheralStartingIndex;
1095  // wrap-around
1096  if (currentIndex <= 13) {
1097  currentIndex -= 3;
1098  } // end of wrap-around
1099  wrappedDDC[i] = ddcOrder[currentIndex];
1100  } // peripheral2
1101  // ------------------
1102  // SECOND-NEAREST NEIGHBOURS
1103  wrappedDDC[9] = ddcOrder[9]; // apex1
1104  wrappedDDC[13] = ddcOrder[13]; // apex2
1105  // ------------------
1106  } // the order of the DDC has to be changed
1107 
1108  // FILL UP THE EIGEN MATRIX
1109  iatomIndex = wrappedDDC[0]; // first point
1110  pointSet(0, 0) = yCloud->pts[iatomIndex].x;
1111  pointSet(0, 1) = yCloud->pts[iatomIndex].y;
1112  pointSet(0, 2) = yCloud->pts[iatomIndex].z;
1113  // Loop through the rest of the equatorial ring points
1114  for (int i = 1; i < 14; i++) {
1115  //
1116  jatomIndex = wrappedDDC[i]; // Atom index to be filled
1117  // Get the distance from iatomIndex
1118  dr = gen::relDist(yCloud, iatomIndex, jatomIndex);
1119 
1120  // basal2
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];
1124  } // end of loop
1125 
1126  return pointSet;
1127 }
std::array< double, 3 > relDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Definition: generic.hpp:196
std::vector< S > pts
Definition: mol_sys.hpp:167
T resize(T... args)

◆ changeHexCageOrder()

Eigen::MatrixXd pntToPnt::changeHexCageOrder ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  basal1,
std::vector< int >  basal2,
int  startingIndex = 0 
)

Fills up an eigen matrix point set using the basal rings basal1 and basal2, changing the order of the point set by filling up from the startingIndex (starting from 0 to 5)

Definition at line 736 of file pntCorrespondence.cpp.

738  {
739  Eigen::MatrixXd pointSet(12, 3);
740  int iatomIndex, jatomIndex; // Current atom index in yCloud, according to
741  // basal1 and basal2 respectively
742  int iPnt; // Current index in the Eigen matrix pointSet
743  int cageSize = 12; // Number of points in the cage
744  std::vector<int> newBasal1, newBasal2;
745  std::array<double, 3> dr; // Components of the distance
746  int iatomOne; // Index of the first atom
747 
748  // Checks and balances
749  //
750  if (startingIndex > 5 || startingIndex < 0) {
751  startingIndex = 0;
752  } // no invalid starting index
753 
754  // Change the order
755  if (startingIndex > 0) {
756  for (int k = 0; k < 6; k++) {
757  iPnt = k + startingIndex;
758  if (iPnt >= 6) {
759  iPnt -= 6;
760  } // wrap-around
761  newBasal1.push_back(basal1[iPnt]);
762  newBasal2.push_back(basal2[iPnt]);
763  } // change the order
764  } // end of filling for startingIndex>0
765  else {
766  newBasal1 = basal1;
767  newBasal2 = basal2;
768  } // end of filling up the reordered basal rings
769 
770  //
771  // FIRST POINT
772  // basal1
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;
777  // basal2
778  jatomIndex = newBasal2[0];
779  // Get the distance from basal1
780  dr = gen::relDist(yCloud, iatomOne, jatomIndex);
781 
782  // basal2
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];
786  //
787  // Loop through the rest of the points
788  for (int i = 1; i < 6; i++) {
789  // basal1
790  iatomIndex = newBasal1[i]; // Atom index to be filled for basal1
791  jatomIndex = newBasal2[i]; // Atom index to be filled for basal2
792  //
793  // Get the distance from the first atom
794  dr = gen::relDist(yCloud, iatomOne, iatomIndex);
795  //
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];
799  //
800  // Get the distance from the first atom
801  dr = gen::relDist(yCloud, iatomOne, jatomIndex);
802  // basal2
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];
806  //
807  } // end of loop
808 
809  return pointSet;
810 }
T push_back(T... args)

◆ createPrismBlock()

Eigen::MatrixXd pntToPnt::createPrismBlock ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
const Eigen::MatrixXd &  refPoints,
int  ringSize,
std::vector< int >  basal1,
std::vector< int >  basal2 
)

Creates an eigen matrix for the points of a prism block, constructed from the points of a perfect polygon of radius 1, given the basal rings and axial dimension

Definition at line 52 of file pntCorrespondence.cpp.

55  {
56  //
57  int nop = ringSize * 2; // Number of particles in the prism block
58  Eigen::MatrixXd pointSet(nop, 3); // Output point set for a regular prism
59  int axialDim; // Has a value of 0, 1 or 2 for x, y or z
60  double avgRadius; // Average radius within which the basal ring points lie
61  double avgHeight; // Get the average height of the prism
62  int iBasalOne,
63  iBasalTwo; // Indices for the basal1 and basal2 points in the point set
64  // --------------------------------------
65  // Get the axial dimension
66  // The axial dimension will have the largest box length
67  // Index -> axial dimension
68  // 0 -> x dim
69  // 1 -> y dim
70  // 2 -> z dim
71  axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
72  yCloud->box.begin();
73  // --------------------------------------
74  // Get the average radius
75  avgRadius = pntToPnt::getRadiusFromRings(yCloud, basal1, basal2);
76  // --------------------------------------
77  // Get the average height of the prism
78  avgHeight = pntToPnt::getAvgHeightPrismBlock(yCloud, basal1, basal2);
79  // --------------------------------------
80  // Fill in the matrix of the point set
81  //
82  // Fill basal1 first, and then basal2
83  for (int i = 0; i < ringSize; i++) {
84  iBasalOne = i; // index in point set for basal1 point
85  iBasalTwo = i + ringSize; // index in point set for basal2 point
86  // Fill up the dimensions
87  for (int k = 0; k < 3; k++) {
88  // For the axial dimension
89  if (k == axialDim) {
90  pointSet(iBasalOne, k) = 0.0; // basal1
91  pointSet(iBasalTwo, k) = avgHeight; // basal2
92  } // end of filling up the axial dimension
93  else {
94  pointSet(iBasalOne, k) = avgRadius * refPoints(i, k); // basal1
95  pointSet(iBasalTwo, k) = avgRadius * refPoints(i, k); // basal2
96  } // fill up the non-axial dimension coordinate dimension
97  } // Fill up all the dimensions
98  } // end of filling in the point set
99  // --------------------------------------
100  return pointSet;
101 } // end of function
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:170
T max_element(T... args)
double getRadiusFromRings(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2)
double getAvgHeightPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2)

◆ fillPointSetPrismBlock()

Eigen::MatrixXd pntToPnt::fillPointSetPrismBlock ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  basal1,
std::vector< int >  basal2,
int  startingIndex 
)

Fill up an Eigen matrix of a prism block (two basal rings) from input vectors for the basal rings

Fill up an Eigen Matrix for a prism block from input vectors of atom indices of the basal rings

Definition at line 504 of file pntCorrespondence.cpp.

506  {
507  //
508  //
509  int dim = 3; // Number of dimensions (3)
510  int ringSize = basal1.size(); // Number of nodes in the ring
511  int nop = ringSize * 2; // Number of particles in prism block
512  Eigen::MatrixXd pointSet(nop, dim); // Output set of 3D points
513  // Indices for the first and second basal rings being filled
514  int currentPosition; // Current position in the basal ring vectors
515  int iatom, jatom; // Current index in the point set
516  int iatomIndex, jatomIndex; // Index in yCloud corresponding to the basal1
517  // and basal2 points
518 
519  // Fill up basal1 points, followed by basal2 points
520 
521  // Check
522  if (startingIndex >= ringSize || startingIndex < 0) {
523  startingIndex = 0;
524  } // end of check
525 
526  // Beginning from the starting index, get the points in the basal ring
527  for (int i = 0; i < ringSize; i++) {
528  //
529  // Getting currentIndex
530  currentPosition = startingIndex + i;
531  // Wrapping around for the ring
532  if (currentPosition >= ringSize) {
533  currentPosition -= ringSize;
534  } // end of wrap-around
535  //
536  // -------------------
537  // Basal1 ring points
538  iatomIndex =
539  basal1[currentPosition]; // Index of the current point in basal1
540  iatom = i; // index in the point set being filled
541  pointSet(iatom, 0) = yCloud->pts[iatomIndex].x; // x coord
542  pointSet(iatom, 1) = yCloud->pts[iatomIndex].y; // y coord
543  pointSet(iatom, 2) = yCloud->pts[iatomIndex].z; // z coord
544  // -------------------
545  // Basal2 ring points
546  jatomIndex =
547  basal2[currentPosition]; // Index of the current point in basal2
548  jatom = i + ringSize; // index in the point set being filled
549  pointSet(jatom, 0) = yCloud->pts[jatomIndex].x; // x coord
550  pointSet(jatom, 1) = yCloud->pts[jatomIndex].y; // y coord
551  pointSet(jatom, 2) = yCloud->pts[jatomIndex].z; // z coord
552  } // end of point filling from the relative ordering vector of vectors
553 
554  // Return the set of points
555  return pointSet;
556 } // end of function
T size(T... args)

◆ fillPointSetPrismRing()

Eigen::MatrixXd pntToPnt::fillPointSetPrismRing ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  basalRing,
int  startingIndex 
)

Fill up an Eigen Matrix of a prism basal ring from an input vector of atom indices

Fill up an Eigen Matrix for a prism basal ring from an input vector of atom indices

Definition at line 460 of file pntCorrespondence.cpp.

462  {
463  //
464  //
465  int dim = 3; // Number of dimensions (3)
466  int ringSize = basalRing.size(); // Number of nodes in the ring
467  int nop = basalRing.size(); // Number of particles in the basal ring
468  Eigen::MatrixXd pointSet(nop, dim); // Output set of 3D points
469  // Indices for the first and second basal rings being filled
470  int currentPosition; // Current index in basalRing
471  int index; // Index in yCloud
472 
473  // Check
474  if (startingIndex >= ringSize || startingIndex < 0) {
475  startingIndex = 0;
476  } // end of check
477 
478  // Beginning from the starting index, get the points in the basal ring
479  for (int i = 0; i < nop; i++) {
480  //
481  // Getting currentIndex
482  currentPosition = startingIndex + i;
483  // Wrapping around for the ring
484  if (currentPosition >= ringSize) {
485  currentPosition -= ringSize;
486  } // end of wrap-around
487  //
488  // -------------------
489  // Basal ring points
490  index = basalRing[currentPosition]; // Index of the current point
491  pointSet(i, 0) = yCloud->pts[index].x; // x coord
492  pointSet(i, 1) = yCloud->pts[index].y; // y coord
493  pointSet(i, 2) = yCloud->pts[index].z; // z coord
494  } // end of point filling from the relative ordering vector of vectors
495 
496  // Return the set of points
497  return pointSet;
498 } // end of function

◆ getAvgHeightPrismBlock()

double pntToPnt::getAvgHeightPrismBlock ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  basal1,
std::vector< int >  basal2 
)

Calculate the average height of the prism block, calculated using the basal rings of the prism and the axial dimension

Definition at line 167 of file pntCorrespondence.cpp.

169  {
170  //
171  double avgHeight = 0.0;
172  double r_ij; // Distance between a point in basal1 and basal2
173  std::vector<double> dist; // Distances
174  int ringSize = basal1.size(); // Number of nodes in the basal rings
175  int iatom; // Atom index in yCloud for the current point in basal1
176  int jatom; // Atom index in yCloud for the current point in basal2
177  // -----------------------
178  // Calculate the distances of each point from the respective centroid
179  for (int i = 0; i < ringSize; i++) {
180  // Get the current point in basal1 and basal2
181  iatom = basal1[i];
182  jatom = basal2[i];
183  // Get the distance between a point in basal1 and the corresponding point in
184  // basal2
185  r_ij = gen::periodicDist(yCloud, iatom, jatom);
186  // Update the distances
187  dist.push_back(r_ij);
188  } // Loop through every point in basal1 and basal2
189  // -----------------------
190  // Calculate the average, excluding the outliers
191  avgHeight = gen::getAverageWithoutOutliers(dist);
192  // -----------------------
193  return avgHeight;
194 } // end of function
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,...
Definition: generic.hpp:104
double getAverageWithoutOutliers(std::vector< double > inpVec)
Get the average, after excluding the outliers, using quartiles.
Definition: generic.cpp:164

◆ getPointSetCage()

Eigen::MatrixXd pntToPnt::getPointSetCage ( ring::strucType  type)

Fills up an eigen matrix point set a reference ring, which is a regular n-gonal polygon, constructed with radius 1 by default; where n is the number of nodes in the ring

REFERENCE POINT SETS FOR BULK ICES Fills up an eigen matrix point set for a reference cage, saved in the templates folder relative to the top-level directory NOT NEEDED MAYBE

Definition at line 564 of file pntCorrespondence.cpp.

564  {
565  //
566  Eigen::MatrixXd pointSet(12, 3); // Reference point set (Eigen matrix)
568  setCloud; // PointCloud for holding the reference point values
569 
570  if (type == ring::HCbasal) {
571  // Read in the XYZ file
572  std::string fileName = "templates/hc.xyz";
573  //
574  sinp::readXYZ(fileName, &setCloud);
575  int n = setCloud.nop; // Number of points in the reference point set
576  Eigen::MatrixXd pointSet(n, 3); // Output point set for a regular polygon
577 
578  // Loop through every particle
579  for (int i = 0; i < n; i++) {
580  pointSet(i, 0) = setCloud.pts[i].x; // x
581  pointSet(i, 1) = setCloud.pts[i].y; // y
582  pointSet(i, 2) = setCloud.pts[i].z; // z
583  } // end of looping through every particle
584 
585  // Return the filled point set
586  return pointSet;
587  } // end of reference point set for HCs
588  // else {
589  // // ddc
590  // Eigen::MatrixXd pointSet(14, 3); // Output point set for a regular
591  // polygon return pointSet;
592  // } // DDCs
593 
594  // Eigen::MatrixXd pointSet(14, 3); // Output point set for a regular polygon
595  // // Never happens
596  return pointSet;
597 
598 } // end of function
int nop
Current frame number.
Definition: mol_sys.hpp:169
@ HCbasal
The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an H...
Definition: ring.hpp:115
int readXYZ(std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Function for reading in atom coordinates from an XYZ file.
Definition: seams_input.cpp:69
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166

◆ getPointSetRefRing()

Eigen::MatrixXd pntToPnt::getPointSetRefRing ( int  n,
int  axialDim 
)

Fills up an eigen matrix point set a reference ring, which is a regular n-gonal polygon; where n is the number of nodes in the ring

Fills up an eigen matrix point set a reference ring, which is a regular n-gonal polygon, constructed with radius 1 by default; where n is the number of nodes in the ring

Definition at line 18 of file pntCorrespondence.cpp.

18  {
19  //
20  Eigen::MatrixXd pointSet(n, 3); // Output point set for a regular polygon
21  std::vector<int> dims; // Apart from the axial dimension
22 
23  // Set the axial dimension to 0, and fill up the vector of the other
24  // dimensions
25  for (int k = 0; k < 3; k++) {
26  if (k != axialDim) {
27  dims.push_back(k);
28  } // other dimensions
29  } // end of setting the axial dimension
30 
31  // To get the vertices of a regular polygon, use the following formula:
32  // x[i] = r * cos(2*pi*i/n)
33  // y[i] = r * sin(2*pi*i/n)
34  // where n is the number of points and the index i is 0<=i<=n
35 
36  // Loop through every particle
37  for (int i = 0; i < n; i++) {
38  pointSet(i, dims[0]) = cos((2.0 * gen::pi * i) / n); // x
39  pointSet(i, dims[1]) = sin((2.0 * gen::pi * i) / n); // y
40  // Set the axial dimension to zero
41  pointSet(i, axialDim) = 0.0;
42  } // end of loop through all the points
43 
44  return pointSet;
45 } // end of function
T cos(T... args)
const double pi
Definition: generic.hpp:50
T sin(T... args)

◆ getRadiusFromRings()

double pntToPnt::getRadiusFromRings ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  basal1,
std::vector< int >  basal2 
)

Calculate the average radial distance for the basal rings, calculated from the centroid of each basal ring

Definition at line 107 of file pntCorrespondence.cpp.

109  {
110  //
111  double avgRadius = 0.0;
112  std::vector<double> centroid1, centroid2;
113  std::vector<double> dist; // Distances
114  int ringSize = basal1.size(); // Number of nodes in the basal rings
115  double r1, r2; // Distance from the centroid
116  int iatom; // Atom index in yCloud for the current point in basal1
117  int jatom; // Atom index in yCloud for the current point in basal2
118  // -----------------------
119  // Calculate the centroid for basal1 and basal2
120  centroid1.resize(3);
121  centroid2.resize(3); // init
122  // Loop through basal1 and basal2
123  for (int i = 0; i < ringSize; i++) {
124  // For basal1
125  centroid1[0] += yCloud->pts[basal1[i]].x;
126  centroid1[1] += yCloud->pts[basal1[i]].y;
127  centroid1[2] += yCloud->pts[basal1[i]].z;
128  //
129  // For basal2
130  centroid2[0] += yCloud->pts[basal2[i]].x;
131  centroid2[1] += yCloud->pts[basal2[i]].y;
132  centroid2[2] += yCloud->pts[basal2[i]].z;
133  //
134  } // end of loop through basal1 and basal2
135  // Normalize by the number of nodes
136  for (int k = 0; k < 3; k++) {
137  centroid1[k] /= ringSize;
138  centroid2[k] /= ringSize;
139  } // end of dividing by the number
140  // Calculated the centroid for basal1 and basal2!
141  // -----------------------
142  // Calculate the distances of each point from the respective centroid
143  for (int i = 0; i < ringSize; i++) {
144  // Get the current point in basal1 and basal2
145  iatom = basal1[i];
146  jatom = basal2[i];
147  // Get the distance from the respective centroid
148  // basal1
149  r1 = gen::unWrappedDistFromPoint(yCloud, iatom, centroid1);
150  // basal2
151  r2 = gen::unWrappedDistFromPoint(yCloud, jatom, centroid2);
152  // Update the distances
153  dist.push_back(r1);
154  dist.push_back(r2);
155  } // Loop through every point in basal1 and basal2
156  // -----------------------
157  // Calculate the average, excluding the outliers
158  avgRadius = gen::getAverageWithoutOutliers(dist);
159  // -----------------------
160  return avgRadius;
161 } // end of function
double unWrappedDistFromPoint(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, std::vector< double > singlePoint)
Definition: generic.hpp:134

◆ relOrderDDC()

std::vector< int > pntToPnt::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.

Reorders the particles of a DDC into a vector, which contains the atom indices of the DDC particles. Here, index is the index in the cageList vector of structs, referring to a specific DDC. The order created is as follows:

  1. The first 6 atoms are the atoms of the equatorial ring, such that adjacent particles are bonded to each other (l1, l2, l3, l4, l5, l6).
  2. The particles (l1, l3, l5) and (l2, l4, l6) are nearest neighbours of atoms in two different sets of peripheral rings. Each set of three peripheral rings have a second-shell neighbour of the equatorial ring particles (called apex1 and apex2 respectively.)
  3. The order is such that the next three atoms in the vector are the first nearest neighbours in the peripheral rings of one triplet (peripheral1), followed by apex1, and three nearest neighbours of the other triplet, followed by apex2. Basically: ((6 equatorial ring atoms), (3 nearest neighbours of [l1,l3,l5] in the peripheral rings), (1 second shell neighbour of [l1,l3,l5]), (first nearest neighbours of [l2,l4,l6]), (second nearest neighbour of [l2,l4,l6]) )

Thus, when you want to change the order of the DDC, this is how the order should be wrapped:

  1. The first 6 atoms should be wrapped around (i.e. the sixth atom should become the first and so on)
  2. The next 3 atoms should be wrapped but only when moved by multiples of 2.
  3. The 9th element should not be changed.
  4. The next three elements should be wrapped around (multiples of 2), since alternate elements of the equatorial ring are bonded.

Definition at line 841 of file pntCorrespondence.cpp.

843  {
844  //
845  std::vector<int> ddcOrder; // Order of the particles in the DDC.
846  int nop = 14; // Number of elements in the DDC
847  int ringSize = 6; // Number of nodes in each ring
848  int iring, jring; // Ring indices of the DDC rings
849  int iatom; // Atom index in the equatorial ring
850  std::vector<int> peripheral1, peripheral2; // Vectors holding the neighbours
851  // of (l1,l3,l5) and (l2,l4,l6)
852  int apex1, apex2;
853  bool neighbourFound; // Neighbour found
854  int nextI, prevI, nextJ, prevJ; // elements around iatom and jatom
855  int jatomIndex, atomIndex;
856 
857  // Add the equatorial ring particles
858  //
859  iring = cageList[index]
860  .rings[0]; // Equatorial ring index in rings vector of vectors
861  //
862  // Loop through all the atoms of the equatorial ring
863  for (int i = 0; i < ringSize; i++) {
864  ddcOrder.push_back(rings[iring][i]);
865  } // end of adding equatorial ring particles
866  //
867  // ------------------------------
868  // Find the neighbouring atoms of (l1, l3 and l5) and (l2,l4,l6) by looping
869  // through the peripheral rings
870  //
871  for (int i = 0; i < ringSize; i++) {
872  // Init
873  iatom = ddcOrder[i]; // element of the equatorial ring
874  // Get the next and previous elements of iatom
875  // next element
876  if (i == ringSize - 1) {
877  nextI = ddcOrder[0];
878  prevI = ddcOrder[i - 1];
879  } // if last element
880  else if (i == 0) {
881  nextI = ddcOrder[i + 1];
882  prevI = ddcOrder[5]; // wrapped
883  } else {
884  nextI = ddcOrder[i + 1];
885  prevI = ddcOrder[i - 1];
886  }
887  // ------------------
888  // Search all the peripheral rings for iatom, get the nearest neighbours and
889  // apex1 and apex2, which are bonded to the nearest neighbours of the
890  // equatorial ring (thus they are second nearest neighbours)
891  for (int j = 1; j <= ringSize; j++) {
892  //
893  jring = cageList[index].rings[j]; // Peripheral ring
894  // Search for iatom
895  //
896  auto it = std::find(rings[jring].begin(), rings[jring].end(), iatom);
897  // If iatom was found in jring peripheral ring
898  if (it != rings[jring].end()) {
899  // ----------------
900  // Atom index for iatom found
901  jatomIndex = std::distance(rings[jring].begin(), it);
902  // ----------------
903  // Get the next and previous element
904  //
905  // Next atom
906  atomIndex = jatomIndex + 1;
907  // wrap-around
908  if (atomIndex >= ringSize) {
909  atomIndex -= ringSize;
910  } // end of wrap-around
911  nextJ = rings[jring][atomIndex];
912  // Previous atom
913  atomIndex = jatomIndex - 1;
914  if (atomIndex < 0) {
915  atomIndex += ringSize;
916  } // end of wrap-around
917  prevJ = rings[jring][atomIndex];
918  // ----------------
919  // Check to see if nextJ or prevJ are different from nextI and prevI
920  //
921  // Checking prevJ
922  if (prevJ != nextI && prevJ != prevI) {
923  // Add to the vector
924  // if even peripheral1
925  if (i % 2 == 0) {
926  peripheral1.push_back(prevJ);
927  // Get apex1 for i=0
928  if (i == 0) {
929  // Go back two elements from jatomIndex
930  atomIndex = jatomIndex - 2;
931  // wrap-around
932  if (atomIndex < 0) {
933  atomIndex += ringSize;
934  } // end of wrap-around
935  apex1 = rings[jring][atomIndex];
936  } // apex1 for i=0
937  } // peripheral1
938  // if odd peripheral2
939  else {
940  peripheral2.push_back(prevJ);
941  // Get apex2 for i=0
942  if (i == 1) {
943  // Go back two elements from jatomIndex
944  atomIndex = jatomIndex - 2;
945  // wrap-around
946  if (atomIndex < 0) {
947  atomIndex += ringSize;
948  } // end of wrap-around
949  apex2 = rings[jring][atomIndex];
950  } // apex2 for i=1
951  } // peripheral 2
952  // Get apex1 or apex2 for i=0 or i=1
953  break;
954  } // check if prevJ is the atom to be added, bonded to iatom
955  // ----------------------
956  //
957  // Checking nextJ
958  if (nextJ != nextI && nextJ != prevI) {
959  // Add to the vector
960  // if even peripheral1
961  if (i % 2 == 0) {
962  peripheral1.push_back(nextJ);
963  // Get apex1 for i=0
964  if (i == 0) {
965  // Go foward two elements from jatomIndex
966  atomIndex = jatomIndex + 2;
967  // wrap-around
968  if (atomIndex >= ringSize) {
969  atomIndex -= ringSize;
970  } // end of wrap-around
971  apex1 = rings[jring][atomIndex];
972  } // apex1 for i=0
973  } // peripheral1
974  // if odd peripheral2
975  else {
976  peripheral2.push_back(prevJ);
977  // Get apex2 for i=0
978  if (i == 1) {
979  // Go foward two elements from jatomIndex
980  atomIndex = jatomIndex + 2;
981  // wrap-around
982  if (atomIndex >= ringSize) {
983  atomIndex -= ringSize;
984  } // end of wrap-around
985  apex2 = rings[jring][atomIndex];
986  } // apex2 for i=1
987  } // peripheral 2
988  // Get apex1 or apex2 for i=0 or i=1
989  break;
990  } // check if prevJ is the atom to be added, bonded to iatom
991  // ----------------
992  } // end of check for iatom in jring
993  } // end of search for the peripheral rings
994  // ------------------
995  } // end of looping through the elements of the equatorial ring
996  // ------------------------------
997  // Update ddcOrder with peripheral1, followed by apex1, peripheral2 and apex2
998  //
999  // peripheral1
1000  for (int i = 0; i < peripheral1.size(); i++) {
1001  ddcOrder.push_back(peripheral1[i]);
1002  } // end of adding peripheral1
1003  //
1004  // apex1
1005  ddcOrder.push_back(apex1);
1006  //
1007  // peripheral2
1008  for (int i = 0; i < peripheral2.size(); i++) {
1009  ddcOrder.push_back(peripheral2[i]);
1010  } // end of adding peripheral2
1011  //
1012  // apex2
1013  ddcOrder.push_back(apex2);
1014  //
1015  return ddcOrder;
1016 } // end of the function
T begin(T... args)
T distance(T... args)
T end(T... args)
T find(T... args)

◆ relOrderHC()

int pntToPnt::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.

Matches the order of the basal rings of an HC or a potential HC. The goal is to find which elements are bonded to which and the relative order

Definition at line 605 of file pntCorrespondence.cpp.

609  {
610  //
611  int l1 = basal1[0]; // First element of basal1
612  int l2 = basal1[1]; // Second element of basal1
613  int ringSize = basal1.size(); // Number of nodes in the basal rings
614  bool neighOne, neighTwo; // Basal2 element is the neighbour of l1 or l2
615  bool neighbourFound; // neighbour found
616  int m_k; // Element of basal2 which is a neighbour of l1 or l2
617  int m_kIndex; // Index of basal2 which is a neighbour of l1 or l2
618  int iatom; // Current element of basal2 being searched for
619  int nextBasal1Element; // Next bonded element of basal1
620  int nextBasal2Element; // Next bonded element of basal2
621  bool isReversedOrder; // basal2 is reversed wrt basal1
622  int index;
623 
624  (*matchedBasal1).resize(ringSize);
625  (*matchedBasal2).resize(ringSize);
626 
627  // Search to see if l1 or l2 is a neighbour
628  // -------------------
629  // Searching for the neighbours of l1 or l2
630  for (int i = 0; i < ringSize; i++) {
631  iatom = basal2[i]; // Element of basal2
632  // Search for the current element in the neighbour list of l1
633  auto it = std::find(nList[l1].begin() + 1, nList[l1].end(), iatom);
634  // If iatom is the neighbour of l1
635  if (it != nList[l1].end()) {
636  neighbourFound = true;
637  neighOne = true;
638  m_k = iatom; // Element found
639  m_kIndex = i; // Index in basal2
640  nextBasal1Element = basal1[2];
641  break;
642  } // iatom is the neighbour of l1
643  // l2 neighbour check
644  else {
645  auto it1 = std::find(nList[l2].begin() + 1, nList[l2].end(), iatom);
646  // If iatom is the neighbour of l2
647  if (it1 != nList[l2].end()) {
648  neighbourFound = true;
649  neighOne = false;
650  neighTwo = true;
651  nextBasal1Element = basal1[3];
652  m_k = iatom; // Element found
653  m_kIndex = i; // Index in basal2
654  break;
655  } // iatom is the neighbour of l2
656  } // Check for the neighbour of l2
657  } // end of search through basal2 elements
658  //
659  // -------------------
660  // If a neighbour was not found, then there is some mistake
661  if (!neighbourFound) {
662  // std::cerr << "This is not an HC\n";
663  return 1;
664  }
665 
666  // ------------------------------
667  //
668  // This element should be the nearest neighbour of the corresponding element
669  // in basal2
670  //
671  // Testing the original order
672  index = m_kIndex + 2;
673  // wrap-around
674  if (index >= ringSize) {
675  index -= ringSize;
676  } // wrap-around
677  nextBasal2Element = basal2[index];
678  // Search for the next basal element
679  auto it = std::find(nList[nextBasal1Element].begin() + 1,
680  nList[nextBasal1Element].end(), nextBasal2Element);
681  // If this element is found, then the original order is correct
682  if (it != nList[nextBasal1Element].end()) {
683  // Fill up the temporary vector with basal2 elements
684  for (int i = 0; i < ringSize; i++) {
685  index = m_kIndex + i; // index in basal2
686  // wrap-around
687  if (index >= ringSize) {
688  index -= ringSize;
689  } // end of wrap-around
690  (*matchedBasal2)[i] = basal2[index]; // fill up values
691  } // end of filling up tempBasal2
692  } // the original order is correct
693  else {
694  //
695  index = m_kIndex - 2;
696  // wrap-around
697  if (index < 0) {
698  index += ringSize;
699  } // wrap-around
700  nextBasal2Element = basal2[index];
701  // Search for the next basal element
702  auto it = std::find(nList[nextBasal1Element].begin() + 1,
703  nList[nextBasal1Element].end(), nextBasal2Element);
704  // If this element is found, then the original order is correct
705  if (it != nList[nextBasal1Element].end()) {
706  // Fill up the temporary vector with basal2 elements
707  for (int i = 0; i < ringSize; i++) {
708  index = m_kIndex - i; // index in basal2
709  // wrap-around
710  if (index < 0) {
711  index += ringSize;
712  } // end of wrap-around
713  (*matchedBasal2)[i] = basal2[index]; // fill up values
714  } // end of filling up tempBasal2
715 
716  }
717  //
718  else {
719  // std::cerr << "not an HC\n";
720  return 1;
721  }
722  //
723  } // the reversed order is correct!
724  // ------------------------------
725  // Fill up basal1
726  (*matchedBasal1) = basal1;
727  return 0;
728 } // end of the function

◆ relOrderPrismBlock() [1/2]

int pntToPnt::relOrderPrismBlock ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< int >  basal1,
std::vector< int >  basal2,
std::vector< int > *  outBasal1,
std::vector< int > *  outBasal2 
)

Get the relative ordering of a pair of basal rings for a deformed prism/perfect prism. Outputs a vector of vectors of indices, such that the first vector is for the first basal ring, and the second vector is for the second basal ring.

Get the relative ordering of a pair of basal rings for a deformed prism/perfect prism. Outputs a vector of vectors of indices, such that the first vector is for the first basal ring, and the second vector is for the second basal ring. The input neighbour list is with respect to indices, not IDs

Definition at line 333 of file pntCorrespondence.cpp.

336  {
337  //
338  int ringSize = basal1.size(); // Number of nodes in basal1 and basal2
339  int nBonds; // Number of bonds between two parallel basal rings
340  bool isNeighbour; // Bool for checking if two atoms are neighbours or not
341  int l_k, m_k; // Elements in basal1 and basal2
342  bool isClock, isAntiClock; // Clockwise and anti-clockwise ordering of basal1
343  // and basal2
344  int iatom, jatom; // Index in the ring to start from, for basal1 and basal2
345  int currentIatom, currentJatom;
346  // ---------------
347  // Find the nearest neighbours of basal1 elements in basal2
348  nBonds = 0;
349  double infHugeNumber = 100000;
350  double leastDist = infHugeNumber;
351  int index = -1; // starting index
352  // For the first element of basal1:
353 
354  l_k = basal1[0]; // This is the atom particle C++ index
355 
356  // Search for the nearest neighbour of l_k in basal2
357  // Loop through basal2 elements
358  for (int m = 0; m < ringSize; m++) {
359  m_k = basal2[m]; // Atom index to find in the neighbour list of iatom
360 
361  // Calculate the distance
362  double dist = gen::periodicDist(yCloud, l_k, m_k);
363 
364  // Update the least distance
365  if (leastDist > dist) {
366  leastDist = dist; // This is the new least distance
367  index = m;
368  } // end of update of the least distance
369 
370  } // found element
371 
372  // If the element has been found, for l1
373  if (leastDist < infHugeNumber) {
374  isNeighbour = true;
375  iatom = 0; // index of basal1
376  jatom = index; // index of basal2
377  } // end of check
378  else {
379  std::cerr << "Something is wrong with your deformed prism.\n";
380  // Error handling
381  return 1;
382  }
383  // ---------------------------------------------------
384  // Find out if the order of basal2 is 'clockwise' or 'anticlockwise'
385  isClock = false; // init
386  isAntiClock = false;
387 
388  // atom index in the ring
389  int tempJfor, tempJback;
390 
391  tempJfor = jatom + 1;
392  tempJback = jatom - 1;
393 
394  if (jatom == ringSize - 1) {
395  tempJfor = 0;
396  tempJback = ringSize - 2;
397  }
398  if (jatom == 0) {
399  tempJfor = 1;
400  tempJback = ringSize - 1;
401  }
402 
403  int forwardJ = basal2[tempJfor];
404  int backwardJ = basal2[tempJback];
405  int currentI = basal1[iatom];
406 
407  // Check clockwise
408  double distClock = gen::periodicDist(yCloud, currentI, forwardJ);
409  double distAntiClock = gen::periodicDist(yCloud, currentI, backwardJ);
410 
411  // Clockwise
412  if (distClock < distAntiClock) {
413  isClock = true;
414  } // end of clockwise check
415  // Anti-clockwise
416  if (distAntiClock < distClock) {
417  isAntiClock = true;
418  } // end of anti-clockwise check
419  // Some error
420  if (isClock == false && isAntiClock == false) {
421  // std::cerr << "The points are equidistant.\n";
422  // Error handling
423  return 1;
424  } // end of error handling
425  // ---------------------------------------------------
426  // Get the order of basal1 and basal2
427  for (int i = 0; i < ringSize; i++) {
428  currentIatom = iatom + i;
429  if (currentIatom >= ringSize) {
430  currentIatom -= ringSize;
431  } // end of basal1 element wrap-around
432 
433  // In clockwise order
434  if (isClock) {
435  currentJatom = jatom + i;
436  if (currentJatom >= ringSize) {
437  currentJatom -= ringSize;
438  } // wrap around
439  } // end of clockwise update
440  else {
441  currentJatom = jatom - i;
442  if (currentJatom < 0) {
443  currentJatom += ringSize;
444  } // wrap around
445  } // end of anti-clockwise update
446 
447  // Add to outBasal1 and outBasal2 now
448  (*outBasal1).push_back(basal1[currentIatom]);
449  (*outBasal2).push_back(basal2[currentJatom]);
450  } //
451  //
452 
453  return 0;
454 } // end of function

◆ relOrderPrismBlock() [2/2]

int pntToPnt::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 
)

Get the relative ordering of a pair of basal rings for a deformed prism/perfect prism. Outputs a vector of vectors of indices, such that the first vector is for the first basal ring, and the second vector is for the second basal ring. The input neighbour list is with respect to indices, not IDs

Definition at line 203 of file pntCorrespondence.cpp.

207  {
208  //
209  int ringSize = basal1.size(); // Number of nodes in basal1 and basal2
210  int nBonds; // Number of bonds between two parallel basal rings
211  bool isNeighbour; // Bool for checking if two atoms are neighbours or not
212  int l_k, m_k; // Elements in basal1 and basal2
213  bool isClock, isAntiClock; // Clockwise and anti-clockwise ordering of basal1
214  // and basal2
215  int iatom, jatom; // Index in the ring to start from, for basal1 and basal2
216  int currentIatom, currentJatom;
217  // ---------------
218  // Find the nearest neighbours of basal1 elements in basal2
219  nBonds = 0;
220  isNeighbour = false;
221  // Loop through every element of basal1
222  for (int l = 0; l < ringSize; l++) {
223  l_k = basal1[l]; // This is the atom particle C++ index
224 
225  // Search for the nearest neighbour of l_k in basal2
226  // Loop through basal2 elements
227  for (int m = 0; m < ringSize; m++) {
228  m_k = basal2[m]; // Atom index to find in the neighbour list of iatom
229 
230  // Find m_k inside l_k neighbour list
231  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
232 
233  // If the element has been found, for l1
234  if (it != nList[l_k].end()) {
235  isNeighbour = true;
236  iatom = l; // index of basal1
237  jatom = m; // index of basal2
238  break;
239  } // found element
240 
241  } // end of loop through all atomIDs in basal2
242 
243  if (isNeighbour) {
244  break;
245  } // nearest neighbour found
246  } // end of loop through all the atomIDs in basal1
247 
248  if (!isNeighbour) {
249  std::cerr << "Something is wrong with your deformed prism.\n";
250  // Error handling
251  return 1;
252  }
253  // ---------------------------------------------------
254  // Find out if the order of basal2 is 'clockwise' or 'anticlockwise'
255  isClock = false; // init
256  isAntiClock = false;
257 
258  // atom index in the ring
259  int tempJfor, tempJback;
260 
261  tempJfor = jatom + 1;
262  tempJback = jatom - 1;
263 
264  if (jatom == ringSize - 1) {
265  tempJfor = 0;
266  tempJback = ringSize - 2;
267  }
268  if (jatom == 0) {
269  tempJfor = 1;
270  tempJback = ringSize - 1;
271  }
272 
273  int forwardJ = basal2[tempJfor];
274  int backwardJ = basal2[tempJback];
275  int currentI = basal1[iatom];
276 
277  // Check clockwise
278  double distClock = gen::periodicDist(yCloud, currentI, forwardJ);
279  double distAntiClock = gen::periodicDist(yCloud, currentI, backwardJ);
280 
281  // Clockwise
282  if (distClock < distAntiClock) {
283  isClock = true;
284  } // end of clockwise check
285  // Anti-clockwise
286  if (distAntiClock < distClock) {
287  isAntiClock = true;
288  } // end of anti-clockwise check
289  // Some error
290  if (isClock == false && isAntiClock == false) {
291  // std::cerr << "The points are equidistant.\n";
292  // Error handling
293  return 1;
294  } // end of error handling
295  // ---------------------------------------------------
296  // Get the order of basal1 and basal2
297  for (int i = 0; i < ringSize; i++) {
298  currentIatom = iatom + i;
299  if (currentIatom >= ringSize) {
300  currentIatom -= ringSize;
301  } // end of basal1 element wrap-around
302 
303  // In clockwise order
304  if (isClock) {
305  currentJatom = jatom + i;
306  if (currentJatom >= ringSize) {
307  currentJatom -= ringSize;
308  } // wrap around
309  } // end of clockwise update
310  else {
311  currentJatom = jatom - i;
312  if (currentJatom < 0) {
313  currentJatom += ringSize;
314  } // wrap around
315  } // end of anti-clockwise update
316 
317  // Add to outBasal1 and outBasal2 now
318  (*outBasal1).push_back(basal1[currentIatom]);
319  (*outBasal2).push_back(basal2[currentJatom]);
320  } //
321  //
322 
323  return 0;
324 } // end of function