Sinp

Functions

std::vector< std::stringsinp::getInpFileList (std::string inputFolder)
 Get file list inside the input folder. More...
 
molSys::PointCloud< molSys::Point< double >, double > sinp::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 > sinp::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 > sinp::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})
 
int sinp::readXYZ (std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 Function for reading in atom coordinates from an XYZ file. More...
 
std::vector< std::vector< int > > sinp::readBonds (std::string filename)
 Reads bonds into a vector of vectors from a file with a specific format. More...
 
bool sinp::atomInSlice (double x, double y, double z, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh)
 

Detailed Description

Function Documentation

◆ atomInSlice()

bool sinp::atomInSlice ( double  x,
double  y,
double  z,
std::array< double, 3 >  coordLow,
std::array< double, 3 >  coordHigh 
)
inline

If this is 3 then the particle is inside the volume slice

Definition at line 85 of file seams_input.hpp.

87  {
88  int flag = 0;
89 
90  if (((x >= coordLow[0]) && (x <= coordHigh[0])) ||
91  coordLow[0] == coordHigh[0]) {
92  flag++;
93  }
94  if (((y >= coordLow[1]) && (y <= coordHigh[1])) ||
95  coordLow[1] == coordHigh[1]) {
96  flag++;
97  }
98  if (((z >= coordLow[2]) && (z <= coordHigh[2])) ||
99  coordLow[2] == coordHigh[2]) {
100  flag++;
101  }
102 
103  if (flag == 3) {
104  return true;
105  } else {
106  return false;
107  }
108 }

◆ getInpFileList()

std::vector<std::string> sinp::getInpFileList ( std::string  inputFolder)

Get file list inside the input folder.

◆ readBonds()

std::vector< std::vector< int > > sinp::readBonds ( std::string  filename)

Reads bonds into a vector of vectors from a file with a specific format.

Get all the ring information, from the R.I.N.G.S. file. Each line contains the IDs of the atoms in the ring. This is saved inside a vector of vectors. Rings which have more than three consecutive water molecules are discarded.

Definition at line 20 of file seams_input.cpp.

20  {
22  inpFile = std::make_unique<std::ifstream>(filename);
24  std::string line; // Current line being read in
25  std::vector<std::string> tokens; // Vector containing word tokens
26  std::vector<int> id; // Vector for the IDs in the ring
27 
28  if (!(gen::file_exists(filename))) {
29  std::cout << "Fatal Error: The bond file does not exist or you gave the "
30  "wrong path.\n";
31  // Throw exception?
32  return bonds;
33  }
34 
35  // Format of the file:
36  // 481 Bonds
37  // 272 214 906 1361 388 1
38  // 388 1361 1042 1548 237 1
39  // 272 1536 1582 1701 1905 1
40 
41  if (inpFile->is_open()) {
42  // ----------------------------------------------------------
43  // At this point we know that the file is open
44  std::getline((*inpFile), line); // Read in bonds
45  // Run this until EOF or you reach the next timestep
46  while (std::getline((*inpFile), line)) {
47  // Read in lines and tokenize them into std::string words
48  tokens = gen::tokenizer(line);
49 
50  id.clear();
51 
52  for (int i = 0; i < tokens.size(); i++) {
53  id.push_back(std::stoi(tokens[i]));
54  } // end of for, assigning values to id
55 
56  bonds.push_back(id);
57 
58  } // end of while, till EOF
59  // ----------------------------------------------------------
60  } // End of if file open statement
61 
62  inpFile->close();
63 
64  return bonds;
65 }
T getline(T... args)
bool file_exists(const std::string &name)
Function for checking if a file exists or not.
Definition: generic.hpp:294
std::vector< std::string > tokenizer(std::string line)
Function for tokenizing line strings into words (strings) delimited by whitespace....
Definition: generic.hpp:255
T push_back(T... args)
T size(T... args)
T stoi(T... args)

◆ readLammpsTrj()

molSys::PointCloud< molSys::Point< double >, double > sinp::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} 
)

Function for reading in a specified frame (frame number and not timestep value)

Function for reading in a lammps file. Reads in a specified frame (frame number and not timestep value).

Parameters
[in]filenameThe name of the lammps trajectory file to be read in
[in]targetFrameThe frame number whose information will be read in
[out]yCloudThe outputted PointCloud
[in]isSliceThis decides whether a slice will be created or not
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 207 of file seams_input.cpp.

