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 264 of file rdf2d.cpp.

264  {
265  //
266  // Sort the vector in descending order
267  std::sort(volumeLengths.begin(), volumeLengths.end());
268  std::reverse(volumeLengths.begin(), volumeLengths.end());
269 
270  return volumeLengths[0] * volumeLengths[1];
271 
272 } // 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 212 of file rdf2d.cpp.

213  {
214  //
215  std::vector<double> lengths; // Volume lengths
217  rMax; // Max of the coordinates {0 is for x, 1 is for y and 2 is for z}
219  rMin; // Min of the coordinates {0 is for x, 1 is for y and 2 is for z}
220  std::vector<double> r_iatom; // Current point coordinates
221  int dim = 3;
222 
223  // Init
224  r_iatom.push_back(yCloud->pts[0].x);
225  r_iatom.push_back(yCloud->pts[0].y);
226  r_iatom.push_back(yCloud->pts[0].z);
227  rMin = r_iatom;
228  rMax = r_iatom;
229 
230  // Loop through the entire PointCloud
231  for (int iatom = 1; iatom < yCloud->nop; iatom++) {
232  r_iatom[0] = yCloud->pts[iatom].x; // x coordinate of iatom
233  r_iatom[1] = yCloud->pts[iatom].y; // y coordinate of iatom
234  r_iatom[2] = yCloud->pts[iatom].z; // z coordinate of iatom
235  // Loop through every dimension
236  for (int k = 0; k < dim; k++) {
237  // Update rMin
238  if (r_iatom[k] < rMin[k]) {
239  rMin[k] = r_iatom[k];
240  } // end of rMin update
241  // Update rMax
242  if (r_iatom[k] > rMax[k]) {
243  rMax[k] = r_iatom[k];
244  } // end of rMax update
245  } // end of looping through every dimension
246  } // end of loop through all the atoms
247 
248  // Get the lengths
249  for (int k = 0; k < dim; k++) {
250  lengths.push_back(rMax[k] - rMin[k]);
251  } // end of updating lengths
252 
253  return lengths;
254 } // end of function
int nop
Current frame number.
Definition: mol_sys.hpp:173
std::vector< S > pts
Definition: mol_sys.hpp:171
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 162 of file rdf2d.cpp.

164  {
165  //
166  auto height = *std::min_element(
167  volumeLengths.begin(),
168  volumeLengths.end()); // The height or the smallest dimension
169  // double planeArea = rdf2::getPlaneArea(
170  // volumeLengths); // The area of the quasi-two-dimensional water
171  double r; // Distance for the current bin
172  double factor; // Factor for accounting for the two-dimensional slab
173  double binVolume; // Volume of the current bin
174  double volumeDensity =
175  nopA / (volumeLengths[0] * volumeLengths[1] * volumeLengths[2]);
176 
177  // Loop through every bin
178  for (int ibin = 0; ibin < nbin; ibin++) {
179  //
180  r = binwidth * (ibin + 0.5);
181  // ----------------------
182  // Factor calculation
183  factor = 1;
184  if (r > height) {
185  // FACTOR = HEIGHT/(2*R)
186  factor = height / (2 * r);
187  } // r>height
188  else {
189  // FACTOR = 1 - R/(2*HEIGHT)
190  factor = 1 - r / (2 * height);
191  } // r<=height
192  // ----------------------
193  // binVolume = 4.0*PI*(DELTA_R**3)*((I_BIN+1)**3-I_BIN**3)/3.0
194  binVolume = 4 * gen::pi * pow(binwidth, 3.0) *
195  (pow(ibin + 1, 3.0) - pow(ibin, 3.0)) / 3.0;
196  // Update the RDF
197  (*rdfValues)[ibin] +=
198  histogram[ibin] / (nIter * binVolume * nopA * volumeDensity * factor);
199  } // end of loop through every bin
200 
201  return 0;
202 
203 } // end of function
const double pi
Definition: generic.hpp:54
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 45 of file rdf2d.cpp.

48  {
49  //
50  int nopA = yCloud->nop; // Number of particles of type A in the current frame.
51  int currentFrame = yCloud->currentFrame; // The current frame
52  int nbin; // Number of bins
53  std::vector<double> volumeLengths; // Lengths of the volume of the
54  // quasi-two-dimensional system
55  std::vector<int> histogram; // Histogram for the RDF, at every step
56  int nIter = finalFrame - firstFrame + 1; // Number of iterations
57  // ----------------------------------------------
58  // INITIALIZATION
59  if (currentFrame == firstFrame) {
60  // -----------------
61  // Checks and balances?
62  // -----------------
63  nbin = 1 + int(cutoff / binwidth); // Number of bins
64  (*rdfValues).resize(nbin); // RDF initialized to 0
65  } // end of initialization
66  // ----------------------------------------------
67  // SAMPLING
68  // Get the volume lengths of the quasi-two-dimensional system
69  volumeLengths = rdf2::getSystemLengths(yCloud);
70  // Update cutoff if required? later
71  nbin = int(cutoff / binwidth); // Number of bins
72  // Sample for the current frame!!
73  histogram = rdf2::sampleRDF_AA(yCloud, cutoff, binwidth, nbin);
74  // ----------------------------------------------
75  // NORMALIZATION
76  rdf2::normalizeRDF(nopA, rdfValues, histogram, binwidth, nbin, volumeLengths,
77  nIter);
78  // ----------------------------------------------
79  // PRINT OUT
80  if (currentFrame == finalFrame) {
81  // Create folder if required
82  sout::makePath(path);
83  std::string outputDirName = path + "topoMonolayer";
84  sout::makePath(outputDirName);
85  //
86  std::string fileName = path + "topoMonolayer/rdf.dat";
87  //
88  // Comment line
89  std::ofstream outputFile; // For the output file
90  outputFile.open(fileName, std::ios_base::app | std::ios_base::out);
91  outputFile << "# r g(r)\n";
92  outputFile.close();
93  //
94  //
95  sout::printRDF(fileName, rdfValues, binwidth, nbin);
96  } // end of print out
97  // ----------------------------------------------
98 
99  return 0;
100 } // end of function
T close(T... args)
int currentFrame
Collection of points.
Definition: mol_sys.hpp:172
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:212
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:115
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:162
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 115 of file rdf2d.cpp.

116  {
117  //
118  std::vector<int> histogram; // Histogram for the RDF
119  double r_ij; // Unwrapped distance between iatom and jatom
120  int ibin; // Index of the bin in which the distance r_ij falls
121 
122  // Init the histogram to 0
123  histogram.resize(nbin);
124 
125  // Loop through pairs of atoms
126  for (int iatom = 0; iatom < yCloud->nop - 1; iatom++) {
127  //
128  // Loop through the next atom in the atom pair
129  for (int jatom = iatom + 1; jatom < yCloud->nop; jatom++) {
130  // Calculate the distance between iatom and jatom
131  r_ij = gen::periodicDist(yCloud, iatom, jatom);
132  // Check if the distance is within the cutoff. Update the histogram if
133  // r_ij is within the cutoff
134  if (r_ij <= cutoff) {
135  ibin = int(r_ij / binwidth); // Bin in which r_ij falls
136  histogram[ibin] += 2; // Update for iatom and jatom both
137  } // end of histogram update
138  } // end of loop through jatom
139  } // end of loop through iatom
140 
141  // Return the histogram
142  return histogram;
143 
144 } // 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:108
T resize(T... args)