Rdf2

Functions

int rdf2::rdf2Danalysis_AA (std::string path, std::vector< double > *rdfValues, molSys::PointCloud< molSys::Point< double >, double > *yCloud, double cutoff, double binwidth, int firstFrame, int finalFrame)
 
std::vector< int > rdf2::sampleRDF_AA (molSys::PointCloud< molSys::Point< double >, double > *yCloud, double cutoff, double binwidth, int nbin)
 Samples the RDF histogram at every step. More...
 
int rdf2::normalizeRDF (int nopA, std::vector< double > *rdfValues, std::vector< int > histogram, double binwidth, int nbin, std::vector< double > volumeLengths, int nIter)
 Normalize the histogram. More...
 
std::vector< double > rdf2::getSystemLengths (molSys::PointCloud< molSys::Point< double >, double > *yCloud)
 Gets the lengths of the volume slice of the quasi-two-dimensional system. More...
 
double rdf2::getPlaneArea (std::vector< double > volumeLengths)
 Gets the plane area from the volume lengths vector. More...
 

Detailed Description

Function Documentation

◆ getPlaneArea()

double rdf2::getPlaneArea ( std::vector< double >  volumeLengths)

Gets the plane area from the volume lengths vector.

Calculates the plane area from the volume lengths. This is the product of the two largest dimensions of the quasi-two-dimensional system.

Parameters
[in]volumeLengthsA vector of the lengths of the volume slice or simulation box
Returns
The plane area of the two significant dimensions

Definition at line 260 of file rdf2d.cpp.

260  {
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 end(T... args)
T reverse(T... args)
T sort(T... args)

◆ getSystemLengths()

std::vector< double > rdf2::getSystemLengths ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud)

Gets the lengths of the volume slice of the quasi-two-dimensional system.

Calculates the lengths of the quasi-two-dimensional system. The smallest length is the 'height'.

Parameters
[in]yCloudThe molSys::PointCloud struct for the system.
Returns
The length (i.e. the 'height') of the smallest dimension of quasi-two-dimensional system.

Definition at line 208 of file rdf2d.cpp.

209  {
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
int nop
Current frame number.
Definition: mol_sys.hpp:169
std::vector< S > pts
Definition: mol_sys.hpp:167
T push_back(T... args)

◆ normalizeRDF()

int rdf2::normalizeRDF ( int  nopA,
std::vector< double > *  rdfValues,
std::vector< int >  histogram,
double  binwidth,
int  nbin,
std::vector< double >  volumeLengths,
int  nIter 
)

Normalize the histogram.

Normalizes the histogram and adds it to the RDF. The normalization requires the plane area and height.

Parameters
[in]nopAThe number of particles of type A.
[in]rdfValuesRadial distribution function values for all the frames, in the form of a vector.
[in]histogramThe histogram for the current frame.
[in]binwidthThe width of each bin for the RDF histogram.
[in]nbinThe number of bins for the RDF.
[in]volumeLengthsThe confining dimensions of the quasi-two-dimensional system, which may be the slice dimensions or the dimensions of the box.
[in]nIterThe number of iterations for which the coordinates will be binned. This is basically equivalent to the number of frames over which the RDF will be calculated.

Definition at line 158 of file rdf2d.cpp.

160  {
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
const double pi
Definition: generic.hpp:50
T min_element(T... args)
T pow(T... args)

◆ rdf2Danalysis_AA()

int rdf2::rdf2Danalysis_AA ( std::string  path,
std::vector< double > *  rdfValues,
molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
double  cutoff,
double  binwidth,
int  firstFrame,
int  finalFrame 
)

Main function for calculating the RDF for the same type of particle: calls other functions for initializing, sampling and normalizing the RDF

Calculates the in-plane RDF for quasi-two-dimensional water, when both the atoms are of the same type. The input PointCloud only has particles of type A in it. This is registered as a Lua function and is accessible to the user. Internally, this function calls the following functions:

  • rdf2::getSystemLengths (Gets the dimensions of the quasi-two-dimensional system).
  • ring::rdf2::sampleRDF_AA (Samples the current frame, binning the coordinates).
  • ring::rdf2::normalizeRDF (Normalizes the RDF).
  • ring::sout::printRDF (Writes out the RDF to the desired output directory, in the form of an ASCII file)
    Parameters
    [in]pathThe file path of the output directory to which output files will be written.
    [in]rdfValuesVector containing the RDF values.
    [in]yCloudThe input PointCloud.
    [in]cutoffCutoff for the RDF. This should not be greater than half the box length.
    [in]binwidthWidth of the bin for histogramming.
    [in]firstFrameThe first frame for RDF binning.
    [in]finalFrameThe final frame for RDF binning.

Definition at line 41 of file rdf2d.cpp.

44  {
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
T close(T... args)
int currentFrame
Collection of points.
Definition: mol_sys.hpp:168
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
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
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)

◆ sampleRDF_AA()

std::vector< int > rdf2::sampleRDF_AA ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
double  cutoff,
double  binwidth,
int  nbin 
)

Samples the RDF histogram at every step.

Samples the RDF for a particular frame The input PointCloud only has particles of type A in it.

  • gen::periodicDist (Periodic distance between a pair of atoms).
    Parameters
    [in]yCloudThe input PointCloud.
    [in]cutoffCutoff for the RDF calculation, which should be less than or equal to half the box length.
    [in]binwidthWidth of the bin.
    [in]nbinNumber of bins.
    Returns
    RDF histogram for the current frame.

Definition at line 111 of file rdf2d.cpp.

112  {
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
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
T resize(T... args)