52 bool onlyTetrahedral) {
71 int numHC, numDDC, mixedRings, prismaticRings, basalRings;
74 bool doShapeMatching =
true;
90 if (onlyTetrahedral) {
97 atomTypes.
resize(yCloud->nop);
99 rmsdPerAtom.
resize(yCloud->nop);
104 for (
int ringSize = initRingSize; ringSize <= maxRingSize; ringSize++) {
110 if (ringsOneType.
size() == 0) {
116 ringsOneType.
size());
122 listHC =
ring::findHC(ringsOneType, &ringType, nList, &cageList);
125 listDDC =
ring::findDDC(ringsOneType, &ringType, listHC, &cageList);
135 &prismaticRings, &basalRings);
139 mixedRings, basalRings, prismaticRings,
204 int totalRingNum = rings.
size();
206 bool cond1, cond2, cond3;
213 notEquatorial.
resize(totalRingNum);
217 for (
int i = 0; i < listHC.
size(); i++) {
218 int currentRingIndex = listHC[i];
220 notEquatorial[currentRingIndex] =
true;
226 for (
int iring = 0; iring < totalRingNum; iring++) {
232 if (notEquatorial[iring]) {
238 peripheralRings.
clear();
243 if (cond1 ==
false) {
252 if (cond2 ==
false) {
264 if (cond3 ==
false) {
269 sort(peripheralRings.
begin(), peripheralRings.
end());
270 peripheralRings.
erase(
271 unique(peripheralRings.
begin(), peripheralRings.
end()),
272 peripheralRings.
end());
274 if (peripheralRings.
size() != 6) {
284 for (
int j = 0; j < peripheralRings.
size(); j++) {
285 jring = peripheralRings[j];
300 notEquatorial[jring] =
true;
333 int noOfCommonRings =
338 for (
int m = 0; m < 6; m++) {
339 index = rings[iring][m];
342 for (
int jring = 0; jring < rings.size(); jring++) {
343 if (iring == jring) {
348 for (
int k = 0; k < 6; k++) {
349 jElement = rings[jring][k];
350 if (jElement == index) {
352 (*peripheralRings).push_back(jring);
363 if (noOfCommonRings < 3) {
374 (*peripheralRings).
size();
375 sort((*peripheralRings).begin(), (*peripheralRings).end());
377 (*peripheralRings).begin() + numElements);
379 (*peripheralRings).resize(
std::distance((*peripheralRings).begin(), ip));
413 for (
int k = 0; k < ringSize; k++) {
416 for (
int i = k; i < k + 3; i++) {
428 for (
int m = 0; m < (*peripheralRings).size(); m++) {
429 jring = (*peripheralRings)[m];
448 (*peripheralRings).swap(newPeripherals);
452 if ((*peripheralRings).size() > 6) {
454 <<
"There are more than 6 peripheral rings. The code will fail. \n";
493 rings[(*peripheralRings)[2]],
494 rings[(*peripheralRings)[4]]);
504 rings[(*peripheralRings)[3]],
505 rings[(*peripheralRings)[5]]);
517 for (
int i = 0; i <= 3; i++) {
520 iring = (*peripheralRings)[i];
521 jring = (*peripheralRings)[i + 2];
525 if (common.
size() < 3) {
570 int totalRingNum = rings.
size();
582 isPrismatic.
resize(totalRingNum);
585 for (
int iring = 0; iring < totalRingNum - 1; iring++) {
588 if (isPrismatic[iring]) {
594 basal1 = rings[iring];
596 for (
int jring = iring + 1; jring < totalRingNum; jring++) {
599 if (isPrismatic[jring]) {
603 basal2 = rings[jring];
617 if (cond2 ==
false) {
626 }
else if ((*ringType)[iring] ==
ring::DDC) {
634 }
else if ((*ringType)[jring] ==
ring::DDC) {
639 prismaticRings.
clear();
643 for (
int k = 0; k < prismaticRings.
size(); k++) {
650 }
else if ((*ringType)[kring] ==
ring::DDC) {
656 isPrismatic[kring] =
true;
702 int l1 = (*basal1)[0];
703 int l2 = (*basal1)[1];
711 int compare1, compare2;
713 bool l1_neighbour, l2_neighbour;
715 bool isNeigh, notNeigh;
720 for (
int k = 0; k < ringSize; k++) {
722 l1_neighbour =
false;
723 l2_neighbour =
false;
728 auto it1 =
std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
730 if (it1 != nList[l1].end()) {
731 compare1 = (*basal1)[2];
732 compare2 = (*basal1)[4];
739 auto it2 =
std::find(nList[l2].begin() + 1, nList[l2].end(), m_k);
741 if (it2 != nList[l2].end()) {
742 compare1 = (*basal1)[3];
743 compare2 = (*basal1)[5];
755 if (l1_neighbour ==
false && l2_neighbour ==
false) {
762 for (
int k = 0; k <= 5; k++) {
763 currentKindex = kIndex + k;
765 if (currentKindex >= ringSize) {
766 currentKindex -= ringSize;
771 evenTriplet.
push_back((*basal2)[currentKindex]);
775 oddTriplet.
push_back((*basal2)[currentKindex]);
826 int needle1 = (*triplet)[1];
827 int needle2 = (*triplet)[2];
828 bool neighbourFound =
false;
829 bool neighOne =
false,
835 std::find(nList[atomOne].begin() + 1, nList[atomOne].end(), needle1);
836 if (it != nList[atomOne].end()) {
837 neighbourFound =
true;
842 it =
std::find(nList[atomTwo].begin() + 1, nList[atomTwo].end(), needle1);
843 if (it != nList[atomTwo].end()) {
844 neighbourFound =
true;
850 if (neighbourFound ==
false) {
859 it =
std::find(nList[atomTwo].begin() + 1, nList[atomTwo].end(), needle2);
861 if (it != nList[atomTwo].end()) {
873 it =
std::find(nList[atomOne].begin() + 1, nList[atomOne].end(), needle2);
875 if (it != nList[atomOne].end()) {
903 for (
int i = 0; i < (*triplet).size(); i++) {
904 iatom = (*triplet)[i];
907 for (
int j = 0; j < (*ring).size(); j++) {
911 it =
std::find(nList[jatom].begin() + 1, nList[jatom].end(), iatom);
913 if (it != nList[jatom].end()) {
940 int ringSize = rings[0].size();
946 for (
int i = 0; i < ringSize; i++) {
951 for (
int m = 0; m < 3; m++) {
953 if (iIndex >= ringSize) {
956 iTriplet.
push_back(rings[iring][iIndex]);
962 for (
int kring = 0; kring < rings.size(); kring++) {
964 if (kring == iring || kring == jring) {
973 if (common.
size() != 3) {
982 for (
int j = 0; j < ringSize; j++) {
983 int katom = rings[kring][j];
987 if (it == iTriplet.
end()) {
997 if (common.
size() == 3) {
999 (*listHC).push_back(kring);
1000 (*prismaticRings).push_back(kring);
1036 int dummyValue = -10;
1040 for (
int iring = 0; iring < (*ringType).size(); iring++) {
1048 std::sort((*listDDC).begin(), (*listDDC).end());
1049 auto iter =
std::find((*listDDC).begin(), (*listDDC).end(), iring);
1050 if (iter != (*listDDC).end()) {
1056 std::sort((*listHC).begin(), (*listHC).end());
1057 auto itr =
std::find((*listHC).begin(), (*listHC).end(), iring);
1058 if (itr != (*listHC).end()) {
1083 int ringSize = rings[0].size();
1086 for (
int iring = 0; iring < ringType.
size(); iring++) {
1124 for (
int i = 0; i < ringSize; i++) {
1125 iatom = rings[iring][i];
1131 if (ringSize == 6) {
1135 (*atomTypes)[iatom] = iRingType;
1138 (*atomTypes)[iatom] = iRingType;
1163 int *numDDC,
int *mixedRings,
int *prismaticRings,
1170 *prismaticRings = 0;
1175 for (
int icage = 0; icage < cageList.
size(); icage++) {
1190 for (
int iring = 0; iring < ringType.
size(); iring++) {
1201 *prismaticRings += 1;
1207 *prismaticRings += 1;
1247 int totalRingNum = rings.size();
1251 int ringSize = rings[0].
size();
1253 Eigen::MatrixXd refPointSet(ringSize, 3);
1260 for (
int iring = 0; iring < totalRingNum - 1; iring++) {
1263 basal1 = rings[iring];
1265 for (
int jring = iring + 1; jring < totalRingNum; jring++) {
1266 basal2 = rings[jring];
1273 if (cond1 ==
true) {
1287 &basal1, &basal2, rmsdPerAtom);
1346 int l1 = (*basal1)[0];
1361 for (
int k = 0; k < ringSize; k++) {
1362 l1_neighbour =
false;
1367 auto it =
std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
1370 if (it != nList[l1].end()) {
1371 l1_neighbour =
true;
1379 if (!l1_neighbour) {
1385 isNeighbour[kIndex] =
true;
1388 for (
int i = 1; i < ringSize; i++) {
1389 lAtomID = (*basal1)[i];
1390 for (
int k = 0; k < ringSize; k++) {
1392 if (isNeighbour[k]) {
1396 kAtomID = (*basal2)[k];
1401 std::find(nList[lAtomID].begin() + 1, nList[lAtomID].end(), kAtomID);
1404 if (it1 != nList[lAtomID].end()) {
1405 isNeighbour[k] =
true;
1413 for (
int k = 0; k < ringSize; k++) {
1415 if (!isNeighbour[k]) {
1436 bool isNeighbour =
false;
1443 for (
int l = 0; l < ringSize; l++) {
1447 for (
int m = 0; m < ringSize; m++) {
1450 auto it =
std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
1453 if (it != nList[l_k].end()) {
1475 int ringSize = basal1.
size();
1477 double infHugeNumber = 100000;
1478 double leastDist = infHugeNumber;
1486 for (
int m = 0; m < ringSize; m++) {
1493 if (leastDist > dist) {
1501 if (leastDist < heightCutoff) {
double periodicDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Inline generic function for obtaining the unwrapped periodic distance between two particles,...
bool notNeighboursOfRing(std::vector< std::vector< int >> nList, std::vector< int > *triplet, std::vector< int > *ring)
bool conditionOneDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
bool basalNeighbours(std::vector< std::vector< int >> nList, std::vector< int > *triplet, int atomOne, int atomTwo)
bool relaxedPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
std::vector< int > findDDC(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > listHC, std::vector< cage::Cage > *cageList)
bool conditionTwoDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings, int iring)
int getAtomTypesTopoBulk(std::vector< std::vector< int >> rings, std::vector< ring::strucType > ringType, std::vector< cage::iceType > *atomTypes)
int getStrucNumbers(std::vector< ring::strucType > ringType, std::vector< cage::Cage > cageList, int *numHC, int *numDDC, int *mixedRings, int *prismaticRings, int *basalRings)
Determines the number of HCs, DDCs, Mixed rings, prismatic and basal rings.
bool basalRingsSeparation(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, double heightCutoff=8)
Check to see that candidate basal prisms are not really far from each other.
bool basalConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
Tests whether two rings are basal rings (true) or not (false)
int findPrismatic(std::vector< std::vector< int >> rings, std::vector< int > *listHC, std::vector< strucType > *ringType, int iring, int jring, std::vector< int > *prismaticRings)
Finds the prismatic rings from basal rings iring and jring.
int topoBulkAnalysis(std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, bool onlyTetrahedral=true)
bool conditionThreeDDC(std::vector< std::vector< int >> rings, std::vector< int > *peripheralRings)
int findBulkPrisms(std::vector< std::vector< int >> rings, std::vector< ring::strucType > *ringType, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, double heightCutoff=8)
Find out which rings are prisms.
std::vector< int > findMixedRings(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< int > *listDDC, std::vector< int > *listHC)
bool basalPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
std::vector< int > findHC(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, std::vector< std::vector< int >> nList, std::vector< cage::Cage > *cageList)
bool commonElementsInThreeRings(std::vector< int > ring1, std::vector< int > ring2, std::vector< int > ring3)
Common elements in 3 rings.
bool findTripletInRing(std::vector< int > ring, std::vector< int > triplet)
Searches a particular ring for a triplet.
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
std::vector< std::vector< int > > getSingleRingSize(std::vector< std::vector< int >> rings, int ringSize)
Returns a vector of vectors of rings of a single size.
std::vector< int > findsCommonElements(std::vector< int > ring1, std::vector< int > ring2)
Returns the common elements of two rings.
@ bothBasal
The ring belongs to both a DDC and HC. It is a 'mixed' ring. The ring is also one of the basal rings ...
@ DDC
The ring belongs to a double-diamond cage (DDC).
@ HCbasal
The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an H...
@ unclassified
The ring is unclassified, which may be either water or a deformed type which cannot be classified by ...
@ bothPrismatic
The ring belongs to both a DDC and HC and is, thus, a 'mixed' ring. The ring is also one of the prism...
@ HCprismatic
The ring belongs only to a hexagonal cage (HC); specifically the ring is purely a prismatic ring of a...
@ Prism
The ring belongs to a prism block, classified according to the prism identification scheme.
bool matchUntetheredPrism(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, std::vector< double > *rmsdPerAtom)
Eigen::MatrixXd getPointSetRefRing(int n, int axialDim)
Topological network criteria functions.
int writeTopoBulkData(std::string path, int currentFrame, int numHC, int numDDC, int mixedRings, int basalRings, int prismaticRings, int firstFrame)
int writeLAMMPSdataTopoBulk(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< cage::iceType > atomTypes, std::string path, bool bondsBetweenDummy=false)
This contains a collection of points; contains information for a particular frame.
This contains per-particle information.
File containing functions used specific to bulk topological network critera.