Clump

Functions

int clump::largestIceCluster (std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int >> nList, std::vector< bool > *isIce, std::vector< int > *clusterID, std::vector< int > *nClusters, std::unordered_map< int, int > *indexNumber, int firstFrame)
 Finds the largest ice cluster. More...
 
int clump::singleClusterLinkedList (molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int >> nList, std::vector< int > *linkedList)
 
int clump::clusterAnalysis (std::string path, molSys::PointCloud< molSys::Point< double >, double > *iceCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< std::vector< int >> &iceNeighbourList, double cutoff, int firstFrame, std::string bopAnalysis="q6")
 
int clump::recenterClusterCloud (molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int >> nList)
 Recenters the coordinates of a pointCloud. More...
 

Detailed Description

Function Documentation

◆ clusterAnalysis()

int clump::clusterAnalysis ( std::string  path,
molSys::PointCloud< molSys::Point< double >, double > *  iceCloud,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< std::vector< int >>  nList,
std::vector< std::vector< int >> &  iceNeighbourList,
double  cutoff,
int  firstFrame,
std::string  bopAnalysis = "q6" 
)

Does the cluster analysis of ice particles in the system. Returns a pointCloud of the largest ice cluster. The neighbour list returned is BY INDEX of the largest ice cluster pointCloud.

Does the cluster analysis of ice particles in the system. Returns a molSys::PointCloud of the largest ice cluster (using the \( q_6 \) parameter by default). Uses the full neighbour list (by ID) according to the full PointCloud yCloud. Returns a neighbour list by index, according to the largest ice cluster.

Parameters
[in]pathOutput directory path, specified by the user
[in,out]iceCloudThe molSys::PointCloud for the largest ice cluster of the ice-like molecules
[in]yCloudThe molSys::PointCloud for all the particles in the frame, regardless of ice type
[in]nListRow-ordered neighbour list by atom ID
[in]iceNeighbourListRow-ordered neighbour list by atom index, not ID, according to the iceCloud atoms
[in]cutoffCutoff for the nearest neighbours
[in]firstFrameFirst frame to be analyzed
[in]bopAnalysisThis determines which method to use for determining the ice-like nature of the particles. This can be "q6" or "chill", for using the \( q_6 \) parameter or CHILL algorithm, respectively

Definition at line 312 of file cluster.cpp.

