24 int ringSize = rings[0].
size();
26 bool padAtoms =
false;
33 if (rings.size() == 0) {
39 const char *path =
"../output";
49 total_size = rings.size() * ringSize;
53 for (
int iring = 0; iring < rings.size(); iring++) {
54 std::move(rings[iring].begin(), rings[iring].end(),
66 if (atoms.
size() != yCloud->nop) {
71 outputFile.
open(
"../output/" + filename);
101 outputFile <<
"Written out by D-SEAMS\n";
103 outputFile << atoms[atoms.
size() - 1] <<
" "
107 outputFile << bonds.size() <<
" bonds"
109 outputFile <<
"0 angles\n0 dihedrals\n0 impropers\n";
112 outputFile <<
"2 atom types\n";
114 outputFile <<
"1 atom types\n";
117 <<
"1 bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
119 outputFile <<
"0 " << yCloud->box[0] <<
" xlo xhi\n";
120 outputFile <<
"0 " << yCloud->box[1] <<
" ylo yhi\n";
121 outputFile <<
"0 " << yCloud->box[2] <<
" zlo zhi\n";
123 outputFile <<
"\nMasses\n\n";
124 outputFile <<
"1 15.999400 # O\n";
126 outputFile <<
"2 1.0 # dummy\n";
129 outputFile <<
"\nAtoms\n\n";
133 for (
int i = 0; i < atoms.
size(); i++) {
134 iatom = atoms[i] - 1;
138 if (atoms[i] != prevAtomID + 1) {
139 dummyAtoms = atoms[i] - prevAtomID - 1;
140 dummyID = prevAtomID;
142 for (
int j = 0; j < dummyAtoms; j++) {
146 outputFile << dummyID <<
" " << yCloud->pts[jatom].molID <<
" 2 0 "
147 << yCloud->pts[jatom].x <<
" " << yCloud->pts[jatom].y <<
" "
148 << yCloud->pts[jatom].z <<
"\n";
154 outputFile << atoms[i] <<
" " << yCloud->pts[iatom].molID <<
" 1 0 "
155 << yCloud->pts[iatom].x <<
" " << yCloud->pts[iatom].y <<
" "
156 << yCloud->pts[iatom].z <<
"\n";
158 prevAtomID = atoms[i];
162 outputFile <<
"\nBonds\n\n";
164 for (
int ibond = 0; ibond < bonds.size(); ibond++) {
165 outputFile << ibond + 1 <<
" 1 " << bonds[ibond][0] <<
" "
166 << bonds[ibond][1] <<
"\n";
183 const char *path =
"../output";
190 outputFile.
open(
"../output/" + filename);
195 for (
int iring = 0; iring < rings.size(); iring++) {
197 for (
int k = 0; k < rings[iring].size(); k++) {
198 outputFile << rings[iring][k] <<
" ";
224 const char *path =
"../output/prisms";
231 outputFile.
open(
"../output/prisms/" + filename);
237 for (
int iring = 0; iring < ringSize; iring++) {
238 iatomIndex = (*basal1)[iring];
240 outputFile << yCloud->pts[iatomIndex].x <<
" ";
241 outputFile << yCloud->pts[iatomIndex].y <<
" ";
242 outputFile << yCloud->pts[iatomIndex].z <<
" ";
248 for (
int iring = 0; iring < ringSize; iring++) {
249 iatomIndex = (*basal2)[iring];
251 outputFile << yCloud->pts[iatomIndex].x <<
" ";
252 outputFile << yCloud->pts[iatomIndex].y <<
" ";
253 outputFile << yCloud->pts[iatomIndex].z <<
" ";
262 outputFile.
open(
"../output/prisms/singleRing.dat");
264 for (
int iring = 0; iring < ringSize; iring++) {
265 iatomIndex = (*basal1)[iring];
267 outputFile << yCloud->pts[iatomIndex].x <<
" ";
268 outputFile << yCloud->pts[iatomIndex].y <<
" ";
269 outputFile << yCloud->pts[iatomIndex].z <<
" ";
291 int totalCages = (*cageList).size();
296 if (totalCages == 0) {
297 std::cerr <<
"There are no cages to print.\n";
307 for (
int icage = 0; icage < totalCages; icage++) {
308 type = (*cageList)[icage].type;
331 std::cerr <<
"The cage is of the wrong type\n";
359 strcpy(cageChar,
"../output/cages/hexCages");
360 actualCageType =
"hexCages";
362 strcpy(cageChar,
"../output/cages/doubleDiaCages");
363 actualCageType =
"doubleDiaCages";
366 std::cerr <<
"The cage is of the wrong type. Exit\n";
372 const char *path =
"../output/cages";
380 const fs::path path1(cageChar);
385 "../output/cages/" + actualCageType +
"/" + filename;
386 outputFile.
open(fileOutNameFull);
395 for (
int i = 0; i < 2; i++) {
396 iring = currentCage[i];
398 for (
int j = 0; j < ringSize; j++) {
399 iatomIndex = rings[iring][j] - 1;
401 outputFile << yCloud->pts[iatomIndex].x <<
" ";
402 outputFile << yCloud->pts[iatomIndex].y <<
" ";
403 outputFile << yCloud->pts[iatomIndex].z <<
" ";
412 for (
int i = 0; i < currentCage.
size(); i++) {
413 iring = currentCage[i];
415 for (
int j = 0; j < ringSize; j++) {
416 iatomIndex = rings[iring][j] - 1;
418 outputFile << yCloud->pts[iatomIndex].x <<
" ";
419 outputFile << yCloud->pts[iatomIndex].y <<
" ";
420 outputFile << yCloud->pts[iatomIndex].z <<
" ";
442 std::string filename =
"basalRings" + number +
".dat";
474 const char *path =
"../output/cages";
482 const fs::path path1(
"../output/cages/hexBasalRings");
486 std::string fileOutNameFull =
"../output/cages/hexBasalRings/" + filename;
487 outputFile.
open(fileOutNameFull);
492 outputFile <<
"# Particle IDs in the two basal rings\n";
496 basal1 = rings[currentCage[0]];
497 unOrderedBasal2 = rings[currentCage[1]];
499 for (
int i = 0; i < basal1.
size(); i++) {
504 for (
int k = 0; k < unOrderedBasal2.
size(); k++) {
512 for (
int n = 1; n < nList[iatom - 1].size(); n++) {
513 natom = nList[iatom - 1][n];
515 if (findAtom == natom) {
539 for (
int i = 0; i < 2; i++) {
540 if (basal2[i] != 0) {
548 for (
int i = 0; i < unOrderedBasal2.
size(); i++) {
549 if (unOrderedBasal2[i] == needle) {
557 nextElement = startHayStack + 2;
558 if (nextElement >= ringSize) {
559 nextElement -= ringSize;
566 if (basal2[startNeedle + 2] == unOrderedBasal2[nextElement]) {
569 for (
int i = 1; i < ringSize; i += 2) {
570 iatom = startNeedle + i;
571 jatom = startHayStack + i;
573 if (iatom >= ringSize) {
576 if (jatom >= ringSize) {
580 basal2[iatom] = unOrderedBasal2[jatom];
592 iatom = startNeedle + 2;
593 jatom = startHayStack - 2;
598 if (basal2[iatom] == unOrderedBasal2[jatom]) {
600 for (
int i = 1; i < ringSize; i += 2) {
601 iatom = startNeedle + i;
602 jatom = startHayStack - i;
604 if (iatom > ringSize) {
611 basal2[iatom] = unOrderedBasal2[jatom];
615 std::cerr <<
"Something is wrong with your HCs.\n";
626 for (
int i = 0; i < basal1.
size(); i++) {
627 outputFile << basal1[i] <<
" ";
632 for (
int i = 0; i < basal2.
size(); i++) {
633 outputFile << basal2[i] <<
" ";
653 std::string filename =
"basalRings" + number +
".dat";
660 int currentIatom, currentJatom;
671 const char *path =
"../output/deformed";
679 const fs::path path1(
"../output/deformed/basalRings");
683 std::string fileOutNameFull =
"../output/deformed/basalRings/" + filename;
684 outputFile.
open(fileOutNameFull);
687 const char *path =
"../output/perfect";
695 const fs::path path1(
"../output/perfect/basalRings");
699 std::string fileOutNameFull =
"../output/perfect/basalRings/" + filename;
700 outputFile.
open(fileOutNameFull);
706 outputFile <<
"# Particle IDs in the two basal rings\n";
714 for (
int l = 0; l < ringSize; l++) {
719 for (
int m = 0; m < ringSize; m++) {
723 auto it =
std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
726 if (it != nList[l_k].end()) {
741 std::cerr <<
"Something is wrong with your deformed prism.\n";
750 int tempJfor, tempJback;
752 tempJfor = jatom + 1;
753 tempJback = jatom - 1;
755 if (jatom == ringSize - 1) {
757 tempJback = ringSize - 2;
761 tempJback = ringSize - 1;
764 int forwardJ = (*basal2)[tempJfor];
765 int backwardJ = (*basal2)[tempJback];
766 int currentI = (*basal1)[iatom];
773 if (distClock < distAntiClock) {
777 if (distAntiClock < distClock) {
781 if (isClock ==
false && isAntiClock ==
false) {
787 for (
int i = 0; i < ringSize; i++) {
788 currentIatom = iatom + i;
789 if (currentIatom >= ringSize) {
790 currentIatom -= ringSize;
795 currentJatom = jatom + i;
796 if (currentJatom >= ringSize) {
797 currentJatom -= ringSize;
801 currentJatom = jatom - i;
802 if (currentJatom < 0) {
803 currentJatom += ringSize;
808 matchedBasal1.
push_back((*basal1)[currentIatom]);
809 matchedBasal2.
push_back((*basal2)[currentJatom]);
817 for (
int i = 0; i < matchedBasal1.
size(); i++) {
818 outputFile << matchedBasal1[i] <<
" ";
823 for (
int i = 0; i < matchedBasal2.
size(); i++) {
824 outputFile << matchedBasal2[i] <<
" ";
837 int largestCluster,
int numOfClusters,
838 int smallestCluster,
double avgClusterSize,
846 outputFile.
open(path +
"clusterStats.dat",
847 std::ios_base::app | std::ios_base::out);
855 if (currentFrame == firstFrame) {
856 outputFile <<
"Frame largestCluster numOfClusters smallestCluster "
861 outputFile << currentFrame <<
" " << largestCluster <<
" " << numOfClusters
862 <<
" " << smallestCluster <<
" " << avgClusterSize <<
"\n";
875 int currentFrame,
int firstFrame) {
885 outputFile.
open(path +
"topoINT/nPrisms.dat",
886 std::ios_base::app | std::ios_base::out);
890 if (currentFrame == firstFrame) {
891 outputFile <<
"Frame RingSize Num_of_prisms Height% RingSize ... Height\n";
898 outputFile << currentFrame <<
" ";
900 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
901 totalPrisms = nPrisms[ringSize - 3] + nDefPrisms[ringSize - 3];
903 outputFile << ringSize <<
" " << totalPrisms <<
" "
904 << nDefPrisms[ringSize - 3] <<
" " << heightPercent[ringSize - 3]
930 std::string outputDirName = path +
"topoMonolayer";
935 outputFileXY.
open(path +
"topoMonolayer/coverageAreaXY.dat",
936 std::ios_base::app | std::ios_base::out);
944 if (currentFrame == firstFrame) {
945 outputFileXY <<
"Frame RingSize Num_of_rings CoverageAreaXY% RingSize ... "
950 outputFileXY << currentFrame <<
" ";
952 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
953 outputFileXY << ringSize <<
" " << nRings[ringSize - 3] <<
" "
954 << coverageAreaXY[ringSize - 3] <<
" ";
957 outputFileXY <<
"\n";
959 outputFileXY.
close();
963 outputFileXZ.
open(path +
"topoMonolayer/coverageAreaXZ.dat",
964 std::ios_base::app | std::ios_base::out);
968 if (currentFrame == firstFrame) {
969 outputFileXZ <<
"Frame RingSize Num_of_rings CoverageAreaXZ% RingSize ... "
978 outputFileXZ << currentFrame <<
" ";
980 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
981 outputFileXZ << ringSize <<
" " << nRings[ringSize - 3] <<
" "
982 << coverageAreaXZ[ringSize - 3] <<
" ";
985 outputFileXZ <<
"\n";
987 outputFileXZ.
close();
991 outputFileYZ.
open(path +
"topoMonolayer/coverageAreaYZ.dat",
992 std::ios_base::app | std::ios_base::out);
996 if (currentFrame == firstFrame) {
997 outputFileYZ <<
"Frame RingSize Num_of_rings CoverageAreaYZ% RingSize ... "
1006 outputFileYZ << currentFrame <<
" ";
1008 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1009 outputFileYZ << ringSize <<
" " << nRings[ringSize - 3] <<
" "
1010 << coverageAreaYZ[ringSize - 3] <<
" ";
1013 outputFileYZ <<
"\n";
1015 outputFileYZ.
close();
1025 double binwidth,
int nbin) {
1031 outputFile.
open(fileName, std::ios_base::app | std::ios_base::out);
1034 for (
int ibin = 0; ibin < nbin; ibin++) {
1036 r = binwidth * (ibin + 0.5);
1037 outputFile << r <<
" " << (*rdfValues)[ibin] <<
"\n";
1050 int numDDC,
int mixedRings,
int basalRings,
1051 int prismaticRings,
int firstFrame) {
1061 outputFile.
open(path +
"bulkTopo/cageData.dat",
1062 std::ios_base::app | std::ios_base::out);
1069 if (currentFrame == firstFrame) {
1070 outputFile <<
"Frame HCnumber DDCnumber MixedRingNumber PrismaticRings "
1074 outputFile << currentFrame <<
" " << numHC <<
" " << numDDC <<
" "
1075 << mixedRings <<
" " << prismaticRings <<
" " << basalRings
1096 std::string outputDirName = path +
"bulkTopo/dumpFiles";
1100 if (yCloud->currentFrame == firstFrame) {
1101 outputFile.
open(path +
"bulkTopo/typeInfo.dat");
1102 outputFile <<
"Atom types in the dump files are:\n";
1103 outputFile <<
" Type 0 (dummy) = unidentified phase\n";
1104 outputFile <<
" Type 1 (hc) = atom belonging to a Hexagonal Cage.\n";
1105 outputFile <<
" Type 2 (ddc) = atom belonging to a Double-Diamond Cage\n";
1106 outputFile <<
" Type 3 (mixed) = atom belonging to a mixed ring shared by "
1109 <<
" Type 4 (pnc) = atom belonging to a pair of pentagonal rings\n";
1110 outputFile <<
" Type 5 (mixed2) = atom belonging to a pentagonal "
1111 "nanochannel, shared by DDCs/HCs\n";
1116 outputFile.
open(path +
"bulkTopo/dumpFiles/" + filename);
1134 outputFile <<
"ITEM: TIMESTEP\n";
1136 outputFile << yCloud->currentFrame <<
"\n";
1138 outputFile <<
"ITEM: NUMBER OF ATOMS\n";
1140 outputFile << yCloud->pts.size() <<
"\n";
1142 outputFile <<
"ITEM: BOX BOUNDS pp pp pp\n";
1144 outputFile << yCloud->boxLow[0] <<
" " << yCloud->boxLow[0] + yCloud->box[0]
1146 outputFile << yCloud->boxLow[1] <<
" " << yCloud->boxLow[1] + yCloud->box[1]
1148 outputFile << yCloud->boxLow[2] <<
" " << yCloud->boxLow[2] + yCloud->box[2]
1151 outputFile <<
"ITEM: ATOMS id mol type x y z rmsd\n";
1158 for (
int i = 0; i < yCloud->pts.size(); i++) {
1160 yCloud->pts[i].atomID;
1162 outputFile << iatom <<
" " << yCloud->pts[i].molID <<
" " << atomTypes[i]
1163 <<
" " << yCloud->pts[i].x <<
" " << yCloud->pts[i].y <<
" "
1164 << yCloud->pts[i].z <<
" " << rmsdPerAtom[i] <<
"\n";
1187 std::string outputDirName = path +
"topoINT/dumpFiles";
1191 outputFile.
open(path +
"topoINT/dumpFiles/" + filename);
1209 outputFile <<
"ITEM: TIMESTEP\n";
1211 outputFile << yCloud->currentFrame <<
"\n";
1213 outputFile <<
"ITEM: NUMBER OF ATOMS\n";
1215 outputFile << yCloud->pts.size() <<
"\n";
1217 outputFile <<
"ITEM: BOX BOUNDS pp pp pp\n";
1219 outputFile << yCloud->boxLow[0] <<
" " << yCloud->boxLow[0] + yCloud->box[0]
1221 outputFile << yCloud->boxLow[1] <<
" " << yCloud->boxLow[1] + yCloud->box[1]
1223 outputFile << yCloud->boxLow[2] <<
" " << yCloud->boxLow[2] + yCloud->box[2]
1226 outputFile <<
"ITEM: ATOMS id mol type x y z rmsd\n";
1233 for (
int i = 0; i < yCloud->pts.size(); i++) {
1235 yCloud->pts[i].atomID;
1237 outputFile << iatom <<
" " << yCloud->pts[i].molID <<
" " << atomTypes[i]
1238 <<
" " << yCloud->pts[i].x <<
" " << yCloud->pts[i].y <<
" "
1239 << yCloud->pts[i].z <<
" " << rmsdPerAtom[i] <<
"\n";
1254 int maxDepth,
std::string path,
bool doShapeMatching) {
1263 "system-prisms-" +
std::to_string(yCloud->currentFrame) +
".data";
1274 outputDirName = path +
"topoINT/dataFiles/";
1278 outputFile.
open(path +
"topoINT/dataFiles/" + filename);
1308 outputFile <<
"Written out by D-SEAMS\n";
1310 outputFile << yCloud->pts.size() <<
" "
1314 outputFile << bonds.
size() <<
" bonds"
1316 outputFile <<
"0 angles\n0 dihedrals\n0 impropers\n";
1318 if (doShapeMatching) {
1319 outputFile << 2 * maxDepth - 2 <<
" atom types\n";
1321 outputFile << maxDepth <<
" atom types\n";
1326 <<
" bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1328 outputFile << yCloud->boxLow[0] <<
" " << yCloud->boxLow[0] + yCloud->box[0]
1330 outputFile << yCloud->boxLow[1] <<
" " << yCloud->boxLow[1] + yCloud->box[1]
1332 outputFile << yCloud->boxLow[2] <<
" " << yCloud->boxLow[2] + yCloud->box[2]
1335 outputFile <<
"\nMasses\n\n";
1336 outputFile <<
"1 15.999400 # dummy\n";
1337 outputFile <<
"2 1.0 # mixedRings \n";
1339 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1340 outputFile << ringSize <<
" 15.999400 # prism" << ringSize <<
"\n";
1343 if (doShapeMatching) {
1344 for (
int ringSize = maxDepth + 1; ringSize <= 2 * maxDepth - 2;
1346 int p = ringSize - maxDepth + 2;
1347 outputFile << ringSize <<
" 15.999400 # deformPrism" << p <<
"\n";
1351 outputFile <<
"\nAtoms\n\n";
1355 for (
int i = 0; i < yCloud->pts.size(); i++) {
1357 yCloud->pts[i].atomID;
1360 outputFile << iatom <<
" " << yCloud->pts[i].molID <<
" " << atomTypes[i]
1361 <<
" 0 " << yCloud->pts[i].x <<
" " << yCloud->pts[i].y <<
" "
1362 << yCloud->pts[i].z <<
"\n";
1367 outputFile <<
"\nBonds\n\n";
1369 for (
int ibond = 0; ibond < bonds.
size(); ibond++) {
1371 outputFile << ibond + 1 <<
" 1 " << bonds[ibond][0] <<
" "
1372 << bonds[ibond][1] <<
"\n";
1397 "system-rings-" +
std::to_string(yCloud->currentFrame) +
".data";
1406 std::string outputDirName = path +
"topoMonolayer";
1408 outputDirName = path +
"topoMonolayer/dataFiles/";
1412 outputFile.
open(path +
"topoMonolayer/dataFiles/" + filename);
1443 outputFile <<
"Written out by D-SEAMS\n";
1445 outputFile << yCloud->pts.size() <<
" "
1449 outputFile << bonds.
size() <<
" bonds"
1451 outputFile <<
"0 angles\n0 dihedrals\n0 impropers\n";
1453 outputFile << maxDepth <<
" atom types\n";
1457 <<
" bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1459 outputFile << yCloud->boxLow[0] <<
" " << yCloud->boxLow[0] + yCloud->box[0]
1461 outputFile << yCloud->boxLow[1] <<
" " << yCloud->boxLow[1] + yCloud->box[1]
1463 outputFile << yCloud->boxLow[2] <<
" " << yCloud->boxLow[2] + yCloud->box[2]
1466 outputFile <<
"\nMasses\n\n";
1467 outputFile <<
"1 15.999400 # dummy\n";
1468 outputFile <<
"2 1.0 # \n";
1470 for (
int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1471 outputFile << ringSize <<
" 15.999400 # ring" << ringSize <<
"\n";
1474 outputFile <<
"\nAtoms\n\n";
1478 for (
int i = 0; i < yCloud->pts.size(); i++) {
1480 yCloud->pts[i].atomID;
1483 outputFile << iatom <<
" " << yCloud->pts[i].molID <<
" " << atomTypes[i]
1484 <<
" 0 " << yCloud->pts[i].x <<
" " << yCloud->pts[i].y <<
" "
1485 << yCloud->pts[i].z <<
"\n";
1490 outputFile <<
"\nBonds\n\n";
1492 for (
int ibond = 0; ibond < bonds.
size(); ibond++) {
1494 outputFile << ibond + 1 <<
" 1 " << bonds[ibond][0] <<
" "
1495 << bonds[ibond][1] <<
"\n";
1514 int ringSize = rings[0].
size();
1516 bool padAtoms =
false;
1525 bool atomOne, atomTwo;
1530 if (listPrism.
size() == 0) {
1547 const char *path =
"../output";
1555 size_t total_size{0};
1557 total_size = listPrism.
size() * ringSize;
1561 for (
int iring = 0; iring < listPrism.
size(); iring++) {
1562 std::move(rings[listPrism[iring]].begin(), rings[listPrism[iring]].end(),
1574 if (atoms.
size() != yCloud->nop) {
1580 outputFile.
open(
"../output/" + filename);
1610 outputFile <<
"Written out by D-SEAMS\n";
1612 outputFile << yCloud->pts.size() <<
" "
1616 outputFile << bonds.
size() <<
" bonds"
1618 outputFile <<
"0 angles\n0 dihedrals\n0 impropers\n";
1621 outputFile <<
"2 atom types\n";
1623 outputFile <<
"1 atom types\n";
1627 <<
" bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1629 outputFile <<
"0 " << yCloud->box[0] <<
" xlo xhi\n";
1630 outputFile <<
"0 " << yCloud->box[1] <<
" ylo yhi\n";
1631 outputFile <<
"0 " << yCloud->box[2] <<
" zlo zhi\n";
1633 outputFile <<
"\nMasses\n\n";
1634 outputFile <<
"1 15.999400 # O\n";
1636 outputFile <<
"2 1.0 # dummy\n";
1639 outputFile <<
"\nAtoms\n\n";
1643 for (
int i = 0; i < atoms.
size(); i++) {
1644 iatom = atoms[i] - 1;
1648 if (atoms[i] != prevAtomID + 1) {
1649 dummyAtoms = atoms[i] - prevAtomID - 1;
1650 dummyID = prevAtomID;
1652 for (
int j = 0; j < dummyAtoms; j++) {
1654 jatom = dummyID - 1;
1656 outputFile << dummyID <<
" " << yCloud->pts[jatom].molID <<
" 2 0 "
1657 << yCloud->pts[jatom].x <<
" " << yCloud->pts[jatom].y <<
" "
1658 << yCloud->pts[jatom].z <<
"\n";
1664 outputFile << atoms[i] <<
" " << yCloud->pts[iatom].molID <<
" 1 0 "
1665 << yCloud->pts[iatom].x <<
" " << yCloud->pts[iatom].y <<
" "
1666 << yCloud->pts[iatom].z <<
"\n";
1668 prevAtomID = atoms[i];
1672 if (atoms[atoms.
size() - 1] != yCloud->nop) {
1674 for (
int id = atoms[atoms.
size() - 1] + 1; id <= yCloud->nop;
id++) {
1676 outputFile <<
id <<
" " << yCloud->pts[jatom].molID <<
" 2 0 "
1677 << yCloud->pts[jatom].x <<
" " << yCloud->pts[jatom].y <<
" "
1678 << yCloud->pts[jatom].z <<
"\n";
1683 outputFile <<
"\nBonds\n\n";
1685 for (
int ibond = 0; ibond < bonds.
size(); ibond++) {
1687 isPrismBond =
false;
1693 if (it != atoms.
end()) {
1695 }
else if (bonds[ibond][0] == atoms[0]) {
1697 }
else if (bonds[ibond][0] == atoms[atoms.
size() - 1]) {
1702 if (it1 != atoms.
end()) {
1704 }
else if (bonds[ibond][1] == atoms[0]) {
1706 }
else if (bonds[ibond][1] == atoms[atoms.
size() - 1]) {
1710 if (atomOne ==
false || atomTwo ==
false) {
1711 isPrismBond =
false;
1718 outputFile << ibond + 1 <<
" 1 " << bonds[ibond][0] <<
" "
1719 << bonds[ibond][1] <<
"\n";
1722 outputFile << ibond + 1 <<
" 2 " << bonds[ibond][0] <<
" "
1723 << bonds[ibond][1] <<
"\n";
1743 int ringSize = rings[0].
size();
1745 bool padAtoms =
false;
1760 if ((*cageList).size() == 0) {
1765 if (numCages == 0) {
1775 const char *path =
"../output";
1783 size_t total_size{0};
1785 total_size = nRings * ringSize;
1790 for (
int icage = 0; icage < (*cageList).size(); icage++) {
1792 if ((*cageList)[icage].type != type) {
1796 for (
int k = 0; k < (*cageList)[icage].rings.size(); k++) {
1797 iring = (*cageList)[icage].rings[k];
1798 std::move(rings[iring].begin(), rings[iring].end(),
1812 if (atoms.
size() != yCloud->nop) {
1818 outputFile.
open(
"../output/" + filename);
1848 outputFile <<
"Written out by D-SEAMS\n";
1850 outputFile << yCloud->pts.size() <<
" "
1854 outputFile << bonds.
size() <<
" bonds"
1856 outputFile <<
"0 angles\n0 dihedrals\n0 impropers\n";
1859 outputFile <<
"2 atom types\n";
1861 outputFile <<
"1 atom types\n";
1865 <<
" bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1867 outputFile <<
"0 " << yCloud->box[0] <<
" xlo xhi\n";
1868 outputFile <<
"0 " << yCloud->box[1] <<
" ylo yhi\n";
1869 outputFile <<
"0 " << yCloud->box[2] <<
" zlo zhi\n";
1871 outputFile <<
"\nMasses\n\n";
1874 actualCageType =
"HC";
1876 actualCageType =
"DDC";
1878 actualCageType =
"error";
1881 outputFile <<
"1 15.999400 # " << actualCageType <<
"\n";
1883 outputFile <<
"2 1.0 # dummy\n";
1886 outputFile <<
"\nAtoms\n\n";
1890 for (
int i = 0; i < atoms.
size(); i++) {
1891 iatom = atoms[i] - 1;
1895 if (atoms[i] != prevAtomID + 1) {
1896 dummyAtoms = atoms[i] - prevAtomID - 1;
1897 dummyID = prevAtomID;
1899 for (
int j = 0; j < dummyAtoms; j++) {
1901 jatom = dummyID - 1;
1903 outputFile << dummyID <<
" " << yCloud->pts[jatom].molID <<
" 2 0 "
1904 << yCloud->pts[jatom].x <<
" " << yCloud->pts[jatom].y <<
" "
1905 << yCloud->pts[jatom].z <<
"\n";
1911 outputFile << atoms[i] <<
" " << yCloud->pts[iatom].molID <<
" 1 0 "
1912 << yCloud->pts[iatom].x <<
" " << yCloud->pts[iatom].y <<
" "
1913 << yCloud->pts[iatom].z <<
"\n";
1915 prevAtomID = atoms[i];
1919 if (atoms[atoms.
size() - 1] != yCloud->nop) {
1921 for (
int id = atoms[atoms.
size() - 1] + 1; id <= yCloud->nop;
id++) {
1923 outputFile <<
id <<
" " << yCloud->pts[jatom].molID <<
" 2 0 "
1924 << yCloud->pts[jatom].x <<
" " << yCloud->pts[jatom].y <<
" "
1925 << yCloud->pts[jatom].z <<
"\n";
1930 outputFile <<
"\nBonds\n\n";
1932 for (
int ibond = 0; ibond < bonds.
size(); ibond++) {
1934 outputFile << ibond + 1 <<
" 1 " << bonds[ibond][0] <<
" "
1935 << bonds[ibond][1] <<
"\n";
1957 outputFile.
open(path + outFile, std::ios_base::app);
1974 outputFile <<
"ITEM: TIMESTEP\n";
1975 outputFile << yCloud->currentFrame <<
"\n";
1976 outputFile <<
"ITEM: NUMBER OF ATOMS\n";
1977 outputFile << yCloud->nop <<
"\n";
1978 outputFile <<
"ITEM: BOX BOUNDS pp pp pp\n";
1979 for (
int k = 0; k < yCloud->boxLow.size(); k++) {
1980 outputFile << yCloud->boxLow[k] <<
" "
1981 << yCloud->boxLow[k] + yCloud->box[k];
1983 if (yCloud->box.size() == 2 * yCloud->boxLow.size()) {
1985 << yCloud->box[k + yCloud->boxLow.size()];
1990 outputFile <<
"ITEM: ATOMS id mol type x y z\n";
1993 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
1994 outputFile << yCloud->pts[iatom].atomID <<
" " << yCloud->pts[iatom].molID;
1998 outputFile <<
" Ic " << yCloud->pts[iatom].x <<
" "
1999 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2003 outputFile <<
" Ih " << yCloud->pts[iatom].x <<
" "
2004 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2008 outputFile <<
" wat " << yCloud->pts[iatom].x <<
" "
2009 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2013 outputFile <<
" intFc " << yCloud->pts[iatom].x <<
" "
2014 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2018 outputFile <<
" clathrate " << yCloud->pts[iatom].x <<
" "
2019 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2023 outputFile <<
" interClathrate " << yCloud->pts[iatom].x <<
" "
2024 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2028 outputFile <<
" unclassified " << yCloud->pts[iatom].x <<
" "
2029 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2033 outputFile <<
" reIc " << yCloud->pts[iatom].x <<
" "
2034 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2038 outputFile <<
" reIh " << yCloud->pts[iatom].x <<
" "
2039 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
2063 cijFile.
open(
"cij.txt", std::ofstream::out | std::ofstream::app);
2064 q3File.
open(
"q3.txt", std::ofstream::out | std::ofstream::app);
2065 q6File.
open(
"q6.txt", std::ofstream::out | std::ofstream::app);
2067 for (
int iatom = 0; iatom < yCloud->nop; iatom++) {
2068 if (yCloud->pts[iatom].type != 1) {
2072 nNumNeighbours = nList[iatom].size() - 1;
2074 for (
int j = 0; j < nNumNeighbours; j++) {
2075 cijFile << yCloud->pts[iatom].c_ij[j].c_value <<
"\n";
2076 avgQ3 += yCloud->pts[iatom].c_ij[j].c_value;
2078 avgQ3 /= nNumNeighbours;
2079 q3File << avgQ3 <<
"\n";
2080 q6File << avgQ6[iatom] <<
"\n";
2099 clusterFile.
open(fileName, std::ofstream::out | std::ofstream::app);
2102 clusterFile.
close();
2118 int currentAtomType;
2119 int numAtomTypes = 6;
2129 if (bondsBetweenDummy) {
2141 outputDirName = path +
"bulkTopo/dataFiles/";
2145 outputFile.
open(path +
"bulkTopo/dataFiles/" + filename);
2175 outputFile <<
"Written out by D-SEAMS\n";
2177 outputFile << yCloud->pts.size() <<
" "
2181 outputFile << bonds.
size() <<
" bonds"
2183 outputFile <<
"0 angles\n0 dihedrals\n0 impropers\n";
2185 outputFile << numAtomTypes <<
" atom types\n";
2189 <<
" bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
2191 outputFile << yCloud->boxLow[0] <<
" " << yCloud->boxLow[0] + yCloud->box[0]
2193 outputFile << yCloud->boxLow[1] <<
" " << yCloud->boxLow[1] + yCloud->box[1]
2195 outputFile << yCloud->boxLow[2] <<
" " << yCloud->boxLow[2] + yCloud->box[2]
2198 outputFile <<
"\nMasses\n\n";
2199 outputFile <<
"1 15.999400 # dummy\n";
2200 outputFile <<
"2 15.999400 # hc \n";
2201 outputFile <<
"3 15.999400 # ddc \n";
2202 outputFile <<
"4 15.999400 # mixed \n";
2203 outputFile <<
"5 15.999400 # pnc \n";
2204 outputFile <<
"6 15.999400 # pncHexaMixed \n";
2206 outputFile <<
"\nAtoms\n\n";
2210 for (
int i = 0; i < yCloud->pts.size(); i++) {
2212 yCloud->pts[i].atomID;
2217 currentAtomType = 2;
2220 currentAtomType = 3;
2224 currentAtomType = 4;
2228 currentAtomType = 5;
2232 currentAtomType = 6;
2236 currentAtomType = 1;
2241 outputFile << iatom <<
" " << yCloud->pts[i].molID <<
" " << currentAtomType
2242 <<
" 0 " << yCloud->pts[i].x <<
" " << yCloud->pts[i].y <<
" "
2243 << yCloud->pts[i].z <<
"\n";
2248 outputFile <<
"\nBonds\n\n";
2250 for (
int ibond = 0; ibond < bonds.
size(); ibond++) {
2252 outputFile << ibond + 1 <<
" 1 " << bonds[ibond][0] <<
" "
2253 << bonds[ibond][1] <<
"\n";
2270 int nAtoms = atoms.
size();
2277 outputDirName = path +
"bulkTopo/clusterXYZ/";
2279 outputDirName = path +
"bulkTopo/clusterXYZ/frame-" +
2284 outputFile.
open(path +
"bulkTopo/clusterXYZ/frame-" +
2292 outputFile << nAtoms <<
"\n";
2294 <<
"Generated by d-SEAMS. 0 type=hc and 1 type =ddc\n";
2297 for (
int i = 0; i < nAtoms; i++) {
2300 outputFile << type <<
" " << yCloud->pts[iatom].x <<
" "
2301 << yCloud->pts[iatom].y <<
" " << yCloud->pts[iatom].z <<
"\n";
T back_inserter(T... args)
std::vector< std::vector< int > > populateBonds(std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
std::vector< std::vector< int > > createBondsFromCages(std::vector< std::vector< int >> rings, std::vector< cage::Cage > *cageList, cage::cageType type, int *nRings)
int largestIceCluster(std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int >> nList, std::vector< bool > *isIce, std::vector< int > *clusterID, std::vector< int > *nClusters, std::unordered_map< int, int > *indexNumber, int firstFrame)
Finds the largest ice cluster.
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,...
@ 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.
@ unclassified
Not classified into any other category.
std::vector< std::vector< int > > readBonds(std::string filename)
Reads bonds into a vector of vectors from a file with a specific format.
int writeBasalRingsPrism(std::vector< int > *basal1, std::vector< int > *basal2, int prismNum, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool isDeformed)
Write out the basal rings for a particular prism.
int writeLAMMPSdata(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> bonds, std::string filename="system-rings.data")
Write a data file for rings.
int writeAllCages(std::string path, std::vector< cage::Cage > *cageList, std::vector< std::vector< int >> rings, std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int currentFrame)
int writePrismNum(std::string path, std::vector< int > nPrisms, std::vector< int > nDefPrisms, std::vector< double > heightPercent, int maxDepth, int currentFrame, int firstFrame)
int writeRingNum(std::string path, int currentFrame, std::vector< int > nRings, std::vector< double > coverageAreaXY, std::vector< double > coverageAreaXZ, std::vector< double > coverageAreaYZ, int maxDepth, int firstFrame)
int makePath(const std::string &path)
int writeLAMMPSdataCages(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> rings, std::vector< cage::Cage > *cageList, cage::cageType type, int numCages, std::string filename="system-cages.data")
int writeTopoBulkData(std::string path, int currentFrame, int numHC, int numDDC, int mixedRings, int basalRings, int prismaticRings, int firstFrame)
int writeDump(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::string outFile)
Generic function for writing out to a dump file.
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 writeBasalRingsHex(std::vector< int > currentCage, int cageNum, std::vector< std::vector< int >> nList, std::vector< std::vector< int >> rings)
Write out the basal rings of a particular Hexagonal cage.
int writePrisms(std::vector< int > *basal1, std::vector< int > *basal2, int prismNum, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Function for writing out each prism.
int writeLAMMPSdataAllRings(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< int > atomTypes, int maxDepth, std::string path)
Write a data file for rings of every type for a monolayer.
int printRDF(std::string fileName, std::vector< double > *rdfValues, double binwidth, int nbin)
Function for printing out the RDF, given the filename.
int writeLAMMPSdataPrisms(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> rings, bool useBondFile, std::string bondFile, std::vector< int > listPrism, std::vector< std::vector< int >> nList, std::string filename="system-prisms.data")
Write a data file for prisms of a single type.
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)
int writeRings(std::vector< std::vector< int >> rings, std::string filename="rings.dat")
Function for printing out ring info, when there is no volume slice.
int writeLAMMPSdumpCages(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > rmsdPerAtom, std::vector< int > atomTypes, std::string path, int firstFrame)
Write out a LAMMPS dump file containing the RMSD per atom for bulk ice.
int writeXYZcluster(std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > atoms, int clusterID, cage::cageType type)
Function for writing out the XYZ files for each cluster.
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.
int writeEachCage(std::vector< int > currentCage, int cageNum, cage::cageType type, std::vector< std::vector< int >> rings, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Write out a particular cage to a file.
int writeCluster(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string fileName="cluster.txt", bool isSlice=false, int largestIceCluster=0)
Function for printing the largest ice cluster.
int writeClusterStats(std::string path, int currentFrame, int largestCluster, int numOfClusters, int smallestCluster, double avgClusterSize, int firstFrame)
Function for writing out cluster statistics.
int writeHisto(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< double > avgQ6)
This contains a collection of points; contains information for a particular frame.
This contains per-particle information.