Ring
+ Collaboration diagram for Ring:

Modules

 Prism3
 

Enumerations

enum class  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...
 
int ring::assignPolygonType (std::vector< std::vector< int >> rings, std::vector< int > *atomTypes, std::vector< int > nRings)
 
molSys::PointCloud< molSys::Point< double >, double > gen::getPointCloudOneAtomType (molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *outCloud, int atomTypeI, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
 
void gen::moleculesInSingleSlice (molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool clearPreviousSliceSelection=true, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
 
void gen::atomsInSingleSlice (molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool clearPreviousSliceSelection=true, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
 
void gen::setAtomsWithSameMolID (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::unordered_multimap< int, int > molIDAtomIDmap, int molID, bool inSliceValue=true)
 
void ring::getEdgeMoleculesInRings (std::vector< std::vector< int >> rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
 
void ring::printSliceGetEdgeMoleculesInRings (std::string path, std::vector< std::vector< int >> rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
 
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)
 

Detailed Description

Enumeration Type Documentation

◆ strucType

enum ring::strucType
strong
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 116 of file ring.hpp.

116  {
117  unclassified,
118  DDC,
119  HCbasal,
120  HCprismatic,
121  bothBasal,
123  Prism,
126 };
@ unclassified
Not classified into any other category.
@ HCbasal
The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an H...
@ DDC
The ring belongs to a double-diamond cage (DDC).
@ bothPrismatic
The ring belongs to both a DDC and HC and is, thus, a 'mixed' ring. The ring is also one of the prism...
@ Prism
The ring belongs to a prism block, classified according to the prism identification scheme.
@ bothBasal
The ring belongs to both a DDC and HC. It is a 'mixed' ring. The ring is also one of the basal rings ...
@ HCprismatic
The ring belongs only to a hexagonal cage (HC); specifically the ring is purely a prismatic ring of a...

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 40 of file ring.cpp.

42  {
43  // Every value in listPrism corresponds to an index in rings.
44  // Every ring contains atom indices, corresponding to the indices (not atom
45  // IDs) in rings
46  int iring; // Index of current ring
47  int iatom; // Index of current atom
48  int ringSize; // Ring size of the current ring
49  int prevRingSize; // Ring size previously assigned to a point
50 
51  // Dummy value corresponds to a value of 1.
52  // If an atom is shared by more than one ring type, it is assigned the
53  // value 2.
54 
55  // Loop through every ring in rings
56  for (int iring = 0; iring < rings.size(); iring++) {
57  ringSize = rings[iring].size();
58  // Loop through every element in iring
59  for (int j = 0; j < ringSize; j++) {
60  iatom = rings[iring][j]; // Atom index
61  // Update the atom type
62  if ((*atomTypes)[iatom] == 1) {
63  (*atomTypes)[iatom] = ringSize;
64  } // The atom is unclassified
65  else {
66  // Only update the ring type if the number is higher
67  prevRingSize = (*atomTypes)[iatom]; // Previously assigned ring size
68  if (ringSize > prevRingSize) {
69  (*atomTypes)[iatom] = ringSize;
70  } // end of assigning the new ring size
71  } // only update if the number is higher
72  } // end of loop through every atom in iring
73  } // end of loop through every ring
74 
75  return 0;
76 } // 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 696 of file topo_one_dim.cpp.

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

◆ 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 801 of file bulkTUM.cpp.

803  {
804  //
805  std::vector<int> atoms; // Contains the atom indices (not IDs) of atoms
806  int ringSize = rings[0].size(); // Number of nodes in each ring
807  int iring; // Index of ring in rings vector of vectors
808  int icage; // Index of the current cage in cageList
809  int iatom; // Atom index
810  //
811  for (int i = 0; i < clusterCages.size(); i++) {
812  //
813  icage = clusterCages[i];
814  // Loop through every ring in the cage
815  for (int j = 0; j < cageList[icage].rings.size(); j++) {
816  iring = cageList[icage].rings[j]; // Current ring index
817  // Loop through every atom in iring
818  for (int k = 0; k < ringSize; k++) {
819  iatom = rings[iring][k];
820  atoms.push_back(iatom); // Add the atom
821  } // Loop through every atom in iring
822  } // loop through every ring in the current cage
823  } // end of loop through all cages
824 
825  // Remove the duplicates in the atoms vector
826  // Duplicate IDs must be removed
827  //
828  sort(atoms.begin(), atoms.end());
829  auto ip = std::unique(atoms.begin(), atoms.end());
830  // Resize peripheral rings to remove undefined terms
831  atoms.resize(std::distance(atoms.begin(), ip));
832  //
833  // Return the vector of atoms
834  return atoms;
835 }
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:644
T push_back(T... args)
T resize(T... args)
T sort(T... args)
T unique(T... args)

◆ atomsInSingleSlice()

void gen::atomsInSingleSlice ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
bool  clearPreviousSliceSelection = true,
std::array< double, 3 >  coordLow = std::array<double, 3>{0, 0, 0},
std::array< double, 3 >  coordHigh = std::array<double, 3>{0, 0, 0} 
)

Given a pointCloud set the inSlice bool for every atom, if the atoms are inside the specified (single) region. Does not handle atoms in molecules straddling the boundary

Function that loops through a given input PointCloud and sets the inSlice bool for every Point according to whether the atom
is in the specified (single) slice or not. Not inclusive of atoms in molecules

Parameters
[in]yCloudThe given input PointCloud
[in]clearPreviousSliceSelectionsets all inSlice bool values to false before adding Points to the slice
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 101 of file selection.cpp.

105  {
106  //
107  bool pointInSlice = true; // If the current point is in the slice, this is true (otherwise this is false)
108 
109  // --------------------
110  // Set inSlice to false for every atom first
111  if (clearPreviousSliceSelection)
112  {
113  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
114  yCloud->pts[iatom].inSlice = false;
115  }
116  } // end of init
117  // --------------------
118 
119  // Loop through every iatom and check if the atom is in the slice or not
120  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
121  // Check if the current atom is in the slice
122  pointInSlice = sinp::atomInSlice(yCloud->pts[iatom].x, yCloud->pts[iatom].y, yCloud->pts[iatom].z,
123  coordLow, coordHigh);
124 
125  yCloud->pts[iatom].inSlice = pointInSlice; // iatom is inside the slice
126  //
127  }
128 
129  return;
130 }
int nop
Current frame number.
Definition: mol_sys.hpp:173
std::vector< S > pts
Definition: mol_sys.hpp:171
bool atomInSlice(double x, double y, double z, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh)
Definition: seams_input.hpp:88

◆ 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 529 of file bulkTUM.cpp.

530  {
531  //
532  int nop = (*rmsdPerAtom).size(); // Number of particles
533 
534  for (int iatom = 0; iatom < nop; iatom++) {
535  //
536  if ((*noOfCommonAtoms)[iatom] == 0) {
537  (*rmsdPerAtom)[iatom] = -1; // Dummy atom
538  } else {
539  (*rmsdPerAtom)[iatom] /= (*noOfCommonAtoms)[iatom];
540  }
541  } // end of averaging
542 
543  return 0;
544 } // 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 351 of file topo_one_dim.cpp.

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

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

◆ 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 377 of file bulkTUM.cpp.

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

22  {
23  //
25 
26  rings.swap(tempEmpty);
27 
28  return 0;
29 }
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 644 of file bulkTUM.cpp.

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

116  {
118  common1; // Vector containing the common elements of the first two rings
120  common2; // Vector containing the common elements of the three rings
121 
122  // Common elements among the first two rings
123  common1 = ring::findsCommonElements(ring1, ring2);
124  if (common1.size() == 0) {
125  return false;
126  } // no common elements in ring1 and ring2
127 
128  // Common elements among all three rings
129  common2 = ring::findsCommonElements(common1, ring3);
130 
131  // If no common elements were found:
132  if (common2.size() == 0) {
133  return false;
134  } // no common elements between ring1, ring2, and ring3
135 
136  return true; // Common elements found!
137 }
std::vector< int > findsCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Returns the common elements of two rings.
Definition: ring.cpp:85

◆ 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 258 of file ring.cpp.

258  {
259  // Sort the rings first
260  sort(ring1.begin(), ring1.end()); // Sort ring1 by ID
261  sort(ring2.begin(), ring2.end()); // Sort ring2 by ID
262  bool result;
263 
264  (ring1 == ring2) ? result = true : result = false;
265 
266  return result;
267 }

◆ 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 759 of file topo_one_dim.cpp.

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

◆ 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 487 of file topo_one_dim.cpp.

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

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

86  {
87  //
88  std::vector<int> common;
89  int iatom; // Index to search for
90 
91  for (int i = 0; i < ring1.size(); i++) {
92  iatom = ring1[i];
93  // Search for iatom in ring2
94 
95  auto it = std::find(ring2.begin(), ring2.end(), iatom);
96 
97  if (it != ring2.end()) {
98  common.push_back(iatom);
99  } // iatom was found!
100  } // end of loop through every element of ring1
101 
102  return common;
103 }

◆ 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 148 of file ring.cpp.

148  {
149  //
150  int ringSize = ring.size(); // should be 6
151  std::vector<int> ringTriplet; // triplet from the ring to be searched
152  int kIndex; // Used for making the triplet
153 
154  // Loop through every possible triplet in the ring to be searched
155  for (int i = 0; i < ringSize; i++) {
156  ringTriplet.clear(); // init
157  // Get the first element of the ring triplet
158  ringTriplet.push_back(ring[i]);
159  //
160  // Get the next two elements
161  for (int k = 1; k < 3; k++) {
162  kIndex = i + k;
163  // Wrap-around
164  if (kIndex >= ringSize) {
165  kIndex -= ringSize;
166  } // end of wrap-around
167  ringTriplet.push_back(ring[kIndex]);
168  } // next two elements of ringTriplet
169  //
170  // Obtained ringTriplet!
171  // Check equality
172  if (triplet == ringTriplet) {
173  return true;
174  } // triplet matches!
175  //
176  // Check the reversed triplet too
177  std::reverse(ringTriplet.begin(), ringTriplet.end());
178 
179  // Check for equality
180  if (triplet == ringTriplet) {
181  return true;
182  } // reversed triplet matches
183  } // first element of ringTriplet
184 
185  return false;
186 }
T clear(T... args)
Topological network criteria functions.
Definition: ring.hpp:64
T reverse(T... args)

◆ getEdgeMoleculesInRings()

void ring::getEdgeMoleculesInRings ( std::vector< std::vector< int >>  rings,
molSys::PointCloud< molSys::Point< double >, double > *  oCloud,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::array< double, 3 >  coordLow,
std::array< double, 3 >  coordHigh,
bool  identicalCloud = false 
)

Select edge molecules and atoms which are part of rings, such that rings formed with even one atom in the slice will be included in the selection Modifies the inSlice bool of a given PointCloud (this may be the same) as the given oxygen atom PointCloud which was used to construct the neighbour list used to construct the rings vector of vectors. We assume that the PointCloud structs have the inSlice bool values set according to the presence of the atom in the slice (this can be done using the gen::moleculesInSingleSlice function.

Function that loops through the PointCloud used to construct the neighbour list (used to generate primitive rings) and sets the inSlice bool values of edge atoms which belong to rings that are formed by atoms in the slice. The output PointCloud may not be the same as the PointCloud used to construct the nList, and the inSlice bool value can be set for both.

Parameters
[in]ringsVector of vectors of the primitive rings (by index) according to oCloud
[in]oCloudPointCloud of O atoms, used to construct the rings vector of vectors
[in]yCloudThe output PointCloud (may contain more than just the O atoms)
[in]identicalCloudbool value; if this is true then oCloud and yCloud are the same
[in]coordLowContains the lower limits of the slice (needed to find all the molecules of oCloud in the slice)
[in]coordHighContains the upper limits of the slice

Definition at line 265 of file selection.cpp.

268  {
269  //
270  // A vector of bool values, such that every ring has a value of true (in the slice) or false (not in the slice)
271  std::vector<bool> ringInSlice(rings.size(), false); // all set to false initially.
272  int jatomIndex, jatomID; // Index and ID in oCloud
273  int jatomIndex1; // Index in yCloud
275  molIDAtomIDmap; // Unordered multimap with molecule IDs of the atoms as the keys and the
276  // atom IDs as the values. More than one atom can have the same molecule ID
277 
278  // --------------------
279  // Sets all O atoms within the slice to an inSlice bool value of true. If a single atom of a molecule is in the
280  // slice, then all atoms in the molecule will also be inside the slice
281  gen::moleculesInSingleSlice(oCloud, true, coordLow, coordHigh);
282  // --------------------
283  // Get the unordered map of the atom IDs (keys) and the molecular IDs
284  // (values)
285  molIDAtomIDmap = molSys::createMolIDAtomIDMultiMap(yCloud);
286  // --------------------
287 
288  // Loop through every iatom present in the slice in oCloud
289  // Check to see if iatom is in any ring. If it is present in a ring, set that ring to true
290  // Change the inSlice bool of every iatom in the ring to true if false. Set also the inSlice bool
291  // of the output cloud if it is not identical.
292 
293  // Loop through every iatom present in the slice in oCloud to determine which rings are in the slice
294  // The indices of oCloud and rings should match (or this will produce unexpected results)
295  for (int iatom = 0; iatom < oCloud->nop; iatom++) {
296  // Skip if iatom is not in the slice
297  if (!oCloud->pts[iatom].inSlice)
298  {
299  continue;
300  } // skip for iatom not in slice
301  //
302  // For iatom in the slice,
303  // loop through all rings to find all the rings it is a part of
304  for (int iring = 0; iring < rings.size(); iring++)
305  {
306  // Skip if iring is in the slice already
307  if (ringInSlice[iring])
308  {
309  continue;
310  } // skip for iring in slice
311  // Check and see if iatom is in iring
312  if(std::find(rings[iring].begin(), rings[iring].end(), iatom)!=rings[iring].end()){
313  // Found iatom; ring is part of the slice
314  ringInSlice[iring] = true; // update the vector of bool values
315  // --------------------------
316  // Change the inSlice bool of every iatom in oCloud
317  // (and optionally, yCloud) in the ring to true
318  // Loop through the elements of the ring
319  for (int j = 0; j < rings[iring].size(); j++)
320  {
321  jatomIndex = rings[iring][j]; // Index of the atom in oCloud
322  // Set this to true in oCloud
323  oCloud->pts[jatomIndex].inSlice = true; // part of slice
324  jatomID = oCloud->pts[jatomIndex].atomID; // Atom ID
325  // Now if oCloud and yCloud are not the same, use the
326  // atom ID to set the inSlice bool value in yCloud
327  if (!identicalCloud)
328  {
329  // Find the index corresponding to the same atom in yCloud
330  auto gotJ = yCloud->idIndexMap.find(jatomID);
331  jatomIndex1 = gotJ->second;
332  // throw if not found ?
333  // Set the jatom inSlice bool to true
334  yCloud->pts[jatomIndex1].inSlice = true; // jatomIndex is inside the slice
335  // set the inSlice value of all atoms in yCloud with the current molecule ID
336  gen::setAtomsWithSameMolID(yCloud, molIDAtomIDmap,
337  yCloud->pts[jatomIndex1].molID, true);
338  } // end of setting values in yCloud
339  } // end of loop through the elements of the current ring
340  // --------------------------
341  } // found iatom in the ring
342  //
343  } // end of loop through all rings searching for iatom
344  //
345  } // end of loop through atoms in oCloud in the slice
346 
347 
348  // Optionally add support for molecule slice update if yCloud and oCloud are not identical?
349  // Some other bool?
350  // IC(ringInSlice);
351 
352  return;
353 }
std::unordered_multimap< int, int > createMolIDAtomIDMultiMap(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: mol_sys.cpp:67
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:176
void setAtomsWithSameMolID(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::unordered_multimap< int, int > molIDAtomIDmap, int molID, bool inSliceValue=true)
Definition: selection.cpp:225
void moleculesInSingleSlice(molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool clearPreviousSliceSelection=true, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
Definition: selection.cpp:147

◆ getPointCloudOneAtomType()

molSys::PointCloud< molSys::Point< double >, double > gen::getPointCloudOneAtomType ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
molSys::PointCloud< molSys::Point< double >, double > *  outCloud,
int  atomTypeI,
bool  isSlice = false,
std::array< double, 3 >  coordLow = std::array<double, 3>{0, 0, 0},
std::array< double, 3 >  coordHigh = std::array<double, 3>{0, 0, 0} 
)

Given a pointCloud containing certain atom types, this returns a pointCloud containing atoms of only the desired type

Function that loops through a given input pointCloud and returns a new pointCloud only containing atoms of a given atom type ID.
This is registered as a Lua function, and is exposed to the user directly.

Parameters
[in]yCloudThe given input PointCloud
[out]outCloudThe output PointCloud
[in]atomTypeIThe type ID of the atoms to save in the output PointCloud
[in]isSliceThis decides whether a slice will be used or not
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 37 of file selection.cpp.

41  {
42  //
43  int natoms = 0; // Number of atoms of the desired type
44  bool pointInSlice = true; // If the current point is in the slice, this is true (otherwise this is false)
45 
46  // --------
47  // Before filling up the PointCloud, if the vectors are filled
48  // empty them
49  *outCloud = molSys::clearPointCloud(outCloud);
50  // --------
51 
52  // Loop through every iatom and check the type
53  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
54  // Skip if the atom is not of the desired type
55  if (yCloud->pts[iatom].type != atomTypeI) {
56  continue;
57  } // end of checking for the type
58  // ----------------
59  // only if a slice has been requested
60  if (isSlice) {
61  pointInSlice = sinp::atomInSlice(yCloud->pts[iatom].x, yCloud->pts[iatom].y, yCloud->pts[iatom].z,
62  coordLow, coordHigh);
63  // Skip if the atom is not part of the slice
64  if (!pointInSlice) {
65  continue;
66  } // skipped for atoms not in the slice
67  } // end of slice handling
68  //
69  // Actually add the atoms to the outCloud
70  natoms++; // Update the number of atoms in outCloud
71  outCloud->pts.push_back(yCloud->pts[iatom]); // Update the pts vector
72  outCloud->idIndexMap[yCloud->pts[iatom].atomID] =
73  outCloud->pts.size() - 1; // array index
74  }
75 
76  // Update the number of particles in the PointCloud
77  outCloud->nop = outCloud->pts.size();
78 
79  // Box and box lengths
80  outCloud->box = yCloud->box;
81  outCloud->boxLow = yCloud->boxLow;
82 
83  // Update the frame number
84  outCloud->currentFrame = yCloud->currentFrame;
85 
86  return *outCloud;
87 }
std::vector< T > boxLow
Periodic box lengths.
Definition: mol_sys.hpp:175
molSys::PointCloud< molSys::Point< double >, double > clearPointCloud(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
//! Function for clearing vectors in PointCloud after multiple usage
Definition: mol_sys.cpp:24

◆ 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 200 of file ring.cpp.

200  {
201  //
202  std::vector<std::vector<int>> ringSingleSize; // rings of one size
203 
204  // rings contains primitive rings of all sizes
205  // Only save rings of a given size (ringSize) to the new
206  // vector of vectors, ringSingleSize
207  for (int iring = 0; iring < rings.size(); iring++) {
208  // Check the size of the current ring
209  // If it is the correct size, save it in ringSingleSize
210  if (rings[iring].size() == ringSize) {
211  ringSingleSize.push_back(rings[iring]);
212  } // End of check of the size of iring
213  } // end of loop through all rings in rings
214 
215  return ringSingleSize;
216 }

◆ 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 226 of file ring.cpp.

226  {
227  std::vector<int> commonElements; // Vector containing common elements
228 
229  // Sort the vectors before finding common elements
230  sort(ring1.begin(), ring1.end());
231  sort(ring2.begin(), ring2.end());
232 
233  // Find intersection of sorted vectors
234  auto it =
235  std::set_intersection(ring1.begin(), ring1.end(), ring2.begin(),
236  ring2.end(), std::back_inserter(commonElements));
237 
238  // If there are no elements in common, then return false
239  if (commonElements.size() == 0) {
240  return false;
241  }
242  // If there are common elements, return true
243  else {
244  return true;
245  }
246 }
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.

◆ moleculesInSingleSlice()

void gen::moleculesInSingleSlice ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
bool  clearPreviousSliceSelection = true,
std::array< double, 3 >  coordLow = std::array<double, 3>{0, 0, 0},
std::array< double, 3 >  coordHigh = std::array<double, 3>{0, 0, 0} 
)

Given a pointCloud set the inSlice bool for every atom, if the molecules are inside the specified (single) region. If even one atom of a molecule is inside the region, then all atoms of that molecule will be inside the region (irrespective of type)

Function that loops through a given input PointCloud and sets the inSlice bool for every Point according to whether the molecule
is in the specified (single) slice or not. If even one atom of a molecule is inside the region, then all atoms belonging to that molecule should be inside the slice as well (therefore, inSlice would be set to true) NOTE: THIS DOES NOT WORK. ERROR

Parameters
[in]yCloudThe given input PointCloud
[in]clearPreviousSliceSelectionsets all inSlice bool values to false before adding Points to the slice
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 147 of file selection.cpp.

151  {
152  //
153  bool pointInSlice = true; // If the current point is in the slice, this is true (otherwise this is false)
155  molIDAtomIDmap; // Unordered multimap with molecule IDs of the atoms as the keys and the
156  // atom IDs as the values. More than one atom can have the same molecule ID
157  int iatomMolID; // molID of the current iatom
158  int jatomID; // atom ID of the current jatom
159  int jatomIndex; // index of jatom
160 
161  // --------------------
162  // Get the unordered map of the atom IDs (keys) and the molecular IDs
163  // (values)
164  molIDAtomIDmap = molSys::createMolIDAtomIDMultiMap(yCloud);
165  // --------------------
166  // Set inSlice to false for every atom first
167  if (clearPreviousSliceSelection)
168  {
169  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
170  yCloud->pts[iatom].inSlice = false;
171  }
172  } // end of init
173  // --------------------
174 
175  // Loop through every iatom and check if the atom is in the slice or not
176  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
177  // Check if the current atom is in the slice
178  pointInSlice = sinp::atomInSlice(yCloud->pts[iatom].x, yCloud->pts[iatom].y, yCloud->pts[iatom].z,
179  coordLow, coordHigh);
180  //
181  // Set the inSlice bool for the particular Point
182  if (pointInSlice)
183  {
184  yCloud->pts[iatom].inSlice = true; // iatom is inside the slice
185  // Find mol ID and atoms with that particular molecular ID
186  iatomMolID = yCloud->pts[iatom].molID; // molecule ID
187  // -----------
188  // Find all atoms with iatomMolID and set the inSlice bool
189  // to true
190  auto range = molIDAtomIDmap.equal_range(iatomMolID);
191  // Loop through all atoms with iatomMolID
192  for (auto it = range.first; it != range.second; it++)
193  {
194  // it->second gives the value (in this case, the atom ID)
195  jatomID = it->second; // Atom ID with molecule ID equal to iatomMolID
196  auto gotJ = yCloud->idIndexMap.find(jatomID);
197  jatomIndex = gotJ->second;
198  // Set the jatom inSlice bool to true
199  yCloud->pts[jatomIndex].inSlice = true; // jatomIndex is inside the slice
200  }
201  // -----------
202  } // the atom is in the slice
203  else{
204  yCloud->pts[iatom].inSlice = false; // iatom is not in the slice
205  } // atom is not in the slice
206  }
207 
208  return;
209 }
T equal_range(T... args)

