Ring
+ Collaboration diagram for Ring:

Modules

 Prism3
 

Enumerations

enum  ring::strucType {
  ring::unclassified , ring::DDC , ring::HCbasal , ring::HCprismatic ,
  ring::bothBasal , ring::bothPrismatic , ring::Prism , ring::deformedPrism ,
  ring::mixedPrismRing
}
 

Functions

int tum3::topoUnitMatchingBulk (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 printClusters, bool onlyTetrahedral)
 
Eigen::MatrixXd tum3::buildRefHC (std::string fileName)
 Build a reference Hexagonal cage, reading in from a template XYZ file. More...
 
Eigen::MatrixXd tum3::buildRefDDC (std::string fileName)
 Build a reference Double-Diamond cage, reading in from a template XYZ file. More...
 
int tum3::shapeMatchHC (molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, cage::Cage cageUnit, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, std::vector< double > *quat, double *rmsd)
 Shape-matching for a target HC. More...
 
int tum3::shapeMatchDDC (molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, std::vector< cage::Cage > cageList, int cageIndex, std::vector< std::vector< int >> rings, std::vector< double > *quat, double *rmsd)
 Shape-matching for a target DDC. More...
 
int tum3::updateRMSDatom (std::vector< std::vector< int >> rings, cage::Cage cageUnit, double rmsd, std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms, std::vector< cage::iceType > atomTypes)
 
int tum3::averageRMSDatom (std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms)
 Average the RMSD per atom. More...
 
std::vector< cage::Cagetum3::topoBulkCriteria (std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, int *numHC, int *numDDC, std::vector< ring::strucType > *ringType)
 
int tum3::clusterCages (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList, int numHC, int numDDC)
 
std::vector< int > tum3::atomsFromCages (std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList, std::vector< int > clusterCages)
 Gets the atoms in the cages of a given cluster. More...
 
std::vector< std::vector< int > > ring::getSingleRingSize (std::vector< std::vector< int >> rings, int ringSize)
 Returns a vector of vectors of rings of a single size. More...
 
bool ring::hasCommonElements (std::vector< int > ring1, std::vector< int > ring2)
 
bool ring::compareRings (std::vector< int > ring1, std::vector< int > ring2)
 
bool ring::findTripletInRing (std::vector< int > ring, std::vector< int > triplet)
 Searches a particular ring for a triplet. More...
 
bool ring::commonElementsInThreeRings (std::vector< int > ring1, std::vector< int > ring2, std::vector< int > ring3)
 Common elements in 3 rings. More...
 
std::vector< int > ring::findsCommonElements (std::vector< int > ring1, std::vector< int > ring2)
 Returns the common elements of two rings. More...
 
int ring::clearRingList (std::vector< std::vector< int >> &rings)
 Erases memory for a vector of vectors for a list of rings. More...
 
std::vector< int > ring::findPrisms (std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, int *nPerfectPrisms, int *nImperfectPrisms, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, bool doShapeMatching=false)
 
