bond.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 // Internal Libraries
12 #include <bond.hpp>
13 #include <generic.hpp>
14 
23  molSys::PointCloud<molSys::Point<double>, double> *yCloud) {
24  //
25  std::vector<std::vector<int>> bonds; // Output vector of vectors
26  std::vector<int> currentBond; // Vector for the current bond
27  int first, second; // Indices for the bonds
28  int iatom, jatom; // Elements of the bond
29  int iatomID, jatomID; // Atom IDs of the elements which are bonded
30 
31  // Error handling
32  if (nList.size() == 0) {
33  // There is some problem!
34  std::cerr << "There are no bonds in the system!\n";
35  return bonds;
36  }
37 
38  // Form of the bonds vector of vectors:
39  // 214 272
40  // 1 2
41  // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
42  // To avoid repitition, discard bonds whose first value is greater than the
43  // second value
44 
45  // Traverse the neighbour list
46 
47  // Loop through every atom in the neighbour list by index
48  for (int i = 0; i < nList.size(); i++) {
49  iatom = nList[i][0]; // Index of the i^th atom
50  // Get the neighbours of iatom
51  for (int j = 1; j < nList[i].size(); j++) {
52  //
53  jatom = nList[iatom][j]; // Index of the neighbour
54  // To avoid duplicates, skip all bonds such
55  // that iatom>jatom
56  if (iatom > jatom) {
57  continue;
58  } // Skip to avoid duplicates
59 
60  // Clear the current bond vector
61  currentBond.clear();
62  // Fill the current bond vector with ATOM IDs
63  iatomID = yCloud->pts[iatom].atomID;
64  jatomID = yCloud->pts[jatom].atomID;
65  currentBond.push_back(iatomID);
66  currentBond.push_back(jatomID);
67  // Add to the bonds vector of vectors
68  bonds.push_back(currentBond);
69  } // end of loop the neighbour list
70  // The last pair is with the last element and the first element
71  // Fill currentBond and update bonds
72 
73  } // end of loop through rings
74 
75  return bonds;
76 }
77 
91  std::vector<cage::iceType> atomTypes) {
92  //
93  std::vector<std::vector<int>> bonds; // Output vector of vectors
94  std::vector<int> currentBond; // Vector for the current bond
95  int first, second; // Indices for the bonds
96  int iatom, jatom; // Elements of the bond
97  int iatomID, jatomID; // Atom IDs of the elements which are bonded
98 
99  // Error handling
100  if (nList.size() == 0) {
101  // There is some problem!
102  std::cerr << "There are no bonds in the system!\n";
103  return bonds;
104  }
105 
106  // Form of the bonds vector of vectors:
107  // 214 272
108  // 1 2
109  // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
110  // To avoid repitition, discard bonds whose first value is greater than the
111  // second value
112 
113  // Traverse the neighbour list
114 
115  // Loop through every atom in the neighbour list by index
116  for (int i = 0; i < nList.size(); i++) {
117  iatom = nList[i][0]; // Index of the i^th atom
118  // Skip for dummy atoms
119  if (atomTypes[iatom] == cage::dummy) {
120  continue;
121  } // Skip for dummy atoms
122  // Get the neighbours of iatom
123  for (int j = 1; j < nList[i].size(); j++) {
124  //
125  jatom = nList[iatom][j]; // Index of the neighbour
126  // Skip for dummy atoms
127  if (atomTypes[jatom] == cage::dummy) {
128  continue;
129  } // Skip for dummy atoms
130  // To avoid duplicates, skip all bonds such
131  // that iatom>jatom
132  if (iatom > jatom) {
133  continue;
134  } // Skip to avoid duplicates
135 
136  // Clear the current bond vector
137  currentBond.clear();
138  // Fill the current bond vector with ATOM IDs
139  iatomID = yCloud->pts[iatom].atomID;
140  jatomID = yCloud->pts[jatom].atomID;
141  currentBond.push_back(iatomID);
142  currentBond.push_back(jatomID);
143  // Add to the bonds vector of vectors
144  bonds.push_back(currentBond);
145  } // end of loop the neighbour list
146  // The last pair is with the last element and the first element
147  // Fill currentBond and update bonds
148 
149  } // end of loop through rings
150 
151  return bonds;
152 }
153 
154 
171  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
172  std::vector<std::vector<int>> nList, int targetFrame,
173  int Htype) {
174  //
176  hBondNet; // Output vector of vectors containing the HBN
178  hCloud; // point Cloud for the hydrogen atoms
180  molIDlist; // Vector of vectors; first element is the molID, and the next
181  // two elements are the hydrogen atom indices
183  idMolIDmap; // Unordered map with atom IDs of oxygens as the keys and the
184  // molecular IDs as the values
185  std::vector<int> currentBondList; // Current bond list for atom
186  int nnumNeighbours; // Number of nearest neighbours for the current atom
187  int iatomID, jatomID; // Atom IDs
188  int iatomIndex, jatomIndex; // Atomic indices of oxygen atoms
189  int hAtomIndex; // Atom index of hydrogen
190  int listIndex; // Index in molIDlist corresponding to a particular molecular
191  // ID
192  int jOxyMolID; // Molecular ID of the jatom oxygen atom
193  double hBondLen; // Length of O-H (between the donor O and acceptor H)
194  double distCutoff = 2.42; // Distance cutoff of O-H, hard-coded
195  double angleCutoff = 30; // Angle cutoff in degrees
196  std::vector<double> ooVec; // Array for the O--O vector
197  std::vector<double> ohVec; // Array for the O-H vector
198 
199  // --------------------
200  // Get all the hydrogen atoms in the frame (no slice)
201  hCloud = sinp::readLammpsTrjreduced(filename, targetFrame, &hCloud, Htype);
202 
203  // Get the unordered map of the oxygen atom IDs (keys) and the molecular IDs
204  // (values)
205  idMolIDmap = molSys::createIDMolIDmap(yCloud);
206 
207  // Get a vector of vectors with the molID in the first column, and the
208  // hydrogen atom indices (not ID) in each row
209  molIDlist = molSys::hAtomMolList(&hCloud, yCloud);
210 
211  // Initialize the vector of vectors hBondNet
212  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
213  hBondNet.push_back(std::vector<int>()); // Empty vector for the index iatom
214  // Fill the first element with the atom ID of iatom itself
215  hBondNet[iatom].push_back(yCloud->pts[iatom].atomID);
216  } // end of init of hBondNet
217 
218  // Loop through the neighbour list
219  for (int iatom = 0; iatom < nList.size(); iatom++) {
220  currentBondList.clear(); // Clear the current bond vector
221  iatomID = nList[iatom][0]; // atom ID corresponding to iatom
222  nnumNeighbours = nList[iatom].size() - 1; // No. of nearest neighbours
223  iatomIndex = iatom; // Atomic index
224  //
225  // Loop through the nearest neighbours
226  for (int j = 1; j <= nnumNeighbours; j++) {
227  jatomID = nList[iatom][j]; // Atom ID of the nearest neighbour
228  // Get the hydrogen atom indices corresponding to the molID of jatomID
229  // Find jOxyMolID
230  auto it = idMolIDmap.find(jatomID);
231  if (it != idMolIDmap.end()) {
232  jOxyMolID = it->second;
233  } // found molecular ID of jatom oxygen atom
234  else {
235  continue;
236  } // not found
237 
238  // Find the index inside the molIDlist corresponding to the molecular ID
239  // to look for
240  listIndex = molSys::searchMolList(molIDlist, jOxyMolID);
241 
242  // Get the atom index of the oxygen atom jatom corresponding jatomID
243  auto gotJ = yCloud->idIndexMap.find(jatomID);
244  if (gotJ != yCloud->idIndexMap.end()) {
245  jatomIndex = gotJ->second;
246  } // found atom index of jatomID
247  else {
248  std::cerr << "Something is wrong with the map.\n";
249  continue;
250  } // index not found
251 
252  // Loop through the hydrogen atoms connected to jatom oxygen atom
253  for (int k = 1; k <= 2; k++) {
254  hAtomIndex = molIDlist[listIndex][k];
255  // --------
256  // Condition One: The O-H length (between the donor hydrogen atom and
257  // the acceptor oxygen atom) should be less than 2.42 Angstrom
258  // (hard-coded)
259  hBondLen =
260  bond::getHbondDistanceOH(yCloud, &hCloud, iatomIndex, hAtomIndex);
261 
262  // If O-H distance is greater than or equal to 2.42 then it is not a
263  // hydrogen bond
264  if (hBondLen >= distCutoff) {
265  continue;
266  } // not a hydrogen bond
267  // --------
268  // Condition Two: The angle between the O-H and O--O vectors is less
269  // than 30 degrees (hard-coded)
270  //
271  ooVec.clear();
272  ohVec.clear();
273  // Get the O--O and O-H vectors
274  // O--O
275  ooVec.push_back(yCloud->pts[iatomIndex].x -
276  yCloud->pts[jatomIndex].x); // dx
277  ooVec.push_back(yCloud->pts[iatomIndex].y -
278  yCloud->pts[jatomIndex].y); // dy
279  ooVec.push_back(yCloud->pts[iatomIndex].z -
280  yCloud->pts[jatomIndex].z); // dz
281  // O-H
282  ohVec.push_back(yCloud->pts[iatomIndex].x -
283  hCloud.pts[hAtomIndex].x); // dx
284  ohVec.push_back(yCloud->pts[iatomIndex].y -
285  hCloud.pts[hAtomIndex].y); // dy
286  ohVec.push_back(yCloud->pts[iatomIndex].z -
287  hCloud.pts[hAtomIndex].z); // dz
288  // Apply PBCs
289  for (int l = 0; l < 3; l++) {
290  ooVec[l] -= yCloud->box[l] * round(ooVec[l] / yCloud->box[l]);
291  ohVec[l] -= yCloud->box[l] * round(ohVec[l] / yCloud->box[l]);
292  } // end of applying PBCs to the O-H and O--O vectors
293  //
294  // Get the angle between the O--O and O-H vectors
295  double eigenAngle = gen::eigenVecAngle(ooVec, ohVec);
296  double eigenAngleDeg = gen::radDeg(eigenAngle);
297 
298  //
299  // A hydrogen bond is formed if the angle is less than 30 degrees
300  if (eigenAngleDeg > angleCutoff) {
301  continue;
302  } // not a hydrogen bond
303 
304  // If you have reached this point, then O and H and indeed
305  // hydrogen-bonded. This means that jatomID should be saved in the new
306  // currentBond
307  hBondNet[iatomIndex].push_back(jatomID);
308  hBondNet[jatomIndex].push_back(iatomID);
309  break; // No need to test the other hydrogen atom if the first has
310  // been tested
311  } // end of loop through hydrogen atoms
312 
313  } // end of loop through the nearest neighbours
314 
315  } // end of loop through the neighbour list
316 
317  // --------------------
318 
319  // Erase all temporary stuff
320  hCloud = molSys::clearPointCloud(&hCloud);
321 
322  return hBondNet;
323 }
324 
335  molSys::PointCloud<molSys::Point<double>, double> *oCloud,
336  molSys::PointCloud<molSys::Point<double>, double> *hCloud, int oAtomIndex,
337  int hAtomIndex) {
338  std::array<double, 3> dr; // relative distance in the X, Y, Z dimensions
339  double r2 = 0.0; // Bond length
340 
341  dr[0] = oCloud->pts[oAtomIndex].x - hCloud->pts[hAtomIndex].x;
342  dr[1] = oCloud->pts[oAtomIndex].y - hCloud->pts[hAtomIndex].y;
343  dr[2] = oCloud->pts[oAtomIndex].z - hCloud->pts[hAtomIndex].z;
344 
345  // Apply the PBCs and get the squared area
346  for (int k = 0; k < 3; k++) {
347  dr[k] -= oCloud->box[k] * round(dr[k] / oCloud->box[k]);
348  r2 += pow(dr[k], 2.0);
349  } // end of applying the PBCs and getting the squared area
350  return sqrt(r2);
351 }
352 
366  std::vector<cage::Cage> *cageList,
367  cage::cageType type, int *nRings) {
368  std::vector<std::vector<int>> bonds; // Output vector of vectors
369  std::vector<int> currentBond; // Vector for the current bond
370  int ringSize = rings[0].size();
371  int currentRing; // (vector) index of the current ring in a particular cage
372 
373  // Error handling
374  if (rings.size() == 0) {
375  // There is some problem!
376  std::cerr << "There are no rings in the system!\n";
377  return bonds;
378  }
379 
380  // Form of the bonds vector of vectors:
381  // 272 214
382  // 1 2
383  // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
384 
385  // Traverse the cageList vector of Cages
386 
387  *nRings = 0; // init
388 
389  // Loop through all the cages
390  for (int icage = 0; icage < (*cageList).size(); icage++) {
391  // Skip if the cage is of a different type
392  if ((*cageList)[icage].type != type) {
393  continue;
394  }
395  *nRings += (*cageList)[icage].rings.size(); // Add to the number of rings
396  //
397  // Now loop through a particular ring inside the i^th cage
398  for (int iring = 0; iring < (*cageList)[icage].rings.size(); iring++) {
399  currentRing = (*cageList)[icage].rings[iring]; // Current ring index
400  // Get the first atom of each pair inside currentRing
401  for (int k = 0; k < rings[currentRing].size() - 1; k++) {
402  // Clear the current bond vector
403  currentBond.clear();
404  // Fill the current bond vector
405  currentBond.push_back(rings[currentRing][k]);
406  currentBond.push_back(rings[currentRing][k + 1]);
407  std::sort(currentBond.begin(), currentBond.end());
408  // Add to the bonds vector of vectors
409  bonds.push_back(currentBond);
410  } // end of loop through ring elements, except the last one
411  // The last pair is with the last element and the first element
412  // Fill currentBond and update bonds
413  currentBond.clear();
414  currentBond.push_back(rings[currentRing][ringSize - 1]);
415  currentBond.push_back(rings[currentRing][0]);
416  std::sort(currentBond.begin(), currentBond.end());
417  bonds.push_back(currentBond);
418  } // end of loop through a particular ring
419  } // end of loop through cages
420 
421  // This may have duplicates, so the duplicates should be removed
422  std::sort(bonds.begin(), bonds.end());
423  bonds.erase(std::unique(bonds.begin(), bonds.end()), bonds.end());
424 
425  return bonds;
426 }
427 
438  reversedBond; // Vector for the current bond in reversed order
439  std::vector<bool> isBondFlag;
440  int temp = 0;
441 
442  std::sort(bonds.begin(), bonds.end());
443  bonds.erase(std::unique(bonds.begin(), bonds.end()), bonds.end());
444 
445  // // Temp uncommented
446  // // Resize the flag vector and init to true
447  // isBondFlag.resize(bonds.size(), true);
448  // // Form of the bonds vector of vectors:
449  // // 272 214
450  // // 1 2
451  // // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
452 
453  // // Traverse the bonds vector of vectors
454 
455  // // Loop through all the bonds
456  // for(int ibond=0; ibond<bonds.size(); ibond++){
457 
458  // // Skip if the bond has been flagged as false
459  // if(isBondFlag[ibond] == false){continue;}
460 
461  // reversedBond.clear();
462  // reversedBond.push_back( bonds[ibond][1] );
463  // reversedBond.push_back( bonds[ibond][0] );
464 
465  // // Compare with other bonds by looping through
466  // // the rest
467  // for(int jbond=0; jbond<bonds.size(); jbond++){
468  // // Skip if jbond has been flagged
469  // if(isBondFlag[jbond] == false){continue;}
470  // // Compare reversed bond and jbond if jbond has not been flagged yet
471  // if(reversedBond==bonds[jbond]){
472  // isBondFlag[jbond] = false;
473  // temp++;
474  // } // end of check of comparison
475  // } // end of looping through all possible bonds
476  // } // end of loop through all possible bonds
477  // // temp
478 
479  return bonds;
480 }
T begin(T... args)
File for bond-related analyses (hydrogen bonds, bonded atoms for data file write-outs etc....
T clear(T... args)
T end(T... args)
T erase(T... args)
T find(T... args)
File for containing generic or common functions.
std::vector< std::vector< int > > populateHbonds(std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int targetFrame, int Htype)
Definition: bond.cpp:170
std::vector< std::vector< int > > trimBonds(std::vector< std::vector< int >> bonds)
Remove duplicate bonds.
Definition: bond.cpp:436
std::vector< std::vector< int > > populateBonds(std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: bond.cpp:22
double getHbondDistanceOH(molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *hCloud, int oAtomIndex, int hAtomIndex)
Definition: bond.cpp:334
std::vector< std::vector< int > > createBondsFromCages(std::vector< std::vector< int >> rings, std::vector< cage::Cage > *cageList, cage::cageType type, int *nRings)
Definition: bond.cpp:365
cageType
Definition: cage.hpp:50
@ dummy
Definition: cage.hpp:71
double eigenVecAngle(std::vector< double > OO, std::vector< double > OH)
Eigen function for getting the angle (in radians) between the O–O and O-H vectors.
Definition: generic.cpp:150
double radDeg(double angle)
Definition: generic.hpp:57
std::vector< std::vector< int > > hAtomMolList(molSys::PointCloud< molSys::Point< double >, double > *hCloud, molSys::PointCloud< molSys::Point< double >, double > *oCloud)
Definition: mol_sys.cpp:62
std::unordered_map< int, int > createIDMolIDmap(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: mol_sys.cpp:40
int searchMolList(std::vector< std::vector< int >> molList, int molIDtoFind)
Definition: mol_sys.cpp:104
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
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrjreduced(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
T push_back(T... args)
T size(T... args)
T sort(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 unique(T... args)