Functions
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 writeRingNumBulk (std::string path, int currentFrame, std::vector< int > nRings, 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 writeMoleculeIDsInSlice (std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 
int writeMoleculeIDsExpressionSelectOVITO (std::string path, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 
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 writeLAMMPSdumpSlice (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path)
 
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, bool isMonolayer=true)
 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 44 of file seams_output.hpp.

44  {
45 #if defined(_WIN32)
46  struct _stat info;
47  if (_stat(path.c_str(), &info) != 0) {
48  return false;
49  }
50  return (info.st_mode & _S_IFDIR) != 0;
51 #else
52  struct stat info;
53  if (stat(path.c_str(), &info) != 0) {
54  return false;
55  }
56  return (info.st_mode & S_IFDIR) != 0;
57 #endif
58 }
std::string::c_str
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 64 of file seams_output.hpp.

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

1227  {
1228  //
1229  std::ofstream outputFile; // For the output file
1230  double r; // Distance for the current bin
1231 
1232  // Append to the file
1233  outputFile.open(fileName, std::ios_base::app | std::ios_base::out);
1234 
1235  // Loop through all the bins
1236  for (int ibin = 0; ibin < nbin; ibin++) {
1237  //
1238  r = binwidth * (ibin + 0.5); // Current distance for ibin
1239  outputFile << r << " " << (*rdfValues)[ibin] << "\n";
1240  } // end of loop through all bins
1241 
1242  outputFile.close();
1243 
1244  return 0;
1245 }
std::ofstream::close
T close(T... args)
std::ofstream::open
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 287 of file seams_output.cpp.

291  {
292  int numDDC; // Number of DDCs
293  int numHC; // Number of HCs
294  int numMC; // Number of MCs
295  int totalCages = (*cageList).size(); // Total number of cages
296  cage::cageType type; // Current cage type
297 
298  // ---------------------------------
299  // Error handling!
300  if (totalCages == 0) {
301  std::cerr << "There are no cages to print.\n";
302  return 1;
303  }
304  // ---------------------------------
305  // Init
306  numDDC = 0;
307  numHC = 0;
308  numMC = 0;
309 
310  // Loop through every cage
311  for (int icage = 0; icage < totalCages; icage++) {
312  type = (*cageList)[icage].type;
313  // ------
314  // Add to the cage type and write out to the appropriate folders
315  // Hexagonal Cages
316  if (type == cage::cageType::HexC) {
317  numHC++;
318  sout::writeEachCage((*cageList)[icage].rings, numHC, type, rings, yCloud);
319  sout::writeBasalRingsHex((*cageList)[icage].rings, numHC, nList, rings);
320  } // end of write out of HCs
321  // Double diamond Cages
322  else if (type == cage::cageType::DoubleDiaC) {
323  numDDC++;
324  sout::writeEachCage((*cageList)[icage].rings, numDDC, type, rings,
325  yCloud);
326  } // end of write out of DDCs
327  // // Mixed Cages
328  // else if (type == cage::Mixed) {
329  // numMC++;
330  // sout::writeEachCage((*cageList)[icage].rings, numMC, type, rings,
331  // yCloud);
332  // } // end of write out of MCs
333  // // Error
334  else {
335  std::cerr << "The cage is of the wrong type\n";
336  continue;
337  } // some error
338  // ------
339  } // end of loop through all cages
340 
341  return 0;
342 }
cage::cageType
cageType
Definition: cage.hpp:54
cage::cageType::HexC
@ HexC
cage::cageType::DoubleDiaC
@ DoubleDiaC
sout::writeBasalRingsHex
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.
Definition: seams_output.cpp:441
sout::writeEachCage
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.
Definition: seams_output.cpp:348

◆ 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 441 of file seams_output.cpp.

443  {
444  std::ofstream outputFile;
445  std::string number = std::to_string(cageNum);
446  std::string filename = "basalRings" + number + ".dat";
447  int ringSize =
448  rings[0].size(); // Size of the ring; each ring contains n elements
449  int iatomIndex; // index of atom coordinate being written out TODO get rid of
450  // this
451  int iring; // Ring index of the current ring
452  // Variables for the hydrogen-bonded 'second' basal ring
453  std::vector<int>
454  basal2; // Elements are bonded to each element of basal1 in order
455  std::vector<int> basal1; // 'First' basal ring
456  std::vector<int>
457  unOrderedBasal2; // Unordered basal2 ring, the ith element is not
458  // necessarily bonded to the ith element of basal1
459  int iatom, jatom; // Atom numbers (starting from 1), not C++ indices; saved
460  // inside rings
461  int natom; // Neighbour list ID for iatom
462  int findAtom; // Atom number of atomID to find in neighbour list
463  bool atomFound; // bool to check if an atomID has been found in the neighbour
464  // list or not
465  // Variables for ordering basal2
466  // After finding the nearest neighbours, elements which are not nearest
467  // neighbours are assigned a value of zero.
468  int needle; // First non-zero element of basal2 after getting the nearest
469  // neighbours
470  int startNeedle; // Index of basal2; the first non-zero element of basal2
471  int startHayStack; // Index of unOrderedBasal2, corresponding to the first
472  // non-zero element of basal2
473  bool isClock; // Original order of unOrderedBasal2
474  int nextElement; // Next element in unOrderedBasal2 after startHayStack
475 
476  // ----------------
477  // Create output dir if it doesn't exist already
478  const char *path = "../output/cages"; // relative to the build directory
479  fs::path dir(path);
480  // if (fs::create_directory(dir)) {
481  // std::cerr << "Output Cage directory created\n";
482  // }
483  // ----------------
484  // Subdirectory
485 
486  const fs::path path1("../output/cages/hexBasalRings");
487  // fs::create_directories(path1);
488 
489  // Write output to file inside the output directory
490  std::string fileOutNameFull = "../output/cages/hexBasalRings/" + filename;
491  outputFile.open(fileOutNameFull);
492 
493  // Format:
494  // Coordinate IDs (starting from 1), ordered according to the input XYZ file
495  // The first line is a comment line
496  outputFile << "# Particle IDs in the two basal rings\n";
497 
498  // ---------------
499  // Find the nearest neighbours of basal1 elements in basal2
500  basal1 = rings[currentCage[0]]; // First basal ring
501  unOrderedBasal2 = rings[currentCage[1]]; // Unordered second basal ring
502 
503  for (int i = 0; i < basal1.size(); i++) {
504  iatom = basal1[i]; // This is the atom particle ID, not the C++ index
505 
506  // Search through unOrderedBasal2 for an element in the neighbourlist of
507  // iatom
508  for (int k = 0; k < unOrderedBasal2.size(); k++) {
509  findAtom =
510  unOrderedBasal2[k]; // Atom ID to find in the neighbour list of iatom
511 
512  atomFound = false; // init
513  jatom = 0;
514 
515  // Search through the neighbour list for findAtom
516  for (int n = 1; n < nList[iatom - 1].size(); n++) {
517  natom = nList[iatom - 1][n]; // Atom ID
518 
519  if (findAtom == natom) {
520  atomFound = true;
521  break;
522  } // Check
523  } // Loop through nList
524 
525  if (atomFound) {
526  jatom = natom;
527  break;
528  } // atom has been found
529  } // end of loop through all atomIDs in unOrderedBasal2
530  basal2.push_back(jatom);
531  } // end of loop through all the atomIDs in basal1
532 
533  // ---------------------------------------------------
534  // Get particles which are not nearest neighbours
535  // 'Alternately' ordered particles
536 
537  // ---------------
538  // Init
539  isClock = false;
540 
541  // ---------------
542  // Get the first non-zero index {needle}
543  for (int i = 0; i < 2; i++) {
544  if (basal2[i] != 0) {
545  needle = basal2[i]; // Set the needle to the first non-zero element
546  startNeedle = i; // Index of basal2
547  break; // Break out of the loop
548  } // end of checking for non-zero index
549  } // end of getting the first non-zero index
550 
551  // Find the index matching needle in unOrderedBasal2
552  for (int i = 0; i < unOrderedBasal2.size(); i++) {
553  if (unOrderedBasal2[i] == needle) {
554  startHayStack = i; // Index at which needle has been found
555  } // end of check for needle
556  } // end of search for needle in unOrderedBasal2
557 
558  // ---------------
559  // Check for 'clockwise' order
560  // Check 'next' element
561  nextElement = startHayStack + 2;
562  if (nextElement >= ringSize) {
563  nextElement -= ringSize;
564  }
565 
566  // Init (indices of basal2 and unOrderedBasal2 respectively)
567  iatom = 0;
568  jatom = 0;
569 
570  if (basal2[startNeedle + 2] == unOrderedBasal2[nextElement]) {
571  isClock = true;
572  // Fill the three elements in increasing order
573  for (int i = 1; i < ringSize; i += 2) {
574  iatom = startNeedle + i;
575  jatom = startHayStack + i;
576  // Make sure the indices are not larger than 6
577  if (iatom >= ringSize) {
578  iatom -= ringSize;
579  }
580  if (jatom >= ringSize) {
581  jatom -= ringSize;
582  }
583 
584  basal2[iatom] = unOrderedBasal2[jatom];
585  } // end of filling next two alternate elements
586  } // check to see if clockwise order is correct
587 
588  // ---------------
589  // Check for 'anticlockwise' order
590  // Check 'next' element
591  if (!isClock) {
592  iatom = 0;
593  jatom = 0; // init
594 
595  // First element
596  iatom = startNeedle + 2;
597  jatom = startHayStack - 2;
598  if (jatom < 0) {
599  jatom += ringSize;
600  }
601 
602  if (basal2[iatom] == unOrderedBasal2[jatom]) {
603  // Fill the three elements in increasing order
604  for (int i = 1; i < ringSize; i += 2) {
605  iatom = startNeedle + i;
606  jatom = startHayStack - i;
607  // Make sure the indices are not larger than 6
608  if (iatom > ringSize) {
609  iatom -= ringSize;
610  }
611  if (jatom < 0) {
612  jatom += ringSize;
613  }
614 
615  basal2[iatom] = unOrderedBasal2[jatom];
616  } // end of filling next two alternate elements
617  } // check to see if anticlockwise order is correct
618  else {
619  std::cerr << "Something is wrong with your HCs.\n";
620  return 1;
621  } // exit with error
622  } // End of check for anticlockwise stuff
623 
624  // ---------------------------------------------------
625  // Print out the ordered rings
626  // For hexagonal cages:
627  // Only print out basal1 and basal2
628 
629  // BASAL1
630  for (int i = 0; i < basal1.size(); i++) {
631  outputFile << basal1[i] << " ";
632  } // end of loop through basal1
633  outputFile << "\n";
634 
635  // BASAL2
636  for (int i = 0; i < basal2.size(); i++) {
637  outputFile << basal2[i] << " ";
638  } // end of loop through basal2
639 
640  // Close the output file
641  outputFile.close();
642 
643  return 0;
644 }
std::vector::push_back
T push_back(T... args)
std::vector::size
T size(T... args)
std::to_string
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 650 of file seams_output.cpp.

654  {
655  std::ofstream outputFile;
656  std::string number = std::to_string(prismNum);
657  std::string filename = "basalRings" + number + ".dat";
658  int ringSize =
659  (*basal1).size(); // Size of the ring; each ring contains n elements
660  int nBonds; // Number of bonds between the two deformed prisms
661  int l_k, m_k; // Atom ID in basal1 and basal2 respectively
662  int iatom,
663  jatom; // Index not ID of basal1 and basal2 respectively which match
664  int currentIatom, currentJatom; // Index not ID for matchedBasal1 and
665  // matchedBasal2 respectively
666  bool isNeighbour; // Basal rings should have at least one nearest neighbour
667  bool isClock; // The order is in the original 'direction' of basal2
668  bool isAntiClock; // The order must be reversed
669  std::vector<int> matchedBasal1, matchedBasal2; // Vectors with revised order
670 
671  // ----------------
672  // Path for deformed prisms
673  if (isDeformed) {
674  // Create output dir if it doesn't exist already
675  const char *path = "../output/deformed"; // relative to the build directory
676  fs::path dir(path);
677  // if (fs::create_directory(dir)) {
678  // std::cerr << "Output Cage directory created\n";
679  // }
680  // ----------------
681  // Subdirectory
682 
683  const fs::path path1("../output/deformed/basalRings");
684  // fs::create_directories(path1);
685 
686  // Write output to file inside the output directory
687  std::string fileOutNameFull = "../output/deformed/basalRings/" + filename;
688  outputFile.open(fileOutNameFull);
689  } else {
690  // Create output dir if it doesn't exist already
691  const char *path = "../output/perfect"; // relative to the build directory
692  fs::path dir(path);
693  // if (fs::create_directory(dir)) {
694  // std::cerr << "Output Cage directory created\n";
695  // }
696  // ----------------
697  // Subdirectory
698 
699  const fs::path path1("../output/perfect/basalRings");
700  // fs::create_directories(path1);
701 
702  // Write output to file inside the output directory
703  std::string fileOutNameFull = "../output/perfect/basalRings/" + filename;
704  outputFile.open(fileOutNameFull);
705  } // end of creating file paths
706 
707  // Format:
708  // Coordinate IDs (starting from 1), ordered according to the input XYZ file
709  // The first line is a comment line
710  outputFile << "# Particle IDs in the two basal rings\n";
711 
712  // ---------------
713  // Find the nearest neighbours of basal1 elements in basal2
714 
715  nBonds = 0;
716  isNeighbour = false;
717  // Loop through every element of basal1
718  for (int l = 0; l < ringSize; l++) {
719  l_k = (*basal1)[l]; // This is the atom particle ID, not the C++ index
720 
721  // Search for the nearest neighbour of l_k in basal2
722  // Loop through basal2 elements
723  for (int m = 0; m < ringSize; m++) {
724  m_k = (*basal2)[m]; // Atom ID to find in the neighbour list of iatom
725 
726  // Find m_k inside l_k neighbour list
727  auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
728 
729  // If the element has been found, for l1
730  if (it != nList[l_k].end()) {
731  isNeighbour = true;
732  iatom = l; // index of basal1
733  jatom = m; // index of basal2
734  break;
735  } // found element
736 
737  } // end of loop through all atomIDs in basal2
738 
739  if (isNeighbour) {
740  break;
741  } // nearest neighbour found
742  } // end of loop through all the atomIDs in basal1
743 
744  if (!isNeighbour) {
745  std::cerr << "Something is wrong with your deformed prism.\n";
746  return 1;
747  }
748  // ---------------------------------------------------
749  // Find out if the order of basal2 is 'clockwise' or 'anticlockwise'
750  isClock = false; // init
751  isAntiClock = false;
752 
753  // atom index (not ID)
754  int tempJfor, tempJback;
755 
756  tempJfor = jatom + 1;
757  tempJback = jatom - 1;
758 
759  if (jatom == ringSize - 1) {
760  tempJfor = 0;
761  tempJback = ringSize - 2;
762  }
763  if (jatom == 0) {
764  tempJfor = 1;
765  tempJback = ringSize - 1;
766  }
767 
768  int forwardJ = (*basal2)[tempJfor];
769  int backwardJ = (*basal2)[tempJback];
770  int currentI = (*basal1)[iatom];
771 
772  // Check clockwise
773  double distClock = gen::periodicDist(yCloud, currentI, forwardJ);
774  double distAntiClock = gen::periodicDist(yCloud, currentI, backwardJ);
775 
776  // Clockwise
777  if (distClock < distAntiClock) {
778  isClock = true;
779  } // end of clockwise check
780  // Anti-clockwise
781  if (distAntiClock < distClock) {
782  isAntiClock = true;
783  } // end of anti-clockwise check
784  // Some error
785  if (isClock == false && isAntiClock == false) {
786  // std::cerr << "The points are equidistant.\n";
787  return 1;
788  } // end of error handling
789  // ---------------------------------------------------
790  // Get the order of basal1 and basal2
791  for (int i = 0; i < ringSize; i++) {
792  currentIatom = iatom + i;
793  if (currentIatom >= ringSize) {
794  currentIatom -= ringSize;
795  } // end of basal1 element wrap-around
796 
797  // In clockwise order
798  if (isClock) {
799  currentJatom = jatom + i;
800  if (currentJatom >= ringSize) {
801  currentJatom -= ringSize;
802  } // wrap around
803  } // end of clockwise update
804  else {
805  currentJatom = jatom - i;
806  if (currentJatom < 0) {
807  currentJatom += ringSize;
808  } // wrap around
809  } // end of anti-clockwise update
810 
811  // Add to matchedBasal1 and matchedBasal2 now
812  matchedBasal1.push_back((*basal1)[currentIatom]);
813  matchedBasal2.push_back((*basal2)[currentJatom]);
814  }
815  // ---------------------------------------------------
816  // Print out the ordered rings
817  // For hexagonal cages:
818  // Only print out basal1 and basal2
819 
820  // BASAL1
821  for (int i = 0; i < matchedBasal1.size(); i++) {
822  outputFile << matchedBasal1[i] << " ";
823  } // end of loop through basal1
824  outputFile << "\n";
825 
826  // BASAL2
827  for (int i = 0; i < matchedBasal2.size(); i++) {
828  outputFile << matchedBasal2[i] << " ";
829  } // end of loop through basal2
830 
831  // Close the output file
832  outputFile.close();
833 
834  return 0;
835 }
std::begin
T begin(T... args)
std::end
T end(T... args)
std::find
T find(T... args)
gen::periodicDist
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:108

◆ 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 2380 of file seams_output.cpp.

2382  {
2383  std::ofstream clusterFile;
2384  // Create a new file in the output directory
2385  clusterFile.open(fileName, std::ofstream::out | std::ofstream::app);
2386  clusterFile << yCloud->currentFrame << " " << largestIceCluster << "\n";
2387  // Close the file
2388  clusterFile.close();
2389  return 0;
2390 }
clump::largestIceCluster
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:40
molSys::PointCloud::currentFrame
int currentFrame
Collection of points.
Definition: mol_sys.hpp:172

◆ 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 840 of file seams_output.cpp.

843  {
844  std::ofstream outputFile;
845  // ----------------
846  // Make the output directory if it doesn't exist
847  sout::makePath(path);
848  // ----------------
849  // Write output to file inside the output directory
850  outputFile.open(path + "clusterStats.dat",
851  std::ios_base::app | std::ios_base::out);
852 
853  // Format:
854  // Comment line
855  // 1 3 0 0 4 35 40 ....
856 
857  // ----------------
858  // Comment line for the first frame
859  if (currentFrame == firstFrame) {
860  outputFile << "Frame largestCluster numOfClusters smallestCluster "
861  "avgClusterSize\n";
862  }
863  // ----------------
864 
865  outputFile << currentFrame << " " << largestCluster << " " << numOfClusters
866  << " " << smallestCluster << " " << avgClusterSize << "\n";
867 
868  outputFile.close();
869 
870  return 0;
871 }

◆ 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 2235 of file seams_output.cpp.

2236  {
2237  std::ofstream outputFile;
2238  // ----------------
2239  // Make the output directory if it doesn't exist
2240  sout::makePath(path);
2241  // ----------------
2242  // Create a new file in the output directory
2243  outputFile.open(path + outFile, std::ios_base::app);
2244 
2245  // Append stuff
2246  // -----------------------
2247  // Header
2248 
2249  // The format of the LAMMPS trajectory file is:
2250  // ITEM: TIMESTEP
2251  // 0
2252  // ITEM: NUMBER OF ATOMS
2253  // 4096
2254  // ITEM: BOX BOUNDS pp pp pp
2255  // -7.9599900000000001e-01 5.0164000000000001e+01
2256  // -7.9599900000000001e-01 5.0164000000000001e+01
2257  // -7.9599900000000001e-01 5.0164000000000001e+01
2258  // ITEM: ATOMS id type x y z
2259  // 1 1 0 0 0 etc
2260  outputFile << "ITEM: TIMESTEP\n";
2261  outputFile << yCloud->currentFrame << "\n";
2262  outputFile << "ITEM: NUMBER OF ATOMS\n";
2263  outputFile << yCloud->nop << "\n";
2264  outputFile << "ITEM: BOX BOUNDS pp pp pp\n";
2265  for (int k = 0; k < yCloud->boxLow.size(); k++) {
2266  outputFile << yCloud->boxLow[k] << " "
2267  << yCloud->boxLow[k] + yCloud->box[k]; // print xlo xhi etc
2268  // print out the tilt factors too if it is a triclinic box
2269  if (yCloud->box.size() == 2 * yCloud->boxLow.size()) {
2270  outputFile << " "
2271  << yCloud->box[k + yCloud->boxLow.size()]; // this would be +2
2272  // for a 2D box
2273  }
2274  outputFile << "\n"; // print end of line
2275  } // end of printing box lengths
2276  outputFile << "ITEM: ATOMS id mol type x y z\n";
2277  // -----------------------
2278  // Atom lines
2279  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
2280  outputFile << yCloud->pts[iatom].atomID << " " << yCloud->pts[iatom].molID;
2281 
2282  // Cubic ice
2283  if (yCloud->pts[iatom].iceType == molSys::atom_state_type::cubic) {
2284  outputFile << " Ic " << yCloud->pts[iatom].x << " "
2285  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2286  } // end of cubic ice
2287  // Hexagonal ice
2288  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::hexagonal) {
2289  outputFile << " Ih " << yCloud->pts[iatom].x << " "
2290  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2291  } // end hexagonal ice
2292  // water
2293  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::water) {
2294  outputFile << " wat " << yCloud->pts[iatom].x << " "
2295  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2296  } // end water
2297  // interfacial
2298  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::interfacial) {
2299  outputFile << " intFc " << yCloud->pts[iatom].x << " "
2300  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2301  } // end interfacial
2302  // clathrate
2303  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::clathrate) {
2304  outputFile << " clathrate " << yCloud->pts[iatom].x << " "
2305  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2306  } // end clathrate
2307  // interClathrate
2308  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::interClathrate) {
2309  outputFile << " interClathrate " << yCloud->pts[iatom].x << " "
2310  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2311  } // end interClathrate
2312  // unclassified
2313  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::unclassified) {
2314  outputFile << " unclassified " << yCloud->pts[iatom].x << " "
2315  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2316  } // end unclassified
2317  // reCubic
2318  else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::reCubic) {
2319  outputFile << " reIc " << yCloud->pts[iatom].x << " "
2320  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2321  } // end reCubic
2322  // reHexagonal
2323  else {
2324  outputFile << " reIh " << yCloud->pts[iatom].x << " "
2325  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2326  } // end reHex
2327 
2328  } // end of loop through all atoms
2329 
2330  // Close the file
2331  outputFile.close();
2332  return 0;
2333 }
molSys::PointCloud::nop
int nop
Current frame number.
Definition: mol_sys.hpp:173
molSys::PointCloud::pts
std::vector< S > pts
Definition: mol_sys.hpp:171
molSys::PointCloud::box
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:174
molSys::PointCloud::boxLow
std::vector< T > boxLow
Periodic box lengths.
Definition: mol_sys.hpp:175
molSys::atom_state_type::hexagonal
@ hexagonal
Ih, or particle type signifying Hexagonal Ice.
molSys::atom_state_type::reCubic
@ reCubic
Reclassified as cubic ice, according to the order parameter.
molSys::atom_state_type::interfacial
@ interfacial
Interfacial ice: ice-like molecules which do not fulfill the strict criteria of the Ic or Ih phases.
molSys::atom_state_type::cubic
@ cubic
Ic, or particle type signifying Cubic Ice.
molSys::atom_state_type::interClathrate
@ interClathrate
Interfacial clathrate ice phase.
molSys::atom_state_type::water
@ water
Liquid/amorphous phase.
molSys::atom_state_type::clathrate
@ clathrate
Clathrate ice phase.
molSys::atom_state_type::unclassified
@ unclassified
Not classified into any other category.