210  {
212  dumpFile = std::make_unique<std::ifstream>(filename);
213  std::string line; // Current line being read in
214  std::vector<std::string> tokens; // Vector containing word tokens
215  std::vector<double> numbers; // Vector containing type double numbers
216  std::vector<double> tilt; // Vector containing tilt factors
217  int currentFrame = 0; // Current frame being read in
218  int nop = -1; // Number of atoms in targetFrame
219  bool foundFrame =
220  false; // Determines whether targetFrame has been found or not
221  bool readNOP = false; // Flag for reading in the number of atoms
222  bool readBox = false; // Flag for reading in the box lengths
223  bool readAtoms = false; // Flag for reading in the atoms
224  int xIndex, yIndex, zIndex,
225  typeIndex; // Indices for x,y,z coordinates, and LAMMPS type ID
226  int molIndex = 0; // Index for molecular ID
227  int atomIndex = 0; // Index for atom ID (Only used if mol ID has not been set)
228  molSys::Point<double> iPoint; // Current point being read in from the file
229  xIndex = yIndex = zIndex = typeIndex = -1; // Default values
230  bool isTriclinic = false; // Flag for an orthogonal or triclinic box
231 
232  if (!(gen::file_exists(filename))) {
233  std::cout
234  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
235  // Throw exception?
236  return *yCloud;
237  }
238 
239  // The format of the LAMMPS trajectory file is:
240  // ITEM: TIMESTEP
241  // 0
242  // ITEM: NUMBER OF ATOMS
243  // 4096
244  // ITEM: BOX BOUNDS pp pp pp
245  // -7.9599900000000001e-01 5.0164000000000001e+01
246  // -7.9599900000000001e-01 5.0164000000000001e+01
247  // -7.9599900000000001e-01 5.0164000000000001e+01
248  // ITEM: ATOMS id type x y z
249  // 1 1 0 0 0 etc
250  if (dumpFile->is_open()) {
251  // ----------------------------------------------------------
252  // At this point we know that the dumpfile is open
253  // This loop searches for targetFrame
254  while (std::getline((*dumpFile), line)) {
255  // Read in lines and tokenize them
256  tokens = gen::tokenizer(line);
257  // Find out which timestep number
258  // you are inside
259  if (tokens[0].compare("ITEM:") == 0) {
260  if (tokens[1].compare("TIMESTEP") == 0) {
261  // Now you are in a new timestep. Update frame number
262  currentFrame++;
263  }
264  }
265 
266  // If targetFrame has been found
267  // break out of the while loop
268  if (currentFrame == targetFrame) {
269  foundFrame = true;
270  break; // Exit the while loop
271  }
272  } // End of while loop searching for targetFrame
273  // ----------------------------------------------------------
274  // Before filling up the PointCloud, if the vectors are filled
275  // empty them
276  *yCloud = molSys::clearPointCloud(yCloud);
277 
278  // ----------------------------------------------------------
279  // If targetFrame has been found, read in the box lengths,
280  // number of atoms and then read in atom positions, type, molID
281  // By default, set molID=1 if not specified
282  if (foundFrame) {
283  // Run this until EOF or you reach the next timestep
284  while (std::getline((*dumpFile), line)) {
285  // Read in lines and tokenize them into std::string words and <double>
286  // numbers
287  tokens = gen::tokenizer(line);
288  numbers = gen::tokenizerDouble(line);
289 
290  // If you've reached the timestep line then you've reached the
291  // next frame. Break out of the while loop
292  if (tokens[0].compare("ITEM:") == 0) {
293  if (tokens[1].compare("TIMESTEP") == 0) {
294  break;
295  }
296  }
297 
298  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
299  // Read number of particles
300  if (readNOP) {
301  nop = std::stoi(line.data());
302  readNOP = false;
303  yCloud->pts.reserve(nop);
304  yCloud->nop = nop;
305  }
306  // Read box lengths
307  if (readBox) {
308  // You've reached the end of box lengths
309  if (tokens[0].compare("ITEM:") == 0) {
310  readBox = false;
311  // If the box is triclinic, get the
312  // orthogonal 'bounding box'
313  if (isTriclinic) {
314  // Update tilt factors
315  for (int k = 0; k < tilt.size(); k++) {
316  yCloud->box.push_back(tilt[k]);
317  }
318  } // end of check for triclinic
319  }
320  // Or else fill up the box lengths
321  else {
322  yCloud->box.push_back(numbers[1] - numbers[0]); // Update box length
323  yCloud->boxLow.push_back(
324  numbers[0]); // Update the lower box coordinate
325  // Do this for a triclinic box only
326  if (numbers.size() == 3) {
327  isTriclinic = true;
328  tilt.push_back(numbers[2]);
329  }
330  }
331  }
332  // Read atoms into yCloud line by line
333  if (readAtoms) {
334  iPoint.type = numbers[typeIndex];
335  iPoint.molID = numbers[molIndex];
336  iPoint.atomID = numbers[atomIndex];
337  iPoint.x = numbers[xIndex];
338  iPoint.y = numbers[yIndex];
339  iPoint.z = numbers[zIndex];
340  // Check if the particle is inside the volume Slice
341  // or not
342  if (isSlice) { // only if a slice has been requested
343  iPoint.inSlice = sinp::atomInSlice(iPoint.x, iPoint.y, iPoint.z,
344  coordLow, coordHigh);
345  }
346  yCloud->pts.push_back(iPoint);
347  yCloud->idIndexMap[iPoint.atomID] = yCloud->pts.size() - 1;
348  }
349  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
350 
351  // Tests for reading in nop, box lengths, and atoms
352  if (tokens[0].compare("ITEM:") == 0) {
353  if (tokens[1].compare("NUMBER") == 0) {
354  readNOP = true;
355  }
356  }
357  if (tokens[0].compare("ITEM:") == 0) {
358  if (tokens[1].compare("BOX") == 0) {
359  readBox = true;
360  }
361  }
362  if (tokens[0].compare("ITEM:") == 0) {
363  if (tokens[1].compare("ATOMS") == 0) {
364  readAtoms = true;
365  // Now find out which index is the coordinate index etc
366  for (int i = 2; i < tokens.size(); i++) {
367  if (tokens[i].compare("type") == 0) {
368  typeIndex = i - 2;
369  }
370  if (tokens[i].compare("x") == 0) {
371  xIndex = i - 2;
372  }
373  if (tokens[i].compare("y") == 0) {
374  yIndex = i - 2;
375  }
376  if (tokens[i].compare("z") == 0) {
377  zIndex = i - 2;
378  }
379  if (tokens[i].compare("mol") == 0) {
380  molIndex = i - 2;
381  }
382  if (tokens[i].compare("id") == 0) {
383  atomIndex = i - 2;
384  }
385  } // End of for loop over tokens
386  if (molIndex == 0) {
387  molIndex = atomIndex;
388  } // Set mol ID=atomID if not given
389  }
390  } // End of nested if loops for checking atom
391 
392  } // End of while
393  } // End of targetFrame found
394  // ----------------------------------------------------------
395  } // End of if file open statement
396 
397  // Check if you filled in the frame correctly
398  if (!(foundFrame)) {
399  std::cout << "You entered a frame that doesn't exist.\n";
400  } // Throw exception
401  if (foundFrame) {
402  if (yCloud->pts.size() != yCloud->nop) {
403  std::cout << "Atoms didn't get filled in properly.\n";
404  }
405  } // Throw exception
406  yCloud->currentFrame = targetFrame;
407 
408  dumpFile->close();
409  return *yCloud;
410 }
T data(T... args)
std::vector< double > tokenizerDouble(std::string line)
Function for tokenizing line strings into a vector of doubles.
Definition: generic.hpp:266
T x
type ID, molID, atomID
Definition: mol_sys.hpp:147
int nop
Current frame number.
Definition: mol_sys.hpp:169
std::vector< S > pts
Definition: mol_sys.hpp:167
bool inSlice
Type of ice/water etc based on cij.
Definition: mol_sys.hpp:151
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:172
std::vector< T > box
Number of atoms.
Definition: mol_sys.hpp:170
std::vector< T > boxLow
Periodic box lengths.
Definition: mol_sys.hpp:171
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:20
int currentFrame
Collection of points.
Definition: mol_sys.hpp:168
bool atomInSlice(double x, double y, double z, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh)
Definition: seams_input.hpp:85
T reserve(T... args)
This contains per-particle information.
Definition: mol_sys.hpp:145

