pntCorrespondence.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS - Deferred Structural Elucidation Analysis for Molecular Simulations
3 //
4 // Copyright (c) 2018--present d-SEAMS core team
5 //
6 // This program is free software: you can redistribute it and/or modify
8 // the Open Source Initiative.
9 //
10 // A copy of the MIT License is included in the LICENSE file of this repository.
11 // You should have received a copy of the MIT License along with this program.
12 // If not, see <https://opensource.org/licenses/MIT>.
13 //-----------------------------------------------------------------------------------
14
15 #include <pntCorrespondence.hpp>
16
22 Eigen::MatrixXd pntToPnt::getPointSetRefRing(int n, int axialDim) {
23  //
24  Eigen::MatrixXd pointSet(n, 3); // Output point set for a regular polygon
25  std::vector<int> dims; // Apart from the axial dimension
26
27  // Set the axial dimension to 0, and fill up the vector of the other
28  // dimensions
29  for (int k = 0; k < 3; k++) {
30  if (k != axialDim) {
31  dims.push_back(k);
32  } // other dimensions
33  } // end of setting the axial dimension
34
35  // To get the vertices of a regular polygon, use the following formula:
36  // x[i] = r * cos(2*pi*i/n)
37  // y[i] = r * sin(2*pi*i/n)
38  // where n is the number of points and the index i is 0<=i<=n
39
40  // Loop through every particle
41  for (int i = 0; i < n; i++) {
42  pointSet(i, dims[0]) = cos((2.0 * gen::pi * i) / n); // x
43  pointSet(i, dims[1]) = sin((2.0 * gen::pi * i) / n); // y
44  // Set the axial dimension to zero
45  pointSet(i, axialDim) = 0.0;
46  } // end of loop through all the points
47
48  return pointSet;
49 } // end of function
50
58  const Eigen::MatrixXd &refPoints, int ringSize, std::vector<int> basal1,
59  std::vector<int> basal2) {
60  //
61  int nop = ringSize * 2; // Number of particles in the prism block
62  Eigen::MatrixXd pointSet(nop, 3); // Output point set for a regular prism
63  int axialDim; // Has a value of 0, 1 or 2 for x, y or z
64  double avgRadius; // Average radius within which the basal ring points lie
65  double avgHeight; // Get the average height of the prism
66  int iBasalOne,
67  iBasalTwo; // Indices for the basal1 and basal2 points in the point set
68  // --------------------------------------
69  // Get the axial dimension
70  // The axial dimension will have the largest box length
71  // Index -> axial dimension
72  // 0 -> x dim
73  // 1 -> y dim
74  // 2 -> z dim
75  axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
76  yCloud->box.begin();
77  // --------------------------------------
78  // Get the average radius
80  // --------------------------------------
81  // Get the average height of the prism
82  avgHeight = pntToPnt::getAvgHeightPrismBlock(yCloud, basal1, basal2);
83  // --------------------------------------
84  // Fill in the matrix of the point set
85  //
86  // Fill basal1 first, and then basal2
87  for (int i = 0; i < ringSize; i++) {
88  iBasalOne = i; // index in point set for basal1 point
89  iBasalTwo = i + ringSize; // index in point set for basal2 point
90  // Fill up the dimensions
91  for (int k = 0; k < 3; k++) {
92  // For the axial dimension
93  if (k == axialDim) {
94  pointSet(iBasalOne, k) = 0.0; // basal1
95  pointSet(iBasalTwo, k) = avgHeight; // basal2
96  } // end of filling up the axial dimension
97  else {
98  pointSet(iBasalOne, k) = avgRadius * refPoints(i, k); // basal1
99  pointSet(iBasalTwo, k) = avgRadius * refPoints(i, k); // basal2
100  } // fill up the non-axial dimension coordinate dimension
101  } // Fill up all the dimensions
102  } // end of filling in the point set
103  // --------------------------------------
104  return pointSet;
105 } // end of function
106
112  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
113  std::vector<int> basal1, std::vector<int> basal2) {
114  //
116  std::vector<double> centroid1, centroid2;
117  std::vector<double> dist; // Distances
118  int ringSize = basal1.size(); // Number of nodes in the basal rings
119  double r1, r2; // Distance from the centroid
120  int iatom; // Atom index in yCloud for the current point in basal1
121  int jatom; // Atom index in yCloud for the current point in basal2
122  // -----------------------
123  // Calculate the centroid for basal1 and basal2
124  centroid1.resize(3);
125  centroid2.resize(3); // init
126  // Loop through basal1 and basal2
127  for (int i = 0; i < ringSize; i++) {
128  // For basal1
129  centroid1[0] += yCloud->pts[basal1[i]].x;
130  centroid1[1] += yCloud->pts[basal1[i]].y;
131  centroid1[2] += yCloud->pts[basal1[i]].z;
132  //
133  // For basal2
134  centroid2[0] += yCloud->pts[basal2[i]].x;
135  centroid2[1] += yCloud->pts[basal2[i]].y;
136  centroid2[2] += yCloud->pts[basal2[i]].z;
137  //
138  } // end of loop through basal1 and basal2
139  // Normalize by the number of nodes
140  for (int k = 0; k < 3; k++) {
141  centroid1[k] /= ringSize;
142  centroid2[k] /= ringSize;
143  } // end of dividing by the number
144  // Calculated the centroid for basal1 and basal2!
145  // -----------------------
146  // Calculate the distances of each point from the respective centroid
147  for (int i = 0; i < ringSize; i++) {
148  // Get the current point in basal1 and basal2
149  iatom = basal1[i];
150  jatom = basal2[i];
151  // Get the distance from the respective centroid
152  // basal1
153  r1 = gen::unWrappedDistFromPoint(yCloud, iatom, centroid1);
154  // basal2
155  r2 = gen::unWrappedDistFromPoint(yCloud, jatom, centroid2);
156  // Update the distances
157  dist.push_back(r1);
158  dist.push_back(r2);
159  } // Loop through every point in basal1 and basal2
160  // -----------------------
161  // Calculate the average, excluding the outliers
163  // -----------------------
165 } // end of function
166
172  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
173  std::vector<int> basal1, std::vector<int> basal2) {
174  //
175  double avgHeight = 0.0;
176  double r_ij; // Distance between a point in basal1 and basal2
177  std::vector<double> dist; // Distances
178  int ringSize = basal1.size(); // Number of nodes in the basal rings
179  int iatom; // Atom index in yCloud for the current point in basal1
180  int jatom; // Atom index in yCloud for the current point in basal2
181  // -----------------------
182  // Calculate the distances of each point from the respective centroid
183  for (int i = 0; i < ringSize; i++) {
184  // Get the current point in basal1 and basal2
185  iatom = basal1[i];
186  jatom = basal2[i];
187  // Get the distance between a point in basal1 and the corresponding point in
188  // basal2
189  r_ij = gen::periodicDist(yCloud, iatom, jatom);
190  // Update the distances
191  dist.push_back(r_ij);
192  } // Loop through every point in basal1 and basal2
193  // -----------------------
194  // Calculate the average, excluding the outliers
195  avgHeight = gen::getAverageWithoutOutliers(dist);
196  // -----------------------
197  return avgHeight;
198 } // end of function
199
208  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
209  std::vector<int> basal1, std::vector<int> basal2,
210  std::vector<std::vector<int>> nList, std::vector<int> *outBasal1,
211  std::vector<int> *outBasal2) {
212  //
213  int ringSize = basal1.size(); // Number of nodes in basal1 and basal2
214  int nBonds; // Number of bonds between two parallel basal rings
215  bool isNeighbour; // Bool for checking if two atoms are neighbours or not
216  int l_k, m_k; // Elements in basal1 and basal2
217  bool isClock, isAntiClock; // Clockwise and anti-clockwise ordering of basal1
218  // and basal2
219  int iatom, jatom; // Index in the ring to start from, for basal1 and basal2
220  int currentIatom, currentJatom;
221  // ---------------
222  // Find the nearest neighbours of basal1 elements in basal2
223  nBonds = 0;
224  isNeighbour = false;
225  // Loop through every element of basal1
226  for (int l = 0; l < ringSize; l++) {
227  l_k = basal1[l]; // This is the atom particle C++ index
228
229  // Search for the nearest neighbour of l_k in basal2
230  // Loop through basal2 elements
231  for (int m = 0; m < ringSize; m++) {
232  m_k = basal2[m]; // Atom index to find in the neighbour list of iatom
233
234  // Find m_k inside l_k neighbour list
235  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
236
237  // If the element has been found, for l1
238  if (it != nList[l_k].end()) {
239  isNeighbour = true;
240  iatom = l; // index of basal1
241  jatom = m; // index of basal2
242  break;
243  } // found element
244
245  } // end of loop through all atomIDs in basal2
246
247  if (isNeighbour) {
248  break;
249  } // nearest neighbour found
250  } // end of loop through all the atomIDs in basal1
251
252  if (!isNeighbour) {
253  std::cerr << "Something is wrong with your deformed prism.\n";
254  // Error handling
255  return 1;
256  }
257  // ---------------------------------------------------
258  // Find out if the order of basal2 is 'clockwise' or 'anticlockwise'
259  isClock = false; // init
260  isAntiClock = false;
261
262  // atom index in the ring
263  int tempJfor, tempJback;
264
265  tempJfor = jatom + 1;
266  tempJback = jatom - 1;
267
268  if (jatom == ringSize - 1) {
269  tempJfor = 0;
270  tempJback = ringSize - 2;
271  }
272  if (jatom == 0) {
273  tempJfor = 1;
274  tempJback = ringSize - 1;
275  }
276
277  int forwardJ = basal2[tempJfor];
278  int backwardJ = basal2[tempJback];
279  int currentI = basal1[iatom];
280
281  // Check clockwise
282  double distClock = gen::periodicDist(yCloud, currentI, forwardJ);
283  double distAntiClock = gen::periodicDist(yCloud, currentI, backwardJ);
284
285  // Clockwise
286  if (distClock < distAntiClock) {
287  isClock = true;
288  } // end of clockwise check
289  // Anti-clockwise
290  if (distAntiClock < distClock) {
291  isAntiClock = true;
292  } // end of anti-clockwise check
293  // Some error
294  if (isClock == false && isAntiClock == false) {
295  // std::cerr << "The points are equidistant.\n";
296  // Error handling
297  return 1;
298  } // end of error handling
299  // ---------------------------------------------------
300  // Get the order of basal1 and basal2
301  for (int i = 0; i < ringSize; i++) {
302  currentIatom = iatom + i;
303  if (currentIatom >= ringSize) {
304  currentIatom -= ringSize;
305  } // end of basal1 element wrap-around
306
307  // In clockwise order
308  if (isClock) {
309  currentJatom = jatom + i;
310  if (currentJatom >= ringSize) {
311  currentJatom -= ringSize;
312  } // wrap around
313  } // end of clockwise update
314  else {
315  currentJatom = jatom - i;
316  if (currentJatom < 0) {
317  currentJatom += ringSize;
318  } // wrap around
319  } // end of anti-clockwise update
320
321  // Add to outBasal1 and outBasal2 now
322  (*outBasal1).push_back(basal1[currentIatom]);
323  (*outBasal2).push_back(basal2[currentJatom]);
324  } //
325  //
326
327  return 0;
328 } // end of function
329
338  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
339  std::vector<int> basal1, std::vector<int> basal2,
340  std::vector<int> *outBasal1, std::vector<int> *outBasal2) {
341  //
342  int ringSize = basal1.size(); // Number of nodes in basal1 and basal2
343  int nBonds; // Number of bonds between two parallel basal rings
344  bool isNeighbour; // Bool for checking if two atoms are neighbours or not
345  int l_k, m_k; // Elements in basal1 and basal2
346  bool isClock, isAntiClock; // Clockwise and anti-clockwise ordering of basal1
347  // and basal2
348  int iatom, jatom; // Index in the ring to start from, for basal1 and basal2
349  int currentIatom, currentJatom;
350  // ---------------
351  // Find the nearest neighbours of basal1 elements in basal2
352  nBonds = 0;
353  double infHugeNumber = 100000;
354  double leastDist = infHugeNumber;
355  int index = -1; // starting index
356  // For the first element of basal1:
357
358  l_k = basal1[0]; // This is the atom particle C++ index
359
360  // Search for the nearest neighbour of l_k in basal2
361  // Loop through basal2 elements
362  for (int m = 0; m < ringSize; m++) {
363  m_k = basal2[m]; // Atom index to find in the neighbour list of iatom
364
365  // Calculate the distance
366  double dist = gen::periodicDist(yCloud, l_k, m_k);
367
368  // Update the least distance
369  if (leastDist > dist) {
370  leastDist = dist; // This is the new least distance
371  index = m;
372  } // end of update of the least distance
373
374  } // found element
375
376  // If the element has been found, for l1
377  if (leastDist < infHugeNumber) {
378  isNeighbour = true;
379  iatom = 0; // index of basal1
380  jatom = index; // index of basal2
381  } // end of check
382  else {
383  std::cerr << "Something is wrong with your deformed prism.\n";
384  // Error handling
385  return 1;
386  }
387  // ---------------------------------------------------
388  // Find out if the order of basal2 is 'clockwise' or 'anticlockwise'
389  isClock = false; // init
390  isAntiClock = false;
391
392  // atom index in the ring
393  int tempJfor, tempJback;
394
395  tempJfor = jatom + 1;
396  tempJback = jatom - 1;
397
398  if (jatom == ringSize - 1) {
399  tempJfor = 0;
400  tempJback = ringSize - 2;
401  }
402  if (jatom == 0) {
403  tempJfor = 1;
404  tempJback = ringSize - 1;
405  }
406
407  int forwardJ = basal2[tempJfor];
408  int backwardJ = basal2[tempJback];
409  int currentI = basal1[iatom];
410
411  // Check clockwise
412  double distClock = gen::periodicDist(yCloud, currentI, forwardJ);
413  double distAntiClock = gen::periodicDist(yCloud, currentI, backwardJ);
414
415  // Clockwise
416  if (distClock < distAntiClock) {
417  isClock = true;
418  } // end of clockwise check
419  // Anti-clockwise
420  if (distAntiClock < distClock) {
421  isAntiClock = true;
422  } // end of anti-clockwise check
423  // Some error
424  if (isClock == false && isAntiClock == false) {
425  // std::cerr << "The points are equidistant.\n";
426  // Error handling
427  return 1;
428  } // end of error handling
429  // ---------------------------------------------------
430  // Get the order of basal1 and basal2
431  for (int i = 0; i < ringSize; i++) {
432  currentIatom = iatom + i;
433  if (currentIatom >= ringSize) {
434  currentIatom -= ringSize;
435  } // end of basal1 element wrap-around
436
437  // In clockwise order
438  if (isClock) {
439  currentJatom = jatom + i;
440  if (currentJatom >= ringSize) {
441  currentJatom -= ringSize;
442  } // wrap around
443  } // end of clockwise update
444  else {
445  currentJatom = jatom - i;
446  if (currentJatom < 0) {
447  currentJatom += ringSize;
448  } // wrap around
449  } // end of anti-clockwise update
450
451  // Add to outBasal1 and outBasal2 now
452  (*outBasal1).push_back(basal1[currentIatom]);
453  (*outBasal2).push_back(basal2[currentJatom]);
454  } //
455  //
456
457  return 0;
458 } // end of function
459
465  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
466  std::vector<int> basalRing, int startingIndex) {
467  //
468  //
469  int dim = 3; // Number of dimensions (3)
470  int ringSize = basalRing.size(); // Number of nodes in the ring
471  int nop = basalRing.size(); // Number of particles in the basal ring
472  Eigen::MatrixXd pointSet(nop, dim); // Output set of 3D points
473  // Indices for the first and second basal rings being filled
474  int currentPosition; // Current index in basalRing
475  int index; // Index in yCloud
476
477  // Check
478  if (startingIndex >= ringSize || startingIndex < 0) {
479  startingIndex = 0;
480  } // end of check
481
482  // Beginning from the starting index, get the points in the basal ring
483  for (int i = 0; i < nop; i++) {
484  //
485  // Getting currentIndex
486  currentPosition = startingIndex + i;
487  // Wrapping around for the ring
488  if (currentPosition >= ringSize) {
489  currentPosition -= ringSize;
490  } // end of wrap-around
491  //
492  // -------------------
493  // Basal ring points
494  index = basalRing[currentPosition]; // Index of the current point
495  pointSet(i, 0) = yCloud->pts[index].x; // x coord
496  pointSet(i, 1) = yCloud->pts[index].y; // y coord
497  pointSet(i, 2) = yCloud->pts[index].z; // z coord
498  } // end of point filling from the relative ordering vector of vectors
499
500  // Return the set of points
501  return pointSet;
502 } // end of function
503
509  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
510  std::vector<int> basal1, std::vector<int> basal2, int startingIndex) {
511  //
512  //
513  int dim = 3; // Number of dimensions (3)
514  int ringSize = basal1.size(); // Number of nodes in the ring
515  int nop = ringSize * 2; // Number of particles in prism block
516  Eigen::MatrixXd pointSet(nop, dim); // Output set of 3D points
517  // Indices for the first and second basal rings being filled
518  int currentPosition; // Current position in the basal ring vectors
519  int iatom, jatom; // Current index in the point set
520  int iatomIndex, jatomIndex; // Index in yCloud corresponding to the basal1
521  // and basal2 points
522
523  // Fill up basal1 points, followed by basal2 points
524
525  // Check
526  if (startingIndex >= ringSize || startingIndex < 0) {
527  startingIndex = 0;
528  } // end of check
529
530  // Beginning from the starting index, get the points in the basal ring
531  for (int i = 0; i < ringSize; i++) {
532  //
533  // Getting currentIndex
534  currentPosition = startingIndex + i;
535  // Wrapping around for the ring
536  if (currentPosition >= ringSize) {
537  currentPosition -= ringSize;
538  } // end of wrap-around
539  //
540  // -------------------
541  // Basal1 ring points
542  iatomIndex =
543  basal1[currentPosition]; // Index of the current point in basal1
544  iatom = i; // index in the point set being filled
545  pointSet(iatom, 0) = yCloud->pts[iatomIndex].x; // x coord
546  pointSet(iatom, 1) = yCloud->pts[iatomIndex].y; // y coord
547  pointSet(iatom, 2) = yCloud->pts[iatomIndex].z; // z coord
548  // -------------------
549  // Basal2 ring points
550  jatomIndex =
551  basal2[currentPosition]; // Index of the current point in basal2
552  jatom = i + ringSize; // index in the point set being filled
553  pointSet(jatom, 0) = yCloud->pts[jatomIndex].x; // x coord
554  pointSet(jatom, 1) = yCloud->pts[jatomIndex].y; // y coord
555  pointSet(jatom, 2) = yCloud->pts[jatomIndex].z; // z coord
556  } // end of point filling from the relative ordering vector of vectors
557
558  // Return the set of points
559  return pointSet;
560 } // end of function
561
562 // ---------------------------------------------------
569  //
570  Eigen::MatrixXd pointSet(12, 3); // Reference point set (Eigen matrix)
572  setCloud; // PointCloud for holding the reference point values
573
574  if (type == ring::HCbasal) {
575  // Read in the XYZ file
576  std::string fileName = "templates/hc.xyz";
577  //
579  int n = setCloud.nop; // Number of points in the reference point set
580  Eigen::MatrixXd pointSet(n, 3); // Output point set for a regular polygon
581
582  // Loop through every particle
583  for (int i = 0; i < n; i++) {
584  pointSet(i, 0) = setCloud.pts[i].x; // x
585  pointSet(i, 1) = setCloud.pts[i].y; // y
586  pointSet(i, 2) = setCloud.pts[i].z; // z
587  } // end of looping through every particle
588
589  // Return the filled point set
590  return pointSet;
591  } // end of reference point set for HCs
592  // else {
593  // // ddc
594  // Eigen::MatrixXd pointSet(14, 3); // Output point set for a regular
595  // polygon return pointSet;
596  // } // DDCs
597
598  // Eigen::MatrixXd pointSet(14, 3); // Output point set for a regular polygon
599  // // Never happens
600  return pointSet;
601
602 } // end of function
603
610  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
611  std::vector<int> basal1, std::vector<int> basal2,
612  std::vector<std::vector<int>> nList, std::vector<int> *matchedBasal1,
613  std::vector<int> *matchedBasal2) {
614  //
615  int l1 = basal1[0]; // First element of basal1
616  int l2 = basal1[1]; // Second element of basal1
617  int ringSize = basal1.size(); // Number of nodes in the basal rings
618  bool neighOne, neighTwo; // Basal2 element is the neighbour of l1 or l2
619  bool neighbourFound; // neighbour found
620  int m_k; // Element of basal2 which is a neighbour of l1 or l2
621  int m_kIndex; // Index of basal2 which is a neighbour of l1 or l2
622  int iatom; // Current element of basal2 being searched for
623  int nextBasal1Element; // Next bonded element of basal1
624  int nextBasal2Element; // Next bonded element of basal2
625  bool isReversedOrder; // basal2 is reversed wrt basal1
626  int index;
627
628  (*matchedBasal1).resize(ringSize);
629  (*matchedBasal2).resize(ringSize);
630
631  // Search to see if l1 or l2 is a neighbour
632  // -------------------
633  // Searching for the neighbours of l1 or l2
634  for (int i = 0; i < ringSize; i++) {
635  iatom = basal2[i]; // Element of basal2
636  // Search for the current element in the neighbour list of l1
637  auto it = std::find(nList[l1].begin() + 1, nList[l1].end(), iatom);
638  // If iatom is the neighbour of l1
639  if (it != nList[l1].end()) {
640  neighbourFound = true;
641  neighOne = true;
642  m_k = iatom; // Element found
643  m_kIndex = i; // Index in basal2
644  nextBasal1Element = basal1[2];
645  break;
646  } // iatom is the neighbour of l1
647  // l2 neighbour check
648  else {
649  auto it1 = std::find(nList[l2].begin() + 1, nList[l2].end(), iatom);
650  // If iatom is the neighbour of l2
651  if (it1 != nList[l2].end()) {
652  neighbourFound = true;
653  neighOne = false;
654  neighTwo = true;
655  nextBasal1Element = basal1[3];
656  m_k = iatom; // Element found
657  m_kIndex = i; // Index in basal2
658  break;
659  } // iatom is the neighbour of l2
660  } // Check for the neighbour of l2
661  } // end of search through basal2 elements
662  //
663  // -------------------
664  // If a neighbour was not found, then there is some mistake
665  if (!neighbourFound) {
666  // std::cerr << "This is not an HC\n";
667  return 1;
668  }
669
670  // ------------------------------
671  //
672  // This element should be the nearest neighbour of the corresponding element
673  // in basal2
674  //
675  // Testing the original order
676  index = m_kIndex + 2;
677  // wrap-around
678  if (index >= ringSize) {
679  index -= ringSize;
680  } // wrap-around
681  nextBasal2Element = basal2[index];
682  // Search for the next basal element
683  auto it = std::find(nList[nextBasal1Element].begin() + 1,
684  nList[nextBasal1Element].end(), nextBasal2Element);
685  // If this element is found, then the original order is correct
686  if (it != nList[nextBasal1Element].end()) {
687  // Fill up the temporary vector with basal2 elements
688  for (int i = 0; i < ringSize; i++) {
689  index = m_kIndex + i; // index in basal2
690  // wrap-around
691  if (index >= ringSize) {
692  index -= ringSize;
693  } // end of wrap-around
694  (*matchedBasal2)[i] = basal2[index]; // fill up values
695  } // end of filling up tempBasal2
696  } // the original order is correct
697  else {
698  //
699  index = m_kIndex - 2;
700  // wrap-around
701  if (index < 0) {
702  index += ringSize;
703  } // wrap-around
704  nextBasal2Element = basal2[index];
705  // Search for the next basal element
706  auto it = std::find(nList[nextBasal1Element].begin() + 1,
707  nList[nextBasal1Element].end(), nextBasal2Element);
708  // If this element is found, then the original order is correct
709  if (it != nList[nextBasal1Element].end()) {
710  // Fill up the temporary vector with basal2 elements
711  for (int i = 0; i < ringSize; i++) {
712  index = m_kIndex - i; // index in basal2
713  // wrap-around
714  if (index < 0) {
715  index += ringSize;
716  } // end of wrap-around
717  (*matchedBasal2)[i] = basal2[index]; // fill up values
718  } // end of filling up tempBasal2
719
720  }
721  //
722  else {
723  // std::cerr << "not an HC\n";
724  return 1;
725  }
726  //
727  } // the reversed order is correct!
728  // ------------------------------
729  // Fill up basal1
730  (*matchedBasal1) = basal1;
731  return 0;
732 } // end of the function
733
741  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
742  std::vector<int> basal1, std::vector<int> basal2, int startingIndex) {
743  Eigen::MatrixXd pointSet(12, 3);
744  int iatomIndex, jatomIndex; // Current atom index in yCloud, according to
745  // basal1 and basal2 respectively
746  int iPnt; // Current index in the Eigen matrix pointSet
747  int cageSize = 12; // Number of points in the cage
748  std::vector<int> newBasal1, newBasal2;
749  std::array<double, 3> dr; // Components of the distance
750  int iatomOne; // Index of the first atom
751
752  // Checks and balances
753  //
754  if (startingIndex > 5 || startingIndex < 0) {
755  startingIndex = 0;
756  } // no invalid starting index
757
758  // Change the order
759  if (startingIndex > 0) {
760  for (int k = 0; k < 6; k++) {
761  iPnt = k + startingIndex;
762  if (iPnt >= 6) {
763  iPnt -= 6;
764  } // wrap-around
765  newBasal1.push_back(basal1[iPnt]);
766  newBasal2.push_back(basal2[iPnt]);
767  } // change the order
768  } // end of filling for startingIndex>0
769  else {
770  newBasal1 = basal1;
771  newBasal2 = basal2;
772  } // end of filling up the reordered basal rings
773
774  //
775  // FIRST POINT
776  // basal1
777  iatomOne = newBasal1[0];
778  pointSet(0, 0) = yCloud->pts[iatomOne].x;
779  pointSet(0, 1) = yCloud->pts[iatomOne].y;
780  pointSet(0, 2) = yCloud->pts[iatomOne].z;
781  // basal2
782  jatomIndex = newBasal2[0];
783  // Get the distance from basal1
784  dr = gen::relDist(yCloud, iatomOne, jatomIndex);
785
786  // basal2
787  pointSet(6, 0) = yCloud->pts[iatomOne].x + dr[0];
788  pointSet(6, 1) = yCloud->pts[iatomOne].y + dr[1];
789  pointSet(6, 2) = yCloud->pts[iatomOne].z + dr[2];
790  //
791  // Loop through the rest of the points
792  for (int i = 1; i < 6; i++) {
793  // basal1
794  iatomIndex = newBasal1[i]; // Atom index to be filled for basal1
795  jatomIndex = newBasal2[i]; // Atom index to be filled for basal2
796  //
797  // Get the distance from the first atom
798  dr = gen::relDist(yCloud, iatomOne, iatomIndex);
799  //
800  pointSet(i, 0) = yCloud->pts[iatomOne].x + dr[0];
801  pointSet(i, 1) = yCloud->pts[iatomOne].y + dr[1];
802  pointSet(i, 2) = yCloud->pts[iatomOne].z + dr[2];
803  //
804  // Get the distance from the first atom
805  dr = gen::relDist(yCloud, iatomOne, jatomIndex);
806  // basal2
807  pointSet(i + 6, 0) = yCloud->pts[iatomOne].x + dr[0];
808  pointSet(i + 6, 1) = yCloud->pts[iatomOne].y + dr[1];
809  pointSet(i + 6, 2) = yCloud->pts[iatomOne].z + dr[2];
810  //
811  } // end of loop
812
813  return pointSet;
814 }
815
847  std::vector<cage::Cage> cageList) {
848  //
849  std::vector<int> ddcOrder; // Order of the particles in the DDC.
850  int nop = 14; // Number of elements in the DDC
851  int ringSize = 6; // Number of nodes in each ring
852  int iring, jring; // Ring indices of the DDC rings
853  int iatom; // Atom index in the equatorial ring
854  std::vector<int> peripheral1, peripheral2; // Vectors holding the neighbours
855  // of (l1,l3,l5) and (l2,l4,l6)
856  int apex1, apex2;
857  bool neighbourFound; // Neighbour found
858  int nextI, prevI, nextJ, prevJ; // elements around iatom and jatom
859  int jatomIndex, atomIndex;
860
861  // Add the equatorial ring particles
862  //
863  iring = cageList[index]
864  .rings[0]; // Equatorial ring index in rings vector of vectors
865  //
866  // Loop through all the atoms of the equatorial ring
867  for (int i = 0; i < ringSize; i++) {
868  ddcOrder.push_back(rings[iring][i]);
869  } // end of adding equatorial ring particles
870  //
871  // ------------------------------
872  // Find the neighbouring atoms of (l1, l3 and l5) and (l2,l4,l6) by looping
873  // through the peripheral rings
874  //
875  for (int i = 0; i < ringSize; i++) {
876  // Init
877  iatom = ddcOrder[i]; // element of the equatorial ring
878  // Get the next and previous elements of iatom
879  // next element
880  if (i == ringSize - 1) {
881  nextI = ddcOrder[0];
882  prevI = ddcOrder[i - 1];
883  } // if last element
884  else if (i == 0) {
885  nextI = ddcOrder[i + 1];
886  prevI = ddcOrder[5]; // wrapped
887  } else {
888  nextI = ddcOrder[i + 1];
889  prevI = ddcOrder[i - 1];
890  }
891  // ------------------
892  // Search all the peripheral rings for iatom, get the nearest neighbours and
893  // apex1 and apex2, which are bonded to the nearest neighbours of the
894  // equatorial ring (thus they are second nearest neighbours)
895  for (int j = 1; j <= ringSize; j++) {
896  //
897  jring = cageList[index].rings[j]; // Peripheral ring
898  // Search for iatom
899  //
900  auto it = std::find(rings[jring].begin(), rings[jring].end(), iatom);
901  // If iatom was found in jring peripheral ring
902  if (it != rings[jring].end()) {
903  // ----------------
904  // Atom index for iatom found
905  jatomIndex = std::distance(rings[jring].begin(), it);
906  // ----------------
907  // Get the next and previous element
908  //
909  // Next atom
910  atomIndex = jatomIndex + 1;
911  // wrap-around
912  if (atomIndex >= ringSize) {
913  atomIndex -= ringSize;
914  } // end of wrap-around
915  nextJ = rings[jring][atomIndex];
916  // Previous atom
917  atomIndex = jatomIndex - 1;
918  if (atomIndex < 0) {
919  atomIndex += ringSize;
920  } // end of wrap-around
921  prevJ = rings[jring][atomIndex];
922  // ----------------
923  // Check to see if nextJ or prevJ are different from nextI and prevI
924  //
925  // Checking prevJ
926  if (prevJ != nextI && prevJ != prevI) {
927  // Add to the vector
928  // if even peripheral1
929  if (i % 2 == 0) {
930  peripheral1.push_back(prevJ);
931  // Get apex1 for i=0
932  if (i == 0) {
933  // Go back two elements from jatomIndex
934  atomIndex = jatomIndex - 2;
935  // wrap-around
936  if (atomIndex < 0) {
937  atomIndex += ringSize;
938  } // end of wrap-around
939  apex1 = rings[jring][atomIndex];
940  } // apex1 for i=0
941  } // peripheral1
942  // if odd peripheral2
943  else {
944  peripheral2.push_back(prevJ);
945  // Get apex2 for i=0
946  if (i == 1) {
947  // Go back two elements from jatomIndex
948  atomIndex = jatomIndex - 2;
949  // wrap-around
950  if (atomIndex < 0) {
951  atomIndex += ringSize;
952  } // end of wrap-around
953  apex2 = rings[jring][atomIndex];
954  } // apex2 for i=1
955  } // peripheral 2
956  // Get apex1 or apex2 for i=0 or i=1
957  break;
958  } // check if prevJ is the atom to be added, bonded to iatom
959  // ----------------------
960  //
961  // Checking nextJ
962  if (nextJ != nextI && nextJ != prevI) {
963  // Add to the vector
964  // if even peripheral1
965  if (i % 2 == 0) {
966  peripheral1.push_back(nextJ);
967  // Get apex1 for i=0
968  if (i == 0) {
969  // Go foward two elements from jatomIndex
970  atomIndex = jatomIndex + 2;
971  // wrap-around
972  if (atomIndex >= ringSize) {
973  atomIndex -= ringSize;
974  } // end of wrap-around
975  apex1 = rings[jring][atomIndex];
976  } // apex1 for i=0
977  } // peripheral1
978  // if odd peripheral2
979  else {
980  peripheral2.push_back(prevJ);
981  // Get apex2 for i=0
982  if (i == 1) {
983  // Go foward two elements from jatomIndex
984  atomIndex = jatomIndex + 2;
985  // wrap-around
986  if (atomIndex >= ringSize) {
987  atomIndex -= ringSize;
988  } // end of wrap-around
989  apex2 = rings[jring][atomIndex];
990  } // apex2 for i=1
991  } // peripheral 2
992  // Get apex1 or apex2 for i=0 or i=1
993  break;
994  } // check if prevJ is the atom to be added, bonded to iatom
995  // ----------------
996  } // end of check for iatom in jring
997  } // end of search for the peripheral rings
998  // ------------------
999  } // end of looping through the elements of the equatorial ring
1000  // ------------------------------
1001  // Update ddcOrder with peripheral1, followed by apex1, peripheral2 and apex2
1002  //
1003  // peripheral1
1004  for (int i = 0; i < peripheral1.size(); i++) {
1005  ddcOrder.push_back(peripheral1[i]);
1006  } // end of adding peripheral1
1007  //
1008  // apex1
1009  ddcOrder.push_back(apex1);
1010  //
1011  // peripheral2
1012  for (int i = 0; i < peripheral2.size(); i++) {
1013  ddcOrder.push_back(peripheral2[i]);
1014  } // end of adding peripheral2
1015  //
1016  // apex2
1017  ddcOrder.push_back(apex2);
1018  //
1019  return ddcOrder;
1020 } // end of the function
1021
1043  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1044  std::vector<int> ddcOrder, int startingIndex) {
1045  int nop = 14; // Number of elements in the DDC
1046  int ringSize = 6; // Six nodes in the rings
1047  Eigen::MatrixXd pointSet(nop, 3); // Output point set
1048  int peripheralStartingIndex; // Index using which the elements of peripheral1
1049  // and peripheral2 will be wrapped
1050  std::vector<int> wrappedDDC; // Changed order of the DDC: should be the same
1051  // size as the original ddcOrder vector
1052  int currentIndex; // Current index
1053  // Variables for filling the point set
1054  int iatomIndex, jatomIndex;
1055  std::array<double, 3> dr; // Components of the distance
1056
1057  if (startingIndex == 0) {
1058  wrappedDDC = ddcOrder;
1059  } // if the order does not have to be changed
1060  else {
1061  wrappedDDC.resize(nop); // Should be the same size as ddcOrder
1062  // ------------------
1063  // EQUATORIAL RING
1064  // Change the order of the equatorial ring
1065  for (int i = 0; i < ringSize; i++) {
1066  currentIndex = startingIndex + i;
1067  // Wrap-around
1068  if (currentIndex >= ringSize) {
1069  currentIndex -= ringSize;
1070  } // end of wrap-around
1071  // Update the equatorial ring
1072  wrappedDDC[i] = ddcOrder[currentIndex];
1073  } // end of wrapping the equatorial ring
1074  // ------------------
1075  // PERIPHERAL RINGS
1076  //
1077  // The index of the peripheral ring starting indices
1078  if (startingIndex <= 1) {
1079  peripheralStartingIndex = 0;
1080  } else if (startingIndex > 1 && startingIndex <= 3) {
1081  peripheralStartingIndex = 1;
1082  } else {
1083  peripheralStartingIndex = 2;
1084  }
1085  //
1086  // Update the portions of the wrappedDDC vector
1087  for (int i = 6; i < 9; i++) {
1088  currentIndex = i + peripheralStartingIndex;
1089  // wrap-around
1090  if (currentIndex <= 9) {
1091  currentIndex -= 3;
1092  } // end of wrap-around
1093  wrappedDDC[i] = ddcOrder[currentIndex];
1094  } // peripheral1
1095  //
1096  // Update the peripheral2 portions of the wrappedDDC vector
1097  for (int i = 10; i < 13; i++) {
1098  currentIndex = i + peripheralStartingIndex;
1099  // wrap-around
1100  if (currentIndex <= 13) {
1101  currentIndex -= 3;
1102  } // end of wrap-around
1103  wrappedDDC[i] = ddcOrder[currentIndex];
1104  } // peripheral2
1105  // ------------------
1106  // SECOND-NEAREST NEIGHBOURS
1107  wrappedDDC[9] = ddcOrder[9]; // apex1
1108  wrappedDDC[13] = ddcOrder[13]; // apex2
1109  // ------------------
1110  } // the order of the DDC has to be changed
1111
1112  // FILL UP THE EIGEN MATRIX
1113  iatomIndex = wrappedDDC[0]; // first point
1114  pointSet(0, 0) = yCloud->pts[iatomIndex].x;
1115  pointSet(0, 1) = yCloud->pts[iatomIndex].y;
1116  pointSet(0, 2) = yCloud->pts[iatomIndex].z;
1117  // Loop through the rest of the equatorial ring points
1118  for (int i = 1; i < 14; i++) {
1119  //
1120  jatomIndex = wrappedDDC[i]; // Atom index to be filled
1121  // Get the distance from iatomIndex
1122  dr = gen::relDist(yCloud, iatomIndex, jatomIndex);
1123
1124  // basal2
1125  pointSet(i, 0) = yCloud->pts[iatomIndex].x + dr[0];
1126  pointSet(i, 1) = yCloud->pts[iatomIndex].y + dr[1];
1127  pointSet(i, 2) = yCloud->pts[iatomIndex].z + dr[2];
1128  } // end of loop
1129
1130  return pointSet;
1131 }
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:200
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:108
double unWrappedDistFromPoint(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, std::vector< double > singlePoint)
Definition: generic.hpp:138
const double pi
Definition: generic.hpp:54
double getAverageWithoutOutliers(std::vector< double > inpVec)
Get the average, after excluding the outliers, using quartiles.
Definition: generic.cpp:168
int nop
Current frame number.
Definition: mol_sys.hpp:173
std::vector< S > pts
Definition: mol_sys.hpp:171
strucType
Definition: ring.hpp:116
@ HCbasal
The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an H...
Definition: ring.hpp:119
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:73
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:170
This contains per-particle information.
Definition: mol_sys.hpp:149