318  {
319  //
320  std::vector<bool> isIce; // For every particle in yCloud, has a value
321  int nTotalIce; // Total number of ice-like molecules
322  std::vector<int> clusterID; // Vector containing the cluster IDs of all the
323  // ice-like molecules.
324  std::vector<int> nClusters; // Number of clusters for each index
326  indexNumber; // Map for cluster index and index in the number vector
327 
328  // -------------------------------------------------------
329  // Init
330  // Clear the largest ice cluster pointCloud.
331  *iceCloud = molSys::clearPointCloud(iceCloud);
332  // Clear the neighbour list by index
333  nneigh::clearNeighbourList(iceNeighbourList);
334  // Init the vector of bools for every particle in yCloud
335  isIce.resize(yCloud->nop);
336  nTotalIce = 0; // Total number of ice-like molecules
337  // -------------------------------------------------------
338  // Use a bond-orientational parameter to find ice-like particles
339  // Q6
340  if (bopAnalysis == "q6") {
341  //
342  std::vector<double> q6Values;
343  // Calculate the Q6 parameters for every point in yCloud
344  q6Values = chill::getq6(yCloud, nList);
345  // Assign values to isIce according to the values
346  // of the q6 parameter. If q6 is greater than 0.5, it is ice-like.
347  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
348  // If q6 is greater than 0.5, it is ice-like
349  if (q6Values[iatom] > 0.5) {
350  isIce[iatom] = true; // is ice-like; by default false
351  nTotalIce++; // Add to the number of ice-like molecules
352  } // ice-like molecule found
353  } // end of loop through every atom in yCloud
354  } // end of getting a vector of bools for ice-like particles
355  //
356  // CHILL
357  // Q6
358  if (bopAnalysis == "chill") {
359  //
360  *yCloud = chill::getCorrel(yCloud, nList, false);
361  // Get the ice types
362  *yCloud = chill::getIceTypeNoPrint(yCloud, nList, false);
363  // Assign values to isIce according to the CHILL algorithm
364  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
365  // If it is an ice-like molecule, add it, otherwise skip
366  if (yCloud->pts[iatom].iceType == molSys::atom_state_type::water) {
367  continue;
368  } // water
369  if (yCloud->pts[iatom].iceType == molSys::atom_state_type::unclassified) {
370  continue;
371  } // unclassified
372  isIce[iatom] = true; // is ice-like; by default false
373  nTotalIce++; // Add to the number of ice-like molecules
374 
375  } // end of loop through every atom in yCloud
376  } // end of getting a vector of bools for ice-like particles
377  // -------------------------------------------------------
378  // Get the largest ice cluster and other data
379  clump::largestIceCluster(path, yCloud, iceCloud, nList, &isIce, &clusterID,
380  &nClusters, &indexNumber, firstFrame);
381 
382  // -------------------------------------------------------
383  // Get the neighbour list by index according to the largest ice cluster
384  // pointCloud
385  iceNeighbourList = nneigh::getNewNeighbourListByIndex(iceCloud, cutoff);
386 
387  return 0;
388 } // end of function
std::vector< double > getq6(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Definition: bop.cpp:868
molSys::PointCloud< molSys::Point< double >, double > getCorrel(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Definition: bop.cpp:299
molSys::PointCloud< molSys::Point< double >, double > getIceTypeNoPrint(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Classifies each atom according to the CHILL algorithm without printing.
Definition: bop.cpp:437
int largestIceCluster(std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int >> nList, std::vector< bool > *isIce, std::vector< int > *clusterID, std::vector< int > *nClusters, std::unordered_map< int, int > *indexNumber, int firstFrame)
Finds the largest ice cluster.
Definition: cluster.cpp:40
int nop
Current frame number.
Definition: mol_sys.hpp:173
std::vector< S > pts
Definition: mol_sys.hpp:171
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
@ water
Liquid/amorphous phase.
@ unclassified
Not classified into any other category.
int clearNeighbourList(std::vector< std::vector< int >> &nList)
Erases memory for a vector of vectors for the neighbour list.
Definition: neighbours.cpp:391
std::vector< std::vector< int > > getNewNeighbourListByIndex(molSys::PointCloud< molSys::Point< double >, double > *yCloud, double cutoff)
Definition: neighbours.cpp:294
T resize(T... args)

◆ largestIceCluster()

int clump::largestIceCluster ( std::string  path,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
molSys::PointCloud< molSys::Point< double >, double > *  iceCloud,
std::vector< std::vector< int >>  nList,
std::vector< bool > *  isIce,
std::vector< int > *  list,
std::vector< int > *  nClusters,
std::unordered_map< int, int > *  indexNumber,
int  firstFrame 
)

Finds the largest ice cluster.

Finds the number of particles in the largest ice cluster, for a given frame, using Stoddard's clustering algorithm (Stoddard J. Comp. Phys., 27, 291, 1977)](https://www.sciencedirect.com/science/article/pii/0021999178900116)

Parameters
[in]pathOutput directory path, specified by the user
[in]yCloudThe molSys::PointCloud for all particles, regardless of type classification
[in,out]iceCloudThe molSys::PointCloud for the largest ice cluster of the ice-like molecules
[in]nListRow-ordered neighbour list by atom ID
[in]isIceHolds a bool value for each particle in yCloud. This is true for ice-like molecules and false otherwise
[in]listLinked list created by the clustering algorithm.
[in]nClustersContains the number of particles in every ice-like cluster found
[in]indexNumberUnordered map for mapping the cluster ID indices to the number in each cluster
[in]firstFrameFirst frame to be analyzed

Definition at line 40 of file cluster.cpp.

45  {
46  //
47  int kAtomID; // Atom ID of the nearest neighbour
48  int iClusterNumber; // Number in the current cluster
49  int nLargestCluster; // Number in the largest cluster
50  std::vector<int> linkedList; // Linked list
51  int j; // Index
52  std::vector<int> startingIndex; // Contains the vector index in list
53  // corresponding to a particular cluster
54  int temp; // For swapping
56  visited; // To make sure you don't go through the same atoms again.
57  molSys::Point<double> iPoint; // Current point
58  int currentIndex; // Current index
59 
60  // -----------------------------------------------------------
61  // INITIALIZATION
62  linkedList.resize(yCloud->nop, -1); // init to dummy value
63  // Initial values of the list. -1 is a dummy value if the molecule is
64  // water or not in the slice
65  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
66  // Skip if the molecule is water or if it is not in the slice
67  if ((*isIce)[iatom] == false || yCloud->pts[iatom].inSlice == false) {
68  continue;
69  } // skip for water or not in slice
70  // Otherwise, assign the index as the ID
71  linkedList[iatom] = iatom;
72  } // init of cluster IDs
73  // -----------------------------------------------------------
74  // Get the linked list
75  for (int i = 0; i < yCloud->nop - 1; i++) {
76  //
77  // Skip if the molecule is water or if it is not in the slice
78  if (linkedList[i] == -1) {
79  continue;
80  } // skip for water or not in slice
81  //
82  // If iatom is already in a cluster, skip it
83  if (linkedList[i] != i) {
84  continue;
85  } // Already part of a cluster
86  //
87  j = i; // Init of j
88  // Execute the next part of the loop while j is not equal to i
89  do {
90  //
91  // Go through the rest of the atoms (KLOOP)
92  for (int k = i + 1; k < yCloud->nop; k++) {
93  // Skip if not ice
94  if ((*isIce)[k] == false) {
95  continue;
96  } // not ice
97  // Skip if already part of a cluster
98  if (linkedList[k] != k) {
99  continue;
100  } // Already part of a cluster
101  //
102  // Check to see if k is a nearest neighbour of j
103  kAtomID = yCloud->pts[k].atomID; // Atom ID
104  auto it = std::find(nList[j].begin() + 1, nList[j].end(), kAtomID);
105  if (it != nList[j].end()) {
106  // Swap!
107  temp = linkedList[j];
108  linkedList[j] = linkedList[k];
109  linkedList[k] = temp;
110  } // j and k are nearest neighbours
111  } // end of loop through k (KLOOP)
112  //
113  j = linkedList[j];
114  } while (j != i); // end of control for j!=i
115  //
116 
117  } // end of looping through every i
118 
119  // -----------------------------------------------------------
120  // Get the starting index (which is the vector index in list) corresponding to
121  // clusters with more than one molecule in them
122  int nextElement; // value in list
123  int index; // starting index value
124  // init
125  visited.resize(yCloud->nop);
126 
127  for (int i = 0; i < yCloud->nop; i++) {
128  //
129  if (visited[i]) {
130  continue;
131  } // Already counted
132  visited[i] = true; // Visited
133  if (linkedList[i] == -1) {
134  continue;
135  } // not ice
136  if (linkedList[i] == i) {
137  continue;
138  } // only one element
139  //
140  currentIndex = i;
141  nextElement = linkedList[currentIndex];
142  index = i;
143  iClusterNumber = 1; // at least one value
144  while (nextElement != index) {
145  iClusterNumber++;
146  currentIndex = nextElement;
147  visited[currentIndex] = true;
148  nextElement = linkedList[currentIndex];
149  } // get number
150  // Update startingIndex with index
151  startingIndex.push_back(index);
152  // Update the number of molecules in the cluster
153  (*nClusters).push_back(iClusterNumber);
154  } // end of loop through
155  // -----------------------------------------------------------
156  // Get the largest cluster and save the atoms to the iceCloud pointCloud
157  nLargestCluster = *std::max_element((*nClusters).begin(), (*nClusters).end());
158  int lClusIndex = distance(
159  (*nClusters).begin(),
160  max_element(
161  (*nClusters).begin(),
162  (*nClusters).end())); // index of the cluster with the largest number
163  // -----------------------------------------------------------
164  // Loop through the linked list from the starting element
165  // startingIndex[lClusIndex].
166  int startCluster =
167  startingIndex[lClusIndex]; // index in linkedList to start from
168 
169  // L[i]->j->L[j]->k->L[k]->i
170  index = startCluster;
171  // Update iceCloud
172  iPoint = yCloud->pts[index];
173  iceCloud->pts.push_back(iPoint);
174  //
175  nextElement = linkedList[startCluster];
176  while (nextElement != index) {
177  currentIndex = nextElement;
178  // update with currentIndex
179  iPoint = yCloud->pts[currentIndex];
180  iceCloud->pts.push_back(iPoint);
181  // end update
182  nextElement = linkedList[currentIndex];
183  } // get the largest cluster
184  // -----------------------------------------------------------
185  // Update other variables in iceCloud
186 
187  iceCloud->currentFrame = yCloud->currentFrame;
188  iceCloud->nop = iceCloud->pts.size();
189  iceCloud->box = yCloud->box;
190  iceCloud->boxLow = yCloud->boxLow;
191 
192  // Update idIndexMap
193  for (int iatom = 0; iatom < iceCloud->nop; iceCloud++) {
194  iceCloud->idIndexMap[iceCloud->pts[iatom].atomID] = iatom;
195  } // end of loop through iceCloud
196 
197  // -----------------------------------------------------------
198  // Write out the cluster statistics
199  int totalClusters = (*nClusters).size(); // Total number of clusters
200  int smallestCluster = nLargestCluster =
201  *std::min_element((*nClusters).begin(),
202  (*nClusters).end()); // Size of the smallest cluster
203  double avgClusterSize = 0.0;
204 
205  // Get the average cluster size
206  for (int iCluster = 0; iCluster < totalClusters; iCluster++) {
207  avgClusterSize += (*nClusters)[iCluster];
208  } // Loop through the clusters
209  // Normalize by the number
210  if (totalClusters == 0) {
211  avgClusterSize = 0;
212  } else {
213  avgClusterSize /= totalClusters;
214  }
215 
216  iceCloud->currentFrame = yCloud->currentFrame;
217  nLargestCluster = (*nClusters)[lClusIndex];
218 
219  // Write out to the file
220  sout::writeClusterStats(path, yCloud->currentFrame, nLargestCluster,
221  totalClusters, smallestCluster, avgClusterSize,
222  firstFrame);
223 
224  // -----------------------------------------------------------
225  return 0;
226 }
T begin(T... args)
T end(T... args)
T find(T... args)
double distance(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Inline generic function for obtaining the wrapped distance between two particles WITHOUT applying PBC...
Definition: generic.hpp:170
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:176
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:174
std::vector< T > boxLow
Periodic box lengths.
Definition: mol_sys.hpp:175
int currentFrame
Collection of points.
Definition: mol_sys.hpp:172
T max_element(T... args)
T min_element(T... args)
int writeClusterStats(std::string path, int currentFrame, int largestCluster, int numOfClusters, int smallestCluster, double avgClusterSize, int firstFrame)
Function for writing out cluster statistics.
T push_back(T... args)
T size(T... args)
This contains per-particle information.
Definition: mol_sys.hpp:149

◆ recenterClusterCloud()

int clump::recenterClusterCloud ( molSys::PointCloud< molSys::Point< double >, double > *  iceCloud,
std::vector< std::vector< int >>  nList 
)

Recenters the coordinates of a pointCloud.

Recenters the largest ice cluster, by applying a transformation on the largest ice cluster coordinates. Requires the neighbour list BY INDEX.

Parameters
[in]iceCloudThe molSys::PointCloud for the largest ice cluster of the ice-like molecules
[in]nListRow-ordered neighbour list by atom index, for the molSys::PointCloud iceCloud

Definition at line 398 of file cluster.cpp.

400  {
401  //
402  int dim = 3; // Dimensions
403  std::vector<double> box = iceCloud->box;
404  std::vector<double> boxLow = iceCloud->boxLow;
405  std::vector<double> boxHigh;
406  double xBoxCenter, yBoxCenter, zBoxCenter; // Centroid of the simulation box
407  double x_centroid, y_centroid, z_centroid; // Centroid of the cluster
408  // Variables for the linked list
409  std::vector<int> linkedList; // Contains the linked list for the cluster
410  std::vector<bool> visited; // Records whether an item has been visited or not
411 
412  // To avoid long confusing lines, fill boxHigh
413  for (int k = 0; k < dim; k++) {
414  boxHigh.push_back(boxLow[k] + box[k]);
415  } // end of filling up boxHigh
416 
417  // --------------------------------------------------------------------------
418  // Get the linked list of the cluster
419  clump::singleClusterLinkedList(iceCloud, nList, &linkedList);
420  // --------------------------------------------------------------------------
421  // Loop through the entire looped list
422  // init
423  visited.resize(iceCloud->nop);
424 
425  // The starting value is the first atom
426  int iatom = 0; // Atom index of the 'starting value'
427  int currentIndex; // Current atom
428  int nextElement; // Next linked atom
429  int index; // Keeps track of the first element in the linked list
430  double x_ij, y_ij, z_ij; // Relative distance between the two atoms
431  double xPBC, yPBC, zPBC; // Actual distance
432 
433  // Loop through the entire linked list
434  for (int i = 0; i < iceCloud->nop; i++) {
435  //
436  if (visited[i]) {
437  continue;
438  } // Already counted
439  visited[i] = true; // Visited
440  // Should never execute
441  if (linkedList[i] == i) {
442  continue;
443  } // only one element
444  //
445  currentIndex = i;
446  nextElement = linkedList[currentIndex];
447  index = i;
448  // Keep looping
449  while (nextElement != index) {
450  currentIndex = nextElement;
451  visited[currentIndex] = true;
452  nextElement = linkedList[currentIndex];
453  // -----------------------------------
454  // Get the relative distance between the central atom (iatom)
455  // and the next element
456  // Coordinates
457  // if (nextElement != index) {
458  x_ij = iceCloud->pts[currentIndex].x - iceCloud->pts[nextElement].x;
459  y_ij = iceCloud->pts[currentIndex].y - iceCloud->pts[nextElement].y;
460  z_ij = iceCloud->pts[currentIndex].z - iceCloud->pts[nextElement].z;
461  // Shift the nextElement if it's on the other side of the box
462  // Shift x
463  if (fabs(x_ij) > 0.5 * box[0]) {
464  // Get the actual distance
465  xPBC = box[0] - fabs(x_ij);
466  if (x_ij < 0) {
467  iceCloud->pts[nextElement].x = iceCloud->pts[currentIndex].x - xPBC;
468  } // To the -x side of currentIndex
469  else {
470  iceCloud->pts[nextElement].x = iceCloud->pts[currentIndex].x + xPBC;
471  } // Add to the + side
472  } // Shift nextElement
473  //
474  // Shift y
475  if (fabs(y_ij) > 0.5 * box[1]) {
476  // Get the actual distance
477  yPBC = box[1] - fabs(y_ij);
478  if (y_ij < 0) {
479  iceCloud->pts[nextElement].y = iceCloud->pts[currentIndex].y - yPBC;
480  } // To the -y side of currentIndex
481  else {
482  iceCloud->pts[nextElement].y = iceCloud->pts[currentIndex].y + yPBC;
483  } // Add to the + side
484  } // Shift nextElement
485  //
486  // Shift z
487  if (fabs(z_ij) > 0.5 * box[2]) {
488  // Get the actual distance
489  zPBC = box[2] - fabs(z_ij);
490  if (z_ij < 0) {
491  iceCloud->pts[nextElement].z = iceCloud->pts[currentIndex].z - zPBC;
492  } // To the -z side of currentIndex
493  else {
494  iceCloud->pts[nextElement].z = iceCloud->pts[currentIndex].z + zPBC;
495  } // Add to the + side
496  } // Shift nextElement
497  // -----------------------------------
498  // } // don't shift the last atom!
499 
500  } // End of going through linked atoms
501  //
502  } // end of loop through atoms
503  // --------------------------------------------------------------------------
504  // Center of the simulation box
505  xBoxCenter = 0.5 * (boxLow[0] + boxHigh[0]);
506  yBoxCenter = 0.5 * (boxLow[1] + boxHigh[1]);
507  zBoxCenter = 0.5 * (boxLow[2] + boxHigh[2]);
508  // --------------------------------------------------------------------------
509  // Get the centroid of the ice cluster.
510  x_centroid = 0.0;
511  y_centroid = 0.0;
512  z_centroid = 0.0;
513 
514  for (int i = 0; i < iceCloud->nop; i++) {
515  x_centroid += iceCloud->pts[i].x;
516  y_centroid += iceCloud->pts[i].y;
517  z_centroid += iceCloud->pts[i].z;
518  } // end of loop through the particles
519 
520  if (iceCloud->nop == 0) {
521  std::cerr << "There are no particles in the cluster.\n";
522  return 1;
523  } // error
524 
525  // Normalize by the number of particles.
526  x_centroid /= iceCloud->nop;
527  y_centroid /= iceCloud->nop;
528  z_centroid /= iceCloud->nop;
529 
530  // --------------------------------------------------------------------------
531  // Get the distance to shift by:
532  double xShift = x_centroid - xBoxCenter;
533  double yShift = y_centroid - yBoxCenter;
534  double zShift = z_centroid - zBoxCenter;
535 
536  // Loop through all atoms and shift them
537  for (int i = 0; i < iceCloud->nop; i++) {
538  iceCloud->pts[i].x -= xShift;
539  iceCloud->pts[i].y -= yShift;
540  iceCloud->pts[i].z -= zShift;
541  } // end of loop through the particles
542 
543  return 0;
544 } // end of function
T fabs(T... args)
int singleClusterLinkedList(molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int >> nList, std::vector< int > *linkedList)
Definition: cluster.cpp:237

◆ singleClusterLinkedList()

int clump::singleClusterLinkedList ( molSys::PointCloud< molSys::Point< double >, double > *  iceCloud,
std::vector< std::vector< int >>  nList,
std::vector< int > *  linkedList 
)

Get the linked list of a cluster, given by iceCloud, for a single cluster. Required for cluster re-centering

Get the linked list of a cluster, given by iceCloud, for a single cluster. Required for cluster re-centering

Parameters
[in]iceCloudThe molSys::PointCloud for the largest ice cluster of the ice-like molecules
[in]nListRow-ordered neighbour list by atom index, not ID
[in,out]linkedListLinked list created by the clustering algorithm for the largest ice cluster

Definition at line 237 of file cluster.cpp.

239  {
240  //
241  int j;
242  int temp; // For swapping indices
243  //
244  // -----------------------------------------------------------
245  // INITIALIZATION
246  (*linkedList).resize(iceCloud->nop);
247  // Initial values of the list.
248  for (int iatom = 0; iatom < iceCloud->nop; iatom++) {
249  // Assign the index as the ID
250  (*linkedList)[iatom] = iatom;
251  } // init of cluster IDs
252  // -----------------------------------------------------------
253  // Get the linked list
254  for (int i = 0; i < iceCloud->nop - 1; i++) {
255  //
256  // If iatom is already in a cluster, skip it
257  if ((*linkedList)[i] != i) {
258  continue;
259  } // Already part of a cluster
260  //
261  j = i; // Init of j
262  // Execute the next part of the loop while j is not equal to i
263  do {
264  //
265  // Go through the rest of the atoms (KLOOP)
266  for (int k = i + 1; k < iceCloud->nop; k++) {
267  // Skip if already part of a cluster
268  if ((*linkedList)[k] != k) {
269  continue;
270  } // Already part of a cluster
271  //
272  // Check to see if k is a nearest neighbour of j
273  auto it = std::find(nList[j].begin() + 1, nList[j].end(), k);
274  if (it != nList[j].end()) {
275  // Swap!
276  temp = (*linkedList)[j];
277  (*linkedList)[j] = (*linkedList)[k];
278  (*linkedList)[k] = temp;
279  } // j and k are nearest neighbours
280  } // end of loop through k (KLOOP)
281  //
282  j = (*linkedList)[j];
283  } while (j != i); // end of control for j!=i
284  //
285 
286  } // end of looping through every i
287 
288  // Return 0
289  return 0;
290 } // end of function