◆ 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 clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
Definition: ring.cpp:22
int assignPolygonType(std::vector< std::vector< int >> rings, std::vector< int > *atomTypes, std::vector< int > nRings)
Definition: ring.cpp:40
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:200
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, bool isMonolayer=true)
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)

◆ printSliceGetEdgeMoleculesInRings()

void ring::printSliceGetEdgeMoleculesInRings ( std::string  path,
std::vector< std::vector< int >>  rings,
molSys::PointCloud< molSys::Point< double >, double > *  oCloud,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::array< double, 3 >  coordLow,
std::array< double, 3 >  coordHigh,
bool  identicalCloud = false 
)

Master function for selecting edge molecules and atoms which are part of rings, such that rings formed with even one atom in the slice will be included in the selection Modifies the inSlice bool of a given PointCloud (this may be the same) as the given oxygen atom PointCloud which was used to construct the neighbour list used to construct the rings vector of vectors (calls ring::getEdgeMoleculesInRings) Prints out molecule IDs individually of molecules in the slice, and also prints out a LAMMPS data file of just the molecules and atoms in the slice

Function that loops through the PointCloud used to construct the neighbour list (used to generate primitive rings) and sets the inSlice bool values of edge atoms which belong to rings that are formed by atoms in the slice. The selected molecule IDs are printed to a separate file, and the molecules and atoms in the slice will be output to a LAMMPS data file (with the original box volume) The output PointCloud may not be the same as the PointCloud used to construct the nList, and the inSlice bool value can be set for both.

