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)
 
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 365 of file bond.cpp.

367  {
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 }
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 334 of file bond.cpp.

337  {
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 }
std::vector< S > pts
Definition: mol_sys.hpp:167
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:170
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 22 of file bond.cpp.

23  {
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 }

◆ 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 89 of file bond.cpp.

91  {
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 }
@ dummy
Definition: cage.hpp:71

◆ 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 170 of file bond.cpp.

173  {
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 }
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:334
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
int nop
Current frame number.
Definition: mol_sys.hpp:169
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::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:172
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})
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166

◆ 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 436 of file bond.cpp.

436  {
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 }