◆ readLammpsTrjO()

molSys::PointCloud< molSys::Point< double >, double > sinp::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} 
)

Function for reading in a specified frame (frame number and not timestep value) / This only reads in oxygen atoms

Function for reading in a lammps file; and saves only the Oxygen atoms. This is an overloaded function. The Oxygen atom ID must be specified.

Parameters
[in]filenameThe name of the lammps trajectory file to be read in
[in]targetFrameThe frame number whose information will be read in
[out]yCloudThe outputted PointCloud
[in]typeOThe type ID of the Oxygen atoms
[in]isSliceThis decides whether a slice will be created or not
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 426 of file seams_input.cpp.

429  {
431  dumpFile = std::make_unique<std::ifstream>(filename);
432  std::string line; // Current line being read in
433  std::vector<std::string> tokens; // Vector containing word tokens
434  std::vector<double> numbers; // Vector containing type double numbers
435  std::vector<double> tilt; // Vector containing tilt factors
436  int currentFrame = 0; // Current frame being read in
437  int nop = -1; // Number of atoms in targetFrame
438  bool foundFrame =
439  false; // Determines whether targetFrame has been found or not
440  bool readNOP = false; // Flag for reading in the number of atoms
441  bool readBox = false; // Flag for reading in the box lengths
442  bool readAtoms = false; // Flag for reading in the atoms
443  int xIndex, yIndex, zIndex,
444  typeIndex; // Indices for x,y,z coordinates, and LAMMPS type ID
445  int molIndex = 0; // Index for molecular ID
446  int atomIndex = 0; // Index for atom ID (Only used if mol ID has not been set)
447  molSys::Point<double> iPoint; // Current point being read in from the file
448  xIndex = yIndex = zIndex = typeIndex = -1; // Default values
449  bool isTriclinic = false; // Flag for an orthogonal or triclinic box
450  int nOxy = 0; // Number of oxygen atoms
451 
452  if (!(gen::file_exists(filename))) {
453  std::cout
454  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
455  // Throw exception?
456  return *yCloud;
457  }
458 
459  // The format of the LAMMPS trajectory file is:
460  // ITEM: TIMESTEP
461  // 0
462  // ITEM: NUMBER OF ATOMS
463  // 4096
464  // ITEM: BOX BOUNDS pp pp pp
465  // -7.9599900000000001e-01 5.0164000000000001e+01
466  // -7.9599900000000001e-01 5.0164000000000001e+01
467  // -7.9599900000000001e-01 5.0164000000000001e+01
468  // ITEM: ATOMS id type x y z
469  // 1 1 0 0 0 etc
470  if (dumpFile->is_open()) {
471  // ----------------------------------------------------------
472  // At this point we know that the dumpfile is open
473  // This loop searches for targetFrame
474  while (std::getline((*dumpFile), line)) {
475  // Read in lines and tokenize them
476  tokens = gen::tokenizer(line);
477  // Find out which timestep number
478  // you are inside
479  if (tokens[0].compare("ITEM:") == 0) {
480  if (tokens[1].compare("TIMESTEP") == 0) {
481  // Now you are in a new timestep. Update frame number
482  currentFrame++;
483  }
484  }
485 
486  // If targetFrame has been found
487  // break out of the while loop
488  if (currentFrame == targetFrame) {
489  foundFrame = true;
490  break; // Exit the while loop
491  }
492  } // End of while loop searching for targetFrame
493  // ----------------------------------------------------------
494  // Before filling up the PointCloud, if the vectors are filled
495  // empty them
496  *yCloud = molSys::clearPointCloud(yCloud);
497 
498  // ----------------------------------------------------------
499  // If targetFrame has been found, read in the box lengths,
500  // number of atoms and then read in atom positions, type, molID
501  // By default, set molID=1 if not specified
502  if (foundFrame) {
503  // Run this until EOF or you reach the next timestep
504  while (std::getline((*dumpFile), line)) {
505  // Read in lines and tokenize them into std::string words and <double>
506  // numbers
507  tokens = gen::tokenizer(line);
508  numbers = gen::tokenizerDouble(line);
509 
510  // If you've reached the timestep line then you've reached the
511  // next frame. Break out of the while loop
512  if (tokens[0].compare("ITEM:") == 0) {
513  if (tokens[1].compare("TIMESTEP") == 0) {
514  break;
515  }
516  }
517 
518  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
519  // Read number of particles
520  if (readNOP) {
521  nop = std::stoi(line.data());
522  readNOP = false;
523  }
524  // Read box lengths
525  if (readBox) {
526  // You've reached the end of box lengths
527  if (tokens[0].compare("ITEM:") == 0) {
528  readBox = false;
529  // If the box is triclinic, get the
530  // orthogonal 'bounding box'
531  if (isTriclinic) {
532  // Update tilt factors
533  for (int k = 0; k < tilt.size(); k++) {
534  yCloud->box.push_back(tilt[k]);
535  }
536  } // end of check for triclinic
537  }
538  // Or else fill up the box lengths
539  else {
540  yCloud->box.push_back(numbers[1] - numbers[0]); // Update box length
541  yCloud->boxLow.push_back(
542  numbers[0]); // Update the lower box coordinate
543  // Do this for a triclinic box only
544  if (numbers.size() == 3) {
545  isTriclinic = true;
546  tilt.push_back(numbers[2]);
547  }
548  }
549  }
550  // Read atoms into yCloud line by line
551  if (readAtoms) {
552  iPoint.type = numbers[typeIndex];
553  iPoint.molID = numbers[molIndex];
554  iPoint.atomID = numbers[atomIndex];
555  iPoint.x = numbers[xIndex];
556  iPoint.y = numbers[yIndex];
557  iPoint.z = numbers[zIndex];
558  // Check if the particle is inside the volume Slice
559  // or not
560  if (isSlice) { // only if a slice has been requested
561  iPoint.inSlice = sinp::atomInSlice(iPoint.x, iPoint.y, iPoint.z,
562  coordLow, coordHigh);
563  }
564  // Save only oxygen atoms
565  if (iPoint.type == typeO) {
566  nOxy++;
567  // yCloud->pts.resize(yCloud->pts.size()+1);
568  yCloud->pts.push_back(iPoint);
569  yCloud->idIndexMap[iPoint.atomID] = yCloud->pts.size() - 1;
570  }
571  }
572  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
573 
574  // Tests for reading in nop, box lengths, and atoms
575  if (tokens[0].compare("ITEM:") == 0) {
576  if (tokens[1].compare("NUMBER") == 0) {
577  readNOP = true;
578  }
579  }
580  if (tokens[0].compare("ITEM:") == 0) {
581  if (tokens[1].compare("BOX") == 0) {
582  readBox = true;
583  }
584  }
585  if (tokens[0].compare("ITEM:") == 0) {
586  if (tokens[1].compare("ATOMS") == 0) {
587  readAtoms = true;
588  // Now find out which index is the coordinate index etc
589  for (int i = 2; i < tokens.size(); i++) {
590  if (tokens[i].compare("type") == 0) {
591  typeIndex = i - 2;
592  }
593  if (tokens[i].compare("x") == 0) {
594  xIndex = i - 2;
595  }
596  if (tokens[i].compare("y") == 0) {
597  yIndex = i - 2;
598  }
599  if (tokens[i].compare("z") == 0) {
600  zIndex = i - 2;
601  }
602  if (tokens[i].compare("mol") == 0) {
603  molIndex = i - 2;
604  }
605  if (tokens[i].compare("id") == 0) {
606  atomIndex = i - 2;
607  }
608  } // End of for loop over tokens
609  if (molIndex == 0) {
610  molIndex = atomIndex;
611  } // Set mol ID=atomID if not given
612  }
613  } // End of nested if loops for checking atom
614 
615  } // End of while
616  } // End of targetFrame found
617  // ----------------------------------------------------------
618  } // End of if file open statement
619 
620  // Check if you filled in the frame correctly
621  if (!(foundFrame)) {
622  std::cout << "You entered a frame that doesn't exist.\n";
623  } // Throw exception
624  if (foundFrame) {
625  yCloud->nop = yCloud->pts.size();
626  if (yCloud->pts.size() != nOxy) {
627  std::cout << "Atoms didn't get filled in properly.\n";
628  }
629  } // Throw exception
630  yCloud->currentFrame = targetFrame;
631 
632  dumpFile->close();
633  return *yCloud;
634 }

