Bond

Functions

std::vector< std::vector< int > > bond::populateHbonds (std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int targetFrame, int Htype)
 
std::vector< std::vector< int > > bond::populateHbondsWithInputClouds (molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *hCloud, std::vector< std::vector< int >> nList)
 
double bond::getHbondDistanceOH (molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *hCloud, int oAtomIndex, int hAtomIndex)
 
std::vector< std::vector< int > > bond::populateBonds (std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 
std::vector< std::vector< int > > bond::populateBonds (std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< cage::iceType > atomTypes)
 
std::vector< std::vector< int > > bond::createBondsFromCages (std::vector< std::vector< int >> rings, std::vector< cage::Cage > *cageList, cage::cageType type, int *nRings)
 
std::vector< std::vector< int > > bond::trimBonds (std::vector< std::vector< int >> bonds)
 Remove duplicate bonds. More...
 

Detailed Description

Function Documentation

◆ createBondsFromCages()

std::vector< std::vector< int > > bond::createBondsFromCages ( std::vector< std::vector< int >>  rings,
std::vector< cage::Cage > *  cageList,
cage::cageType  type,
int *  nRings 
)

Creates a vector of vectors containing bond connectivity information from the rings vector of vectors and cage information

Create a vector of vectors containing bond information (bonded atom IDs, not vector or array indices!) from the ring vector of vectors and cageList

Parameters
[in]ringsRow-ordered vector of vectors atom indices of ring information. Each row is a ring, containing the indices of the particles in that ring
[in]cageListA vector of cage::Cage containing a list of HCs or DDCs
[in]typeThe type of cage to get bonds for
[in,out]nRingsThe total number of rings for all the cages, for the particular cage type

Definition at line 529 of file bond.cpp.

531  {
532  std::vector<std::vector<int>> bonds; // Output vector of vectors
533  std::vector<int> currentBond; // Vector for the current bond
534  int ringSize = rings[0].size();
535  int currentRing; // (vector) index of the current ring in a particular cage
536 
537  // Error handling
538  if (rings.size() == 0) {
539  // There is some problem!
540  std::cerr << "There are no rings in the system!\n";
541  return bonds;
542  }
543 
544  // Form of the bonds vector of vectors:
545  // 272 214
546  // 1 2
547  // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
548 
549  // Traverse the cageList vector of Cages
550 
551  *nRings = 0; // init
552 
553  // Loop through all the cages
554  for (int icage = 0; icage < (*cageList).size(); icage++) {
555  // Skip if the cage is of a different type
556  if ((*cageList)[icage].type != type) {
557  continue;
558  }
559  *nRings += (*cageList)[icage].rings.size(); // Add to the number of rings
560  //
561  // Now loop through a particular ring inside the i^th cage
562  for (int iring = 0; iring < (*cageList)[icage].rings.size(); iring++) {
563  currentRing = (*cageList)[icage].rings[iring]; // Current ring index
564  // Get the first atom of each pair inside currentRing
565  for (int k = 0; k < rings[currentRing].size() - 1; k++) {
566  // Clear the current bond vector
567  currentBond.clear();
568  // Fill the current bond vector
569  currentBond.push_back(rings[currentRing][k]);
570  currentBond.push_back(rings[currentRing][k + 1]);
571  std::sort(currentBond.begin(), currentBond.end());
572  // Add to the bonds vector of vectors
573  bonds.push_back(currentBond);
574  } // end of loop through ring elements, except the last one
575  // The last pair is with the last element and the first element
576  // Fill currentBond and update bonds
577  currentBond.clear();
578  currentBond.push_back(rings[currentRing][ringSize - 1]);
579  currentBond.push_back(rings[currentRing][0]);
580  std::sort(currentBond.begin(), currentBond.end());
581  bonds.push_back(currentBond);
582  } // end of loop through a particular ring
583  } // end of loop through cages
584 
585  // This may have duplicates, so the duplicates should be removed
586  std::sort(bonds.begin(), bonds.end());
587  bonds.erase(std::unique(bonds.begin(), bonds.end()), bonds.end());
588 
589  return bonds;
590 }
T begin(T... args)
T clear(T... args)
T end(T... args)
T erase(T... args)
T push_back(T... args)
T size(T... args)
T sort(T... args)
T unique(T... args)

◆ getHbondDistanceOH()

double bond::getHbondDistanceOH ( molSys::PointCloud< molSys::Point< double >, double > *  oCloud,
molSys::PointCloud< molSys::Point< double >, double > *  hCloud,
int  oAtomIndex,
int  hAtomIndex 
)

Calculates the distance of the hydrogen bond between O and H (of different atoms), given the respective pointClouds and the indices to each atom

Calculates the bond length between a Hydrogen and Oxygen atom of two different atoms, given their respective pointClouds and the indices to each atom.

Parameters
[in]oCloudThe molSys::PointCloud for the oxygen atoms only
[in]hCloudThe molSys::PointCloud for the hydrogen atoms only
[in]oAtomIndexThe index (in the oCloud) of the oxygen atom
[in]hAtomIndexThe index (in the hCloud) of the hydrogen atom

Definition at line 498 of file bond.cpp.

501  {
502  std::array<double, 3> dr; // relative distance in the X, Y, Z dimensions
503  double r2 = 0.0; // Bond length
504 
505  dr[0] = oCloud->pts[oAtomIndex].x - hCloud->pts[hAtomIndex].x;
506  dr[1] = oCloud->pts[oAtomIndex].y - hCloud->pts[hAtomIndex].y;
507  dr[2] = oCloud->pts[oAtomIndex].z - hCloud->pts[hAtomIndex].z;
508 
509  // Apply the PBCs and get the squared area
510  for (int k = 0; k < 3; k++) {
511  dr[k] -= oCloud->box[k] * round(dr[k] / oCloud->box[k]);
512  r2 += pow(dr[k], 2.0);
513  } // end of applying the PBCs and getting the squared area
514  return sqrt(r2);
515 }
std::vector< S > pts
Definition: mol_sys.hpp:171
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:174
T pow(T... args)
T round(T... args)
T sqrt(T... args)

◆ populateBonds() [1/2]

std::vector< std::vector< int > > bond::populateBonds ( std::vector< std::vector< int >>  nList,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud 
)

Create a vector of vectors containing bond connectivity information. May contain duplicates! Gets the bond information from the vector of vectors containing the neighbour list by index

Create a vector of vectors containing bond information (outputs bonded atom IDs, not indices!) from the neighbour list vector of vectors (which contains atom INDICES). Moreover, the first element corresponds to the atom whose neighbours have been found.

Definition at line 26 of file bond.cpp.

27  {
28  //
29  std::vector<std::vector<int>> bonds; // Output vector of vectors
30  std::vector<int> currentBond; // Vector for the current bond
31  int first, second; // Indices for the bonds
32  int iatom, jatom; // Elements of the bond
33  int iatomID, jatomID; // Atom IDs of the elements which are bonded
34 
35  // Error handling
36  if (nList.size() == 0) {
37  // There is some problem!
38  std::cerr << "There are no bonds in the system!\n";
39  return bonds;
40  }
41 
42  // Form of the bonds vector of vectors:
43  // 214 272
44  // 1 2
45  // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
46  // To avoid repitition, discard bonds whose first value is greater than the
47  // second value
48 
49  // Traverse the neighbour list
50 
51  // Loop through every atom in the neighbour list by index
52  for (int i = 0; i < nList.size(); i++) {
53  iatom = nList[i][0]; // Index of the i^th atom
54  // Get the neighbours of iatom
55  for (int j = 1; j < nList[i].size(); j++) {
56  //
57  jatom = nList[iatom][j]; // Index of the neighbour
58  // To avoid duplicates, skip all bonds such
59  // that iatom>jatom
60  if (iatom > jatom) {
61  continue;
62  } // Skip to avoid duplicates
63 
64  // Clear the current bond vector
65  currentBond.clear();
66  // Fill the current bond vector with ATOM IDs
67  iatomID = yCloud->pts[iatom].atomID;
68  jatomID = yCloud->pts[jatom].atomID;
69  currentBond.push_back(iatomID);
70  currentBond.push_back(jatomID);
71  // Add to the bonds vector of vectors
72  bonds.push_back(currentBond);
73  } // end of loop the neighbour list
74  // The last pair is with the last element and the first element
75  // Fill currentBond and update bonds
76 
77  } // end of loop through rings
78 
79  return bonds;
80 }

◆ populateBonds() [2/2]

std::vector< std::vector< int > > bond::populateBonds ( std::vector< std::vector< int >>  nList,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< cage::iceType atomTypes 
)

Create a vector of vectors containing bond connectivity information Gets the bond information from the vector of vectors containing the neighbour list by index. Bonds between dummy atoms are not filled.

Create a vector of vectors containing bond information (outputs bonded atom IDs, not indices!) from the neighbour list vector of vectors (which contains atom INDICES). Bonds between dummy atoms, and between dummy and ice atoms are not added. Moreover, the first element corresponds to the atom whose neighbours have been found.

Parameters
[in]nListRow-ordered neighbour list by ID
[in]yCloudThe input molSys::PointCloud
[in]atomTypesContains an atom type for each particle in yCloud

Definition at line 93 of file bond.cpp.

95  {
96  //
97  std::vector<std::vector<int>> bonds; // Output vector of vectors
98  std::vector<int> currentBond; // Vector for the current bond
99  int first, second; // Indices for the bonds
100  int iatom, jatom; // Elements of the bond
101  int iatomID, jatomID; // Atom IDs of the elements which are bonded
102 
103  // Error handling
104  if (nList.size() == 0) {
105  // There is some problem!
106  std::cerr << "There are no bonds in the system!\n";
107  return bonds;
108  }
109 
110  // Form of the bonds vector of vectors:
111  // 214 272
112  // 1 2
113  // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
114  // To avoid repitition, discard bonds whose first value is greater than the
115  // second value
116 
117  // Traverse the neighbour list
118 
119  // Loop through every atom in the neighbour list by index
120  for (int i = 0; i < nList.size(); i++) {
121  iatom = nList[i][0]; // Index of the i^th atom
122  // Skip for dummy atoms
123  if (atomTypes[iatom] == cage::iceType::dummy) {
124  continue;
125  } // Skip for dummy atoms
126  // Get the neighbours of iatom
127  for (int j = 1; j < nList[i].size(); j++) {
128  //
129  jatom = nList[iatom][j]; // Index of the neighbour
130  // Skip for dummy atoms
131  if (atomTypes[jatom] == cage::iceType::dummy) {
132  continue;
133  } // Skip for dummy atoms
134  // To avoid duplicates, skip all bonds such
135  // that iatom>jatom
136  if (iatom > jatom) {
137  continue;
138  } // Skip to avoid duplicates
139 
140  // Clear the current bond vector
141  currentBond.clear();
142  // Fill the current bond vector with ATOM IDs
143  iatomID = yCloud->pts[iatom].atomID;
144  jatomID = yCloud->pts[jatom].atomID;
145  currentBond.push_back(iatomID);
146  currentBond.push_back(jatomID);
147  // Add to the bonds vector of vectors
148  bonds.push_back(currentBond);
149  } // end of loop the neighbour list
150  // The last pair is with the last element and the first element
151  // Fill currentBond and update bonds
152 
153  } // end of loop through rings
154 
155  return bonds;
156 }

◆ populateHbonds()

std::vector< std::vector< int > > bond::populateHbonds ( std::string  filename,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< std::vector< int >>  nList,
int  targetFrame,
int  Htype 
)

Create a vector of vectors (similar to the neighbour list conventions) containing the hydrogen bond connectivity information. Decides the existence of the hydrogen bond depending on the O–O and O–H vectors from the neighbour list already constructed

Create a vector of vectors (similar to the neighbour list conventions). The output vector of vectors is row-ordered. The first element is the atom ID of the particle for which the neighbours are enumerated (the central atom), followed by the central atom's neighbour's IDs (not indices). Decides the existence of the hydrogen bond depending on the O–O and O–H vectors from the neighbour list (by ID) already constructed.

Parameters
[in]filenameFilename of the trajectory, with the hydrogen and oxygen coordinates
[in]yCloudThe input molSys::PointCloud for the oxygen atoms only
[in]nListRow-ordered neighbour list by atom ID
[in]targetFrameThe target or current frame number (starts from 1) and is not the timestep value
[in]HtypeThe type ID of the hydrogen atoms

Definition at line 174 of file bond.cpp.

177  {
178  //
180  hBondNet; // Output vector of vectors containing the HBN
182  hCloud; // point Cloud for the hydrogen atoms
184  molIDlist; // Vector of vectors; first element is the molID, and the next
185  // two elements are the hydrogen atom indices
187  idMolIDmap; // Unordered map with atom IDs of oxygens as the keys and the
188  // molecular IDs as the values
189  std::vector<int> currentBondList; // Current bond list for atom
190  int nnumNeighbours; // Number of nearest neighbours for the current atom
191  int iatomID, jatomID; // Atom IDs
192  int iatomIndex, jatomIndex; // Atomic indices of oxygen atoms
193  int hAtomIndex; // Atom index of hydrogen
194  int listIndex; // Index in molIDlist corresponding to a particular molecular
195  // ID
196  int jOxyMolID; // Molecular ID of the jatom oxygen atom
197  double hBondLen; // Length of O-H (between the donor O and acceptor H)
198  double distCutoff = 2.42; // Distance cutoff of O-H, hard-coded
199  double angleCutoff = 30; // Angle cutoff in degrees
200  std::vector<double> ooVec; // Array for the O--O vector
201  std::vector<double> ohVec; // Array for the O-H vector
202 
203  // --------------------
204  // Get all the hydrogen atoms in the frame (no slice)
205  hCloud = sinp::readLammpsTrjreduced(filename, targetFrame, &hCloud, Htype);
206 
207  // Get the unordered map of the oxygen atom IDs (keys) and the molecular IDs
208  // (values)
209  idMolIDmap = molSys::createIDMolIDmap(yCloud);
210 
211  // Get a vector of vectors with the molID in the first column, and the
212  // hydrogen atom indices (not ID) in each row
213  molIDlist = molSys::hAtomMolList(&hCloud, yCloud);
214 
215  // Initialize the vector of vectors hBondNet
216  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
217  hBondNet.push_back(std::vector<int>()); // Empty vector for the index iatom
218  // Fill the first element with the atom ID of iatom itself
219  hBondNet[iatom].push_back(yCloud->pts[iatom].atomID);
220  } // end of init of hBondNet
221 
222  // Loop through the neighbour list
223  for (int iatom = 0; iatom < nList.size(); iatom++) {
224  currentBondList.clear(); // Clear the current bond vector
225  iatomID = nList[iatom][0]; // atom ID corresponding to iatom
226  nnumNeighbours = nList[iatom].size() - 1; // No. of nearest neighbours
227  iatomIndex = iatom; // Atomic index
228  //
229  // Loop through the nearest neighbours
230  for (int j = 1; j <= nnumNeighbours; j++) {
231  jatomID = nList[iatom][j]; // Atom ID of the nearest neighbour
232  // Get the hydrogen atom indices corresponding to the molID of jatomID
233  // Find jOxyMolID
234  auto it = idMolIDmap.find(jatomID);
235  if (it != idMolIDmap.end()) {
236  jOxyMolID = it->second;
237  } // found molecular ID of jatom oxygen atom
238  else {
239  continue;
240  } // not found
241 
242  // Find the index inside the molIDlist corresponding to the molecular ID
243  // to look for
244  listIndex = molSys::searchMolList(molIDlist, jOxyMolID);
245 
246  // Get the atom index of the oxygen atom jatom corresponding jatomID
247  auto gotJ = yCloud->idIndexMap.find(jatomID);
248  if (gotJ != yCloud->idIndexMap.end()) {
249  jatomIndex = gotJ->second;
250  } // found atom index of jatomID
251  else {
252  std::cerr << "Something is wrong with the map.\n";
253  continue;
254  } // index not found
255 
256  // Loop through the hydrogen atoms connected to jatom oxygen atom
257  for (int k = 1; k <= 2; k++) {
258  hAtomIndex = molIDlist[listIndex][k];
259  // --------
260  // Condition One: The O-H length (between the donor hydrogen atom and
261  // the acceptor oxygen atom) should be less than 2.42 Angstrom
262  // (hard-coded)
263  hBondLen =
264  bond::getHbondDistanceOH(yCloud, &hCloud, iatomIndex, hAtomIndex);
265 
266  // If O-H distance is greater than or equal to 2.42 then it is not a
267  // hydrogen bond
268  if (hBondLen >= distCutoff) {
269  continue;
270  } // not a hydrogen bond
271  // --------
272  // Condition Two: The angle between the O-H and O--O vectors is less
273  // than 30 degrees (hard-coded)
274  //
275  ooVec.clear();
276  ohVec.clear();
277  // Get the O--O and O-H vectors
278  // O--O
279  ooVec.push_back(yCloud->pts[iatomIndex].x -
280  yCloud->pts[jatomIndex].x); // dx
281  ooVec.push_back(yCloud->pts[iatomIndex].y -
282  yCloud->pts[jatomIndex].y); // dy
283  ooVec.push_back(yCloud->pts[iatomIndex].z -
284  yCloud->pts[jatomIndex].z); // dz
285  // O-H
286  ohVec.push_back(yCloud->pts[iatomIndex].x -
287  hCloud.pts[hAtomIndex].x); // dx
288  ohVec.push_back(yCloud->pts[iatomIndex].y -
289  hCloud.pts[hAtomIndex].y); // dy
290  ohVec.push_back(yCloud->pts[iatomIndex].z -
291  hCloud.pts[hAtomIndex].z); // dz
292  // Apply PBCs
293  for (int l = 0; l < 3; l++) {
294  ooVec[l] -= yCloud->box[l] * round(ooVec[l] / yCloud->box[l]);
295  ohVec[l] -= yCloud->box[l] * round(ohVec[l] / yCloud->box[l]);
296  } // end of applying PBCs to the O-H and O--O vectors
297  //
298  // Get the angle between the O--O and O-H vectors
299  double eigenAngle = gen::eigenVecAngle(ooVec, ohVec);
300  double eigenAngleDeg = gen::radDeg(eigenAngle);
301 
302  //
303  // A hydrogen bond is formed if the angle is less than 30 degrees
304  if (eigenAngleDeg > angleCutoff) {
305  continue;
306  } // not a hydrogen bond
307 
308  // If you have reached this point, then O and H and indeed
309  // hydrogen-bonded. This means that jatomID should be saved in the new
310  // currentBond
311  hBondNet[iatomIndex].push_back(jatomID);
312  hBondNet[jatomIndex].push_back(iatomID);
313  break; // No need to test the other hydrogen atom if the first has
314  // been tested
315  } // end of loop through hydrogen atoms
316 
317  } // end of loop through the nearest neighbours
318 
319  } // end of loop through the neighbour list
320 
321  // --------------------
322 
323  // Erase all temporary stuff
324  hCloud = molSys::clearPointCloud(&hCloud);
325 
326  return hBondNet;
327 }
T find(T... args)
double getHbondDistanceOH(molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *hCloud, int oAtomIndex, int hAtomIndex)
Definition: bond.cpp:498
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:155
double radDeg(double angle)
Definition: generic.hpp:61
int nop
Current frame number.
Definition: mol_sys.hpp:173
std::vector< std::vector< int > > hAtomMolList(molSys::PointCloud< molSys::Point< double >, double > *hCloud, molSys::PointCloud< molSys::Point< double >, double > *oCloud)
Definition: mol_sys.cpp:89
std::unordered_map< int, int > createIDMolIDmap(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: mol_sys.cpp:44
int searchMolList(std::vector< std::vector< int >> molList, int molIDtoFind)
Definition: mol_sys.cpp:131
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:176
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
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})
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:170

