seams_output.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 <seams_input.hpp>
12 #include <seams_output.hpp>
13 
21  std::string filename) {
22  std::ofstream outputFile;
23  std::vector<int> atoms; // Holds all atom IDs to print
24  int ringSize = rings[0].size(); // Ring size of each ring in rings
25  int iatom; // Index, not atom ID
26  bool padAtoms = false; // Add extra atoms if the atom IDs are skipped
27  int prevAtomID = 0; // Check for previous atom ID
28  int dummyAtoms = 0; // Number of dummy atoms to fill
29  int dummyID;
30  int jatom; // Array index is 1 less than the ID (index for dummy atom)
31  // ----------------
32  // Return if there are no rings
33  if (rings.size() == 0) {
34  return 1;
35  }
36  // ----------------
37  // Otherwise create file
38  // Create output dir if it doesn't exist already
39  const char *path = "../output"; // relative to the build directory
40  fs::path dir(path);
41  // if (fs::create_directory(dir)) {
42  // std::cerr << "Output directory created\n";
43  // }
44  // ----------------
45  // Get all the unique atomIDs of the atoms in the rings of this type
46  // Put all atom IDs into one 1-D vector
47  size_t total_size{0};
48  // Get the total number of atoms (repeated)
49  total_size = rings.size() * ringSize;
50  // Reserve this size inside atoms
51  atoms.reserve(total_size);
52  // Fill up all these atom IDs
53  for (int iring = 0; iring < rings.size(); iring++) {
54  std::move(rings[iring].begin(), rings[iring].end(),
55  std::back_inserter(atoms));
56  } // end of loop through all rings in the list
57 
58  // Sort the array according to atom ID, which will be needed to get the
59  // unique IDs and to remove duplicates
60  sort(atoms.begin(), atoms.end());
61  // Get the unique atom IDs
62  auto ip = std::unique(atoms.begin(), atoms.end());
63  // Resize the vector to remove undefined terms
64  atoms.resize(std::distance(atoms.begin(), ip));
65  // If the number of atoms is less than the total nop, add dummy atoms
66  if (atoms.size() != yCloud->nop) {
67  padAtoms = true;
68  }
69  // ----------------
70  // Write output to file inside the output directory
71  outputFile.open("../output/" + filename);
72  // FORMAT:
73  // Comment Line
74  // 4 atoms
75  // 4 bonds
76  // 0 angles
77  // 0 dihedrals
78  // 0 impropers
79  // 1 atom types
80  // 1 bond types
81  // 0 angle types
82  // 0 dihedral types
83  // 0 improper types
84  // -1.124000 52.845002 xlo xhi
85  // 0.000000 54.528999 ylo yhi
86  // 1.830501 53.087501 zlo zhi
87 
88  // Masses
89 
90  // 1 15.999400 # O
91 
92  // Atoms
93 
94  // 1 1 1 0 20.239 1.298 6.873 # O
95  // 2 1 1 0 0 5.193 6.873 # O
96  // 3 1 1 0 2.249 1.298 6.873 # O
97 
98  // -------
99  // Write the header
100  // Write comment line
101  outputFile << "Written out by D-SEAMS\n";
102  // Write out the number of atoms
103  outputFile << atoms[atoms.size() - 1] << " "
104  << "atoms"
105  << "\n";
106  // Number of bonds
107  outputFile << bonds.size() << " bonds"
108  << "\n";
109  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
110  // If padded atoms are required, two atom types will be required
111  if (padAtoms) {
112  outputFile << "2 atom types\n";
113  } else {
114  outputFile << "1 atom types\n";
115  } // end of atom types
116  outputFile
117  << "1 bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
118  // Box lengths
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";
122  // Masses
123  outputFile << "\nMasses\n\n";
124  outputFile << "1 15.999400 # O\n";
125  if (padAtoms) {
126  outputFile << "2 1.0 # dummy\n";
127  }
128  // Atoms
129  outputFile << "\nAtoms\n\n";
130  // -------
131  // Write out the atom coordinates
132  // Loop through atoms
133  for (int i = 0; i < atoms.size(); i++) {
134  iatom = atoms[i] - 1; // The actual index is one less than the ID
135  // -----------
136  // Pad out
137  // Fill in dummy atoms if some have been skipped
138  if (atoms[i] != prevAtomID + 1) {
139  dummyAtoms = atoms[i] - prevAtomID - 1;
140  dummyID = prevAtomID;
141  // Loop to write out dummy atoms
142  for (int j = 0; j < dummyAtoms; j++) {
143  dummyID++;
144  jatom = dummyID - 1;
145  // 1 molecule-tag atom-type q x y z
146  outputFile << dummyID << " " << yCloud->pts[jatom].molID << " 2 0 "
147  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
148  << yCloud->pts[jatom].z << "\n";
149  } // end of dummy atom write-out
150  } // end of check for dummy atom printing
151  // -----------
152  // Write out coordinates
153  // 1 molecule-tag atom-type q x y z
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";
157  // update the previous atom ID
158  prevAtomID = atoms[i];
159  } // end of loop through all atoms in atomID
160 
161  // Print the bonds now!
162  outputFile << "\nBonds\n\n";
163  // Loop through all bonds
164  for (int ibond = 0; ibond < bonds.size(); ibond++) {
165  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
166  << bonds[ibond][1] << "\n";
167  }
168 
169  // Once the datafile has been printed, exit
170  return 0;
171 }
172 
179  std::string filename) {
180  std::ofstream outputFile;
181  // ----------------
182  // Create output dir if it doesn't exist already
183  const char *path = "../output"; // relative to the build directory
184  fs::path dir(path);
185  // if (fs::create_directory(dir)) {
186  // std::cerr << "Output directory created\n";
187  // }
188  // ----------------
189  // Write output to file inside the output directory
190  outputFile.open("../output/" + filename);
191 
192  // Format:
193  // 272 214 906 1361 388 1
194 
195  for (int iring = 0; iring < rings.size(); iring++) {
196  // Otherwise, write out to the file
197  for (int k = 0; k < rings[iring].size(); k++) {
198  outputFile << rings[iring][k] << " ";
199  } // end of loop through ring elements
200  outputFile << "\n";
201  } // end of loop through rings
202 
203  outputFile.close();
204 
205  return 0;
206 }
207 
213  std::vector<int> *basal1, std::vector<int> *basal2, int prismNum,
214  molSys::PointCloud<molSys::Point<double>, double> *yCloud) {
215  std::ofstream outputFile;
216  std::string number = std::to_string(prismNum);
217  std::string filename = "prism" + number + ".dat";
218  int ringSize =
219  (*basal1).size(); // Size of the ring; each ring contains n elements
220  int iatomIndex; // index of atom coordinate being written out
221 
222  // ----------------
223  // Create output dir if it doesn't exist already
224  const char *path = "../output/prisms"; // relative to the build directory
225  fs::path dir(path);
226  // if (fs::create_directory(dir)) {
227  // std::cerr << "Output Prism directory created\n";
228  // }
229  // ----------------
230  // Write output to file inside the output directory
231  outputFile.open("../output/prisms/" + filename);
232 
233  // Format:
234  // x y z coordinates of each node
235 
236  // For basal 1
237  for (int iring = 0; iring < ringSize; iring++) {
238  iatomIndex = (*basal1)[iring]; // C++ indices are one less
239  // Write the coordinates out to the file
240  outputFile << yCloud->pts[iatomIndex].x << " ";
241  outputFile << yCloud->pts[iatomIndex].y << " ";
242  outputFile << yCloud->pts[iatomIndex].z << " ";
243 
244  outputFile << "\n";
245  } // end of loop through basal1
246 
247  // For basal 2
248  for (int iring = 0; iring < ringSize; iring++) {
249  iatomIndex = (*basal2)[iring]; // C++ indices are one less
250  // Write the coordinates out to the file
251  outputFile << yCloud->pts[iatomIndex].x << " ";
252  outputFile << yCloud->pts[iatomIndex].y << " ";
253  outputFile << yCloud->pts[iatomIndex].z << " ";
254 
255  outputFile << "\n";
256  } // end of loop through basal1
257 
258  outputFile.close();
259 
260  // ---- Print out all the coordinates of a single ring, for prismNum=1 only
261  if (prismNum == 1) {
262  outputFile.open("../output/prisms/singleRing.dat");
263  // For basal 1
264  for (int iring = 0; iring < ringSize; iring++) {
265  iatomIndex = (*basal1)[iring]; // C++ indices are one less
266  // Write the coordinates out to the file
267  outputFile << yCloud->pts[iatomIndex].x << " ";
268  outputFile << yCloud->pts[iatomIndex].y << " ";
269  outputFile << yCloud->pts[iatomIndex].z << " ";
270 
271  outputFile << "\n";
272  } // end of loop through basal1
273  outputFile.close();
274  }
275 
276  return 0;
277 }
278 
284  std::string path, std::vector<cage::Cage> *cageList,
286  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
287  int currentFrame) {
288  int numDDC; // Number of DDCs
289  int numHC; // Number of HCs
290  int numMC; // Number of MCs
291  int totalCages = (*cageList).size(); // Total number of cages
292  cage::cageType type; // Current cage type
293 
294  // ---------------------------------
295  // Error handling!
296  if (totalCages == 0) {
297  std::cerr << "There are no cages to print.\n";
298  return 1;
299  }
300  // ---------------------------------
301  // Init
302  numDDC = 0;
303  numHC = 0;
304  numMC = 0;
305 
306  // Loop through every cage
307  for (int icage = 0; icage < totalCages; icage++) {
308  type = (*cageList)[icage].type;
309  // ------
310  // Add to the cage type and write out to the appropriate folders
311  // Hexagonal Cages
312  if (type == cage::HexC) {
313  numHC++;
314  sout::writeEachCage((*cageList)[icage].rings, numHC, type, rings, yCloud);
315  sout::writeBasalRingsHex((*cageList)[icage].rings, numHC, nList, rings);
316  } // end of write out of HCs
317  // Double diamond Cages
318  else if (type == cage::DoubleDiaC) {
319  numDDC++;
320  sout::writeEachCage((*cageList)[icage].rings, numDDC, type, rings,
321  yCloud);
322  } // end of write out of DDCs
323  // // Mixed Cages
324  // else if (type == cage::Mixed) {
325  // numMC++;
326  // sout::writeEachCage((*cageList)[icage].rings, numMC, type, rings,
327  // yCloud);
328  // } // end of write out of MCs
329  // // Error
330  else {
331  std::cerr << "The cage is of the wrong type\n";
332  continue;
333  } // some error
334  // ------
335  } // end of loop through all cages
336 
337  return 0;
338 }
339 
345  std::vector<int> currentCage, int cageNum, cage::cageType type,
347  molSys::PointCloud<molSys::Point<double>, double> *yCloud) {
348  std::ofstream outputFile;
349  std::string number = std::to_string(cageNum);
350  std::string filename = "cage" + number + ".dat";
351  int ringSize =
352  rings[0].size(); // Size of the ring; each ring contains n elements
353  int iatomIndex; // index of atom coordinate being written out
354  std::string actualCageType; // is icage a DDC, HC or MC?
355  char cageChar[100]; // is icage a DDC, HC or MC?
356  int iring; // Ring index of the current ring
357 
358  if (type == cage::HexC) {
359  strcpy(cageChar, "../output/cages/hexCages");
360  actualCageType = "hexCages";
361  } else if (type == cage::DoubleDiaC) {
362  strcpy(cageChar, "../output/cages/doubleDiaCages");
363  actualCageType = "doubleDiaCages";
364  } else {
365  // throw error
366  std::cerr << "The cage is of the wrong type. Exit\n";
367  return 1;
368  } //
369 
370  // ----------------
371  // Create output dir if it doesn't exist already
372  const char *path = "../output/cages"; // relative to the build directory
373  fs::path dir(path);
374  // if (fs::create_directory(dir)) {
375  // std::cerr << "Output Cage directory created\n";
376  // }
377  // ----------------
378  // Subdirectory
379 
380  const fs::path path1(cageChar);
381  // fs::create_directories(path1);
382 
383  // Write output to file inside the output directory
384  std::string fileOutNameFull =
385  "../output/cages/" + actualCageType + "/" + filename;
386  outputFile.open(fileOutNameFull);
387 
388  // Format:
389  // x y z coordinates of each node
390 
391  // For hexagonal cages:
392  if (type == cage::HexC) {
393  // Print out only basal ring atoms, since they describe the outer structure
394  // The first two rings are basal rings
395  for (int i = 0; i < 2; i++) {
396  iring = currentCage[i]; // Current iring
397  // Get every node of iring
398  for (int j = 0; j < ringSize; j++) {
399  iatomIndex = rings[iring][j] - 1; // C++ indices are one less
400  // Write out the coordinates to the file
401  outputFile << yCloud->pts[iatomIndex].x << " ";
402  outputFile << yCloud->pts[iatomIndex].y << " ";
403  outputFile << yCloud->pts[iatomIndex].z << " ";
404 
405  outputFile << "\n";
406  } // end of loop through iring
407  } // end of getting every ring in the current cage
408  } // end of printing basal ring atoms for hexagonal cages
409  // For currentCage
410  else {
411  // Loop through all cages (could be duplicates) TODO: remove duplicates
412  for (int i = 0; i < currentCage.size(); i++) {
413  iring = currentCage[i]; // Current iring
414  // Get every node of iring
415  for (int j = 0; j < ringSize; j++) {
416  iatomIndex = rings[iring][j] - 1; // C++ indices are one less
417  // Write out the coordinates to the file
418  outputFile << yCloud->pts[iatomIndex].x << " ";
419  outputFile << yCloud->pts[iatomIndex].y << " ";
420  outputFile << yCloud->pts[iatomIndex].z << " ";
421 
422  outputFile << "\n";
423  } // end of loop through iring
424  } // end of getting every ring in the current cage
425  } // end of cage printing (has duplicates)
426 
427  // Close the output file
428  outputFile.close();
429 
430  return 0;
431 }
432 
437 int sout::writeBasalRingsHex(std::vector<int> currentCage, int cageNum,
439  std::vector<std::vector<int>> rings) {
440  std::ofstream outputFile;
441  std::string number = std::to_string(cageNum);
442  std::string filename = "basalRings" + number + ".dat";
443  int ringSize =
444  rings[0].size(); // Size of the ring; each ring contains n elements
445  int iatomIndex; // index of atom coordinate being written out TODO get rid of
446  // this
447  int iring; // Ring index of the current ring
448  // Variables for the hydrogen-bonded 'second' basal ring
450  basal2; // Elements are bonded to each element of basal1 in order
451  std::vector<int> basal1; // 'First' basal ring
453  unOrderedBasal2; // Unordered basal2 ring, the ith element is not
454  // necessarily bonded to the ith element of basal1
455  int iatom, jatom; // Atom numbers (starting from 1), not C++ indices; saved
456  // inside rings
457  int natom; // Neighbour list ID for iatom
458  int findAtom; // Atom number of atomID to find in neighbour list
459  bool atomFound; // bool to check if an atomID has been found in the neighbour
460  // list or not
461  // Variables for ordering basal2
462  // After finding the nearest neighbours, elements which are not nearest
463  // neighbours are assigned a value of zero.
464  int needle; // First non-zero element of basal2 after getting the nearest
465  // neighbours
466  int startNeedle; // Index of basal2; the first non-zero element of basal2
467  int startHayStack; // Index of unOrderedBasal2, corresponding to the first
468  // non-zero element of basal2
469  bool isClock; // Original order of unOrderedBasal2
470  int nextElement; // Next element in unOrderedBasal2 after startHayStack
471 
472  // ----------------
473  // Create output dir if it doesn't exist already
474  const char *path = "../output/cages"; // relative to the build directory
475  fs::path dir(path);
476  // if (fs::create_directory(dir)) {
477  // std::cerr << "Output Cage directory created\n";
478  // }
479  // ----------------
480  // Subdirectory
481 
482  const fs::path path1("../output/cages/hexBasalRings");
483  // fs::create_directories(path1);
484 
485  // Write output to file inside the output directory
486  std::string fileOutNameFull = "../output/cages/hexBasalRings/" + filename;
487  outputFile.open(fileOutNameFull);
488 
489  // Format:
490  // Coordinate IDs (starting from 1), ordered according to the input XYZ file
491  // The first line is a comment line
492  outputFile << "# Particle IDs in the two basal rings\n";
493 
494  // ---------------
495  // Find the nearest neighbours of basal1 elements in basal2
496  basal1 = rings[currentCage[0]]; // First basal ring
497  unOrderedBasal2 = rings[currentCage[1]]; // Unordered second basal ring
498 
499  for (int i = 0; i < basal1.size(); i++) {
500  iatom = basal1[i]; // This is the atom particle ID, not the C++ index
501 
502  // Search through unOrderedBasal2 for an element in the neighbourlist of
503  // iatom
504  for (int k = 0; k < unOrderedBasal2.size(); k++) {
505  findAtom =
506  unOrderedBasal2[k]; // Atom ID to find in the neighbour list of iatom
507 
508  atomFound = false; // init
509  jatom = 0;
510 
511  // Search through the neighbour list for findAtom
512  for (int n = 1; n < nList[iatom - 1].size(); n++) {
513  natom = nList[iatom - 1][n]; // Atom ID
514 
515  if (findAtom == natom) {
516  atomFound = true;
517  break;
518  } // Check
519  } // Loop through nList
520 
521  if (atomFound) {
522  jatom = natom;
523  break;
524  } // atom has been found
525  } // end of loop through all atomIDs in unOrderedBasal2
526  basal2.push_back(jatom);
527  } // end of loop through all the atomIDs in basal1
528 
529  // ---------------------------------------------------
530  // Get particles which are not nearest neighbours
531  // 'Alternately' ordered particles
532 
533  // ---------------
534  // Init
535  isClock = false;
536 
537  // ---------------
538  // Get the first non-zero index {needle}
539  for (int i = 0; i < 2; i++) {
540  if (basal2[i] != 0) {
541  needle = basal2[i]; // Set the needle to the first non-zero element
542  startNeedle = i; // Index of basal2
543  break; // Break out of the loop
544  } // end of checking for non-zero index
545  } // end of getting the first non-zero index
546 
547  // Find the index matching needle in unOrderedBasal2
548  for (int i = 0; i < unOrderedBasal2.size(); i++) {
549  if (unOrderedBasal2[i] == needle) {
550  startHayStack = i; // Index at which needle has been found
551  } // end of check for needle
552  } // end of search for needle in unOrderedBasal2
553 
554  // ---------------
555  // Check for 'clockwise' order
556  // Check 'next' element
557  nextElement = startHayStack + 2;
558  if (nextElement >= ringSize) {
559  nextElement -= ringSize;
560  }
561 
562  // Init (indices of basal2 and unOrderedBasal2 respectively)
563  iatom = 0;
564  jatom = 0;
565 
566  if (basal2[startNeedle + 2] == unOrderedBasal2[nextElement]) {
567  isClock = true;
568  // Fill the three elements in increasing order
569  for (int i = 1; i < ringSize; i += 2) {
570  iatom = startNeedle + i;
571  jatom = startHayStack + i;
572  // Make sure the indices are not larger than 6
573  if (iatom >= ringSize) {
574  iatom -= ringSize;
575  }
576  if (jatom >= ringSize) {
577  jatom -= ringSize;
578  }
579 
580  basal2[iatom] = unOrderedBasal2[jatom];
581  } // end of filling next two alternate elements
582  } // check to see if clockwise order is correct
583 
584  // ---------------
585  // Check for 'anticlockwise' order
586  // Check 'next' element
587  if (!isClock) {
588  iatom = 0;
589  jatom = 0; // init
590 
591  // First element
592  iatom = startNeedle + 2;
593  jatom = startHayStack - 2;
594  if (jatom < 0) {
595  jatom += ringSize;
596  }
597 
598  if (basal2[iatom] == unOrderedBasal2[jatom]) {
599  // Fill the three elements in increasing order
600  for (int i = 1; i < ringSize; i += 2) {
601  iatom = startNeedle + i;
602  jatom = startHayStack - i;
603  // Make sure the indices are not larger than 6
604  if (iatom > ringSize) {
605  iatom -= ringSize;
606  }
607  if (jatom < 0) {
608  jatom += ringSize;
609  }
610 
611  basal2[iatom] = unOrderedBasal2[jatom];
612  } // end of filling next two alternate elements
613  } // check to see if anticlockwise order is correct
614  else {
615  std::cerr << "Something is wrong with your HCs.\n";
616  return 1;
617  } // exit with error
618  } // End of check for anticlockwise stuff
619 
620  // ---------------------------------------------------
621  // Print out the ordered rings
622  // For hexagonal cages:
623  // Only print out basal1 and basal2
624 
625  // BASAL1
626  for (int i = 0; i < basal1.size(); i++) {
627  outputFile << basal1[i] << " ";
628  } // end of loop through basal1
629  outputFile << "\n";
630 
631  // BASAL2
632  for (int i = 0; i < basal2.size(); i++) {
633  outputFile << basal2[i] << " ";
634  } // end of loop through basal2
635 
636  // Close the output file
637  outputFile.close();
638 
639  return 0;
640 }
641 
647  std::vector<int> *basal1, std::vector<int> *basal2, int prismNum,
649  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
650  bool isDeformed) {
651  std::ofstream outputFile;
652  std::string number = std::to_string(prismNum);
653  std::string filename = "basalRings" + number + ".dat";
654  int ringSize =
655  (*basal1).size(); // Size of the ring; each ring contains n elements
656  int nBonds; // Number of bonds between the two deformed prisms
657  int l_k, m_k; // Atom ID in basal1 and basal2 respectively
658  int iatom,
659  jatom; // Index not ID of basal1 and basal2 respectively which match
660  int currentIatom, currentJatom; // Index not ID for matchedBasal1 and
661  // matchedBasal2 respectively
662  bool isNeighbour; // Basal rings should have at least one nearest neighbour
663  bool isClock; // The order is in the original 'direction' of basal2
664  bool isAntiClock; // The order must be reversed
665  std::vector<int> matchedBasal1, matchedBasal2; // Vectors with revised order
666 
667  // ----------------
668  // Path for deformed prisms
669  if (isDeformed) {
670  // Create output dir if it doesn't exist already
671  const char *path = "../output/deformed"; // relative to the build directory
672  fs::path dir(path);
673  // if (fs::create_directory(dir)) {
674  // std::cerr << "Output Cage directory created\n";
675  // }
676  // ----------------
677  // Subdirectory
678 
679  const fs::path path1("../output/deformed/basalRings");
680  // fs::create_directories(path1);
681 
682  // Write output to file inside the output directory
683  std::string fileOutNameFull = "../output/deformed/basalRings/" + filename;
684  outputFile.open(fileOutNameFull);
685  } else {
686  // Create output dir if it doesn't exist already
687  const char *path = "../output/perfect"; // relative to the build directory
688  fs::path dir(path);
689  // if (fs::create_directory(dir)) {
690  // std::cerr << "Output Cage directory created\n";
691  // }
692  // ----------------
693  // Subdirectory
694 
695  const fs::path path1("../output/perfect/basalRings");
696  // fs::create_directories(path1);
697 
698  // Write output to file inside the output directory
699  std::string fileOutNameFull = "../output/perfect/basalRings/" + filename;
700  outputFile.open(fileOutNameFull);
701  } // end of creating file paths
702 
703  // Format:
704  // Coordinate IDs (starting from 1), ordered according to the input XYZ file
705  // The first line is a comment line
706  outputFile << "# Particle IDs in the two basal rings\n";
707 
708  // ---------------
709  // Find the nearest neighbours of basal1 elements in basal2
710 
711  nBonds = 0;
712  isNeighbour = false;
713  // Loop through every element of basal1
714  for (int l = 0; l < ringSize; l++) {
715  l_k = (*basal1)[l]; // This is the atom particle ID, not the C++ index
716 
717  // Search for the nearest neighbour of l_k in basal2
718  // Loop through basal2 elements
719  for (int m = 0; m < ringSize; m++) {
720  m_k = (*basal2)[m]; // Atom ID to find in the neighbour list of iatom
721 
722  // Find m_k inside l_k neighbour list
723  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
724 
725  // If the element has been found, for l1
726  if (it != nList[l_k].end()) {
727  isNeighbour = true;
728  iatom = l; // index of basal1
729  jatom = m; // index of basal2
730  break;
731  } // found element
732 
733  } // end of loop through all atomIDs in basal2
734 
735  if (isNeighbour) {
736  break;
737  } // nearest neighbour found
738  } // end of loop through all the atomIDs in basal1
739 
740  if (!isNeighbour) {
741  std::cerr << "Something is wrong with your deformed prism.\n";
742  return 1;
743  }
744  // ---------------------------------------------------
745  // Find out if the order of basal2 is 'clockwise' or 'anticlockwise'
746  isClock = false; // init
747  isAntiClock = false;
748 
749  // atom index (not ID)
750  int tempJfor, tempJback;
751 
752  tempJfor = jatom + 1;
753  tempJback = jatom - 1;
754 
755  if (jatom == ringSize - 1) {
756  tempJfor = 0;
757  tempJback = ringSize - 2;
758  }
759  if (jatom == 0) {
760  tempJfor = 1;
761  tempJback = ringSize - 1;
762  }
763 
764  int forwardJ = (*basal2)[tempJfor];
765  int backwardJ = (*basal2)[tempJback];
766  int currentI = (*basal1)[iatom];
767 
768  // Check clockwise
769  double distClock = gen::periodicDist(yCloud, currentI, forwardJ);
770  double distAntiClock = gen::periodicDist(yCloud, currentI, backwardJ);
771 
772  // Clockwise
773  if (distClock < distAntiClock) {
774  isClock = true;
775  } // end of clockwise check
776  // Anti-clockwise
777  if (distAntiClock < distClock) {
778  isAntiClock = true;
779  } // end of anti-clockwise check
780  // Some error
781  if (isClock == false && isAntiClock == false) {
782  // std::cerr << "The points are equidistant.\n";
783  return 1;
784  } // end of error handling
785  // ---------------------------------------------------
786  // Get the order of basal1 and basal2
787  for (int i = 0; i < ringSize; i++) {
788  currentIatom = iatom + i;
789  if (currentIatom >= ringSize) {
790  currentIatom -= ringSize;
791  } // end of basal1 element wrap-around
792 
793  // In clockwise order
794  if (isClock) {
795  currentJatom = jatom + i;
796  if (currentJatom >= ringSize) {
797  currentJatom -= ringSize;
798  } // wrap around
799  } // end of clockwise update
800  else {
801  currentJatom = jatom - i;
802  if (currentJatom < 0) {
803  currentJatom += ringSize;
804  } // wrap around
805  } // end of anti-clockwise update
806 
807  // Add to matchedBasal1 and matchedBasal2 now
808  matchedBasal1.push_back((*basal1)[currentIatom]);
809  matchedBasal2.push_back((*basal2)[currentJatom]);
810  }
811  // ---------------------------------------------------
812  // Print out the ordered rings
813  // For hexagonal cages:
814  // Only print out basal1 and basal2
815 
816  // BASAL1
817  for (int i = 0; i < matchedBasal1.size(); i++) {
818  outputFile << matchedBasal1[i] << " ";
819  } // end of loop through basal1
820  outputFile << "\n";
821 
822  // BASAL2
823  for (int i = 0; i < matchedBasal2.size(); i++) {
824  outputFile << matchedBasal2[i] << " ";
825  } // end of loop through basal2
826 
827  // Close the output file
828  outputFile.close();
829 
830  return 0;
831 }
832 
836 int sout::writeClusterStats(std::string path, int currentFrame,
837  int largestCluster, int numOfClusters,
838  int smallestCluster, double avgClusterSize,
839  int firstFrame) {
840  std::ofstream outputFile;
841  // ----------------
842  // Make the output directory if it doesn't exist
843  sout::makePath(path);
844  // ----------------
845  // Write output to file inside the output directory
846  outputFile.open(path + "clusterStats.dat",
847  std::ios_base::app | std::ios_base::out);
848 
849  // Format:
850  // Comment line
851  // 1 3 0 0 4 35 40 ....
852 
853  // ----------------
854  // Comment line for the first frame
855  if (currentFrame == firstFrame) {
856  outputFile << "Frame largestCluster numOfClusters smallestCluster "
857  "avgClusterSize\n";
858  }
859  // ----------------
860 
861  outputFile << currentFrame << " " << largestCluster << " " << numOfClusters
862  << " " << smallestCluster << " " << avgClusterSize << "\n";
863 
864  outputFile.close();
865 
866  return 0;
867 }
868 
873  std::vector<int> nDefPrisms,
874  std::vector<double> heightPercent, int maxDepth,
875  int currentFrame, int firstFrame) {
876  std::ofstream outputFile;
877  int totalPrisms; // Number of total prisms
878  // ----------------
879  // Make the output directory if it doesn't exist
880  sout::makePath(path);
881  std::string outputDirName = path + "topoINT";
882  sout::makePath(outputDirName);
883  // ----------------
884  // Write output to file inside the output directory
885  outputFile.open(path + "topoINT/nPrisms.dat",
886  std::ios_base::app | std::ios_base::out);
887 
888  // ----------------
889  // Write the comment line if the first frame is being written out
890  if (currentFrame == firstFrame) {
891  outputFile << "Frame RingSize Num_of_prisms Height% RingSize ... Height\n";
892  }
893  // ----------------
894  // Format:
895  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
896  // 1 3 0 0 4 35 40 ....
897 
898  outputFile << currentFrame << " ";
899 
900  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
901  totalPrisms = nPrisms[ringSize - 3] + nDefPrisms[ringSize - 3];
902  // Write out
903  outputFile << ringSize << " " << totalPrisms << " "
904  << nDefPrisms[ringSize - 3] << " " << heightPercent[ringSize - 3]
905  << " ";
906  }
907 
908  outputFile << "\n";
909 
910  outputFile.close();
911 
912  return 0;
913 }
914 
918 int sout::writeRingNum(std::string path, int currentFrame,
919  std::vector<int> nRings,
920  std::vector<double> coverageAreaXY,
921  std::vector<double> coverageAreaXZ,
922  std::vector<double> coverageAreaYZ, int maxDepth,
923  int firstFrame) {
924  std::ofstream outputFileXY;
925  std::ofstream outputFileXZ;
926  std::ofstream outputFileYZ;
927  // ----------------
928  // Make the output directory if it doesn't exist
929  sout::makePath(path);
930  std::string outputDirName = path + "topoMonolayer";
931  sout::makePath(outputDirName);
932  // ----------------
933  // Coverage Area of XY
934  // Write output to file inside the output directory
935  outputFileXY.open(path + "topoMonolayer/coverageAreaXY.dat",
936  std::ios_base::app | std::ios_base::out);
937 
938  // Format:
939  // Comment line
940  // 1 3 0 0 4 35 40 ....
941 
942  // ----------------
943  // Add comment for the first frame
944  if (currentFrame == firstFrame) {
945  outputFileXY << "Frame RingSize Num_of_rings CoverageAreaXY% RingSize ... "
946  "CoverageAreaXY%\n";
947  }
948  // ----------------
949 
950  outputFileXY << currentFrame << " ";
951 
952  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
953  outputFileXY << ringSize << " " << nRings[ringSize - 3] << " "
954  << coverageAreaXY[ringSize - 3] << " ";
955  }
956 
957  outputFileXY << "\n";
958 
959  outputFileXY.close();
960  // ----------------
961  // Coverage Area of XZ
962  // Write output to file inside the output directory
963  outputFileXZ.open(path + "topoMonolayer/coverageAreaXZ.dat",
964  std::ios_base::app | std::ios_base::out);
965 
966  // ----------------
967  // Add comment for the first frame
968  if (currentFrame == firstFrame) {
969  outputFileXZ << "Frame RingSize Num_of_rings CoverageAreaXZ% RingSize ... "
970  "CoverageAreaXZ%\n";
971  }
972  // ----------------
973 
974  // Format:
975  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
976  // 1 3 0 0 4 35 40 ....
977 
978  outputFileXZ << currentFrame << " ";
979 
980  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
981  outputFileXZ << ringSize << " " << nRings[ringSize - 3] << " "
982  << coverageAreaXZ[ringSize - 3] << " ";
983  }
984 
985  outputFileXZ << "\n";
986 
987  outputFileXZ.close();
988  // ----------------
989  // Coverage Area of YZ
990  // Write output to file inside the output directory
991  outputFileYZ.open(path + "topoMonolayer/coverageAreaYZ.dat",
992  std::ios_base::app | std::ios_base::out);
993 
994  // ----------------
995  // Add comment for the first frame
996  if (currentFrame == firstFrame) {
997  outputFileYZ << "Frame RingSize Num_of_rings CoverageAreaYZ% RingSize ... "
998  "CoverageAreaYZ%\n";
999  }
1000  // ----------------
1001 
1002  // Format:
1003  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
1004  // 1 3 0 0 4 35 40 ....
1005 
1006  outputFileYZ << currentFrame << " ";
1007 
1008  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1009  outputFileYZ << ringSize << " " << nRings[ringSize - 3] << " "
1010  << coverageAreaYZ[ringSize - 3] << " ";
1011  }
1012 
1013  outputFileYZ << "\n";
1014 
1015  outputFileYZ.close();
1016 
1017  return 0;
1018 }
1019 
1025  double binwidth, int nbin) {
1026  //
1027  std::ofstream outputFile; // For the output file
1028  double r; // Distance for the current bin
1029 
1030  // Append to the file
1031  outputFile.open(fileName, std::ios_base::app | std::ios_base::out);
1032 
1033  // Loop through all the bins
1034  for (int ibin = 0; ibin < nbin; ibin++) {
1035  //
1036  r = binwidth * (ibin + 0.5); // Current distance for ibin
1037  outputFile << r << " " << (*rdfValues)[ibin] << "\n";
1038  } // end of loop through all bins
1039 
1040  outputFile.close();
1041 
1042  return 0;
1043 }
1044 
1049 int sout::writeTopoBulkData(std::string path, int currentFrame, int numHC,
1050  int numDDC, int mixedRings, int basalRings,
1051  int prismaticRings, int firstFrame) {
1052  //
1053  std::ofstream outputFile;
1054  // ----------------
1055  // Make the output directory if it doesn't exist
1056  sout::makePath(path);
1057  std::string outputDirName = path + "bulkTopo";
1058  sout::makePath(outputDirName);
1059  // ----------------
1060  // Write output to file inside the output directory
1061  outputFile.open(path + "bulkTopo/cageData.dat",
1062  std::ios_base::app | std::ios_base::out);
1063 
1064  // Format:
1065  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
1066  // 1 3 0 0 4 35 40 ....
1067  // -------------------
1068  // If first frame then write the comment line
1069  if (currentFrame == firstFrame) {
1070  outputFile << "Frame HCnumber DDCnumber MixedRingNumber PrismaticRings "
1071  "basalRings\n";
1072  }
1073  // -------------------
1074  outputFile << currentFrame << " " << numHC << " " << numDDC << " "
1075  << mixedRings << " " << prismaticRings << " " << basalRings
1076  << "\n";
1077 
1078  outputFile.close();
1079  return 0;
1080 } // end of function
1081 
1087  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1088  std::vector<double> rmsdPerAtom, std::vector<int> atomTypes,
1089  std::string path, int firstFrame) {
1090  std::ofstream outputFile;
1091  int iatom; // Index, not atom ID
1092  std::string filename =
1093  "dump-" + std::to_string(yCloud->currentFrame) + ".lammpstrj";
1094  // ----------------
1095  // Make the output directory if it doesn't exist
1096  std::string outputDirName = path + "bulkTopo/dumpFiles";
1097  sout::makePath(outputDirName);
1098  // ----------------
1099  // Write out information about the data types
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 "
1107  "a DDC and HC\n";
1108  outputFile
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";
1112  outputFile.close();
1113  } // end of writing out information
1114  // ----------------
1115  // Write output to file inside the output directory
1116  outputFile.open(path + "bulkTopo/dumpFiles/" + filename);
1117  // ----------------------------------------------------
1118  // Header Format
1119 
1120  // ITEM: TIMESTEP
1121  // 0
1122  // ITEM: NUMBER OF ATOMS
1123  // 500
1124  // ITEM: BOX BOUNDS pp pp pp
1125  // -9.0400100000000005e-01 1.7170999999999999e+01
1126  // -9.0400100000000005e-01 1.7170999999999999e+01
1127  // -9.0400100000000005e-01 1.7170999999999999e+01
1128  // ITEM: ATOMS id mol type x y z rmsd
1129 
1130  // -----------------
1131  // -------
1132  // Write the header
1133  // ITEM: TIMESTEP
1134  outputFile << "ITEM: TIMESTEP\n";
1135  // Write out frame number
1136  outputFile << yCloud->currentFrame << "\n";
1137  // ITEM: NUMBER OF ATOMS
1138  outputFile << "ITEM: NUMBER OF ATOMS\n";
1139  // Number of atoms
1140  outputFile << yCloud->pts.size() << "\n";
1141  // ITEM: BOX BOUNDS pp pp pp
1142  outputFile << "ITEM: BOX BOUNDS pp pp pp\n";
1143  // Box lengths
1144  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1145  << "\n";
1146  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1147  << "\n";
1148  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1149  << "\n";
1150  // ITEM: ATOMS id mol type x y z rmsd
1151  outputFile << "ITEM: ATOMS id mol type x y z rmsd\n";
1152  // -------
1153  // Write out the atom coordinates
1154  // Format
1155  // ITEM: ATOMS id mol type x y z rmsd
1156  //
1157  // Loop through atoms
1158  for (int i = 0; i < yCloud->pts.size(); i++) {
1159  iatom =
1160  yCloud->pts[i].atomID; // The actual ID can be different from the index
1161  // Write out coordinates
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";
1165 
1166  } // end of loop through all atoms in pointCloud
1167  // -----------------------------------------------------
1168  outputFile.close(); // Close the file
1169  return 0;
1170 } // end of function
1171 
1177  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1178  std::vector<double> rmsdPerAtom, std::vector<int> atomTypes, int maxDepth,
1179  std::string path) {
1180  //
1181  std::ofstream outputFile;
1182  int iatom; // Index, not atom ID
1183  std::string filename =
1184  "dump-" + std::to_string(yCloud->currentFrame) + ".lammpstrj";
1185  // ----------------
1186  // Make the output directory if it doesn't exist
1187  std::string outputDirName = path + "topoINT/dumpFiles";
1188  sout::makePath(outputDirName);
1189  // ----------------
1190  // Write output to file inside the output directory
1191  outputFile.open(path + "topoINT/dumpFiles/" + filename);
1192  // ----------------------------------------------------
1193  // Header Format
1194 
1195  // ITEM: TIMESTEP
1196  // 0
1197  // ITEM: NUMBER OF ATOMS
1198  // 500
1199  // ITEM: BOX BOUNDS pp pp pp
1200  // -9.0400100000000005e-01 1.7170999999999999e+01
1201  // -9.0400100000000005e-01 1.7170999999999999e+01
1202  // -9.0400100000000005e-01 1.7170999999999999e+01
1203  // ITEM: ATOMS id mol type x y z rmsd
1204 
1205  // -----------------
1206  // -------
1207  // Write the header
1208  // ITEM: TIMESTEP
1209  outputFile << "ITEM: TIMESTEP\n";
1210  // Write out frame number
1211  outputFile << yCloud->currentFrame << "\n";
1212  // ITEM: NUMBER OF ATOMS
1213  outputFile << "ITEM: NUMBER OF ATOMS\n";
1214  // Number of atoms
1215  outputFile << yCloud->pts.size() << "\n";
1216  // ITEM: BOX BOUNDS pp pp pp
1217  outputFile << "ITEM: BOX BOUNDS pp pp pp\n";
1218  // Box lengths
1219  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1220  << "\n";
1221  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1222  << "\n";
1223  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1224  << "\n";
1225  // ITEM: ATOMS id mol type x y z rmsd
1226  outputFile << "ITEM: ATOMS id mol type x y z rmsd\n";
1227  // -------
1228  // Write out the atom coordinates
1229  // Format
1230  // ITEM: ATOMS id mol type x y z rmsd
1231  //
1232  // Loop through atoms
1233  for (int i = 0; i < yCloud->pts.size(); i++) {
1234  iatom =
1235  yCloud->pts[i].atomID; // The actual ID can be different from the index
1236  // Write out coordinates
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";
1240 
1241  } // end of loop through all atoms in pointCloud
1242  // -----------------------------------------------------
1243  return 0;
1244 } // end of function
1245 
1252  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1253  std::vector<std::vector<int>> nList, std::vector<int> atomTypes,
1254  int maxDepth, std::string path, bool doShapeMatching) {
1255  //
1256  std::ofstream outputFile;
1257  int iatom; // Index, not atom ID
1258  int bondTypes = 1;
1259  // Bond stuff
1260  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1261  // containing the atom IDs of each bond
1262  std::string filename =
1263  "system-prisms-" + std::to_string(yCloud->currentFrame) + ".data";
1264 
1265  // ---------------
1266  // Get the bonds
1267  bonds = bond::populateBonds(nList, yCloud);
1268  //
1269  // ----------------
1270  // Make the output directory if it doesn't exist
1271  sout::makePath(path);
1272  std::string outputDirName = path + "topoINT";
1273  sout::makePath(outputDirName);
1274  outputDirName = path + "topoINT/dataFiles/";
1275  sout::makePath(outputDirName);
1276  // ----------------
1277  // Write output to file inside the output directory
1278  outputFile.open(path + "topoINT/dataFiles/" + filename);
1279  // FORMAT:
1280  // Comment Line
1281  // 4 atoms
1282  // 4 bonds
1283  // 0 angles
1284  // 0 dihedrals
1285  // 0 impropers
1286  // 1 atom types
1287  // 1 bond types
1288  // 0 angle types
1289  // 0 dihedral types
1290  // 0 improper types
1291  // -1.124000 52.845002 xlo xhi
1292  // 0.000000 54.528999 ylo yhi
1293  // 1.830501 53.087501 zlo zhi
1294 
1295  // Masses
1296 
1297  // 1 15.999400 # O
1298 
1299  // Atoms
1300 
1301  // 1 1 1 0 20.239 1.298 6.873 # O
1302  // 2 1 1 0 0 5.193 6.873 # O
1303  // 3 1 1 0 2.249 1.298 6.873 # O
1304 
1305  // -------
1306  // Write the header
1307  // Write comment line
1308  outputFile << "Written out by D-SEAMS\n";
1309  // Write out the number of atoms
1310  outputFile << yCloud->pts.size() << " "
1311  << "atoms"
1312  << "\n";
1313  // Number of bonds
1314  outputFile << bonds.size() << " bonds"
1315  << "\n";
1316  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1317  // There are maxDepth-2 total types of prisms + dummy
1318  if (doShapeMatching) {
1319  outputFile << 2 * maxDepth - 2 << " atom types\n";
1320  } else {
1321  outputFile << maxDepth << " atom types\n";
1322  }
1323  // Bond types
1324  outputFile
1325  << bondTypes
1326  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1327  // Box lengths
1328  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1329  << " xlo xhi\n";
1330  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1331  << " ylo yhi\n";
1332  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1333  << " zlo zhi\n";
1334  // Masses
1335  outputFile << "\nMasses\n\n";
1336  outputFile << "1 15.999400 # dummy\n";
1337  outputFile << "2 1.0 # mixedRings \n";
1338  // There are maxDepth-2 other prism types
1339  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1340  outputFile << ringSize << " 15.999400 # prism" << ringSize << "\n";
1341  } // end of writing out perfect atom types
1342  // Write out the types for the deformed prism blocks
1343  if (doShapeMatching) {
1344  for (int ringSize = maxDepth + 1; ringSize <= 2 * maxDepth - 2;
1345  ringSize++) {
1346  int p = ringSize - maxDepth + 2;
1347  outputFile << ringSize << " 15.999400 # deformPrism" << p << "\n";
1348  } // end of writing out perfect atom types
1349  } // Deformed prism types
1350  // Atoms
1351  outputFile << "\nAtoms\n\n";
1352  // -------
1353  // Write out the atom coordinates
1354  // Loop through atoms
1355  for (int i = 0; i < yCloud->pts.size(); i++) {
1356  iatom =
1357  yCloud->pts[i].atomID; // The actual ID can be different from the index
1358  // Write out coordinates
1359  // atomID molecule-tag atom-type q x y z
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";
1363 
1364  } // end of loop through all atoms in pointCloud
1365 
1366  // Print the bonds now!
1367  outputFile << "\nBonds\n\n";
1368  // Loop through all bonds
1369  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1370  //
1371  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
1372  << bonds[ibond][1] << "\n";
1373  } // end of for loop for bonds
1374 
1375  outputFile.close();
1376  // Once the datafile has been printed, exit
1377  return 0;
1378 }
1379 
1386  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1387  std::vector<std::vector<int>> nList, std::vector<int> atomTypes,
1388  int maxDepth, std::string path) {
1389  //
1390  std::ofstream outputFile;
1391  int iatom; // Index, not atom ID
1392  int bondTypes = 1;
1393  // Bond stuff
1394  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1395  // containing the atom IDs of each bond
1396  std::string filename =
1397  "system-rings-" + std::to_string(yCloud->currentFrame) + ".data";
1398 
1399  // ---------------
1400  // Get the bonds
1401  bonds = bond::populateBonds(nList, yCloud);
1402  //
1403  // ----------------
1404  // Make the output directory if it doesn't exist
1405  sout::makePath(path);
1406  std::string outputDirName = path + "topoMonolayer";
1407  sout::makePath(outputDirName);
1408  outputDirName = path + "topoMonolayer/dataFiles/";
1409  sout::makePath(outputDirName);
1410 
1411  // Write output to file inside the output directory
1412  outputFile.open(path + "topoMonolayer/dataFiles/" + filename);
1413 
1414  // FORMAT:
1415  // Comment Line
1416  // 4 atoms
1417  // 4 bonds
1418  // 0 angles
1419  // 0 dihedrals
1420  // 0 impropers
1421  // 1 atom types
1422  // 1 bond types
1423  // 0 angle types
1424  // 0 dihedral types
1425  // 0 improper types
1426  // -1.124000 52.845002 xlo xhi
1427  // 0.000000 54.528999 ylo yhi
1428  // 1.830501 53.087501 zlo zhi
1429 
1430  // Masses
1431 
1432  // 1 15.999400 # O
1433 
1434  // Atoms
1435 
1436  // 1 1 1 0 20.239 1.298 6.873 # O
1437  // 2 1 1 0 0 5.193 6.873 # O
1438  // 3 1 1 0 2.249 1.298 6.873 # O
1439 
1440  // -------
1441  // Write the header
1442  // Write comment line
1443  outputFile << "Written out by D-SEAMS\n";
1444  // Write out the number of atoms
1445  outputFile << yCloud->pts.size() << " "
1446  << "atoms"
1447  << "\n";
1448  // Number of bonds
1449  outputFile << bonds.size() << " bonds"
1450  << "\n";
1451  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1452  // There are maxDepth-2 total types of prisms + dummy
1453  outputFile << maxDepth << " atom types\n";
1454  // Bond types
1455  outputFile
1456  << bondTypes
1457  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1458  // Box lengths
1459  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1460  << " xlo xhi\n";
1461  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1462  << " ylo yhi\n";
1463  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1464  << " zlo zhi\n";
1465  // Masses
1466  outputFile << "\nMasses\n\n";
1467  outputFile << "1 15.999400 # dummy\n";
1468  outputFile << "2 1.0 # \n";
1469  // There are maxDepth-2 other prism types
1470  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1471  outputFile << ringSize << " 15.999400 # ring" << ringSize << "\n";
1472  } // end of writing out atom types
1473  // Atoms
1474  outputFile << "\nAtoms\n\n";
1475  // -------
1476  // Write out the atom coordinates
1477  // Loop through atoms
1478  for (int i = 0; i < yCloud->pts.size(); i++) {
1479  iatom =
1480  yCloud->pts[i].atomID; // The actual ID can be different from the index
1481  // Write out coordinates
1482  // atomID molecule-tag atom-type q x y z
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";
1486 
1487  } // end of loop through all atoms in pointCloud
1488 
1489  // Print the bonds now!
1490  outputFile << "\nBonds\n\n";
1491  // Loop through all bonds
1492  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1493  //
1494  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
1495  << bonds[ibond][1] << "\n";
1496  } // end of for loop for bonds
1497 
1498  outputFile.close();
1499  // Once the datafile has been printed, exit
1500  return 0;
1501 }
1502 
1508  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1509  std::vector<std::vector<int>> rings, bool useBondFile, std::string bondFile,
1510  std::vector<int> listPrism, std::vector<std::vector<int>> nList,
1511  std::string filename) {
1512  std::ofstream outputFile;
1513  std::vector<int> atoms; // Holds all atom IDs to print
1514  int ringSize = rings[0].size(); // Ring size of each ring in rings
1515  int iatom; // Index, not atom ID
1516  bool padAtoms = false; // Add extra atoms if the atom IDs are skipped
1517  int prevAtomID = 0; // Check for previous atom ID
1518  int dummyAtoms = 0; // Number of dummy atoms to fill
1519  int dummyID;
1520  int jatom; // Array index is 1 less than the ID (index for dummy atom)
1521  int bondTypes = 1;
1522  // Bond stuff
1523  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1524  // containing the atom IDs of each bond
1525  bool atomOne, atomTwo; // If bond atoms are in the prism or not
1526  bool isPrismBond;
1527 
1528  // ----------------
1529  // Return if there are no prisms
1530  if (listPrism.size() == 0) {
1531  return 1;
1532  }
1533 
1534  // ---------------
1535  // Get the bonds
1536  if (useBondFile) {
1537  // Bonds from file
1538  bonds = sinp::readBonds(bondFile);
1539  } // get bonds from file
1540  else {
1541  bonds = bond::populateBonds(nList, yCloud);
1542  } // Bonds from rings
1543  //
1544  // ----------------
1545  // Otherwise create file
1546  // Create output dir if it doesn't exist already
1547  const char *path = "../output"; // relative to the build directory
1548  fs::path dir(path);
1549  // if (fs::create_directory(dir)) {
1550  // std::cerr << "Output directory created\n";
1551  // }
1552  // ----------------
1553  // Get all the unique atomIDs of the atoms in the rings of this type
1554  // Put all atom IDs into one 1-D vector
1555  size_t total_size{0};
1556  // Get the total number of atoms (repeated)
1557  total_size = listPrism.size() * ringSize;
1558  // Reserve this size inside atoms
1559  atoms.reserve(total_size);
1560  // Fill up all these atom IDs
1561  for (int iring = 0; iring < listPrism.size(); iring++) {
1562  std::move(rings[listPrism[iring]].begin(), rings[listPrism[iring]].end(),
1563  std::back_inserter(atoms));
1564  } // end of loop through all rings in the list
1565 
1566  // Sort the array according to atom ID, which will be needed to get the
1567  // unique IDs and to remove duplicates
1568  sort(atoms.begin(), atoms.end());
1569  // Get the unique atom IDs
1570  auto ip = std::unique(atoms.begin(), atoms.end());
1571  // Resize the vector to remove undefined terms
1572  atoms.resize(std::distance(atoms.begin(), ip));
1573  // If the number of atoms is less than the total nop, add dummy atoms
1574  if (atoms.size() != yCloud->nop) {
1575  padAtoms = true;
1576  bondTypes = 2;
1577  }
1578  // ----------------
1579  // Write output to file inside the output directory
1580  outputFile.open("../output/" + filename);
1581  // FORMAT:
1582  // Comment Line
1583  // 4 atoms
1584  // 4 bonds
1585  // 0 angles
1586  // 0 dihedrals
1587  // 0 impropers
1588  // 1 atom types
1589  // 1 bond types
1590  // 0 angle types
1591  // 0 dihedral types
1592  // 0 improper types
1593  // -1.124000 52.845002 xlo xhi
1594  // 0.000000 54.528999 ylo yhi
1595  // 1.830501 53.087501 zlo zhi
1596 
1597  // Masses
1598 
1599  // 1 15.999400 # O
1600 
1601  // Atoms
1602 
1603  // 1 1 1 0 20.239 1.298 6.873 # O
1604  // 2 1 1 0 0 5.193 6.873 # O
1605  // 3 1 1 0 2.249 1.298 6.873 # O
1606 
1607  // -------
1608  // Write the header
1609  // Write comment line
1610  outputFile << "Written out by D-SEAMS\n";
1611  // Write out the number of atoms
1612  outputFile << yCloud->pts.size() << " "
1613  << "atoms"
1614  << "\n";
1615  // Number of bonds
1616  outputFile << bonds.size() << " bonds"
1617  << "\n";
1618  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1619  // If padded atoms are required, two atom types will be required
1620  if (padAtoms) {
1621  outputFile << "2 atom types\n";
1622  } else {
1623  outputFile << "1 atom types\n";
1624  } // end of atom types
1625  outputFile
1626  << bondTypes
1627  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1628  // Box lengths
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";
1632  // Masses
1633  outputFile << "\nMasses\n\n";
1634  outputFile << "1 15.999400 # O\n";
1635  if (padAtoms) {
1636  outputFile << "2 1.0 # dummy\n";
1637  }
1638  // Atoms
1639  outputFile << "\nAtoms\n\n";
1640  // -------
1641  // Write out the atom coordinates
1642  // Loop through atoms
1643  for (int i = 0; i < atoms.size(); i++) {
1644  iatom = atoms[i] - 1; // The actual index is one less than the ID
1645  // -----------
1646  // Pad out
1647  // Fill in dummy atoms if some have been skipped
1648  if (atoms[i] != prevAtomID + 1) {
1649  dummyAtoms = atoms[i] - prevAtomID - 1;
1650  dummyID = prevAtomID;
1651  // Loop to write out dummy atoms
1652  for (int j = 0; j < dummyAtoms; j++) {
1653  dummyID++;
1654  jatom = dummyID - 1;
1655  // 1 molecule-tag atom-type q x y z
1656  outputFile << dummyID << " " << yCloud->pts[jatom].molID << " 2 0 "
1657  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
1658  << yCloud->pts[jatom].z << "\n";
1659  } // end of dummy atom write-out
1660  } // end of check for dummy atom printing
1661  // -----------
1662  // Write out coordinates
1663  // 1 molecule-tag atom-type q x y z
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";
1667  // update the previous atom ID
1668  prevAtomID = atoms[i];
1669  } // end of loop through all atoms in atomID
1670 
1671  // Fill in the rest of the dummy atoms
1672  if (atoms[atoms.size() - 1] != yCloud->nop) {
1673  //
1674  for (int id = atoms[atoms.size() - 1] + 1; id <= yCloud->nop; id++) {
1675  jatom = id - 1;
1676  outputFile << id << " " << yCloud->pts[jatom].molID << " 2 0 "
1677  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
1678  << yCloud->pts[jatom].z << "\n";
1679  } // end of printing out dummy atoms
1680  }
1681 
1682  // Print the bonds now!
1683  outputFile << "\nBonds\n\n";
1684  // Loop through all bonds
1685  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1686  // Init
1687  isPrismBond = false;
1688  atomOne = false;
1689  atomTwo = false;
1690  // --------
1691  // Check if the bond is in the prism or not
1692  auto it = std::find(atoms.begin() + 1, atoms.end(), bonds[ibond][0]);
1693  if (it != atoms.end()) {
1694  atomOne = true;
1695  } else if (bonds[ibond][0] == atoms[0]) {
1696  atomOne = true;
1697  } else if (bonds[ibond][0] == atoms[atoms.size() - 1]) {
1698  atomOne = true;
1699  }
1700 
1701  auto it1 = std::find(atoms.begin() + 1, atoms.end(), bonds[ibond][1]);
1702  if (it1 != atoms.end()) {
1703  atomTwo = true;
1704  } else if (bonds[ibond][1] == atoms[0]) {
1705  atomTwo = true;
1706  } else if (bonds[ibond][1] == atoms[atoms.size() - 1]) {
1707  atomTwo = true;
1708  }
1709 
1710  if (atomOne == false || atomTwo == false) {
1711  isPrismBond = false;
1712  } else {
1713  isPrismBond = true;
1714  }
1715  // --------
1716  if (isPrismBond) {
1717  // is inside the prism (type 1)
1718  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
1719  << bonds[ibond][1] << "\n";
1720  } else {
1721  // not inside the prism (type 2)
1722  outputFile << ibond + 1 << " 2 " << bonds[ibond][0] << " "
1723  << bonds[ibond][1] << "\n";
1724  }
1725 
1726  } // end of for loop for bonds
1727 
1728  outputFile.close();
1729  // Once the datafile has been printed, exit
1730  return 0;
1731 }
1732 
1738  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
1740  cage::cageType type, int numCages, std::string filename) {
1741  std::ofstream outputFile;
1742  std::vector<int> atoms; // Holds all atom IDs to print
1743  int ringSize = rings[0].size(); // Ring size of each ring in rings
1744  int iatom; // Index, not atom ID
1745  bool padAtoms = false; // Add extra atoms if the atom IDs are skipped
1746  int prevAtomID = 0; // Check for previous atom ID
1747  int dummyAtoms = 0; // Number of dummy atoms to fill
1748  int dummyID;
1749  int iring; // Current ring index (vector index)
1750  int jatom; // Array index is 1 less than the ID (index for dummy atom)
1751  int bondTypes = 1;
1752  std::string actualCageType; // The actual name of the cage types
1753  // Bond stuff
1754  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1755  // containing the atom IDs of each bond
1756  int nRings = 0; // Number of rings
1757 
1758  // ----------------
1759  // Return if there are no cages at all
1760  if ((*cageList).size() == 0) {
1761  return 1;
1762  }
1763 
1764  // Return if there are no cages of the required type
1765  if (numCages == 0) {
1766  return 1;
1767  }
1768  // ---------------
1769  // Get the bonds
1770  bonds = bond::createBondsFromCages(rings, cageList, type, &nRings);
1771  //
1772  // ----------------
1773  // Otherwise create file
1774  // Create output dir if it doesn't exist already
1775  const char *path = "../output"; // relative to the build directory
1776  fs::path dir(path);
1777  // if (fs::create_directory(dir)) {
1778  // std::cerr << "Output directory created\n";
1779  // }
1780  // ----------------
1781  // Get all the unique atomIDs of the atoms in the rings of this type
1782  // Put all atom IDs into one 1-D vector
1783  size_t total_size{0};
1784  // Get the total number of atoms (repeated)
1785  total_size = nRings * ringSize;
1786  // Reserve this size inside atoms
1787  atoms.reserve(total_size);
1788  // Fill up all these atom IDs
1789  // Loop through every cage in cageList
1790  for (int icage = 0; icage < (*cageList).size(); icage++) {
1791  // Skip if the cage is of a different type
1792  if ((*cageList)[icage].type != type) {
1793  continue;
1794  }
1795  // Loop through every ring inside Cage
1796  for (int k = 0; k < (*cageList)[icage].rings.size(); k++) {
1797  iring = (*cageList)[icage].rings[k]; // Current ring index
1798  std::move(rings[iring].begin(), rings[iring].end(),
1799  std::back_inserter(atoms));
1800  } // end of loop through every ring in icage
1801 
1802  } // end of loop through all cages in cageList
1803 
1804  // Sort the array according to atom ID, which will be needed to get the
1805  // unique IDs and to remove duplicates
1806  sort(atoms.begin(), atoms.end());
1807  // Get the unique atom IDs
1808  auto ip = std::unique(atoms.begin(), atoms.end());
1809  // Resize the vector to remove undefined terms
1810  atoms.resize(std::distance(atoms.begin(), ip));
1811  // If the number of atoms is less than the total nop, add dummy atoms
1812  if (atoms.size() != yCloud->nop) {
1813  padAtoms = true;
1814  bondTypes = 1;
1815  }
1816  // ----------------
1817  // Write output to file inside the output directory
1818  outputFile.open("../output/" + filename);
1819  // FORMAT:
1820  // Comment Line
1821  // 4 atoms
1822  // 4 bonds
1823  // 0 angles
1824  // 0 dihedrals
1825  // 0 impropers
1826  // 1 atom types
1827  // 1 bond types
1828  // 0 angle types
1829  // 0 dihedral types
1830  // 0 improper types
1831  // -1.124000 52.845002 xlo xhi
1832  // 0.000000 54.528999 ylo yhi
1833  // 1.830501 53.087501 zlo zhi
1834 
1835  // Masses
1836 
1837  // 1 15.999400 # O
1838 
1839  // Atoms
1840 
1841  // 1 1 1 0 20.239 1.298 6.873 # O
1842  // 2 1 1 0 0 5.193 6.873 # O
1843  // 3 1 1 0 2.249 1.298 6.873 # O
1844 
1845  // -------
1846  // Write the header
1847  // Write comment line
1848  outputFile << "Written out by D-SEAMS\n";
1849  // Write out the number of atoms
1850  outputFile << yCloud->pts.size() << " "
1851  << "atoms"
1852  << "\n";
1853  // Number of bonds
1854  outputFile << bonds.size() << " bonds"
1855  << "\n";
1856  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1857  // If padded atoms are required, two atom types will be required
1858  if (padAtoms) {
1859  outputFile << "2 atom types\n";
1860  } else {
1861  outputFile << "1 atom types\n";
1862  } // end of atom types
1863  outputFile
1864  << bondTypes
1865  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1866  // Box lengths
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";
1870  // Masses
1871  outputFile << "\nMasses\n\n";
1872  // For DDCs and HCs
1873  if (type == cage::HexC) {
1874  actualCageType = "HC";
1875  } else if (type == cage::DoubleDiaC) {
1876  actualCageType = "DDC";
1877  } else {
1878  actualCageType = "error";
1879  }
1880  //
1881  outputFile << "1 15.999400 # " << actualCageType << "\n";
1882  if (padAtoms) {
1883  outputFile << "2 1.0 # dummy\n";
1884  }
1885  // Atoms
1886  outputFile << "\nAtoms\n\n";
1887  // -------
1888  // Write out the atom coordinates
1889  // Loop through atoms
1890  for (int i = 0; i < atoms.size(); i++) {
1891  iatom = atoms[i] - 1; // The actual index is one less than the ID
1892  // -----------
1893  // Pad out
1894  // Fill in dummy atoms if some have been skipped
1895  if (atoms[i] != prevAtomID + 1) {
1896  dummyAtoms = atoms[i] - prevAtomID - 1;
1897  dummyID = prevAtomID;
1898  // Loop to write out dummy atoms
1899  for (int j = 0; j < dummyAtoms; j++) {
1900  dummyID++;
1901  jatom = dummyID - 1;
1902  // 1 molecule-tag atom-type q x y z
1903  outputFile << dummyID << " " << yCloud->pts[jatom].molID << " 2 0 "
1904  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
1905  << yCloud->pts[jatom].z << "\n";
1906  } // end of dummy atom write-out
1907  } // end of check for dummy atom printing
1908  // -----------
1909  // Write out coordinates
1910  // 1 molecule-tag atom-type q x y z
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";
1914  // update the previous atom ID
1915  prevAtomID = atoms[i];
1916  } // end of loop through all atoms in atomID
1917 
1918  // Fill in the rest of the dummy atoms
1919  if (atoms[atoms.size() - 1] != yCloud->nop) {
1920  //
1921  for (int id = atoms[atoms.size() - 1] + 1; id <= yCloud->nop; id++) {
1922  jatom = id - 1;
1923  outputFile << id << " " << yCloud->pts[jatom].molID << " 2 0 "
1924  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
1925  << yCloud->pts[jatom].z << "\n";
1926  } // end of printing out dummy atoms
1927  }
1928 
1929  // Print the bonds now!
1930  outputFile << "\nBonds\n\n";
1931  // Loop through all bonds
1932  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1933  // write out the bond
1934  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
1935  << bonds[ibond][1] << "\n";
1936 
1937  } // end of for loop for bonds
1938 
1939  outputFile.close();
1940  // Once the datafile has been printed, exit
1941  return 0;
1942 }
1943 
1944 // Legacy
1945 
1950  std::string path, std::string outFile) {
1951  std::ofstream outputFile;
1952  // ----------------
1953  // Make the output directory if it doesn't exist
1954  sout::makePath(path);
1955  // ----------------
1956  // Create a new file in the output directory
1957  outputFile.open(path + outFile, std::ios_base::app);
1958 
1959  // Append stuff
1960  // -----------------------
1961  // Header
1962 
1963  // The format of the LAMMPS trajectory file is:
1964  // ITEM: TIMESTEP
1965  // 0
1966  // ITEM: NUMBER OF ATOMS
1967  // 4096
1968  // ITEM: BOX BOUNDS pp pp pp
1969  // -7.9599900000000001e-01 5.0164000000000001e+01
1970  // -7.9599900000000001e-01 5.0164000000000001e+01
1971  // -7.9599900000000001e-01 5.0164000000000001e+01
1972  // ITEM: ATOMS id type x y z
1973  // 1 1 0 0 0 etc
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]; // print xlo xhi etc
1982  // print out the tilt factors too if it is a triclinic box
1983  if (yCloud->box.size() == 2 * yCloud->boxLow.size()) {
1984  outputFile << " "
1985  << yCloud->box[k + yCloud->boxLow.size()]; // this would be +2
1986  // for a 2D box
1987  }
1988  outputFile << "\n"; // print end of line
1989  } // end of printing box lengths
1990  outputFile << "ITEM: ATOMS id mol type x y z\n";
1991  // -----------------------
1992  // Atom lines
1993  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
1994  outputFile << yCloud->pts[iatom].atomID << " " << yCloud->pts[iatom].molID;
1995 
1996  // Cubic ice
1997  if (yCloud->pts[iatom].iceType == molSys::cubic) {
1998  outputFile << " Ic " << yCloud->pts[iatom].x << " "
1999  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2000  } // end of cubic ice
2001  // Hexagonal ice
2002  else if (yCloud->pts[iatom].iceType == molSys::hexagonal) {
2003  outputFile << " Ih " << yCloud->pts[iatom].x << " "
2004  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2005  } // end hexagonal ice
2006  // water
2007  else if (yCloud->pts[iatom].iceType == molSys::water) {
2008  outputFile << " wat " << yCloud->pts[iatom].x << " "
2009  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2010  } // end water
2011  // interfacial
2012  else if (yCloud->pts[iatom].iceType == molSys::interfacial) {
2013  outputFile << " intFc " << yCloud->pts[iatom].x << " "
2014  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2015  } // end interfacial
2016  // clathrate
2017  else if (yCloud->pts[iatom].iceType == molSys::clathrate) {
2018  outputFile << " clathrate " << yCloud->pts[iatom].x << " "
2019  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2020  } // end clathrate
2021  // interClathrate
2022  else if (yCloud->pts[iatom].iceType == molSys::interClathrate) {
2023  outputFile << " interClathrate " << yCloud->pts[iatom].x << " "
2024  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2025  } // end interClathrate
2026  // unclassified
2027  else if (yCloud->pts[iatom].iceType == molSys::unclassified) {
2028  outputFile << " unclassified " << yCloud->pts[iatom].x << " "
2029  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2030  } // end unclassified
2031  // reCubic
2032  else if (yCloud->pts[iatom].iceType == molSys::reCubic) {
2033  outputFile << " reIc " << yCloud->pts[iatom].x << " "
2034  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2035  } // end reCubic
2036  // reHexagonal
2037  else {
2038  outputFile << " reIh " << yCloud->pts[iatom].x << " "
2039  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2040  } // end reHex
2041 
2042  } // end of loop through all atoms
2043 
2044  // Close the file
2045  outputFile.close();
2046  return 0;
2047 }
2048 
2055  std::vector<double> avgQ6) {
2056  std::ofstream cijFile;
2057  std::ofstream q3File;
2058  std::ofstream q6File;
2059  // Create a new file in the output directory
2060  int nNumNeighbours;
2061  double avgQ3;
2062 
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);
2066 
2067  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
2068  if (yCloud->pts[iatom].type != 1) {
2069  continue;
2070  }
2071  // Check for slice later
2072  nNumNeighbours = nList[iatom].size() - 1;
2073  avgQ3 = 0.0;
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;
2077  } // Loop through neighbours
2078  avgQ3 /= nNumNeighbours;
2079  q3File << avgQ3 << "\n";
2080  q6File << avgQ6[iatom] << "\n";
2081  } // loop through all atoms
2082 
2083  // Close the file
2084  cijFile.close();
2085  q3File.close();
2086  q6File.close();
2087 
2088  return 0;
2089 }
2090 
2095  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
2096  std::string fileName, bool isSlice, int largestIceCluster) {
2097  std::ofstream clusterFile;
2098  // Create a new file in the output directory
2099  clusterFile.open(fileName, std::ofstream::out | std::ofstream::app);
2100  clusterFile << yCloud->currentFrame << " " << largestIceCluster << "\n";
2101  // Close the file
2102  clusterFile.close();
2103  return 0;
2104 }
2105 
2112  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
2114  std::string path, bool bondsBetweenDummy) {
2115  //
2116  std::ofstream outputFile;
2117  int iatom; // Index, not atom ID
2118  int currentAtomType; // Current atom type: a value from 1 to 4
2119  int numAtomTypes = 6; // DDC, HC, Mixed, dummy, mixed2 and pnc
2120  int bondTypes = 1;
2121  // Bond stuff
2122  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
2123  // containing the atom IDs of each bond
2124  std::string filename =
2125  "system-" + std::to_string(yCloud->currentFrame) + ".data";
2126 
2127  // ---------------
2128  // Get the bonds
2129  if (bondsBetweenDummy) {
2130  bonds = bond::populateBonds(nList, yCloud);
2131  } // create bonds between dummy atoms
2132  else {
2133  bonds = bond::populateBonds(nList, yCloud, atomTypes);
2134  } // only create bonds between non-dummy atoms
2135  //
2136  // ----------------
2137  // Make the output directory if it doesn't exist
2138  sout::makePath(path);
2139  std::string outputDirName = path + "bulkTopo";
2140  sout::makePath(outputDirName);
2141  outputDirName = path + "bulkTopo/dataFiles/";
2142  sout::makePath(outputDirName);
2143  // ----------------
2144  // Write output to file inside the output directory
2145  outputFile.open(path + "bulkTopo/dataFiles/" + filename);
2146  // FORMAT:
2147  // Comment Line
2148  // 4 atoms
2149  // 4 bonds
2150  // 0 angles
2151  // 0 dihedrals
2152  // 0 impropers
2153  // 1 atom types
2154  // 1 bond types
2155  // 0 angle types
2156  // 0 dihedral types
2157  // 0 improper types
2158  // -1.124000 52.845002 xlo xhi
2159  // 0.000000 54.528999 ylo yhi
2160  // 1.830501 53.087501 zlo zhi
2161 
2162  // Masses
2163 
2164  // 1 15.999400 # O
2165 
2166  // Atoms
2167 
2168  // 1 1 1 0 20.239 1.298 6.873 # O
2169  // 2 1 1 0 0 5.193 6.873 # O
2170  // 3 1 1 0 2.249 1.298 6.873 # O
2171 
2172  // -------
2173  // Write the header
2174  // Write comment line
2175  outputFile << "Written out by D-SEAMS\n";
2176  // Write out the number of atoms
2177  outputFile << yCloud->pts.size() << " "
2178  << "atoms"
2179  << "\n";
2180  // Number of bonds
2181  outputFile << bonds.size() << " bonds"
2182  << "\n";
2183  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
2184  // There are maxDepth-2 total types of prisms + dummy
2185  outputFile << numAtomTypes << " atom types\n";
2186  // Bond types
2187  outputFile
2188  << bondTypes
2189  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
2190  // Box lengths
2191  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
2192  << " xlo xhi\n";
2193  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
2194  << " ylo yhi\n";
2195  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
2196  << " zlo zhi\n";
2197  // Masses
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";
2205  // Atoms
2206  outputFile << "\nAtoms\n\n";
2207  // -------
2208  // Write out the atom coordinates
2209  // Loop through atoms
2210  for (int i = 0; i < yCloud->pts.size(); i++) {
2211  iatom =
2212  yCloud->pts[i].atomID; // The actual ID can be different from the index
2213  //
2214  // Get the atom type
2215  // hc atom type
2216  if (atomTypes[i] == cage::hc) {
2217  currentAtomType = 2;
2218  } // hc
2219  else if (atomTypes[i] == cage::ddc) {
2220  currentAtomType = 3;
2221  } // ddc
2222  // mixed
2223  else if (atomTypes[i] == cage::mixed) {
2224  currentAtomType = 4;
2225  } // mixed
2226  // pnc
2227  else if (atomTypes[i] == cage::pnc) {
2228  currentAtomType = 5;
2229  } // mixed
2230  // pnc and DDCs/HCs mixed
2231  else if (atomTypes[i] == cage::mixed2) {
2232  currentAtomType = 6;
2233  } // mixed
2234  // dummy
2235  else {
2236  currentAtomType = 1;
2237  } // dummy
2238  //
2239  // Write out coordinates
2240  // atomID molecule-tag atom-type q x y z
2241  outputFile << iatom << " " << yCloud->pts[i].molID << " " << currentAtomType
2242  << " 0 " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
2243  << yCloud->pts[i].z << "\n";
2244 
2245  } // end of loop through all atoms in pointCloud
2246 
2247  // Print the bonds now!
2248  outputFile << "\nBonds\n\n";
2249  // Loop through all bonds
2250  for (int ibond = 0; ibond < bonds.size(); ibond++) {
2251  //
2252  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
2253  << bonds[ibond][1] << "\n";
2254  } // end of for loop for bonds
2255 
2256  // Once the datafile has been printed, exit
2257  return 0;
2258 }
2259 
2265  std::string path, molSys::PointCloud<molSys::Point<double>, double> *yCloud,
2266  std::vector<int> atoms, int clusterID, cage::cageType type) {
2267  //
2268  std::ofstream outputFile;
2269  std::string filename = "cluster-" + std::to_string(clusterID) + ".xyz";
2270  int nAtoms = atoms.size(); // Number of atoms
2271  int iatom; // Current atom
2272  // ----------------
2273  // Make the output directory if it doesn't exist
2274  sout::makePath(path);
2275  std::string outputDirName = path + "bulkTopo";
2276  sout::makePath(outputDirName);
2277  outputDirName = path + "bulkTopo/clusterXYZ/";
2278  sout::makePath(outputDirName);
2279  outputDirName = path + "bulkTopo/clusterXYZ/frame-" +
2280  std::to_string(yCloud->currentFrame);
2281  sout::makePath(outputDirName);
2282  // ----------------
2283  // Write output to file inside the output directory
2284  outputFile.open(path + "bulkTopo/clusterXYZ/frame-" +
2285  std::to_string(yCloud->currentFrame) + "/" + filename);
2286 
2287  // Format of an XYZ file:
2288  // 1970
2289  // generated by VMD
2290  // O 43.603500 16.926201 15.215700
2291  // O 39.912601 14.775100 19.379200
2292  outputFile << nAtoms << "\n"; // Number of atoms
2293  outputFile
2294  << "Generated by d-SEAMS. 0 type=hc and 1 type =ddc\n"; // Comment line
2295  //
2296  // Write out all the atom coordinates
2297  for (int i = 0; i < nAtoms; i++) {
2298  iatom = atoms[i]; // Atom index to be printed out
2299  //
2300  outputFile << type << " " << yCloud->pts[iatom].x << " "
2301  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2302  } // end of loop through all atoms
2303 
2304  outputFile.close();
2305 
2306  return 0;
2307 }
T back_inserter(T... args)
T begin(T... args)
T close(T... args)
T distance(T... args)
T end(T... args)
T find(T... args)
std::vector< std::vector< int > > populateBonds(std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: bond.cpp:22
std::vector< std::vector< int > > createBondsFromCages(std::vector< std::vector< int >> rings, std::vector< cage::Cage > *cageList, cage::cageType type, int *nRings)
Definition: bond.cpp:365
cageType
Definition: cage.hpp:50
@ HexC
Definition: cage.hpp:50
@ DoubleDiaC
Definition: cage.hpp:50
@ pnc
Definition: cage.hpp:71
@ ddc
Definition: cage.hpp:71
@ mixed
Definition: cage.hpp:71
@ hc
Definition: cage.hpp:71
@ mixed2
Definition: cage.hpp:71
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.
Definition: cluster.cpp:36
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,...
Definition: generic.hpp:104
@ 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
@ unclassified
Not classified into any other category.
Definition: mol_sys.hpp:116
std::vector< std::vector< int > > readBonds(std::string filename)
Reads bonds into a vector of vectors from a file with a specific format.
Definition: seams_input.cpp:20
T move(T... args)
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)
T open(T... args)
T push_back(T... args)
T reserve(T... args)
T resize(T... args)
File for functions that read in files).
T size(T... args)
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
T to_string(T... args)
T unique(T... args)