◆ readLammpsTrjreduced()

molSys::PointCloud< molSys::Point< double >, double > sinp::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} 
)

Function that reads in only atoms pf the desired type and ignores all atoms which are not in the slice as well

Function for reading in a lammps file; and saves only the atoms of the desired type. Atoms which are not inside the slice or not of type I are not saved at all This is an overloaded function. The type atom ID must be specified.

Parameters
[in]filenameThe name of the lammps trajectory file to be read in
[in]targetFrameThe frame number whose information will be read in
[out]yCloudThe outputted PointCloud
[in]typeIThe type ID of the desired type of atoms
[in]isSliceThis decides whether a slice will be created or not
[in]coordLowContains the lower limits of the slice, if a slice is to be created
[in]coordHighContains the upper limits of the slice, if a slice is to be created

Definition at line 651 of file seams_input.cpp.

655  {
657  dumpFile = std::make_unique<std::ifstream>(filename);
658  std::string line; // Current line being read in
659  std::vector<std::string> tokens; // Vector containing word tokens
660  std::vector<double> numbers; // Vector containing type double numbers
661  std::vector<double> tilt; // Vector containing tilt factors
662  int currentFrame = 0; // Current frame being read in
663  int nop = -1; // Number of atoms in targetFrame
664  bool foundFrame =
665  false; // Determines whether targetFrame has been found or not
666  bool readNOP = false; // Flag for reading in the number of atoms
667  bool readBox = false; // Flag for reading in the box lengths
668  bool readAtoms = false; // Flag for reading in the atoms
669  int xIndex, yIndex, zIndex,
670  typeIndex; // Indices for x,y,z coordinates, and LAMMPS type ID
671  int molIndex = 0; // Index for molecular ID
672  int atomIndex = 0; // Index for atom ID (Only used if mol ID has not been set)
673  molSys::Point<double> iPoint; // Current point being read in from the file
674  xIndex = yIndex = zIndex = typeIndex = -1; // Default values
675  bool isTriclinic = false; // Flag for an orthogonal or triclinic box
676  int nOxy = 0; // Number of oxygen atoms
677 
678  if (!(gen::file_exists(filename))) {
679  std::cout
680  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
681  // Throw exception?
682  return *yCloud;
683  }
684 
685  // The format of the LAMMPS trajectory file is:
686  // ITEM: TIMESTEP
687  // 0
688  // ITEM: NUMBER OF ATOMS
689  // 4096
690  // ITEM: BOX BOUNDS pp pp pp
691  // -7.9599900000000001e-01 5.0164000000000001e+01
692  // -7.9599900000000001e-01 5.0164000000000001e+01
693  // -7.9599900000000001e-01 5.0164000000000001e+01
694  // ITEM: ATOMS id type x y z
695  // 1 1 0 0 0 etc
696  if (dumpFile->is_open()) {
697  // ----------------------------------------------------------
698  // At this point we know that the dumpfile is open
699  // This loop searches for targetFrame
700  while (std::getline((*dumpFile), line)) {
701  // Read in lines and tokenize them
702  tokens = gen::tokenizer(line);
703  // Find out which timestep number
704  // you are inside
705  if (tokens[0].compare("ITEM:") == 0) {
706  if (tokens[1].compare("TIMESTEP") == 0) {
707  // Now you are in a new timestep. Update frame number
708  currentFrame++;
709  }
710  }
711 
712  // If targetFrame has been found
713  // break out of the while loop
714  if (currentFrame == targetFrame) {
715  foundFrame = true;
716  break; // Exit the while loop
717  }
718  } // End of while loop searching for targetFrame
719  // ----------------------------------------------------------
720  // Before filling up the PointCloud, if the vectors are filled
721  // empty them
722  *yCloud = molSys::clearPointCloud(yCloud);
723 
724  // ----------------------------------------------------------
725  // If targetFrame has been found, read in the box lengths,
726  // number of atoms and then read in atom positions, type, molID
727  // By default, set molID=1 if not specified
728  if (foundFrame) {
729  // Run this until EOF or you reach the next timestep
730  while (std::getline((*dumpFile), line)) {
731  // Read in lines and tokenize them into std::string words and <double>
732  // numbers
733  tokens = gen::tokenizer(line);
734  numbers = gen::tokenizerDouble(line);
735 
736  // If you've reached the timestep line then you've reached the
737  // next frame. Break out of the while loop
738  if (tokens[0].compare("ITEM:") == 0) {
739  if (tokens[1].compare("TIMESTEP") == 0) {
740  break;
741  }
742  }
743 
744  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
745  // Read number of particles
746  if (readNOP) {
747  nop = std::stoi(line.data());
748  readNOP = false;
749  }
750  // Read box lengths
751  if (readBox) {
752  // You've reached the end of box lengths
753  if (tokens[0].compare("ITEM:") == 0) {
754  readBox = false;
755  // If the box is triclinic, get the
756  // orthogonal 'bounding box'
757  if (isTriclinic) {
758  // Update tilt factors
759  for (int k = 0; k < tilt.size(); k++) {
760  yCloud->box.push_back(tilt[k]);
761  }
762  } // end of check for triclinic
763  }
764  // Or else fill up the box lengths
765  else {
766  yCloud->box.push_back(numbers[1] - numbers[0]); // Update box length
767  yCloud->boxLow.push_back(
768  numbers[0]); // Update the lower box coordinate
769  // Do this for a triclinic box only
770  if (numbers.size() == 3) {
771  isTriclinic = true;
772  tilt.push_back(numbers[2]);
773  }
774  }
775  }
776  // Read atoms into yCloud line by line
777  if (readAtoms) {
778  iPoint.type = numbers[typeIndex];
779  iPoint.molID = numbers[molIndex];
780  iPoint.atomID = numbers[atomIndex];
781  iPoint.x = numbers[xIndex];
782  iPoint.y = numbers[yIndex];
783  iPoint.z = numbers[zIndex];
784  // Check if the particle is inside the volume Slice
785  // or not
786  if (isSlice) { // only if a slice has been requested
787  iPoint.inSlice = sinp::atomInSlice(iPoint.x, iPoint.y, iPoint.z,
788  coordLow, coordHigh);
789  // Skip if the atom is not part of the slice
790  if (!iPoint.inSlice) {
791  continue;
792  } // do not save if not inside the slice
793  }
794  // Save only atoms of the desired type
795  if (iPoint.type == typeI) {
796  nOxy++;
797  // yCloud->pts.resize(yCloud->pts.size()+1);
798  yCloud->pts.push_back(iPoint);
799  yCloud->idIndexMap[iPoint.atomID] =
800  yCloud->pts.size() - 1; // array index
801  }
802  }
803  // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
804 
805  // Tests for reading in nop, box lengths, and atoms
806  if (tokens[0].compare("ITEM:") == 0) {
807  if (tokens[1].compare("NUMBER") == 0) {
808  readNOP = true;
809  }
810  }
811  if (tokens[0].compare("ITEM:") == 0) {
812  if (tokens[1].compare("BOX") == 0) {
813  readBox = true;
814  }
815  }
816  if (tokens[0].compare("ITEM:") == 0) {
817  if (tokens[1].compare("ATOMS") == 0) {
818  readAtoms = true;
819  // Now find out which index is the coordinate index etc
820  for (int i = 2; i < tokens.size(); i++) {
821  if (tokens[i].compare("type") == 0) {
822  typeIndex = i - 2;
823  }
824  if (tokens[i].compare("x") == 0) {
825  xIndex = i - 2;
826  }
827  if (tokens[i].compare("y") == 0) {
828  yIndex = i - 2;
829  }
830  if (tokens[i].compare("z") == 0) {
831  zIndex = i - 2;
832  }
833  if (tokens[i].compare("mol") == 0) {
834  molIndex = i - 2;
835  }
836  if (tokens[i].compare("id") == 0) {
837  atomIndex = i - 2;
838  }
839  } // End of for loop over tokens
840  if (molIndex == 0) {
841  molIndex = atomIndex;
842  } // Set mol ID=atomID if not given
843  }
844  } // End of nested if loops for checking atom
845 
846  } // End of while
847  } // End of targetFrame found
848  // ----------------------------------------------------------
849  } // End of if file open statement
850 
851  // Check if you filled in the frame correctly
852  if (!(foundFrame)) {
853  std::cout << "You entered a frame that doesn't exist.\n";
854  } // Throw exception
855 
856  // Update the number of particles
857  yCloud->nop = yCloud->pts.size();
858 
859  // Update the frame number
860  yCloud->currentFrame = targetFrame;
861 
862  dumpFile->close();
863  return *yCloud;
864 }