◆ 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 348 of file seams_output.cpp.

351  {
352  std::ofstream outputFile;
353  std::string number = std::to_string(cageNum);
354  std::string filename = "cage" + number + ".dat";
355  int ringSize =
356  rings[0].size(); // Size of the ring; each ring contains n elements
357  int iatomIndex; // index of atom coordinate being written out
358  std::string actualCageType; // is icage a DDC, HC or MC?
359  char cageChar[100]; // is icage a DDC, HC or MC?
360  int iring; // Ring index of the current ring
361 
362  if (type == cage::cageType::HexC) {
363  strcpy(cageChar, "../output/cages/hexCages");
364  actualCageType = "hexCages";
365  } else if (type == cage::cageType::DoubleDiaC) {
366  strcpy(cageChar, "../output/cages/doubleDiaCages");
367  actualCageType = "doubleDiaCages";
368  } else {
369  // throw error
370  std::cerr << "The cage is of the wrong type. Exit\n";
371  return 1;
372  } //
373 
374  // ----------------
375  // Create output dir if it doesn't exist already
376  const char *path = "../output/cages"; // relative to the build directory
377  fs::path dir(path);
378  // if (fs::create_directory(dir)) {
379  // std::cerr << "Output Cage directory created\n";
380  // }
381  // ----------------
382  // Subdirectory
383 
384  const fs::path path1(cageChar);
385  // fs::create_directories(path1);
386 
387  // Write output to file inside the output directory
388  std::string fileOutNameFull =
389  "../output/cages/" + actualCageType + "/" + filename;
390  outputFile.open(fileOutNameFull);
391 
392  // Format:
393  // x y z coordinates of each node
394 
395  // For hexagonal cages:
396  if (type == cage::cageType::HexC) {
397  // Print out only basal ring atoms, since they describe the outer structure
398  // The first two rings are basal rings
399  for (int i = 0; i < 2; i++) {
400  iring = currentCage[i]; // Current iring
401  // Get every node of iring
402  for (int j = 0; j < ringSize; j++) {
403  iatomIndex = rings[iring][j] - 1; // C++ indices are one less
404  // Write out the coordinates to the file
405  outputFile << yCloud->pts[iatomIndex].x << " ";
406  outputFile << yCloud->pts[iatomIndex].y << " ";
407  outputFile << yCloud->pts[iatomIndex].z << " ";
408 
409  outputFile << "\n";
410  } // end of loop through iring
411  } // end of getting every ring in the current cage
412  } // end of printing basal ring atoms for hexagonal cages
413  // For currentCage
414  else {
415  // Loop through all cages (could be duplicates) TODO: remove duplicates
416  for (int i = 0; i < currentCage.size(); i++) {
417  iring = currentCage[i]; // Current iring
418  // Get every node of iring
419  for (int j = 0; j < ringSize; j++) {
420  iatomIndex = rings[iring][j] - 1; // C++ indices are one less
421  // Write out the coordinates to the file
422  outputFile << yCloud->pts[iatomIndex].x << " ";
423  outputFile << yCloud->pts[iatomIndex].y << " ";
424  outputFile << yCloud->pts[iatomIndex].z << " ";
425 
426  outputFile << "\n";
427  } // end of loop through iring
428  } // end of getting every ring in the current cage
429  } // end of cage printing (has duplicates)
430 
431  // Close the output file
432  outputFile.close();
433 
434  return 0;
435 }
std::strcpy
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 2339 of file seams_output.cpp.