bool ring::basalPrismConditions (std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
 
bool ring::relaxedPrismConditions (std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
 
bool ring::discardExtraTetragonBlocks (std::vector< int > *basal1, std::vector< int > *basal2, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 
std::vector< std::vector< int > > ring::keepAxialRingsOnly (std::vector< std::vector< int >> rings, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 Saves only axial rings out of all possible rings. More...
 
int ring::prismAnalysis (std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, int *atomID, int firstFrame, int currentFrame, bool doShapeMatching=false)
 
int ring::assignPrismType (std::vector< std::vector< int >> rings, std::vector< int > listPrism, int ringSize, std::vector< ring::strucType > ringType, std::vector< int > *atomTypes, std::vector< ring::strucType > *atomState)
 
int ring::deformedPrismTypes (std::vector< ring::strucType > atomState, std::vector< int > *atomTypes, int maxDepth)
 Get the atom type values for deformed prisms. More...
 
int ring::rmAxialTranslations (molSys::PointCloud< molSys::Point< double >, double > *yCloud, int *atomID, int firstFrame, int currentFrame)
 Shift the entire ice nanotube and remove axial translations. More...
 
int ring::polygonRingAnalysis (std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, double sheetArea, int firstFrame)
 
int ring::assignPolygonType (std::vector< std::vector< int >> rings, std::vector< int > *atomTypes, std::vector< int > nRings)
 

Detailed Description

Enumeration Type Documentation

◆ strucType

Enumerator
unclassified 

The ring is unclassified, which may be either water or a deformed type which cannot be classified by the criteria.

DDC 

The ring belongs to a double-diamond cage (DDC).

HCbasal 

The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an HC.

HCprismatic 

The ring belongs only to a hexagonal cage (HC); specifically the ring is purely a prismatic ring of an HC. It is not shared by a DDC.

bothBasal 

The ring belongs to both a DDC and HC. It is a 'mixed' ring. The ring is also one of the basal rings of the HC of which it is part. A mixed ring must be a peripheral ring of the DDC of which it is part by definition.

bothPrismatic 

The ring belongs to both a DDC and HC and is, thus, a 'mixed' ring. The ring is also one of the prismatic rings of the HC of which it is part. A mixed ring must be a peripheral ring of the DDC of which it is part by definition (can never be an equatorial ring of a DDC and also be part of an HC).

Prism 

The ring belongs to a prism block, classified according to the prism identification scheme.

deformedPrism 
mixedPrismRing 

Definition at line 112 of file ring.hpp.

112  {
113  unclassified,
114  DDC,
115  HCbasal,
116  HCprismatic,
117  bothBasal,
119  Prism,
122 };
@ 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

Function Documentation

◆ assignPolygonType()

int ring::assignPolygonType ( std::vector< std::vector< int >>  rings,
std::vector< int > *  atomTypes,
std::vector< int >  nRings 
)

Assign an atomType (equal to the number of nodes in the ring) given n-membered rings.

Assign an atomType (equal to the number of nodes in the ring) given the rings vector.

Parameters
[in]ringsThe vector of vectors containing the primitive rings, of a particular ring size.
[in,out]atomTypesA vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.
[in]nRingsNumber of rings.

Definition at line 125 of file topo_two_dim.cpp.

127  {
128  // Every value in listPrism corresponds to an index in rings.
129  // Every ring contains atom indices, corresponding to the indices (not atom
130  // IDs) in rings
131  int iring; // Index of current ring
132  int iatom; // Index of current atom
133  int ringSize; // Ring size of the current ring
134  int prevRingSize; // Ring size previously assigned to a point
135 
136  // Dummy value corresponds to a value of 1.
137  // If an atom is shared by more than one ring type, it is assigned the
138  // value 2.
139 
140  // Loop through every ring in rings
141  for (int iring = 0; iring < rings.size(); iring++) {
142  ringSize = rings[iring].size();
143  // Loop through every element in iring
144  for (int j = 0; j < ringSize; j++) {
145  iatom = rings[iring][j]; // Atom index
146  // Update the atom type
147  if ((*atomTypes)[iatom] == 1) {
148  (*atomTypes)[iatom] = ringSize;
149  } // The atom is unclassified
150  else {
151  // Only update the ring type if the number is higher
152  prevRingSize = (*atomTypes)[iatom]; // Previously assigned ring size
153  if (ringSize > prevRingSize) {
154  (*atomTypes)[iatom] = ringSize;
155  } // end of assigning the new ring size
156  } // only update if the number is higher
157  } // end of loop through every atom in iring
158  } // end of loop through every ring
159 
160  return 0;
161 } // end of function
T size(T... args)

◆ assignPrismType()

int ring::assignPrismType ( std::vector< std::vector< int >>  rings,
std::vector< int >  listPrism,
int  ringSize,
std::vector< ring::strucType ringType,
std::vector< int > *  atomTypes,
std::vector< ring::strucType > *  atomState 
)

Assign an atomType (equal to the number of nodes in the ring) given a vector with a list of indices of rings comprising the prisms

Assign an atomType (equal to the number of nodes in the ring) given a vector with a list of indices of rings comprising the prisms. Note that the ring indices in the listPrism vector correspond to the indices in the input rings vector of vectors.

Parameters
[in]ringsThe vector of vectors containing the primitive rings, of a particular ring size.
[in]listPrismThe list of prism blocks found.
[in]ringSizeThe current ring size or number of nodes in each ring.
[in,out]atomTypesA vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.

Definition at line 692 of file topo_one_dim.cpp.

696  {
697  // Every value in listPrism corresponds to an index in rings.
698  // Every ring contains atom indices, corresponding to the indices (not atom
699  // IDs) in rings
700  int iring; // Index of current ring
701  int iatom; // Index of current atom
702  ring::strucType currentState; // Current state of the ring being filled
703 
704  // Dummy value corresponds to a value of 1.
705  // Each value is initialized to the value of 1.
706 
707  // Loop through every ring in rings
708  for (int i = 0; i < listPrism.size(); i++) {
709  iring = listPrism[i];
710  // Get the state of all atoms in the ring
711  currentState = ringType[iring];
712  // Loop through every element in iring
713  for (int j = 0; j < ringSize; j++) {
714  iatom = rings[iring][j]; // Atom index
715  // Update the atom type
716  (*atomTypes)[iatom] = ringSize;
717  //
718  // Update the state of the atom
719  if ((*atomState)[iatom] == ring::unclassified) {
720  (*atomState)[iatom] = currentState;
721  } // Update the unclassified atom
722  else {
723  if ((*atomState)[iatom] != currentState) {
724  // For mixed, there is a preference
725  if (currentState == ring::mixedPrismRing) {
726  (*atomState)[iatom] = currentState;
727  } // fill
728  else if ((*atomState)[iatom] == ring::deformedPrism &&
729  currentState == ring::Prism) {
730  (*atomState)[iatom] = ring::mixedPrismRing;
731  } else if ((*atomState)[iatom] == ring::Prism &&
732  currentState == ring::deformedPrism) {
733  (*atomState)[iatom] = ring::mixedPrismRing;
734  }
735  } //
736  } // already filled?
737 
738  } // end of loop through every atom in iring
739  } // end of loop through every ring
740 
741  return 0;
742 } // end of function
strucType
Definition: ring.hpp:112

◆ atomsFromCages()

std::vector< int > tum3::atomsFromCages ( std::vector< std::vector< int >>  rings,
std::vector< cage::Cage cageList,
std::vector< int >  clusterCages 
)

Gets the atoms in the cages of a given cluster.

Returns a vector containing the atom indices of all the atoms in a cluster of cages, according to the input cageList vector of Cages

Definition at line 797 of file bulkTUM.cpp.

799  {
800  //
801  std::vector<int> atoms; // Contains the atom indices (not IDs) of atoms
802  int ringSize = rings[0].size(); // Number of nodes in each ring
803  int iring; // Index of ring in rings vector of vectors
804  int icage; // Index of the current cage in cageList
805  int iatom; // Atom index
806  //
807  for (int i = 0; i < clusterCages.size(); i++) {
808  //
809  icage = clusterCages[i];
810  // Loop through every ring in the cage
811  for (int j = 0; j < cageList[icage].rings.size(); j++) {
812  iring = cageList[icage].rings[j]; // Current ring index
813  // Loop through every atom in iring
814  for (int k = 0; k < ringSize; k++) {
815  iatom = rings[iring][k];
816  atoms.push_back(iatom); // Add the atom
817  } // Loop through every atom in iring
818  } // loop through every ring in the current cage
819  } // end of loop through all cages
820 
821  // Remove the duplicates in the atoms vector
822  // Duplicate IDs must be removed
823  //
824  sort(atoms.begin(), atoms.end());
825  auto ip = std::unique(atoms.begin(), atoms.end());
826  // Resize peripheral rings to remove undefined terms
827  atoms.resize(std::distance(atoms.begin(), ip));
828  //
829  // Return the vector of atoms
830  return atoms;
831 }
T begin(T... args)
T distance(T... args)
T end(T... args)
int clusterCages(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList, int numHC, int numDDC)
Definition: bulkTUM.cpp:640
T push_back(T... args)
T resize(T... args)
T sort(T... args)
T unique(T... args)

◆ averageRMSDatom()

int tum3::averageRMSDatom ( std::vector< double > *  rmsdPerAtom,
std::vector< int > *  noOfCommonAtoms 
)

Average the RMSD per atom.

Average the RMSD per atom, by the number of common elements

Definition at line 525 of file bulkTUM.cpp.

526  {
527  //
528  int nop = (*rmsdPerAtom).size(); // Number of particles
529 
530  for (int iatom = 0; iatom < nop; iatom++) {
531  //
532  if ((*noOfCommonAtoms)[iatom] == 0) {
533  (*rmsdPerAtom)[iatom] = -1; // Dummy atom
534  } else {
535  (*rmsdPerAtom)[iatom] /= (*noOfCommonAtoms)[iatom];
536  }
537  } // end of averaging
538 
539  return 0;
540 } // end of function

◆ basalPrismConditions()

bool ring::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 347 of file topo_one_dim.cpp.

349  {
350  int l1 = (*basal1)[0]; // first element of basal1 ring
351  int ringSize =
352  (*basal1).size(); // Size of the ring; each ring contains n elements
353  int m_k; // Atom ID of element in basal2
354  bool l1_neighbour; // m_k is a neighbour of l1(true) or not (false)
355 
356  // isNeighbour is initialized to false for all basal2 elements; indication if
357  // basal2 elements are neighbours of basal1
358  std::vector<bool> isNeighbour(ringSize, false);
359  int kIndex; // m_k index
360  int lAtomID; // atomID of the current element of basal1
361  int kAtomID; // atomID of the current element of basal2
362 
363  // ---------------------------------------------
364  // COMPARISON OF basal2 ELEMENTS WITH l1
365  for (int k = 0; k < ringSize; k++) {
366  l1_neighbour = false;
367  m_k = (*basal2)[k];
368  // =================================
369  // Checking to seee if m_k is be a neighbour of l1
370  // Find m_k inside l1 neighbour list
371  auto it = std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
372 
373  // If the element has been found, for l1
374  if (it != nList[l1].end()) {
375  l1_neighbour = true;
376  kIndex = k;
377  break;
378  }
379  } // l1 is a neighbour of m_k
380 
381  // If there is no nearest neighbour, then the two rings are not part of the
382  // prism
383  if (!l1_neighbour) {
384  return false;
385  }
386 
387  // ---------------------------------------------
388  // NEIGHBOURS of basal1 in basal2
389  isNeighbour[kIndex] = true;
390 
391  // All elements of basal1 must be neighbours of basal2
392  for (int i = 1; i < ringSize; i++) {
393  lAtomID = (*basal1)[i]; // element of basal1 ring
394  for (int k = 0; k < ringSize; k++) {
395  // Skip if already a neighbour
396  if (isNeighbour[k]) {
397  continue;
398  }
399  // Get the comparison basal2 element
400  kAtomID = (*basal2)[k]; // element of basal2 ring;
401 
402  // Checking to see if kAtomID is a neighbour of lAtomID
403  // Find kAtomID inside lAtomID neighbour list
404  auto it1 =
405  std::find(nList[lAtomID].begin() + 1, nList[lAtomID].end(), kAtomID);
406 
407  // If the element has been found, for l1
408  if (it1 != nList[lAtomID].end()) {
409  isNeighbour[k] = true;
410  }
411  } // Loop through basal2
412  } // end of check for neighbours of basal1
413 
414  // ---------------------------------------------
415 
416  // They should all be neighbours
417  for (int k = 0; k < ringSize; k++) {
418  // Check to see if any element is false
419  if (!isNeighbour[k]) {
420  return false;
421  }
422  }
423 
424  // Everything works out!
425  return true;
426 }
T find(T... args)

◆ buildRefDDC()

Eigen::MatrixXd tum3::buildRefDDC ( std::string  fileName)

Build a reference Double-Diamond cage, reading in from a template XYZ file.

Build a reference Double-Diamond cage, reading it in from a templates directory

Definition at line 431 of file bulkTUM.cpp.

431  {
432  //
433  Eigen::MatrixXd refPnts(14, 3); // Reference point set (Eigen matrix)
434  // Get the reference HC point set
436  setCloud; // PointCloud for holding the reference point values
437  // Variables for rings
438  std::vector<std::vector<int>> nList; // Neighbour list
439  std::vector<std::vector<int>> rings; // Rings
441  ringType; // This vector will have a value for each ring inside
442  std::vector<int> listDDC,
443  listHC; // Contains atom indices of atoms making up DDCs and HCs
444  std::vector<int> ddcOrder; // Atom indices of particles in the DDC
445  // Make a list of all the DDCs and HCs
446  std::vector<cage::Cage> cageList;
447  int iring, jring;
448  //
449  // read in the XYZ file into the pointCloud setCloud
450  //
451  sinp::readXYZ(fileName, &setCloud);
452 
453  // box lengths
454  for (int i = 0; i < 3; i++) {
455  setCloud.box[i] = 50;
456  } // end of setting box lengths
457  //
458 
459  nList = nneigh::neighListO(3.5, &setCloud, 1);
460  // Neighbour list by index
461  nList = nneigh::neighbourListByIndex(&setCloud, nList);
462  // Find the vector of vector of rings
463  rings = primitive::ringNetwork(nList, 6);
464  // init the ringType vector
465  ringType.resize(rings.size());
466  // Find the DDCs
467  listDDC = ring::findDDC(rings, &ringType, listHC, &cageList);
468  // Save the order of the DDC in a vector
469  ddcOrder = pntToPnt::relOrderDDC(0, rings, cageList);
470  // Get the reference point set
471  refPnts = pntToPnt::changeDiaCageOrder(&setCloud, ddcOrder, 0);
472  //
473  return refPnts;
474 }
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:170
std::vector< std::vector< int > > neighListO(double rcutoff, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI)
Definition: neighbours.cpp:114
std::vector< std::vector< int > > neighbourListByIndex(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList)
Definition: neighbours.cpp:337
std::vector< std::vector< int > > ringNetwork(std::vector< std::vector< int >> nList, int maxDepth)
Definition: franzblau.cpp:28
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 readXYZ(std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Function for reading in atom coordinates from an XYZ file.
Definition: seams_input.cpp:69
Eigen::MatrixXd changeDiaCageOrder(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > ddcOrder, int startingIndex=0)
std::vector< int > relOrderDDC(int index, std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList)
Matches the order of the basal rings of an DDC or a potential HC.
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166

◆ buildRefHC()

Eigen::MatrixXd tum3::buildRefHC ( std::string  fileName)

Build a reference Hexagonal cage, reading in from a template XYZ file.

Build a reference Hexagonal cage, reading it in from a templates directory

Definition at line 373 of file bulkTUM.cpp.

373  {
374  //
375  Eigen::MatrixXd refPnts(12, 3); // Reference point set (Eigen matrix)
376  // Get the reference HC point set
378  setCloud; // PointCloud for holding the reference point values
379  // Variables for rings
380  std::vector<std::vector<int>> nList; // Neighbour list
381  std::vector<std::vector<int>> rings; // Rings
383  ringType; // This vector will have a value for each ring inside
384  std::vector<int> listHC; // Contains atom indices of atoms making up HCs
385  // Make a list of all the DDCs and HCs
386  std::vector<cage::Cage> cageList;
387  int iring, jring;
388  //
389  // read in the XYZ file into the pointCloud setCloud
390  //
391  sinp::readXYZ(fileName, &setCloud);
392  // box lengths
393  for (int i = 0; i < 3; i++) {
394  setCloud.box[i] = 50;
395  } // end of setting box lengths
396  //
397 
398  nList = nneigh::neighListO(3.5, &setCloud, 1);
399  // Neighbour list by index
400  nList = nneigh::neighbourListByIndex(&setCloud, nList);
401  // Find the vector of vector of rings
402  rings = primitive::ringNetwork(nList, 6);
403  // init the ringType vector
404  ringType.resize(rings.size());
405  // Find the HCs
406  listHC = ring::findHC(rings, &ringType, nList, &cageList);
407  // Get the basal rings from cageList
408  iring = cageList[0].rings[0];
409  jring = cageList[0].rings[1];
410  //
411  std::vector<int> matchedBasal1,
412  matchedBasal2; // Re-ordered basal rings 1 and 2
413  // Reordered basal rings
414  // Getting the target Eigen vectors
415  // Get the re-ordered matched basal rings, ordered with respect to each
416  // other
417 
418  pntToPnt::relOrderHC(&setCloud, rings[iring], rings[jring], nList,
419  &matchedBasal1, &matchedBasal2);
420  // Get the reference point set
421  refPnts =
422  pntToPnt::changeHexCageOrder(&setCloud, matchedBasal1, matchedBasal2, 0);
423  //
424  return refPnts;
425 }
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 relOrderHC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, std::vector< std::vector< int >> nList, std::vector< int > *matchedBasal1, std::vector< int > *matchedBasal2)
Matches the order of the basal rings of an HC or a potential HC.
Eigen::MatrixXd changeHexCageOrder(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex=0)

◆ clearRingList()

int ring::clearRingList ( std::vector< std::vector< int >> &  rings)

Erases memory for a vector of vectors for a list of rings.

Deletes the memory of a vector of vectors (specifically for rings).

Parameters
[in,out]ringsThe vector of vectors to be cleared.

Definition at line 18 of file ring.cpp.

18  {
19  //
21 
22  rings.swap(tempEmpty);
23 
24  return 0;
25 }
T swap(T... args)

◆ clusterCages()

int tum3::clusterCages ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::string  path,
std::vector< std::vector< int >>  rings,
std::vector< cage::Cage cageList,
int  numHC,
int  numDDC 
)

Clustering Clusters cages using the Stillinger algorithm and prints out individual XYZ files of clusters.

Groups cages into clusters and prints out each cluster into an XYZ file

Definition at line 640 of file bulkTUM.cpp.

643  {
644  //
645  std::vector<int> linkedList; // Linked list for the clustered cages
646  int nCages = numHC + numDDC; // Number of cages in total
647  int j, temp; // Variables used inside the loops
648  bool isNeighbour; // true if cages are adjacent to each other, false otherwise
650  visited; // To make sure you don't go through the same cages again.
651  int singleDDCs, singleHCs; // Number of single DDCs and HCs
652  // Units are in Angstrom^3
653  double volDDC = 72.7; // The calculated alpha volume of a single DDC
654  double volHC = 40.63; // The calculated alpha volume of a single HC
655  int nextElement; // value in list
656  int index; // starting index value
657  int currentIndex; // Current cage index
658  int iClusterNumber; // Number in the current cluster
659  std::vector<int> nClusters; // Number of cages in each cluster
660  std::vector<int> clusterCages; // Indices of each cage in the cluster
661  std::vector<int> atoms; // Vector containing the atom indices of all the
662  // atoms in the cages
663  cage::cageType type; // Type of the cages in the cluster
664  // -----------------------------------------------------------
665  // INITIALIZATION
666  linkedList.resize(nCages); // init to dummy value
667  // Initial values of the list.
668  for (int icage = 0; icage < nCages; icage++) {
669  // Assign the index as the ID
670  linkedList[icage] = icage; // Index
671  } // init of cluster IDs
672  // -----------------------------------------------------------
673  // GETTING THE LINKED LIST
674  //
675  // HCs first
676  for (int i = 0; i < nCages - 1; i++) {
677  // If iatom is already in a cluster, skip it
678  if (linkedList[i] != i) {
679  continue;
680  } // Already part of a cluster
681  //
682  j = i; // Init of j
683  // Execute the next part of the loop while j is not equal to i
684  do {
685  //
686  // Go through the rest of the atoms (KLOOP)
687  for (int k = i + 1; k < nCages; k++) {
688  // Skip if already part of a cluster
689  if (linkedList[k] != k) {
690  continue;
691  } // Already part of a cluster
692  //
693  // k and j must be of the same type
694  if (cageList[k].type != cageList[j].type) {
695  continue;
696  } // skip if k and j are not of the same type
697  // Check to see if k is a nearest neighbour of j
698  isNeighbour =
699  ring::hasCommonElements(cageList[k].rings, cageList[j].rings);
700  if (isNeighbour) {
701  // Swap!
702  temp = linkedList[j];
703  linkedList[j] = linkedList[k];
704  linkedList[k] = temp;
705  } // j and k are neighbouring cages
706  } // end of loop through k (KLOOP)
707  //
708  j = linkedList[j];
709  } while (j != i); // end of control for j!=i
710  } // end of getting the linked list (looping through every i)
711  // -----------------------------------------------------------
712  // WRITE-OUTS
713  //
714  // init
715  visited.resize(nCages);
716  singleDDCs = 0;
717  singleHCs = 0;
718 
719  for (int i = 0; i < nCages; i++) {
720  //
721  if (visited[i]) {
722  continue;
723  } // Already counted
724  // Now that the cage has been visited, set it to true
725  visited[i] = true; // Visited
726  // If only one cage is in the cluster
727  if (linkedList[i] == i) {
728  // Add to the number of single DDCs
729  if (cageList[i].type == cage::DoubleDiaC) {
730  singleDDCs++;
731  } // add to the single DDCs
732  else {
733  singleHCs++;
734  } // single HCs
735  continue;
736  } // cluster has only one cage
737  //
738  // init
739  clusterCages.resize(0);
740  //
741  type = cageList[i].type;
742  currentIndex = i;
743  nextElement = linkedList[currentIndex];
744  index = i;
745  iClusterNumber = 1; // at least one value
746  clusterCages.push_back(i); // Add the first cage
747  // Get the other cages in the cluster
748  while (nextElement != index) {
749  iClusterNumber++;
750  currentIndex = nextElement;
751  visited[currentIndex] = true;
752  clusterCages.push_back(currentIndex); // Add the first cage
753  nextElement = linkedList[currentIndex];
754  } // get number
755  // Update the number of molecules in the cluster
756  nClusters.push_back(iClusterNumber);
757  // --------------
758  // PRINT OUT
759  atoms = tum3::atomsFromCages(rings, cageList, clusterCages);
760  int clusterID = nClusters.size();
761  // Write out the cluster
762  sout::writeXYZcluster(path, yCloud, atoms, clusterID, type);
763  // Print out the cages into an XYZ file
764  // --------------
765  } // end of loop through
766  // -----------------------------------------------------------
767  // Write out the stuff for single cages
768  // ----------------
769  std::ofstream outputFile;
770  // Make the output directory if it doesn't exist
771  sout::makePath(path);
772  std::string outputDirName = path + "bulkTopo";
773  sout::makePath(outputDirName);
774  outputDirName = path + "bulkTopo/clusterXYZ/";
775  sout::makePath(outputDirName);
776  outputDirName = path + "bulkTopo/clusterXYZ/frame-" +
777  std::to_string(yCloud->currentFrame);
778  sout::makePath(outputDirName);
779  // ----------------
780  // Write output to file inside the output directory
781  outputFile.open(path + "bulkTopo/clusterXYZ/frame-" +
782  std::to_string(yCloud->currentFrame) + "/info.dat");
783  outputFile << "# volDDC volHC in Angstrom^3\n";
784  outputFile << singleDDCs * volDDC << " " << singleHCs * volHC << "\n";
785  outputFile << "There are " << singleDDCs << " single DDCs\n";
786  outputFile << "There are " << singleHCs << " single HCs\n";
787  outputFile.close(); // Close the file
788  // -----------------------------------------------------------
789  return 0;
790 } // end of the function
T close(T... args)
cageType
Definition: cage.hpp:50
@ DoubleDiaC
Definition: cage.hpp:50
int currentFrame
Collection of points.
Definition: mol_sys.hpp:168
std::vector< int > atomsFromCages(std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList, std::vector< int > clusterCages)
Gets the atoms in the cages of a given cluster.
Definition: bulkTUM.cpp:797
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Definition: ring.cpp:175
int makePath(const std::string &path)
int writeXYZcluster(std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > atoms, int clusterID, cage::cageType type)
Function for writing out the XYZ files for each cluster.
T open(T... args)
T to_string(T... args)

◆ commonElementsInThreeRings()

bool ring::commonElementsInThreeRings ( std::vector< int >  ring1,
std::vector< int >  ring2,
std::vector< int >  ring3 
)

Common elements in 3 rings.

Function that finds the common elements in three input rings

Parameters
[in]ring1The first ring.
[in]ring2The second ring.
[in]ring3The third ring.
Returns
A value which is true if the three rings have at least one common element, and false if the three rings have no elements in common.

Definition at line 63 of file ring.cpp.

65  {
67  common1; // Vector containing the common elements of the first two rings
69  common2; // Vector containing the common elements of the three rings
70 
71  // Common elements among the first two rings
72  common1 = ring::findsCommonElements(ring1, ring2);
73  if (common1.size() == 0) {
74  return false;
75  } // no common elements in ring1 and ring2
76 
77  // Common elements among all three rings
78  common2 = ring::findsCommonElements(common1, ring3);
79 
80  // If no common elements were found:
81  if (common2.size() == 0) {
82  return false;
83  } // no common elements between ring1, ring2, and ring3
84 
85  return true; // Common elements found!
86 }
std::vector< int > findsCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Returns the common elements of two rings.
Definition: ring.cpp:34

◆ compareRings()

bool ring::compareRings ( std::vector< int >  ring1,
std::vector< int >  ring2 
)

Compares two disordered vectors and checks to see if they contain the same elements

Checks to see if two vectors (ring1, ring2) have the same elements (which may or may not be disordered). The order or sequence of elements is not important, so the rings are sorted. Returns true if the rings have the same elements

Parameters
[in]ring1The first ring.
[in]ring2The second ring.
Returns
A bool; true if the rings contain the same elements (not necessarily in the same sequence) and false if they do not have the same elements.

Definition at line 207 of file ring.cpp.

207  {
208  // Sort the rings first
209  sort(ring1.begin(), ring1.end()); // Sort ring1 by ID
210  sort(ring2.begin(), ring2.end()); // Sort ring2 by ID
211  bool result;
212 
213  (ring1 == ring2) ? result = true : result = false;
214 
215  return result;
216 }

◆ deformedPrismTypes()

int ring::deformedPrismTypes ( std::vector< ring::strucType atomState,
std::vector< int > *  atomTypes,
int  maxDepth 
)

Get the atom type values for deformed prisms.

Assign an atomType value for atoms belonging to deformed prisms.

Parameters
[in]ringsThe vector of vectors containing the primitive rings, of a particular ring size.
[in]listPrismThe list of prism blocks found.
[in]ringSizeThe current ring size or number of nodes in each ring.
[in,out]atomTypesA vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.

Definition at line 755 of file topo_one_dim.cpp.

756  {
757  //
758  int nop = atomState.size(); // Number of particles
759 
760  // Loop through all particles
761  for (int iatom = 0; iatom < nop; iatom++) {
762  // Check the atom state
763  // Deformed
764  if (atomState[iatom] == ring::deformedPrism) {
765  (*atomTypes)[iatom] += maxDepth - 2;
766  } // type for a deformed prism atom
767  else if (atomState[iatom] == ring::mixedPrismRing) {
768  (*atomTypes)[iatom] = 2;
769  } // type for a mixed prism ring
770  } // end of reassignation
771 
772  //
773  return 0;
774 }

◆ discardExtraTetragonBlocks()

bool ring::discardExtraTetragonBlocks ( std::vector< int > *  basal1,
std::vector< int > *  basal2,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud 
)

Checks whether two 4-membered rings are parallel in one dimension or not to prevent overcounting

A function that discards basal tetragonal rings which are not oriented perpendicular to the axial direction. This is will only work for the XY, YZ or XZ planes. If true is returned, then the rings are axial. This function is only needed for tetragonal prism blocks because in the other cases, the basal rings and lateral planes have a different number of nodes.

Parameters
[in]basal1The first basal ring.
[in]basal2The other candidate basal ring.
[in]yCloudThe input PointCloud.
Returns
A bool value which is true if the candidate basal rings are axial, and is false if the planes are lateral planes.

Definition at line 483 of file topo_one_dim.cpp.

485  {
486  int ringSize =
487  (*basal1).size(); // Size of the ring; each ring contains n elements
488  int iatomIndex,
489  jatomIndex; // Indices of the elements in basal1 and basal2 respectively
490  double r_i, r_j; // Coordinates in the axial dimension of iatom and jatom of
491  // basal1 and basal2 respectively
492  int axialDim; // 0 for x, 1 for y and 2 for z dimensions respectively
493  // Variables for getting the projected area
494  bool axialBasal1, axialBasal2; // bools for checking if basal1 and basal2 are
495  // axial (true) respectively
496  double areaXY, areaXZ,
497  areaYZ; // Projected area on the XY, XZ and YZ planes respectively
498  // ----------------------------------------
499  // Find the axial dimension for a quasi-one-dimensional ice nanotube
500  // The axial dimension will have the largest box length
501  // Index -> axial dimension
502  // 0 -> x dim
503  // 1 -> y dim
504  // 2 -> z dim
505  axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
506  yCloud->box.begin();
507  // ----------------------------------------
508  // Calculate projected area onto the XY, YZ and XZ planes for basal1
509  axialBasal1 = false; // Init to false
510  axialBasal2 = false; // Init
511 
512  // Init the projected area
513  areaXY = 0.0;
514  areaXZ = 0.0;
515  areaYZ = 0.0;
516 
517  jatomIndex = (*basal1)[0];
518 
519  // All points except the first pair
520  for (int k = 1; k < ringSize; k++) {
521  iatomIndex = (*basal1)[k]; // Current vertex
522 
523  // Add to the polygon area
524  // ------
525  // XY plane
526  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
527  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
528  // ------
529  // XZ plane
530  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
531  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
532  // ------
533  // YZ plane
534  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
535  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
536  // ------
537  jatomIndex = iatomIndex;
538  }
539 
540  // Closure point
541  iatomIndex = (*basal1)[0];
542  // ------
543  // XY plane
544  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
545  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
546  // ------
547  // XZ plane
548  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
549  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
550  // ------
551  // YZ plane
552  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
553  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
554  // ------
555  // The actual projected area is half of this
556  areaXY *= 0.5;
557  areaXZ *= 0.5;
558  areaYZ *= 0.5;
559  // Get the absolute value
560  areaXY = fabs(areaXY);
561  areaXZ = fabs(areaXZ);
562  areaYZ = fabs(areaYZ);
563 
564  // If the axial dimension is x, y, or z:
565  // then the maximum basal area should be in the YZ, XZ and XY dimensions
566  // respectively
567  // x dim
568  if (axialDim == 0) {
569  if (areaYZ > areaXY && areaYZ > areaXZ) {
570  axialBasal1 = true;
571  } // end of check for axial ring for basal1
572  } // x dim
573  // y dim
574  else if (axialDim == 1) {
575  if (areaXZ > areaXY && areaXZ > areaYZ) {
576  axialBasal1 = true;
577  } // end of check for axial ring for basal1
578  } // x dim
579  // z dim
580  else if (axialDim == 2) {
581  if (areaXY > areaXZ && areaXY > areaYZ) {
582  axialBasal1 = true;
583  } // end of check for axial ring for basal1
584  } // x dim
585  else {
586  std::cerr << "Could not find the axial dimension.\n";
587  return false;
588  }
589  // ----------------------------------------
590  // Calculate projected area onto the XY, YZ and XZ planes for basal2
591 
592  // Init the projected area
593  areaXY = 0.0;
594  areaXZ = 0.0;
595  areaYZ = 0.0;
596 
597  jatomIndex = (*basal2)[0];
598 
599  // All points except the first pair
600  for (int k = 1; k < ringSize; k++) {
601  iatomIndex = (*basal2)[k]; // Current vertex
602 
603  // Add to the polygon area
604  // ------
605  // XY plane
606  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
607  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
608  // ------
609  // XZ plane
610  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
611  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
612  // ------
613  // YZ plane
614  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
615  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
616  // ------
617  jatomIndex = iatomIndex;
618  }
619 
620  // Closure point
621  iatomIndex = (*basal2)[0];
622  // ------
623  // XY plane
624  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
625  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
626  // ------
627  // XZ plane
628  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
629  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
630  // ------
631  // YZ plane
632  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
633  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
634  // ------
635  // The actual projected area is half of this
636  areaXY *= 0.5;
637  areaXZ *= 0.5;
638  areaYZ *= 0.5;
639  // Get the absolute value
640  areaXY = fabs(areaXY);
641  areaXZ = fabs(areaXZ);
642  areaYZ = fabs(areaYZ);
643 
644  // Check if xy projected area is the greatest
645  // If the axial dimension is x, y, or z:
646  // then the maximum basal area should be in the YZ, XZ and XY dimensions
647  // respectively
648  // x dim
649  if (axialDim == 0) {
650  if (areaYZ > areaXY && areaYZ > areaXZ) {
651  axialBasal2 = true;
652  } // end of check for axial ring for basal1
653  } // x dim
654  // y dim
655  else if (axialDim == 1) {
656  if (areaXZ > areaXY && areaXZ > areaYZ) {
657  axialBasal2 = true;
658  } // end of check for axial ring for basal1
659  } // x dim
660  // z dim
661  else if (axialDim == 2) {
662  if (areaXY > areaXZ && areaXY > areaYZ) {
663  axialBasal2 = true;
664  } // end of check for axial ring for basal1
665  } // x dim
666  else {
667  std::cerr << "Could not find the axial dimension.\n";
668  return false;
669  }
670  // ----------------------------------------
671 
672  // Now check if basal1 and basal2 are axial or not
673  if (axialBasal1 == true && axialBasal2 == true) {
674  return true;
675  } else {
676  return false;
677  } // Check for basal1 and basal2
678 }
T fabs(T... args)
std::vector< S > pts
Definition: mol_sys.hpp:167
T max_element(T... args)

◆ findPrisms()

std::vector< int > ring::findPrisms ( std::vector< std::vector< int >>  rings,
std::vector< strucType > *  ringType,
int *  nPerfectPrisms,
int *  nImperfectPrisms,
std::vector< std::vector< int >>  nList,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< double > *  rmsdPerAtom,
bool  doShapeMatching = false 
)

Find out which rings are prisms. Returns a vector containing all the ring IDs which 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 189 of file topo_one_dim.cpp.

193  {
194  std::vector<int> listPrism;
195  int totalRingNum = rings.size(); // Total number of rings
196  std::vector<int> basal1; // First basal ring
197  std::vector<int> basal2; // Second basal ring
198  bool cond1, cond2; // Conditions for rings to be basal (true) or not (false)
199  bool relaxedCond; // Condition so that at least one bond exists between the
200  // two basal rings
201  bool isAxialPair; // Basal rings should be parallel in one dimension to
202  // prevent overcounting
203  int ringSize = rings[0].size(); // Number of nodes in each ring
204  *nImperfectPrisms = 0; // Number of undeformed prisms
205  *nPerfectPrisms = 0; // Number of undeformed prisms
206  // Matrix for the reference ring for a given ringSize.
207  Eigen::MatrixXd refPointSet(ringSize, 3);
208 
209  // Get the reference ring point set for a given ring size.
210  // Get the axial dimension
211  int axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
212  yCloud->box.begin();
213  refPointSet = pntToPnt::getPointSetRefRing(ringSize, axialDim);
214  //
215 
216  // Two loops through all the rings are required to find pairs of basal rings
217  for (int iring = 0; iring < totalRingNum - 1; iring++) {
218  cond1 = false;
219  cond2 = false;
220  basal1 = rings[iring]; // Assign iring to basal1
221  // Loop through the other rings to get a pair
222  for (int jring = iring + 1; jring < totalRingNum; jring++) {
223  basal2 = rings[jring]; // Assign jring to basal2
224  // ------------
225  // Put extra check for axial basal rings if shapeMatching is being done
226  if (doShapeMatching == true || ringSize == 4) {
227  isAxialPair = false; // init
228  isAxialPair =
229  ring::discardExtraTetragonBlocks(&basal1, &basal2, yCloud);
230  if (isAxialPair == false) {
231  continue;
232  }
233  } // end of check for tetragonal prism blocks
234  // ------------
235  // Step one: Check to see if basal1 and basal2 have common
236  // elements or not. If they don't, then they cannot be basal rings
237  cond1 = ring::hasCommonElements(basal1, basal2);
238  if (cond1 == true) {
239  continue;
240  }
241  // -----------
242  // Step two and three: One of the elements of basal2 must be the nearest
243  // neighbour of the first (index0; l1) If m_k is the nearest neighbour of
244  // l1, m_{k+1} ... m_{k+(n-1)} must be neighbours of l_i+1 etc or l_i-1
245  cond2 = ring::basalPrismConditions(nList, &basal1, &basal2);
246  // If cond2 is false, the strict criteria for prisms has not been met
247  if (cond2 == false) {
248  // Skip if shape-matching is not desired
249  if (!doShapeMatching) {
250  continue;
251  } // shape-matching not desired
252  //
253  // If shape-matching is to be done:
254  // Check for the reduced criteria fulfilment
255  relaxedCond = ring::relaxedPrismConditions(nList, &basal1, &basal2);
256  // Skip if relaxed criteria are not met
257  if (relaxedCond == false) {
258  continue;
259  } // end of skipping if the prisms do not fulfil relaxed criteria
260 
261  // Do shape matching here
262  bool isDeformedPrism = match::matchPrism(
263  yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, false);
264 
265  // Success! The rings are basal rings of a deformed prism!
266  if (isDeformedPrism) {
267  //
268  // // Update the number of prism blocks
269  *nImperfectPrisms += 1;
270  // Update iring
271  if ((*ringType)[iring] == ring::unclassified) {
272  (*ringType)[iring] = ring::deformedPrism;
273  listPrism.push_back(iring);
274  } else if ((*ringType)[iring] == ring::Prism) {
275  (*ringType)[iring] = ring::mixedPrismRing;
276  } // if it is deformed
277  // Update jring
278  if ((*ringType)[jring] == ring::unclassified) {
279  (*ringType)[jring] = ring::deformedPrism;
280  listPrism.push_back(jring);
281  } else if ((*ringType)[jring] == ring::Prism) {
282  (*ringType)[jring] = ring::mixedPrismRing;
283  } // if it is deformed
284  } // end of update of ring types
285 
286  // // Write outs
287  // // Now write out axial basal rings
288  // sout::writeBasalRingsPrism(&basal1, &basal2, nDeformedPrisms, nList,
289  // yCloud, true);
290  } // end of reduced criteria
291  // Strict criteria
292  else {
293  // Update the number of prism blocks
294  *nPerfectPrisms += 1;
295  // Update iring
296  if ((*ringType)[iring] == ring::unclassified) {
297  (*ringType)[iring] = ring::Prism;
298  listPrism.push_back(iring);
299  } else if ((*ringType)[iring] == ring::deformedPrism) {
300  (*ringType)[iring] = ring::mixedPrismRing;
301  } // if it is deformed
302  // Update jring
303  if ((*ringType)[jring] == ring::unclassified) {
304  (*ringType)[jring] = ring::Prism;
305  listPrism.push_back(jring);
306  } else if ((*ringType)[jring] == ring::deformedPrism) {
307  (*ringType)[jring] = ring::mixedPrismRing;
308  } // if it is deformed
309  //
310  // Shape-matching to get the RMSD (if shape-matching is desired)
311  if (doShapeMatching) {
312  bool isKnownPrism = match::matchPrism(
313  yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, true);
314  } // end of shape-matching to get rmsd
315  //
316  // // Now write out axial basal rings for convex hull calculations
317  // sout::writePrisms(&basal1, &basal2, *nPrisms, yCloud);
318  // // Write out prisms for shape-matching
319  // sout::writeBasalRingsPrism(&basal1, &basal2, *nPrisms, nList, yCloud,
320  // false);
321  // -----------
322  } // end of strict criteria
323 
324  } // end of loop through rest of the rings to get the second basal ring
325  } // end of loop through all rings for first basal ring
326 
327  sort(listPrism.begin(), listPrism.end());
328  auto ip = std::unique(listPrism.begin(), listPrism.end());
329  // Resize peripheral rings to remove undefined terms
330  listPrism.resize(std::distance(listPrism.begin(), ip));
331 
332  return listPrism;
333 }
bool relaxedPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
bool discardExtraTetragonBlocks(std::vector< int > *basal1, std::vector< int > *basal2, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
bool basalPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
bool matchPrism(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, bool isPerfect=true)
Definition: shapeMatch.cpp:17
Eigen::MatrixXd getPointSetRefRing(int n, int axialDim)

◆ findsCommonElements()

std::vector< int > ring::findsCommonElements ( std::vector< int >  ring1,
std::vector< int >  ring2 
)

Returns the common elements of two rings.

Function that finds and returns a vector containing the elements shared by two input rings (each ring is a vector).

Parameters
[in]ring1The first ring.
[in]ring2The second ring.
Returns
A vector containing the common elements between the input rings.

Definition at line 34 of file ring.cpp.

35  {
36  //
37  std::vector<int> common;
38  int iatom; // Index to search for
39 
40  for (int i = 0; i < ring1.size(); i++) {
41  iatom = ring1[i];
42  // Search for iatom in ring2
43 
44  auto it = std::find(ring2.begin(), ring2.end(), iatom);
45 
46  if (it != ring2.end()) {
47  common.push_back(iatom);
48  } // iatom was found!
49  } // end of loop through every element of ring1
50 
51  return common;
52 }

◆ findTripletInRing()

bool ring::findTripletInRing ( std::vector< int >  ring,
std::vector< int >  triplet 
)

Searches a particular ring for a triplet.

Function that finds out if a given triplet is is present (in the same order or in reversed order) in a given ring.

Parameters
[in]ringThe input ring containing the indices of atoms.
[in]tripletVector containing the triplet, for whose presence the input ring vector will be checked.
Returns
A bool value which is true if the triplet is present in the ring, and is false if the triplet is not in the ring.

Definition at line 97 of file ring.cpp.

97  {
98  //
99  int ringSize = ring.size(); // should be 6
100  std::vector<int> ringTriplet; // triplet from the ring to be searched
101  int kIndex; // Used for making the triplet
102 
103  // Loop through every possible triplet in the ring to be searched
104  for (int i = 0; i < ringSize; i++) {
105  ringTriplet.clear(); // init
106  // Get the first element of the ring triplet
107  ringTriplet.push_back(ring[i]);
108  //
109  // Get the next two elements
110  for (int k = 1; k < 3; k++) {
111  kIndex = i + k;
112  // Wrap-around
113  if (kIndex >= ringSize) {
114  kIndex -= ringSize;
115  } // end of wrap-around
116  ringTriplet.push_back(ring[kIndex]);
117  } // next two elements of ringTriplet
118  //
119  // Obtained ringTriplet!
120  // Check equality
121  if (triplet == ringTriplet) {
122  return true;
123  } // triplet matches!
124  //
125  // Check the reversed triplet too
126  std::reverse(ringTriplet.begin(), ringTriplet.end());
127 
128  // Check for equality
129  if (triplet == ringTriplet) {
130  return true;
131  } // reversed triplet matches
132  } // first element of ringTriplet
133 
134  return false;
135 }
T clear(T... args)
Topological network criteria functions.
Definition: ring.hpp:60
T reverse(T... args)

◆ getSingleRingSize()

std::vector< std::vector< int > > ring::getSingleRingSize ( std::vector< std::vector< int >>  rings,
int  ringSize 
)

Returns a vector of vectors of rings of a single size.

Function that gets rings of a single ring size (i.e. a particular number of nodes) from all primitive rings, and returns a vector of vectors containing the rings of the specified size.

Parameters
[in]ringsThe vector of vectors containing the primitive rings of all sizes.
[in]ringSizeThe desired ring size or number of nodes in each ring to be saved.
Returns
A vector of vectors containing primitive rings of one ring size or length.

Definition at line 149 of file ring.cpp.

149  {
150  //
151  std::vector<std::vector<int>> ringSingleSize; // rings of one size
152 
153  // rings contains primitive rings of all sizes
154  // Only save rings of a given size (ringSize) to the new
155  // vector of vectors, ringSingleSize
156  for (int iring = 0; iring < rings.size(); iring++) {
157  // Check the size of the current ring
158  // If it is the correct size, save it in ringSingleSize
159  if (rings[iring].size() == ringSize) {
160  ringSingleSize.push_back(rings[iring]);
161  } // End of check of the size of iring
162  } // end of loop through all rings in rings
163 
164  return ringSingleSize;
165 }

◆ hasCommonElements()

bool ring::hasCommonElements ( std::vector< int >  ring1,
std::vector< int >  ring2 
)

Check to see if two vectors have common elements or not True, if common elements are present and false if there are no common elements

For two vectors, checks to see if there are common elements (true) or not (false).

Parameters
[in]ring1The vector of the first ring.
[in]ring2The vector of the second ring.
Returns
A bool which is true if the input vectors have at least one common element, and false if there are no common elements.

Definition at line 175 of file ring.cpp.

175  {
176  std::vector<int> commonElements; // Vector containing common elements
177 
178  // Sort the vectors before finding common elements
179  sort(ring1.begin(), ring1.end());
180  sort(ring2.begin(), ring2.end());
181 
182  // Find intersection of sorted vectors
183  auto it =
184  std::set_intersection(ring1.begin(), ring1.end(), ring2.begin(),
185  ring2.end(), std::back_inserter(commonElements));
186 
187  // If there are no elements in common, then return false
188  if (commonElements.size() == 0) {
189  return false;
190  }
191  // If there are common elements, return true
192  else {
193  return true;
194  }
195 }
T back_inserter(T... args)
T set_intersection(T... args)

◆ keepAxialRingsOnly()

std::vector<std::vector<int> > ring::keepAxialRingsOnly ( std::vector< std::vector< int >>  rings,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud 
)

Saves only axial rings out of all possible rings.

◆ polygonRingAnalysis()

int ring::polygonRingAnalysis ( std::string  path,
std::vector< std::vector< int >>  rings,
std::vector< std::vector< int >>  nList,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
int  maxDepth,
double  sheetArea,
int  firstFrame 
)

Find out which rings are prisms, looping through all ring sizes upto the maxDepth The input ringsAllSizes array has rings of every size.

Function that loops through the primitive rings (which is a vector of vectors) of all sizes, upto maxDepth (the largest ring size). The function returns a vector which contains the ring indices of all the rings which constitute prism blocks. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. This is registered as a Lua function, and is exposed to the user directly. The function calls the following functions internally:

  • ring::clearRingList (Clears the vector of vectors for rings of a single type, to prevent excessive memory being blocked).
  • ring::getSingleRingSize (Fill a vector of vectors for rings of a particular ring size).
  • ring::findPrisms (Now that rings of a particular size have been obtained, prism blocks are found and saved).
  • topoparam::normHeightPercent (Gets the height% for the prism blocks).
  • ring::assignPrismType (Assigns a prism type to each atom type).
  • sout::writePrismNum (Write out the prism information for the current frame).
  • sout::writeLAMMPSdataAllPrisms (Writes out a LAMMPS data file for the current frame, which can be visualized in OVITO).
    Parameters
    [in]pathThe string to the output directory, in which files will be written out.
    [in]ringsRow-ordered vector of vectors for rings of a single type.
    [in]nListRow-ordered neighbour list by index.
    [in]yCloudThe input PointCloud.
    [in]maxDepthThe maximum possible size of the primitive rings.
    [in]sheetAreaArea calculated using the two significant dimensions of the quasi-two-dimensional sheet.
    [in]firstFrameThe first frame to be analyzed

Definition at line 37 of file topo_two_dim.cpp.

41  {
42  //
44  ringsOneType; // Vector of vectors of rings of a single size
45  int nRings; // Number of rings of the current type
46  std::vector<int> nRingList; // Vector of the values of the number of prisms
47  // for a particular frame
48  std::vector<double> coverageAreaXY; // Coverage area percent on the XY plane
49  // for a particular n and frame
50  std::vector<double> coverageAreaXZ; // Coverage area percent on the XZ plane
51  // for a particular n and frame
52  std::vector<double> coverageAreaYZ; // Coverage area percent on the YZ plane
53  // for a particular n and frame
54  std::vector<double> singleAreas; // Coverage area on the XY, XZ and YZ planes
55  // for a single ring
57  atomTypes; // contains int values for each prism type considered
58  // -------------------------------------------------------------------------------
59  // Init
60  nRingList.resize(
61  maxDepth -
62  2); // Has a value for every value of ringSize from 3, upto maxDepth
63  coverageAreaXY.resize(maxDepth - 2);
64  coverageAreaXZ.resize(maxDepth - 2);
65  coverageAreaYZ.resize(maxDepth - 2);
66  // The atomTypes vector is the same size as the pointCloud atoms
67  atomTypes.resize(yCloud->nop, 1); // The dummy or unclassified value is 1
68  // -------------------------------------------------------------------------------
69  // Run this loop for rings of sizes upto maxDepth
70  // The smallest possible ring is of size 3
71  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
72  // Clear ringsOneType
73  ring::clearRingList(ringsOneType);
74  // Get rings of the current ring size
75  ringsOneType = ring::getSingleRingSize(rings, ringSize);
76  //
77  // Continue if there are zero rings of ringSize
78  if (ringsOneType.size() == 0) {
79  nRingList[ringSize - 3] = 0; // Update the number of prisms
80  // Update the coverage area%
81  coverageAreaXY[ringSize - 3] = 0.0;
82  coverageAreaXZ[ringSize - 3] = 0.0;
83  coverageAreaYZ[ringSize - 3] = 0.0;
84  continue;
85  } // skip if there are no rings
86  //
87  // -------------
88  // Number of rings with n nodes
89  nRings = ringsOneType.size();
90  // -------------
91  // Now that you have rings of a certain size:
92  nRingList[ringSize - 3] = nRings; // Update the number of n-membered rings
93  // Update the coverage area% for the n-membered ring phase.
94  singleAreas = topoparam::calcCoverageArea(yCloud, ringsOneType, sheetArea);
95  // XY plane
96  coverageAreaXY[ringSize - 3] = singleAreas[0];
97  // XZ plane
98  coverageAreaXZ[ringSize - 3] = singleAreas[1];
99  // YZ plane
100  coverageAreaYZ[ringSize - 3] = singleAreas[2];
101  // -------------
102  } // end of loop through every possible ringSize
103 
104  // Get the atom types for all the ring types
105  ring::assignPolygonType(rings, &atomTypes, nRingList);
106 
107  // Write out the ring information
108  sout::writeRingNum(path, yCloud->currentFrame, nRingList, coverageAreaXY,
109  coverageAreaXZ, coverageAreaYZ, maxDepth, firstFrame);
110  // Write out the lammps data file for the particular frame
111  sout::writeLAMMPSdataAllRings(yCloud, nList, atomTypes, maxDepth, path);
112 
113  return 0;
114 }
int nop
Current frame number.
Definition: mol_sys.hpp:169
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
Definition: ring.cpp:18
int assignPolygonType(std::vector< std::vector< int >> rings, std::vector< int > *atomTypes, std::vector< int > nRings)
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 writeRingNum(std::string path, int currentFrame, std::vector< int > nRings, std::vector< double > coverageAreaXY, std::vector< double > coverageAreaXZ, std::vector< double > coverageAreaYZ, int maxDepth, int firstFrame)
int writeLAMMPSdataAllRings(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< int > atomTypes, int maxDepth, std::string path)
Write a data file for rings of every type for a monolayer.
std::vector< double > calcCoverageArea(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> rings, double sheetArea)

◆ prismAnalysis()

int ring::prismAnalysis ( std::string  path,
std::vector< std::vector< int >>  rings,
std::vector< std::vector< int >>  nList,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
int  maxDepth,
int *  atomID,
int  firstFrame,
int  currentFrame,
bool  doShapeMatching = false 
)

Find out which rings are prisms, looping through all ring sizes upto the maxDepth The input ringsAllSizes array has rings of every size.

Function that loops through the primitive rings (which is a vector of vectors) of all sizes, upto maxDepth (the largest ring size). The function returns a vector which contains the ring indices of all the rings which constitute prism blocks. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. This is registered as a Lua function, and is exposed to the user directly. The function calls the following functions internally:

  • ring::clearRingList (Clears the vector of vectors for rings of a single type, to prevent excessive memory being blocked).
  • ring::getSingleRingSize (Fill a vector of vectors for rings of a particular ring size).
  • ring::findPrisms (Now that rings of a particular size have been obtained, prism blocks are found and saved).
  • topoparam::normHeightPercent (Gets the height% for the prism blocks).
  • ring::assignPrismType (Assigns a prism type to each atom type).
  • sout::writePrismNum (Write out the prism information for the current frame).
  • sout::writeLAMMPSdataAllPrisms (Writes out a LAMMPS data file for the current frame, which can be visualized in OVITO).
    Parameters
    [in]pathThe string to the output directory, in which files will be written out.
    [in]ringsRow-ordered vector of vectors for rings of a single type.
    [in]nListRow-ordered neighbour list by index.
    [in]yCloudThe input PointCloud.
    [in]maxDepthThe maximum possible size of the primitive rings.
    [in]maxDepthThe first frame.

Definition at line 45 of file topo_one_dim.cpp.

49  {
50  //
52  ringsOneType; // Vector of vectors of rings of a single size
53  std::vector<int> listPrism; // Vector for ring indices of n-sided prism
55  ringType; // This vector will have a value for each ring inside
57  atomState; // This vector will have a value for each atom, depending on
58  // the ring type
59  int nPerfectPrisms; // Number of perfect prisms of each type
60  int nImperfectPrisms; // Number of deformed prisms of each type
61  std::vector<int> nPrismList; // Vector of the values of the number of perfect
62  // prisms for a particular frame
63  std::vector<int> nDefPrismList; // Vector of the values of the number of
64  // deformed prisms for a particular frame
66  heightPercent; // Height percent for a particular n and frame
68  atomTypes; // contains int values for each prism type considered
69  double avgPrismHeight = 2.845; // A value of 2.7-2.85 Angstrom is reasonable
70  // Qualifier for the RMSD per atom:
71  std::vector<double> rmsdPerAtom;
72  // -------------------------------------------------------------------------------
73  // Init
74  nPrismList.resize(
75  maxDepth -
76  2); // Has a value for every value of ringSize from 3, upto maxDepth
77  nDefPrismList.resize(maxDepth - 2);
78  heightPercent.resize(maxDepth - 2);
79  // The atomTypes vector is the same size as the pointCloud atoms
80  atomTypes.resize(yCloud->nop, 1); // The dummy or unclassified value is 1
81  // The rmsdPerAtom vector is the same size as the pointCloud atoms and has an
82  // RMSD value for every atom
83  rmsdPerAtom.resize(yCloud->nop, -1);
84  // Resize the atom state vector
85  atomState.resize(yCloud->nop); // Dummy or unclassified
86  // -------------------------------------------------------------------------------
87  // Run this loop for rings of sizes upto maxDepth
88  // The smallest possible ring is of size 3
89  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
90  // Clear ringsOneType
91  ring::clearRingList(ringsOneType);
92  // Get rings of the current ring size
93  ringsOneType = ring::getSingleRingSize(rings, ringSize);
94  //
95  // Continue if there are zero rings of ringSize
96  if (ringsOneType.size() == 0) {
97  nPrismList[ringSize - 3] = 0; // Update the number of prisms
98  nDefPrismList[ringSize - 3] = 0; // Update the number of deformed prisms
99  heightPercent[ringSize - 3] = 0.0; // Update the height%
100  continue;
101  } // skip if there are no rings
102  //
103  // -------------
104  // Init of variables specific to ringSize prisms
105  listPrism.resize(0);
106  ringType.resize(0);
107  nPerfectPrisms = 0;
108  nImperfectPrisms = 0;
109  ringType.resize(
110  ringsOneType.size()); // Has a value for each ring. init to zero.
111  // -------------
112  // Now that you have rings of a certain size:
113  // Find prisms, saving the ring indices to listPrism
114  listPrism = ring::findPrisms(ringsOneType, &ringType, &nPerfectPrisms,
115  &nImperfectPrisms, nList, yCloud, &rmsdPerAtom,
116  doShapeMatching);
117  // -------------
118  nPrismList[ringSize - 3] =
119  nPerfectPrisms; // Update the number of perfect prisms
120  nDefPrismList[ringSize - 3] =
121  nImperfectPrisms; // Update the number of defromed prisms
122  int totalPrisms = nPerfectPrisms + nImperfectPrisms;
123  // Update the height% for the phase
124  heightPercent[ringSize - 3] =
125  topoparam::normHeightPercent(yCloud, totalPrisms, avgPrismHeight);
126  // Continue if there are no prism units
127  if (nPerfectPrisms + nImperfectPrisms == 0) {
128  continue;
129  } // skip for no prisms
130  // Do a bunch of write-outs and calculations
131  // TODO: Write out each individual prism as data files (maybe with an
132  // option)
133  // Get the atom types for a particular prism type
134  ring::assignPrismType(ringsOneType, listPrism, ringSize, ringType,
135  &atomTypes, &atomState);
136  // -------------
137  } // end of loop through every possible ringSize
138 
139  // Calculate the height%
140 
141  // Write out the prism information
142  sout::writePrismNum(path, nPrismList, nDefPrismList, heightPercent, maxDepth,
143  yCloud->currentFrame, firstFrame);
144  // Reassign the prism block types for the deformed prisms
145  if (doShapeMatching) {
146  ring::deformedPrismTypes(atomState, &atomTypes, maxDepth);
147  } // reassign prism block types for deformed prisms
148 
149  // Get rid of translations along the axial dimension
150  ring::rmAxialTranslations(yCloud, atomID, firstFrame, currentFrame);
151 
152  // Write out dump files with the RMSD per atom for the shape-matching
153  // criterion
154  if (doShapeMatching) {
155  sout::writeLAMMPSdumpINT(yCloud, rmsdPerAtom, atomTypes, maxDepth, path);
156  } // reassign prism block types for deformed prisms
157 
158  // Write out the lammps data file for the particular frame
159  sout::writeLAMMPSdataAllPrisms(yCloud, nList, atomTypes, maxDepth, path,
160  doShapeMatching);
161 
162  return 0;
163 }
std::vector< int > findPrisms(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, int *nPerfectPrisms, int *nImperfectPrisms, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, bool doShapeMatching=false)
int deformedPrismTypes(std::vector< ring::strucType > atomState, std::vector< int > *atomTypes, int maxDepth)
Get the atom type values for deformed prisms.
int rmAxialTranslations(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int *atomID, int firstFrame, int currentFrame)
Shift the entire ice nanotube and remove axial translations.
int assignPrismType(std::vector< std::vector< int >> rings, std::vector< int > listPrism, int ringSize, std::vector< ring::strucType > ringType, std::vector< int > *atomTypes, std::vector< ring::strucType > *atomState)
int writePrismNum(std::string path, std::vector< int > nPrisms, std::vector< int > nDefPrisms, std::vector< double > heightPercent, int maxDepth, int currentFrame, int firstFrame)
int writeLAMMPSdumpINT(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > rmsdPerAtom, std::vector< int > atomTypes, int maxDepth, std::string path)
Write out a LAMMPS dump file containing the RMSD per atom.
int writeLAMMPSdataAllPrisms(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< int > atomTypes, int maxDepth, std::string path, bool doShapeMatching=false)
Write a data file for prisms of every type.
double normHeightPercent(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int nPrisms, double avgPrismHeight)

◆ relaxedPrismConditions()

bool ring::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 432 of file topo_one_dim.cpp.

434  {
435  int ringSize =
436  (*basal1).size(); // Size of the ring; each ring contains n elements
437  int m_k; // Atom ID of element in basal2
438  bool isNeighbour = false; // This is true if there is at least one bond
439  // between the basal rings
440  int l_k; // Atom ID of element in basal1
441 
442  // ---------------------------------------------
443  // COMPARISON OF basal2 ELEMENTS (m_k) WITH basal1 ELEMENTS (l_k)
444  // Loop through all the elements of basal1
445  for (int l = 0; l < ringSize; l++) {
446  l_k = (*basal1)[l];
447  // Search for the nearest neighbour of l_k in basal2
448  // Loop through basal2 elements
449  for (int m = 0; m < ringSize; m++) {
450  m_k = (*basal2)[m];
451  // Find m_k inside l_k neighbour list
452  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
453 
454  // If the element has been found, for l1
455  if (it != nList[l_k].end()) {
456  isNeighbour = true;
457  break;
458  } // found element
459  } // end of loop through all the elements of basal2
460 
461  // If a neighbour has been found then
462  if (isNeighbour) {
463  return true;
464  }
465  } // end of loop through all the elements of basal1
466 
467  // If a neighbour has not been found, return false
468  return false;
469 }

◆ rmAxialTranslations()

int ring::rmAxialTranslations ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
int *  atomID,
int  firstFrame,
int  currentFrame 
)

Shift the entire ice nanotube and remove axial translations.

Assign an atomType value for atoms belonging to deformed prisms.

Parameters
[in]ringsThe vector of vectors containing the primitive rings, of a particular ring size.
[in]listPrismThe list of prism blocks found.
[in]ringSizeThe current ring size or number of nodes in each ring.
[in,out]atomTypesA vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.

Definition at line 787 of file topo_one_dim.cpp.

789  {
790  //
791  int atomIndex; // Index of the atom to be shifted first
792  double shiftDistance; // Value by which all atoms will be shifted
793  double distFrmBoundary; // Distance from either the upper or lower boundary
794  // along the axial dimension
795  // Get the axial dimension
796  int axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
797  yCloud->box.begin();
798  // Check: (default z)
799  if (axialDim < 0 || axialDim > 2) {
800  axialDim = 2;
801  } // end of check to set the axial dimension
802  // Lower and higher limits of the box in the axial dimension
803  double boxLowAxial = yCloud->boxLow[axialDim];
804  double boxHighAxial = boxLowAxial + yCloud->box[axialDim];
805  //
806  // -----------------------------------
807  // Update the atomID of the 'first' atom in yCloud if the current frame is the
808  // first frame.
809  // Get the atom index of the first atom to be shifted down or up
810  if (currentFrame == firstFrame) {
811  *atomID = yCloud->pts[0].atomID;
812  atomIndex = 0;
813  } // end of update of the atom ID to be shifted for the first frame
814  else {
815  //
816  int iatomID = *atomID; // Atom ID of the first particle to be moved down
817  // Find the index given the atom ID
818  auto it = yCloud->idIndexMap.find(iatomID);
819 
820  if (it != yCloud->idIndexMap.end()) {
821  atomIndex = it->second;
822  } // found jatom
823  else {
824  std::cerr << "Lost atoms.\n";
825  return 1;
826  } // error handling
827  //
828  } // Update of atomIndex for all other frames
829  // -----------------------------------
830  // Calculate the distance by which the atoms must be shifted (negative value)
831  if (axialDim == 0) {
832  shiftDistance = boxLowAxial - yCloud->pts[atomIndex].x;
833  } // x coordinate
834  else if (axialDim == 1) {
835  shiftDistance = boxLowAxial - yCloud->pts[atomIndex].y;
836  } // y coordinate
837  else {
838  shiftDistance = boxLowAxial - yCloud->pts[atomIndex].z;
839  } // z coordinate
840  // -----------------------------------
841  // Shift all the particles
842  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
843  // Shift the particles along the axial dimension only
844  if (axialDim == 0) {
845  yCloud->pts[iatom].x += shiftDistance;
846  // Wrap-around
847  if (yCloud->pts[iatom].x < boxLowAxial) {
848  distFrmBoundary = boxLowAxial - yCloud->pts[iatom].x; // positive value
849  yCloud->pts[iatom].x = boxHighAxial - distFrmBoundary;
850  } // end of wrap-around
851  } // x coordinate
852  else if (axialDim == 1) {
853  yCloud->pts[iatom].y += shiftDistance;
854  // Wrap-around
855  if (yCloud->pts[iatom].y < boxLowAxial) {
856  distFrmBoundary = boxLowAxial - yCloud->pts[iatom].y; // positive value
857  yCloud->pts[iatom].y = boxHighAxial - distFrmBoundary;
858  } // end of wrap-around
859  } // y coordinate
860  else {
861  yCloud->pts[iatom].z += shiftDistance;
862  // Wrap-around
863  if (yCloud->pts[iatom].z < boxLowAxial) {
864  distFrmBoundary = boxLowAxial - yCloud->pts[iatom].z; // positive value
865  yCloud->pts[iatom].z = boxHighAxial - distFrmBoundary;
866  } // end of wrap-around
867  } // z coordinate
868  } // end of shifting all paritcles downward
869  // -----------------------------------
870  return 0;
871 } // end of function
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:172
std::vector< T > boxLow
Periodic box lengths.
Definition: mol_sys.hpp:171

