bulkTUM.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <bulkTUM.hpp>
12 
13 // -----------------------------------------------------------------------------------------------------
14 // TOPOLOGICAL UNIT MATCHING ALGORITHMS
15 // -----------------------------------------------------------------------------------------------------
16 
51  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int firstFrame,
52  bool printClusters, bool onlyTetrahedral) {
53  // The input rings vector has rings of all sizes
54 
55  // ringType has a value for rings of a particular size
57  ringType; // This vector will have a value for each ring inside
58  // ringList, of type enum strucType in gen.hpp
59  // Make a list of all the DDCs and HCs
60  std::vector<cage::Cage> cageList;
62  ringsOneType; // Vector of vectors of rings of a single size
63  int initRingSize; // Todo or not: calculate the PNCs or not
64  int maxRingSize = 6; // DDCs and HCs are for 6-membered rings
66  atomTypes; // This vector will have a value for every atom
67  // Number of types
68  int numHC, numDDC, mixedRings, prismaticRings, basalRings;
69  // Shape-matching variables ----
70  double rmsd; // RMSD value for a particular cage type
71  std::vector<double> quat; // Quaternion obtained from shape-matching
72  std::vector<std::vector<double>> quatList; // List of quaternions
73  // Reference points
74  Eigen::MatrixXd refPntsDDC(14,
75  3); // Reference point set (Eigen matrix) for a DDC
76  Eigen::MatrixXd refPntsHC(12,
77  3); // Reference point set (Eigen matrix) for a DDC
78 
79  // Vector for the RMSD per atom and the RMSD per ring:
80  std::vector<double> rmsdPerAtom;
82  noOfCommonElements; // An atom can belong to more than one ring or shape
83  //
84 
85  if (onlyTetrahedral) {
86  initRingSize = 6;
87  } else {
88  initRingSize = 5;
89  }
90 
91  // Init the atom type vector
92  atomTypes.resize(yCloud->nop); // Has a value for each atom
93  // Init the rmsd per atom
94  rmsdPerAtom.resize(yCloud->nop, -1); // Has a value for each atom
95  noOfCommonElements.resize(yCloud->nop); // Has a value for each atom
96 
97  // ----------------------------------------------
98  // Init
99  // Get rings of size 5 or 6.
100  for (int ringSize = initRingSize; ringSize <= maxRingSize; ringSize++) {
101  // Clear ringsOneType
102  ring::clearRingList(ringsOneType);
103  // Get rings of the current ring size
104  ringsOneType = ring::getSingleRingSize(rings, ringSize);
105  // Skip for zero rings
106  if (ringsOneType.size() == 0) {
107  continue;
108  } // skip for no rings of ringSize
109  //
110  // Init the ringType vector
111  ringType.resize(
112  ringsOneType.size()); // Has a value for each ring. init to zero.
113  // ----------------------------------------------
114  if (ringSize == 6) {
115  // Get the cages
116 
117  // Find the DDCs and HCs
118  cageList = tum3::topoBulkCriteria(path, ringsOneType, nList, yCloud,
119  firstFrame, &numHC, &numDDC, &ringType);
120 
121  // Gets the atom type for every atom, to be used for printing out the ice
122  // types found
123  ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
124  }
125  // Finding the 5-membered prismatic blocks
126  else {
127  // Get the prism block classifications
128  prism3::findBulkPrisms(ringsOneType, &ringType, nList, yCloud,
129  &rmsdPerAtom);
130  // Gets the atom type for every atom, to be used for printing out the ice
131  // types found
132  ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
133  }
134  }
135 
136  // --------------------------------------------------
137  // SHAPE-MATCHING
138  //
139  // Get the reference point sets
140  //
141  std::string filePathHC = "templates/hc.xyz";
142  std::string filePathDDC = "templates/ddc.xyz";
143  refPntsHC = tum3::buildRefHC(filePathHC); // HC
144  refPntsDDC = tum3::buildRefDDC(filePathDDC); // DDC
145  //
146  // Init
147  //
148  // Loop through all the atoms and set commonElements to 1 for PNC atoms
149  for (int i = 0; i < yCloud->nop; i++) {
150  if (rmsdPerAtom[i] != -1) {
151  noOfCommonElements[i] = 1;
152  } // for atoms which have been updated
153  } // end of updating commonElements
154  //
155  // Loop through the entire cageList vector of cages to match the HCs and DDCs
156  //
157  // HCs are first, followed by DDCs.
158  //
159  // Go through all the HCs
160  for (int icage = 0; icage < numHC; icage++) {
161  // Match against a perfect HC
162  tum3::shapeMatchHC(yCloud, refPntsHC, cageList[icage], ringsOneType, nList,
163  &quat, &rmsd);
164  // Update the vector of quaternions
165  quatList.push_back(quat);
166  // Update the RMSD per ring
167  tum3::updateRMSDatom(ringsOneType, cageList[icage], rmsd, &rmsdPerAtom,
168  &noOfCommonElements, atomTypes);
169  } // end of looping through all HCs
170  // --------------------------------------------------
171  // Go through all the DDCs
172  for (int icage = numHC; icage < cageList.size(); icage++) {
173  // Match against a perfect DDC
174  tum3::shapeMatchDDC(yCloud, refPntsDDC, cageList, icage, ringsOneType,
175  &quat, &rmsd);
176  // Update the vector of quaternions
177  quatList.push_back(quat);
178  // Update the RMSD per ring
179  tum3::updateRMSDatom(ringsOneType, cageList[icage], rmsd, &rmsdPerAtom,
180  &noOfCommonElements, atomTypes);
181  } // end of looping through all HCs
182 
183  // --------------------------------------------------
184  // Getting the RMSD per atom
185  // Average the RMSD per atom
186  tum3::averageRMSDatom(&rmsdPerAtom, &noOfCommonElements);
187  // --------------------------------------------------
188 
189  // Print out the lammps data file with the bonds and types
190  sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path);
191  // To output the bonds between dummy atoms, uncomment the following line
192  // sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path, true);
193 
194  // --------------------------------------------------
195  // Writes out the dumpfile
196  //
197  std::vector<int> dumpAtomTypes; // Must be typecast to int
198  dumpAtomTypes.resize(atomTypes.size());
199  // Change from enum to int
200  for (int i = 0; i < atomTypes.size(); i++) {
201  if (atomTypes[i] == cage::hc) {
202  dumpAtomTypes[i] = 1;
203  } // HC
204  else if (atomTypes[i] == cage::ddc) {
205  dumpAtomTypes[i] = 2;
206  } // DDC
207  else if (atomTypes[i] == cage::mixed) {
208  dumpAtomTypes[i] = 3;
209  } // mixed rings
210  else if (atomTypes[i] == cage::pnc) {
211  dumpAtomTypes[i] = 4;
212  } // pnc
213  else if (atomTypes[i] == cage::mixed2) {
214  dumpAtomTypes[i] = 5;
215  } // shared by pnc and ddc/hc
216  else {
217  dumpAtomTypes[i] = 0;
218  } // dummy atoms
219  } // end of changing the type to int
220  sout::writeLAMMPSdumpCages(yCloud, rmsdPerAtom, dumpAtomTypes, path,
221  firstFrame);
222  // --------------------------------------------------
223  //
224  // PRINT THE CLUSTERS
225  //
226  if (printClusters) {
227  // Print the clusters
228  tum3::clusterCages(yCloud, path, ringsOneType, cageList, numHC, numDDC);
229  } // end of printing individual clusters
230  // --------------------------------------------------
231  return 0;
232 }
233 
234 // Shape-matching for an HC
248  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
249  const Eigen::MatrixXd &refPoints, cage::Cage cageUnit,
251  std::vector<double> *quat, double *rmsd) {
252  //
253  int iring,
254  jring; // Indices in the ring vector of vectors for the basal rings
255  std::vector<int> basal1,
256  basal2; // Re-ordered basal rings matched to each other
257  int ringSize = 6; // Each ring has 6 nodes
258  Eigen::MatrixXd targetPointSet(12, 3); // Target point set (Eigen matrix)
259  //
260  std::vector<double> rmsdList; // List of RMSD per atom
261  // Variables for looping through possible permutations
262  //
263  std::vector<double> currentQuat; // quaternion rotation
264  double currentRmsd; // least RMSD
265  double currentScale; // scale
266  double scale; // Final scale
267  int index; // Int for describing which permutation matches the best
268 
269  // Init
270  iring = cageUnit.rings[0]; // Index of basal1
271  jring = cageUnit.rings[1]; // Index of basal2
272  rmsdList.resize(yCloud->nop); // Not actually updated here
273  //
274  // ----------------
275  // Re-order the basal rings so that they are matched
276  pntToPnt::relOrderHC(yCloud, rings[iring], rings[jring], nList, &basal1,
277  &basal2);
278  // ----------------
279  // Loop through all possible permutations
280  //
281  for (int i = 0; i < ringSize; i++) {
282  // Change the order of the target points somehow!
283  //
284  targetPointSet = pntToPnt::changeHexCageOrder(yCloud, basal1, basal2, i);
285  // Shape-matching
286  absor::hornAbsOrientation(refPoints, targetPointSet, &currentQuat,
287  &currentRmsd, &rmsdList, &currentScale);
288  if (i == 0) {
289  *quat = currentQuat;
290  *rmsd = currentRmsd;
291  scale = currentScale;
292  index = 0;
293  } else {
294  if (currentRmsd < *rmsd) {
295  *quat = currentQuat;
296  *rmsd = currentRmsd;
297  scale = currentScale;
298  index = i;
299  } // update
300  } // Update if this is a better match
301  } // Loop through possible permutations
302  // ----------------
303  return 0;
304 }
305 
306 // Shape-matching for a DDC
320  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
321  const Eigen::MatrixXd &refPoints, std::vector<cage::Cage> cageList,
322  int cageIndex, std::vector<std::vector<int>> rings,
323  std::vector<double> *quat, double *rmsd) {
324  //
325  std::vector<int> ddcOrder; // Connectivity of the DDC
326  int ringSize = 6; // Each ring has 6 nodes
327  Eigen::MatrixXd targetPointSet(14, 3); // Target point set (Eigen matrix)
328  //
329  std::vector<double> rmsdList; // List of RMSD per atom
330  // Variables for looping through possible permutations
331  //
332  std::vector<double> currentQuat; // quaternion rotation
333  double currentRmsd; // least RMSD
334  double currentScale; // scale
335  double scale; // Final scale
336  int index; // Int for describing which permutation matches the best
337 
338  // ----------------
339  // Save the order of the DDC in a vector
340  ddcOrder = pntToPnt::relOrderDDC(cageIndex, rings, cageList);
341  // ----------------
342  // Loop through all possible permutations
343  //
344  for (int i = 0; i < ringSize; i++) {
345  // Change the order of the target points somehow!
346  //
347  targetPointSet = pntToPnt::changeDiaCageOrder(yCloud, ddcOrder, i);
348  // Shape-matching
349  absor::hornAbsOrientation(refPoints, targetPointSet, &currentQuat,
350  &currentRmsd, &rmsdList, &currentScale);
351  if (i == 0) {
352  *quat = currentQuat;
353  *rmsd = currentRmsd;
354  scale = currentScale;
355  index = 0;
356  } else {
357  if (currentRmsd < *rmsd) {
358  *quat = currentQuat;
359  *rmsd = currentRmsd;
360  scale = currentScale;
361  index = i;
362  } // update
363  } // Update if this is a better match
364  } // Loop through possible permutations
365  // ----------------
366  return 0;
367 }
368 
373 Eigen::MatrixXd tum3::buildRefHC(std::string fileName) {
374  //
375  Eigen::MatrixXd refPnts(12, 3); // Reference point set (Eigen matrix)
376  // Get the reference HC point set
378  setCloud; // PointCloud for holding the reference point values
379  // Variables for rings
380  std::vector<std::vector<int>> nList; // Neighbour list
381  std::vector<std::vector<int>> rings; // Rings
383  ringType; // This vector will have a value for each ring inside
384  std::vector<int> listHC; // Contains atom indices of atoms making up HCs
385  // Make a list of all the DDCs and HCs
386  std::vector<cage::Cage> cageList;
387  int iring, jring;
388  //
389  // read in the XYZ file into the pointCloud setCloud
390  //
391  sinp::readXYZ(fileName, &setCloud);
392  // box lengths
393  for (int i = 0; i < 3; i++) {
394  setCloud.box[i] = 50;
395  } // end of setting box lengths
396  //
397 
398  nList = nneigh::neighListO(3.5, &setCloud, 1);
399  // Neighbour list by index
400  nList = nneigh::neighbourListByIndex(&setCloud, nList);
401  // Find the vector of vector of rings
402  rings = primitive::ringNetwork(nList, 6);
403  // init the ringType vector
404  ringType.resize(rings.size());
405  // Find the HCs
406  listHC = ring::findHC(rings, &ringType, nList, &cageList);
407  // Get the basal rings from cageList
408  iring = cageList[0].rings[0];
409  jring = cageList[0].rings[1];
410  //
411  std::vector<int> matchedBasal1,
412  matchedBasal2; // Re-ordered basal rings 1 and 2
413  // Reordered basal rings
414  // Getting the target Eigen vectors
415  // Get the re-ordered matched basal rings, ordered with respect to each
416  // other
417 
418  pntToPnt::relOrderHC(&setCloud, rings[iring], rings[jring], nList,
419  &matchedBasal1, &matchedBasal2);
420  // Get the reference point set
421  refPnts =
422  pntToPnt::changeHexCageOrder(&setCloud, matchedBasal1, matchedBasal2, 0);
423  //
424  return refPnts;
425 }
426 
431 Eigen::MatrixXd tum3::buildRefDDC(std::string fileName) {
432  //
433  Eigen::MatrixXd refPnts(14, 3); // Reference point set (Eigen matrix)
434  // Get the reference HC point set
436  setCloud; // PointCloud for holding the reference point values
437  // Variables for rings
438  std::vector<std::vector<int>> nList; // Neighbour list
439  std::vector<std::vector<int>> rings; // Rings
441  ringType; // This vector will have a value for each ring inside
442  std::vector<int> listDDC,
443  listHC; // Contains atom indices of atoms making up DDCs and HCs
444  std::vector<int> ddcOrder; // Atom indices of particles in the DDC
445  // Make a list of all the DDCs and HCs
446  std::vector<cage::Cage> cageList;
447  int iring, jring;
448  //
449  // read in the XYZ file into the pointCloud setCloud
450  //
451  sinp::readXYZ(fileName, &setCloud);
452 
453  // box lengths
454  for (int i = 0; i < 3; i++) {
455  setCloud.box[i] = 50;
456  } // end of setting box lengths
457  //
458 
459  nList = nneigh::neighListO(3.5, &setCloud, 1);
460  // Neighbour list by index
461  nList = nneigh::neighbourListByIndex(&setCloud, nList);
462  // Find the vector of vector of rings
463  rings = primitive::ringNetwork(nList, 6);
464  // init the ringType vector
465  ringType.resize(rings.size());
466  // Find the DDCs
467  listDDC = ring::findDDC(rings, &ringType, listHC, &cageList);
468  // Save the order of the DDC in a vector
469  ddcOrder = pntToPnt::relOrderDDC(0, rings, cageList);
470  // Get the reference point set
471  refPnts = pntToPnt::changeDiaCageOrder(&setCloud, ddcOrder, 0);
472  //
473  return refPnts;
474 }
475 
483  cage::Cage cageUnit, double rmsd,
484  std::vector<double> *rmsdPerAtom,
485  std::vector<int> *noOfCommonAtoms,
486  std::vector<cage::iceType> atomTypes) {
487  //
488  int nRings = cageUnit.rings.size(); // Number of rings in the current cage
489  int iring; // Index according to the rings vector of vector, for the current
490  // ring inside the cage
491  int iatom; // Current atom index
492  int ringSize = rings[0].size(); // Number of nodes in each ring (6)
493 
494  // Loop through the rings in each cage
495  for (int i = 0; i < nRings; i++) {
496  iring = cageUnit.rings[i]; // Current ring
497 
498  // Loop through all the atoms in the ring
499  for (int j = 0; j < ringSize; j++) {
500  iatom = rings[iring][j]; // Current atom index
501 
502  // Skip for PNC atoms
503  if (atomTypes[iatom] == cage::pnc || atomTypes[iatom] == cage::mixed2) {
504  continue;
505  } // Do not update if the atom is a PNC
506  //
507  // UPDATE
508  if ((*rmsdPerAtom)[iatom] == -1) {
509  (*rmsdPerAtom)[iatom] = rmsd;
510  (*noOfCommonAtoms)[iatom] = 1;
511  } else {
512  (*rmsdPerAtom)[iatom] += rmsd;
513  (*noOfCommonAtoms)[iatom] += 1;
514  }
515  } // end of looping through all the atoms in the ring
516 
517  } // end of looping through the rings in the cage
518 
519  return 0;
520 }
521 
526  std::vector<int> *noOfCommonAtoms) {
527  //
528  int nop = (*rmsdPerAtom).size(); // Number of particles
529 
530  for (int iatom = 0; iatom < nop; iatom++) {
531  //
532  if ((*noOfCommonAtoms)[iatom] == 0) {
533  (*rmsdPerAtom)[iatom] = -1; // Dummy atom
534  } else {
535  (*rmsdPerAtom)[iatom] /= (*noOfCommonAtoms)[iatom];
536  }
537  } // end of averaging
538 
539  return 0;
540 } // end of function
541 
542 // -----------------------------------------------------------------
543 //
544 // TOPOLOGICAL NETWORK ALGORITHMS
545 //
546 // -----------------------------------------------------------------
580  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int firstFrame,
581  int *numHC, int *numDDC, std::vector<ring::strucType> *ringType) {
582  //
583  // Ring IDs of each type will be saved in these vectors
584  std::vector<int> listDDC; // Vector for ring indices of DDC
585  std::vector<int> listHC; // Vector for ring indices of HC
587  listMixed; // Vector for ring indices of rings that are both DDC and HC
588  // ringList, of type enum strucType in gen.hpp
589  // Make a list of all the DDCs and HCs
590  std::vector<cage::Cage> cageList;
591  // Number of types
592  int mixedRings, prismaticRings, basalRings;
593 
594  // ----------------------------------------------
595  // Init
596  //
597  *numHC = 0; // Number of hexagonal cages
598  *numDDC = 0; // Init the number of DDCs
599  // Quit the function for zero rings
600  if (rings.size() == 0) {
601  return cageList;
602  } // skip for no rings of ringSize
603  //
604  // Init the ringType vector
605  (*ringType).resize(rings.size()); // Has a value for each ring. init to zero.
606  // ----------------------------------------------
607  // Get the cages
608 
609  // Find HC rings, saving the ring IDs (starting from 0) to listHC
610  listHC = ring::findHC(rings, ringType, nList, &cageList);
611 
612  // Find DDC rings, saving the IDs to listDDC
613  listDDC = ring::findDDC(rings, ringType, listHC, &cageList);
614 
615  // Find rings which are both DDCs and HCs (mixed)
616  // A dummy value of -10 in the listDDC and listHC vectors for mixed rings
617  listMixed = ring::findMixedRings(rings, ringType, &listDDC, &listHC);
618 
619  // Get the number of structures (DDCs, HCs, mixed rings, basal rings,
620  // prismatic rings)
621  ring::getStrucNumbers(*ringType, cageList, numHC, numDDC, &mixedRings,
622  &prismaticRings, &basalRings);
623 
624  // Write out to a file
625  sout::writeTopoBulkData(path, yCloud->currentFrame, *numHC, *numDDC,
626  mixedRings, basalRings, prismaticRings, firstFrame);
627 
628  return cageList;
629 }
630 
631 // -----------------------------------------------------------------
632 //
633 // CLUSTERING ALGORITHMS
634 //
635 // -----------------------------------------------------------------
643  int numHC, int numDDC) {
644  //
645  std::vector<int> linkedList; // Linked list for the clustered cages
646  int nCages = numHC + numDDC; // Number of cages in total
647  int j, temp; // Variables used inside the loops
648  bool isNeighbour; // true if cages are adjacent to each other, false otherwise
650  visited; // To make sure you don't go through the same cages again.
651  int singleDDCs, singleHCs; // Number of single DDCs and HCs
652  // Units are in Angstrom^3
653  double volDDC = 72.7; // The calculated alpha volume of a single DDC
654  double volHC = 40.63; // The calculated alpha volume of a single HC
655  int nextElement; // value in list
656  int index; // starting index value
657  int currentIndex; // Current cage index
658  int iClusterNumber; // Number in the current cluster
659  std::vector<int> nClusters; // Number of cages in each cluster
660  std::vector<int> clusterCages; // Indices of each cage in the cluster
661  std::vector<int> atoms; // Vector containing the atom indices of all the
662  // atoms in the cages
663  cage::cageType type; // Type of the cages in the cluster
664  // -----------------------------------------------------------
665  // INITIALIZATION
666  linkedList.resize(nCages); // init to dummy value
667  // Initial values of the list.
668  for (int icage = 0; icage < nCages; icage++) {
669  // Assign the index as the ID
670  linkedList[icage] = icage; // Index
671  } // init of cluster IDs
672  // -----------------------------------------------------------
673  // GETTING THE LINKED LIST
674  //
675  // HCs first
676  for (int i = 0; i < nCages - 1; i++) {
677  // If iatom is already in a cluster, skip it
678  if (linkedList[i] != i) {
679  continue;
680  } // Already part of a cluster
681  //
682  j = i; // Init of j
683  // Execute the next part of the loop while j is not equal to i
684  do {
685  //
686  // Go through the rest of the atoms (KLOOP)
687  for (int k = i + 1; k < nCages; k++) {
688  // Skip if already part of a cluster
689  if (linkedList[k] != k) {
690  continue;
691  } // Already part of a cluster
692  //
693  // k and j must be of the same type
694  if (cageList[k].type != cageList[j].type) {
695  continue;
696  } // skip if k and j are not of the same type
697  // Check to see if k is a nearest neighbour of j
698  isNeighbour =
699  ring::hasCommonElements(cageList[k].rings, cageList[j].rings);
700  if (isNeighbour) {
701  // Swap!
702  temp = linkedList[j];
703  linkedList[j] = linkedList[k];
704  linkedList[k] = temp;
705  } // j and k are neighbouring cages
706  } // end of loop through k (KLOOP)
707  //
708  j = linkedList[j];
709  } while (j != i); // end of control for j!=i
710  } // end of getting the linked list (looping through every i)
711  // -----------------------------------------------------------
712  // WRITE-OUTS
713  //
714  // init
715  visited.resize(nCages);
716  singleDDCs = 0;
717  singleHCs = 0;
718 
719  for (int i = 0; i < nCages; i++) {
720  //
721  if (visited[i]) {
722  continue;
723  } // Already counted
724  // Now that the cage has been visited, set it to true
725  visited[i] = true; // Visited
726  // If only one cage is in the cluster
727  if (linkedList[i] == i) {
728  // Add to the number of single DDCs
729  if (cageList[i].type == cage::DoubleDiaC) {
730  singleDDCs++;
731  } // add to the single DDCs
732  else {
733  singleHCs++;
734  } // single HCs
735  continue;
736  } // cluster has only one cage
737  //
738  // init
739  clusterCages.resize(0);
740  //
741  type = cageList[i].type;
742  currentIndex = i;
743  nextElement = linkedList[currentIndex];
744  index = i;
745  iClusterNumber = 1; // at least one value
746  clusterCages.push_back(i); // Add the first cage
747  // Get the other cages in the cluster
748  while (nextElement != index) {
749  iClusterNumber++;
750  currentIndex = nextElement;
751  visited[currentIndex] = true;
752  clusterCages.push_back(currentIndex); // Add the first cage
753  nextElement = linkedList[currentIndex];
754  } // get number
755  // Update the number of molecules in the cluster
756  nClusters.push_back(iClusterNumber);
757  // --------------
758  // PRINT OUT
759  atoms = tum3::atomsFromCages(rings, cageList, clusterCages);
760  int clusterID = nClusters.size();
761  // Write out the cluster
762  sout::writeXYZcluster(path, yCloud, atoms, clusterID, type);
763  // Print out the cages into an XYZ file
764  // --------------
765  } // end of loop through
766  // -----------------------------------------------------------
767  // Write out the stuff for single cages
768  // ----------------
769  std::ofstream outputFile;
770  // Make the output directory if it doesn't exist
771  sout::makePath(path);
772  std::string outputDirName = path + "bulkTopo";
773  sout::makePath(outputDirName);
774  outputDirName = path + "bulkTopo/clusterXYZ/";
775  sout::makePath(outputDirName);
776  outputDirName = path + "bulkTopo/clusterXYZ/frame-" +
777  std::to_string(yCloud->currentFrame);
778  sout::makePath(outputDirName);
779  // ----------------
780  // Write output to file inside the output directory
781  outputFile.open(path + "bulkTopo/clusterXYZ/frame-" +
782  std::to_string(yCloud->currentFrame) + "/info.dat");
783  outputFile << "# volDDC volHC in Angstrom^3\n";
784  outputFile << singleDDCs * volDDC << " " << singleHCs * volHC << "\n";
785  outputFile << "There are " << singleDDCs << " single DDCs\n";
786  outputFile << "There are " << singleHCs << " single HCs\n";
787  outputFile.close(); // Close the file
788  // -----------------------------------------------------------
789  return 0;
790 } // end of the function
791 
792 
798  std::vector<cage::Cage> cageList,
800  //
801  std::vector<int> atoms; // Contains the atom indices (not IDs) of atoms
802  int ringSize = rings[0].size(); // Number of nodes in each ring
803  int iring; // Index of ring in rings vector of vectors
804  int icage; // Index of the current cage in cageList
805  int iatom; // Atom index
806  //
807  for (int i = 0; i < clusterCages.size(); i++) {
808  //
809  icage = clusterCages[i];
810  // Loop through every ring in the cage
811  for (int j = 0; j < cageList[icage].rings.size(); j++) {
812  iring = cageList[icage].rings[j]; // Current ring index
813  // Loop through every atom in iring
814  for (int k = 0; k < ringSize; k++) {
815  iatom = rings[iring][k];
816  atoms.push_back(iatom); // Add the atom
817  } // Loop through every atom in iring
818  } // loop through every ring in the current cage
819  } // end of loop through all cages
820 
821  // Remove the duplicates in the atoms vector
822  // Duplicate IDs must be removed
823  //
824  sort(atoms.begin(), atoms.end());
825  auto ip = std::unique(atoms.begin(), atoms.end());
826  // Resize peripheral rings to remove undefined terms
827  atoms.resize(std::distance(atoms.begin(), ip));
828  //
829  // Return the vector of atoms
830  return atoms;
831 }
T begin(T... args)
File containing functions used specific to bulk topological unit matching (TUM) criterion.
T close(T... args)
T distance(T... args)
T end(T... args)
cageType
Definition: cage.hpp:50
std::vector< int > rings
type of the cage : can be DDC or HC
Definition: cage.hpp:84
@ DoubleDiaC
Definition: cage.hpp:50
@ pnc
Definition: cage.hpp:71
@ ddc
Definition: cage.hpp:71
@ mixed
Definition: cage.hpp:71
@ hc
Definition: cage.hpp:71
@ mixed2
Definition: cage.hpp:71
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:170
std::vector< std::vector< int > > neighListO(double rcutoff, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI)
Definition: neighbours.cpp:114
std::vector< std::vector< int > > neighbourListByIndex(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList)
Definition: neighbours.cpp:337
std::vector< std::vector< int > > ringNetwork(std::vector< std::vector< int >> nList, int maxDepth)
Definition: franzblau.cpp:28
std::vector< int > findDDC(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > listHC, std::vector< cage::Cage > *cageList)
Definition: topo_bulk.cpp:199
int getAtomTypesTopoBulk(std::vector< std::vector< int >> rings, std::vector< ring::strucType > ringType, std::vector< cage::iceType > *atomTypes)
Definition: topo_bulk.cpp:1078
int getStrucNumbers(std::vector< ring::strucType > ringType, std::vector< cage::Cage > cageList, int *numHC, int *numDDC, int *mixedRings, int *prismaticRings, int *basalRings)
Determines the number of HCs, DDCs, Mixed rings, prismatic and basal rings.
Definition: topo_bulk.cpp:1161
int findBulkPrisms(std::vector< std::vector< int >> rings, std::vector< ring::strucType > *ringType, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, double heightCutoff=8)
Find out which rings are prisms.
Definition: topo_bulk.cpp:1242
std::vector< int > findMixedRings(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > *listDDC, std::vector< int > *listHC)
Definition: topo_bulk.cpp:1031
std::vector< int > findHC(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< std::vector< int >> nList, std::vector< cage::Cage > *cageList)
Definition: topo_bulk.cpp:565
int shapeMatchHC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, cage::Cage cageUnit, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, std::vector< double > *quat, double *rmsd)
Shape-matching for a target HC.
Definition: bulkTUM.cpp:247
Eigen::MatrixXd buildRefHC(std::string fileName)
Build a reference Hexagonal cage, reading in from a template XYZ file.
Definition: bulkTUM.cpp:373
std::vector< int > atomsFromCages(std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList, std::vector< int > clusterCages)
Gets the atoms in the cages of a given cluster.
Definition: bulkTUM.cpp:797
int clusterCages(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::vector< std::vector< int >> rings, std::vector< cage::Cage > cageList, int numHC, int numDDC)
Definition: bulkTUM.cpp:640
int averageRMSDatom(std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms)
Average the RMSD per atom.
Definition: bulkTUM.cpp:525
int shapeMatchDDC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, std::vector< cage::Cage > cageList, int cageIndex, std::vector< std::vector< int >> rings, std::vector< double > *quat, double *rmsd)
Shape-matching for a target DDC.
Definition: bulkTUM.cpp:319
int 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)
Definition: bulkTUM.cpp:48
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
Definition: ring.cpp:18
int updateRMSDatom(std::vector< std::vector< int >> rings, cage::Cage cageUnit, double rmsd, std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms, std::vector< cage::iceType > atomTypes)
Definition: bulkTUM.cpp:482
std::vector< cage::Cage > topoBulkCriteria(std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, int *numHC, int *numDDC, std::vector< ring::strucType > *ringType)
Definition: bulkTUM.cpp:577
Eigen::MatrixXd buildRefDDC(std::string fileName)
Build a reference Double-Diamond cage, reading in from a template XYZ file.
Definition: bulkTUM.cpp:431
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Definition: ring.cpp:175
std::vector< std::vector< int > > getSingleRingSize(std::vector< std::vector< int >> rings, int ringSize)
Returns a vector of vectors of rings of a single size.
Definition: ring.cpp:149
int readXYZ(std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Function for reading in atom coordinates from an XYZ file.
Definition: seams_input.cpp:69
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)
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.
int relOrderHC(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, std::vector< std::vector< int >> nList, std::vector< int > *matchedBasal1, std::vector< int > *matchedBasal2)
Matches the order of the basal rings of an HC or a potential HC.
Eigen::MatrixXd changeHexCageOrder(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex=0)
int makePath(const std::string &path)
int writeTopoBulkData(std::string path, int currentFrame, int numHC, int numDDC, int mixedRings, int basalRings, int prismaticRings, int firstFrame)
int writeLAMMPSdataTopoBulk(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< cage::iceType > atomTypes, std::string path, bool bondsBetweenDummy=false)
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.
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 push_back(T... args)
T resize(T... args)
T size(T... args)
This contains a cage, with the constituent rings.
Definition: cage.hpp:82
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166
This contains per-particle information.
Definition: mol_sys.hpp:145
T to_string(T... args)
T unique(T... args)