2341  {
2342  std::ofstream cijFile;
2343  std::ofstream q3File;
2344  std::ofstream q6File;
2345  // Create a new file in the output directory
2346  int nNumNeighbours;
2347  double avgQ3;
2348 
2349  cijFile.open("cij.txt", std::ofstream::out | std::ofstream::app);
2350  q3File.open("q3.txt", std::ofstream::out | std::ofstream::app);
2351  q6File.open("q6.txt", std::ofstream::out | std::ofstream::app);
2352 
2353  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
2354  if (yCloud->pts[iatom].type != 1) {
2355  continue;
2356  }
2357  // Check for slice later
2358  nNumNeighbours = nList[iatom].size() - 1;
2359  avgQ3 = 0.0;
2360  for (int j = 0; j < nNumNeighbours; j++) {
2361  cijFile << yCloud->pts[iatom].c_ij[j].c_value << "\n";
2362  avgQ3 += yCloud->pts[iatom].c_ij[j].c_value;
2363  } // Loop through neighbours
2364  avgQ3 /= nNumNeighbours;
2365  q3File << avgQ3 << "\n";
2366  q6File << avgQ6[iatom] << "\n";
2367  } // loop through all atoms
2368 
2369  // Close the file
2370  cijFile.close();
2371  q3File.close();
2372  q6File.close();
2373 
2374  return 0;
2375 }

◆ 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 22 of file seams_output.cpp.

25  {
26  std::ofstream outputFile;
27  std::vector<int> atoms; // Holds all atom IDs to print
28  int ringSize = rings[0].size(); // Ring size of each ring in rings
29  int iatom; // Index, not atom ID
30  bool padAtoms = false; // Add extra atoms if the atom IDs are skipped
31  int prevAtomID = 0; // Check for previous atom ID
32  int dummyAtoms = 0; // Number of dummy atoms to fill
33  int dummyID;
34  int jatom; // Array index is 1 less than the ID (index for dummy atom)
35  // ----------------
36  // Return if there are no rings
37  if (rings.size() == 0) {
38  return 1;
39  }
40  // ----------------
41  // Otherwise create file
42  // Create output dir if it doesn't exist already
43  const char *path = "../output"; // relative to the build directory
44  fs::path dir(path);
45  // if (fs::create_directory(dir)) {
46  // std::cerr << "Output directory created\n";
47  // }
48  // ----------------
49  // Get all the unique atomIDs of the atoms in the rings of this type
50  // Put all atom IDs into one 1-D vector
51  size_t total_size{0};
52  // Get the total number of atoms (repeated)
53  total_size = rings.size() * ringSize;
54  // Reserve this size inside atoms
55  atoms.reserve(total_size);
56  // Fill up all these atom IDs
57  for (int iring = 0; iring < rings.size(); iring++) {
58  std::move(rings[iring].begin(), rings[iring].end(),
59  std::back_inserter(atoms));
60  } // end of loop through all rings in the list
61 
62  // Sort the array according to atom ID, which will be needed to get the
63  // unique IDs and to remove duplicates
64  sort(atoms.begin(), atoms.end());
65  // Get the unique atom IDs
66  auto ip = std::unique(atoms.begin(), atoms.end());
67  // Resize the vector to remove undefined terms
68  atoms.resize(std::distance(atoms.begin(), ip));
69  // If the number of atoms is less than the total nop, add dummy atoms
70  if (atoms.size() != yCloud->nop) {
71  padAtoms = true;
72  }
73  // ----------------
74  // Write output to file inside the output directory
75  outputFile.open("../output/" + filename);
76  // FORMAT:
77  // Comment Line
78  // 4 atoms
79  // 4 bonds
80  // 0 angles
81  // 0 dihedrals
82  // 0 impropers
83  // 1 atom types
84  // 1 bond types
85  // 0 angle types
86  // 0 dihedral types
87  // 0 improper types
88  // -1.124000 52.845002 xlo xhi
89  // 0.000000 54.528999 ylo yhi
90  // 1.830501 53.087501 zlo zhi
91 
92  // Masses
93 
94  // 1 15.999400 # O
95 
96  // Atoms
97 
98  // 1 1 1 0 20.239 1.298 6.873 # O
99  // 2 1 1 0 0 5.193 6.873 # O
100  // 3 1 1 0 2.249 1.298 6.873 # O
101 
102  // -------
103  // Write the header
104  // Write comment line
105  outputFile << "Written out by D-SEAMS\n";
106  // Write out the number of atoms
107  outputFile << atoms[atoms.size() - 1] << " "
108  << "atoms"
109  << "\n";
110  // Number of bonds
111  outputFile << bonds.size() << " bonds"
112  << "\n";
113  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
114  // If padded atoms are required, two atom types will be required
115  if (padAtoms) {
116  outputFile << "2 atom types\n";
117  } else {
118  outputFile << "1 atom types\n";
119  } // end of atom types
120  outputFile
121  << "1 bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
122  // Box lengths
123  outputFile << "0 " << yCloud->box[0] << " xlo xhi\n";
124  outputFile << "0 " << yCloud->box[1] << " ylo yhi\n";
125  outputFile << "0 " << yCloud->box[2] << " zlo zhi\n";
126  // Masses
127  outputFile << "\nMasses\n\n";
128  outputFile << "1 15.999400 # O\n";
129  if (padAtoms) {
130  outputFile << "2 1.0 # dummy\n";
131  }
132  // Atoms
133  outputFile << "\nAtoms\n\n";
134  // -------
135  // Write out the atom coordinates
136  // Loop through atoms
137  for (int i = 0; i < atoms.size(); i++) {
138  iatom = atoms[i] - 1; // The actual index is one less than the ID
139  // -----------
140  // Pad out
141  // Fill in dummy atoms if some have been skipped
142  if (atoms[i] != prevAtomID + 1) {
143  dummyAtoms = atoms[i] - prevAtomID - 1;
144  dummyID = prevAtomID;
145  // Loop to write out dummy atoms
146  for (int j = 0; j < dummyAtoms; j++) {
147  dummyID++;
148  jatom = dummyID - 1;
149  // 1 molecule-tag atom-type q x y z
150  outputFile << dummyID << " " << yCloud->pts[jatom].molID << " 2 0 "
151  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
152  << yCloud->pts[jatom].z << "\n";
153  } // end of dummy atom write-out
154  } // end of check for dummy atom printing
155  // -----------
156  // Write out coordinates
157  // 1 molecule-tag atom-type q x y z
158  outputFile << atoms[i] << " " << yCloud->pts[iatom].molID << " 1 0 "
159  << yCloud->pts[iatom].x << " " << yCloud->pts[iatom].y << " "
160  << yCloud->pts[iatom].z << "\n";
161  // update the previous atom ID
162  prevAtomID = atoms[i];
163  } // end of loop through all atoms in atomID
164 
165  // Print the bonds now!
166  outputFile << "\nBonds\n\n";
167  // Loop through all bonds
168  for (int ibond = 0; ibond < bonds.size(); ibond++) {
169  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
170  << bonds[ibond][1] << "\n";
171  }
172 
173  // Once the datafile has been printed, exit
174  return 0;
175 }
std::back_inserter
T back_inserter(T... args)
std::distance
T distance(T... args)
std::move
T move(T... args)
std::vector::reserve
T reserve(T... args)
std::vector::resize
T resize(T... args)
std::sort
T sort(T... args)
std::unique
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 1527 of file seams_output.cpp.