◆ shapeMatchDDC()

int tum3::shapeMatchDDC ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
const Eigen::MatrixXd &  refPoints,
std::vector< cage::Cage cageList,
int  cageIndex,
std::vector< std::vector< int >>  rings,
std::vector< double > *  quat,
double *  rmsd 
)

Shape-matching for a target DDC.

Match the input cage with a perfect DDC. The quaternion for the rotation and the RMSD is outputted.

Parameters
[in]refPointsThe point set of the reference system (or right system). This is a \( (n \times 3) \) Eigen matrix. Here, \( n \) is the number of particles.
[in]targetPoints\( (n \times 3) \) Eigen matrix of the candidate/test system (or left system).
[in,out]quatThe quaternion for the optimum rotation of the left system into the right system.
[in,out]scaleThe scale factor obtained from Horn's algorithm.

Definition at line 319 of file bulkTUM.cpp.

323  {
324  //
325  std::vector<int> ddcOrder; // Connectivity of the DDC
326  int ringSize = 6; // Each ring has 6 nodes
327  Eigen::MatrixXd targetPointSet(14, 3); // Target point set (Eigen matrix)
328  //
329  std::vector<double> rmsdList; // List of RMSD per atom
330  // Variables for looping through possible permutations
331  //
332  std::vector<double> currentQuat; // quaternion rotation
333  double currentRmsd; // least RMSD
334  double currentScale; // scale
335  double scale; // Final scale
336  int index; // Int for describing which permutation matches the best
337 
338  // ----------------
339  // Save the order of the DDC in a vector
340  ddcOrder = pntToPnt::relOrderDDC(cageIndex, rings, cageList);
341  // ----------------
342  // Loop through all possible permutations
343  //
344  for (int i = 0; i < ringSize; i++) {
345  // Change the order of the target points somehow!
346  //
347  targetPointSet = pntToPnt::changeDiaCageOrder(yCloud, ddcOrder, i);
348  // Shape-matching
349  absor::hornAbsOrientation(refPoints, targetPointSet, &currentQuat,
350  &currentRmsd, &rmsdList, &currentScale);
351  if (i == 0) {
352  *quat = currentQuat;
353  *rmsd = currentRmsd;
354  scale = currentScale;
355  index = 0;
356  } else {
357  if (currentRmsd < *rmsd) {
358  *quat = currentQuat;
359  *rmsd = currentRmsd;
360  scale = currentScale;
361  index = i;
362  } // update
363  } // Update if this is a better match
364  } // Loop through possible permutations
365  // ----------------
366  return 0;
367 }
int hornAbsOrientation(const Eigen::MatrixXd &refPoints, const Eigen::MatrixXd &targetPoints, std::vector< double > *quat, double *rmsd, std::vector< double > *rmsdList, double *scale)
Get the absolute orientation using Horn's algorithm (with quaternions)

