topo_bulk.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 <topo_bulk.hpp>
12 
13 // -----------------------------------------------------------------------------------------------------
14 // DDC / HC ALGORITHMS
15 // -----------------------------------------------------------------------------------------------------
16 
51  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int firstFrame,
52  bool onlyTetrahedral) {
53  //
54  // Ring IDs of each type will be saved in these vectors
55  std::vector<int> listDDC; // Vector for ring indices of DDC
56  std::vector<int> listHC; // Vector for ring indices of HC
58  listMixed; // Vector for ring indices of rings that are both DDC and HC
60  ringType; // This vector will have a value for each ring inside
61  // ringList, of type enum strucType in gen.hpp
62  // Make a list of all the DDCs and HCs
63  std::vector<cage::Cage> cageList;
65  ringsOneType; // Vector of vectors of rings of a single size
66  int initRingSize; // Todo or not: calculate the PNCs or not
67  int maxRingSize = 6; // DDCs and HCs are for 6-membered rings
69  atomTypes; // This vector will have a value for every atom
70  // Number of types
71  int numHC, numDDC, mixedRings, prismaticRings, basalRings;
72 
73  // TODO: handle shapeMatching
74  bool doShapeMatching = true;
75  // Qualifier for the RMSD per atom:
76  std::vector<double> rmsdPerAtom;
77  //
78 
79  // test
80  // std::cout << "rings"
81  // << "\n";
82  // for (int i = 0; i < rings.size(); i++) {
83  // for (int j = 0; j < rings[i].size(); j++) {
84  // std::cout << rings[i][j] << " ";
85  // }
86  // std::cout << "\n";
87  // }
88  // std::cout << "\n";
89 
90  if (onlyTetrahedral) {
91  initRingSize = 6;
92  } else {
93  initRingSize = 5;
94  }
95 
96  // Init the atom type vector
97  atomTypes.resize(yCloud->nop); // Has a value for each atom
98  // Init the rmsd per atom (not used yet)
99  rmsdPerAtom.resize(yCloud->nop); // Has a value for each atom
100 
101  // ----------------------------------------------
102  // Init
103  // Get rings of size 5 or 6.
104  for (int ringSize = initRingSize; ringSize <= maxRingSize; ringSize++) {
105  // Clear ringsOneType
106  ring::clearRingList(ringsOneType);
107  // Get rings of the current ring size
108  ringsOneType = ring::getSingleRingSize(rings, ringSize);
109  // Skip for zero rings
110  if (ringsOneType.size() == 0) {
111  continue;
112  } // skip for no rings of ringSize
113  //
114  // Init the ringType vector
115  ringType.resize(
116  ringsOneType.size()); // Has a value for each ring. init to zero.
117  // ----------------------------------------------
118  if (ringSize == 6) {
119  // Get the cages
120 
121  // Find HC rings, saving the ring IDs (starting from 0) to listHC
122  listHC = ring::findHC(ringsOneType, &ringType, nList, &cageList);
123 
124  // Find DDC rings, saving the IDs to listDDC
125  listDDC = ring::findDDC(ringsOneType, &ringType, listHC, &cageList);
126 
127  // Find rings which are both DDCs and HCs (mixed)
128  // A dummy value of -10 in the listDDC and listHC vectors for mixed rings
129  listMixed =
130  ring::findMixedRings(ringsOneType, &ringType, &listDDC, &listHC);
131 
132  // Get the number of structures (DDCs, HCs, mixed rings, basal rings,
133  // prismatic rings)
134  ring::getStrucNumbers(ringType, cageList, &numHC, &numDDC, &mixedRings,
135  &prismaticRings, &basalRings);
136 
137  // Write out to a file
138  sout::writeTopoBulkData(path, yCloud->currentFrame, numHC, numDDC,
139  mixedRings, basalRings, prismaticRings,
140  firstFrame);
141  // Gets the atom type for every atom, to be used for printing out the ice
142  // types found
143  ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
144  }
145  // Finding prismatic blocks
146  else {
147  // Get the prism block classifications
148  prism3::findBulkPrisms(ringsOneType, &ringType, nList, yCloud,
149  &rmsdPerAtom);
150  // Gets the atom type for every atom, to be used for printing out the ice
151  // types found
152  ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
153  }
154  }
155 
156  // Print out the lammps data file with the bonds
157  sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path);
158  // To output the bonds between dummy atoms, uncomment the following line
159  // sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path, true);
160 
161  return 0;
162 }
163 
201  std::vector<int> listHC,
202  std::vector<cage::Cage> *cageList) {
203  std::vector<int> listDDC;
204  int totalRingNum = rings.size(); // Total number of hexagonal rings
205  std::vector<int> peripheralRings; // Indices which may be peripheral rings
206  bool cond1, cond2, cond3; // Conditions for DDC
207  int jring; // Index for peripheral ring being added
208  std::vector<int> DDCRings; // Indices of rings which constitute a single DDC,
209  // with the equatorial ring first
210  // Vector for checking if a ring is basal, prismatic or peripheral
212  notEquatorial; // true if the ring is prismatic, basal or peripheral.
213  notEquatorial.resize(totalRingNum); // Initialized to false
214 
215  // --------
216  // Set all basal and prismatic rings to true (which are part of an HC)
217  for (int i = 0; i < listHC.size(); i++) {
218  int currentRingIndex = listHC[i];
219  // Set this to true
220  notEquatorial[currentRingIndex] = true;
221  } // end of update of notEquatorial
222  // --------
223 
224  // To search for equatorial rings, loop through all
225  // the hexagonal rings
226  for (int iring = 0; iring < totalRingNum; iring++) {
227  // ------------
228  // Step zero: If the ring has been classified as a basal or prismatic ring
229  // in an HC or is a peripheral ring, then it cannot be the equatiorial ring
230  // in a DDC
231  //
232  if (notEquatorial[iring]) {
233  continue;
234  } // skip for rings which are not equatorial
235  //
236  // ------------
237  // Init
238  peripheralRings.clear();
239  // ------------
240  // Step one: Find all rings which contain each index (m_k) of the equatorial
241  // ring, iring, in at least three other rings
242  cond1 = ring::conditionOneDDC(rings, &peripheralRings, iring);
243  if (cond1 == false) {
244  continue;
245  }
246  // ------------
247  // Step two: For every triplet in iring, there is at least one
248  // hexagonal ring other than iring that passes through the triplet.
249  // The peripheral rings are stored in order of the starting element
250  // of each triplet.
251  cond2 = ring::conditionTwoDDC(rings, &peripheralRings, iring);
252  if (cond2 == false) {
253  continue;
254  }
255  // ------------
256  // Step three: For every triplet in the equatorial ring, there is at least
257  // one hexagonal ring other than iring that passes through the triplet.
258  // Rings corresponding to triplets need not be searched again since
259  // peripheralRings are stored in that order. Rings corresponding to T1, T3,
260  // T5 must have a common element. Similarly rings corresponding to T2, T4,
261  // T6 must have at least one common element. Alternating rings corresponding
262  // to triplets must have at least three common elements
263  cond3 = ring::conditionThreeDDC(rings, &peripheralRings);
264  if (cond3 == false) {
265  continue;
266  }
267  // ------------
268  // If the peripheral rings are duplicates, skip everything
269  sort(peripheralRings.begin(), peripheralRings.end());
270  peripheralRings.erase(
271  unique(peripheralRings.begin(), peripheralRings.end()),
272  peripheralRings.end());
273  // There should be 6 unique peripheral rings
274  if (peripheralRings.size() != 6) {
275  continue;
276  }
277  // ------------
278  // If iring is an equatorial ring, add it to the listDDC vector
279  if ((*ringType)[iring] == ring::unclassified) {
280  (*ringType)[iring] = ring::DDC;
281  listDDC.push_back(iring);
282  }
283  // Add the peripheral ring IDs too
284  for (int j = 0; j < peripheralRings.size(); j++) {
285  jring = peripheralRings[j];
286  if ((*ringType)[jring] == ring::unclassified) {
287  (*ringType)[jring] = ring::DDC;
288  listDDC.push_back(jring);
289  } else if ((*ringType)[jring] == ring::HCbasal) {
290  (*ringType)[jring] = ring::bothBasal;
291  listDDC.push_back(jring);
292  } // end of update
293  // never true
294  else if ((*ringType)[jring] == ring::HCprismatic) {
295  (*ringType)[jring] = ring::bothPrismatic;
296  listDDC.push_back(jring);
297  }
298  //
299  // Update the notEquatorial vector
300  notEquatorial[jring] = true;
301  } // end of update for peripheral rings
302  // Add rings to the cageList vector of struct Cages.
303  DDCRings.clear(); // init
304  DDCRings.push_back(iring); // Add the equatorial ring first
305  DDCRings.insert(std::end(DDCRings), std::begin(peripheralRings),
306  std::end(peripheralRings));
307  (*cageList).push_back({cage::DoubleDiaC, DDCRings});
308  // ------------
309  } // end of loop through all hexagonal rings
310 
311  return listDDC;
312 }
313 
331  std::vector<int> *peripheralRings, int iring) {
332  int index; // Atom ID to be compared
333  int noOfCommonRings =
334  0; // No of rings in which the element to be matched has been found
335  int jElement; // Atom ID being compared to index
336 
337  // Loop through each element of iring for finding matches
338  for (int m = 0; m < 6; m++) {
339  index = rings[iring][m]; // Atom Index to be compared and matched with
340  noOfCommonRings = 0; // init to zero.
341  // Loop through every ring except iring
342  for (int jring = 0; jring < rings.size(); jring++) {
343  if (iring == jring) {
344  continue;
345  } // Skip for iring
346  // -------
347  // Search every element of jring
348  for (int k = 0; k < 6; k++) {
349  jElement = rings[jring][k];
350  if (jElement == index) {
351  noOfCommonRings++;
352  (*peripheralRings).push_back(jring);
353  break;
354  } // if index is found inside jring
355  else {
356  continue;
357  }
358  } // end of loop through every element of jring
359  // -------
360  } // end of loop through all rings except iring
361  // If less than 3 rings have been found for each element, then this is
362  // not an equatorial ring
363  if (noOfCommonRings < 3) {
364  return false;
365  } // end of check for common ring number per element in iring
366  } // end of loop through each element of iring
367 
368  // iring is an equatorial ring. The duplicate ring IDs inside
369  // peripheralRings should be removed
370  std::vector<int>::iterator ip; // Iterator to find the last element upto
371  // which unique elements are present
372  // Duplicate IDs must be removed
373  int numElements =
374  (*peripheralRings).size(); // number of elements in peripheralRings
375  sort((*peripheralRings).begin(), (*peripheralRings).end());
376  ip = std::unique((*peripheralRings).begin(),
377  (*peripheralRings).begin() + numElements);
378  // Resize peripheral rings to remove undefined terms
379  (*peripheralRings).resize(std::distance((*peripheralRings).begin(), ip));
380 
381  return true;
382 }
383 
403  std::vector<int> *peripheralRings, int iring) {
404  std::vector<int> triplet; // Triplet formed from iring
405  int ringSize = 6; // Here, all the rings are hexagons
406  int j; // Used for making the triplet
407  int jring; // Peripheral ring ID to be searched
408  int count; // Number of rings found that match the triplet
410  newPeripherals; // Vector in which the new peripheral ring IDs are saved.
411  // This will be swapped with peripheralRings later
412 
413  for (int k = 0; k < ringSize; k++) {
414  triplet.clear(); // Clear the triplet
415  // Get a triplet
416  for (int i = k; i < k + 3; i++) {
417  j = i;
418  if (i >= ringSize) {
419  j = i - ringSize;
420  }
421  triplet.push_back(rings[iring][j]);
422  } // end of getting a triplet from k
423  // -------------
424  // Compare the triplet with every possible peripheral
425  // ring inside peripheralRings.
426  count = 0; // init to zero
427  // Loop through all possible peripheral rings
428  for (int m = 0; m < (*peripheralRings).size(); m++) {
429  jring = (*peripheralRings)[m]; // Ring ID of ring to be searched
430  // Search inside the ring with index jring for the triplet
431  bool foundTriplet = ring::findTripletInRing(rings[jring], triplet);
432 
433  // If the triplet has been found inside jring
434  if (foundTriplet) {
435  newPeripherals.push_back(jring); // Update new peripheral vector
436  count++;
437  break;
438  } // end of ring found
439  } // end of loop through all possible peripheral rings
440  // If count is 0, then the triplet was not found in any peripheral ring
441  if (count == 0) {
442  return false;
443  } // Return false since the triplet was not found
444  // -------------
445  } // end of looping through 0-6 to get triplets
446 
447  // Swap the old peripheral rings vector with the new one
448  (*peripheralRings).swap(newPeripherals);
449 
450  // If there are more than 6 peripheral rings, the code will fail
451  // Comment this out if you want
452  if ((*peripheralRings).size() > 6) {
453  std::cerr
454  << "There are more than 6 peripheral rings. The code will fail. \n";
455  return false;
456  } // end of check for more than 6 peripherals
457 
458  return true;
459 }
460 
484  std::vector<int> *peripheralRings) {
485  // New
486  std::vector<int> common; // Vector containing common elements
487  bool hasCommon; // true if the rings have a common element
488  int iring, jring; // Pairs of peripheral rings
489  // ----------------------------------------------------------------------------
490  // CONDITION 1: Rings corresponding to T1, T3, T5 should have at least one
491  // common element.
492  hasCommon = ring::commonElementsInThreeRings(rings[(*peripheralRings)[0]],
493  rings[(*peripheralRings)[2]],
494  rings[(*peripheralRings)[4]]);
495 
496  // If T1, T3, T5 don't have a common element, return false
497  if (!hasCommon) {
498  return false;
499  } // not a DDC
500  // ----------------------------------------------------------------------------
501  // CONDITION 2: Rings corresponding to T1, T3, T5 should have at least one
502  // common element.
503  hasCommon = ring::commonElementsInThreeRings(rings[(*peripheralRings)[1]],
504  rings[(*peripheralRings)[3]],
505  rings[(*peripheralRings)[5]]);
506 
507  // If T1, T3, T5 don't have a common element, return false
508  if (!hasCommon) {
509  return false;
510  } // not a DDC
511  // ----------------------------------------------------------------------------
512  // CONDITION 3: Rings corresponding to {T1, T3}, {T2, T4}, {T3, T5}, {T4,
513  // T6}
514  // must have three elements in common amongst them
515 
516  // Loops to get pairs of rings corresponding to the right triplets
517  for (int i = 0; i <= 3; i++) {
518  common.clear(); // init
519  // Pairs of rings corresponding to triplets.
520  iring = (*peripheralRings)[i];
521  jring = (*peripheralRings)[i + 2];
522  // Get the common elements
523  common = ring::findsCommonElements(rings[iring], rings[jring]);
524  // There should be at least three elements
525  if (common.size() < 3) {
526  return false;
527  } // not a DDC
528  } // end of getting iring and jring
529  // ----------------------------------------------------------------------------
530 
531  // iring is an equatorial ring and peripheralRings has the 6 peripheral rings
532  return true;
533 }
534 
568  std::vector<cage::Cage> *cageList) {
569  std::vector<int> listHC;
570  int totalRingNum = rings.size(); // Total number of hexagonal rings
571  std::vector<int> basal1; // First basal ring
572  std::vector<int> basal2; // Second basal ring
573  bool cond1, cond2; // Conditions for rings to be basal (true) or not (false)
575  HCRings; // Indices of rings which constitute a single HC, with the basal
576  // rings first, followed by prismatic rings
577  std::vector<int> prismaticRings; // Ring indices of prismatic rings
578  int kring; // Ring index of the prismatic rings
580  isPrismatic; // Flag for checking if the ring is prismatic (true) or not
581  // (false), since the basal rings are checked
582  isPrismatic.resize(totalRingNum); // Initialized to false
583 
584  // Two loops through all the rings are required to find pairs of basal rings
585  for (int iring = 0; iring < totalRingNum - 1; iring++) {
586  // -----------------------
587  // Skip if iring is prismatic
588  if (isPrismatic[iring]) {
589  continue;
590  } // Skip if prismatic
591  // -----------------------
592  cond1 = false;
593  cond2 = false;
594  basal1 = rings[iring]; // Assign iring to basal1
595  // Loop through the other rings to get a pair
596  for (int jring = iring + 1; jring < totalRingNum; jring++) {
597  // -----------------------
598  // Skip if iring is prismatic
599  if (isPrismatic[jring]) {
600  continue;
601  } // Skip if prismatic
602  // -----------------------
603  basal2 = rings[jring]; // Assign jring to basal2
604  // ------------
605  // Step one: Check to see if basal1 and basal2 have common
606  // elements or not. If they don't, then they cannot be basal rings
607  cond1 = ring::hasCommonElements(basal1, basal2);
608  if (cond1 == true) {
609  continue;
610  }
611  // -----------
612  // Step two and three: One of the elements of basal2 must be the nearest
613  // neighbour of either the first (index0; l1) or second (index1; l2)
614  // element of basal1. If m_k is the nearest neighbour of l1, m_{k+2} and
615  // m_{k+4} must be neighbours of l3 and l5(l5 or l3). Modify for l2.
616  cond2 = ring::basalConditions(nList, &basal1, &basal2);
617  if (cond2 == false) {
618  continue;
619  }
620  // -----------
621  // iring and jring are basal rings!
622  // Update iring
623  if ((*ringType)[iring] == ring::unclassified) {
624  (*ringType)[iring] = ring::HCbasal;
625  listHC.push_back(iring);
626  } else if ((*ringType)[iring] == ring::DDC) {
627  (*ringType)[iring] = ring::bothBasal;
628  listHC.push_back(iring);
629  }
630  // Update jring
631  if ((*ringType)[jring] == ring::unclassified) {
632  (*ringType)[jring] = ring::HCbasal;
633  listHC.push_back(jring);
634  } else if ((*ringType)[jring] == ring::DDC) {
635  (*ringType)[jring] = ring::bothBasal;
636  listHC.push_back(jring);
637  }
638  // Find the prismatic rings
639  prismaticRings.clear(); // Clear the prismatic ring vector first
640  ring::findPrismatic(rings, &listHC, ringType, iring, jring,
641  &prismaticRings);
642  // Update the prismatic rings
643  for (int k = 0; k < prismaticRings.size(); k++) {
644  kring =
645  prismaticRings[k]; // Current ring index of the (3) prismatic rings
646  // Update kring
647  if ((*ringType)[kring] == ring::unclassified) {
648  (*ringType)[kring] = ring::HCprismatic;
649  listHC.push_back(kring);
650  } else if ((*ringType)[kring] == ring::DDC) {
651  (*ringType)[kring] = ring::bothPrismatic;
652  listHC.push_back(kring);
653  }
654  //
655  // Update the isPrismatic vector
656  isPrismatic[kring] = true;
657  //
658  } // End of update of prismatic rings in listHC
659  // -----------
660  // Update the cageList vector of Cages
661  // Update the basal rings
662  HCRings.clear();
663  HCRings.push_back(iring);
664  HCRings.push_back(jring);
665  // Add the prismaticRings
666  HCRings.insert(std::end(HCRings), std::begin(prismaticRings),
667  std::end(prismaticRings));
668  (*cageList).push_back({cage::HexC, HCRings});
669  // -----------
670  } // end of loop through rest of the rings to get the second basal ring
671  } // end of loop through all rings for first basal ring
672 
673  sort(listHC.begin(), listHC.end());
674  auto ip = std::unique(listHC.begin(), listHC.end());
675  // Resize peripheral rings to remove undefined terms
676  listHC.resize(std::distance(listHC.begin(), ip));
677 
678  return listHC;
679 }
680 
701  std::vector<int> *basal1, std::vector<int> *basal2) {
702  int l1 = (*basal1)[0]; // first element of basal1 ring
703  int l2 = (*basal1)[1]; // second element of basal1 ring
704  int ringSize = 6; // Size of the ring; each ring contains 6 elements
705  int m_k; // Atom Index (in pointCloud) of element in basal2
706  int kIndex; // Index of m_k in basal2, corresponding to m_k
707  int currentKindex; // Current k index when finding alternating elements of
708  // basal2
709  std::vector<int> evenTriplet; // contains m_k, m_{k+2}, m_{k+4}
710  std::vector<int> oddTriplet; // contains m_{k+1}, m_{k+3}, m_{k+5}
711  int compare1, compare2; // l3 and l5 OR l4 and l6
712  int index;
713  bool l1_neighbour, l2_neighbour; // m_k is a neighbour of l1(true) or not
714  // (false); m_k is a neighbour of l2(true)
715  bool isNeigh, notNeigh; // Used to check if the rings are basal or not
716 
717  // ---------------------------------------------
718  // SEARCH FOR L1_NEIGHBOUR OR L2_NEIGHBOUR
719  // Search for whether an element of basal2 is a neighbour of l1 or l2
720  for (int k = 0; k < ringSize; k++) {
721  // init
722  l1_neighbour = false;
723  l2_neighbour = false;
724  m_k = (*basal2)[k];
725 
726  // ---------------
727  // CHECK IF M_K MATCHES L1 NEIGHBOURS
728  auto it1 = std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
729  // If m_k was found in l1's nList
730  if (it1 != nList[l1].end()) {
731  compare1 = (*basal1)[2]; // l3
732  compare2 = (*basal1)[4]; // l5
733  kIndex = k; // Saving the array index of m_k
734  l1_neighbour = true;
735  break;
736  } // m_k found in l1's nList
737  // ---------------
738  // CHECK IF M_K MATCHES L2 NEIGHBOURS
739  auto it2 = std::find(nList[l2].begin() + 1, nList[l2].end(), m_k);
740  // If m_k was found in l1's nList
741  if (it2 != nList[l2].end()) {
742  compare1 = (*basal1)[3]; // l4
743  compare2 = (*basal1)[5]; // l6
744  kIndex = k; // Saving the array index of m_k
745  l2_neighbour = true;
746  break;
747  } // m_k found in l1's nList
748  // ---------------
749  } // End of search for basal2 elements in l1 or l2's nList
750  // ---------------------------------------------
751 
752  // Return false if neither l1 nor l2 have any neighbours
753  // in basal2
754 
755  if (l1_neighbour == false && l2_neighbour == false) {
756  return false;
757  } // basal conditions not fulfilled
758 
759  // Get the alternating elements starting with kIndex.
760  // 'evenTriplet': m_k, m_{k+2}, m_{k+4} - neighbours of compare1 and compare2.
761  // 'oddTriplet': m_{k+1}, m_{k+3}, m_{k+5}- cannot be neighbours of basal1
762  for (int k = 0; k <= 5; k++) {
763  currentKindex = kIndex + k; // k
764  // Wrap-around
765  if (currentKindex >= ringSize) {
766  currentKindex -= ringSize;
767  } // end of wrap-around of k
768  //
769  // Update 'evenTriplet'
770  if (k % 2 == 0) {
771  evenTriplet.push_back((*basal2)[currentKindex]);
772  } // end of update of evenTriplet
773  // Update 'oddTriplet'
774  else {
775  oddTriplet.push_back((*basal2)[currentKindex]);
776  } // end of update of oddTriplet
777  } // End of getting alternating triplets
778 
779  // ---------------------------------------------
780  // CONDITION1: m_{k+2} and m_{k+4} must be bonded to l3 and l5 (if l1 is a
781  // neighbour) or m_{k+2} and m_{k+4} must be bonded to l4 and l6 (if l2 is a
782  // neighbour) Basically, this boils down to checking whether compare1 and
783  // compare2 are in the neighbour lists of the last two elements of evenTriplet
784 
785  isNeigh = ring::basalNeighbours(nList, &evenTriplet, compare1, compare2);
786 
787  // If condition1 is not true, then the candidate
788  // rings are not part of an HC
789  if (!isNeigh) {
790  return false;
791  } // Not an HC
792 
793  // ---------------------------------------------
794  // CONDITION2: m_{k+1}, m_{k+3} and m_{k+5} must NOT be bonded to any element
795  // in basal1.
796  // Basically, this boils down to checking whether the elements of oddTriplet
797  // are in the neighbour lists of all the elements of basal1.
798 
799  // condition 2. This must be true for an HC
800  notNeigh = ring::notNeighboursOfRing(nList, &oddTriplet, basal1);
801 
802  // If condition2 is not true, the the candidate rings
803  // are not part of an HC
804  if (!notNeigh) {
805  return false;
806  } // Not an HC
807 
808  // Otherwise, all the conditions are true and this is an HC
809  return true;
810 
811 } // end of function
812 
823  std::vector<int> *triplet, int atomOne,
824  int atomTwo) {
825  // Search for needles in a haystack :)
826  int needle1 = (*triplet)[1];
827  int needle2 = (*triplet)[2];
828  bool neighbourFound = false;
829  bool neighOne = false,
830  neighTwo = false; // Neighbour for atomOne, or neighbour for atomTwo
831  // ----------------------------
832  // For first element needle1, which must belong to either atomOne's or
833  // atomTwo's neighbour list Search atomOne's neighbours
834  auto it =
835  std::find(nList[atomOne].begin() + 1, nList[atomOne].end(), needle1);
836  if (it != nList[atomOne].end()) {
837  neighbourFound = true;
838  neighOne = true;
839  } // atomOne's neighbour
840  // If it is not atomOne's neighbour, it might be atomTwo's neighbour
841  if (!neighOne) {
842  it = std::find(nList[atomTwo].begin() + 1, nList[atomTwo].end(), needle1);
843  if (it != nList[atomTwo].end()) {
844  neighbourFound = true;
845  neighTwo = true;
846  } // end of check to see if neighbour was found
847  } // End of check to see if needle1 is atomTwo's neighbour
848  // ----------------------------
849  // If needle1 is not a neighbour of atomOne or atomTwo, return false
850  if (neighbourFound == false) {
851  return false;
852  }
853 
854  // Check to see if needle2 is a neighbour of either atomOne or atomTwo
855  // ===============
856  // if atomOne was a neighbour of needle1, needle2 must be a neighbour of
857  // atomTwo
858  if (neighOne) {
859  it = std::find(nList[atomTwo].begin() + 1, nList[atomTwo].end(), needle2);
860  // It is a neighbour
861  if (it != nList[atomTwo].end()) {
862  return true;
863  }
864  // It is not a neighbour
865  else {
866  return false;
867  }
868  } // End of check for neighbour of atomTwo
869  // ===============
870  // if atomTwo was a neighbour of needle1, needle2 must be a neighbour of
871  // atomOne
872  else {
873  it = std::find(nList[atomOne].begin() + 1, nList[atomOne].end(), needle2);
874  // It is a neighbour
875  if (it != nList[atomOne].end()) {
876  return true;
877  }
878  // It is not a neighbour
879  else {
880  return false;
881  }
882  }
883  // ===============
884 }
885 
896  std::vector<int> *triplet,
898  int iatom; // AtomID of the atom to be searched for inside the neighbour
899  // lists
900  int jatom; // AtomID of in whose neighbour list iatom will be searched for
902 
903  for (int i = 0; i < (*triplet).size(); i++) {
904  iatom = (*triplet)[i]; // AtomID to be searched for
905  // iatom must be searched for in the neighbour lists of all elements
906  // of the ring vector
907  for (int j = 0; j < (*ring).size(); j++) {
908  jatom = (*ring)[j];
909  // ------------------
910  // Search for iatom in the neighbour list of jatom
911  it = std::find(nList[jatom].begin() + 1, nList[jatom].end(), iatom);
912  // It is a neighbour!
913  if (it != nList[jatom].end()) {
914  return false;
915  }
916  // ------------------
917  } // end of loop through every element of ring
918  } // end of loop through all triplet elements
919 
920  return true;
921 }
922 
936  std::vector<int> *listHC,
937  std::vector<ring::strucType> *ringType, int iring,
938  int jring, std::vector<int> *prismaticRings) {
939  int iIndex, jIndex; // Used for making rings to be searched
940  int ringSize = rings[0].size(); // This is 6 for hexagons
941  std::vector<int> iTriplet; // triplet formed from iring
942  std::vector<int> jTriplet; // triplet formed from jring
943  std::vector<int> common; // Common elements
944 
945  // Make all possible triplets out of iring
946  for (int i = 0; i < ringSize; i++) {
947  // init
948  iTriplet.clear();
949  // ------
950  // Get a triplet
951  for (int m = 0; m < 3; m++) {
952  iIndex = i + m;
953  if (iIndex >= ringSize) {
954  iIndex -= ringSize;
955  }
956  iTriplet.push_back(rings[iring][iIndex]);
957  } // end of getting a triplet from iring
958 
959  // -------------------------------------------
960  // Now that a triplet has been found, find all rings with that triplet in
961  // it!
962  for (int kring = 0; kring < rings.size(); kring++) {
963  // If this is the same as iring or jring, skip
964  if (kring == iring || kring == jring) {
965  continue;
966  } // is not prismatic
967  //
968  // Now find out whether kring has the triplet or not!
969  common = ring::findsCommonElements(iTriplet, rings[kring]);
970 
971  // If this triplet is not shared by kring
972  // skip
973  if (common.size() != 3) {
974  continue;
975  } // kring does not have iTriplet
976 
977  // -----------------
978  // If kring does have the triplet, check to see if at least three other
979  // elements of kring are shared by jring
980  jTriplet.clear(); // init
981  // Make jTriplet
982  for (int j = 0; j < ringSize; j++) {
983  int katom = rings[kring][j];
984  auto it = std::find(iTriplet.begin(), iTriplet.end(), katom);
985 
986  // If not found, add it to jTriplet
987  if (it == iTriplet.end()) {
988  jTriplet.push_back(katom);
989  } // update jTriplet
990  } // end of making jTriplet out of kring
991  // -----------------
992 
993  // Now search for jTriplet inside jring
994  common = ring::findsCommonElements(jTriplet, rings[jring]);
995 
996  // Update the prismatic rings
997  if (common.size() == 3) {
998  //
999  (*listHC).push_back(kring); // Update listHC vector
1000  (*prismaticRings).push_back(kring); // Update prismatic rings
1001  // Update the type inside ringType
1002  // If the ring is already a DDC ring, it is a mixed ring
1003  if ((*ringType)[kring] == ring::DDC) {
1004  (*ringType)[kring] = ring::bothPrismatic;
1005  }
1006  // If it is unclassified, it is just a prismatic ring
1007  if ((*ringType)[kring] == ring::unclassified) {
1008  (*ringType)[kring] = ring::HCprismatic;
1009  } // end ring update
1010  } // add kring to the list of prismatic rings
1011  } // end of searching through rings for kring
1012  // -------------------------------------------
1013  } // end of getting pairs of triplets
1014 
1015  return 0;
1016 }
1017 
1032  std::vector<strucType> *ringType,
1033  std::vector<int> *listDDC,
1034  std::vector<int> *listHC) {
1035  std::vector<int> listMixed;
1036  int dummyValue = -10;
1037 
1038  // Loop through all rings in the ringType and
1039  // adds the ring Indices of all rings which are both DDCs and HCs
1040  for (int iring = 0; iring < (*ringType).size(); iring++) {
1041  // If iring is of mixed type, add it to the listMixed vector
1042  if ((*ringType)[iring] == ring::bothBasal ||
1043  (*ringType)[iring] == ring::bothPrismatic) {
1044  listMixed.push_back(iring);
1045 
1046  //-----------------
1047  // Search for iring in listDDC
1048  std::sort((*listDDC).begin(), (*listDDC).end());
1049  auto iter = std::find((*listDDC).begin(), (*listDDC).end(), iring);
1050  if (iter != (*listDDC).end()) {
1051  *iter = dummyValue; // Assign dummy value to the mixed ring
1052  } // found in listDDC
1053  //-----------------
1054  //-----------------
1055  // Search for iring in listHC
1056  std::sort((*listHC).begin(), (*listHC).end());
1057  auto itr = std::find((*listHC).begin(), (*listHC).end(), iring);
1058  if (itr != (*listHC).end()) {
1059  *itr = dummyValue; // Assign dummy value to the mixed ring
1060  } // found in listHC
1061  //-----------------
1062 
1063  } // end of check for type
1064  } // end of loop through all every ring
1065 
1066  return listMixed;
1067 }
1068 
1080  std::vector<cage::iceType> *atomTypes) {
1081  cage::iceType iRingType; // Type of the current ring
1082  int iatom; // Current ring
1083  int ringSize = rings[0].size(); // Size of the ring
1084 
1085  // Loop through every ring in ringType
1086  for (int iring = 0; iring < ringType.size(); iring++) {
1087  //
1088  // Skip if the ring is unclassified
1089  if (ringType[iring] == ring::unclassified) {
1090  continue;
1091  } // skip for unclassified rings
1092 
1093  // ------------
1094  // Get the current ring type
1095  // DDC
1096  if (ringType[iring] == ring::DDC) {
1097  iRingType = cage::ddc;
1098  } // DDC atoms
1099  //
1100  // HC
1101  else if (ringType[iring] == ring::HCbasal ||
1102  ringType[iring] == ring::HCprismatic) {
1103  iRingType = cage::hc;
1104  } // HC atoms
1105  //
1106  // Mixed
1107  else if (ringType[iring] == ring::bothBasal ||
1108  ringType[iring] == ring::bothPrismatic) {
1109  iRingType = cage::mixed;
1110  } // HC atoms
1111  // Prism
1112  else if (ringType[iring] == ring::Prism ||
1113  ringType[iring] == ring::deformedPrism ||
1114  ringType[iring] == ring::mixedPrismRing) {
1115  iRingType = cage::pnc; // 5 membered pnc
1116  } // prism
1117  // Should never go here
1118  else {
1119  continue;
1120  } //
1121  // ------------
1122  // Otherwise, loop through every inside the ring and assign atomTypes the
1123  // iRingType
1124  for (int i = 0; i < ringSize; i++) {
1125  iatom = rings[iring][i]; // Atom index in ring
1126  if ((*atomTypes)[iatom] == cage::mixed ||
1127  (*atomTypes)[iatom] == cage::mixed2) {
1128  continue;
1129  } // Don't reassign
1130  // For atoms shared by PNCs and DDCs/HCs
1131  if (ringSize == 6) {
1132  if ((*atomTypes)[iatom] == cage::pnc) {
1133  (*atomTypes)[iatom] = cage::mixed2;
1134  } else {
1135  (*atomTypes)[iatom] = iRingType;
1136  }
1137  } else {
1138  (*atomTypes)[iatom] = iRingType;
1139  }
1140  } // end of loop thorugh the current ring
1141  } // end of loop through every ring
1142 
1143  return 0;
1144 }
1145 
1162  std::vector<cage::Cage> cageList, int *numHC,
1163  int *numDDC, int *mixedRings, int *prismaticRings,
1164  int *basalRings) {
1165  // Determines the number of HCs, DDCs, Mixed rings, prismatic and basal rings
1166  // Init
1167  *numHC = 0;
1168  *numDDC = 0;
1169  *mixedRings = 0;
1170  *prismaticRings = 0;
1171  *basalRings = 0;
1172  // ------------------------------------
1173  // GETTING THE CAGES (DDCs and HCs)
1174  // Loop through cages
1175  for (int icage = 0; icage < cageList.size(); icage++) {
1176  //
1177  // HC
1178  if (cageList[icage].type == cage::HexC) {
1179  *numHC += 1;
1180  } // end of updating HC number
1181  //
1182  // DDC
1183  if (cageList[icage].type == cage::DoubleDiaC) {
1184  *numDDC += 1;
1185  } // end of updating DDC number
1186  } // end of loop through cages
1187  // ------------------------------------
1188  // GETTING THE RINGSS (Mixed, Prismatic and Basal rings)
1189  // Loop through the rings
1190  for (int iring = 0; iring < ringType.size(); iring++) {
1191  // Mixed
1192  if (ringType[iring] == ring::bothBasal ||
1193  ringType[iring] == ring::bothPrismatic) {
1194  *mixedRings += 1;
1195  // Also update basal rings
1196  if (ringType[iring] == ring::bothBasal) {
1197  *basalRings += 1;
1198  } // mixed basal rings
1199  // Also update prismatic rings
1200  if (ringType[iring] == ring::bothPrismatic) {
1201  *prismaticRings += 1;
1202  } // mixed prismatic rings
1203  } // end of updating mixed
1204  //
1205  // HCs
1206  if (ringType[iring] == ring::HCprismatic) {
1207  *prismaticRings += 1;
1208  } // HC prismatic
1209  // basal HCs
1210  if (ringType[iring] == ring::HCbasal) {
1211  *basalRings += 1;
1212  } // HC basal
1213  } // end of loop through every ring
1214  // ------------------------------------
1215 
1216  return 0;
1217 } // end of function
1218 
1245  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1246  std::vector<double> *rmsdPerAtom, double heightCutoff) {
1247  int totalRingNum = rings.size(); // Total number of rings
1248  std::vector<int> basal1; // First basal ring
1249  std::vector<int> basal2; // Second basal ring
1250  bool cond1, cond2; // Conditions for rings to be basal (true) or not (false)
1251  int ringSize = rings[0].size(); // Number of nodes in each ring
1252  // Matrix for the reference ring for a given ringSize.
1253  Eigen::MatrixXd refPointSet(ringSize, 3);
1254 
1255  int axialDim = 2; // Default=z
1256  refPointSet = pntToPnt::getPointSetRefRing(ringSize, axialDim);
1257  //
1258 
1259  // Two loops through all the rings are required to find pairs of basal rings
1260  for (int iring = 0; iring < totalRingNum - 1; iring++) {
1261  cond1 = false;
1262  cond2 = false;
1263  basal1 = rings[iring]; // Assign iring to basal1
1264  // Loop through the other rings to get a pair
1265  for (int jring = iring + 1; jring < totalRingNum; jring++) {
1266  basal2 = rings[jring]; // Assign jring to basal2
1267  // ------------
1268  // Put extra check for axial basal rings if shapeMatching is being done
1269  // ------------
1270  // Step one: Check to see if basal1 and basal2 have common
1271  // elements or not. If they don't, then they cannot be basal rings
1272  cond1 = ring::hasCommonElements(basal1, basal2);
1273  if (cond1 == true) {
1274  continue;
1275  }
1276 
1277  // ------------
1278  bool smallDist =
1279  prism3::basalRingsSeparation(yCloud, basal1, basal2, heightCutoff);
1280  if (!smallDist) {
1281  continue;
1282  } // the basal rings are too far apart
1283 
1284  // Otherwise
1285  // Do shape matching here
1286  bool isPrism = match::matchUntetheredPrism(yCloud, nList, refPointSet,
1287  &basal1, &basal2, rmsdPerAtom);
1288 
1289  // Success! The rings are basal rings of a prism!
1290  if (isPrism) {
1291  //
1292  // Update iring
1293  if ((*ringType)[iring] == ring::unclassified) {
1294  (*ringType)[iring] = ring::Prism;
1295  }
1296  // Update jring
1297  if ((*ringType)[jring] == ring::unclassified) {
1298  (*ringType)[jring] = ring::Prism;
1299  }
1300  } // end of reduced criteria
1301  // Strict criteria
1302  else {
1303  cond2 = prism3::basalPrismConditions(nList, &basal1, &basal2);
1304  // If the condition is false then the strict criterion has not been met
1305  if (!cond2) {
1306  continue;
1307  }
1308  // Update iring
1309  if ((*ringType)[iring] == ring::unclassified) {
1310  (*ringType)[iring] = ring::Prism;
1311  }
1312  // Update jring
1313  if ((*ringType)[jring] == ring::unclassified) {
1314  (*ringType)[jring] = ring::Prism;
1315  }
1316  //
1317  // Shape-matching to get the RMSD (if shape-matching is desired)
1318 
1319  // bool isKnownPrism = match::matchPrism(
1320  // yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, true);
1321 
1322  // -----------
1323  } // end of strict criteria
1324 
1325  } // end of loop through rest of the rings to get the second basal ring
1326  } // end of loop through all rings for first basal ring
1327 
1328  return 0;
1329 }
1330 
1344  std::vector<int> *basal1,
1345  std::vector<int> *basal2) {
1346  int l1 = (*basal1)[0]; // first element of basal1 ring
1347  int ringSize =
1348  (*basal1).size(); // Size of the ring; each ring contains n elements
1349  int m_k; // Atom ID of element in basal2
1350  bool l1_neighbour; // m_k is a neighbour of l1(true) or not (false)
1351 
1352  // isNeighbour is initialized to false for all basal2 elements; indication if
1353  // basal2 elements are neighbours of basal1
1354  std::vector<bool> isNeighbour(ringSize, false);
1355  int kIndex; // m_k index
1356  int lAtomID; // atomID of the current element of basal1
1357  int kAtomID; // atomID of the current element of basal2
1358 
1359  // ---------------------------------------------
1360  // COMPARISON OF basal2 ELEMENTS WITH l1
1361  for (int k = 0; k < ringSize; k++) {
1362  l1_neighbour = false;
1363  m_k = (*basal2)[k];
1364  // =================================
1365  // Checking to seee if m_k is be a neighbour of l1
1366  // Find m_k inside l1 neighbour list
1367  auto it = std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
1368 
1369  // If the element has been found, for l1
1370  if (it != nList[l1].end()) {
1371  l1_neighbour = true;
1372  kIndex = k;
1373  break;
1374  }
1375  } // l1 is a neighbour of m_k
1376 
1377  // If there is no nearest neighbour, then the two rings are not part of the
1378  // prism
1379  if (!l1_neighbour) {
1380  return false;
1381  }
1382 
1383  // ---------------------------------------------
1384  // NEIGHBOURS of basal1 in basal2
1385  isNeighbour[kIndex] = true;
1386 
1387  // All elements of basal1 must be neighbours of basal2
1388  for (int i = 1; i < ringSize; i++) {
1389  lAtomID = (*basal1)[i]; // element of basal1 ring
1390  for (int k = 0; k < ringSize; k++) {
1391  // Skip if already a neighbour
1392  if (isNeighbour[k]) {
1393  continue;
1394  }
1395  // Get the comparison basal2 element
1396  kAtomID = (*basal2)[k]; // element of basal2 ring;
1397 
1398  // Checking to see if kAtomID is a neighbour of lAtomID
1399  // Find kAtomID inside lAtomID neighbour list
1400  auto it1 =
1401  std::find(nList[lAtomID].begin() + 1, nList[lAtomID].end(), kAtomID);
1402 
1403  // If the element has been found, for l1
1404  if (it1 != nList[lAtomID].end()) {
1405  isNeighbour[k] = true;
1406  }
1407  } // Loop through basal2
1408  } // end of check for neighbours of basal1
1409 
1410  // ---------------------------------------------
1411 
1412  // They should all be neighbours
1413  for (int k = 0; k < ringSize; k++) {
1414  // Check to see if any element is false
1415  if (!isNeighbour[k]) {
1416  return false;
1417  }
1418  }
1419 
1420  // Everything works out!
1421  return true;
1422 }
1423 
1431  std::vector<int> *basal1,
1432  std::vector<int> *basal2) {
1433  int ringSize =
1434  (*basal1).size(); // Size of the ring; each ring contains n elements
1435  int m_k; // Atom ID of element in basal2
1436  bool isNeighbour = false; // This is true if there is at least one bond
1437  // between the basal rings
1438  int l_k; // Atom ID of element in basal1
1439 
1440  // ---------------------------------------------
1441  // COMPARISON OF basal2 ELEMENTS (m_k) WITH basal1 ELEMENTS (l_k)
1442  // Loop through all the elements of basal1
1443  for (int l = 0; l < ringSize; l++) {
1444  l_k = (*basal1)[l];
1445  // Search for the nearest neighbour of l_k in basal2
1446  // Loop through basal2 elements
1447  for (int m = 0; m < ringSize; m++) {
1448  m_k = (*basal2)[m];
1449  // Find m_k inside l_k neighbour list
1450  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
1451 
1452  // If the element has been found, for l1
1453  if (it != nList[l_k].end()) {
1454  isNeighbour = true;
1455  break;
1456  } // found element
1457  } // end of loop through all the elements of basal2
1458 
1459  // If a neighbour has been found then
1460  if (isNeighbour) {
1461  return true;
1462  }
1463  } // end of loop through all the elements of basal1
1464 
1465  // If a neighbour has not been found, return false
1466  return false;
1467 }
1468 
1472  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1473  std::vector<int> basal1, std::vector<int> basal2, double heightCutoff) {
1474  //
1475  int ringSize = basal1.size();
1476  int l_k, m_k; // Atom indices
1477  double infHugeNumber = 100000;
1478  double leastDist = infHugeNumber;
1479  int index = -1; // starting index
1480  // For the first element of basal1:
1481 
1482  l_k = basal1[0]; // This is the atom particle C++ index
1483 
1484  // Search for the nearest neighbour of l_k in basal2
1485  // Loop through basal2 elements
1486  for (int m = 0; m < ringSize; m++) {
1487  m_k = basal2[m]; // Atom index to find in the neighbour list of iatom
1488 
1489  // Calculate the distance
1490  double dist = gen::periodicDist(yCloud, l_k, m_k);
1491 
1492  // Update the least distance
1493  if (leastDist > dist) {
1494  leastDist = dist; // This is the new least distance
1495  index = m;
1496  } // end of update of the least distance
1497 
1498  } // found element
1499 
1500  // If the element has been found, for l1
1501  if (leastDist < heightCutoff) {
1502  return true;
1503  } // end of check
1504  else {
1505  return false;
1506  }
1507 }
T begin(T... args)
T clear(T... args)
T distance(T... args)
T end(T... args)
T erase(T... args)
T find(T... args)
iceType
Definition: cage.hpp:71
@ HexC
Definition: cage.hpp:50
@ DoubleDiaC
Definition: cage.hpp:50
@ pnc
Definition: cage.hpp:71
@ ddc
Definition: cage.hpp:71
@ mixed
Definition: cage.hpp:71
@ hc
Definition: cage.hpp:71
@ mixed2
Definition: cage.hpp:71
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
bool notNeighboursOfRing(std::vector< std::vector< int >> nList, std::vector< int > *triplet, std::vector< int > *ring)
Definition: topo_bulk.cpp:895
bool conditionOneDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
Definition: topo_bulk.cpp:330
bool basalNeighbours(std::vector< std::vector< int >> nList, std::vector< int > *triplet, int atomOne, int atomTwo)
Definition: topo_bulk.cpp:822
bool relaxedPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
Definition: topo_bulk.cpp:1430
std::vector< int > findDDC(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > listHC, std::vector< cage::Cage > *cageList)
Definition: topo_bulk.cpp:199
bool conditionTwoDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
Definition: topo_bulk.cpp:402
int getAtomTypesTopoBulk(std::vector< std::vector< int >> rings, std::vector< ring::strucType > ringType, std::vector< cage::iceType > *atomTypes)
Definition: topo_bulk.cpp:1078
int getStrucNumbers(std::vector< ring::strucType > ringType, std::vector< cage::Cage > cageList, int *numHC, int *numDDC, int *mixedRings, int *prismaticRings, int *basalRings)
Determines the number of HCs, DDCs, Mixed rings, prismatic and basal rings.
Definition: topo_bulk.cpp:1161
bool basalRingsSeparation(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, double heightCutoff=8)
Check to see that candidate basal prisms are not really far from each other.
Definition: topo_bulk.cpp:1471
bool basalConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
Tests whether two rings are basal rings (true) or not (false)
Definition: topo_bulk.cpp:700
int findPrismatic(std::vector< std::vector< int >> rings, std::vector< int > *listHC, std::vector< strucType > *ringType, int iring, int jring, std::vector< int > *prismaticRings)
Finds the prismatic rings from basal rings iring and jring.
Definition: topo_bulk.cpp:935
int topoBulkAnalysis(std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, bool onlyTetrahedral=true)
Definition: topo_bulk.cpp:48
bool conditionThreeDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings)
Definition: topo_bulk.cpp:483
int findBulkPrisms(std::vector< std::vector< int >> rings, std::vector< ring::strucType > *ringType, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, double heightCutoff=8)
Find out which rings are prisms.
Definition: topo_bulk.cpp:1242
std::vector< int > findMixedRings(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > *listDDC, std::vector< int > *listHC)
Definition: topo_bulk.cpp:1031
bool basalPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
Definition: topo_bulk.cpp:1343
std::vector< int > findHC(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< std::vector< int >> nList, std::vector< cage::Cage > *cageList)
Definition: topo_bulk.cpp:565
bool commonElementsInThreeRings(std::vector< int > ring1, std::vector< int > ring2, std::vector< int > ring3)
Common elements in 3 rings.
Definition: ring.cpp:63
bool findTripletInRing(std::vector< int > ring, std::vector< int > triplet)
Searches a particular ring for a triplet.
Definition: ring.cpp:97
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
Definition: ring.cpp:18
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Definition: ring.cpp:175
std::vector< std::vector< int > > getSingleRingSize(std::vector< std::vector< int >> rings, int ringSize)
Returns a vector of vectors of rings of a single size.
Definition: ring.cpp:149
std::vector< int > findsCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Returns the common elements of two rings.
Definition: ring.cpp:34
@ bothBasal
The ring belongs to both a DDC and HC. It is a 'mixed' ring. The ring is also one of the basal rings ...
Definition: ring.hpp:117
@ DDC
The ring belongs to a double-diamond cage (DDC).
Definition: ring.hpp:114
@ 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
@ deformedPrism
Definition: ring.hpp:120
@ unclassified
The ring is unclassified, which may be either water or a deformed type which cannot be classified by ...
Definition: ring.hpp:113
@ bothPrismatic
The ring belongs to both a DDC and HC and is, thus, a 'mixed' ring. The ring is also one of the prism...
Definition: ring.hpp:118
@ HCprismatic
The ring belongs only to a hexagonal cage (HC); specifically the ring is purely a prismatic ring of a...
Definition: ring.hpp:116
@ mixedPrismRing
Definition: ring.hpp:121
@ Prism
The ring belongs to a prism block, classified according to the prism identification scheme.
Definition: ring.hpp:119
T insert(T... args)
bool matchUntetheredPrism(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, std::vector< double > *rmsdPerAtom)
Definition: shapeMatch.cpp:106
Eigen::MatrixXd getPointSetRefRing(int n, int axialDim)
Topological network criteria functions.
Definition: ring.hpp:60
int writeTopoBulkData(std::string path, int currentFrame, int numHC, int numDDC, int mixedRings, int basalRings, int prismaticRings, int firstFrame)
int writeLAMMPSdataTopoBulk(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< cage::iceType > atomTypes, std::string path, bool bondsBetweenDummy=false)
T push_back(T... args)
T resize(T... args)
T size(T... args)
T sort(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
File containing functions used specific to bulk topological network critera.
T unique(T... args)