1530  {
1531  //
1532  std::ofstream outputFile;
1533  int iatom; // Index, not atom ID
1534  int bondTypes = 1;
1535  // Bond stuff
1536  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1537  // containing the atom IDs of each bond
1538  std::string filename =
1539  "system-prisms-" + std::to_string(yCloud->currentFrame) + ".data";
1540 
1541  // ---------------
1542  // Get the bonds
1543  bonds = bond::populateBonds(nList, yCloud);
1544  //
1545  // ----------------
1546  // Make the output directory if it doesn't exist
1547  sout::makePath(path);
1548  std::string outputDirName = path + "topoINT";
1549  sout::makePath(outputDirName);
1550  outputDirName = path + "topoINT/dataFiles/";
1551  sout::makePath(outputDirName);
1552  // ----------------
1553  // Write output to file inside the output directory
1554  outputFile.open(path + "topoINT/dataFiles/" + filename);
1555  // FORMAT:
1556  // Comment Line
1557  // 4 atoms
1558  // 4 bonds
1559  // 0 angles
1560  // 0 dihedrals
1561  // 0 impropers
1562  // 1 atom types
1563  // 1 bond types
1564  // 0 angle types
1565  // 0 dihedral types
1566  // 0 improper types
1567  // -1.124000 52.845002 xlo xhi
1568  // 0.000000 54.528999 ylo yhi
1569  // 1.830501 53.087501 zlo zhi
1570 
1571  // Masses
1572 
1573  // 1 15.999400 # O
1574 
1575  // Atoms
1576 
1577  // 1 1 1 0 20.239 1.298 6.873 # O
1578  // 2 1 1 0 0 5.193 6.873 # O
1579  // 3 1 1 0 2.249 1.298 6.873 # O
1580 
1581  // -------
1582  // Write the header
1583  // Write comment line
1584  outputFile << "Written out by D-SEAMS\n";
1585  // Write out the number of atoms
1586  outputFile << yCloud->pts.size() << " "
1587  << "atoms"
1588  << "\n";
1589  // Number of bonds
1590  outputFile << bonds.size() << " bonds"
1591  << "\n";
1592  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1593  // There are maxDepth-2 total types of prisms + dummy
1594  if (doShapeMatching) {
1595  outputFile << 2 * maxDepth - 2 << " atom types\n";
1596  } else {
1597  outputFile << maxDepth << " atom types\n";
1598  }
1599  // Bond types
1600  outputFile
1601  << bondTypes
1602  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1603  // Box lengths
1604  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1605  << " xlo xhi\n";
1606  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1607  << " ylo yhi\n";
1608  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1609  << " zlo zhi\n";
1610  // Masses
1611  outputFile << "\nMasses\n\n";
1612  outputFile << "1 15.999400 # dummy\n";
1613  outputFile << "2 1.0 # mixedRings \n";
1614  // There are maxDepth-2 other prism types
1615  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1616  outputFile << ringSize << " 15.999400 # prism" << ringSize << "\n";
1617  } // end of writing out perfect atom types
1618  // Write out the types for the deformed prism blocks
1619  if (doShapeMatching) {
1620  for (int ringSize = maxDepth + 1; ringSize <= 2 * maxDepth - 2;
1621  ringSize++) {
1622  int p = ringSize - maxDepth + 2;
1623  outputFile << ringSize << " 15.999400 # deformPrism" << p << "\n";
1624  } // end of writing out perfect atom types
1625  } // Deformed prism types
1626  // Atoms
1627  outputFile << "\nAtoms\n\n";
1628  // -------
1629  // Write out the atom coordinates
1630  // Loop through atoms
1631  for (int i = 0; i < yCloud->pts.size(); i++) {
1632  iatom =
1633  yCloud->pts[i].atomID; // The actual ID can be different from the index
1634  // Write out coordinates
1635  // atomID molecule-tag atom-type q x y z
1636  outputFile << iatom << " " << yCloud->pts[i].molID << " " << atomTypes[i]
1637  << " 0 " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
1638  << yCloud->pts[i].z << "\n";
1639 
1640  } // end of loop through all atoms in pointCloud
1641 
1642  // Print the bonds now!
1643  outputFile << "\nBonds\n\n";
1644  // Loop through all bonds
1645  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1646  //
1647  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
1648  << bonds[ibond][1] << "\n";
1649  } // end of for loop for bonds
1650 
1651  outputFile.close();
1652  // Once the datafile has been printed, exit
1653  return 0;
1654 }
bond::populateBonds
std::vector< std::vector< int > > populateBonds(std::vector< std::vector< int >> nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Definition: bond.cpp:26

◆ 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,
bool  isMonolayer = true 
)

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 1661 of file seams_output.cpp.

1664  {
1665  //
1666  std::ofstream outputFile;
1667  int iatom; // Index, not atom ID
1668  int bondTypes = 1;
1669  // Bond stuff
1670  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1671  // containing the atom IDs of each bond
1672  std::string filename =
1673  "system-rings-" + std::to_string(yCloud->currentFrame) + ".data";
1674  std::string pathName, pathFolder;
1675 
1676  // ---------------
1677  // Get the bonds
1678  bonds = bond::populateBonds(nList, yCloud);
1679  //
1680  // ----------------
1681  // Make the output directory if it doesn't exist
1682  if (isMonolayer)
1683  {
1684  pathFolder = "topoMonolayer";
1685  pathName = "topoMonolayer/dataFiles/";
1686  } else{
1687  pathFolder = "bulkTopo";
1688  pathName = "bulkTopo/dataFiles/";
1689  }
1690 
1691  sout::makePath(path);
1692  std::string outputDirName = path + pathFolder;
1693  sout::makePath(outputDirName);
1694  outputDirName = path + pathName;
1695  sout::makePath(outputDirName);
1696 
1697  // Write output to file inside the output directory
1698  outputFile.open(path + pathName + filename);
1699 
1700  // FORMAT:
1701  // Comment Line
1702  // 4 atoms
1703  // 4 bonds
1704  // 0 angles
1705  // 0 dihedrals
1706  // 0 impropers
1707  // 1 atom types
1708  // 1 bond types
1709  // 0 angle types
1710  // 0 dihedral types
1711  // 0 improper types
1712  // -1.124000 52.845002 xlo xhi
1713  // 0.000000 54.528999 ylo yhi
1714  // 1.830501 53.087501 zlo zhi
1715 
1716  // Masses
1717 
1718  // 1 15.999400 # O
1719 
1720  // Atoms
1721 
1722  // 1 1 1 0 20.239 1.298 6.873 # O
1723  // 2 1 1 0 0 5.193 6.873 # O
1724  // 3 1 1 0 2.249 1.298 6.873 # O
1725 
1726  // -------
1727  // Write the header
1728  // Write comment line
1729  outputFile << "Written out by D-SEAMS\n";
1730  // Write out the number of atoms
1731  outputFile << yCloud->pts.size() << " "
1732  << "atoms"
1733  << "\n";
1734  // Number of bonds
1735  outputFile << bonds.size() << " bonds"
1736  << "\n";
1737  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1738  // There are maxDepth-2 total types of prisms + dummy
1739  outputFile << maxDepth << " atom types\n";
1740  // Bond types
1741  outputFile
1742  << bondTypes
1743  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1744  // Box lengths
1745  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1746  << " xlo xhi\n";
1747  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1748  << " ylo yhi\n";
1749  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1750  << " zlo zhi\n";
1751  // Masses
1752  outputFile << "\nMasses\n\n";
1753  outputFile << "1 15.999400 # dummy\n";
1754  outputFile << "2 1.0 # \n";
1755  // There are maxDepth-2 other prism types
1756  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1757  outputFile << ringSize << " 15.999400 # ring" << ringSize << "\n";
1758  } // end of writing out atom types
1759  // Atoms
1760  outputFile << "\nAtoms\n\n";
1761  // -------
1762  // Write out the atom coordinates
1763  // Loop through atoms
1764  for (int i = 0; i < yCloud->pts.size(); i++) {
1765  iatom =
1766  yCloud->pts[i].atomID; // The actual ID can be different from the index
1767  // Write out coordinates
1768  // atomID molecule-tag atom-type q x y z
1769  outputFile << iatom << " " << yCloud->pts[i].molID << " " << atomTypes[i]
1770  << " 0 " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
1771  << yCloud->pts[i].z << "\n";
1772 
1773  } // end of loop through all atoms in pointCloud
1774 
1775  // Print the bonds now!
1776  outputFile << "\nBonds\n\n";
1777  // Loop through all bonds
1778  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1779  //
1780  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
1781  << bonds[ibond][1] << "\n";
1782  } // end of for loop for bonds
1783 
1784  outputFile.close();
1785  // Once the datafile has been printed, exit
1786  return 0;
1787 }

◆ 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 2023 of file seams_output.cpp.

