sout Namespace Reference

Functions

bool isDirExist (const std::string &path)
 
int makePath (const std::string &path)
 
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. More...
 
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 printRDF (std::string fileName, std::vector< double > *rdfValues, double binwidth, int nbin)
 Function for printing out the RDF, given the filename. More...
 
int writeTopoBulkData (std::string path, int currentFrame, int numHC, int numDDC, int mixedRings, int basalRings, int prismaticRings, int firstFrame)
 
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. More...
 
int writeClusterStats (std::string path, int currentFrame, int largestCluster, int numOfClusters, int smallestCluster, double avgClusterSize, int firstFrame)
 Function for writing out cluster statistics. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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 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. More...
 
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 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 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. More...
 
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. More...
 
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. More...
 
int writeDump (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::string outFile)
 Generic function for writing out to a dump file. More...
 
int writeHisto (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::vector< double > avgQ6)
 
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. More...
 
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. More...
 

Function Documentation

◆ isDirExist()

bool sout::isDirExist ( const std::string path)
inline

Inline function for checking if the directory exists or not

Parameters
[in]pathThe path of the directory
Returns
True or false

Definition at line 40 of file seams_output.hpp.

40  {
41 #if defined(_WIN32)
42  struct _stat info;
43  if (_stat(path.c_str(), &info) != 0) {
44  return false;
45  }
46  return (info.st_mode & _S_IFDIR) != 0;
47 #else
48  struct stat info;
49  if (stat(path.c_str(), &info) != 0) {
50  return false;
51  }
52  return (info.st_mode & S_IFDIR) != 0;
53 #endif
54 }
T c_str(T... args)

◆ makePath()

int sout::makePath ( const std::string path)
inline

Inline function for creating the desried directory.

Parameters
[in]pathThe path of the directory

Definition at line 60 of file seams_output.hpp.

60  {
61 #if defined(_WIN32)
62  int ret = _mkdir(path.c_str());
63 #else
64  mode_t mode = 0755;
65  int ret = mkdir(path.c_str(), mode);
66 #endif
67  if (ret == 0)
68  return 0;
69 
70  switch (errno) {
71  case ENOENT:
72  // parent didn't exist, try to create it
73  {
74  int pos = path.find_last_of('/');
75  if (pos == std::string::npos)
76 #if defined(_WIN32)
77  pos = path.find_last_of('\\');
78  if (pos == std::string::npos)
79 #endif
80  return 1;
81  if (!makePath(path.substr(0, pos)))
82  return 1;
83  }
84 // now, try to create again
85 #if defined(_WIN32)
86  return 0 == _mkdir(path.c_str());
87 #else
88  return 0 == mkdir(path.c_str(), mode);
89 #endif
90 
91  case EEXIST:
92  // done!
93  if (isDirExist(path)) {
94  return 0;
95  } else {
96  return 1;
97  }
98 
99  default:
100  return 1;
101  }
102 }
T find_last_of(T... args)
bool isDirExist(const std::string &path)
int makePath(const std::string &path)
T substr(T... args)

◆ printRDF()

int sout::printRDF ( std::string  fileName,
std::vector< double > *  rdfValues,
double  binwidth,
int  nbin 
)

Function for printing out the RDF, given the filename.

Function for printing out the RDF to a file, given that the file already exists and given the filename.

Definition at line 1024 of file seams_output.cpp.

1025  {
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 }
T close(T... args)
T open(T... args)

◆ writeAllCages()

int sout::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 
)

Write out all cages of all types into a folder called cages inside the output directory

Function for writing out all the cages of all types to a folder in the output

Definition at line 283 of file seams_output.cpp.

287  {
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 }
cageType
Definition: cage.hpp:50
@ HexC
Definition: cage.hpp:50
@ DoubleDiaC
Definition: cage.hpp:50
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 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.

◆ writeBasalRingsHex()

int sout::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.

Function for printing out the basal rings only of the hexagonal cage described by the number cageNum Uses Boost!

Definition at line 437 of file seams_output.cpp.

439  {
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 }
T push_back(T... args)
T size(T... args)
T to_string(T... args)

◆ writeBasalRingsPrism()

int sout::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.

Function for printing out the basal rings of the prism blocks described by the number prismNum. Uses Boost!

Definition at line 646 of file seams_output.cpp.

650  {
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 }
T begin(T... args)
T end(T... args)
T find(T... args)
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

◆ writeCluster()

int sout::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.

Function to print out the largest ice cluster

Definition at line 2094 of file seams_output.cpp.

2096  {
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 }
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
int currentFrame
Collection of points.
Definition: mol_sys.hpp:168

◆ writeClusterStats()

int sout::writeClusterStats ( std::string  path,
int  currentFrame,
int  largestCluster,
int  numOfClusters,
int  smallestCluster,
double  avgClusterSize,
int  firstFrame 
)

Function for writing out cluster statistics.

Function for printing out cluster statistics

Definition at line 836 of file seams_output.cpp.

839  {
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 }

◆ writeDump()

int sout::writeDump ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::string  path,
std::string  outFile 
)

Generic function for writing out to a dump file.

Function for printing out info in PairCorrel struct

Definition at line 1949 of file seams_output.cpp.