◆ populateHbondsWithInputClouds()

std::vector< std::vector< int > > bond::populateHbondsWithInputClouds ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
molSys::PointCloud< molSys::Point< double >, double > *  hCloud,
std::vector< std::vector< int >>  nList 
)

Create a vector of vectors (similar to the neighbour list conventions) containing the hydrogen bond connectivity information. Decides the existence of the hydrogen bond depending on the O–O and O–H vectors from the neighbour list already constructed, taking a PointCloud for the H atoms as input

Create a vector of vectors (similar to the neighbour list conventions). The output vector of vectors is row-ordered. The first element is the atom ID of the particle for which the neighbours are enumerated (the central atom), followed by the central atom's neighbour's IDs (not indices). Decides the existence of the hydrogen bond depending on the O–O and O–H vectors from the neighbour list (by ID) already constructed. The only difference between this function and the similar populateHbonds function is that this takes in the H atom PointCloud as input

Parameters
[in]yCloudThe input molSys::PointCloud for the hydrogen atoms only, for the entire system (regardless of whether there is a slice or not)
[in]hCloudThe input molSys::PointCloud for the oxygen atoms only
[in]nListRow-ordered neighbour list by atom ID

Definition at line 343 of file bond.cpp.

345  {
346  //
348  hBondNet; // Output vector of vectors containing the HBN
350  molIDlist; // Vector of vectors; first element is the molID, and the next
351  // two elements are the hydrogen atom indices
353  idMolIDmap; // Unordered map with atom IDs of oxygens as the keys and the
354  // molecular IDs as the values
355  std::vector<int> currentBondList; // Current bond list for atom
356  int nnumNeighbours; // Number of nearest neighbours for the current atom
357  int iatomID, jatomID; // Atom IDs
358  int iatomIndex, jatomIndex; // Atomic indices of oxygen atoms
359  int hAtomIndex; // Atom index of hydrogen
360  int listIndex; // Index in molIDlist corresponding to a particular molecular
361  // ID
362  int jOxyMolID; // Molecular ID of the jatom oxygen atom
363  double hBondLen; // Length of O-H (between the donor O and acceptor H)
364  double distCutoff = 2.42; // Distance cutoff of O-H, hard-coded
365  double angleCutoff = 30; // Angle cutoff in degrees
366  std::vector<double> ooVec; // Array for the O--O vector
367  std::vector<double> ohVec; // Array for the O-H vector
368 
369  // --------------------
370  // Get the unordered map of the oxygen atom IDs (keys) and the molecular IDs
371  // (values)
372  idMolIDmap = molSys::createIDMolIDmap(yCloud);
373 
374  // Get a vector of vectors with the molID in the first column, and the
375  // hydrogen atom indices (not ID) in each row
376  molIDlist = molSys::hAtomMolList(hCloud, yCloud);
377 
378  // Initialize the vector of vectors hBondNet
379  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
380  hBondNet.push_back(std::vector<int>()); // Empty vector for the index iatom
381  // Fill the first element with the atom ID of iatom itself
382  hBondNet[iatom].push_back(yCloud->pts[iatom].atomID);
383  } // end of init of hBondNet
384 
385  // Loop through the neighbour list
386  for (int iatom = 0; iatom < nList.size(); iatom++) {
387  currentBondList.clear(); // Clear the current bond vector
388  iatomID = nList[iatom][0]; // atom ID corresponding to iatom
389  nnumNeighbours = nList[iatom].size() - 1; // No. of nearest neighbours
390  iatomIndex = iatom; // Atomic index
391  //
392  // Loop through the nearest neighbours
393  for (int j = 1; j <= nnumNeighbours; j++) {
394  jatomID = nList[iatom][j]; // Atom ID of the nearest neighbour
395  // Get the hydrogen atom indices corresponding to the molID of jatomID
396  // Find jOxyMolID
397  auto it = idMolIDmap.find(jatomID);
398  if (it != idMolIDmap.end()) {
399  jOxyMolID = it->second;
400  } // found molecular ID of jatom oxygen atom
401  else {
402  continue;
403  } // not found
404 
405  // Find the index inside the molIDlist corresponding to the molecular ID
406  // to look for
407  listIndex = molSys::searchMolList(molIDlist, jOxyMolID);
408 
409  // Get the atom index of the oxygen atom jatom corresponding jatomID
410  auto gotJ = yCloud->idIndexMap.find(jatomID);
411  if (gotJ != yCloud->idIndexMap.end()) {
412  jatomIndex = gotJ->second;
413  } // found atom index of jatomID
414  else {
415  std::cerr << "Something is wrong with the map.\n";
416  continue;
417  } // index not found
418 
419  // Loop through the hydrogen atoms connected to jatom oxygen atom
420  for (int k = 1; k <= 2; k++) {
421  hAtomIndex = molIDlist[listIndex][k];
422  // --------
423  // Condition One: The O-H length (between the donor hydrogen atom and
424  // the acceptor oxygen atom) should be less than 2.42 Angstrom
425  // (hard-coded)
426  hBondLen =
427  bond::getHbondDistanceOH(yCloud, hCloud, iatomIndex, hAtomIndex);
428 
429  // If O-H distance is greater than or equal to 2.42 then it is not a
430  // hydrogen bond
431  if (hBondLen >= distCutoff) {
432  continue;
433  } // not a hydrogen bond
434  // --------
435  // Condition Two: The angle between the O-H and O--O vectors is less
436  // than 30 degrees (hard-coded)
437  //
438  ooVec.clear();
439  ohVec.clear();
440  // Get the O--O and O-H vectors
441  // O--O
442  ooVec.push_back(yCloud->pts[iatomIndex].x -
443  yCloud->pts[jatomIndex].x); // dx
444  ooVec.push_back(yCloud->pts[iatomIndex].y -
445  yCloud->pts[jatomIndex].y); // dy
446  ooVec.push_back(yCloud->pts[iatomIndex].z -
447  yCloud->pts[jatomIndex].z); // dz
448  // O-H
449  ohVec.push_back(yCloud->pts[iatomIndex].x -
450  hCloud->pts[hAtomIndex].x); // dx
451  ohVec.push_back(yCloud->pts[iatomIndex].y -
452  hCloud->pts[hAtomIndex].y); // dy
453  ohVec.push_back(yCloud->pts[iatomIndex].z -
454  hCloud->pts[hAtomIndex].z); // dz
455  // Apply PBCs
456  for (int l = 0; l < 3; l++) {
457  ooVec[l] -= yCloud->box[l] * round(ooVec[l] / yCloud->box[l]);
458  ohVec[l] -= yCloud->box[l] * round(ohVec[l] / yCloud->box[l]);
459  } // end of applying PBCs to the O-H and O--O vectors
460  //
461  // Get the angle between the O--O and O-H vectors
462  double eigenAngle = gen::eigenVecAngle(ooVec, ohVec);
463  double eigenAngleDeg = gen::radDeg(eigenAngle);
464 
465  //
466  // A hydrogen bond is formed if the angle is less than 30 degrees
467  if (eigenAngleDeg > angleCutoff) {
468  continue;
469  } // not a hydrogen bond
470 
471  // If you have reached this point, then O and H and indeed
472  // hydrogen-bonded. This means that jatomID should be saved in the new
473  // currentBond
474  hBondNet[iatomIndex].push_back(jatomID);
475  hBondNet[jatomIndex].push_back(iatomID);
476  break; // No need to test the other hydrogen atom if the first has
477  // been tested
478  } // end of loop through hydrogen atoms
479 
480  } // end of loop through the nearest neighbours
481 
482  } // end of loop through the neighbour list
483 
484  // --------------------
485 
486  return hBondNet;
487 }

