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 308 of file cluster.cpp.

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

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

◆ 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 394 of file cluster.cpp.

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

◆ 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 233 of file cluster.cpp.

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