2026  {
2027  std::ofstream outputFile;
2028  std::vector<int> atoms; // Holds all atom IDs to print
2029  int ringSize = rings[0].size(); // Ring size of each ring in rings
2030  int iatom; // Index, not atom ID
2031  bool padAtoms = false; // Add extra atoms if the atom IDs are skipped
2032  int prevAtomID = 0; // Check for previous atom ID
2033  int dummyAtoms = 0; // Number of dummy atoms to fill
2034  int dummyID;
2035  int iring; // Current ring index (vector index)
2036  int jatom; // Array index is 1 less than the ID (index for dummy atom)
2037  int bondTypes = 1;
2038  std::string actualCageType; // The actual name of the cage types
2039  // Bond stuff
2040  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
2041  // containing the atom IDs of each bond
2042  int nRings = 0; // Number of rings
2043 
2044  // ----------------
2045  // Return if there are no cages at all
2046  if ((*cageList).size() == 0) {
2047  return 1;
2048  }
2049 
2050  // Return if there are no cages of the required type
2051  if (numCages == 0) {
2052  return 1;
2053  }
2054  // ---------------
2055  // Get the bonds
2056  bonds = bond::createBondsFromCages(rings, cageList, type, &nRings);
2057  //
2058  // ----------------
2059  // Otherwise create file
2060  // Create output dir if it doesn't exist already
2061  const char *path = "../output"; // relative to the build directory
2062  fs::path dir(path);
2063  // if (fs::create_directory(dir)) {
2064  // std::cerr << "Output directory created\n";
2065  // }
2066  // ----------------
2067  // Get all the unique atomIDs of the atoms in the rings of this type
2068  // Put all atom IDs into one 1-D vector
2069  size_t total_size{0};
2070  // Get the total number of atoms (repeated)
2071  total_size = nRings * ringSize;
2072  // Reserve this size inside atoms
2073  atoms.reserve(total_size);
2074  // Fill up all these atom IDs
2075  // Loop through every cage in cageList
2076  for (int icage = 0; icage < (*cageList).size(); icage++) {
2077  // Skip if the cage is of a different type
2078  if ((*cageList)[icage].type != type) {
2079  continue;
2080  }
2081  // Loop through every ring inside Cage
2082  for (int k = 0; k < (*cageList)[icage].rings.size(); k++) {
2083  iring = (*cageList)[icage].rings[k]; // Current ring index
2084  std::move(rings[iring].begin(), rings[iring].end(),
2085  std::back_inserter(atoms));
2086  } // end of loop through every ring in icage
2087 
2088  } // end of loop through all cages in cageList
2089 
2090  // Sort the array according to atom ID, which will be needed to get the
2091  // unique IDs and to remove duplicates
2092  sort(atoms.begin(), atoms.end());
2093  // Get the unique atom IDs
2094  auto ip = std::unique(atoms.begin(), atoms.end());
2095  // Resize the vector to remove undefined terms
2096  atoms.resize(std::distance(atoms.begin(), ip));
2097  // If the number of atoms is less than the total nop, add dummy atoms
2098  if (atoms.size() != yCloud->nop) {
2099  padAtoms = true;
2100  bondTypes = 1;
2101  }
2102  // ----------------
2103  // Write output to file inside the output directory
2104  outputFile.open("../output/" + filename);
2105  // FORMAT:
2106  // Comment Line
2107  // 4 atoms
2108  // 4 bonds
2109  // 0 angles
2110  // 0 dihedrals
2111  // 0 impropers
2112  // 1 atom types
2113  // 1 bond types
2114  // 0 angle types
2115  // 0 dihedral types
2116  // 0 improper types
2117  // -1.124000 52.845002 xlo xhi
2118  // 0.000000 54.528999 ylo yhi
2119  // 1.830501 53.087501 zlo zhi
2120 
2121  // Masses
2122 
2123  // 1 15.999400 # O
2124 
2125  // Atoms
2126 
2127  // 1 1 1 0 20.239 1.298 6.873 # O
2128  // 2 1 1 0 0 5.193 6.873 # O
2129  // 3 1 1 0 2.249 1.298 6.873 # O
2130 
2131  // -------
2132  // Write the header
2133  // Write comment line
2134  outputFile << "Written out by D-SEAMS\n";
2135  // Write out the number of atoms
2136  outputFile << yCloud->pts.size() << " "
2137  << "atoms"
2138  << "\n";
2139  // Number of bonds
2140  outputFile << bonds.size() << " bonds"
2141  << "\n";
2142  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
2143  // If padded atoms are required, two atom types will be required
2144  if (padAtoms) {
2145  outputFile << "2 atom types\n";
2146  } else {
2147  outputFile << "1 atom types\n";
2148  } // end of atom types
2149  outputFile
2150  << bondTypes
2151  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
2152  // Box lengths
2153  outputFile << "0 " << yCloud->box[0] << " xlo xhi\n";
2154  outputFile << "0 " << yCloud->box[1] << " ylo yhi\n";
2155  outputFile << "0 " << yCloud->box[2] << " zlo zhi\n";
2156  // Masses
2157  outputFile << "\nMasses\n\n";
2158  // For DDCs and HCs
2159  if (type == cage::cageType::HexC) {
2160  actualCageType = "HC";
2161  } else if (type == cage::cageType::DoubleDiaC) {
2162  actualCageType = "DDC";
2163  } else {
2164  actualCageType = "error";
2165  }
2166  //
2167  outputFile << "1 15.999400 # " << actualCageType << "\n";
2168  if (padAtoms) {
2169  outputFile << "2 1.0 # dummy\n";
2170  }
2171  // Atoms
2172  outputFile << "\nAtoms\n\n";
2173  // -------
2174  // Write out the atom coordinates
2175  // Loop through atoms
2176  for (int i = 0; i < atoms.size(); i++) {
2177  iatom = atoms[i] - 1; // The actual index is one less than the ID
2178  // -----------
2179  // Pad out
2180  // Fill in dummy atoms if some have been skipped
2181  if (atoms[i] != prevAtomID + 1) {
2182  dummyAtoms = atoms[i] - prevAtomID - 1;
2183  dummyID = prevAtomID;
2184  // Loop to write out dummy atoms
2185  for (int j = 0; j < dummyAtoms; j++) {
2186  dummyID++;
2187  jatom = dummyID - 1;
2188  // 1 molecule-tag atom-type q x y z
2189  outputFile << dummyID << " " << yCloud->pts[jatom].molID << " 2 0 "
2190  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
2191  << yCloud->pts[jatom].z << "\n";
2192  } // end of dummy atom write-out
2193  } // end of check for dummy atom printing
2194  // -----------
2195  // Write out coordinates
2196  // 1 molecule-tag atom-type q x y z
2197  outputFile << atoms[i] << " " << yCloud->pts[iatom].molID << " 1 0 "
2198  << yCloud->pts[iatom].x << " " << yCloud->pts[iatom].y << " "
2199  << yCloud->pts[iatom].z << "\n";
2200  // update the previous atom ID
2201  prevAtomID = atoms[i];
2202  } // end of loop through all atoms in atomID
2203 
2204  // Fill in the rest of the dummy atoms
2205  if (atoms[atoms.size() - 1] != yCloud->nop) {
2206  //
2207  for (int id = atoms[atoms.size() - 1] + 1; id <= yCloud->nop; id++) {
2208  jatom = id - 1;
2209  outputFile << id << " " << yCloud->pts[jatom].molID << " 2 0 "
2210  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
2211  << yCloud->pts[jatom].z << "\n";
2212  } // end of printing out dummy atoms
2213  }
2214 
2215  // Print the bonds now!
2216  outputFile << "\nBonds\n\n";
2217  // Loop through all bonds
2218  for (int ibond = 0; ibond < bonds.size(); ibond++) {
2219  // write out the bond
2220  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
2221  << bonds[ibond][1] << "\n";
2222 
2223  } // end of for loop for bonds
2224 
2225  outputFile.close();
2226  // Once the datafile has been printed, exit
2227  return 0;
2228 }
bond::createBondsFromCages
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:529

◆ 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 1793 of file seams_output.cpp.

1797  {
1798  std::ofstream outputFile;
1799  std::vector<int> atoms; // Holds all atom IDs to print
1800  int ringSize = rings[0].size(); // Ring size of each ring in rings
1801  int iatom; // Index, not atom ID
1802  bool padAtoms = false; // Add extra atoms if the atom IDs are skipped
1803  int prevAtomID = 0; // Check for previous atom ID
1804  int dummyAtoms = 0; // Number of dummy atoms to fill
1805  int dummyID;
1806  int jatom; // Array index is 1 less than the ID (index for dummy atom)
1807  int bondTypes = 1;
1808  // Bond stuff
1809  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
1810  // containing the atom IDs of each bond
1811  bool atomOne, atomTwo; // If bond atoms are in the prism or not
1812  bool isPrismBond;
1813 
1814  // ----------------
1815  // Return if there are no prisms
1816  if (listPrism.size() == 0) {
1817  return 1;
1818  }
1819 
1820  // ---------------
1821  // Get the bonds
1822  if (useBondFile) {
1823  // Bonds from file
1824  bonds = sinp::readBonds(bondFile);
1825  } // get bonds from file
1826  else {
1827  bonds = bond::populateBonds(nList, yCloud);
1828  } // Bonds from rings
1829  //
1830  // ----------------
1831  // Otherwise create file
1832  // Create output dir if it doesn't exist already
1833  const char *path = "../output"; // relative to the build directory
1834  fs::path dir(path);
1835  // if (fs::create_directory(dir)) {
1836  // std::cerr << "Output directory created\n";
1837  // }
1838  // ----------------
1839  // Get all the unique atomIDs of the atoms in the rings of this type
1840  // Put all atom IDs into one 1-D vector
1841  size_t total_size{0};
1842  // Get the total number of atoms (repeated)
1843  total_size = listPrism.size() * ringSize;
1844  // Reserve this size inside atoms
1845  atoms.reserve(total_size);
1846  // Fill up all these atom IDs
1847  for (int iring = 0; iring < listPrism.size(); iring++) {
1848  std::move(rings[listPrism[iring]].begin(), rings[listPrism[iring]].end(),
1849  std::back_inserter(atoms));
1850  } // end of loop through all rings in the list
1851 
1852  // Sort the array according to atom ID, which will be needed to get the
1853  // unique IDs and to remove duplicates
1854  sort(atoms.begin(), atoms.end());
1855  // Get the unique atom IDs
1856  auto ip = std::unique(atoms.begin(), atoms.end());
1857  // Resize the vector to remove undefined terms
1858  atoms.resize(std::distance(atoms.begin(), ip));
1859  // If the number of atoms is less than the total nop, add dummy atoms
1860  if (atoms.size() != yCloud->nop) {
1861  padAtoms = true;
1862  bondTypes = 2;
1863  }
1864  // ----------------
1865  // Write output to file inside the output directory
1866  outputFile.open("../output/" + filename);
1867  // FORMAT:
1868  // Comment Line
1869  // 4 atoms
1870  // 4 bonds
1871  // 0 angles
1872  // 0 dihedrals
1873  // 0 impropers
1874  // 1 atom types
1875  // 1 bond types
1876  // 0 angle types
1877  // 0 dihedral types
1878  // 0 improper types
1879  // -1.124000 52.845002 xlo xhi
1880  // 0.000000 54.528999 ylo yhi
1881  // 1.830501 53.087501 zlo zhi
1882 
1883  // Masses
1884 
1885  // 1 15.999400 # O
1886 
1887  // Atoms
1888 
1889  // 1 1 1 0 20.239 1.298 6.873 # O
1890  // 2 1 1 0 0 5.193 6.873 # O
1891  // 3 1 1 0 2.249 1.298 6.873 # O
1892 
1893  // -------
1894  // Write the header
1895  // Write comment line
1896  outputFile << "Written out by D-SEAMS\n";
1897  // Write out the number of atoms
1898  outputFile << yCloud->pts.size() << " "
1899  << "atoms"
1900  << "\n";
1901  // Number of bonds
1902  outputFile << bonds.size() << " bonds"
1903  << "\n";
1904  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
1905  // If padded atoms are required, two atom types will be required
1906  if (padAtoms) {
1907  outputFile << "2 atom types\n";
1908  } else {
1909  outputFile << "1 atom types\n";
1910  } // end of atom types
1911  outputFile
1912  << bondTypes
1913  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
1914  // Box lengths
1915  outputFile << "0 " << yCloud->box[0] << " xlo xhi\n";
1916  outputFile << "0 " << yCloud->box[1] << " ylo yhi\n";
1917  outputFile << "0 " << yCloud->box[2] << " zlo zhi\n";
1918  // Masses
1919  outputFile << "\nMasses\n\n";
1920  outputFile << "1 15.999400 # O\n";
1921  if (padAtoms) {
1922  outputFile << "2 1.0 # dummy\n";
1923  }
1924  // Atoms
1925  outputFile << "\nAtoms\n\n";
1926  // -------
1927  // Write out the atom coordinates
1928  // Loop through atoms
1929  for (int i = 0; i < atoms.size(); i++) {
1930  iatom = atoms[i] - 1; // The actual index is one less than the ID
1931  // -----------
1932  // Pad out
1933  // Fill in dummy atoms if some have been skipped
1934  if (atoms[i] != prevAtomID + 1) {
1935  dummyAtoms = atoms[i] - prevAtomID - 1;
1936  dummyID = prevAtomID;
1937  // Loop to write out dummy atoms
1938  for (int j = 0; j < dummyAtoms; j++) {
1939  dummyID++;
1940  jatom = dummyID - 1;
1941  // 1 molecule-tag atom-type q x y z
1942  outputFile << dummyID << " " << yCloud->pts[jatom].molID << " 2 0 "
1943  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
1944  << yCloud->pts[jatom].z << "\n";
1945  } // end of dummy atom write-out
1946  } // end of check for dummy atom printing
1947  // -----------
1948  // Write out coordinates
1949  // 1 molecule-tag atom-type q x y z
1950  outputFile << atoms[i] << " " << yCloud->pts[iatom].molID << " 1 0 "
1951  << yCloud->pts[iatom].x << " " << yCloud->pts[iatom].y << " "
1952  << yCloud->pts[iatom].z << "\n";
1953  // update the previous atom ID
1954  prevAtomID = atoms[i];
1955  } // end of loop through all atoms in atomID
1956 
1957  // Fill in the rest of the dummy atoms
1958  if (atoms[atoms.size() - 1] != yCloud->nop) {
1959  //
1960  for (int id = atoms[atoms.size() - 1] + 1; id <= yCloud->nop; id++) {
1961  jatom = id - 1;
1962  outputFile << id << " " << yCloud->pts[jatom].molID << " 2 0 "
1963  << yCloud->pts[jatom].x << " " << yCloud->pts[jatom].y << " "
1964  << yCloud->pts[jatom].z << "\n";
1965  } // end of printing out dummy atoms
1966  }
1967 
1968  // Print the bonds now!
1969  outputFile << "\nBonds\n\n";
1970  // Loop through all bonds
1971  for (int ibond = 0; ibond < bonds.size(); ibond++) {
1972  // Init
1973  isPrismBond = false;
1974  atomOne = false;
1975  atomTwo = false;
1976  // --------
1977  // Check if the bond is in the prism or not
1978  auto it = std::find(atoms.begin() + 1, atoms.end(), bonds[ibond][0]);
1979  if (it != atoms.end()) {
1980  atomOne = true;
1981  } else if (bonds[ibond][0] == atoms[0]) {
1982  atomOne = true;
1983  } else if (bonds[ibond][0] == atoms[atoms.size() - 1]) {
1984  atomOne = true;
1985  }
1986 
1987  auto it1 = std::find(atoms.begin() + 1, atoms.end(), bonds[ibond][1]);
1988  if (it1 != atoms.end()) {
1989  atomTwo = true;
1990  } else if (bonds[ibond][1] == atoms[0]) {
1991  atomTwo = true;
1992  } else if (bonds[ibond][1] == atoms[atoms.size() - 1]) {
1993  atomTwo = true;
1994  }
1995 
1996  if (atomOne == false || atomTwo == false) {
1997  isPrismBond = false;
1998  } else {
1999  isPrismBond = true;
2000  }
2001  // --------
2002  if (isPrismBond) {
2003  // is inside the prism (type 1)
2004  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
2005  << bonds[ibond][1] << "\n";
2006  } else {
2007  // not inside the prism (type 2)
2008  outputFile << ibond + 1 << " 2 " << bonds[ibond][0] << " "
2009  << bonds[ibond][1] << "\n";
2010  }
2011 
2012  } // end of for loop for bonds
2013 
2014  outputFile.close();
2015  // Once the datafile has been printed, exit
2016  return 0;
2017 }
sinp::readBonds
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:24

