seams_input.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS - Deferred Structural Elucidation Analysis for Molecular Simulations
3 //
4 // Copyright (c) 2018--present d-SEAMS core team
5 //
6 // This program is free software: you can redistribute it and/or modify
7 // it under the terms of the MIT License as published by
8 // the Open Source Initiative.
9 //
10 // A copy of the MIT License is included in the LICENSE file of this repository.
11 // You should have received a copy of the MIT License along with this program.
12 // If not, see <https://opensource.org/licenses/MIT>.
13 //-----------------------------------------------------------------------------------
14 
15 #include <generic.hpp>
16 #include <seams_input.hpp>
17 
26  inpFile = std::make_unique<std::ifstream>(filename);
28  std::string line; // Current line being read in
29  std::vector<std::string> tokens; // Vector containing word tokens
30  std::vector<int> id; // Vector for the IDs in the ring
31 
32  if (!(gen::file_exists(filename))) {
33  std::cout << "Fatal Error: The bond file does not exist or you gave the "
34  "wrong path.\n";
35  // Throw exception?
36  return bonds;
37  }
38 
39  // Format of the file:
40  // 481 Bonds
41  // 272 214 906 1361 388 1
42  // 388 1361 1042 1548 237 1
43  // 272 1536 1582 1701 1905 1
44 
45  if (inpFile->is_open()) {
46  // ----------------------------------------------------------
47  // At this point we know that the file is open
48  std::getline((*inpFile), line); // Read in bonds
49  // Run this until EOF or you reach the next timestep
50  while (std::getline((*inpFile), line)) {
51  // Read in lines and tokenize them into std::string words
52  tokens = gen::tokenizer(line);
53 
54  id.clear();
55 
56  for (int i = 0; i < tokens.size(); i++) {
57  id.push_back(std::stoi(tokens[i]));
58  } // end of for, assigning values to id
59 
60  bonds.push_back(id);
61 
62  } // end of while, till EOF
63  // ----------------------------------------------------------
64  } // End of if file open statement
65 
66  inpFile->close();
67 
68  return bonds;
69 }
75  xyzFile = std::make_unique<std::ifstream>(filename);
77  yCloud;
78  std::string line; // Current line being read in
79  std::vector<std::string> tokens; // Vector containing word tokens
80  std::vector<double> numbers; // Vector containing type double numbers
81  int nop = -1; // Number of atoms in targetFrame
82  int iatom = 1; // Current atom being filled into the PointCloud
83  molSys::Point<double> iPoint; // Current point being read in from the file
84  double xLo, xHi, yLo, yHi, zLo,
85  zHi; // Box lengths extrapolated from the least and greatest coordinates
86 
87  if (!(gen::file_exists(filename))) {
88  std::cout
89  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
90  // Throw exception?
91  throw std::runtime_error("Wrong filepath");
92  }
93 
94  // --------
95  // Before filling up the PointCloud, if the vectors are filled
96  // empty them
97  //*yCloud = molSys::clearPointCloud(yCloud);
98  // --------
99 
100  // Format of an XYZ file:
101  // 1970
102  // generated by VMD
103  // O 43.603500 16.926201 15.215700
104  // O 39.912601 14.775100 19.379200
105  if (xyzFile->is_open()) {
106 
107  // ----------------------------------------------------------
108  // At this point we know that the XYZ file is open
109 
110  // The first line contains the number of atoms
111  std::getline((*xyzFile), line);
112  nop = std::stoi(line);
113  // Skip the comment line
114  std::getline((*xyzFile), line);
115  // Reserve memory for atom coordinates using nop read in
116  yCloud.pts.reserve(nop);
117  yCloud.nop = nop;
118  // Run this until EOF or you reach the next timestep
119  while (std::getline((*xyzFile), line)) {
120 
121  // Read in lines and tokenize them into std::string words and <double>
122  // numbers
123  tokens = gen::tokenizer(line);
124  numbers = gen::tokenizerDouble(line);
125 
126  // Skip whitespace
127  if (tokens.size() == 0) {
128  continue;
129  }
130 
131  // Put logic for checking atom type here later
132  iPoint.type = 1; // Oxygen type; hard-coded!
133  iPoint.molID = iatom;
134  iPoint.atomID = iatom;
135  iPoint.x = std::stod(tokens[1]);
136  iPoint.y = std::stod(tokens[2]);
137  iPoint.z = std::stod(tokens[3]);
138 
139  yCloud.pts.push_back(iPoint);
140  yCloud.idIndexMap[iPoint.atomID] = yCloud.pts.size() - 1;
141  iatom++; // Increase index
142 
143  // First point
144  if (iatom == 1) {
145  // Loop through the dimensions
146  xLo = iPoint.x;
147  xHi = iPoint.x;
148  yLo = iPoint.y;
149  yHi = iPoint.y;
150  zLo = iPoint.z;
151  zHi = iPoint.z;
152  } // first point
153  else {
154  if (iPoint.x < xLo) {
155  xLo = iPoint.x;
156  } // xLo
157  else if (iPoint.x > xHi) {
158  xHi = iPoint.x;
159  } // xHi
160  else if (iPoint.y < yLo) {
161  yLo = iPoint.y;
162  } // yLo
163  else if (iPoint.y > yHi) {
164  yHi = iPoint.y;
165  } // yHi
166  else if (iPoint.z < zLo) {
167  zLo = iPoint.z;
168  } // zLo
169  else if (iPoint.z > zHi) {
170  zHi = iPoint.z;
171  } // zHi
172  } // update
173  } // end of while, looping through lines till EOF
174  // ----------------------------------------------------------
175  } // End of if file open statement
176 
177  xyzFile->close();
178 
179  if (yCloud.pts.size() == 1) {
180  xHi = xLo + 10;
181  yHi = yLo + 10;
182  zHi = zLo + 10;
183  } // for a single point in the system (never happens)
184 
185  // Fill up the box lengths
186  yCloud.box.push_back(xHi - xLo);
187  yCloud.box.push_back(yHi - yLo);
188  yCloud.box.push_back(zHi - zLo);
189 
190  if (yCloud.pts.size() != yCloud.nop) {
191  std::cout << "Atoms didn't get filled in properly.\n";
192  }
193 
194  return yCloud;
195 }
196 
197 // External Libraries
198 
212 sinp::readLammpsTrj(std::string filename, int targetFrame,
213  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
214  bool isSlice, std::array<double, 3> coordLow,
215  std::array<double, 3> coordHigh) {
217  dumpFile = std::make_unique<std::ifstream>(filename);
218  std::string line; // Current line being read in
219  std::vector<std::string> tokens; // Vector containing word tokens
220  std::vector<double> numbers; // Vector containing type double numbers
221  std::vector<double> tilt; // Vector containing tilt factors
222  int currentFrame = 0; // Current frame being read in
223  int nop = -1; // Number of atoms in targetFrame
224  bool foundFrame =
225  false; // Determines whether targetFrame has been found or not
226  bool readNOP = false; // Flag for reading in the number of atoms
227  bool readBox = false; // Flag for reading in the box lengths
228  bool readAtoms = false; // Flag for reading in the atoms
229  int xIndex, yIndex, zIndex,
230  typeIndex; // Indices for x,y,z coordinates, and LAMMPS type ID
231  int molIndex = 0; // Index for molecular ID
232  int atomIndex = 0; // Index for atom ID (Only used if mol ID has not been set)
233  molSys::Point<double> iPoint; // Current point being read in from the file
234  xIndex = yIndex = zIndex = typeIndex = -1; // Default values
235  bool isTriclinic = false; // Flag for an orthogonal or triclinic box
236 
237  if (!(gen::file_exists(filename))) {
238  std::cout
239  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
240  // Throw exception?
241  return *yCloud;
242  }
243 
244  // The format of the LAMMPS trajectory file is:
245  // ITEM: TIMESTEP
246  // 0
247  // ITEM: NUMBER OF ATOMS
248  // 4096
249  // ITEM: BOX BOUNDS pp pp pp
250  // -7.9599900000000001e-01 5.0164000000000001e+01
251  // -7.9599900000000001e-01 5.0164000000000001e+01
252  // -7.9599900000000001e-01 5.0164000000000001e+01
253  // ITEM: ATOMS id type x y z
254  // 1 1 0 0 0 etc
255  if (dumpFile->is_open()) {
256  // ----------------------------------------------------------
257  // At this point we know that the dumpfile is open
258  // This loop searches for targetFrame
259  while (std::getline((*dumpFile), line)) {
260  // Read in lines and tokenize them
261  tokens = gen::tokenizer(line);
262  // Find out which timestep number
263  // you are inside
264  if (tokens[0].compare("ITEM:") == 0) {
265  if (tokens[1].compare("TIMESTEP") == 0) {
266  // Now you are in a new timestep. Update frame number
267  currentFrame++;
268  }
269  }
270 
271  // If targetFrame has been found
272  // break out of the while loop
273  if (currentFrame == targetFrame) {
274  foundFrame = true;
275  break; // Exit the while loop
276  }
277  } // End of while loop searching for targetFrame
278  // ----------------------------------------------------------
279  // Before filling up the PointCloud, if the vectors are filled
280  // empty them
281  *yCloud = molSys::clearPointCloud(yCloud);
282 
283  // ----------------------------------------------------------
284  // If targetFrame has been found, read in the box lengths,
285  // number of atoms and then read in atom positions, type, molID
286  // By default, set molID=1 if not specified
287  if (foundFrame) {
288  // Run this until EOF or you reach the next timestep
289  while (std::getline((*dumpFile), line)) {
290  // Read in lines and tokenize them into std::string words and <double>
291  // numbers
292  tokens = gen::tokenizer(line);
293  numbers = gen::tokenizerDouble(line);
294 
295  // If you've reached the timestep line then you've reached the
296  // next frame. Break out of the while loop
297  if (tokens[0].compare("ITEM:") == 0) {
298  if (tokens[1].compare("TIMESTEP") == 0) {
299  break;
300  }
301  }
302 
303  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
304  // Read number of particles
305  if (readNOP) {
306  nop = std::stoi(line.data());
307  readNOP = false;
308  yCloud->pts.reserve(nop);
309  yCloud->nop = nop;
310  }
311  // Read box lengths
312  if (readBox) {
313  // You've reached the end of box lengths
314  if (tokens[0].compare("ITEM:") == 0) {
315  readBox = false;
316  // If the box is triclinic, get the
317  // orthogonal 'bounding box'
318  if (isTriclinic) {
319  // Update tilt factors
320  for (int k = 0; k < tilt.size(); k++) {
321  yCloud->box.push_back(tilt[k]);
322  }
323  } // end of check for triclinic
324  }
325  // Or else fill up the box lengths
326  else {
327  yCloud->box.push_back(numbers[1] - numbers[0]); // Update box length
328  yCloud->boxLow.push_back(
329  numbers[0]); // Update the lower box coordinate
330  // Do this for a triclinic box only
331  if (numbers.size() == 3) {
332  isTriclinic = true;
333  tilt.push_back(numbers[2]);
334  }
335  }
336  }
337  // Read atoms into yCloud line by line
338  if (readAtoms) {
339  iPoint.type = numbers[typeIndex];
340  iPoint.molID = numbers[molIndex];
341  iPoint.atomID = numbers[atomIndex];
342  iPoint.x = numbers[xIndex];
343  iPoint.y = numbers[yIndex];
344  iPoint.z = numbers[zIndex];
345  // Check if the particle is inside the volume Slice
346  // or not
347  if (isSlice) { // only if a slice has been requested
348  iPoint.inSlice = sinp::atomInSlice(iPoint.x, iPoint.y, iPoint.z,
349  coordLow, coordHigh);
350  }
351  yCloud->pts.push_back(iPoint);
352  yCloud->idIndexMap[iPoint.atomID] = yCloud->pts.size() - 1;
353  }
354  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
355 
356  // Tests for reading in nop, box lengths, and atoms
357  if (tokens[0].compare("ITEM:") == 0) {
358  if (tokens[1].compare("NUMBER") == 0) {
359  readNOP = true;
360  }
361  }
362  if (tokens[0].compare("ITEM:") == 0) {
363  if (tokens[1].compare("BOX") == 0) {
364  readBox = true;
365  }
366  }
367  if (tokens[0].compare("ITEM:") == 0) {
368  if (tokens[1].compare("ATOMS") == 0) {
369  readAtoms = true;
370  // Now find out which index is the coordinate index etc
371  for (int i = 2; i < tokens.size(); i++) {
372  if (tokens[i].compare("type") == 0) {
373  typeIndex = i - 2;
374  }
375  if (tokens[i].compare("x") == 0) {
376  xIndex = i - 2;
377  }
378  if (tokens[i].compare("y") == 0) {
379  yIndex = i - 2;
380  }
381  if (tokens[i].compare("z") == 0) {
382  zIndex = i - 2;
383  }
384  if (tokens[i].compare("mol") == 0) {
385  molIndex = i - 2;
386  }
387  if (tokens[i].compare("id") == 0) {
388  atomIndex = i - 2;
389  }
390  } // End of for loop over tokens
391  if (molIndex == 0) {
392  molIndex = atomIndex;
393  } // Set mol ID=atomID if not given
394  }
395  } // End of nested if loops for checking atom
396 
397  } // End of while
398  } // End of targetFrame found
399  // ----------------------------------------------------------
400  } // End of if file open statement
401 
402  // Check if you filled in the frame correctly
403  if (!(foundFrame)) {
404  std::cout << "You entered a frame that doesn't exist.\n";
405  } // Throw exception
406  if (foundFrame) {
407  if (yCloud->pts.size() != yCloud->nop) {
408  std::cout << "Atoms didn't get filled in properly.\n";
409  }
410  } // Throw exception
411  yCloud->currentFrame = targetFrame;
412 
413  dumpFile->close();
414  return *yCloud;
415 }
416 
431 sinp::readLammpsTrjO(std::string filename, int targetFrame,
432  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
433  int typeO, bool isSlice, std::array<double, 3> coordLow,
434  std::array<double, 3> coordHigh) {
436  dumpFile = std::make_unique<std::ifstream>(filename);
437  std::string line; // Current line being read in
438  std::vector<std::string> tokens; // Vector containing word tokens
439  std::vector<double> numbers; // Vector containing type double numbers
440  std::vector<double> tilt; // Vector containing tilt factors
441  int currentFrame = 0; // Current frame being read in
442  int nop = -1; // Number of atoms in targetFrame
443  bool foundFrame =
444  false; // Determines whether targetFrame has been found or not
445  bool readNOP = false; // Flag for reading in the number of atoms
446  bool readBox = false; // Flag for reading in the box lengths
447  bool readAtoms = false; // Flag for reading in the atoms
448  int xIndex, yIndex, zIndex,
449  typeIndex; // Indices for x,y,z coordinates, and LAMMPS type ID
450  int molIndex = 0; // Index for molecular ID
451  int atomIndex = 0; // Index for atom ID (Only used if mol ID has not been set)
452  molSys::Point<double> iPoint; // Current point being read in from the file
453  xIndex = yIndex = zIndex = typeIndex = -1; // Default values
454  bool isTriclinic = false; // Flag for an orthogonal or triclinic box
455  int nOxy = 0; // Number of oxygen atoms
456 
457  if (!(gen::file_exists(filename))) {
458  std::cout
459  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
460  // Throw exception?
461  return *yCloud;
462  }
463 
464  // The format of the LAMMPS trajectory file is:
465  // ITEM: TIMESTEP
466  // 0
467  // ITEM: NUMBER OF ATOMS
468  // 4096
469  // ITEM: BOX BOUNDS pp pp pp
470  // -7.9599900000000001e-01 5.0164000000000001e+01
471  // -7.9599900000000001e-01 5.0164000000000001e+01
472  // -7.9599900000000001e-01 5.0164000000000001e+01
473  // ITEM: ATOMS id type x y z
474  // 1 1 0 0 0 etc
475  if (dumpFile->is_open()) {
476  // ----------------------------------------------------------
477  // At this point we know that the dumpfile is open
478  // This loop searches for targetFrame
479  while (std::getline((*dumpFile), line)) {
480  // Read in lines and tokenize them
481  tokens = gen::tokenizer(line);
482  // Find out which timestep number
483  // you are inside
484  if (tokens[0].compare("ITEM:") == 0) {
485  if (tokens[1].compare("TIMESTEP") == 0) {
486  // Now you are in a new timestep. Update frame number
487  currentFrame++;
488  }
489  }
490 
491  // If targetFrame has been found
492  // break out of the while loop
493  if (currentFrame == targetFrame) {
494  foundFrame = true;
495  break; // Exit the while loop
496  }
497  } // End of while loop searching for targetFrame
498  // ----------------------------------------------------------
499  // Before filling up the PointCloud, if the vectors are filled
500  // empty them
501  *yCloud = molSys::clearPointCloud(yCloud);
502 
503  // ----------------------------------------------------------
504  // If targetFrame has been found, read in the box lengths,
505  // number of atoms and then read in atom positions, type, molID
506  // By default, set molID=1 if not specified
507  if (foundFrame) {
508  // Run this until EOF or you reach the next timestep
509  while (std::getline((*dumpFile), line)) {
510  // Read in lines and tokenize them into std::string words and <double>
511  // numbers
512  tokens = gen::tokenizer(line);
513  numbers = gen::tokenizerDouble(line);
514 
515  // If you've reached the timestep line then you've reached the
516  // next frame. Break out of the while loop
517  if (tokens[0].compare("ITEM:") == 0) {
518  if (tokens[1].compare("TIMESTEP") == 0) {
519  break;
520  }
521  }
522 
523  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
524  // Read number of particles
525  if (readNOP) {
526  nop = std::stoi(line.data());
527  readNOP = false;
528  }
529  // Read box lengths
530  if (readBox) {
531  // You've reached the end of box lengths
532  if (tokens[0].compare("ITEM:") == 0) {
533  readBox = false;
534  // If the box is triclinic, get the
535  // orthogonal 'bounding box'
536  if (isTriclinic) {
537  // Update tilt factors
538  for (int k = 0; k < tilt.size(); k++) {
539  yCloud->box.push_back(tilt[k]);
540  }
541  } // end of check for triclinic
542  }
543  // Or else fill up the box lengths
544  else {
545  yCloud->box.push_back(numbers[1] - numbers[0]); // Update box length
546  yCloud->boxLow.push_back(
547  numbers[0]); // Update the lower box coordinate
548  // Do this for a triclinic box only
549  if (numbers.size() == 3) {
550  isTriclinic = true;
551  tilt.push_back(numbers[2]);
552  }
553  }
554  }
555  // Read atoms into yCloud line by line
556  if (readAtoms) {
557  iPoint.type = numbers[typeIndex];
558  iPoint.molID = numbers[molIndex];
559  iPoint.atomID = numbers[atomIndex];
560  iPoint.x = numbers[xIndex];
561  iPoint.y = numbers[yIndex];
562  iPoint.z = numbers[zIndex];
563  // Check if the particle is inside the volume Slice
564  // or not
565  if (isSlice) { // only if a slice has been requested
566  iPoint.inSlice = sinp::atomInSlice(iPoint.x, iPoint.y, iPoint.z,
567  coordLow, coordHigh);
568  }
569  // Save only oxygen atoms
570  if (iPoint.type == typeO) {
571  nOxy++;
572  // yCloud->pts.resize(yCloud->pts.size()+1);
573  yCloud->pts.push_back(iPoint);
574  yCloud->idIndexMap[iPoint.atomID] = yCloud->pts.size() - 1;
575  }
576  }
577  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
578 
579  // Tests for reading in nop, box lengths, and atoms
580  if (tokens[0].compare("ITEM:") == 0) {
581  if (tokens[1].compare("NUMBER") == 0) {
582  readNOP = true;
583  }
584  }
585  if (tokens[0].compare("ITEM:") == 0) {
586  if (tokens[1].compare("BOX") == 0) {
587  readBox = true;
588  }
589  }
590  if (tokens[0].compare("ITEM:") == 0) {
591  if (tokens[1].compare("ATOMS") == 0) {
592  readAtoms = true;
593  // Now find out which index is the coordinate index etc
594  for (int i = 2; i < tokens.size(); i++) {
595  if (tokens[i].compare("type") == 0) {
596  typeIndex = i - 2;
597  }
598  if (tokens[i].compare("x") == 0) {
599  xIndex = i - 2;
600  }
601  if (tokens[i].compare("y") == 0) {
602  yIndex = i - 2;
603  }
604  if (tokens[i].compare("z") == 0) {
605  zIndex = i - 2;
606  }
607  if (tokens[i].compare("mol") == 0) {
608  molIndex = i - 2;
609  }
610  if (tokens[i].compare("id") == 0) {
611  atomIndex = i - 2;
612  }
613  } // End of for loop over tokens
614  if (molIndex == 0) {
615  molIndex = atomIndex;
616  } // Set mol ID=atomID if not given
617  }
618  } // End of nested if loops for checking atom
619 
620  } // End of while
621  } // End of targetFrame found
622  // ----------------------------------------------------------
623  } // End of if file open statement
624 
625  // Check if you filled in the frame correctly
626  if (!(foundFrame)) {
627  std::cout << "You entered a frame that doesn't exist.\n";
628  } // Throw exception
629  if (foundFrame) {
630  yCloud->nop = yCloud->pts.size();
631  if (yCloud->pts.size() != nOxy) {
632  std::cout << "Atoms didn't get filled in properly.\n";
633  }
634  } // Throw exception
635  yCloud->currentFrame = targetFrame;
636 
637  dumpFile->close();
638  return *yCloud;
639 }
640 
657  std::string filename, int targetFrame,
658  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int typeI,
659  bool isSlice, std::array<double, 3> coordLow,
660  std::array<double, 3> coordHigh) {
662  dumpFile = std::make_unique<std::ifstream>(filename);
663  std::string line; // Current line being read in
664  std::vector<std::string> tokens; // Vector containing word tokens
665  std::vector<double> numbers; // Vector containing type double numbers
666  std::vector<double> tilt; // Vector containing tilt factors
667  int currentFrame = 0; // Current frame being read in
668  int nop = -1; // Number of atoms in targetFrame
669  bool foundFrame =
670  false; // Determines whether targetFrame has been found or not
671  bool readNOP = false; // Flag for reading in the number of atoms
672  bool readBox = false; // Flag for reading in the box lengths
673  bool readAtoms = false; // Flag for reading in the atoms
674  int xIndex, yIndex, zIndex,
675  typeIndex; // Indices for x,y,z coordinates, and LAMMPS type ID
676  int molIndex = 0; // Index for molecular ID
677  int atomIndex = 0; // Index for atom ID (Only used if mol ID has not been set)
678  molSys::Point<double> iPoint; // Current point being read in from the file
679  xIndex = yIndex = zIndex = typeIndex = -1; // Default values
680  bool isTriclinic = false; // Flag for an orthogonal or triclinic box
681  int nOxy = 0; // Number of oxygen atoms
682 
683  if (!(gen::file_exists(filename))) {
684  std::cout
685  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
686  // Throw exception?
687  return *yCloud;
688  }
689 
690  // The format of the LAMMPS trajectory file is:
691  // ITEM: TIMESTEP
692  // 0
693  // ITEM: NUMBER OF ATOMS
694  // 4096
695  // ITEM: BOX BOUNDS pp pp pp
696  // -7.9599900000000001e-01 5.0164000000000001e+01
697  // -7.9599900000000001e-01 5.0164000000000001e+01
698  // -7.9599900000000001e-01 5.0164000000000001e+01
699  // ITEM: ATOMS id type x y z
700  // 1 1 0 0 0 etc
701  if (dumpFile->is_open()) {
702  // ----------------------------------------------------------
703  // At this point we know that the dumpfile is open
704  // This loop searches for targetFrame
705  while (std::getline((*dumpFile), line)) {
706  // Read in lines and tokenize them
707  tokens = gen::tokenizer(line);
708  // Find out which timestep number
709  // you are inside
710  if (tokens[0].compare("ITEM:") == 0) {
711  if (tokens[1].compare("TIMESTEP") == 0) {
712  // Now you are in a new timestep. Update frame number
713  currentFrame++;
714  }
715  }
716 
717  // If targetFrame has been found
718  // break out of the while loop
719  if (currentFrame == targetFrame) {
720  foundFrame = true;
721  break; // Exit the while loop
722  }
723  } // End of while loop searching for targetFrame
724  // ----------------------------------------------------------
725  // Before filling up the PointCloud, if the vectors are filled
726  // empty them
727  *yCloud = molSys::clearPointCloud(yCloud);
728 
729  // ----------------------------------------------------------
730  // If targetFrame has been found, read in the box lengths,
731  // number of atoms and then read in atom positions, type, molID
732  // By default, set molID=1 if not specified
733  if (foundFrame) {
734  // Run this until EOF or you reach the next timestep
735  while (std::getline((*dumpFile), line)) {
736  // Read in lines and tokenize them into std::string words and <double>
737  // numbers
738  tokens = gen::tokenizer(line);
739  numbers = gen::tokenizerDouble(line);
740 
741  // If you've reached the timestep line then you've reached the
742  // next frame. Break out of the while loop
743  if (tokens[0].compare("ITEM:") == 0) {
744  if (tokens[1].compare("TIMESTEP") == 0) {
745  break;
746  }
747  }
748 
749  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
750  // Read number of particles
751  if (readNOP) {
752  nop = std::stoi(line.data());
753  readNOP = false;
754  }
755  // Read box lengths
756  if (readBox) {
757  // You've reached the end of box lengths
758  if (tokens[0].compare("ITEM:") == 0) {
759  readBox = false;
760  // If the box is triclinic, get the
761  // orthogonal 'bounding box'
762  if (isTriclinic) {
763  // Update tilt factors
764  for (int k = 0; k < tilt.size(); k++) {
765  yCloud->box.push_back(tilt[k]);
766  }
767  } // end of check for triclinic
768  }
769  // Or else fill up the box lengths
770  else {
771  yCloud->box.push_back(numbers[1] - numbers[0]); // Update box length
772  yCloud->boxLow.push_back(
773  numbers[0]); // Update the lower box coordinate
774  // Do this for a triclinic box only
775  if (numbers.size() == 3) {
776  isTriclinic = true;
777  tilt.push_back(numbers[2]);
778  }
779  }
780  }
781  // Read atoms into yCloud line by line
782  if (readAtoms) {
783  iPoint.type = numbers[typeIndex];
784  iPoint.molID = numbers[molIndex];
785  iPoint.atomID = numbers[atomIndex];
786  iPoint.x = numbers[xIndex];
787  iPoint.y = numbers[yIndex];
788  iPoint.z = numbers[zIndex];
789  // Check if the particle is inside the volume Slice
790  // or not
791  if (isSlice) { // only if a slice has been requested
792  iPoint.inSlice = sinp::atomInSlice(iPoint.x, iPoint.y, iPoint.z,
793  coordLow, coordHigh);
794  // Skip if the atom is not part of the slice
795  if (!iPoint.inSlice) {
796  continue;
797  } // do not save if not inside the slice
798  }
799  // Save only atoms of the desired type
800  if (iPoint.type == typeI) {
801  nOxy++;
802  // yCloud->pts.resize(yCloud->pts.size()+1);
803  yCloud->pts.push_back(iPoint);
804  yCloud->idIndexMap[iPoint.atomID] =
805  yCloud->pts.size() - 1; // array index
806  }
807  }
808  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
809 
810  // Tests for reading in nop, box lengths, and atoms
811  if (tokens[0].compare("ITEM:") == 0) {
812  if (tokens[1].compare("NUMBER") == 0) {
813  readNOP = true;
814  }
815  }
816  if (tokens[0].compare("ITEM:") == 0) {
817  if (tokens[1].compare("BOX") == 0) {
818  readBox = true;
819  }
820  }
821  if (tokens[0].compare("ITEM:") == 0) {
822  if (tokens[1].compare("ATOMS") == 0) {
823  readAtoms = true;
824  // Now find out which index is the coordinate index etc
825  for (int i = 2; i < tokens.size(); i++) {
826  if (tokens[i].compare("type") == 0) {
827  typeIndex = i - 2;
828  }
829  if (tokens[i].compare("x") == 0) {
830  xIndex = i - 2;
831  }
832  if (tokens[i].compare("y") == 0) {
833  yIndex = i - 2;
834  }
835  if (tokens[i].compare("z") == 0) {
836  zIndex = i - 2;
837  }
838  if (tokens[i].compare("mol") == 0) {
839  molIndex = i - 2;
840  }
841  if (tokens[i].compare("id") == 0) {
842  atomIndex = i - 2;
843  }
844  } // End of for loop over tokens
845  if (molIndex == 0) {
846  molIndex = atomIndex;
847  } // Set mol ID=atomID if not given
848  }
849  } // End of nested if loops for checking atom
850 
851  } // End of while
852  } // End of targetFrame found
853  // ----------------------------------------------------------
854  } // End of if file open statement
855 
856  // Check if you filled in the frame correctly
857  if (!(foundFrame)) {
858  std::cout << "You entered a frame that doesn't exist.\n";
859  } // Throw exception
860 
861  // Update the number of particles
862  yCloud->nop = yCloud->pts.size();
863 
864  // Update the frame number
865  yCloud->currentFrame = targetFrame;
866 
867  dumpFile->close();
868  return *yCloud;
869 }
T data(T... args)
File for containing generic or common functions.
T getline(T... args)
bool file_exists(const std::string &name)
Function for checking if a file exists or not.
Definition: generic.hpp:298
std::vector< std::string > tokenizer(std::string line)
Function for tokenizing line strings into words (strings) delimited by whitespace....
Definition: generic.hpp:259
std::vector< double > tokenizerDouble(std::string line)
Function for tokenizing line strings into a vector of doubles.
Definition: generic.hpp:270
T x
type ID, molID, atomID
Definition: mol_sys.hpp:151
int nop
Current frame number.
Definition: mol_sys.hpp:173
std::vector< S > pts
Definition: mol_sys.hpp:171
bool inSlice
Type of ice/water etc based on cij.
Definition: mol_sys.hpp:155
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:176
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:174
molSys::PointCloud< molSys::Point< double >, double > clearPointCloud(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
//! Function for clearing vectors in PointCloud after multiple usage
Definition: mol_sys.cpp:24
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrj(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrjreduced(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrjO(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeO, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
molSys::PointCloud< molSys::Point< double >, double > readXYZ(std::string filename)
Function for reading in atom coordinates from an XYZ file.
Definition: seams_input.cpp:73
bool atomInSlice(double x, double y, double z, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh)
Definition: seams_input.hpp:88
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
T push_back(T... args)
T reserve(T... args)
File for functions that read in files).
T size(T... args)
T stod(T... args)
T stoi(T... args)
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:170
This contains per-particle information.
Definition: mol_sys.hpp:149