pntCorrespondence.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <pntCorrespondence.hpp>
12 
18 Eigen::MatrixXd pntToPnt::getPointSetRefRing(int n, int axialDim) {
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
46 
54  const Eigen::MatrixXd &refPoints, int ringSize, std::vector<int> basal1,
55  std::vector<int> basal2) {
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
102 
108  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
109  std::vector<int> basal1, std::vector<int> basal2) {
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
162 
168  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
169  std::vector<int> basal1, std::vector<int> basal2) {
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
195 
204  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
205  std::vector<int> basal1, std::vector<int> basal2,
206  std::vector<std::vector<int>> nList, std::vector<int> *outBasal1,
207  std::vector<int> *outBasal2) {
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
325 
334  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
335  std::vector<int> basal1, std::vector<int> basal2,
336  std::vector<int> *outBasal1, std::vector<int> *outBasal2) {
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
455 
461  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
462  std::vector<int> basalRing, int startingIndex) {
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
499 
505  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
506  std::vector<int> basal1, std::vector<int> basal2, int startingIndex) {
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
557 
558 // ---------------------------------------------------
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
599 
606  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
607  std::vector<int> basal1, std::vector<int> basal2,
608  std::vector<std::vector<int>> nList, std::vector<int> *matchedBasal1,
609  std::vector<int> *matchedBasal2) {
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
729 
737  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
738  std::vector<int> basal1, std::vector<int> basal2, int startingIndex) {
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 }
811 
843  std::vector<cage::Cage> cageList) {
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
1017 
1039  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1040  std::vector<int> ddcOrder, int startingIndex) {
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 }
T distance(T... args)
T find(T... args)
std::array< double, 3 > relDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Definition: generic.hpp:196
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 unWrappedDistFromPoint(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, std::vector< double > singlePoint)
Definition: generic.hpp:134
const double pi
Definition: generic.hpp:50
double getAverageWithoutOutliers(std::vector< double > inpVec)
Get the average, after excluding the outliers, using quartiles.
Definition: generic.cpp:164
int nop
Current frame number.
Definition: mol_sys.hpp:169
std::vector< S > pts
Definition: mol_sys.hpp:167
strucType
Definition: ring.hpp:112
@ 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
T max_element(T... args)
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)
T push_back(T... args)
T resize(T... args)
T size(T... args)
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166
This contains per-particle information.
Definition: mol_sys.hpp:145