14 namespace bg = boost::geometry;
59 result.
resize(2 * orderL + 1);
61 for (
int k = 0; k < 2 * orderL + 1; k++) {
65 result[k] = boost::math::spherical_harmonic(orderL, m, theta, phi);
81 bg::model::point<long double, 3, bg::cs::cartesian> cartesianPoint;
82 bg::model::point<long double, 3, bg::cs::spherical<bg::radian>> azuPoint;
84 bg::set<0>(cartesianPoint, cartCoord[0]);
85 bg::set<1>(cartesianPoint, cartCoord[1]);
86 bg::set<2>(cartesianPoint, cartCoord[2]);
88 bg::transform(cartesianPoint, azuPoint);
89 result[0] = bg::get<0>(azuPoint);
90 result[1] = bg::get<1>(azuPoint);
106 double theta = angles[1];
107 double phi = angles[0];
111 for (
int m = 0; m < 7; m++) {
134 double theta = angles[1];
135 double phi = angles[0];
148 result = constant *
std::exp(-1.0 * i * phi) * sin(theta) *
149 (5.0 *
std::pow(cos(theta), 2.0) - 1);
153 result = constant * (5.0 *
std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
157 result = constant *
std::exp(i * phi) * sin(theta) *
158 (5.0 *
std::pow(cos(theta), 2.0) - 1);
188 double theta = angles[1];
189 double phi = angles[0];
193 for (
int m = 0; m < 13; m++) {
216 double theta = angles[1];
217 double phi = angles[0];
231 (11.0 *
std::pow(cos(theta), 2.0) - 1);
236 (11.0 *
std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
242 18.0 *
std::pow(cos(theta), 2.0) + 1.0);
246 result = constant *
std::exp(-1.0 * i * phi) * sin(theta) *
248 30.0 *
std::pow(cos(theta), 3.0) + 5.0 * cos(theta));
252 result = constant * (231.0 *
std::pow(cos(theta), 6.0) -
254 105.0 *
std::pow(cos(theta), 2.0) - 5.0);
258 result = constant *
std::exp(i * phi) * sin(theta) *
260 30.0 *
std::pow(cos(theta), 3.0) + 5.0 * cos(theta));
266 18.0 *
std::pow(cos(theta), 2.0) + 1.0);
271 (11.0 *
std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
273 }
else if (m == 10) {
276 (11.0 *
std::pow(cos(theta), 2.0) - 1);
278 }
else if (m == 11) {
283 }
else if (m == 12) {
319 QlmTotal.
ptq.resize(yCloud->nop);
322 for (
int iatom = 0; iatom < nList.size(); iatom++) {
326 nList[iatomIndex][0];
327 nnumNeighbours = nList[iatomIndex].size() - 1;
329 for (
int j = 1; j <= nnumNeighbours; j++) {
331 jatomID = nList[iatomIndex][j];
334 auto it = yCloud->idIndexMap.find(jatomID);
336 if (it != yCloud->idIndexMap.end()) {
337 jatomIndex = it->second;
340 std::cerr <<
"Something is wrong with the ID and index map.\n";
348 angles[1] = acos(delta[2] / r);
349 angles[0] = atan2(delta[0], delta[1]);
359 for (
int m = 0; m < 2 * l + 1; m++) {
360 QlmTotal.
ptq[iatomIndex].ylm[m] += yl[m];
366 QlmTotal.
ptq[iatomIndex].ylm =
373 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
378 if (yCloud->pts[iatom].inSlice ==
false) {
384 nnumNeighbours = nList[iatomIndex].size() - 1;
385 yCloud->pts[iatomIndex].c_ij.reserve(nnumNeighbours);
387 for (
int j = 1; j <= nnumNeighbours; j++) {
389 dot_product = {0, 0};
393 jatomID = nList[iatomIndex][j];
395 auto it = yCloud->idIndexMap.find(jatomID);
397 if (it != yCloud->idIndexMap.end()) {
398 jatomIndex = it->second;
401 std::cerr <<
"Something is wrong with the ID and index map.\n";
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));
409 Inorm = Inorm + (qI * std::conj(qI));
410 Jnorm = Jnorm + (qJ * std::conj(qJ));
413 complexDenominator =
std::sqrt(Inorm * Jnorm);
414 complexCij = dot_product / complexDenominator;
416 cij_real = complexCij.
real();
418 if (cij_real < -0.8) {
420 }
else if (cij_real > -0.2 && cij_real < -0.05) {
425 yCloud->pts[iatomIndex].c_ij.push_back(temp_cij);
436 int ih, ic,
water, interIce, unknown, total;
437 ih = ic =
water = unknown = interIce = total = 0;
438 int num_staggrd, num_eclipsd, na;
442 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
447 if (yCloud->pts[iatom].inSlice ==
false) {
453 num_staggrd = num_eclipsd = na =
456 nnumNeighbours = nList[iatom].size() - 1;
458 for (
int j = 0; j < nnumNeighbours; j++) {
459 bondType = yCloud->pts[iatom].c_ij[j].classifier;
476 if (num_staggrd >= 4) {
481 else if (num_eclipsd == 1 && num_staggrd == 3) {
486 else if (
isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
503 int firstFrame,
bool isSlice,
std::string outputFileName) {
504 int ih, ic,
water, interIce, unknown, total;
505 ih = ic =
water = unknown = interIce = total = 0;
506 int num_staggrd, num_eclipsd, na;
510 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
515 if (yCloud->pts[iatom].inSlice ==
false) {
521 num_staggrd = num_eclipsd = na =
524 nnumNeighbours = nList[iatom].size() - 1;
526 for (
int j = 0; j < nnumNeighbours; j++) {
527 bondType = yCloud->pts[iatom].c_ij[j].classifier;
544 if (num_staggrd >= 4) {
549 else if (num_eclipsd == 1 && num_staggrd == 3) {
554 else if (
isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
575 outputFile.
open(path +
"bop/" + outputFileName, std::ios_base::app);
578 if (yCloud->currentFrame == firstFrame) {
579 outputFile <<
"Frame Ic Ih Interfacial Water Total \n";
582 outputFile << yCloud->currentFrame <<
" " << ic <<
" " << ih <<
" "
583 << interIce <<
" " <<
water <<
" " << total <<
"\n";
622 QlmTotal.
ptq.resize(yCloud->nop);
625 for (
int iatom = 0; iatom < nList.size(); iatom++) {
629 nList[iatomIndex][0];
630 nnumNeighbours = nList[iatomIndex].size() - 1;
632 for (
int j = 1; j <= nnumNeighbours; j++) {
634 jatomID = nList[iatomIndex][j];
637 auto it = yCloud->idIndexMap.find(jatomID);
639 if (it != yCloud->idIndexMap.end()) {
640 jatomIndex = it->second;
643 std::cerr <<
"Something is wrong with the ID and index map.\n";
651 angles[1] = acos(delta[2] / r);
652 angles[0] = atan2(delta[0], delta[1]);
662 for (
int m = 0; m < 2 * l + 1; m++) {
663 QlmTotal.
ptq[iatomIndex].ylm[m] += yl[m];
669 QlmTotal.
ptq[iatomIndex].ylm =
676 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
681 if (yCloud->pts[iatom].inSlice ==
false) {
687 nnumNeighbours = nList[iatomIndex].size() - 1;
688 yCloud->pts[iatomIndex].c_ij.reserve(nnumNeighbours);
690 for (
int j = 1; j <= nnumNeighbours; j++) {
692 dot_product = {0, 0};
696 jatomID = nList[iatomIndex][j];
698 auto it = yCloud->idIndexMap.find(jatomID);
700 if (it != yCloud->idIndexMap.end()) {
701 jatomIndex = it->second;
704 std::cerr <<
"Something is wrong with the ID and index map.\n";
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));
712 Inorm = Inorm + (qI * std::conj(qI));
713 Jnorm = Jnorm + (qJ * std::conj(qJ));
716 complexDenominator =
std::sqrt(Inorm * Jnorm);
717 complexCij = dot_product / complexDenominator;
719 cij_real = complexCij.
real();
721 if (cij_real <= -0.8) {
723 }
else if (cij_real >= -0.35 && cij_real <= 0.25) {
728 yCloud->pts[iatomIndex].c_ij.push_back(temp_cij);
753 int firstFrame,
bool isSlice,
755 int ih, ic, interIce,
water, unknown, clath, interClath,
757 ih = ic =
water = unknown = interIce = total = 0;
758 clath = interClath = 0;
759 int num_staggrd, num_eclipsd, na;
763 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
768 if (yCloud->pts[iatom].inSlice ==
false) {
775 nList[iatom].size() - 1;
776 num_staggrd = num_eclipsd = na =
779 for (
int j = 0; j < nnumNeighbours; j++) {
780 bondType = yCloud->pts[iatom].c_ij[j].classifier;
792 if (nnumNeighbours == 4) {
794 if (num_eclipsd == 0 && num_staggrd == 4) {
799 else if (num_eclipsd == 1 && num_staggrd == 3) {
804 else if (
isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
809 else if (num_eclipsd == 4 && num_staggrd == 0) {
814 else if (num_eclipsd == 3) {
840 outputFile.
open(path +
"bop/" + outputFileName, std::ios_base::app);
843 if (yCloud->currentFrame == firstFrame) {
844 outputFile <<
"Frame Ic Ih Interfacial Clath InterClath Water Total\n";
847 outputFile << yCloud->currentFrame <<
" " << ic <<
" " << ih <<
" "
848 << interIce <<
" " << clath <<
" " << interClath <<
" " <<
water
849 <<
" " << total <<
"\n";
885 double q_value = 0.0;
887 QlmTotal.
ptq.resize(yCloud->nop);
888 resultQ.
resize(yCloud->nop);
891 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
894 nnumNeighbours = nList[iatom].size() - 1;
896 for (
int j = 1; j <= nnumNeighbours; j++) {
898 jatomID = nList[iatom][j];
901 auto it = yCloud->idIndexMap.find(jatomID);
903 if (it != yCloud->idIndexMap.end()) {
904 jatomIndex = it->second;
906 std::cerr <<
"Your map must be wrong.\n";
916 angles[1] = acos(delta[2] / r);
917 angles[0] = atan2(delta[0], delta[1]);
927 for (
int m = 0; m < 2 * l + 1; m++) {
928 QlmTotal.
ptq[iatom].ylm[m] += yl[m];
934 QlmTotal.
ptq[iatom].ylm =
941 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
949 nnumNeighbours = nList[iatom].size() - 1;
953 for (
int j = 1; j <= nnumNeighbours; j++) {
955 dot_product = {0, 0};
959 jatomID = nList[iatom][j];
962 auto it = yCloud->idIndexMap.find(jatomID);
964 if (it != yCloud->idIndexMap.end()) {
965 jatomIndex = it->second;
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));
972 Inorm = Inorm + (qI * std::conj(qI));
973 Jnorm = Jnorm + (qJ * std::conj(qJ));
976 complexDenominator =
std::sqrt(Inorm * Jnorm);
977 complexCij = dot_product / complexDenominator;
979 cij_real = complexCij.
real();
986 q_value /= (double)nnumNeighbours;
988 resultQ[iatom] = q_value;
1016 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
1019 if ((*q6)[iatom] > 0.5) {
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;
1026 avgQ3 /= (double)nnumNeighbours;
1029 if (avgQ3 <= -0.75) {
1030 if (avgQ3 < -0.85) {
1058 int firstFrame,
bool isSlice,
std::string outputFileName) {
1059 int ih, ic, interIce,
water, unknown, clath, interClath,
1061 ih = ic =
water = unknown = interIce = total = 0;
1062 clath = interClath = 0;
1064 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
1068 if (yCloud->pts[iatom].inSlice ==
false) {
1098 outputFile.
open(path +
"bop/" + outputFileName, std::ios_base::app);
1101 if (yCloud->currentFrame == firstFrame) {
1102 outputFile <<
"Frame Ic Ih Interfacial Clath InterClath Water Total\n";
1105 outputFile << yCloud->currentFrame <<
" " << ic <<
" " << ih <<
" "
1106 << interIce <<
" " << clath <<
" " << interClath <<
" " <<
water
1107 <<
" " << total <<
"\n";
1129 int nnumNeighbours =
1130 nList[iatom].size() - 1;
1131 int neighStaggered =
1139 if (num_staggrd == 2) {
1141 for (
int j = 1; j <= nnumNeighbours; j++) {
1143 jatomID = nList[iatom][j];
1145 auto it = yCloud->idIndexMap.find(jatomID);
1147 if (it != yCloud->idIndexMap.end()) {
1148 jatomIndex = it->second;
1150 std::cerr <<
"Something is gravely wrong with your map.\n";
1155 if (neighStaggered > 2) {
1162 if (num_staggrd == 3 && num_eclipsd == 0) {
1164 for (
int j = 1; j <= nnumNeighbours; j++) {
1167 jatomID = nList[iatom][j];
1169 auto it = yCloud->idIndexMap.find(jatomID);
1171 if (it != yCloud->idIndexMap.end()) {
1172 jatomIndex = it->second;
1174 std::cerr <<
"Something is gravely wrong with your map.\n";
1179 if (neighStaggered == 2) {
1200 int num_staggrd = 0;
1203 int nnumNeighbours =
1204 nList[jatom].size() - 1;
1207 for (
int i = 0; i < nnumNeighbours; i++) {
1208 bondType = yCloud->pts[jatom].c_ij[i].classifier;
File for the bond order parameter analysis.
std::vector< std::complex< double > > lookupTableQ6Vec(std::array< double, 2 > angles)
Lookup table for Q6.
std::vector< double > getq6(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
std::vector< std::complex< double > > lookupTableQ3Vec(std::array< double, 2 > angles)
Lookup table for Q3.
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.
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.
std::vector< YlmAtom > ptq
molSys::PointCloud< molSys::Point< double >, double > getCorrel(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
std::array< double, 2 > radialCoord(std::array< double, 3 > cartCoord)
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.
std::complex< double > lookupTableQ6(int m, std::array< double, 2 > angles)
Lookup table for Q6 (m=0 to m=12)
bool isInterfacial(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int iatom, int num_staggrd, int num_eclipsd)
std::vector< std::complex< double > > spheriHarmo(int orderL, std::array< double, 2 > radialCoord)
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.
std::complex< double > lookupTableQ3(int m, std::array< double, 2 > angles)
Lookup table for Q3 (m=0 to m=6)
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.
molSys::PointCloud< molSys::Point< double >, double > reclassifyWater(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *q6)
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.
std::vector< std::complex< double > > avgVector(std::vector< std::complex< double >> v, int l, int neigh)
std::array< double, 3 > relDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
double c_value
Classifier according to CHILL, CHILL+ etc.
@ cubic
Ic, or particle type signifying Cubic Ice.
@ hexagonal
Ih, or particle type signifying Hexagonal Ice.
@ interClathrate
Interfacial clathrate ice phase.
@ interfacial
Interfacial ice: ice-like molecules which do not fulfill the strict criteria of the Ic or Ih phases.
@ reCubic
Reclassified as cubic ice, according to the order parameter.
@ clathrate
Clathrate ice phase.
@ water
Liquid/amorphous phase.
@ reHex
Reclassified as hexagonal ice, according to the order parameter.
@ unclassified
Not classified into any other category.
@ eclipsed
The bond is an eclipsed bond.
@ out_of_range
The bond cannot be classified as either staggered or eclipsed.
@ staggered
The bond is a staggered bond, according to the or value.
int makePath(const std::string &path)
This is the local orientational bond order parameter , of length .
This contains a collection of points; contains information for a particular frame.
This contains per-particle information.
This contains the bond classifier of enum type bond_type, and the bond correlation factor.