49 int *atomID,
int firstFrame,
int currentFrame,
bool doShapeMatching) {
69 double avgPrismHeight = 2.845;
77 nDefPrismList.
resize(maxDepth - 2);
78 heightPercent.
resize(maxDepth - 2);
80 atomTypes.
resize(yCloud->nop, 1);
83 rmsdPerAtom.
resize(yCloud->nop, -1);
85 atomState.
resize(yCloud->nop);
89 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
96 if (ringsOneType.
size() == 0) {
97 nPrismList[ringSize - 3] = 0;
98 nDefPrismList[ringSize - 3] = 0;
99 heightPercent[ringSize - 3] = 0.0;
108 nImperfectPrisms = 0;
110 ringsOneType.
size());
115 &nImperfectPrisms, nList, yCloud, &rmsdPerAtom,
118 nPrismList[ringSize - 3] =
120 nDefPrismList[ringSize - 3] =
122 int totalPrisms = nPerfectPrisms + nImperfectPrisms;
124 heightPercent[ringSize - 3] =
127 if (nPerfectPrisms + nImperfectPrisms == 0) {
135 &atomTypes, &atomState);
143 yCloud->currentFrame, firstFrame);
145 if (doShapeMatching) {
154 if (doShapeMatching) {
195 int totalRingNum = rings.
size();
203 int ringSize = rings[0].
size();
204 *nImperfectPrisms = 0;
207 Eigen::MatrixXd refPointSet(ringSize, 3);
217 for (
int iring = 0; iring < totalRingNum - 1; iring++) {
220 basal1 = rings[iring];
222 for (
int jring = iring + 1; jring < totalRingNum; jring++) {
223 basal2 = rings[jring];
226 if (doShapeMatching ==
true || ringSize == 4) {
230 if (isAxialPair ==
false) {
247 if (cond2 ==
false) {
249 if (!doShapeMatching) {
257 if (relaxedCond ==
false) {
263 yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom,
false);
266 if (isDeformedPrism) {
269 *nImperfectPrisms += 1;
294 *nPerfectPrisms += 1;
311 if (doShapeMatching) {
313 yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom,
true);
327 sort(listPrism.
begin(), listPrism.
end());
350 int l1 = (*basal1)[0];
365 for (
int k = 0; k < ringSize; k++) {
366 l1_neighbour =
false;
371 auto it =
std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
374 if (it != nList[l1].end()) {
389 isNeighbour[kIndex] =
true;
392 for (
int i = 1; i < ringSize; i++) {
393 lAtomID = (*basal1)[i];
394 for (
int k = 0; k < ringSize; k++) {
396 if (isNeighbour[k]) {
400 kAtomID = (*basal2)[k];
405 std::find(nList[lAtomID].begin() + 1, nList[lAtomID].end(), kAtomID);
408 if (it1 != nList[lAtomID].end()) {
409 isNeighbour[k] =
true;
417 for (
int k = 0; k < ringSize; k++) {
419 if (!isNeighbour[k]) {
438 bool isNeighbour =
false;
445 for (
int l = 0; l < ringSize; l++) {
449 for (
int m = 0; m < ringSize; m++) {
452 auto it =
std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
455 if (it != nList[l_k].end()) {
494 bool axialBasal1, axialBasal2;
496 double areaXY, areaXZ,
517 jatomIndex = (*basal1)[0];
520 for (
int k = 1; k < ringSize; k++) {
521 iatomIndex = (*basal1)[k];
526 areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
527 (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
530 areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
531 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
534 areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
535 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
537 jatomIndex = iatomIndex;
541 iatomIndex = (*basal1)[0];
544 areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
545 (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
548 areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
549 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
552 areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
553 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
560 areaXY = fabs(areaXY);
561 areaXZ = fabs(areaXZ);
562 areaYZ = fabs(areaYZ);
569 if (areaYZ > areaXY && areaYZ > areaXZ) {
574 else if (axialDim == 1) {
575 if (areaXZ > areaXY && areaXZ > areaYZ) {
580 else if (axialDim == 2) {
581 if (areaXY > areaXZ && areaXY > areaYZ) {
586 std::cerr <<
"Could not find the axial dimension.\n";
597 jatomIndex = (*basal2)[0];
600 for (
int k = 1; k < ringSize; k++) {
601 iatomIndex = (*basal2)[k];
606 areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
607 (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
610 areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
611 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
614 areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
615 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
617 jatomIndex = iatomIndex;
621 iatomIndex = (*basal2)[0];
624 areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
625 (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
628 areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
629 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
632 areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
633 (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
640 areaXY = fabs(areaXY);
641 areaXZ = fabs(areaXZ);
642 areaYZ = fabs(areaYZ);
650 if (areaYZ > areaXY && areaYZ > areaXZ) {
655 else if (axialDim == 1) {
656 if (areaXZ > areaXY && areaXZ > areaYZ) {
661 else if (axialDim == 2) {
662 if (areaXY > areaXZ && areaXY > areaYZ) {
667 std::cerr <<
"Could not find the axial dimension.\n";
673 if (axialBasal1 ==
true && axialBasal2 ==
true) {
708 for (
int i = 0; i < listPrism.
size(); i++) {
709 iring = listPrism[i];
711 currentState = ringType[iring];
713 for (
int j = 0; j < ringSize; j++) {
714 iatom = rings[iring][j];
716 (*atomTypes)[iatom] = ringSize;
720 (*atomState)[iatom] = currentState;
723 if ((*atomState)[iatom] != currentState) {
726 (*atomState)[iatom] = currentState;
758 int nop = atomState.
size();
761 for (
int iatom = 0; iatom < nop; iatom++) {
765 (*atomTypes)[iatom] += maxDepth - 2;
768 (*atomTypes)[iatom] = 2;
789 int firstFrame,
int currentFrame) {
792 double shiftDistance;
793 double distFrmBoundary;
799 if (axialDim < 0 || axialDim > 2) {
803 double boxLowAxial = yCloud->boxLow[axialDim];
804 double boxHighAxial = boxLowAxial + yCloud->box[axialDim];
810 if (currentFrame == firstFrame) {
811 *atomID = yCloud->pts[0].atomID;
816 int iatomID = *atomID;
818 auto it = yCloud->idIndexMap.find(iatomID);
820 if (it != yCloud->idIndexMap.end()) {
821 atomIndex = it->second;
832 shiftDistance = boxLowAxial - yCloud->pts[atomIndex].x;
834 else if (axialDim == 1) {
835 shiftDistance = boxLowAxial - yCloud->pts[atomIndex].y;
838 shiftDistance = boxLowAxial - yCloud->pts[atomIndex].z;
842 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
845 yCloud->pts[iatom].x += shiftDistance;
847 if (yCloud->pts[iatom].x < boxLowAxial) {
848 distFrmBoundary = boxLowAxial - yCloud->pts[iatom].x;
849 yCloud->pts[iatom].x = boxHighAxial - distFrmBoundary;
852 else if (axialDim == 1) {
853 yCloud->pts[iatom].y += shiftDistance;
855 if (yCloud->pts[iatom].y < boxLowAxial) {
856 distFrmBoundary = boxLowAxial - yCloud->pts[iatom].y;
857 yCloud->pts[iatom].y = boxHighAxial - distFrmBoundary;
861 yCloud->pts[iatom].z += shiftDistance;
863 if (yCloud->pts[iatom].z < boxLowAxial) {
864 distFrmBoundary = boxLowAxial - yCloud->pts[iatom].z;
865 yCloud->pts[iatom].z = boxHighAxial - distFrmBoundary;
bool relaxedPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
int prismAnalysis(std::string path, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, int *atomID, int firstFrame, int currentFrame, bool doShapeMatching=false)
std::vector< int > findPrisms(std::vector< std::vector< int >> rings, std::vector< strucType > *ringType, int *nPerfectPrisms, int *nImperfectPrisms, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, bool doShapeMatching=false)
int deformedPrismTypes(std::vector< ring::strucType > atomState, std::vector< int > *atomTypes, int maxDepth)
Get the atom type values for deformed prisms.
int clearRingList(std::vector< std::vector< int >> &rings)
Erases memory for a vector of vectors for a list of rings.
bool discardExtraTetragonBlocks(std::vector< int > *basal1, std::vector< int > *basal2, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
int rmAxialTranslations(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int *atomID, int firstFrame, int currentFrame)
Shift the entire ice nanotube and remove axial translations.
int assignPrismType(std::vector< std::vector< int >> rings, std::vector< int > listPrism, int ringSize, std::vector< ring::strucType > ringType, std::vector< int > *atomTypes, std::vector< ring::strucType > *atomState)
bool hasCommonElements(std::vector< int > ring1, std::vector< int > ring2)
bool basalPrismConditions(std::vector< std::vector< int >> nList, std::vector< int > *basal1, std::vector< int > *basal2)
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.
@ unclassified
The ring is unclassified, which may be either water or a deformed type which cannot be classified by ...
@ Prism
The ring belongs to a prism block, classified according to the prism identification scheme.
bool matchPrism(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, bool isPerfect=true)
Eigen::MatrixXd getPointSetRefRing(int n, int axialDim)
int writePrismNum(std::string path, std::vector< int > nPrisms, std::vector< int > nDefPrisms, std::vector< double > heightPercent, int maxDepth, int currentFrame, int firstFrame)
int writeLAMMPSdumpINT(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > rmsdPerAtom, std::vector< int > atomTypes, int maxDepth, std::string path)
Write out a LAMMPS dump file containing the RMSD per atom.
int writeLAMMPSdataAllPrisms(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< int > atomTypes, int maxDepth, std::string path, bool doShapeMatching=false)
Write a data file for prisms of every type.
double normHeightPercent(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int nPrisms, double avgPrismHeight)
This contains a collection of points; contains information for a particular frame.
This contains per-particle information.
File containing functions used specific to quasi-one-dimensional topological network critera (the pri...