◆ 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 2397 of file seams_output.cpp.

2400  {
2401  //
2402  std::ofstream outputFile;
2403  int iatom; // Index, not atom ID
2404  int currentAtomType; // Current atom type: a value from 1 to 4
2405  int numAtomTypes = 6; // DDC, HC, Mixed, dummy, mixed2 and pnc
2406  int bondTypes = 1;
2407  // Bond stuff
2408  std::vector<std::vector<int>> bonds; // Vector of vector, with each row
2409  // containing the atom IDs of each bond
2410  std::string filename =
2411  "system-" + std::to_string(yCloud->currentFrame) + ".data";
2412 
2413  // ---------------
2414  // Get the bonds
2415  if (bondsBetweenDummy) {
2416  bonds = bond::populateBonds(nList, yCloud);
2417  } // create bonds between dummy atoms
2418  else {
2419  bonds = bond::populateBonds(nList, yCloud, atomTypes);
2420  } // only create bonds between non-dummy atoms
2421  //
2422  // ----------------
2423  // Make the output directory if it doesn't exist
2424  sout::makePath(path);
2425  std::string outputDirName = path + "bulkTopo";
2426  sout::makePath(outputDirName);
2427  outputDirName = path + "bulkTopo/dataFiles/";
2428  sout::makePath(outputDirName);
2429  // ----------------
2430  // Write output to file inside the output directory
2431  outputFile.open(path + "bulkTopo/dataFiles/" + filename);
2432  // FORMAT:
2433  // Comment Line
2434  // 4 atoms
2435  // 4 bonds
2436  // 0 angles
2437  // 0 dihedrals
2438  // 0 impropers
2439  // 1 atom types
2440  // 1 bond types
2441  // 0 angle types
2442  // 0 dihedral types
2443  // 0 improper types
2444  // -1.124000 52.845002 xlo xhi
2445  // 0.000000 54.528999 ylo yhi
2446  // 1.830501 53.087501 zlo zhi
2447 
2448  // Masses
2449 
2450  // 1 15.999400 # O
2451 
2452  // Atoms
2453 
2454  // 1 1 1 0 20.239 1.298 6.873 # O
2455  // 2 1 1 0 0 5.193 6.873 # O
2456  // 3 1 1 0 2.249 1.298 6.873 # O
2457 
2458  // -------
2459  // Write the header
2460  // Write comment line
2461  outputFile << "Written out by D-SEAMS\n";
2462  // Write out the number of atoms
2463  outputFile << yCloud->pts.size() << " "
2464  << "atoms"
2465  << "\n";
2466  // Number of bonds
2467  outputFile << bonds.size() << " bonds"
2468  << "\n";
2469  outputFile << "0 angles\n0 dihedrals\n0 impropers\n";
2470  // There are maxDepth-2 total types of prisms + dummy
2471  outputFile << numAtomTypes << " atom types\n";
2472  // Bond types
2473  outputFile
2474  << bondTypes
2475  << " bond types\n0 angle types\n0 dihedral types\n0 improper types\n";
2476  // Box lengths
2477  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
2478  << " xlo xhi\n";
2479  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
2480  << " ylo yhi\n";
2481  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
2482  << " zlo zhi\n";
2483  // Masses
2484  outputFile << "\nMasses\n\n";
2485  outputFile << "1 15.999400 # dummy\n";
2486  outputFile << "2 15.999400 # hc \n";
2487  outputFile << "3 15.999400 # ddc \n";
2488  outputFile << "4 15.999400 # mixed \n";
2489  outputFile << "5 15.999400 # pnc \n";
2490  outputFile << "6 15.999400 # pncHexaMixed \n";
2491  // Atoms
2492  outputFile << "\nAtoms\n\n";
2493  // -------
2494  // Write out the atom coordinates
2495  // Loop through atoms
2496  for (int i = 0; i < yCloud->pts.size(); i++) {
2497  iatom =
2498  yCloud->pts[i].atomID; // The actual ID can be different from the index
2499  //
2500  // Get the atom type
2501  // hc atom type
2502  if (atomTypes[i] == cage::iceType::hc) {
2503  currentAtomType = 2;
2504  } // hc
2505  else if (atomTypes[i] == cage::iceType::ddc) {
2506  currentAtomType = 3;
2507  } // ddc
2508  // mixed
2509  else if (atomTypes[i] == cage::iceType::mixed) {
2510  currentAtomType = 4;
2511  } // mixed
2512  // pnc
2513  else if (atomTypes[i] == cage::iceType::pnc) {
2514  currentAtomType = 5;
2515  } // mixed
2516  // pnc and DDCs/HCs mixed
2517  else if (atomTypes[i] == cage::iceType::mixed2) {
2518  currentAtomType = 6;
2519  } // mixed
2520  // dummy
2521  else {
2522  currentAtomType = 1;
2523  } // dummy
2524  //
2525  // Write out coordinates
2526  // atomID molecule-tag atom-type q x y z
2527  outputFile << iatom << " " << yCloud->pts[i].molID << " " << currentAtomType
2528  << " 0 " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
2529  << yCloud->pts[i].z << "\n";
2530 
2531  } // end of loop through all atoms in pointCloud
2532 
2533  // Print the bonds now!
2534  outputFile << "\nBonds\n\n";
2535  // Loop through all bonds
2536  for (int ibond = 0; ibond < bonds.size(); ibond++) {
2537  //
2538  outputFile << ibond + 1 << " 1 " << bonds[ibond][0] << " "
2539  << bonds[ibond][1] << "\n";
2540  } // end of for loop for bonds
2541 
2542  // Once the datafile has been printed, exit
2543  return 0;
2544 }
cage::iceType::mixed
@ mixed
cage::iceType::mixed2
@ mixed2
cage::iceType::hc
@ hc
cage::iceType::pnc
@ pnc
cage::iceType::ddc
@ ddc

◆ 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 1288 of file seams_output.cpp.

1291  {
1292  std::ofstream outputFile;
1293  int iatom; // Index, not atom ID
1294  std::string filename =
1295  "dump-" + std::to_string(yCloud->currentFrame) + ".lammpstrj";
1296  // ----------------
1297  // Make the output directory if it doesn't exist
1298  std::string outputDirName = path + "bulkTopo/dumpFiles";
1299  sout::makePath(outputDirName);
1300  // ----------------
1301  // Write out information about the data types
1302  if (yCloud->currentFrame == firstFrame) {
1303  outputFile.open(path + "bulkTopo/typeInfo.dat");
1304  outputFile << "Atom types in the dump files are:\n";
1305  outputFile << " Type 0 (dummy) = unidentified phase\n";
1306  outputFile << " Type 1 (hc) = atom belonging to a Hexagonal Cage.\n";
1307  outputFile << " Type 2 (ddc) = atom belonging to a Double-Diamond Cage\n";
1308  outputFile << " Type 3 (mixed) = atom belonging to a mixed ring shared by "
1309  "a DDC and HC\n";
1310  outputFile
1311  << " Type 4 (pnc) = atom belonging to a pair of pentagonal rings\n";
1312  outputFile << " Type 5 (mixed2) = atom belonging to a pentagonal "
1313  "nanochannel, shared by DDCs/HCs\n";
1314  outputFile.close();
1315  } // end of writing out information
1316  // ----------------
1317  // Write output to file inside the output directory
1318  outputFile.open(path + "bulkTopo/dumpFiles/" + filename);
1319  // ----------------------------------------------------
1320  // Header Format
1321 
1322  // ITEM: TIMESTEP
1323  // 0
1324  // ITEM: NUMBER OF ATOMS
1325  // 500
1326  // ITEM: BOX BOUNDS pp pp pp
1327  // -9.0400100000000005e-01 1.7170999999999999e+01
1328  // -9.0400100000000005e-01 1.7170999999999999e+01
1329  // -9.0400100000000005e-01 1.7170999999999999e+01
1330  // ITEM: ATOMS id mol type x y z rmsd
1331 
1332  // -----------------
1333  // -------
1334  // Write the header
1335  // ITEM: TIMESTEP
1336  outputFile << "ITEM: TIMESTEP\n";
1337  // Write out frame number
1338  outputFile << yCloud->currentFrame << "\n";
1339  // ITEM: NUMBER OF ATOMS
1340  outputFile << "ITEM: NUMBER OF ATOMS\n";
1341  // Number of atoms
1342  outputFile << yCloud->pts.size() << "\n";
1343  // ITEM: BOX BOUNDS pp pp pp
1344  outputFile << "ITEM: BOX BOUNDS pp pp pp\n";
1345  // Box lengths
1346  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1347  << "\n";
1348  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1349  << "\n";
1350  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1351  << "\n";
1352  // ITEM: ATOMS id mol type x y z rmsd
1353  outputFile << "ITEM: ATOMS id mol type x y z rmsd\n";
1354  // -------
1355  // Write out the atom coordinates
1356  // Format
1357  // ITEM: ATOMS id mol type x y z rmsd
1358  //
1359  // Loop through atoms
1360  for (int i = 0; i < yCloud->pts.size(); i++) {
1361  iatom =
1362  yCloud->pts[i].atomID; // The actual ID can be different from the index
1363  // Write out coordinates
1364  outputFile << iatom << " " << yCloud->pts[i].molID << " " << atomTypes[i]
1365  << " " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
1366  << yCloud->pts[i].z << " " << rmsdPerAtom[i] << "\n";
1367 
1368  } // end of loop through all atoms in pointCloud
1369  // -----------------------------------------------------
1370  outputFile.close(); // Close the file
1371  return 0;
1372 } // 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 1378 of file seams_output.cpp.