◆ shapeMatchHC()

int tum3::shapeMatchHC ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
const Eigen::MatrixXd &  refPoints,
cage::Cage  cageUnit,
std::vector< std::vector< int >>  rings,
std::vector< std::vector< int >>  nList,
std::vector< double > *  quat,
double *  rmsd 
)

Shape-matching for a target HC.

Match the input cage with a perfect HC. The quaternion for the rotation and the RMSD is outputted.

Parameters
[in]refPointsThe point set of the reference system (or right system). This is a \( (n \times 3) \) Eigen matrix. Here, \( n \) is the number of particles.
[in]targetPoints\( (n \times 3) \) Eigen matrix of the candidate/test system (or left system).
[in,out]quatThe quaternion for the optimum rotation of the left system into the right system.
[in,out]scaleThe scale factor obtained from Horn's algorithm.

Definition at line 247 of file bulkTUM.cpp.

251  {
252  //
253  int iring,
254  jring; // Indices in the ring vector of vectors for the basal rings
255  std::vector<int> basal1,
256  basal2; // Re-ordered basal rings matched to each other
257  int ringSize = 6; // Each ring has 6 nodes
258  Eigen::MatrixXd targetPointSet(12, 3); // Target point set (Eigen matrix)
259  //
260  std::vector<double> rmsdList; // List of RMSD per atom
261  // Variables for looping through possible permutations
262  //
263  std::vector<double> currentQuat; // quaternion rotation
264  double currentRmsd; // least RMSD
265  double currentScale; // scale
266  double scale; // Final scale
267  int index; // Int for describing which permutation matches the best
268 
269  // Init
270  iring = cageUnit.rings[0]; // Index of basal1
271  jring = cageUnit.rings[1]; // Index of basal2
272  rmsdList.resize(yCloud->nop); // Not actually updated here
273  //
274  // ----------------
275  // Re-order the basal rings so that they are matched
276  pntToPnt::relOrderHC(yCloud, rings[iring], rings[jring], nList, &basal1,
277  &basal2);
278  // ----------------
279  // Loop through all possible permutations
280  //
281  for (int i = 0; i < ringSize; i++) {
282  // Change the order of the target points somehow!
283  //
284  targetPointSet = pntToPnt::changeHexCageOrder(yCloud, basal1, basal2, i);
285  // Shape-matching
286  absor::hornAbsOrientation(refPoints, targetPointSet, &currentQuat,
287  &currentRmsd, &rmsdList, &currentScale);
288  if (i == 0) {
289  *quat = currentQuat;
290  *rmsd = currentRmsd;
291  scale = currentScale;
292  index = 0;
293  } else {
294  if (currentRmsd < *rmsd) {
295  *quat = currentQuat;
296  *rmsd = currentRmsd;
297  scale = currentScale;
298  index = i;
299  } // update
300  } // Update if this is a better match
301  } // Loop through possible permutations
302  // ----------------
303  return 0;
304 }
std::vector< int > rings
type of the cage : can be DDC or HC
Definition: cage.hpp:84