◆ trimBonds()

std::vector< std::vector< int > > bond::trimBonds ( std::vector< std::vector< int >>  bonds)

Remove duplicate bonds.

The bond 1 2 and 2 1 are the same. To prevent multiple bonds between the same atoms, remove all bonds which are duplicates of the reversed vectors (denoting individual bonds) within the bonds vector of vectors

Parameters
[in,out]bondsRow-ordered vector of vectors of the bond matrix Each row is a ring, containing the indices of the particles in that ring

Definition at line 600 of file bond.cpp.

600  {
602  reversedBond; // Vector for the current bond in reversed order
603  std::vector<bool> isBondFlag;
604  int temp = 0;
605 
606  std::sort(bonds.begin(), bonds.end());
607  bonds.erase(std::unique(bonds.begin(), bonds.end()), bonds.end());
608 
609  // // Temp uncommented
610  // // Resize the flag vector and init to true
611  // isBondFlag.resize(bonds.size(), true);
612  // // Form of the bonds vector of vectors:
613  // // 272 214
614  // // 1 2
615  // // Meaning that 272 and 214 are bonded; similarly 1 and 2 are bonded
616 
617  // // Traverse the bonds vector of vectors
618 
619  // // Loop through all the bonds
620  // for(int ibond=0; ibond<bonds.size(); ibond++){
621 
622  // // Skip if the bond has been flagged as false
623  // if(isBondFlag[ibond] == false){continue;}
624 
625  // reversedBond.clear();
626  // reversedBond.push_back( bonds[ibond][1] );
627  // reversedBond.push_back( bonds[ibond][0] );
628 
629  // // Compare with other bonds by looping through
630  // // the rest
631  // for(int jbond=0; jbond<bonds.size(); jbond++){
632  // // Skip if jbond has been flagged
633  // if(isBondFlag[jbond] == false){continue;}
634  // // Compare reversed bond and jbond if jbond has not been flagged yet
635  // if(reversedBond==bonds[jbond]){
636  // isBondFlag[jbond] = false;
637  // temp++;
638  // } // end of check of comparison
639  // } // end of looping through all possible bonds
640  // } // end of loop through all possible bonds
641  // // temp
642 
643  return bonds;
644 }