1950  {
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 }
int nop
Current frame number.
Definition: mol_sys.hpp:169
std::vector< S > pts
Definition: mol_sys.hpp:167
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:170
std::vector< T > boxLow
Periodic box lengths.
Definition: mol_sys.hpp:171
@ 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

◆ writeEachCage()

int sout::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.

Function for printing out each cage's coordinates, when there is no volume slice Uses Boost!

Definition at line 344 of file seams_output.cpp.

347  {
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 }
T strcpy(T... args)

◆ writeHisto()

int sout::writeHisto ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< std::vector< int >>  nList,
std::vector< double >  avgQ6 
)

Function for printing out Q6, Cij and averaged Q3 values as single columns to text files The file names are cij, q6, q3

Function for printing out values of averaged Q6, averaged Q3 and Cij values

Definition at line 2053 of file seams_output.cpp.

2055  {
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 }

◆ writeLAMMPSdata()

int sout::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.

Prints out a LAMMPS data file, with some default options. Only Oxygen atoms are printed out

Definition at line 18 of file seams_output.cpp.

21  {
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 }
T back_inserter(T... args)
T distance(T... args)
T move(T... args)
T reserve(T... args)
T resize(T... args)
T sort(T... args)
T unique(T... args)

◆ writeLAMMPSdataAllPrisms()

int sout::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.

Prints out a LAMMPS data file for all the prisms for a single frame, with some default options. Only Oxygen atoms are printed out. Bonds are inferred from the rings vector of vectors

Definition at line 1251 of file seams_output.cpp.

1254  {
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 }
std::vector< std::vector< int > > populateBonds(std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: bond.cpp:22

◆ writeLAMMPSdataAllRings()

int sout::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.

Prints out a LAMMPS data file for all the rings for a single frame, with some default options. Only Oxygen atoms are printed out. Bonds are inferred from the rings vector of vectors

Definition at line 1385 of file seams_output.cpp.

1388  {
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 }

◆ writeLAMMPSdataCages()

int sout::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" 
)

Write out a lammps data file for DDCs or HCs, assuming that there is no slice

Prints out a LAMMPS data file for the either the DDCs or HCs, with some default options. Only Oxygen atoms are printed out

Definition at line 1737 of file seams_output.cpp.

1740  {
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 }
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

◆ writeLAMMPSdataPrisms()

int sout::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.

Prints out a LAMMPS data file for the prisms, with some default options. Only Oxygen atoms are printed out

Definition at line 1507 of file seams_output.cpp.

1511  {
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 }
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

◆ writeLAMMPSdataTopoBulk()

int sout::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 
)

Write a data file for a particular frame, writing out topological bulk ice structures (DDCs/HCs)

Prints out a LAMMPS data file for all the topological bulk ice for a single frame, with some default options. Only Oxygen atoms are printed out. Bonds are inferred from the neighbour list

Definition at line 2111 of file seams_output.cpp.

2114  {
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 }
@ 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

◆ writeLAMMPSdumpCages()

int sout::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.

Prints out a LAMMPS dump file for all the cages in bulk ice, for every frame, printing the RMSD per atom as well

Definition at line 1086 of file seams_output.cpp.

1089  {
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

◆ writeLAMMPSdumpINT()

int sout::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.

Prints out a LAMMPS dump file for all the prisms, for every frame, printing the RMSD per atom as well

Definition at line 1176 of file seams_output.cpp.

1179  {
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

◆ writePrismNum()

int sout::writePrismNum ( std::string  path,
std::vector< int >  nPrisms,
std::vector< int >  nDefPrisms,
std::vector< double >  heightPercent,
int  maxDepth,
int  currentFrame,
int  firstFrame 
)

Function for printing out the number of prism blocks, with or without slices. Be careful when using slices!

Function for printing out prism info, when there is no volume slice

Definition at line 872 of file seams_output.cpp.

875  {
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 }

◆ writePrisms()

int sout::writePrisms ( std::vector< int > *  basal1,
std::vector< int > *  basal2,
int  prismNum,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud 
)

Function for writing out each prism.

Function for printing out each info, when there is no volume slice Uses Boost!

Definition at line 212 of file seams_output.cpp.

214  {
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 }

◆ writeRingNum()

int sout::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 
)

Function for printing out the coverage area and the number of rings of each type

Function for printing out ring info, for a monolayer

Definition at line 918 of file seams_output.cpp.

923  {
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 }

◆ writeRings()

int sout::writeRings ( std::vector< std::vector< int >>  rings,
std::string  filename = "rings.dat" 
)

Function for printing out ring info, when there is no volume slice.

Function for printing out ring info, when there is no volume slice Uses Boost!

Definition at line 178 of file seams_output.cpp.

179  {
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 }

◆ writeTopoBulkData()

int sout::writeTopoBulkData ( std::string  path,
int  currentFrame,
int  numHC,
int  numDDC,
int  mixedRings,
int  basalRings,
int  prismaticRings,
int  firstFrame 
)

Function for printing out the number of DDCs, HCs, mixed rings, basal and prismatic rings

Function for out the number of DDCs, HCs, mixed rings, basal and prismatic rings.

Definition at line 1049 of file seams_output.cpp.

1051  {
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

◆ writeXYZcluster()

int sout::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.

Function for writing out the XYZ file of a particular cluster. The vector atoms contains the atom indices of the atoms to be written out

Definition at line 2264 of file seams_output.cpp.

2266  {
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 }