1381  {
1382  //
1383  std::ofstream outputFile;
1384  int iatom; // Index, not atom ID
1385  std::string filename =
1386  "dump-" + std::to_string(yCloud->currentFrame) + ".lammpstrj";
1387  // ----------------
1388  // Make the output directory if it doesn't exist
1389  std::string outputDirName = path + "topoINT/dumpFiles";
1390  sout::makePath(outputDirName);
1391  // ----------------
1392  // Write output to file inside the output directory
1393  outputFile.open(path + "topoINT/dumpFiles/" + filename);
1394  // ----------------------------------------------------
1395  // Header Format
1396 
1397  // ITEM: TIMESTEP
1398  // 0
1399  // ITEM: NUMBER OF ATOMS
1400  // 500
1401  // ITEM: BOX BOUNDS pp pp pp
1402  // -9.0400100000000005e-01 1.7170999999999999e+01
1403  // -9.0400100000000005e-01 1.7170999999999999e+01
1404  // -9.0400100000000005e-01 1.7170999999999999e+01
1405  // ITEM: ATOMS id mol type x y z rmsd
1406 
1407  // -----------------
1408  // -------
1409  // Write the header
1410  // ITEM: TIMESTEP
1411  outputFile << "ITEM: TIMESTEP\n";
1412  // Write out frame number
1413  outputFile << yCloud->currentFrame << "\n";
1414  // ITEM: NUMBER OF ATOMS
1415  outputFile << "ITEM: NUMBER OF ATOMS\n";
1416  // Number of atoms
1417  outputFile << yCloud->pts.size() << "\n";
1418  // ITEM: BOX BOUNDS pp pp pp
1419  outputFile << "ITEM: BOX BOUNDS pp pp pp\n";
1420  // Box lengths
1421  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1422  << "\n";
1423  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1424  << "\n";
1425  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1426  << "\n";
1427  // ITEM: ATOMS id mol type x y z rmsd
1428  outputFile << "ITEM: ATOMS id mol type x y z rmsd\n";
1429  // -------
1430  // Write out the atom coordinates
1431  // Format
1432  // ITEM: ATOMS id mol type x y z rmsd
1433  //
1434  // Loop through atoms
1435  for (int i = 0; i < yCloud->pts.size(); i++) {
1436  iatom =
1437  yCloud->pts[i].atomID; // The actual ID can be different from the index
1438  // Write out coordinates
1439  outputFile << iatom << " " << yCloud->pts[i].molID << " " << atomTypes[i]
1440  << " " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
1441  << yCloud->pts[i].z << " " << rmsdPerAtom[i] << "\n";
1442 
1443  } // end of loop through all atoms in pointCloud
1444  // -----------------------------------------------------
1445  return 0;
1446 } // end of function

◆ writeLAMMPSdumpSlice()

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

Write out a LAMMPS dump file containing the inSlice value for each atom for a user-defined slice

Prints out a LAMMPS dump file for all atoms, for every frame, printing the inSlice attribute for a user-defined slice in a separate column

Definition at line 1452 of file seams_output.cpp.

1454  {
1455  //
1456  std::ofstream outputFile;
1457  int iatom; // Index, not atom ID
1458  std::string filename =
1459  "dump-" + std::to_string(yCloud->currentFrame) + ".lammpstrj";
1460  // ----------------
1461  // Make the output directory if it doesn't exist
1462  sout::makePath(path+"selection");
1463  std::string outputDirName = path + "selection/dumpFiles";
1464  sout::makePath(outputDirName);
1465  // ----------------
1466  // Write output to file inside the output directory
1467  outputFile.open(path + "selection/dumpFiles/" + filename);
1468  // ----------------------------------------------------
1469  // Header Format
1470 
1471  // ITEM: TIMESTEP
1472  // 0
1473  // ITEM: NUMBER OF ATOMS
1474  // 500
1475  // ITEM: BOX BOUNDS pp pp pp
1476  // -9.0400100000000005e-01 1.7170999999999999e+01
1477  // -9.0400100000000005e-01 1.7170999999999999e+01
1478  // -9.0400100000000005e-01 1.7170999999999999e+01
1479  // ITEM: ATOMS id mol type x y z rmsd
1480 
1481  // -----------------
1482  // -------
1483  // Write the header
1484  // ITEM: TIMESTEP
1485  outputFile << "ITEM: TIMESTEP\n";
1486  // Write out frame number
1487  outputFile << yCloud->currentFrame << "\n";
1488  // ITEM: NUMBER OF ATOMS
1489  outputFile << "ITEM: NUMBER OF ATOMS\n";
1490  // Number of atoms
1491  outputFile << yCloud->pts.size() << "\n";
1492  // ITEM: BOX BOUNDS pp pp pp
1493  outputFile << "ITEM: BOX BOUNDS pp pp pp\n";
1494  // Box lengths
1495  outputFile << yCloud->boxLow[0] << " " << yCloud->boxLow[0] + yCloud->box[0]
1496  << "\n";
1497  outputFile << yCloud->boxLow[1] << " " << yCloud->boxLow[1] + yCloud->box[1]
1498  << "\n";
1499  outputFile << yCloud->boxLow[2] << " " << yCloud->boxLow[2] + yCloud->box[2]
1500  << "\n";
1501  // ITEM: ATOMS id mol type x y z rmsd
1502  outputFile << "ITEM: ATOMS id mol type x y z inSlice\n";
1503  // -------
1504  // Write out the atom coordinates
1505  // Format
1506  // ITEM: ATOMS id mol type x y z rmsd
1507  //
1508  // Loop through atoms
1509  for (int i = 0; i < yCloud->pts.size(); i++) {
1510  iatom =
1511  yCloud->pts[i].atomID; // The actual ID can be different from the index
1512  // Write out coordinates
1513  outputFile << iatom << " " << yCloud->pts[i].molID << " " << yCloud->pts[i].type
1514  << " " << yCloud->pts[i].x << " " << yCloud->pts[i].y << " "
1515  << yCloud->pts[i].z << " " << yCloud->pts[i].inSlice << "\n";
1516 
1517  } // end of loop through all atoms in pointCloud
1518  // -----------------------------------------------------
1519  return 0;
1520 } // end of function

◆ writeMoleculeIDsExpressionSelectOVITO()

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

Function for printing out the molecule IDs present in the slice (compatible with the OVITO Expression Select command

Function for printing out the molecule IDs present in the slice. The format should be compatible with the group command in LAMMPS

Definition at line 964 of file seams_output.cpp.

965  {
966  std::ofstream outputFile;
967  std::string filename =
968  "ovito-molIDSelect-" + std::to_string(yCloud->currentFrame) + ".dat";
969  std::vector<int> idVec; // Vector which will contain all molecule IDs present in the slice
970  // Eventually we will sort this and only keep unique molecule IDs.
971  int prevElem, currentElem; // previous and current mol ID in the sequence
972  int lastPrintedElemSeries; // Last element printed
973  // ----------------
974  // Make the output directory if it doesn't exist
975  sout::makePath(path);
976  sout::makePath(path+"selection");
977  sout::makePath(path+"selection/IDovitoFiles");
978  // ----------------
979  // Write output to file inside the output directory
980  outputFile.open(path + "selection/IDovitoFiles/" + filename);
981 
982  // ----------------
983  // Find all molecule IDs from yCloud
984  // Loop through all iatom in yCloud
985  for (int iatom = 0; iatom < yCloud->nop; iatom++)
986  {
987  // If iatom is in the slice, add the molecule ID to idVec
988  if (yCloud->pts[iatom].inSlice)
989  {
990  idVec.push_back(yCloud->pts[iatom].molID);
991  } // end of adding molecule ID to vector for slice
992  } // end of loop through iatom
993  // ----------------
994  // Sort in ascending order
995  std::sort(idVec.begin(), idVec.end());
996  auto it = std::unique(idVec.begin(), idVec.end());
997  // Get rid of undefined elements
998  idVec.resize(std::distance(idVec.begin(), it));
999  // ----------------
1000 
1001  // Format:
1002  // Comment line
1003  // 1 2 3 4 6 7
1004 
1005  // ----------------
1006  // Comment line
1007  outputFile << "# Molecule IDs in slice\n";
1008  outputFile << "# OVITO Expression select command \n";
1009  // ----------------
1010 
1011  // Print other molecule IDs to the file
1012  for (int i=0; i<idVec.size()-1; i++)
1013  {
1014  currentElem = idVec[i]; // current mol ID
1015 
1016  outputFile << "MoleculeIdentifier == " << currentElem << " || ";
1017 
1018  } // end of printing all elements to file except the last one
1019 
1020  // Print the last element
1021  outputFile << "MoleculeIdentifier == " << idVec.back();
1022 
1023  outputFile.close();
1024 
1025  return 0;
1026 }
std::vector::back
T back(T... args)

◆ writeMoleculeIDsInSlice()

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

Function for printing out the molecule IDs present in the slice (compatible with the LAMMPS group command

Function for printing out the molecule IDs present in the slice. The format should be compatible with the group command in LAMMPS

Definition at line 877 of file seams_output.cpp.

878  {
879  std::ofstream outputFile;
880  std::string filename =
881  "molID-" + std::to_string(yCloud->currentFrame) + ".dat";
882  std::vector<int> idVec; // Vector which will contain all molecule IDs present in the slice
883  // Eventually we will sort this and only keep unique molecule IDs.
884  int prevElem, currentElem; // previous and current mol ID in the sequence
885  int lastPrintedElemSeries; // Last element printed
886  // ----------------
887  // Make the output directory if it doesn't exist
888  sout::makePath(path);
889  sout::makePath(path+"selection");
890  sout::makePath(path+"selection/IDtextFiles");
891  // ----------------
892  // Write output to file inside the output directory
893  outputFile.open(path + "selection/IDtextFiles/" + filename);
894 
895  // ----------------
896  // Find all molecule IDs from yCloud
897  // Loop through all iatom in yCloud
898  for (int iatom = 0; iatom < yCloud->nop; iatom++)
899  {
900  // If iatom is in the slice, add the molecule ID to idVec
901  if (yCloud->pts[iatom].inSlice)
902  {
903  idVec.push_back(yCloud->pts[iatom].molID);
904  } // end of adding molecule ID to vector for slice
905  } // end of loop through iatom
906  // ----------------
907  // Sort in ascending order
908  std::sort(idVec.begin(), idVec.end());
909  auto it = std::unique(idVec.begin(), idVec.end());
910  // Get rid of undefined elements
911  idVec.resize(std::distance(idVec.begin(), it));
912  // ----------------
913 
914  // Format:
915  // Comment line
916  // 1 2 3 4 6 7
917 
918  // ----------------
919  // Comment line
920  outputFile << "# Molecule IDs in slice\n";
921  outputFile << "# LAMMPS command : group groupName molecule 100:10000 \n";
922  // Format
923  // In LAMMPS, groups can be assigned by ID using the following command
924  // group groupName molecule 100:10000
925  // ----------------
926  // First element
927  outputFile << idVec[0];
928  prevElem = idVec[0];
929  lastPrintedElemSeries = idVec[0];
930  // ----------------
931 
932  // Print other molecule IDs to the file
933  for (int i=1; i<idVec.size(); i++)
934  {
935  currentElem = idVec[i]; // current mol ID
936  // prevElem is the previous mol ID
937  // if the currentElem-prevElem>1
938  if (currentElem-prevElem>1 || i==idVec.size()-1)
939  {
940  if (lastPrintedElemSeries!=prevElem)
941  {
942  outputFile << ":" << prevElem << " " << currentElem;
943  lastPrintedElemSeries = currentElem;
944  } // the previous element has not been printed
945  else{
946  outputFile << " " << currentElem;
947  lastPrintedElemSeries = currentElem;
948  } // the previous element has already been printed
949  //
950  } // print currentElem
951  //
952  prevElem = currentElem;
953  } // end of printing all elements to file
954 
955  outputFile.close();
956 
957  return 0;
958 }

◆ 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 1031 of file seams_output.cpp.

1034  {
1035  std::ofstream outputFile;
1036  int totalPrisms; // Number of total prisms
1037  // ----------------
1038  // Make the output directory if it doesn't exist
1039  sout::makePath(path);
1040  std::string outputDirName = path + "topoINT";
1041  sout::makePath(outputDirName);
1042  // ----------------
1043  // Write output to file inside the output directory
1044  outputFile.open(path + "topoINT/nPrisms.dat",
1045  std::ios_base::app | std::ios_base::out);
1046 
1047  // ----------------
1048  // Write the comment line if the first frame is being written out
1049  if (currentFrame == firstFrame) {
1050  outputFile << "Frame RingSize Num_of_prisms Height% RingSize ... Height\n";
1051  }
1052  // ----------------
1053  // Format:
1054  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
1055  // 1 3 0 0 4 35 40 ....
1056 
1057  outputFile << currentFrame << " ";
1058 
1059  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1060  totalPrisms = nPrisms[ringSize - 3] + nDefPrisms[ringSize - 3];
1061  // Write out
1062  outputFile << ringSize << " " << totalPrisms << " "
1063  << nDefPrisms[ringSize - 3] << " " << heightPercent[ringSize - 3]
1064  << " ";
1065  }
1066 
1067  outputFile << "\n";
1068 
1069  outputFile.close();
1070 
1071  return 0;
1072 }

◆ 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 216 of file seams_output.cpp.

218  {
219  std::ofstream outputFile;
220  std::string number = std::to_string(prismNum);
221  std::string filename = "prism" + number + ".dat";
222  int ringSize =
223  (*basal1).size(); // Size of the ring; each ring contains n elements
224  int iatomIndex; // index of atom coordinate being written out
225 
226  // ----------------
227  // Create output dir if it doesn't exist already
228  const char *path = "../output/prisms"; // relative to the build directory
229  fs::path dir(path);
230  // if (fs::create_directory(dir)) {
231  // std::cerr << "Output Prism directory created\n";
232  // }
233  // ----------------
234  // Write output to file inside the output directory
235  outputFile.open("../output/prisms/" + filename);
236 
237  // Format:
238  // x y z coordinates of each node
239 
240  // For basal 1
241  for (int iring = 0; iring < ringSize; iring++) {
242  iatomIndex = (*basal1)[iring]; // C++ indices are one less
243  // Write the coordinates out to the file
244  outputFile << yCloud->pts[iatomIndex].x << " ";
245  outputFile << yCloud->pts[iatomIndex].y << " ";
246  outputFile << yCloud->pts[iatomIndex].z << " ";
247 
248  outputFile << "\n";
249  } // end of loop through basal1
250 
251  // For basal 2
252  for (int iring = 0; iring < ringSize; iring++) {
253  iatomIndex = (*basal2)[iring]; // C++ indices are one less
254  // Write the coordinates out to the file
255  outputFile << yCloud->pts[iatomIndex].x << " ";
256  outputFile << yCloud->pts[iatomIndex].y << " ";
257  outputFile << yCloud->pts[iatomIndex].z << " ";
258 
259  outputFile << "\n";
260  } // end of loop through basal1
261 
262  outputFile.close();
263 
264  // ---- Print out all the coordinates of a single ring, for prismNum=1 only
265  if (prismNum == 1) {
266  outputFile.open("../output/prisms/singleRing.dat");
267  // For basal 1
268  for (int iring = 0; iring < ringSize; iring++) {
269  iatomIndex = (*basal1)[iring]; // C++ indices are one less
270  // Write the coordinates out to the file
271  outputFile << yCloud->pts[iatomIndex].x << " ";
272  outputFile << yCloud->pts[iatomIndex].y << " ";
273  outputFile << yCloud->pts[iatomIndex].z << " ";
274 
275  outputFile << "\n";
276  } // end of loop through basal1
277  outputFile.close();
278  }
279 
280  return 0;
281 }

◆ 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 1077 of file seams_output.cpp.

1082  {
1083  std::ofstream outputFileXY;
1084  std::ofstream outputFileXZ;
1085  std::ofstream outputFileYZ;
1086  // ----------------
1087  // Make the output directory if it doesn't exist
1088  sout::makePath(path);
1089  std::string outputDirName = path + "topoMonolayer";
1090  sout::makePath(outputDirName);
1091  // ----------------
1092  // Coverage Area of XY
1093  // Write output to file inside the output directory
1094  outputFileXY.open(path + "topoMonolayer/coverageAreaXY.dat",
1095  std::ios_base::app | std::ios_base::out);
1096 
1097  // Format:
1098  // Comment line
1099  // 1 3 0 0 4 35 40 ....
1100 
1101  // ----------------
1102  // Add comment for the first frame
1103  if (currentFrame == firstFrame) {
1104  outputFileXY << "Frame RingSize Num_of_rings CoverageAreaXY% RingSize ... "
1105  "CoverageAreaXY%\n";
1106  }
1107  // ----------------
1108 
1109  outputFileXY << currentFrame << " ";
1110 
1111  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1112  outputFileXY << ringSize << " " << nRings[ringSize - 3] << " "
1113  << coverageAreaXY[ringSize - 3] << " ";
1114  }
1115 
1116  outputFileXY << "\n";
1117 
1118  outputFileXY.close();
1119  // ----------------
1120  // Coverage Area of XZ
1121  // Write output to file inside the output directory
1122  outputFileXZ.open(path + "topoMonolayer/coverageAreaXZ.dat",
1123  std::ios_base::app | std::ios_base::out);
1124 
1125  // ----------------
1126  // Add comment for the first frame
1127  if (currentFrame == firstFrame) {
1128  outputFileXZ << "Frame RingSize Num_of_rings CoverageAreaXZ% RingSize ... "
1129  "CoverageAreaXZ%\n";
1130  }
1131  // ----------------
1132 
1133  // Format:
1134  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
1135  // 1 3 0 0 4 35 40 ....
1136 
1137  outputFileXZ << currentFrame << " ";
1138 
1139  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1140  outputFileXZ << ringSize << " " << nRings[ringSize - 3] << " "
1141  << coverageAreaXZ[ringSize - 3] << " ";
1142  }
1143 
1144  outputFileXZ << "\n";
1145 
1146  outputFileXZ.close();
1147  // ----------------
1148  // Coverage Area of YZ
1149  // Write output to file inside the output directory
1150  outputFileYZ.open(path + "topoMonolayer/coverageAreaYZ.dat",
1151  std::ios_base::app | std::ios_base::out);
1152 
1153  // ----------------
1154  // Add comment for the first frame
1155  if (currentFrame == firstFrame) {
1156  outputFileYZ << "Frame RingSize Num_of_rings CoverageAreaYZ% RingSize ... "
1157  "CoverageAreaYZ%\n";
1158  }
1159  // ----------------
1160 
1161  // Format:
1162  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
1163  // 1 3 0 0 4 35 40 ....
1164 
1165  outputFileYZ << currentFrame << " ";
1166 
1167  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1168  outputFileYZ << ringSize << " " << nRings[ringSize - 3] << " "
1169  << coverageAreaYZ[ringSize - 3] << " ";
1170  }
1171 
1172  outputFileYZ << "\n";
1173 
1174  outputFileYZ.close();
1175 
1176  return 0;
1177 }

