•All Data Structures Namespaces Files Functions Variables Enumerations Enumerator Modules Pages
Namespaces | Functions
Bond

Namespaces

namespace  bond
 Functions for bond-related analyses.
 

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.
 

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.

Code
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}

◆ 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.

Code
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}
molSys::PointCloud::pts
std::vector< S > pts
Definition mol_sys.hpp:171
molSys::PointCloud::box
std::vector< T > box
Number of atoms.
Definition mol_sys.hpp:174

◆ 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.

Code
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.

Code
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}
cage::iceType::dummy
@ dummy

◆ 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.

Code
177 {
178 //
179 std::vector<std::vector<int>>
180 hBondNet; // Output vector of vectors containing the HBN
181 molSys::PointCloud<molSys::Point<double>, double>
182 hCloud; // point Cloud for the hydrogen atoms
183 std::vector<std::vector<int>>
184 molIDlist; // Vector of vectors; first element is the molID, and the next
185 // two elements are the hydrogen atom indices
186 std::unordered_map<int, int>
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}
bond::getHbondDistanceOH
double getHbondDistanceOH(molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *hCloud, int oAtomIndex, int hAtomIndex)
Definition bond.cpp:498
gen::eigenVecAngle
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
gen::radDeg
double radDeg(double angle)
Definition generic.hpp:61
molSys::PointCloud::nop
int nop
Current frame number.
Definition mol_sys.hpp:173
molSys::hAtomMolList
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
molSys::createIDMolIDmap
std::unordered_map< int, int > createIDMolIDmap(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition mol_sys.cpp:44
molSys::searchMolList
int searchMolList(std::vector< std::vector< int > > molList, int molIDtoFind)
Definition mol_sys.cpp:131
molSys::PointCloud::idIndexMap
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition mol_sys.hpp:176
molSys::clearPointCloud
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
sinp::readLammpsTrjreduced
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})
Definition seams_input.cpp:656
molSys::PointCloud
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.

Code
345 {
346 //
347 std::vector<std::vector<int>>
348 hBondNet; // Output vector of vectors containing the HBN
349 std::vector<std::vector<int>>
350 molIDlist; // Vector of vectors; first element is the molID, and the next
351 // two elements are the hydrogen atom indices
352 std::unordered_map<int, int>
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.

Code
600 {
601 std::vector<int>
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}
Generated by Doxygen 1.9.8 and the doxyYoda theme