Parameters
[in]ringsVector of vectors of the primitive rings (by index) according to oCloud
[in]oCloudPointCloud of O atoms, used to construct the rings vector of vectors
[in]yCloudThe output PointCloud (may contain more than just the O atoms)
[in]identicalCloudbool value; if this is true then oCloud and yCloud are the same
[in]coordLowContains the lower limits of the slice (needed to find all the molecules of oCloud in the slice)
[in]coordHighContains the upper limits of the slice

Definition at line 375 of file selection.cpp.

379  {
380  //
381 
382  //Given the full yCloud PointCloud, set the inSlice bool for every atom,
383  // if the molecules are inside the specified (single) region.
384  // gen::atomsInSingleSlice(yCloud, true, coordLow, coordHigh);
385  gen::moleculesInSingleSlice(yCloud, true, coordLow, coordHigh);
386 
387  // Make sure that molecules which participate in the rings inside the slice are also
388  // in the selection
389  ring::getEdgeMoleculesInRings(rings, oCloud, yCloud, coordLow, coordHigh, identicalCloud);
390 
391  // Print out the molecule IDs of all the atoms in the slice
392  sout::writeMoleculeIDsInSlice(path, yCloud);
393  // Print out a command that could be used for an expression select command in OVITO
395 
396  // Print out the dump of all atoms and molecules, with an inSlice value printed in a separate column
397  // H atoms not included in the slice (TODO: fix)
398  sout::writeLAMMPSdumpSlice(yCloud, path);
399 
400  IC(rings[0]);
401 
402  return;
403 }
void getEdgeMoleculesInRings(std::vector< std::vector< int >> rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
Definition: selection.cpp:265
int writeMoleculeIDsInSlice(std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
int writeLAMMPSdumpSlice(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path)
int writeMoleculeIDsExpressionSelectOVITO(std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud)

◆ 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 49 of file topo_one_dim.cpp.

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

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

◆ 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 791 of file topo_one_dim.cpp.

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

◆ setAtomsWithSameMolID()

void gen::setAtomsWithSameMolID ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::unordered_multimap< int, int >  molIDAtomIDmap,
int  molID,
bool  inSliceValue = true 
)

Given a particular molecule ID and a pointCloud set the inSlice bool for all atoms, with that molecule ID

Function that loops through a given input PointCloud and sets the inSlice bool for every Point according to whether the molecule
is in the specified (single) slice or not. If even one atom of a molecule is inside the region, then all atoms belonging to that molecule should be inside the slice as well (therefore, inSlice would be set to true)

Parameters
[in]yCloudThe given input PointCloud
[in]clearPreviousSliceSelectionsets all inSlice bool values to false before adding Points to the slice
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 225 of file selection.cpp.

228  {
229  //
230  int jatomID; // atom ID of the current jatom
231  int jatomIndex; // index of jatom
232 
233  // Find all atoms with molID and set the inSlice bool to inSliceValue
234  auto range = molIDAtomIDmap.equal_range(molID);
235 
236  // Loop through all atoms with molID
237  for (auto it = range.first; it != range.second; it++){
238  // it->second gives the value (in this case, the atom ID)
239  jatomID = it->second; // Atom ID with molecule ID equal to iatomMolID
240  auto gotJ = yCloud->idIndexMap.find(jatomID);
241  jatomIndex = gotJ->second;
242  // Set the jatom inSlice bool to true
243  yCloud->pts[jatomIndex].inSlice = inSliceValue; // jatomIndex is assigned inSliceValue
244  } // end of loop through all atoms with molID
245 
246  // IC(yCloud->pts[jatomIndex]);
247 
248  return;
249 }

◆ 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 323 of file bulkTUM.cpp.

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

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

◆ 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 581 of file bulkTUM.cpp.

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

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

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