◆ writeRingNumBulk()

int sout::writeRingNumBulk ( std::string  path,
int  currentFrame,
std::vector< int >  nRings,
int  maxDepth,
int  firstFrame 
)

Function for printing out the number of rings of each type in a bulk system

Function for printing out ring info, for a bulk system

Definition at line 1182 of file seams_output.cpp.

1185  {
1186  std::ofstream outputFile;
1187  // ----------------
1188  // Make the output directory if it doesn't exist
1189  sout::makePath(path);
1190  std::string outputDirName = path + "bulkTopo";
1191  sout::makePath(outputDirName);
1192  // ----------------
1193  // Ring output file
1194  // Write output to file inside the output directory
1195  outputFile.open(path + "bulkTopo/num_rings.dat",
1196  std::ios_base::app | std::ios_base::out);
1197 
1198  // Format:
1199  // Comment line
1200  // 1 3 0 4 35 ....
1201 
1202  // ----------------
1203  // Add comment for the first frame
1204  if (currentFrame == firstFrame) {
1205  outputFile << "Frame RingSize Num_of_rings RingSize Num_of_rings...\n";
1206  }
1207  // ----------------
1208 
1209  outputFile << currentFrame << " ";
1210 
1211  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
1212  outputFile << ringSize << " " << nRings[ringSize - 3] << " ";
1213  }
1214 
1215  outputFile << "\n";
1216 
1217  outputFile.close();
1218 
1219  return 0;
1220 }

◆ 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 182 of file seams_output.cpp.

183  {
184  std::ofstream outputFile;
185  // ----------------
186  // Create output dir if it doesn't exist already
187  const char *path = "../output"; // relative to the build directory
188  fs::path dir(path);
189  // if (fs::create_directory(dir)) {
190  // std::cerr << "Output directory created\n";
191  // }
192  // ----------------
193  // Write output to file inside the output directory
194  outputFile.open("../output/" + filename);
195 
196  // Format:
197  // 272 214 906 1361 388 1
198 
199  for (int iring = 0; iring < rings.size(); iring++) {
200  // Otherwise, write out to the file
201  for (int k = 0; k < rings[iring].size(); k++) {
202  outputFile << rings[iring][k] << " ";
203  } // end of loop through ring elements
204  outputFile << "\n";
205  } // end of loop through rings
206 
207  outputFile.close();
208 
209  return 0;
210 }

◆ 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 1251 of file seams_output.cpp.

1253  {
1254  //
1255  std::ofstream outputFile;
1256  // ----------------
1257  // Make the output directory if it doesn't exist
1258  sout::makePath(path);
1259  std::string outputDirName = path + "bulkTopo";
1260  sout::makePath(outputDirName);
1261  // ----------------
1262  // Write output to file inside the output directory
1263  outputFile.open(path + "bulkTopo/cageData.dat",
1264  std::ios_base::app | std::ios_base::out);
1265 
1266  // Format:
1267  // Frame RingSize Num_of_prisms Height% RingSize ... Height%
1268  // 1 3 0 0 4 35 40 ....
1269  // -------------------
1270  // If first frame then write the comment line
1271  if (currentFrame == firstFrame) {
1272  outputFile << "Frame HCnumber DDCnumber MixedRingNumber PrismaticRings "
1273  "basalRings\n";
1274  }
1275  // -------------------
1276  outputFile << currentFrame << " " << numHC << " " << numDDC << " "
1277  << mixedRings << " " << prismaticRings << " " << basalRings
1278  << "\n";
1279 
1280  outputFile.close();
1281  return 0;
1282 } // 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 2550 of file seams_output.cpp.

2552  {
2553  //
2554  std::ofstream outputFile;
2555  std::string filename = "cluster-" + std::to_string(clusterID) + ".xyz";
2556  int nAtoms = atoms.size(); // Number of atoms
2557  int iatom; // Current atom
2558  // ----------------
2559  // Make the output directory if it doesn't exist
2560  sout::makePath(path);
2561  std::string outputDirName = path + "bulkTopo";
2562  sout::makePath(outputDirName);
2563  outputDirName = path + "bulkTopo/clusterXYZ/";
2564  sout::makePath(outputDirName);
2565  outputDirName = path + "bulkTopo/clusterXYZ/frame-" +
2566  std::to_string(yCloud->currentFrame);
2567  sout::makePath(outputDirName);
2568  // ----------------
2569  // Write output to file inside the output directory
2570  outputFile.open(path + "bulkTopo/clusterXYZ/frame-" +
2571  std::to_string(yCloud->currentFrame) + "/" + filename);
2572 
2573  // Format of an XYZ file:
2574  // 1970
2575  // generated by VMD
2576  // O 43.603500 16.926201 15.215700
2577  // O 39.912601 14.775100 19.379200
2578  outputFile << nAtoms << "\n"; // Number of atoms
2579  outputFile
2580  << "Generated by d-SEAMS. 0 type=hc and 1 type =ddc\n"; // Comment line
2581  //
2582  // Write out all the atom coordinates
2583  for (int i = 0; i < nAtoms; i++) {
2584  iatom = atoms[i]; // Atom index to be printed out
2585  // TODO: Should print string representation
2586  outputFile << static_cast<int>(type) << " " << yCloud->pts[iatom].x << " "
2587  << yCloud->pts[iatom].y << " " << yCloud->pts[iatom].z << "\n";
2588  } // end of loop through all atoms
2589 
2590  outputFile.close();
2591 
2592  return 0;
2593 }
Generated by Doxygen 1.9.1 and the doxyYoda theme