◆ topoBulkCriteria()

std::vector< cage::Cage > tum3::topoBulkCriteria ( std::string  path,
std::vector< std::vector< int >>  rings,
std::vector< std::vector< int >>  nList,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
int  firstFrame,
int *  numHC,
int *  numDDC,
std::vector< ring::strucType > *  ringType 
)

Topological network methods Finds the HCs and DDCs for the system

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 only for rings of size 6! 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 should contain rings of only size 6.
    [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 577 of file bulkTUM.cpp.

581  {
582  //
583  // Ring IDs of each type will be saved in these vectors
584  std::vector<int> listDDC; // Vector for ring indices of DDC
585  std::vector<int> listHC; // Vector for ring indices of HC
587  listMixed; // Vector for ring indices of rings that are both DDC and HC
588  // ringList, of type enum strucType in gen.hpp
589  // Make a list of all the DDCs and HCs
590  std::vector<cage::Cage> cageList;
591  // Number of types
592  int mixedRings, prismaticRings, basalRings;
593 
594  // ----------------------------------------------
595  // Init
596  //
597  *numHC = 0; // Number of hexagonal cages
598  *numDDC = 0; // Init the number of DDCs
599  // Quit the function for zero rings
600  if (rings.size() == 0) {
601  return cageList;
602  } // skip for no rings of ringSize
603  //
604  // Init the ringType vector
605  (*ringType).resize(rings.size()); // Has a value for each ring. init to zero.
606  // ----------------------------------------------
607  // Get the cages
608 
609  // Find HC rings, saving the ring IDs (starting from 0) to listHC
610  listHC = ring::findHC(rings, ringType, nList, &cageList);
611 
612  // Find DDC rings, saving the IDs to listDDC
613  listDDC = ring::findDDC(rings, ringType, listHC, &cageList);
614 
615  // Find rings which are both DDCs and HCs (mixed)
616  // A dummy value of -10 in the listDDC and listHC vectors for mixed rings
617  listMixed = ring::findMixedRings(rings, ringType, &listDDC, &listHC);
618 
619  // Get the number of structures (DDCs, HCs, mixed rings, basal rings,
620  // prismatic rings)
621  ring::getStrucNumbers(*ringType, cageList, numHC, numDDC, &mixedRings,
622  &prismaticRings, &basalRings);
623 
624  // Write out to a file
625  sout::writeTopoBulkData(path, yCloud->currentFrame, *numHC, *numDDC,
626  mixedRings, basalRings, prismaticRings, firstFrame);
627 
628  return cageList;
629 }
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
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
int writeTopoBulkData(std::string path, int currentFrame, int numHC, int numDDC, int mixedRings, int basalRings, int prismaticRings, int firstFrame)

◆ topoUnitMatchingBulk()

int tum3::topoUnitMatchingBulk ( 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  printClusters,
bool  onlyTetrahedral 
)

Topological unit matching for bulk water. If printClusters is true, individual clusters of connected cages are printed.

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 bulkTUM.cpp.

52  {
53  // The input rings vector has rings of all sizes
54 
55  // ringType has a value for rings of a particular size
57  ringType; // This vector will have a value for each ring inside
58  // ringList, of type enum strucType in gen.hpp
59  // Make a list of all the DDCs and HCs
60  std::vector<cage::Cage> cageList;
62  ringsOneType; // Vector of vectors of rings of a single size
63  int initRingSize; // Todo or not: calculate the PNCs or not
64  int maxRingSize = 6; // DDCs and HCs are for 6-membered rings
66  atomTypes; // This vector will have a value for every atom
67  // Number of types
68  int numHC, numDDC, mixedRings, prismaticRings, basalRings;
69  // Shape-matching variables ----
70  double rmsd; // RMSD value for a particular cage type
71  std::vector<double> quat; // Quaternion obtained from shape-matching
72  std::vector<std::vector<double>> quatList; // List of quaternions
73  // Reference points
74  Eigen::MatrixXd refPntsDDC(14,
75  3); // Reference point set (Eigen matrix) for a DDC
76  Eigen::MatrixXd refPntsHC(12,
77  3); // Reference point set (Eigen matrix) for a DDC
78 
79  // Vector for the RMSD per atom and the RMSD per ring:
80  std::vector<double> rmsdPerAtom;
82  noOfCommonElements; // An atom can belong to more than one ring or shape
83  //
84 
85  if (onlyTetrahedral) {
86  initRingSize = 6;
87  } else {
88  initRingSize = 5;
89  }
90 
91  // Init the atom type vector
92  atomTypes.resize(yCloud->nop); // Has a value for each atom
93  // Init the rmsd per atom
94  rmsdPerAtom.resize(yCloud->nop, -1); // Has a value for each atom
95  noOfCommonElements.resize(yCloud->nop); // Has a value for each atom
96 
97  // ----------------------------------------------
98  // Init
99  // Get rings of size 5 or 6.
100  for (int ringSize = initRingSize; ringSize <= maxRingSize; ringSize++) {
101  // Clear ringsOneType
102  ring::clearRingList(ringsOneType);
103  // Get rings of the current ring size
104  ringsOneType = ring::getSingleRingSize(rings, ringSize);
105  // Skip for zero rings
106  if (ringsOneType.size() == 0) {
107  continue;
108  } // skip for no rings of ringSize
109  //
110  // Init the ringType vector
111  ringType.resize(
112  ringsOneType.size()); // Has a value for each ring. init to zero.
113  // ----------------------------------------------
114  if (ringSize == 6) {
115  // Get the cages
116 
117  // Find the DDCs and HCs
118  cageList = tum3::topoBulkCriteria(path, ringsOneType, nList, yCloud,
119  firstFrame, &numHC, &numDDC, &ringType);
120 
121  // Gets the atom type for every atom, to be used for printing out the ice
122  // types found
123  ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
124  }
125  // Finding the 5-membered prismatic blocks
126  else {
127  // Get the prism block classifications
128  prism3::findBulkPrisms(ringsOneType, &ringType, nList, yCloud,
129  &rmsdPerAtom);
130  // Gets the atom type for every atom, to be used for printing out the ice
131  // types found
132  ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
133  }
134  }
135 
136  // --------------------------------------------------
137  // SHAPE-MATCHING
138  //
139  // Get the reference point sets
140  //
141  std::string filePathHC = "templates/hc.xyz";
142  std::string filePathDDC = "templates/ddc.xyz";
143  refPntsHC = tum3::buildRefHC(filePathHC); // HC
144  refPntsDDC = tum3::buildRefDDC(filePathDDC); // DDC
145  //
146  // Init
147  //
148  // Loop through all the atoms and set commonElements to 1 for PNC atoms
149  for (int i = 0; i < yCloud->nop; i++) {
150  if (rmsdPerAtom[i] != -1) {
151  noOfCommonElements[i] = 1;
152  } // for atoms which have been updated
153  } // end of updating commonElements
154  //
155  // Loop through the entire cageList vector of cages to match the HCs and DDCs
156  //
157  // HCs are first, followed by DDCs.
158  //
159  // Go through all the HCs
160  for (int icage = 0; icage < numHC; icage++) {
161  // Match against a perfect HC
162  tum3::shapeMatchHC(yCloud, refPntsHC, cageList[icage], ringsOneType, nList,
163  &quat, &rmsd);
164  // Update the vector of quaternions
165  quatList.push_back(quat);
166  // Update the RMSD per ring
167  tum3::updateRMSDatom(ringsOneType, cageList[icage], rmsd, &rmsdPerAtom,
168  &noOfCommonElements, atomTypes);
169  } // end of looping through all HCs
170  // --------------------------------------------------
171  // Go through all the DDCs
172  for (int icage = numHC; icage < cageList.size(); icage++) {
173  // Match against a perfect DDC
174  tum3::shapeMatchDDC(yCloud, refPntsDDC, cageList, icage, ringsOneType,
175  &quat, &rmsd);
176  // Update the vector of quaternions
177  quatList.push_back(quat);
178  // Update the RMSD per ring
179  tum3::updateRMSDatom(ringsOneType, cageList[icage], rmsd, &rmsdPerAtom,
180  &noOfCommonElements, atomTypes);
181  } // end of looping through all HCs
182 
183  // --------------------------------------------------
184  // Getting the RMSD per atom
185  // Average the RMSD per atom
186  tum3::averageRMSDatom(&rmsdPerAtom, &noOfCommonElements);
187  // --------------------------------------------------
188 
189  // Print out the lammps data file with the bonds and types
190  sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path);
191  // To output the bonds between dummy atoms, uncomment the following line
192  // sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path, true);
193 
194  // --------------------------------------------------
195  // Writes out the dumpfile
196  //
197  std::vector<int> dumpAtomTypes; // Must be typecast to int
198  dumpAtomTypes.resize(atomTypes.size());
199  // Change from enum to int
200  for (int i = 0; i < atomTypes.size(); i++) {
201  if (atomTypes[i] == cage::hc) {
202  dumpAtomTypes[i] = 1;
203  } // HC
204  else if (atomTypes[i] == cage::ddc) {
205  dumpAtomTypes[i] = 2;
206  } // DDC
207  else if (atomTypes[i] == cage::mixed) {
208  dumpAtomTypes[i] = 3;
209  } // mixed rings
210  else if (atomTypes[i] == cage::pnc) {
211  dumpAtomTypes[i] = 4;
212  } // pnc
213  else if (atomTypes[i] == cage::mixed2) {
214  dumpAtomTypes[i] = 5;
215  } // shared by pnc and ddc/hc
216  else {
217  dumpAtomTypes[i] = 0;
218  } // dummy atoms
219  } // end of changing the type to int
220  sout::writeLAMMPSdumpCages(yCloud, rmsdPerAtom, dumpAtomTypes, path,
221  firstFrame);
222  // --------------------------------------------------
223  //
224  // PRINT THE CLUSTERS
225  //
226  if (printClusters) {
227  // Print the clusters
228  tum3::clusterCages(yCloud, path, ringsOneType, cageList, numHC, numDDC);
229  } // end of printing individual clusters
230  // --------------------------------------------------
231  return 0;
232 }
@ 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
int getAtomTypesTopoBulk(std::vector< std::vector< int >> rings, std::vector< ring::strucType > ringType, std::vector< cage::iceType > *atomTypes)
Definition: topo_bulk.cpp:1078
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
int shapeMatchHC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, cage::Cage cageUnit, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, std::vector< double > *quat, double *rmsd)
Shape-matching for a target HC.
Definition: bulkTUM.cpp:247
Eigen::MatrixXd buildRefHC(std::string fileName)
Build a reference Hexagonal cage, reading in from a template XYZ file.
Definition: bulkTUM.cpp:373
int averageRMSDatom(std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms)
Average the RMSD per atom.
Definition: bulkTUM.cpp:525
int shapeMatchDDC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, std::vector< cage::Cage > cageList, int cageIndex, std::vector< std::vector< int >> rings, std::vector< double > *quat, double *rmsd)
Shape-matching for a target DDC.
Definition: bulkTUM.cpp:319
int updateRMSDatom(std::vector< std::vector< int >> rings, cage::Cage cageUnit, double rmsd, std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms, std::vector< cage::iceType > atomTypes)
Definition: bulkTUM.cpp:482
std::vector< cage::Cage > topoBulkCriteria(std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, int *numHC, int *numDDC, std::vector< ring::strucType > *ringType)
Definition: bulkTUM.cpp:577
Eigen::MatrixXd buildRefDDC(std::string fileName)
Build a reference Double-Diamond cage, reading in from a template XYZ file.
Definition: bulkTUM.cpp:431
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)
int writeLAMMPSdumpCages(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > rmsdPerAtom, std::vector< int > atomTypes, std::string path, int firstFrame)
Write out a LAMMPS dump file containing the RMSD per atom for bulk ice.