◆ readXYZ()

int sinp::readXYZ ( std::string  filename,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud 
)

Function for reading in atom coordinates from an XYZ file.

Function for reading in an XYZ file

Definition at line 69 of file seams_input.cpp.

70  {
72  xyzFile = std::make_unique<std::ifstream>(filename);
73  std::string line; // Current line being read in
74  std::vector<std::string> tokens; // Vector containing word tokens
75  std::vector<double> numbers; // Vector containing type double numbers
76  int nop = -1; // Number of atoms in targetFrame
77  int iatom = 1; // Current atom being filled into the PointCloud
78  molSys::Point<double> iPoint; // Current point being read in from the file
79  double xLo, xHi, yLo, yHi, zLo,
80  zHi; // Box lengths extrapolated from the least and greatest coordinates
81 
82  if (!(gen::file_exists(filename))) {
83  std::cout
84  << "Fatal Error: The file does not exist or you gave the wrong path.\n";
85  // Throw exception?
86  return 1;
87  }
88 
89  // --------
90  // Before filling up the PointCloud, if the vectors are filled
91  // empty them
92  *yCloud = molSys::clearPointCloud(yCloud);
93  // --------
94 
95  // Format of an XYZ file:
96  // 1970
97  // generated by VMD
98  // O 43.603500 16.926201 15.215700
99  // O 39.912601 14.775100 19.379200
100  if (xyzFile->is_open()) {
101  // ----------------------------------------------------------
102  // At this point we know that the XYZ file is open
103 
104  // The first line contains the number of atoms
105  std::getline((*xyzFile), line);
106  nop = std::stoi(line);
107  // Skip the comment line
108  std::getline((*xyzFile), line);
109  // Reserve memory for atom coordinates using nop read in
110  yCloud->pts.reserve(nop);
111  yCloud->nop = nop;
112 
113  // Run this until EOF or you reach the next timestep
114  while (std::getline((*xyzFile), line)) {
115  // Read in lines and tokenize them into std::string words and <double>
116  // numbers
117  tokens = gen::tokenizer(line);
118  numbers = gen::tokenizerDouble(line);
119 
120  // Skip whitespace
121  if (tokens.size() == 0) {
122  continue;
123  }
124 
125  // Put logic for checking atom type here later
126  iPoint.type = 1; // Oxygen type; hard-coded!
127  iPoint.molID = iatom;
128  iPoint.atomID = iatom;
129  iPoint.x = std::stod(tokens[1]);
130  iPoint.y = std::stod(tokens[2]);
131  iPoint.z = std::stod(tokens[3]);
132 
133  yCloud->pts.push_back(iPoint);
134  yCloud->idIndexMap[iPoint.atomID] = yCloud->pts.size() - 1;
135  iatom++; // Increase index
136 
137  // First point
138  if (iatom == 1) {
139  // Loop through the dimensions
140  xLo = iPoint.x;
141  xHi = iPoint.x;
142  yLo = iPoint.y;
143  yHi = iPoint.y;
144  zLo = iPoint.z;
145  zHi = iPoint.z;
146  } // first point
147  else {
148  if (iPoint.x < xLo) {
149  xLo = iPoint.x;
150  } // xLo
151  else if (iPoint.x > xHi) {
152  xHi = iPoint.x;
153  } // xHi
154  else if (iPoint.y < yLo) {
155  yLo = iPoint.y;
156  } // yLo
157  else if (iPoint.y > yHi) {
158  yHi = iPoint.y;
159  } // yHi
160  else if (iPoint.z < zLo) {
161  zLo = iPoint.z;
162  } // zLo
163  else if (iPoint.z > zHi) {
164  zHi = iPoint.z;
165  } // zHi
166  } // update
167 
168  } // end of while, looping through lines till EOF
169  // ----------------------------------------------------------
170  } // End of if file open statement
171 
172  xyzFile->close();
173 
174  if (yCloud->pts.size() == 1) {
175  xHi = xLo + 10;
176  yHi = yLo + 10;
177  zHi = zLo + 10;
178  } // for a single point in the system (never happens)
179 
180  // Fill up the box lengths
181  yCloud->box.push_back(xHi - xLo);
182  yCloud->box.push_back(yHi - yLo);
183  yCloud->box.push_back(zHi - zLo);
184 
185  if (yCloud->pts.size() != yCloud->nop) {
186  std::cout << "Atoms didn't get filled in properly.\n";
187  }
188 
189  return 0;
190 }
T stod(T... args)