neighbours.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 <iostream>
12 #include <math.h>
13 #include <neighbours.hpp>
14 
26 nneigh::neighList(double rcutoff,
28  int typeI, int typeJ) {
29  std::vector<std::vector<int>> nList; // Vector of vector of ints
30  int jatomIndex; // Atom ID corresponding to jatom
31  int iatomIndex; // Atom ID corresponding to iatom
32  double r_ij; // cutoff
33 
34  // Initialize with nop (irrespective of type)
35  // Initialize and fill the first element with the current atom ID whose
36  // neighbour list will be filled
37  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
38  // Find the atom ID (key) given the index or iatom (value)
39  auto itr = std::find_if(
40  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
41  [&iatom](const std::pair<int, int> &p) { return p.second == iatom; });
42  // If found:
43  if (itr == yCloud->idIndexMap.end()) {
44  std::cerr << "Something is wrong with your idIndexMap!\n";
45  continue;
46  } else {
47  iatomIndex = itr->first;
48  } // End of finding the atom ID to fill as the first element in the
49  // neighbour list
50  nList.push_back(std::vector<int>()); // Empty vector for the index iatom
51  // Fill the first element with the atom ID of iatom itself
52  nList[iatom].push_back(iatomIndex);
53  } // end of init
54 
55  // pairs of atoms of type I and J
56  // Loop through every iatom and find nearest neighbours within rcutoff
57  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
58  if (yCloud->pts[iatom].type != typeI) {
59  continue;
60  }
61  // Loop through the other atoms
62  for (int jatom = 0; jatom < yCloud->nop; jatom++) {
63  if (yCloud->pts[jatom].type != typeJ) {
64  continue;
65  }
66  // If the distance is greater than rcutoff, continue
67  r_ij = gen::periodicDist(yCloud, iatom, jatom);
68  if (r_ij > rcutoff) {
69  continue;
70  }
71 
72  // Get the atom IDs for iatom and jatom
73  auto gotI = std::find_if(
74  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
75  [&iatom](const std::pair<int, int> &p) { return p.second == iatom; });
76  if (gotI == yCloud->idIndexMap.end()) {
77  std::cerr << "Something is wrong with your idIndexMap!\n";
78  return nList;
79  } else {
80  iatomIndex = gotI->first;
81  } // End of finding the atom ID for iatom
82  // Find the atom ID of jatom
83  auto gotJ = std::find_if(
84  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
85  [&jatom](const std::pair<int, int> &p) { return p.second == jatom; });
86  if (gotJ == yCloud->idIndexMap.end()) {
87  std::cerr << "Something is wrong with your idIndexMap!\n";
88  return nList;
89  } else {
90  jatomIndex = gotJ->first;
91  } // End of finding the atom ID for jatom
92  // Update the neighbour indices with atom IDs for iatom and jatom both
93  // (full list)
94  nList[iatom].push_back(jatomIndex);
95  nList[jatom].push_back(iatomIndex);
96 
97  } // End of loop through jatom
98  } // End of loop for iatom
99 
100  return nList;
101 }
102 
114 nneigh::neighListO(double rcutoff,
115  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
116  int typeI) {
118  nList; // Vector of vectors of the neighbour list
119  double r_ij; // Distance between iatom and jatom
120  int iatomIndex; // Atomic ID of the atom with index iatom
121  int jatomIndex; // Atomic ID of the atom with index jatom
122  int indexYay;
123  std::vector<int> tempListIatom;
124 
125  // Initialize and fill the first element with the current atom ID whose
126  // neighbour list will be filled
127  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
128  // Find the atom ID (key) given the index or iatom (value)
129  auto itr = std::find_if(
130  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
131  [&iatom](const std::pair<int, int> &p) { return p.second == iatom; });
132  // If found:
133  if (itr == yCloud->idIndexMap.end()) {
134  std::cerr << "Something is wrong with your idIndexMap!\n";
135  continue;
136  } else {
137  iatomIndex = itr->first;
138  } // End of finding the atom ID to fill as the first element in the
139  // neighbour list
140  nList.push_back(std::vector<int>()); // Empty vector for the index iatom
141  // Fill the first element with the atom ID of iatom itself
142  nList[iatom].push_back(iatomIndex);
143  } // end of init
144 
145  // Loop through every iatom and find nearest neighbours within rcutoff
146  for (int iatom = 0; iatom < yCloud->nop - 1; iatom++) {
147  if (yCloud->pts[iatom].type != typeI) {
148  continue;
149  }
150  // Loop through the other atoms
151  for (int jatom = iatom + 1; jatom < yCloud->nop; jatom++) {
152  if (yCloud->pts[jatom].type != typeI) {
153  continue;
154  }
155  // If the distance is greater than rcutoff, continue
156  r_ij = gen::periodicDist(yCloud, iatom, jatom);
157  if (r_ij > rcutoff) {
158  continue;
159  }
160 
161  // Get the atom IDs for iatom and jatom
162  auto gotI = std::find_if(
163  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
164  [&iatom](const std::pair<int, int> &p) { return p.second == iatom; });
165  if (gotI == yCloud->idIndexMap.end()) {
166  std::cerr << "Something is wrong with your idIndexMap!\n";
167  return nList;
168  } else {
169  iatomIndex = gotI->first;
170  } // End of finding the atom ID for iatom
171  // Find the atom ID of jatom
172  auto gotJ = std::find_if(
173  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
174  [&jatom](const std::pair<int, int> &p) { return p.second == jatom; });
175  if (gotJ == yCloud->idIndexMap.end()) {
176  std::cerr << "Something is wrong with your idIndexMap!\n";
177  return nList;
178  } else {
179  jatomIndex = gotJ->first;
180  } // End of finding the atom ID for jatom
181  // Update the neighbour indices with atom IDs for iatom and jatom both
182  // (full list)
183  nList[iatom].push_back(jatomIndex);
184  nList[jatom].push_back(iatomIndex);
185 
186  } // End of loop through jatom
187  } // End of loop for iatom
188 
189  return nList;
190 }
191 
203 nneigh::halfNeighList(double rcutoff,
204  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
205  int typeI) {
207  nList; // Vector of vectors of the neighbour list
208  double r_ij; // Distance between iatom and jatom
209  int iatomIndex; // Atomic ID of the atom with index iatom
210  int jatomIndex; // Atomic ID of the atom with index jatom
211  int indexYay;
212  std::vector<int> tempListIatom;
213 
214  // Initialize and fill the first element with the current atom ID whose
215  // neighbour list will be filled
216  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
217  // Find the atom ID (key) given the index or iatom (value)
218  auto itr = std::find_if(
219  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
220  [&iatom](const std::pair<int, int> &p) { return p.second == iatom; });
221  // If found:
222  if (itr == yCloud->idIndexMap.end()) {
223  std::cerr << "Something is wrong with your idIndexMap!\n";
224  continue;
225  } else {
226  iatomIndex = itr->first;
227  } // End of finding the atom ID to fill as the first element in the
228  // neighbour list
229  nList.push_back(std::vector<int>()); // Empty vector for the index iatom
230  // Fill the first element with the atom ID of iatom itself
231  nList[iatom].push_back(iatomIndex);
232  } // end of init
233 
234  // Loop through every iatom and find nearest neighbours within rcutoff
235  for (int iatom = 0; iatom < yCloud->nop - 1; iatom++) {
236  if (yCloud->pts[iatom].type != typeI) {
237  continue;
238  }
239  // Loop through the other atoms
240  for (int jatom = iatom + 1; jatom < yCloud->nop; jatom++) {
241  if (yCloud->pts[jatom].type != typeI) {
242  continue;
243  }
244  // If the distance is greater than rcutoff, continue
245  r_ij = gen::periodicDist(yCloud, iatom, jatom);
246  if (r_ij > rcutoff) {
247  continue;
248  }
249 
250  // Get the atom IDs for iatom and jatom
251  auto gotI = std::find_if(
252  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
253  [&iatom](const std::pair<int, int> &p) { return p.second == iatom; });
254  if (gotI == yCloud->idIndexMap.end()) {
255  std::cerr << "Something is wrong with your idIndexMap!\n";
256  return nList;
257  } else {
258  iatomIndex = gotI->first;
259  } // End of finding the atom ID for iatom
260  // Find the atom ID of jatom
261  auto gotJ = std::find_if(
262  yCloud->idIndexMap.begin(), yCloud->idIndexMap.end(),
263  [&jatom](const std::pair<int, int> &p) { return p.second == jatom; });
264  if (gotJ == yCloud->idIndexMap.end()) {
265  std::cerr << "Something is wrong with your idIndexMap!\n";
266  return nList;
267  } else {
268  jatomIndex = gotJ->first;
269  } // End of finding the atom ID for jatom
270  // Update the neighbour indices with atom IDs for iatom and jatom both
271  // (full list)
272  nList[iatom].push_back(jatomIndex);
273 
274  } // End of loop through jatom
275  } // End of loop for iatom
276 
277  return nList;
278 }
279 
291  molSys::PointCloud<molSys::Point<double>, double> *yCloud, double cutoff) {
292  //
294  double r_ij; // Distance between iatom and jatom
295  std::vector<int> tempListIatom;
296 
297  // Initialize and fill the first element with the current atom ID whose
298  // neighbour list will be filled
299  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
300  //
301  nList.push_back(std::vector<int>()); // Empty vector for the index iatom
302  // Fill the first element with the atom ID of iatom itself
303  nList[iatom].push_back(iatom);
304  } // end of init
305  // -------------------------------------------------------
306  // Loop through every iatom and find nearest neighbours within rcutoff
307  for (int iatom = 0; iatom < yCloud->nop - 1; iatom++) {
308  // Loop through the other atoms
309  for (int jatom = iatom + 1; jatom < yCloud->nop; jatom++) {
310  // If the distance is greater than rcutoff, continue
311  r_ij = gen::periodicDist(yCloud, iatom, jatom);
312  if (r_ij > cutoff) {
313  continue;
314  }
315 
316  // Update the neighbour indices with atom IDs for iatom and jatom both
317  // (full list)
318  nList[iatom].push_back(jatom);
319  nList[jatom].push_back(iatom);
320 
321  } // End of loop through jatom
322  } // End of loop for iatom
323 
324  return nList;
325 } // end of function
326 
338  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
339  std::vector<std::vector<int>> nList) {
340  //
341  std::vector<std::vector<int>> indexNlist; // Desired neighbour list of indices
342  int iatomID, jatomID; // Atom IDs
343  int iatomIndex, jatomIndex; // Indices of iatom and jatom
344  int nnumNeighbours; // Number of nearest neighbours
345 
346  // Loop through every atom whose neighbours are contained in the neighbour
347  // list
348  for (int iatom = 0; iatom < nList.size(); iatom++) {
349  iatomID = nList[iatom][0]; // Atom ID
350  // Get the index of iatom
351  auto gotI = yCloud->idIndexMap.find(iatomID);
352  if (gotI != yCloud->idIndexMap.end()) {
353  iatomIndex = gotI->second;
354  } // found iatomIndex
355  //
356  nnumNeighbours = nList[iatomIndex].size() - 1;
357  // Update the new neighbour list
358  indexNlist.push_back(
359  std::vector<int>()); // Empty vector for the index iatom
360  // Fill the first element with the atom ID of iatom itself
361  indexNlist[iatom].push_back(iatomIndex);
362  //
363  // Loop through the neighbours of iatom
364  for (int jatom = 1; jatom <= nnumNeighbours; jatom++) {
365  jatomID = nList[iatomIndex][jatom]; // Atom ID of neighbour
366  //
367  // Get the index of the j^th atom
368  auto gotJ = yCloud->idIndexMap.find(jatomID);
369  if (gotJ != yCloud->idIndexMap.end()) {
370  jatomIndex = gotJ->second;
371  } // found jatomIndex
372  // Add to the neighbour list
373  indexNlist[iatom].push_back(jatomIndex);
374  } // end of loop through neighbours
375  } // end of loop through every atom
376 
377  // Return the new neighbour list
378  return indexNlist;
379 }
380 
388  //
389  std::vector<std::vector<int>> tempEmpty;
390 
391  nList.swap(tempEmpty);
392 
393  return 0;
394 }
T find_if(T... args)
double periodicDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Inline generic function for obtaining the unwrapped periodic distance between two particles,...
Definition: generic.hpp:104
std::vector< std::vector< int > > neighList(double rcutoff, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI, int typeJ)
All these functions use atom IDs and not indices.
Definition: neighbours.cpp:26
std::vector< std::vector< int > > halfNeighList(double rcutoff, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI=1)
Definition: neighbours.cpp:203
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
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
Header file for neighbour list generation.
T push_back(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
T swap(T... args)