bop.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <bop.hpp>
12 #include <iostream>
13 
14 namespace bg = boost::geometry;
15 
16 // //SPHERICAL HARMONIC FUNCTIONS
17 // /********************************************/ /**
18 // * Spherical harmonics using boost
19 // ***********************************************/
20 // std::vector<std::complex<double>>
21 // sph::spheriHarmo(int orderL, std::array<double, 2> radialCoord) {
22 // // Iterate over the index of order
23 // std::vector<int> v = {-3, -2, -1, 0, 1, 2, 3};
24 // // For keeping track of the index of the output vector
25 // int i(0);
26 // std::vector<std::complex<double>> result;
27 // for (auto n : v) {
28 // auto theta = radialCoord[1];
29 // auto phi = radialCoord[0];
30 // result.resize(result.size() + 1);
31 // // This is for l=3
32 // std::complex<double> b =
33 // boost::math::spherical_harmonic(orderL, n, theta, phi);
34 // result[i] = b;
35 // // Update the index
36 // i++;
37 // }
38 // return result;
39 // }
40 
54  // For keeping track of the index of the output vector
56  std::complex<double> b; // Boost temp value
57  int m;
58 
59  result.resize(2 * orderL + 1);
60 
61  for (int k = 0; k < 2 * orderL + 1; k++) {
62  double theta = radialCoord[1];
63  double phi = radialCoord[0];
64  m = k - orderL;
65  result[k] = boost::math::spherical_harmonic(orderL, m, theta, phi);
66  }
67 
68  return result;
69 }
70 
78  // The output
79  std::array<double, 2> result;
80  // Point Definitions
81  bg::model::point<long double, 3, bg::cs::cartesian> cartesianPoint;
82  bg::model::point<long double, 3, bg::cs::spherical<bg::radian>> azuPoint;
83  // Set Value (TODO: Recurse this)
84  bg::set<0>(cartesianPoint, cartCoord[0]);
85  bg::set<1>(cartesianPoint, cartCoord[1]);
86  bg::set<2>(cartesianPoint, cartCoord[2]);
87  // Transform
88  bg::transform(cartesianPoint, azuPoint);
89  result[0] = bg::get<0>(azuPoint);
90  result[1] = bg::get<1>(azuPoint);
91  return result;
92 }
93 
104  // For keeping track of the index of the output vector
106  double theta = angles[1];
107  double phi = angles[0];
108 
109  result.resize(7);
110 
111  for (int m = 0; m < 7; m++) {
112  result[m] = sph::lookupTableQ3(m, angles);
113  }
114 
115  return result;
116 }
117 
130  std::complex<double> result(0.0, 0.0);
131  const double pi = std::acos(-1);
132  const std::complex<double> i(0.0, 1.0);
133  double constant;
134  double theta = angles[1];
135  double phi = angles[0];
136 
137  if (m == 0) {
138  constant = 0.125 * std::sqrt(35 / pi);
139  result = constant * std::exp(-3.0 * i * phi) * std::pow(sin(theta), 3.0);
140  return result;
141  } else if (m == 1) {
142  constant = 0.25 * std::sqrt(105 / (2 * pi));
143  result = constant * std::exp(-2.0 * i * phi) * std::pow(sin(theta), 2.0) *
144  cos(theta);
145  return result;
146  } else if (m == 2) {
147  constant = 0.125 * std::sqrt(21 / pi);
148  result = constant * std::exp(-1.0 * i * phi) * sin(theta) *
149  (5.0 * std::pow(cos(theta), 2.0) - 1);
150  return result;
151  } else if (m == 3) {
152  constant = 0.25 * std::sqrt(7 / pi);
153  result = constant * (5.0 * std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
154  return result;
155  } else if (m == 4) {
156  constant = -0.125 * std::sqrt(21 / pi);
157  result = constant * std::exp(i * phi) * sin(theta) *
158  (5.0 * std::pow(cos(theta), 2.0) - 1);
159  return result;
160  } else if (m == 5) {
161  constant = 0.25 * std::sqrt(105 / (2 * pi));
162  result = constant * std::exp(2.0 * i * phi) * std::pow(sin(theta), 2.0) *
163  cos(theta);
164  return result;
165  } else if (m == 6) {
166  constant = -0.125 * std::sqrt(35 / pi);
167  result = constant * std::exp(3.0 * i * phi) * std::pow(sin(theta), 3.0);
168  return result;
169  }
170 
171  return result;
172 }
173 
186  // For keeping track of the index of the output vector
188  double theta = angles[1];
189  double phi = angles[0];
190 
191  result.resize(13);
192 
193  for (int m = 0; m < 13; m++) {
194  result[m] = sph::lookupTableQ6(m, angles);
195  }
196 
197  return result;
198 }
199 
212  std::complex<double> result(0.0, 0.0);
213  const double pi = std::acos(-1);
214  const std::complex<double> i(0.0, 1.0);
215  double constant;
216  double theta = angles[1];
217  double phi = angles[0];
218 
219  if (m == 0) {
220  constant = 0.015625 * std::sqrt(3003 / pi);
221  result = constant * std::exp(-6.0 * i * phi) * std::pow(sin(theta), 6.0);
222  return result;
223  } else if (m == 1) {
224  constant = 0.09375 * std::sqrt(1001 / pi);
225  result = constant * std::exp(-5.0 * i * phi) * std::pow(sin(theta), 5.0) *
226  cos(theta);
227  return result;
228  } else if (m == 2) {
229  constant = 0.09375 * std::sqrt(91 / (2 * pi));
230  result = constant * std::exp(-4.0 * i * phi) * std::pow(sin(theta), 4.0) *
231  (11.0 * std::pow(cos(theta), 2.0) - 1);
232  return result;
233  } else if (m == 3) {
234  constant = 0.03125 * std::sqrt(1365 / pi);
235  result = constant * std::exp(-3.0 * i * phi) * std::pow(sin(theta), 3.0) *
236  (11.0 * std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
237  return result;
238  } else if (m == 4) {
239  constant = 0.015625 * std::sqrt(1365 / pi);
240  result = constant * std::exp(-2.0 * i * phi) * std::pow(sin(theta), 2.0) *
241  (33.0 * std::pow(cos(theta), 4.0) -
242  18.0 * std::pow(cos(theta), 2.0) + 1.0);
243  return result;
244  } else if (m == 5) {
245  constant = 0.0625 * std::sqrt(273 / (2 * pi));
246  result = constant * std::exp(-1.0 * i * phi) * sin(theta) *
247  (33.0 * std::pow(cos(theta), 5.0) -
248  30.0 * std::pow(cos(theta), 3.0) + 5.0 * cos(theta));
249  return result;
250  } else if (m == 6) {
251  constant = 0.03125 * std::sqrt(13 / pi);
252  result = constant * (231.0 * std::pow(cos(theta), 6.0) -
253  315.0 * std::pow(cos(theta), 4.0) +
254  105.0 * std::pow(cos(theta), 2.0) - 5.0);
255  return result;
256  } else if (m == 7) {
257  constant = -0.0625 * std::sqrt(273 / (2 * pi));
258  result = constant * std::exp(i * phi) * sin(theta) *
259  (33.0 * std::pow(cos(theta), 5.0) -
260  30.0 * std::pow(cos(theta), 3.0) + 5.0 * cos(theta));
261  return result;
262  } else if (m == 8) {
263  constant = 0.015625 * std::sqrt(1365 / pi);
264  result = constant * std::exp(2.0 * i * phi) * std::pow(sin(theta), 2.0) *
265  (33.0 * std::pow(cos(theta), 4.0) -
266  18.0 * std::pow(cos(theta), 2.0) + 1.0);
267  return result;
268  } else if (m == 9) {
269  constant = -0.03125 * std::sqrt(1365 / pi);
270  result = constant * std::exp(3.0 * i * phi) * std::pow(sin(theta), 3.0) *
271  (11.0 * std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
272  return result;
273  } else if (m == 10) {
274  constant = 0.09375 * std::sqrt(91 / (2 * pi));
275  result = constant * std::exp(4.0 * i * phi) * std::pow(sin(theta), 4.0) *
276  (11.0 * std::pow(cos(theta), 2.0) - 1);
277  return result;
278  } else if (m == 11) {
279  constant = -0.09375 * std::sqrt(1001 / pi);
280  result = constant * std::exp(5.0 * i * phi) * std::pow(sin(theta), 5.0) *
281  cos(theta);
282  return result;
283  } else if (m == 12) {
284  constant = 0.015625 * std::sqrt(3003 / pi);
285  result = constant * std::exp(6.0 * i * phi) * std::pow(sin(theta), 6.0);
286  return result;
287  }
288 
289  return result;
290 }
291 
296  std::vector<std::vector<int>> nList, bool isSlice) {
297  //
298  int l = 3; // TODO: Don't hard-code this; change later
299  int iatomID; // Atom ID (key) of iatom
300  int iatomIndex; // Index (value) of iatom
301  int jatomID; // Atom ID (key) of jatom
302  int jatomIndex; // Index (value) of nearest neighbour
303  std::array<double, 3> delta;
304  std::array<double, 2> angles;
305  chill::QlmAtom QlmTotal; // Qlm for each iatom
307  yl; // temp q_lm for each pair of iatom and jatom
308  std::complex<double> dot_product = {0, 0};
309  std::complex<double> qI = {0, 0};
310  std::complex<double> qJ = {0, 0};
311  std::complex<double> Inorm = {0, 0};
312  std::complex<double> Jnorm = {0, 0};
313  std::complex<double> complexDenominator = {0, 0};
314  std::complex<double> complexCij = {0, 0};
315  molSys::Result temp_cij; // Holds the c_ij value
316  double cij_real;
317  int nnumNeighbours; // Number of nearest neighbours for iatom
318 
319  QlmTotal.ptq.resize(yCloud->nop);
320 
321  // Loop through the neighbour list
322  for (int iatom = 0; iatom < nList.size(); iatom++) {
323  // The atom index is iatom
324  iatomIndex = iatom;
325  iatomID =
326  nList[iatomIndex][0]; // The first element in nList is the ID of iatom
327  nnumNeighbours = nList[iatomIndex].size() - 1;
328  // Now loop over the first four neighbours
329  for (int j = 1; j <= nnumNeighbours; j++) {
330  // Get the ID of jatom saved in the neighbour list
331  jatomID = nList[iatomIndex][j];
332 
333  // Get the index of jatom
334  auto it = yCloud->idIndexMap.find(jatomID);
335 
336  if (it != yCloud->idIndexMap.end()) {
337  jatomIndex = it->second;
338  } // found jatom
339  else {
340  std::cerr << "Something is wrong with the ID and index map.\n";
341  return *yCloud;
342  } // error handling
343 
344  // Get the relative distance now that the index values are known
345  delta = gen::relDist(yCloud, iatomIndex, jatomIndex);
346  double r = std::sqrt(std::pow(delta[0], 2.0) + std::pow(delta[1], 2.0) +
347  std::pow(delta[2], 2.0));
348  angles[1] = acos(delta[2] / r); // theta
349  angles[0] = atan2(delta[0], delta[1]); // phi
350 
351  // Now add over all nearest neighbours
352  if (j == 1) {
353  QlmTotal.ptq[iatomIndex].ylm = sph::spheriHarmo(3, angles);
354  // QlmTotal.ptq[iatom].ylm = sph::lookupTableQ3Vec(angles);
355  continue;
356  }
357  // Not for the first jatom
358  yl = sph::spheriHarmo(3, angles);
359  for (int m = 0; m < 2 * l + 1; m++) {
360  QlmTotal.ptq[iatomIndex].ylm[m] += yl[m];
361  // QlmTotal.ptq[iatom].ylm[m] += sph::lookupTableQ3(m, angles);
362  }
363  } // End of loop over 4 nearest neighbours
364 
365  // Divide by 4
366  QlmTotal.ptq[iatomIndex].ylm =
367  gen::avgVector(QlmTotal.ptq[iatom].ylm, l, nnumNeighbours);
368  } // End of looping over all iatom
369 
370  // ------------------------------------------------
371  // Now that you have all qlm for the particles,
372  // find c_ij
373  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
374  // if(yCloud->pts[iatom].type!=typeO){continue;}
375  // if this is a slice and the particle is not in the slice
376  // then skip
377  if (isSlice) {
378  if (yCloud->pts[iatom].inSlice == false) {
379  continue;
380  }
381  }
382  // The index is what we are looping through
383  iatomIndex = iatom;
384  nnumNeighbours = nList[iatomIndex].size() - 1;
385  yCloud->pts[iatomIndex].c_ij.reserve(nnumNeighbours);
386  // loop over the 4 nearest neighbours
387  for (int j = 1; j <= nnumNeighbours; j++) {
388  // Init to zero
389  dot_product = {0, 0};
390  Inorm = {0, 0};
391  Jnorm = {0, 0};
392  // Get ID of the nearest neighbour
393  jatomID = nList[iatomIndex][j];
394  // Get the index (value) from the ID (key)
395  auto it = yCloud->idIndexMap.find(jatomID);
396 
397  if (it != yCloud->idIndexMap.end()) {
398  jatomIndex = it->second;
399  } // found jatom
400  else {
401  std::cerr << "Something is wrong with the ID and index map.\n";
402  return *yCloud;
403  } // error handling
404  // Spherical harmonics
405  for (int m = 0; m < 2 * l + 1; m++) {
406  qI = QlmTotal.ptq[iatomIndex].ylm[m];
407  qJ = QlmTotal.ptq[jatomIndex].ylm[m];
408  dot_product = dot_product + (qI * std::conj(qJ)); // unnormalized
409  Inorm = Inorm + (qI * std::conj(qI));
410  Jnorm = Jnorm + (qJ * std::conj(qJ));
411  } // end loop over m components
412  // Get the denominator
413  complexDenominator = std::sqrt(Inorm * Jnorm);
414  complexCij = dot_product / complexDenominator;
415  // Update c_ij and type
416  cij_real = complexCij.real();
417  temp_cij.c_value = cij_real;
418  if (cij_real < -0.8) {
419  temp_cij.classifier = molSys::staggered;
420  } else if (cij_real > -0.2 && cij_real < -0.05) {
421  temp_cij.classifier = molSys::eclipsed;
422  } else {
423  temp_cij.classifier = molSys::out_of_range;
424  }
425  yCloud->pts[iatomIndex].c_ij.push_back(temp_cij);
426  } // end loop over nearest neighbours
427  }
428 
429  return *yCloud;
430 }
431 
434  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
435  std::vector<std::vector<int>> nList, bool isSlice) {
436  int ih, ic, water, interIce, unknown, total; // No. of particles of each type
437  ih = ic = water = unknown = interIce = total = 0;
438  int num_staggrd, num_eclipsd, na;
439  molSys::bond_type bondType;
440  int nnumNeighbours; // Number of nearest neighbours
441 
442  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
443  // if(yCloud->pts[iatom].type!=typeO){continue;}
444  // if this is a slice and the particle is not in the slice
445  // then skip
446  if (isSlice) {
447  if (yCloud->pts[iatom].inSlice == false) {
448  continue;
449  }
450  }
451  total++; // Update the total number of atoms considered. Change this to
452  // check for slices
453  num_staggrd = num_eclipsd = na =
454  0; // init to zero before loop through neighbours
455 
456  nnumNeighbours = nList[iatom].size() - 1;
457  // Loop through the bond cij and get the number of staggered, eclipsed bonds
458  for (int j = 0; j < nnumNeighbours; j++) {
459  bondType = yCloud->pts[iatom].c_ij[j].classifier;
460  if (bondType == molSys::eclipsed) {
461  num_eclipsd++;
462  } else if (bondType == molSys::staggered) {
463  num_staggrd++;
464  } else {
465  na++;
466  }
467  } // End of loop through neighbours
468 
469  // Add more tests later
470  yCloud->pts[iatom].iceType = molSys::unclassified; // default
471  // Cubic ice
472  // if (num_eclipsd==0 && num_staggrd==4){
473  // yCloud->pts[iatom].iceType = molSys::cubic;
474  // ic++;
475  // }
476  if (num_staggrd >= 4) {
477  yCloud->pts[iatom].iceType = molSys::cubic;
478  ic++;
479  }
480  // Hexagonal
481  else if (num_eclipsd == 1 && num_staggrd == 3) {
482  yCloud->pts[iatom].iceType = molSys::hexagonal;
483  ih++;
484  }
485  // Interfacial
486  else if (isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
487  yCloud->pts[iatom].iceType = molSys::interfacial;
488  interIce++;
489  } else {
490  yCloud->pts[iatom].iceType = molSys::water;
491  water++;
492  }
493 
494  } // End of loop through every iatom
495 
496  return *yCloud;
497 }
498 
503  int firstFrame, bool isSlice, std::string outputFileName) {
504  int ih, ic, water, interIce, unknown, total; // No. of particles of each type
505  ih = ic = water = unknown = interIce = total = 0;
506  int num_staggrd, num_eclipsd, na;
507  molSys::bond_type bondType;
508  int nnumNeighbours; // Number of nearest neighbours
509 
510  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
511  // if(yCloud->pts[iatom].type!=typeO){continue;}
512  // if this is a slice and the particle is not in the slice
513  // then skip
514  if (isSlice) {
515  if (yCloud->pts[iatom].inSlice == false) {
516  continue;
517  }
518  }
519  total++; // Update the total number of atoms considered. Change this to
520  // check for slices
521  num_staggrd = num_eclipsd = na =
522  0; // init to zero before loop through neighbours
523 
524  nnumNeighbours = nList[iatom].size() - 1;
525  // Loop through the bond cij and get the number of staggered, eclipsed bonds
526  for (int j = 0; j < nnumNeighbours; j++) {
527  bondType = yCloud->pts[iatom].c_ij[j].classifier;
528  if (bondType == molSys::eclipsed) {
529  num_eclipsd++;
530  } else if (bondType == molSys::staggered) {
531  num_staggrd++;
532  } else {
533  na++;
534  }
535  } // End of loop through neighbours
536 
537  // Add more tests later
538  yCloud->pts[iatom].iceType = molSys::unclassified; // default
539  // Cubic ice
540  // if (num_eclipsd==0 && num_staggrd==4){
541  // yCloud->pts[iatom].iceType = molSys::cubic;
542  // ic++;
543  // }
544  if (num_staggrd >= 4) {
545  yCloud->pts[iatom].iceType = molSys::cubic;
546  ic++;
547  }
548  // Hexagonal
549  else if (num_eclipsd == 1 && num_staggrd == 3) {
550  yCloud->pts[iatom].iceType = molSys::hexagonal;
551  ih++;
552  }
553  // Interfacial
554  else if (isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
555  yCloud->pts[iatom].iceType = molSys::interfacial;
556  interIce++;
557  } else {
558  yCloud->pts[iatom].iceType = molSys::water;
559  water++;
560  }
561 
562  } // End of loop through every iatom
563 
564  // water = total - ic -ih;
565 
566  // --------------------
567  // Create the directories if needed
568  sout::makePath(path);
569  std::string outputDirName = path + "bop";
570  sout::makePath(outputDirName);
571  // --------------------
572 
573  // Print to file
574  std::ofstream outputFile;
575  outputFile.open(path + "bop/" + outputFileName, std::ios_base::app);
576  // --------------------
577  // Write out the comment line for the first frame
578  if (yCloud->currentFrame == firstFrame) {
579  outputFile << "Frame Ic Ih Interfacial Water Total \n";
580  }
581  // --------------------
582  outputFile << yCloud->currentFrame << " " << ic << " " << ih << " "
583  << interIce << " " << water << " " << total << "\n";
584  outputFile.close();
585 
586  return *yCloud;
587 }
588 
599  std::vector<std::vector<int>> nList, bool isSlice) {
600  //
601  int l = 3; // TODO: Don't hard-code this; change later
602  int iatomID; // Atom ID (key) of iatom
603  int iatomIndex; // Index (value) of iatom
604  int jatomID; // Atom ID (key) of jatom
605  int jatomIndex; // Index (value) of nearest neighbour
606  std::array<double, 3> delta;
607  std::array<double, 2> angles;
608  chill::QlmAtom QlmTotal; // Qlm for each iatom
610  yl; // temp q_lm for each pair of iatom and jatom
611  std::complex<double> dot_product = {0, 0};
612  std::complex<double> qI = {0, 0};
613  std::complex<double> qJ = {0, 0};
614  std::complex<double> Inorm = {0, 0};
615  std::complex<double> Jnorm = {0, 0};
616  std::complex<double> complexDenominator = {0, 0};
617  std::complex<double> complexCij = {0, 0};
618  molSys::Result temp_cij; // Holds the c_ij value
619  double cij_real;
620  int nnumNeighbours; // Number of nearest neighbours for iatom
621 
622  QlmTotal.ptq.resize(yCloud->nop);
623 
624  // Loop through the neighbour list
625  for (int iatom = 0; iatom < nList.size(); iatom++) {
626  // The atom index is iatom
627  iatomIndex = iatom;
628  iatomID =
629  nList[iatomIndex][0]; // The first element in nList is the ID of iatom
630  nnumNeighbours = nList[iatomIndex].size() - 1;
631  // Now loop over the first four neighbours
632  for (int j = 1; j <= nnumNeighbours; j++) {
633  // Get the ID of jatom saved in the neighbour list
634  jatomID = nList[iatomIndex][j];
635 
636  // Get the index of jatom
637  auto it = yCloud->idIndexMap.find(jatomID);
638 
639  if (it != yCloud->idIndexMap.end()) {
640  jatomIndex = it->second;
641  } // found jatom
642  else {
643  std::cerr << "Something is wrong with the ID and index map.\n";
644  return *yCloud;
645  } // error handling
646 
647  // Get the relative distance now that the index values are known
648  delta = gen::relDist(yCloud, iatomIndex, jatomIndex);
649  double r = std::sqrt(std::pow(delta[0], 2.0) + std::pow(delta[1], 2.0) +
650  std::pow(delta[2], 2.0));
651  angles[1] = acos(delta[2] / r); // theta
652  angles[0] = atan2(delta[0], delta[1]); // phi
653 
654  // Now add over all nearest neighbours
655  if (j == 1) {
656  QlmTotal.ptq[iatomIndex].ylm = sph::spheriHarmo(3, angles);
657  // QlmTotal.ptq[iatom].ylm = sph::lookupTableQ3Vec(angles);
658  continue;
659  }
660  // Not for the first jatom
661  yl = sph::spheriHarmo(3, angles);
662  for (int m = 0; m < 2 * l + 1; m++) {
663  QlmTotal.ptq[iatomIndex].ylm[m] += yl[m];
664  // QlmTotal.ptq[iatom].ylm[m] += sph::lookupTableQ3(m, angles);
665  }
666  } // End of loop over 4 nearest neighbours
667 
668  // Divide by 4
669  QlmTotal.ptq[iatomIndex].ylm =
670  gen::avgVector(QlmTotal.ptq[iatom].ylm, l, nnumNeighbours);
671  } // End of looping over all iatom
672 
673  // ------------------------------------------------
674  // Now that you have all qlm for the particles,
675  // find c_ij
676  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
677  // if(yCloud->pts[iatom].type!=typeO){continue;}
678  // if this is a slice and the particle is not in the slice
679  // then skip
680  if (isSlice) {
681  if (yCloud->pts[iatom].inSlice == false) {
682  continue;
683  }
684  }
685  // The index is what we are looping through
686  iatomIndex = iatom;
687  nnumNeighbours = nList[iatomIndex].size() - 1;
688  yCloud->pts[iatomIndex].c_ij.reserve(nnumNeighbours);
689  // loop over the 4 nearest neighbours
690  for (int j = 1; j <= nnumNeighbours; j++) {
691  // Init to zero
692  dot_product = {0, 0};
693  Inorm = {0, 0};
694  Jnorm = {0, 0};
695  // Get ID of the nearest neighbour
696  jatomID = nList[iatomIndex][j];
697  // Get the index (value) from the ID (key)
698  auto it = yCloud->idIndexMap.find(jatomID);
699 
700  if (it != yCloud->idIndexMap.end()) {
701  jatomIndex = it->second;
702  } // found jatom
703  else {
704  std::cerr << "Something is wrong with the ID and index map.\n";
705  return *yCloud;
706  } // error handling
707  // Spherical harmonics
708  for (int m = 0; m < 2 * l + 1; m++) {
709  qI = QlmTotal.ptq[iatomIndex].ylm[m];
710  qJ = QlmTotal.ptq[jatomIndex].ylm[m];
711  dot_product = dot_product + (qI * std::conj(qJ)); // unnormalized
712  Inorm = Inorm + (qI * std::conj(qI));
713  Jnorm = Jnorm + (qJ * std::conj(qJ));
714  } // end loop over m components
715  // Get the denominator
716  complexDenominator = std::sqrt(Inorm * Jnorm);
717  complexCij = dot_product / complexDenominator;
718  // Update c_ij and type
719  cij_real = complexCij.real();
720  temp_cij.c_value = cij_real;
721  if (cij_real <= -0.8) {
722  temp_cij.classifier = molSys::staggered;
723  } else if (cij_real >= -0.35 && cij_real <= 0.25) {
724  temp_cij.classifier = molSys::eclipsed;
725  } else {
726  temp_cij.classifier = molSys::out_of_range;
727  }
728  yCloud->pts[iatomIndex].c_ij.push_back(temp_cij);
729  } // end loop over nearest neighbours
730  }
731 
732  // ------------------------------------------------
733 
734  return *yCloud;
735 }
736 
753  int firstFrame, bool isSlice,
754  std::string outputFileName) {
755  int ih, ic, interIce, water, unknown, clath, interClath,
756  total; // No. of particles of each type
757  ih = ic = water = unknown = interIce = total = 0;
758  clath = interClath = 0;
759  int num_staggrd, num_eclipsd, na;
760  molSys::bond_type bondType;
761  int nnumNeighbours; // number of nearest neighbours
762 
763  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
764  // if(yCloud->pts[iatom].type!=typeO){continue;}
765  // if this is a slice and the particle is not in the slice
766  // then skip
767  if (isSlice) {
768  if (yCloud->pts[iatom].inSlice == false) {
769  continue;
770  }
771  }
772  total++; // Update the total number of atoms considered. Change this to a
773  // check for slices
774  nnumNeighbours =
775  nList[iatom].size() - 1; // number of nearest neighbours (total -1)
776  num_staggrd = num_eclipsd = na =
777  0; // init to zero before loop through neighbours
778  // Loop through the bond cij and get the number of staggered, eclipsed bonds
779  for (int j = 0; j < nnumNeighbours; j++) {
780  bondType = yCloud->pts[iatom].c_ij[j].classifier;
781  if (bondType == molSys::eclipsed) {
782  num_eclipsd++;
783  } else if (bondType == molSys::staggered) {
784  num_staggrd++;
785  } else {
786  na++;
787  }
788  } // End of loop through neighbours
789 
790  // Add more tests later
791  yCloud->pts[iatom].iceType = molSys::unclassified; // default
792  if (nnumNeighbours == 4) {
793  // Cubic ice
794  if (num_eclipsd == 0 && num_staggrd == 4) {
795  yCloud->pts[iatom].iceType = molSys::cubic;
796  ic++;
797  }
798  // Hexagonal
799  else if (num_eclipsd == 1 && num_staggrd == 3) {
800  yCloud->pts[iatom].iceType = molSys::hexagonal;
801  ih++;
802  }
803  // Interfacial
804  else if (isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
805  yCloud->pts[iatom].iceType = molSys::interfacial;
806  interIce++;
807  }
808  // Clathrate
809  else if (num_eclipsd == 4 && num_staggrd == 0) {
810  yCloud->pts[iatom].iceType = molSys::clathrate;
811  clath++;
812  }
813  // Interfacial clathrate
814  else if (num_eclipsd == 3) {
815  yCloud->pts[iatom].iceType = molSys::interClathrate;
816  interClath++;
817  }
818  // Water
819  else {
820  yCloud->pts[iatom].iceType = molSys::water;
821  water++;
822  }
823  } else {
824  yCloud->pts[iatom].iceType = molSys::water;
825  water++;
826  }
827 
828  } // End of loop through every iatom
829 
830  // water = total - ic -ih;
831 
832  // --------------------
833  // Create the directories if needed
834  sout::makePath(path);
835  std::string outputDirName = path + "bop";
836  sout::makePath(outputDirName);
837  // --------------------
838 
839  std::ofstream outputFile;
840  outputFile.open(path + "bop/" + outputFileName, std::ios_base::app);
841  // --------------------
842  // Comment line for the first line
843  if (yCloud->currentFrame == firstFrame) {
844  outputFile << "Frame Ic Ih Interfacial Clath InterClath Water Total\n";
845  }
846  // --------------------
847  outputFile << yCloud->currentFrame << " " << ic << " " << ih << " "
848  << interIce << " " << clath << " " << interClath << " " << water
849  << " " << total << "\n";
850  outputFile.close();
851 
852  return *yCloud;
853 }
854 
855 // TODO: Add code for slices!
865  std::vector<std::vector<int>> nList, bool isSlice) {
866  //
867  int l = 6; // We're using q6 here
868  int jatomID; // Atom ID of the nearest neighbour
869  int jatomIndex; // Index of nearest neighbour
870  std::array<double, 3> delta;
871  std::array<double, 2> angles;
872  chill::QlmAtom QlmTotal; // Qlm for each iatom
874  yl; // temp q_lm for each pair of iatom and jatom
875  std::complex<double> dot_product = {0, 0};
876  std::complex<double> qI = {0, 0};
877  std::complex<double> qJ = {0, 0};
878  std::complex<double> Inorm = {0, 0};
879  std::complex<double> Jnorm = {0, 0};
880  std::complex<double> complexDenominator = {0, 0};
881  std::complex<double> complexCij = {0, 0};
882  double cij_real;
883  int nnumNeighbours;
884  std::vector<double> resultQ; // Vector with averaged q values
885  double q_value = 0.0; // Averaged q value per neighbour pair
886 
887  QlmTotal.ptq.resize(yCloud->nop);
888  resultQ.resize(yCloud->nop);
889 
890  // Loop through every index in yCloud
891  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
892  // if(yCloud->pts[iatom].type!=typeO){continue;}
893 
894  nnumNeighbours = nList[iatom].size() - 1; // One less than the actual length
895  // Now loop over the first four neighbours
896  for (int j = 1; j <= nnumNeighbours; j++) {
897  // Get the atom ID
898  jatomID = nList[iatom][j]; // Atom ID (key)
899 
900  // Get the atom index (value) from the atom ID (key)
901  auto it = yCloud->idIndexMap.find(jatomID);
902 
903  if (it != yCloud->idIndexMap.end()) {
904  jatomIndex = it->second;
905  } else {
906  std::cerr << "Your map must be wrong.\n";
907  return resultQ; // return with error
908  }
909 
910  // Get the relative distances
911  delta = gen::relDist(yCloud, iatom, jatomIndex);
912 
913  // angles = sph::radialCoord(delta);
914  double r = std::sqrt(std::pow(delta[0], 2.0) + std::pow(delta[1], 2.0) +
915  std::pow(delta[2], 2.0));
916  angles[1] = acos(delta[2] / r); // theta
917  angles[0] = atan2(delta[0], delta[1]); // phi
918 
919  // Now add over all nearest neighbours
920  if (j == 1) {
921  // QlmTotal.ptq[iatom].ylm = sph::spheriHarmo(l, angles);
922  QlmTotal.ptq[iatom].ylm = sph::lookupTableQ6Vec(angles);
923  continue;
924  }
925  // Not for the first jatom
926  yl = sph::spheriHarmo(l, angles);
927  for (int m = 0; m < 2 * l + 1; m++) {
928  QlmTotal.ptq[iatom].ylm[m] += yl[m];
929  // QlmTotal.ptq[iatom].ylm[m] += sph::lookupTableQ6(m, angles);
930  }
931  } // End of loop over 4 nearest neighbours
932 
933  // Divide by 4
934  QlmTotal.ptq[iatom].ylm =
935  gen::avgVector(QlmTotal.ptq[iatom].ylm, l, nnumNeighbours);
936  } // End of looping over all iatom
937 
938  // ------------------------------------------------
939  // Now that you have all qlm for the particles,
940  // find c_ij
941  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
942  // if(yCloud->pts[iatom].type!=typeO){continue;}
943  // if this is a slice and the particle is not in the slice
944  // then skip TODO: UNCOMMENT
945  // if(isSlice){
946  // if(yCloud->pts[iatom].inSlice==false){continue;}
947  // }
948 
949  nnumNeighbours = nList[iatom].size() - 1; // Number of nearest neighbours
950  q_value = 0.0; // initialize to zero
951  // yCloud->pts[iatom].c_ij.reserve(nnumNeighbours);
952  // loop over the 4 nearest neighbours
953  for (int j = 1; j <= nnumNeighbours; j++) {
954  // Init to zero
955  dot_product = {0, 0};
956  Inorm = {0, 0};
957  Jnorm = {0, 0};
958  // Get index of the nearest neighbour!
959  jatomID = nList[iatom][j]; // Atom ID (key)
960 
961  // Get the index jatomIndex
962  auto it = yCloud->idIndexMap.find(jatomID);
963 
964  if (it != yCloud->idIndexMap.end()) {
965  jatomIndex = it->second;
966  } // end of getting the index of jatom
967 
968  for (int m = 0; m < 2 * l + 1; m++) {
969  qI = QlmTotal.ptq[iatom].ylm[m];
970  qJ = QlmTotal.ptq[jatomIndex].ylm[m];
971  dot_product = dot_product + (qI * std::conj(qJ)); // unnormalized
972  Inorm = Inorm + (qI * std::conj(qI));
973  Jnorm = Jnorm + (qJ * std::conj(qJ));
974  } // end loop over m components
975  // Get the denominator
976  complexDenominator = std::sqrt(Inorm * Jnorm);
977  complexCij = dot_product / complexDenominator;
978  // Update c_ij and type
979  cij_real = complexCij.real();
980 
981  q_value += cij_real;
982 
983  } // end loop over nearest neighbours
984 
985  // Average q_value over all nearest neighbours
986  q_value /= (double)nnumNeighbours;
987 
988  resultQ[iatom] = q_value; // Update the vector of averaged q6
989  }
990 
991  // ------------------------------------------------
992 
993  return resultQ;
994 }
995 
1008  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1009  std::vector<double> *q6) {
1010  // If averaged q6 > 0.5, then consider it to be ice
1011  // If averaged q3 < -0.75 then it is ih or ic. If q3 < -0.85 then it is cubic,
1012  // otherwise it is hexagonal
1013  double avgQ3 = 0.0;
1014  int nnumNeighbours;
1015 
1016  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
1017  // Check if it has been classified as water
1018  if (yCloud->pts[iatom].iceType == molSys::water) {
1019  if ((*q6)[iatom] > 0.5) {
1020  avgQ3 = 0.0; // init to zero
1021  // Loop through all c_ij
1022  nnumNeighbours = yCloud->pts[iatom].c_ij.size();
1023  for (int j = 0; j < nnumNeighbours; j++) {
1024  avgQ3 += yCloud->pts[iatom].c_ij[j].c_value;
1025  }
1026  avgQ3 /= (double)nnumNeighbours;
1027 
1028  // If averaged q3 < -0.75, then reclassify
1029  if (avgQ3 <= -0.75) {
1030  if (avgQ3 < -0.85) {
1031  yCloud->pts[iatom].iceType = molSys::reCubic;
1032  } // molSys::cubic
1033  else {
1034  yCloud->pts[iatom].iceType = molSys::reHex;
1035  } // molSys::hexagonal
1036  } // end of reclassification
1037  } // check for solid atom!
1038  } // end of check for water
1039  } // End loop through every iatom
1040 
1041  return *yCloud;
1042 }
1043 
1057  molSys::PointCloud<molSys::Point<double>, double> *yCloud, std::string path,
1058  int firstFrame, bool isSlice, std::string outputFileName) {
1059  int ih, ic, interIce, water, unknown, clath, interClath,
1060  total; // No. of particles of each type
1061  ih = ic = water = unknown = interIce = total = 0;
1062  clath = interClath = 0;
1063 
1064  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
1065  // if(yCloud->pts[iatom].type != typeO){continue;}
1066  // check for slice
1067  if (isSlice) {
1068  if (yCloud->pts[iatom].inSlice == false) {
1069  continue;
1070  }
1071  }
1072  total++;
1073  if (yCloud->pts[iatom].iceType == molSys::cubic) {
1074  ic++;
1075  } else if (yCloud->pts[iatom].iceType == molSys::hexagonal) {
1076  ih++;
1077  } else if (yCloud->pts[iatom].iceType == molSys::water) {
1078  water++;
1079  } else if (yCloud->pts[iatom].iceType == molSys::interfacial) {
1080  interIce++;
1081  } else if (yCloud->pts[iatom].iceType == molSys::clathrate) {
1082  clath++;
1083  } else if (yCloud->pts[iatom].iceType == molSys::interClathrate) {
1084  interClath++;
1085  } else {
1086  unknown++;
1087  }
1088  }
1089 
1090  // --------------------
1091  // Create the directories if needed
1092  sout::makePath(path);
1093  std::string outputDirName = path + "bop";
1094  sout::makePath(outputDirName);
1095  // --------------------
1096  // Print to file
1097  std::ofstream outputFile;
1098  outputFile.open(path + "bop/" + outputFileName, std::ios_base::app);
1099  // --------------------
1100  // Write out the comment line
1101  if (yCloud->currentFrame == firstFrame) {
1102  outputFile << "Frame Ic Ih Interfacial Clath InterClath Water Total\n";
1103  }
1104  // --------------------
1105  outputFile << yCloud->currentFrame << " " << ic << " " << ih << " "
1106  << interIce << " " << clath << " " << interClath << " " << water
1107  << " " << total << "\n";
1108  outputFile.close();
1109 
1110  return 0;
1111 }
1112 
1126  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1127  std::vector<std::vector<int>> nList, int iatom, int num_staggrd,
1128  int num_eclipsd) {
1129  int nnumNeighbours =
1130  nList[iatom].size() - 1; // number of nearest neighbours of iatom
1131  int neighStaggered =
1132  0; // number of staggered bonds in the neighbours of iatom
1133  int jatomID; // ID of the nearest neighbour
1134  int jatomIndex; // Index (value) of nearest neighbour
1135 
1136  // INTERFACIAL
1137  // Condition 1 : only two staggered bonds and at least
1138  // one neighbor with more than two staggered bonds
1139  if (num_staggrd == 2) {
1140  // Loop over the nearest neighbours
1141  for (int j = 1; j <= nnumNeighbours; j++) {
1142  // Get index of the nearest neighbour
1143  jatomID = nList[iatom][j];
1144  // Get the index (value) from the ID (key)
1145  auto it = yCloud->idIndexMap.find(jatomID);
1146 
1147  if (it != yCloud->idIndexMap.end()) {
1148  jatomIndex = it->second;
1149  } else {
1150  std::cerr << "Something is gravely wrong with your map.\n";
1151  return false;
1152  }
1153  //
1154  neighStaggered = chill::numStaggered(yCloud, nList, jatomIndex);
1155  if (neighStaggered > 2) {
1156  return true;
1157  }
1158  } // End loop over nearest neighbours
1159  } // end condition 1
1160  // Condition 2 : three staggered bonds, no eclipsed bond,
1161  // and at least one neighbor with two staggered bonds
1162  if (num_staggrd == 3 && num_eclipsd == 0) {
1163  // Loop over the nearest neighbours
1164  for (int j = 1; j <= nnumNeighbours; j++) {
1165  // Get index of the nearest neighbour
1166  // ID of the nearest neighbour
1167  jatomID = nList[iatom][j];
1168  // Get the index (value) from the ID (key)
1169  auto it = yCloud->idIndexMap.find(jatomID);
1170 
1171  if (it != yCloud->idIndexMap.end()) {
1172  jatomIndex = it->second;
1173  } else {
1174  std::cerr << "Something is gravely wrong with your map.\n";
1175  return false;
1176  }
1177  //
1178  neighStaggered = chill::numStaggered(yCloud, nList, jatomIndex);
1179  if (neighStaggered == 2) {
1180  return true;
1181  }
1182  }
1183  }
1184  // not interfacial
1185  return false;
1186 }
1187 
1198  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1199  std::vector<std::vector<int>> nList, int jatom) {
1200  int num_staggrd = 0; // Number of staggered bonds
1201  molSys::bond_type bondType; // Bond type
1202  int num_bonds; // No. of bonds of the jatom
1203  int nnumNeighbours =
1204  nList[jatom].size() - 1; // No. of nearest neighbours of index jatom
1205 
1206  // Loop over all bonds
1207  for (int i = 0; i < nnumNeighbours; i++) {
1208  bondType = yCloud->pts[jatom].c_ij[i].classifier;
1209  // If the bond is staggered increment the number of staggered bonds
1210  if (bondType == molSys::staggered) {
1211  num_staggrd++;
1212  }
1213  } // end of loop over c_ij
1214 
1215  return num_staggrd;
1216 }
T acos(T... args)
File for the bond order parameter analysis.
T close(T... args)
T exp(T... args)
std::vector< std::complex< double > > lookupTableQ6Vec(std::array< double, 2 > angles)
Lookup table for Q6.
Definition: bop.cpp:185
std::vector< double > getq6(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Definition: bop.cpp:864
std::vector< std::complex< double > > lookupTableQ3Vec(std::array< double, 2 > angles)
Lookup table for Q3.
Definition: bop.cpp:103
molSys::PointCloud< molSys::Point< double >, double > getCorrelPlus(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Gets c_ij and then classifies bond types according to the CHILL+ algorithm.
Definition: bop.cpp:598
molSys::PointCloud< molSys::Point< double >, double > getIceTypePlus(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="chillPlus.txt")
Classifies each atom according to the CHILL+ algorithm.
Definition: bop.cpp:751
std::vector< YlmAtom > ptq
Definition: bop.hpp:166
molSys::PointCloud< molSys::Point< double >, double > getCorrel(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Definition: bop.cpp:295
std::array< double, 2 > radialCoord(std::array< double, 3 > cartCoord)
Definition: bop.cpp:77
int printIceType(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="superChill.txt")
Prints out the iceType for a particular frame onto the terminal.
Definition: bop.cpp:1056
std::complex< double > lookupTableQ6(int m, std::array< double, 2 > angles)
Lookup table for Q6 (m=0 to m=12)
Definition: bop.cpp:211
bool isInterfacial(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int iatom, int num_staggrd, int num_eclipsd)
Definition: bop.cpp:1125
std::vector< std::complex< double > > spheriHarmo(int orderL, std::array< double, 2 > radialCoord)
Definition: bop.cpp:53
int numStaggered(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int jatom)
Finds the number of staggered bonds for a given atom of index jatom.
Definition: bop.cpp:1197
std::complex< double > lookupTableQ3(int m, std::array< double, 2 > angles)
Lookup table for Q3 (m=0 to m=6)
Definition: bop.cpp:129
molSys::PointCloud< molSys::Point< double >, double > getIceTypeNoPrint(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
Classifies each atom according to the CHILL algorithm without printing.
Definition: bop.cpp:433
molSys::PointCloud< molSys::Point< double >, double > reclassifyWater(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *q6)
Definition: bop.cpp:1007
molSys::PointCloud< molSys::Point< double >, double > getIceType(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="chill.txt")
Classifies each atom according to the CHILL algorithm.
Definition: bop.cpp:501
std::vector< std::complex< double > > avgVector(std::vector< std::complex< double >> v, int l, int neigh)
Definition: generic.hpp:309
std::array< double, 3 > relDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Definition: generic.hpp:196
const double pi
Definition: generic.hpp:50
bond_type
Definition: mol_sys.hpp:71
bond_type classifier
Definition: mol_sys.hpp:129
double c_value
Classifier according to CHILL, CHILL+ etc.
Definition: mol_sys.hpp:130
@ cubic
Ic, or particle type signifying Cubic Ice.
Definition: mol_sys.hpp:110
@ hexagonal
Ih, or particle type signifying Hexagonal Ice.
Definition: mol_sys.hpp:111
@ interClathrate
Interfacial clathrate ice phase.
Definition: mol_sys.hpp:115
@ interfacial
Interfacial ice: ice-like molecules which do not fulfill the strict criteria of the Ic or Ih phases.
Definition: mol_sys.hpp:113
@ reCubic
Reclassified as cubic ice, according to the order parameter.
Definition: mol_sys.hpp:117
@ clathrate
Clathrate ice phase.
Definition: mol_sys.hpp:114
@ water
Liquid/amorphous phase.
Definition: mol_sys.hpp:112
@ reHex
Reclassified as hexagonal ice, according to the order parameter.
Definition: mol_sys.hpp:118
@ unclassified
Not classified into any other category.
Definition: mol_sys.hpp:116
@ eclipsed
The bond is an eclipsed bond.
Definition: mol_sys.hpp:71
@ out_of_range
The bond cannot be classified as either staggered or eclipsed.
Definition: mol_sys.hpp:71
@ staggered
The bond is a staggered bond, according to the or value.
Definition: mol_sys.hpp:71
int makePath(const std::string &path)
T open(T... args)
T pow(T... args)
T real(T... args)
T resize(T... args)
T sqrt(T... args)
This is the local orientational bond order parameter , of length .
Definition: bop.hpp:165
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
This contains the bond classifier of enum type bond_type, and the bond correlation factor.
Definition: mol_sys.hpp:128