◆ updateRMSDatom()

int tum3::updateRMSDatom ( std::vector< std::vector< int >>  rings,
cage::Cage  cageUnit,
double  rmsd,
std::vector< double > *  rmsdPerAtom,
std::vector< int > *  noOfCommonAtoms,
std::vector< cage::iceType atomTypes 
)

Calulate the RMSD for each ring, using RMSD values (rmsd) obtained from the shape-matching of each cage

Update the calculated RMSD per ring using the RMSD values of each cage, and also update the values in the noOfCommonRings vector, which will be used for averaging the RMSD per atom depending on the number of cages that share that particular ring.

Definition at line 482 of file bulkTUM.cpp.

486  {
487  //
488  int nRings = cageUnit.rings.size(); // Number of rings in the current cage
489  int iring; // Index according to the rings vector of vector, for the current
490  // ring inside the cage
491  int iatom; // Current atom index
492  int ringSize = rings[0].size(); // Number of nodes in each ring (6)
493 
494  // Loop through the rings in each cage
495  for (int i = 0; i < nRings; i++) {
496  iring = cageUnit.rings[i]; // Current ring
497 
498  // Loop through all the atoms in the ring
499  for (int j = 0; j < ringSize; j++) {
500  iatom = rings[iring][j]; // Current atom index
501 
502  // Skip for PNC atoms
503  if (atomTypes[iatom] == cage::pnc || atomTypes[iatom] == cage::mixed2) {
504  continue;
505  } // Do not update if the atom is a PNC
506  //
507  // UPDATE
508  if ((*rmsdPerAtom)[iatom] == -1) {
509  (*rmsdPerAtom)[iatom] = rmsd;
510  (*noOfCommonAtoms)[iatom] = 1;
511  } else {
512  (*rmsdPerAtom)[iatom] += rmsd;
513  (*noOfCommonAtoms)[iatom] += 1;
514  }
515  } // end of looping through all the atoms in the ring
516 
517  } // end of looping through the rings in the cage
518 
519  return 0;
520 }