Prism3
+ Collaboration diagram for Prism3:

Functions

int ring::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)
 
std::vector< int > ring::findDDC (std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > listHC, std::vector< cage::Cage > *cageList)
 
std::vector< int > ring::findMixedRings (std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > *listDDC, std::vector< int > *listHC)
 
std::vector< int > ring::findHC (std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< std::vector< int >> nList, std::vector< cage::Cage > *cageList)
 
bool ring::conditionOneDDC (std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
 
bool ring::conditionTwoDDC (std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
 
bool ring::conditionThreeDDC (std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings)
 
bool ring::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) More...
 
bool ring::basalNeighbours (std::vector< std::vector< int >> nList, std::vector< int > *triplet, int atomOne, int atomTwo)
 
bool ring::notNeighboursOfRing (std::vector< std::vector< int >> nList, std::vector< int > *triplet, std::vector< int > *ring)
 
int ring::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. More...
 
int ring::getAtomTypesTopoBulk (std::vector< std::vector< int >> rings, std::vector< ring::strucType > ringType, std::vector< cage::iceType > *atomTypes)
 
int ring::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. More...
 
int prism3::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. More...
 
bool prism3::basalPrismConditions (std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
 
bool prism3::relaxedPrismConditions (std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
 
bool prism3::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. More...
 

Detailed Description

Function Documentation

◆ basalConditions()

bool ring::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)

Check to see if two basal rings are HCs or not, using the neighbour list information. The neighbour list nList is a vector of vectors, containing atom indices (not atom IDs!). The first element of each subvector in nList is the atom index of the particle for which the other elements are the nearest neighbours. Internally, the following functions are called:

  • ring::basalNeighbours (CONDITION1: \( m_{k+2} \) and \( m_{k+4} \) must be bonded to \( l_3 \) and \( l_5 \) (if \( l_1 \) is a neighbour) or \( m_{k+2} \) and \( m_{k+4} \) must be bonded to \( l_4 \) and \( l_6 \) (if \( l_2 \) is a neighbour)).
  • ring::notNeighboursOfRing (CONDITION2: \( m_{k+1} \), \( m_{k+3} \) and \( m_{k+5}\) must NOT be bonded to any element in basal1).
    Parameters
    [in]nListVector of vectors of the neighbour list (by index).
    [in]basal1Vector containing the first candidate basal ring.
    [in]basal2Vector containing the second candidate basal ring.
    Returns
    A bool; true if the basal rings being tested fulfill this condition for being the basal rings of an HC.

Definition at line 700 of file topo_bulk.cpp.

701  {
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
T begin(T... args)
T end(T... args)
T find(T... args)
bool notNeighboursOfRing(std::vector< std::vector< int >> nList, std::vector< int > *triplet, std::vector< int > *ring)
Definition: topo_bulk.cpp:895
bool basalNeighbours(std::vector< std::vector< int >> nList, std::vector< int > *triplet, int atomOne, int atomTwo)
Definition: topo_bulk.cpp:822
T push_back(T... args)

◆ basalNeighbours()

bool ring::basalNeighbours ( std::vector< std::vector< int >>  nList,
std::vector< int > *  triplet,
int  atomOne,
int  atomTwo 
)

Tests whether the last two elements of a triplet are neighbours of two atom IDs passed in

Tests whether the last two elements of a triplet are neighbours of two atom indices which have been passed in as inputs.

Parameters
[in]nListVector of vectors of the neighbour list (by index).
[in]tripletVector containing the current triplet being tested.
[in]atomOneIndex of the first atom.
[in]atomTwoIndex of the second atom.
Returns
A bool; true if the condition is met and false otherwise.

Definition at line 822 of file topo_bulk.cpp.

824  {
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 }

◆ basalPrismConditions()

bool prism3::basalPrismConditions ( 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) for a prism (strict criterion)

A function that checks to see if two basal rings are basal rings of a prism block or not, using the neighbour list information. The neighbour list nList is a row-ordered vector of vectors, containing atom indices (not atom IDs!). The first element of each subvector in nList is the atom index of the particle for which the other elements are the nearest neighbours.

Parameters
[in]nListRow-ordered neighbour list by atom index.
[in]basal1The vector for one of the basal rings.
[in]basal2The vector for the other basal ring.
Returns
A value that is true if the basal rings constitute a prism block, and false if they do not make up a prism block.

Definition at line 1343 of file topo_bulk.cpp.

1345  {
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 }

◆ basalRingsSeparation()

bool prism3::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.

Check to see that candidate basal prisms are not really far from each other Return true if the basal rings are within the heightCutoff

Definition at line 1471 of file topo_bulk.cpp.

1473  {
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 }
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
T size(T... args)

◆ conditionOneDDC()

bool ring::conditionOneDDC ( std::vector< std::vector< int >>  rings,
std::vector< int > *  peripheralRings,
int  iring 
)

First condition for the DDC: There must be at least 3 other rings in which each element of the equatorial ring is present

For a given ring, which is being tested as the equatorial ring, this function tests if each element of the ring \( (m_k) \) is present in at least three other rings or not. Returns false if this is not true. The ring indices of rings that have the common element are saved inside the periperal ring vector (peripheralRings) as potential peripheral rings, which is then passed as an input to the subsequent functions testing conditions.

Parameters
[in]ringsVector of vectors containing the 6-membered primitive rings.
[in]peripheralRingsVector containing the indices of rings which are potential peripheral rings.
[in]iringIndex of the ring being tested as equatorial.
Returns
A bool; true if ring being tested as equatorial, iring, satisfies the above condition for being an equatorial ring, and false otherwise.

Definition at line 330 of file topo_bulk.cpp.

331  {
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 }
T distance(T... args)
T sort(T... args)
T unique(T... args)

◆ conditionThreeDDC()

bool ring::conditionThreeDDC ( std::vector< std::vector< int >>  rings,
std::vector< int > *  peripheralRings 
)

Third condition for the DDC: Even (by vector index) numbered index triplets and odd triplets must have at least one element in common

For a given ring, which is being tested as the equatorial ring, this function tests the following, given peripheralRings stored in increasing order of the triplet starting element:

  1. Rings corresponding to T1, T3, T5 should have at least one common element.
  2. Rings corresponding to T2, T4, T6 should have at least one common element.
  3. The following rings should have at least three common elements- {T1, T3}, {T2, T4}, {T3, T5}, {T4, T6}. Internally calls the following functions:
  • ring::commonElementsInThreeRings (CONDITION 1: Rings corresponding to T1, T3, T5 should have at least one common element).
  • ring::commonElementsInThreeRings (CONDITION 2: Rings corresponding to T1, T3, T5 should have at least one common element).
  • ring::findsCommonElements (Gets the common elements between a pair of rings).
    Parameters
    [in]ringsVector of vectors containing the 6-membered primitive rings.
    [in]peripheralRingsVector containing the indices of rings which are potential peripheral rings.
    Returns
    A bool; true if ring being tested as equatorial satisfies the above condition for being an equatorial ring, and false otherwise.

Definition at line 483 of file topo_bulk.cpp.

484  {
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 }
T clear(T... args)
bool commonElementsInThreeRings(std::vector< int > ring1, std::vector< int > ring2, std::vector< int > ring3)
Common elements in 3 rings.
Definition: ring.cpp:63
std::vector< int > findsCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Returns the common elements of two rings.
Definition: ring.cpp:34

◆ conditionTwoDDC()

bool ring::conditionTwoDDC ( std::vector< std::vector< int >>  rings,
std::vector< int > *  peripheralRings,
int  iring 
)

Second condition for the DDC: There must be at least 1 other ring for every triplet in the equatorial ring

For a given ring, which is being tested as the equatorial ring, this function tests if each triplet that can be formed from the ring is common to at least one other ring or not. Returns false if this is not true. The ring indices of rings that have the common triplet are ultimately saved inside the periperal ring vector as potential peripheral rings, which is passed as an input to the subsequent condition-testing functions. The function calls the following function internally:

  • ring::findTripletInRing (Searches inside another ring with index jring for the current triplet, which was obtained from iring).
    Parameters
    [in]ringsVector of vectors containing the 6-membered primitive rings.
    [in]peripheralRingsVector containing the indices of rings which are potential peripheral rings.
    [in]iringIndex of the ring being tested as equatorial.
    Returns
    A bool; true if ring being tested as equatorial, iring, satisfies the above condition for being an equatorial ring, and false otherwise.

Definition at line 402 of file topo_bulk.cpp.

403  {
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 }
T count(T... args)
bool findTripletInRing(std::vector< int > ring, std::vector< int > triplet)
Searches a particular ring for a triplet.
Definition: ring.cpp:97

◆ findBulkPrisms()

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

Determines which rings are n-sided prisms. This function returns a vector which contains the ring indices of all the rings which are prisms. The ring indices correspond to the index of the rings inside the vector of vector rings, starting from 0. Prism rings can be found using a three-step procedure, in which first two basal rings are found. Prismatic rings are simply rings which share every face made by upper and lower triplets of the basal rings The neighbour list is also required as an input, which is a vector of vectors, containing atom IDs. The first element of the neighbour list is the atom index of the atom for which the other elements are nearest neighbours.\

Parameters
[in]ringsThe input vector of vectors containing the primitive rings of a single ring size (number of nodes).
[in]ringTypeA vector containing a ring::strucType value (a classification type) for each ring.
[in]nPrismsThe number of prism blocks identified for the number of nodes.
[in]nListThe row-ordered neighbour list (by atom index).
[in]yCloudThe input PointCloud.
Returns
A vector containing the ring indices of all the rings which have been classified as prisms. The indices are with respect to the input rings vector of vectors.

Definition at line 1242 of file topo_bulk.cpp.

1246  {
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 }
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 basalPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
Definition: topo_bulk.cpp:1343
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Definition: ring.cpp:175
@ unclassified
The ring is unclassified, which may be either water or a deformed type which cannot be classified by ...
Definition: ring.hpp:113
@ Prism
The ring belongs to a prism block, classified according to the prism identification scheme.
Definition: ring.hpp:119
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)

◆ findDDC()

std::vector< int > ring::findDDC ( std::vector< std::vector< int >>  rings,
std::vector< strucType > *  ringType,
std::vector< int >  listHC,
std::vector< cage::Cage > *  cageList 
)

Find out which hexagonal rings are DDC (Double Diamond Cages) rings. Returns a vector containing all the ring IDs which are DDC rings

Determines which hexagonal rings are DDC rings. This function returns a vector which contains the ring IDs of all the rings which are DDC rings. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. DDC rings can be found using a three-step procedure, in which equatorial rings and their corresponding rings can be found. Peripheral rings can be shared, so first all the equatorial rings must be found. Whenever an equatorial ring and the corresponding peripheral rings are found, their values must be updated to DDC enum type inside the ringType vector, which has been passed as an input to this function. At the end, the output vector can be updated to avoid adding the same ring index more than once (this may happen for peripheral rings, which can be shared). Internally, the function calls the following:

  • ring::conditionOneDDC (Step one: Finds all rings which contain each index (m_k) of the equatorial ring, iring, in at least three other rings).
  • ring::conditionTwoDDC (Step two: For every triplet in iring, there is at least one hexagonal ring other than iring that passes through the triplet. The peripheral rings are stored in order of the starting element of each triplet.)
  • ring::conditionThreeDDC (Step three: For every triplet in the equatorial ring, there is at least one hexagonal ring other than iring that passes through the triplet. Rings corresponding to triplets need not be searched again since peripheralRings are stored in that order. Rings corresponding to T1, T3, T5 must have a common element. Similarly rings corresponding to T2, T4, T6 must have at least one common element. Alternating rings corresponding to triplets must have at least three common elements).
Parameters
[in]ringsVector of vectors containing 6-membered primitive rings (wherein each ring contains the atom indices of the particles which constitute the ring).
[in]ringTypeVector containing a ring::strucType (structure type) for each ring.
[in]cageListVector in which every cage is saved.
Returns
A vector of all the ring indices which constitute DDCs.

Definition at line 199 of file topo_bulk.cpp.

202  {
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 }
T erase(T... args)
@ DoubleDiaC
Definition: cage.hpp:50
bool conditionOneDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
Definition: topo_bulk.cpp:330
bool conditionTwoDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
Definition: topo_bulk.cpp:402
bool conditionThreeDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings)
Definition: topo_bulk.cpp:483
@ 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
@ 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
T insert(T... args)
T resize(T... args)

◆ findHC()

std::vector< int > ring::findHC ( std::vector< std::vector< int >>  rings,
std::vector< strucType > *  ringType,
std::vector< std::vector< int >>  nList,
std::vector< cage::Cage > *  cageList 
)

Find out which hexagonal rings are HC rings. Returns a vector containing all the ring IDs which are HC rings

Determines which hexagonal rings are HCs. This function returns a vector which contains the ring IDs of all the rings which are HC rings. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. HC rings can be found using a three-step procedure, in which first two basal rings are found. Prismatic rings are simply rings which share every face made by upper and lower triplets of the basal rings The neighbour list is also required as an input, which is a vector of vectors, containing atom IDs. The first element of the neighbour list is the atomID of the atom for which the other elements are nearest neighbours. The function calls the following functions internally:

  • ring::hasCommonElements (Step one: Checks to see if basal1 and basal2, two candidate basal rings, have common elements or not. If they don't, then they cannot be basal rings).
  • ring::basalConditions (Step two and three: One of the elements of basal2 must be the nearest neighbour of either the first (index0; \( l_1 \)) or second (index1; \( l_2 \)) element of basal1. If \( m_k \) is the nearest neighbour of \( l_1 \), \( m_{k+2} \) and \( m_{k+4} \) must be neighbours of \( l_3 \) and \( l_5 \)( \( l_5 \) or \( l_3 \)). Similarly modified for \( l_2 \)).
    Parameters
    [in]ringsVector of vectors containing 6-membered primitive rings (wherein each ring contains the atom indices of the particles which constitute the ring).
    [in]ringTypeVector containing a ring::strucType (structure type) for each ring.
    [in]nListRow-ordered vector of vectors of the neighbour list (by index).
    [in]cageListVector in which every cage is saved.
    Returns
    A vector of all the ring indices which constitute HCs.

Definition at line 565 of file topo_bulk.cpp.

568  {
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 }
@ HexC
Definition: cage.hpp:50
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

◆ findMixedRings()

std::vector< int > ring::findMixedRings ( std::vector< std::vector< int >>  rings,
std::vector< strucType > *  ringType,
std::vector< int > *  listDDC,
std::vector< int > *  listHC 
)

Find out which hexagonal rings are both DDCs (Double Diamond Cages) and HCs (Hexagonal Cages). Returns a vector containing all the ring IDs which are of this type

Determines which hexagonal rings are both DDCs and HCs. This function returns a vector which contains the ring IDs of all the rings which are both. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. Rings which are both are of enum type bothBasal OR bothPrismatic. Reassign rings which are mixed in listDDC and listHC as the dummy value -10.

Parameters
[in]ringsThe 6-membered primitive rings.
[in]ringTypeStructure type for every ring.
[in]listDDCContains a vector of ring indices of DDCs.
[in]listHCContains a vector of ring indices of HCs.
Returns
A vector containing the indices of rings which are mixed.

Definition at line 1031 of file topo_bulk.cpp.

1034  {
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 }

◆ findPrismatic()

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

Finding prismatic rings when passed the information in the ringType input vector.

Parameters
[in]ringsThe 6-membered primitive rings.
[in]listHCVector containing the ring indices of rings which are part of HCs.
[in]ringTypeContains a structure type for each ring.
[in]iringIndex of the \( i^{th} \) ring.
[in]jringIndex of the \( j^{th} \) ring.
[in]prismaticRingsVector containing the indices of rings which are prismatic.

Definition at line 935 of file topo_bulk.cpp.

938  {
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 }

◆ getAtomTypesTopoBulk()

int ring::getAtomTypesTopoBulk ( std::vector< std::vector< int >>  rings,
std::vector< ring::strucType ringType,
std::vector< cage::iceType > *  atomTypes 
)

Assigns a type of enum iceType, to every atom, using information from ringType, which has the information of every ring

Assigns a type (cage::iceType) to each atom, according to the previous classification of every ring in ringType. Each subring or vector inside the vector of vector rings, is by index, meaning that the atoms are saved by their indices starting from 0 in the PointCloud.

Parameters
[in]ringsVector of vectors of 6-membered rings.
[in]ringTypeVector containing the structural type for each ring.
[in]atomTypesStructural type for each atom.

Definition at line 1078 of file topo_bulk.cpp.

1080  {
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 }
iceType
Definition: cage.hpp:71
@ 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
@ deformedPrism
Definition: ring.hpp:120
@ mixedPrismRing
Definition: ring.hpp:121

◆ getStrucNumbers()

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

Determines the number of HCs, DDCs from the cageList vector, containing a list of cages. The number of mixed rings, prismatic rings and basal rings are obtained from the ringType vector.

Parameters
[in]ringTypeVector containing the structural type for each ring.
[in]cageListContains all the cages (DDCs and HCs).
[in]numHCThe number of HCs.
[in]numDDCThe number of DDCs.
[in]mixedRingsThe number of mixed rings.
[in]prismaticRingsThe number of prismatic rings (of HCs/mixed rings).
[in]basalRingsTThe number of basal rings (of HCs/mixed rings).

Definition at line 1161 of file topo_bulk.cpp.

1164  {
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

◆ notNeighboursOfRing()

bool ring::notNeighboursOfRing ( std::vector< std::vector< int >>  nList,
std::vector< int > *  triplet,
std::vector< int > *  ring 
)

Tests to check that elements of a triplet are not neighbours of a ring (vector) passed

Checks to make sure that the elements of the triplet are NOT neighbours of any elements inside a vector (ring) passed in (false) If any of them are neighbours, this function returns false.

Parameters
[in]nListVector of vectors of the neighbour list (by index).
[in]tripletVector containing the current triplet being tested.
[in]ringRing passed in.
Returns
A bool; true if the condition is met and false otherwise.

Definition at line 895 of file topo_bulk.cpp.

897  {
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 }

◆ relaxedPrismConditions()

bool prism3::relaxedPrismConditions ( std::vector< std::vector< int >>  nList,
std::vector< int > *  basal1,
std::vector< int > *  basal2 
)

Reduced criterion: Two candidate basal rings of a prism block should have at least one bond between them

Relaxed criteria for deformed prism blocks: at least one bond should exist between the basal rings.

Definition at line 1430 of file topo_bulk.cpp.

1432  {
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 }

◆ topoBulkAnalysis()

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

Find out which rings are DDCs or HCs, which are comprised of 6-membered primitive rings. Start with a neighbour list (by index) and a vector of vectors of rings (also by index). TODO: try 'square' ice and ice0

Finds out if rings constitute double-diamond cages or hexagonal cages. Requires a neighbour list (by index) and a vector of vectors of primitive rings which should also be by index. This is registered as a Lua function and is accessible to the user. Internally, this function calls the following functions:

  • ring::getSingleRingSize (Saves rings of a single ring size into a new vector of vectors, which is subsequently used for finding DDCs, HCs etc).
  • ring::findDDC (Finds the DDCs).
  • ring::findHC (Finds the HCs).
  • ring::findMixedRings (Finds the mixed rings, which are shared by DDCs and HCs both).
  • ring::getStrucNumbers (Gets the number of structures (DDCs, HCs, mixed rings, basal rings, prismatic rings, to be used for write-outs).
  • sout::writeTopoBulkData (Writes out the numbers and data obtained for the current frame).
  • ring::getAtomTypesTopoBulk (Gets the atom type for every atom, to be used for printing out the ice types found).
  • sout::writeLAMMPSdataTopoBulk (Writes out the atoms, with the classified types, into a LAMMPS data file, which can be visualized in OVITO).
    Parameters
    [in]pathThe file path of the output directory to which output files will be written.
    [in]ringsVector of vectors containing the primitive rings. This contains rings of all sizes.
    [in]nListRow-ordered neighbour list, by index.
    [in]yCloudThe input PointCloud, with respect to which the indices in the rings and nList vector of vectors have been saved.
    [in]firstFrameFirst frame to be analyzed
    [in]onlyTetrahedralFlag for only finding DDCs and HCs (true) or also finding PNCs (false)

Definition at line 48 of file topo_bulk.cpp.

52  {
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 }
int nop
Current frame number.
Definition: mol_sys.hpp:169
int currentFrame
Collection of points.
Definition: mol_sys.hpp:168
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
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
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
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
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
Definition: ring.cpp:18
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
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)