rdf2d.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <rdf2d.hpp>
12 
13 // -----------------------------------------------------------------------------------------------------
14 // IN-PLANE RDF
15 // -----------------------------------------------------------------------------------------------------
16 
42  std::string path, std::vector<double> *rdfValues,
43  molSys::PointCloud<molSys::Point<double>, double> *yCloud, double cutoff,
44  double binwidth, int firstFrame, int finalFrame) {
45  //
46  int nopA = yCloud->nop; // Number of particles of type A in the current frame.
47  int currentFrame = yCloud->currentFrame; // The current frame
48  int nbin; // Number of bins
49  std::vector<double> volumeLengths; // Lengths of the volume of the
50  // quasi-two-dimensional system
51  std::vector<int> histogram; // Histogram for the RDF, at every step
52  int nIter = finalFrame - firstFrame + 1; // Number of iterations
53  // ----------------------------------------------
54  // INITIALIZATION
55  if (currentFrame == firstFrame) {
56  // -----------------
57  // Checks and balances?
58  // -----------------
59  nbin = 1 + int(cutoff / binwidth); // Number of bins
60  (*rdfValues).resize(nbin); // RDF initialized to 0
61  } // end of initialization
62  // ----------------------------------------------
63  // SAMPLING
64  // Get the volume lengths of the quasi-two-dimensional system
65  volumeLengths = rdf2::getSystemLengths(yCloud);
66  // Update cutoff if required? later
67  nbin = int(cutoff / binwidth); // Number of bins
68  // Sample for the current frame!!
69  histogram = rdf2::sampleRDF_AA(yCloud, cutoff, binwidth, nbin);
70  // ----------------------------------------------
71  // NORMALIZATION
72  rdf2::normalizeRDF(nopA, rdfValues, histogram, binwidth, nbin, volumeLengths,
73  nIter);
74  // ----------------------------------------------
75  // PRINT OUT
76  if (currentFrame == finalFrame) {
77  // Create folder if required
78  sout::makePath(path);
79  std::string outputDirName = path + "topoMonolayer";
80  sout::makePath(outputDirName);
81  //
82  std::string fileName = path + "topoMonolayer/rdf.dat";
83  //
84  // Comment line
85  std::ofstream outputFile; // For the output file
86  outputFile.open(fileName, std::ios_base::app | std::ios_base::out);
87  outputFile << "# r g(r)\n";
88  outputFile.close();
89  //
90  //
91  sout::printRDF(fileName, rdfValues, binwidth, nbin);
92  } // end of print out
93  // ----------------------------------------------
94 
95  return 0;
96 } // end of function
97 
112  double cutoff, double binwidth, int nbin) {
113  //
114  std::vector<int> histogram; // Histogram for the RDF
115  double r_ij; // Unwrapped distance between iatom and jatom
116  int ibin; // Index of the bin in which the distance r_ij falls
117 
118  // Init the histogram to 0
119  histogram.resize(nbin);
120 
121  // Loop through pairs of atoms
122  for (int iatom = 0; iatom < yCloud->nop - 1; iatom++) {
123  //
124  // Loop through the next atom in the atom pair
125  for (int jatom = iatom + 1; jatom < yCloud->nop; jatom++) {
126  // Calculate the distance between iatom and jatom
127  r_ij = gen::periodicDist(yCloud, iatom, jatom);
128  // Check if the distance is within the cutoff. Update the histogram if
129  // r_ij is within the cutoff
130  if (r_ij <= cutoff) {
131  ibin = int(r_ij / binwidth); // Bin in which r_ij falls
132  histogram[ibin] += 2; // Update for iatom and jatom both
133  } // end of histogram update
134  } // end of loop through jatom
135  } // end of loop through iatom
136 
137  // Return the histogram
138  return histogram;
139 
140 } // end of function
141 
158 int rdf2::normalizeRDF(int nopA, std::vector<double> *rdfValues,
159  std::vector<int> histogram, double binwidth, int nbin,
160  std::vector<double> volumeLengths, int nIter) {
161  //
162  auto height = *std::min_element(
163  volumeLengths.begin(),
164  volumeLengths.end()); // The height or the smallest dimension
165  // double planeArea = rdf2::getPlaneArea(
166  // volumeLengths); // The area of the quasi-two-dimensional water
167  double r; // Distance for the current bin
168  double factor; // Factor for accounting for the two-dimensional slab
169  double binVolume; // Volume of the current bin
170  double volumeDensity =
171  nopA / (volumeLengths[0] * volumeLengths[1] * volumeLengths[2]);
172 
173  // Loop through every bin
174  for (int ibin = 0; ibin < nbin; ibin++) {
175  //
176  r = binwidth * (ibin + 0.5);
177  // ----------------------
178  // Factor calculation
179  factor = 1;
180  if (r > height) {
181  // FACTOR = HEIGHT/(2*R)
182  factor = height / (2 * r);
183  } // r>height
184  else {
185  // FACTOR = 1 - R/(2*HEIGHT)
186  factor = 1 - r / (2 * height);
187  } // r<=height
188  // ----------------------
189  // binVolume = 4.0*PI*(DELTA_R**3)*((I_BIN+1)**3-I_BIN**3)/3.0
190  binVolume = 4 * gen::pi * pow(binwidth, 3.0) *
191  (pow(ibin + 1, 3.0) - pow(ibin, 3.0)) / 3.0;
192  // Update the RDF
193  (*rdfValues)[ibin] +=
194  histogram[ibin] / (nIter * binVolume * nopA * volumeDensity * factor);
195  } // end of loop through every bin
196 
197  return 0;
198 
199 } // end of function
200 
209  molSys::PointCloud<molSys::Point<double>, double> *yCloud) {
210  //
211  std::vector<double> lengths; // Volume lengths
213  rMax; // Max of the coordinates {0 is for x, 1 is for y and 2 is for z}
215  rMin; // Min of the coordinates {0 is for x, 1 is for y and 2 is for z}
216  std::vector<double> r_iatom; // Current point coordinates
217  int dim = 3;
218 
219  // Init
220  r_iatom.push_back(yCloud->pts[0].x);
221  r_iatom.push_back(yCloud->pts[0].y);
222  r_iatom.push_back(yCloud->pts[0].z);
223  rMin = r_iatom;
224  rMax = r_iatom;
225 
226  // Loop through the entire PointCloud
227  for (int iatom = 1; iatom < yCloud->nop; iatom++) {
228  r_iatom[0] = yCloud->pts[iatom].x; // x coordinate of iatom
229  r_iatom[1] = yCloud->pts[iatom].y; // y coordinate of iatom
230  r_iatom[2] = yCloud->pts[iatom].z; // z coordinate of iatom
231  // Loop through every dimension
232  for (int k = 0; k < dim; k++) {
233  // Update rMin
234  if (r_iatom[k] < rMin[k]) {
235  rMin[k] = r_iatom[k];
236  } // end of rMin update
237  // Update rMax
238  if (r_iatom[k] > rMax[k]) {
239  rMax[k] = r_iatom[k];
240  } // end of rMax update
241  } // end of looping through every dimension
242  } // end of loop through all the atoms
243 
244  // Get the lengths
245  for (int k = 0; k < dim; k++) {
246  lengths.push_back(rMax[k] - rMin[k]);
247  } // end of updating lengths
248 
249  return lengths;
250 } // end of function
251 
261  //
262  // Sort the vector in descending order
263  std::sort(volumeLengths.begin(), volumeLengths.end());
264  std::reverse(volumeLengths.begin(), volumeLengths.end());
265 
266  return volumeLengths[0] * volumeLengths[1];
267 
268 } // end of function
T begin(T... args)
T close(T... args)
T end(T... args)
double periodicDist(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatom, int jatom)
Inline generic function for obtaining the unwrapped periodic distance between two particles,...
Definition: generic.hpp:104
const double pi
Definition: generic.hpp:50
std::vector< double > getSystemLengths(molSys::PointCloud< molSys::Point< double >, double > *yCloud)
Gets the lengths of the volume slice of the quasi-two-dimensional system.
Definition: rdf2d.cpp:208
std::vector< int > sampleRDF_AA(molSys::PointCloud< molSys::Point< double >, double > *yCloud, double cutoff, double binwidth, int nbin)
Samples the RDF histogram at every step.
Definition: rdf2d.cpp:111
double getPlaneArea(std::vector< double > volumeLengths)
Gets the plane area from the volume lengths vector.
Definition: rdf2d.cpp:260
int rdf2Danalysis_AA(std::string path, std::vector< double > *rdfValues, molSys::PointCloud< molSys::Point< double >, double > *yCloud, double cutoff, double binwidth, int firstFrame, int finalFrame)
Definition: rdf2d.cpp:41
int normalizeRDF(int nopA, std::vector< double > *rdfValues, std::vector< int > histogram, double binwidth, int nbin, std::vector< double > volumeLengths, int nIter)
Normalize the histogram.
Definition: rdf2d.cpp:158
T min_element(T... args)
int makePath(const std::string &path)
int printRDF(std::string fileName, std::vector< double > *rdfValues, double binwidth, int nbin)
Function for printing out the RDF, given the filename.
T open(T... args)
T push_back(T... args)
File containing functions used to calculate the in-plane radial distribution functions.
T resize(T... args)
T reverse(T... args)
T sort(T... args)
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166
This contains per-particle information.
Definition: mol_sys.hpp:145