Functions
Sinp

Functions

std::vector< std::string > sinp::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})
 
molSys::PointCloud< molSys::Point< double >, double > sinp::readXYZ (std::string filename)
 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 88 of file seams_input.hpp.

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

◆ 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 24 of file seams_input.cpp.

24  {
25  std::unique_ptr<std::ifstream> inpFile;
26  inpFile = std::make_unique<std::ifstream>(filename);
27  std::vector<std::vector<int>> bonds;
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 }
std::getline
T getline(T... args)
gen::file_exists
bool file_exists(const std::string &name)
Function for checking if a file exists or not.
Definition: generic.hpp:298
gen::tokenizer
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::push_back
T push_back(T... args)
std::vector::size
T size(T... args)
std::stoi
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 212 of file seams_input.cpp.

215  {
216  std::unique_ptr<std::ifstream> dumpFile;
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 }
std::string::data
T data(T... args)
gen::tokenizerDouble
std::vector< double > tokenizerDouble(std::string line)
Function for tokenizing line strings into a vector of doubles.
Definition: generic.hpp:270
molSys::Point::x
T x
type ID, molID, atomID
Definition: mol_sys.hpp:151
molSys::PointCloud::nop
int nop
Current frame number.
Definition: mol_sys.hpp:173
molSys::Point::y
T y
Definition: mol_sys.hpp:151
molSys::Point::z
T z
Definition: mol_sys.hpp:151
molSys::Point::type
int type
Definition: mol_sys.hpp:150
molSys::PointCloud::pts
std::vector< S > pts
Definition: mol_sys.hpp:171
molSys::Point::inSlice
bool inSlice
Type of ice/water etc based on cij.
Definition: mol_sys.hpp:155
molSys::PointCloud::idIndexMap
std::unordered_map< int, int > idIndexMap
xlo, ylo, zlo
Definition: mol_sys.hpp:176
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::clearPointCloud
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::Point::molID
int molID
Definition: mol_sys.hpp:150
molSys::Point::atomID
int atomID
Definition: mol_sys.hpp:150
molSys::PointCloud::currentFrame
int currentFrame
Collection of points.
Definition: mol_sys.hpp:172
sinp::atomInSlice
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::reserve
T reserve(T... args)
molSys::Point
This contains per-particle information.
Definition: mol_sys.hpp:149

◆ 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 431 of file seams_input.cpp.

434  {
435  std::unique_ptr<std::ifstream> dumpFile;
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 }

◆ 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 656 of file seams_input.cpp.

660  {
661  std::unique_ptr<std::ifstream> dumpFile;
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 }

◆ readXYZ()

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

Function for reading in atom coordinates from an XYZ file.

Function for reading in an XYZ file

Definition at line 73 of file seams_input.cpp.

73  {
74  std::unique_ptr<std::ifstream> xyzFile;
75  xyzFile = std::make_unique<std::ifstream>(filename);
76  molSys::PointCloud<molSys::Point<double>, double>
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 }
std::stod
T stod(T... args)
molSys::PointCloud
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:170
Generated by Doxygen 1.9.1 and the doxyYoda theme