topo_one_dim.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <topo_one_dim.hpp>
12 
13 // -----------------------------------------------------------------------------------------------------
14 // PRISM ALGORITHMS
15 // -----------------------------------------------------------------------------------------------------
16 
48  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int maxDepth,
49  int *atomID, int firstFrame, int currentFrame, bool doShapeMatching) {
50  //
52  ringsOneType; // Vector of vectors of rings of a single size
53  std::vector<int> listPrism; // Vector for ring indices of n-sided prism
55  ringType; // This vector will have a value for each ring inside
57  atomState; // This vector will have a value for each atom, depending on
58  // the ring type
59  int nPerfectPrisms; // Number of perfect prisms of each type
60  int nImperfectPrisms; // Number of deformed prisms of each type
61  std::vector<int> nPrismList; // Vector of the values of the number of perfect
62  // prisms for a particular frame
63  std::vector<int> nDefPrismList; // Vector of the values of the number of
64  // deformed prisms for a particular frame
66  heightPercent; // Height percent for a particular n and frame
68  atomTypes; // contains int values for each prism type considered
69  double avgPrismHeight = 2.845; // A value of 2.7-2.85 Angstrom is reasonable
70  // Qualifier for the RMSD per atom:
71  std::vector<double> rmsdPerAtom;
72  // -------------------------------------------------------------------------------
73  // Init
74  nPrismList.resize(
75  maxDepth -
76  2); // Has a value for every value of ringSize from 3, upto maxDepth
77  nDefPrismList.resize(maxDepth - 2);
78  heightPercent.resize(maxDepth - 2);
79  // The atomTypes vector is the same size as the pointCloud atoms
80  atomTypes.resize(yCloud->nop, 1); // The dummy or unclassified value is 1
81  // The rmsdPerAtom vector is the same size as the pointCloud atoms and has an
82  // RMSD value for every atom
83  rmsdPerAtom.resize(yCloud->nop, -1);
84  // Resize the atom state vector
85  atomState.resize(yCloud->nop); // Dummy or unclassified
86  // -------------------------------------------------------------------------------
87  // Run this loop for rings of sizes upto maxDepth
88  // The smallest possible ring is of size 3
89  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
90  // Clear ringsOneType
91  ring::clearRingList(ringsOneType);
92  // Get rings of the current ring size
93  ringsOneType = ring::getSingleRingSize(rings, ringSize);
94  //
95  // Continue if there are zero rings of ringSize
96  if (ringsOneType.size() == 0) {
97  nPrismList[ringSize - 3] = 0; // Update the number of prisms
98  nDefPrismList[ringSize - 3] = 0; // Update the number of deformed prisms
99  heightPercent[ringSize - 3] = 0.0; // Update the height%
100  continue;
101  } // skip if there are no rings
102  //
103  // -------------
104  // Init of variables specific to ringSize prisms
105  listPrism.resize(0);
106  ringType.resize(0);
107  nPerfectPrisms = 0;
108  nImperfectPrisms = 0;
109  ringType.resize(
110  ringsOneType.size()); // Has a value for each ring. init to zero.
111  // -------------
112  // Now that you have rings of a certain size:
113  // Find prisms, saving the ring indices to listPrism
114  listPrism = ring::findPrisms(ringsOneType, &ringType, &nPerfectPrisms,
115  &nImperfectPrisms, nList, yCloud, &rmsdPerAtom,
116  doShapeMatching);
117  // -------------
118  nPrismList[ringSize - 3] =
119  nPerfectPrisms; // Update the number of perfect prisms
120  nDefPrismList[ringSize - 3] =
121  nImperfectPrisms; // Update the number of defromed prisms
122  int totalPrisms = nPerfectPrisms + nImperfectPrisms;
123  // Update the height% for the phase
124  heightPercent[ringSize - 3] =
125  topoparam::normHeightPercent(yCloud, totalPrisms, avgPrismHeight);
126  // Continue if there are no prism units
127  if (nPerfectPrisms + nImperfectPrisms == 0) {
128  continue;
129  } // skip for no prisms
130  // Do a bunch of write-outs and calculations
131  // TODO: Write out each individual prism as data files (maybe with an
132  // option)
133  // Get the atom types for a particular prism type
134  ring::assignPrismType(ringsOneType, listPrism, ringSize, ringType,
135  &atomTypes, &atomState);
136  // -------------
137  } // end of loop through every possible ringSize
138 
139  // Calculate the height%
140 
141  // Write out the prism information
142  sout::writePrismNum(path, nPrismList, nDefPrismList, heightPercent, maxDepth,
143  yCloud->currentFrame, firstFrame);
144  // Reassign the prism block types for the deformed prisms
145  if (doShapeMatching) {
146  ring::deformedPrismTypes(atomState, &atomTypes, maxDepth);
147  } // reassign prism block types for deformed prisms
148 
149  // Get rid of translations along the axial dimension
150  ring::rmAxialTranslations(yCloud, atomID, firstFrame, currentFrame);
151 
152  // Write out dump files with the RMSD per atom for the shape-matching
153  // criterion
154  if (doShapeMatching) {
155  sout::writeLAMMPSdumpINT(yCloud, rmsdPerAtom, atomTypes, maxDepth, path);
156  } // reassign prism block types for deformed prisms
157 
158  // Write out the lammps data file for the particular frame
159  sout::writeLAMMPSdataAllPrisms(yCloud, nList, atomTypes, maxDepth, path,
160  doShapeMatching);
161 
162  return 0;
163 }
164 
190  std::vector<ring::strucType> *ringType, int *nPerfectPrisms,
191  int *nImperfectPrisms, std::vector<std::vector<int>> nList,
192  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
193  std::vector<double> *rmsdPerAtom, bool doShapeMatching) {
194  std::vector<int> listPrism;
195  int totalRingNum = rings.size(); // Total number of rings
196  std::vector<int> basal1; // First basal ring
197  std::vector<int> basal2; // Second basal ring
198  bool cond1, cond2; // Conditions for rings to be basal (true) or not (false)
199  bool relaxedCond; // Condition so that at least one bond exists between the
200  // two basal rings
201  bool isAxialPair; // Basal rings should be parallel in one dimension to
202  // prevent overcounting
203  int ringSize = rings[0].size(); // Number of nodes in each ring
204  *nImperfectPrisms = 0; // Number of undeformed prisms
205  *nPerfectPrisms = 0; // Number of undeformed prisms
206  // Matrix for the reference ring for a given ringSize.
207  Eigen::MatrixXd refPointSet(ringSize, 3);
208 
209  // Get the reference ring point set for a given ring size.
210  // Get the axial dimension
211  int axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
212  yCloud->box.begin();
213  refPointSet = pntToPnt::getPointSetRefRing(ringSize, axialDim);
214  //
215 
216  // Two loops through all the rings are required to find pairs of basal rings
217  for (int iring = 0; iring < totalRingNum - 1; iring++) {
218  cond1 = false;
219  cond2 = false;
220  basal1 = rings[iring]; // Assign iring to basal1
221  // Loop through the other rings to get a pair
222  for (int jring = iring + 1; jring < totalRingNum; jring++) {
223  basal2 = rings[jring]; // Assign jring to basal2
224  // ------------
225  // Put extra check for axial basal rings if shapeMatching is being done
226  if (doShapeMatching == true || ringSize == 4) {
227  isAxialPair = false; // init
228  isAxialPair =
229  ring::discardExtraTetragonBlocks(&basal1, &basal2, yCloud);
230  if (isAxialPair == false) {
231  continue;
232  }
233  } // end of check for tetragonal prism blocks
234  // ------------
235  // Step one: Check to see if basal1 and basal2 have common
236  // elements or not. If they don't, then they cannot be basal rings
237  cond1 = ring::hasCommonElements(basal1, basal2);
238  if (cond1 == true) {
239  continue;
240  }
241  // -----------
242  // Step two and three: One of the elements of basal2 must be the nearest
243  // neighbour of the first (index0; l1) If m_k is the nearest neighbour of
244  // l1, m_{k+1} ... m_{k+(n-1)} must be neighbours of l_i+1 etc or l_i-1
245  cond2 = ring::basalPrismConditions(nList, &basal1, &basal2);
246  // If cond2 is false, the strict criteria for prisms has not been met
247  if (cond2 == false) {
248  // Skip if shape-matching is not desired
249  if (!doShapeMatching) {
250  continue;
251  } // shape-matching not desired
252  //
253  // If shape-matching is to be done:
254  // Check for the reduced criteria fulfilment
255  relaxedCond = ring::relaxedPrismConditions(nList, &basal1, &basal2);
256  // Skip if relaxed criteria are not met
257  if (relaxedCond == false) {
258  continue;
259  } // end of skipping if the prisms do not fulfil relaxed criteria
260 
261  // Do shape matching here
262  bool isDeformedPrism = match::matchPrism(
263  yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, false);
264 
265  // Success! The rings are basal rings of a deformed prism!
266  if (isDeformedPrism) {
267  //
268  // // Update the number of prism blocks
269  *nImperfectPrisms += 1;
270  // Update iring
271  if ((*ringType)[iring] == ring::unclassified) {
272  (*ringType)[iring] = ring::deformedPrism;
273  listPrism.push_back(iring);
274  } else if ((*ringType)[iring] == ring::Prism) {
275  (*ringType)[iring] = ring::mixedPrismRing;
276  } // if it is deformed
277  // Update jring
278  if ((*ringType)[jring] == ring::unclassified) {
279  (*ringType)[jring] = ring::deformedPrism;
280  listPrism.push_back(jring);
281  } else if ((*ringType)[jring] == ring::Prism) {
282  (*ringType)[jring] = ring::mixedPrismRing;
283  } // if it is deformed
284  } // end of update of ring types
285 
286  // // Write outs
287  // // Now write out axial basal rings
288  // sout::writeBasalRingsPrism(&basal1, &basal2, nDeformedPrisms, nList,
289  // yCloud, true);
290  } // end of reduced criteria
291  // Strict criteria
292  else {
293  // Update the number of prism blocks
294  *nPerfectPrisms += 1;
295  // Update iring
296  if ((*ringType)[iring] == ring::unclassified) {
297  (*ringType)[iring] = ring::Prism;
298  listPrism.push_back(iring);
299  } else if ((*ringType)[iring] == ring::deformedPrism) {
300  (*ringType)[iring] = ring::mixedPrismRing;
301  } // if it is deformed
302  // Update jring
303  if ((*ringType)[jring] == ring::unclassified) {
304  (*ringType)[jring] = ring::Prism;
305  listPrism.push_back(jring);
306  } else if ((*ringType)[jring] == ring::deformedPrism) {
307  (*ringType)[jring] = ring::mixedPrismRing;
308  } // if it is deformed
309  //
310  // Shape-matching to get the RMSD (if shape-matching is desired)
311  if (doShapeMatching) {
312  bool isKnownPrism = match::matchPrism(
313  yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, true);
314  } // end of shape-matching to get rmsd
315  //
316  // // Now write out axial basal rings for convex hull calculations
317  // sout::writePrisms(&basal1, &basal2, *nPrisms, yCloud);
318  // // Write out prisms for shape-matching
319  // sout::writeBasalRingsPrism(&basal1, &basal2, *nPrisms, nList, yCloud,
320  // false);
321  // -----------
322  } // end of strict criteria
323 
324  } // end of loop through rest of the rings to get the second basal ring
325  } // end of loop through all rings for first basal ring
326 
327  sort(listPrism.begin(), listPrism.end());
328  auto ip = std::unique(listPrism.begin(), listPrism.end());
329  // Resize peripheral rings to remove undefined terms
330  listPrism.resize(std::distance(listPrism.begin(), ip));
331 
332  return listPrism;
333 }
334 
348  std::vector<int> *basal1,
349  std::vector<int> *basal2) {
350  int l1 = (*basal1)[0]; // first element of basal1 ring
351  int ringSize =
352  (*basal1).size(); // Size of the ring; each ring contains n elements
353  int m_k; // Atom ID of element in basal2
354  bool l1_neighbour; // m_k is a neighbour of l1(true) or not (false)
355 
356  // isNeighbour is initialized to false for all basal2 elements; indication if
357  // basal2 elements are neighbours of basal1
358  std::vector<bool> isNeighbour(ringSize, false);
359  int kIndex; // m_k index
360  int lAtomID; // atomID of the current element of basal1
361  int kAtomID; // atomID of the current element of basal2
362 
363  // ---------------------------------------------
364  // COMPARISON OF basal2 ELEMENTS WITH l1
365  for (int k = 0; k < ringSize; k++) {
366  l1_neighbour = false;
367  m_k = (*basal2)[k];
368  // =================================
369  // Checking to seee if m_k is be a neighbour of l1
370  // Find m_k inside l1 neighbour list
371  auto it = std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
372 
373  // If the element has been found, for l1
374  if (it != nList[l1].end()) {
375  l1_neighbour = true;
376  kIndex = k;
377  break;
378  }
379  } // l1 is a neighbour of m_k
380 
381  // If there is no nearest neighbour, then the two rings are not part of the
382  // prism
383  if (!l1_neighbour) {
384  return false;
385  }
386 
387  // ---------------------------------------------
388  // NEIGHBOURS of basal1 in basal2
389  isNeighbour[kIndex] = true;
390 
391  // All elements of basal1 must be neighbours of basal2
392  for (int i = 1; i < ringSize; i++) {
393  lAtomID = (*basal1)[i]; // element of basal1 ring
394  for (int k = 0; k < ringSize; k++) {
395  // Skip if already a neighbour
396  if (isNeighbour[k]) {
397  continue;
398  }
399  // Get the comparison basal2 element
400  kAtomID = (*basal2)[k]; // element of basal2 ring;
401 
402  // Checking to see if kAtomID is a neighbour of lAtomID
403  // Find kAtomID inside lAtomID neighbour list
404  auto it1 =
405  std::find(nList[lAtomID].begin() + 1, nList[lAtomID].end(), kAtomID);
406 
407  // If the element has been found, for l1
408  if (it1 != nList[lAtomID].end()) {
409  isNeighbour[k] = true;
410  }
411  } // Loop through basal2
412  } // end of check for neighbours of basal1
413 
414  // ---------------------------------------------
415 
416  // They should all be neighbours
417  for (int k = 0; k < ringSize; k++) {
418  // Check to see if any element is false
419  if (!isNeighbour[k]) {
420  return false;
421  }
422  }
423 
424  // Everything works out!
425  return true;
426 }
427 
433  std::vector<int> *basal1,
434  std::vector<int> *basal2) {
435  int ringSize =
436  (*basal1).size(); // Size of the ring; each ring contains n elements
437  int m_k; // Atom ID of element in basal2
438  bool isNeighbour = false; // This is true if there is at least one bond
439  // between the basal rings
440  int l_k; // Atom ID of element in basal1
441 
442  // ---------------------------------------------
443  // COMPARISON OF basal2 ELEMENTS (m_k) WITH basal1 ELEMENTS (l_k)
444  // Loop through all the elements of basal1
445  for (int l = 0; l < ringSize; l++) {
446  l_k = (*basal1)[l];
447  // Search for the nearest neighbour of l_k in basal2
448  // Loop through basal2 elements
449  for (int m = 0; m < ringSize; m++) {
450  m_k = (*basal2)[m];
451  // Find m_k inside l_k neighbour list
452  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
453 
454  // If the element has been found, for l1
455  if (it != nList[l_k].end()) {
456  isNeighbour = true;
457  break;
458  } // found element
459  } // end of loop through all the elements of basal2
460 
461  // If a neighbour has been found then
462  if (isNeighbour) {
463  return true;
464  }
465  } // end of loop through all the elements of basal1
466 
467  // If a neighbour has not been found, return false
468  return false;
469 }
470 
484  std::vector<int> *basal1, std::vector<int> *basal2,
485  molSys::PointCloud<molSys::Point<double>, double> *yCloud) {
486  int ringSize =
487  (*basal1).size(); // Size of the ring; each ring contains n elements
488  int iatomIndex,
489  jatomIndex; // Indices of the elements in basal1 and basal2 respectively
490  double r_i, r_j; // Coordinates in the axial dimension of iatom and jatom of
491  // basal1 and basal2 respectively
492  int axialDim; // 0 for x, 1 for y and 2 for z dimensions respectively
493  // Variables for getting the projected area
494  bool axialBasal1, axialBasal2; // bools for checking if basal1 and basal2 are
495  // axial (true) respectively
496  double areaXY, areaXZ,
497  areaYZ; // Projected area on the XY, XZ and YZ planes respectively
498  // ----------------------------------------
499  // Find the axial dimension for a quasi-one-dimensional ice nanotube
500  // The axial dimension will have the largest box length
501  // Index -> axial dimension
502  // 0 -> x dim
503  // 1 -> y dim
504  // 2 -> z dim
505  axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
506  yCloud->box.begin();
507  // ----------------------------------------
508  // Calculate projected area onto the XY, YZ and XZ planes for basal1
509  axialBasal1 = false; // Init to false
510  axialBasal2 = false; // Init
511 
512  // Init the projected area
513  areaXY = 0.0;
514  areaXZ = 0.0;
515  areaYZ = 0.0;
516 
517  jatomIndex = (*basal1)[0];
518 
519  // All points except the first pair
520  for (int k = 1; k < ringSize; k++) {
521  iatomIndex = (*basal1)[k]; // Current vertex
522 
523  // Add to the polygon area
524  // ------
525  // XY plane
526  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
527  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
528  // ------
529  // XZ plane
530  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
531  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
532  // ------
533  // YZ plane
534  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
535  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
536  // ------
537  jatomIndex = iatomIndex;
538  }
539 
540  // Closure point
541  iatomIndex = (*basal1)[0];
542  // ------
543  // XY plane
544  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
545  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
546  // ------
547  // XZ plane
548  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
549  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
550  // ------
551  // YZ plane
552  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
553  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
554  // ------
555  // The actual projected area is half of this
556  areaXY *= 0.5;
557  areaXZ *= 0.5;
558  areaYZ *= 0.5;
559  // Get the absolute value
560  areaXY = fabs(areaXY);
561  areaXZ = fabs(areaXZ);
562  areaYZ = fabs(areaYZ);
563 
564  // If the axial dimension is x, y, or z:
565  // then the maximum basal area should be in the YZ, XZ and XY dimensions
566  // respectively
567  // x dim
568  if (axialDim == 0) {
569  if (areaYZ > areaXY && areaYZ > areaXZ) {
570  axialBasal1 = true;
571  } // end of check for axial ring for basal1
572  } // x dim
573  // y dim
574  else if (axialDim == 1) {
575  if (areaXZ > areaXY && areaXZ > areaYZ) {
576  axialBasal1 = true;
577  } // end of check for axial ring for basal1
578  } // x dim
579  // z dim
580  else if (axialDim == 2) {
581  if (areaXY > areaXZ && areaXY > areaYZ) {
582  axialBasal1 = true;
583  } // end of check for axial ring for basal1
584  } // x dim
585  else {
586  std::cerr << "Could not find the axial dimension.\n";
587  return false;
588  }
589  // ----------------------------------------
590  // Calculate projected area onto the XY, YZ and XZ planes for basal2
591 
592  // Init the projected area
593  areaXY = 0.0;
594  areaXZ = 0.0;
595  areaYZ = 0.0;
596 
597  jatomIndex = (*basal2)[0];
598 
599  // All points except the first pair
600  for (int k = 1; k < ringSize; k++) {
601  iatomIndex = (*basal2)[k]; // Current vertex
602 
603  // Add to the polygon area
604  // ------
605  // XY plane
606  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
607  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
608  // ------
609  // XZ plane
610  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
611  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
612  // ------
613  // YZ plane
614  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
615  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
616  // ------
617  jatomIndex = iatomIndex;
618  }
619 
620  // Closure point
621  iatomIndex = (*basal2)[0];
622  // ------
623  // XY plane
624  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
625  (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
626  // ------
627  // XZ plane
628  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
629  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
630  // ------
631  // YZ plane
632  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
633  (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
634  // ------
635  // The actual projected area is half of this
636  areaXY *= 0.5;
637  areaXZ *= 0.5;
638  areaYZ *= 0.5;
639  // Get the absolute value
640  areaXY = fabs(areaXY);
641  areaXZ = fabs(areaXZ);
642  areaYZ = fabs(areaYZ);
643 
644  // Check if xy projected area is the greatest
645  // If the axial dimension is x, y, or z:
646  // then the maximum basal area should be in the YZ, XZ and XY dimensions
647  // respectively
648  // x dim
649  if (axialDim == 0) {
650  if (areaYZ > areaXY && areaYZ > areaXZ) {
651  axialBasal2 = true;
652  } // end of check for axial ring for basal1
653  } // x dim
654  // y dim
655  else if (axialDim == 1) {
656  if (areaXZ > areaXY && areaXZ > areaYZ) {
657  axialBasal2 = true;
658  } // end of check for axial ring for basal1
659  } // x dim
660  // z dim
661  else if (axialDim == 2) {
662  if (areaXY > areaXZ && areaXY > areaYZ) {
663  axialBasal2 = true;
664  } // end of check for axial ring for basal1
665  } // x dim
666  else {
667  std::cerr << "Could not find the axial dimension.\n";
668  return false;
669  }
670  // ----------------------------------------
671 
672  // Now check if basal1 and basal2 are axial or not
673  if (axialBasal1 == true && axialBasal2 == true) {
674  return true;
675  } else {
676  return false;
677  } // Check for basal1 and basal2
678 }
679 
693  std::vector<int> listPrism, int ringSize,
695  std::vector<int> *atomTypes,
696  std::vector<ring::strucType> *atomState) {
697  // Every value in listPrism corresponds to an index in rings.
698  // Every ring contains atom indices, corresponding to the indices (not atom
699  // IDs) in rings
700  int iring; // Index of current ring
701  int iatom; // Index of current atom
702  ring::strucType currentState; // Current state of the ring being filled
703 
704  // Dummy value corresponds to a value of 1.
705  // Each value is initialized to the value of 1.
706 
707  // Loop through every ring in rings
708  for (int i = 0; i < listPrism.size(); i++) {
709  iring = listPrism[i];
710  // Get the state of all atoms in the ring
711  currentState = ringType[iring];
712  // Loop through every element in iring
713  for (int j = 0; j < ringSize; j++) {
714  iatom = rings[iring][j]; // Atom index
715  // Update the atom type
716  (*atomTypes)[iatom] = ringSize;
717  //
718  // Update the state of the atom
719  if ((*atomState)[iatom] == ring::unclassified) {
720  (*atomState)[iatom] = currentState;
721  } // Update the unclassified atom
722  else {
723  if ((*atomState)[iatom] != currentState) {
724  // For mixed, there is a preference
725  if (currentState == ring::mixedPrismRing) {
726  (*atomState)[iatom] = currentState;
727  } // fill
728  else if ((*atomState)[iatom] == ring::deformedPrism &&
729  currentState == ring::Prism) {
730  (*atomState)[iatom] = ring::mixedPrismRing;
731  } else if ((*atomState)[iatom] == ring::Prism &&
732  currentState == ring::deformedPrism) {
733  (*atomState)[iatom] = ring::mixedPrismRing;
734  }
735  } //
736  } // already filled?
737 
738  } // end of loop through every atom in iring
739  } // end of loop through every ring
740 
741  return 0;
742 } // end of function
743 
744 // Get the atom type values for deformed prisms
756  std::vector<int> *atomTypes, int maxDepth) {
757  //
758  int nop = atomState.size(); // Number of particles
759 
760  // Loop through all particles
761  for (int iatom = 0; iatom < nop; iatom++) {
762  // Check the atom state
763  // Deformed
764  if (atomState[iatom] == ring::deformedPrism) {
765  (*atomTypes)[iatom] += maxDepth - 2;
766  } // type for a deformed prism atom
767  else if (atomState[iatom] == ring::mixedPrismRing) {
768  (*atomTypes)[iatom] = 2;
769  } // type for a mixed prism ring
770  } // end of reassignation
771 
772  //
773  return 0;
774 }
775 
776 // Shift the entire ice nanotube and remove axial translations
788  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int *atomID,
789  int firstFrame, int currentFrame) {
790  //
791  int atomIndex; // Index of the atom to be shifted first
792  double shiftDistance; // Value by which all atoms will be shifted
793  double distFrmBoundary; // Distance from either the upper or lower boundary
794  // along the axial dimension
795  // Get the axial dimension
796  int axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
797  yCloud->box.begin();
798  // Check: (default z)
799  if (axialDim < 0 || axialDim > 2) {
800  axialDim = 2;
801  } // end of check to set the axial dimension
802  // Lower and higher limits of the box in the axial dimension
803  double boxLowAxial = yCloud->boxLow[axialDim];
804  double boxHighAxial = boxLowAxial + yCloud->box[axialDim];
805  //
806  // -----------------------------------
807  // Update the atomID of the 'first' atom in yCloud if the current frame is the
808  // first frame.
809  // Get the atom index of the first atom to be shifted down or up
810  if (currentFrame == firstFrame) {
811  *atomID = yCloud->pts[0].atomID;
812  atomIndex = 0;
813  } // end of update of the atom ID to be shifted for the first frame
814  else {
815  //
816  int iatomID = *atomID; // Atom ID of the first particle to be moved down
817  // Find the index given the atom ID
818  auto it = yCloud->idIndexMap.find(iatomID);
819 
820  if (it != yCloud->idIndexMap.end()) {
821  atomIndex = it->second;
822  } // found jatom
823  else {
824  std::cerr << "Lost atoms.\n";
825  return 1;
826  } // error handling
827  //
828  } // Update of atomIndex for all other frames
829  // -----------------------------------
830  // Calculate the distance by which the atoms must be shifted (negative value)
831  if (axialDim == 0) {
832  shiftDistance = boxLowAxial - yCloud->pts[atomIndex].x;
833  } // x coordinate
834  else if (axialDim == 1) {
835  shiftDistance = boxLowAxial - yCloud->pts[atomIndex].y;
836  } // y coordinate
837  else {
838  shiftDistance = boxLowAxial - yCloud->pts[atomIndex].z;
839  } // z coordinate
840  // -----------------------------------
841  // Shift all the particles
842  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
843  // Shift the particles along the axial dimension only
844  if (axialDim == 0) {
845  yCloud->pts[iatom].x += shiftDistance;
846  // Wrap-around
847  if (yCloud->pts[iatom].x < boxLowAxial) {
848  distFrmBoundary = boxLowAxial - yCloud->pts[iatom].x; // positive value
849  yCloud->pts[iatom].x = boxHighAxial - distFrmBoundary;
850  } // end of wrap-around
851  } // x coordinate
852  else if (axialDim == 1) {
853  yCloud->pts[iatom].y += shiftDistance;
854  // Wrap-around
855  if (yCloud->pts[iatom].y < boxLowAxial) {
856  distFrmBoundary = boxLowAxial - yCloud->pts[iatom].y; // positive value
857  yCloud->pts[iatom].y = boxHighAxial - distFrmBoundary;
858  } // end of wrap-around
859  } // y coordinate
860  else {
861  yCloud->pts[iatom].z += shiftDistance;
862  // Wrap-around
863  if (yCloud->pts[iatom].z < boxLowAxial) {
864  distFrmBoundary = boxLowAxial - yCloud->pts[iatom].z; // positive value
865  yCloud->pts[iatom].z = boxHighAxial - distFrmBoundary;
866  } // end of wrap-around
867  } // z coordinate
868  } // end of shifting all paritcles downward
869  // -----------------------------------
870  return 0;
871 } // end of function
T begin(T... args)
T distance(T... args)
T end(T... args)
T find(T... args)
bool relaxedPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
int 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)
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)
strucType
Definition: ring.hpp:112
int deformedPrismTypes(std::vector< ring::strucType > atomState, std::vector< int > *atomTypes, int maxDepth)
Get the atom type values for deformed prisms.
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
Definition: ring.cpp:18
bool discardExtraTetragonBlocks(std::vector< int > *basal1, std::vector< int > *basal2, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
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)
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Definition: ring.cpp:175
bool basalPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
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
@ deformedPrism
Definition: ring.hpp:120
@ unclassified
The ring is unclassified, which may be either water or a deformed type which cannot be classified by ...
Definition: ring.hpp:113
@ mixedPrismRing
Definition: ring.hpp:121
@ Prism
The ring belongs to a prism block, classified according to the prism identification scheme.
Definition: ring.hpp:119
T max_element(T... args)
bool matchPrism(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, std::vector< double > *rmsdPerAtom, bool isPerfect=true)
Definition: shapeMatch.cpp:17
Eigen::MatrixXd getPointSetRefRing(int n, int axialDim)
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)
T push_back(T... args)
T resize(T... args)
T size(T... args)
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166
This contains per-particle information.
Definition: mol_sys.hpp:145
File containing functions used specific to quasi-one-dimensional topological network critera (the pri...
T unique(T... args)