order_parameter.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 <order_parameter.hpp>
12 
29  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int nPrisms,
30  double avgPrismHeight) {
31  //
32  double hPercent; // Normalized height percent
33  double nanoTubeHeight; // Height of the SWCT
34  double numberMax; // Maximum number possible, given the average prism height
35 
36  // ---------------------------------------
37  // Calculate the height of the SWCT
38  // This is the longest dimension of the simulation box
39  nanoTubeHeight = *max_element(yCloud->box.begin(), yCloud->box.end());
40  // ---------------------------------------
41  // Calculate the maximum possible height, given the average prism height
42  // and the height of the nanotube
43  numberMax = nanoTubeHeight / avgPrismHeight;
44  // ---------------------------------------
45  // Calculate the normalized height percentage
46  hPercent = nPrisms / numberMax * 100.0;
47 
48  return hPercent;
49 }
50 
68  std::vector<std::vector<int>> rings, double sheetArea) {
69  //
70  double areaXY, areaXZ, areaYZ; // Total coverage area
71  std::vector<double> singleAreas; // Area of single rings
72 
73  // ---------------------------------------
74  // Initialization
75  areaXY = 0.0;
76  areaXZ = 0.0;
77  areaYZ = 0.0;
78  // ---------------------------------------
79  // Loop through all the rings
80  for (int iring = 0; iring < rings.size(); iring++) {
81  // Get the coverage area for the current ring
82  singleAreas = topoparam::projAreaSingleRing(yCloud, rings[iring]);
83  // Add these to the total coverage area
84  areaXY += singleAreas[0];
85  areaXZ += singleAreas[1];
86  areaYZ += singleAreas[2];
87  } // end of loop through all the rings
88  // ---------------------------------------
89  // Normalize the coverage area by the sheet area
90  areaXY = areaXY / sheetArea * 100.0;
91  areaXZ = areaXZ / sheetArea * 100.0;
92  areaYZ = areaYZ / sheetArea * 100.0;
93 
94  return {areaXY, areaXZ, areaYZ};
95 }
96 
102  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
104  //
105  int iatomIndex, jatomIndex; // Atom indices of the i^th and j^th atoms
106  int ringSize = ring.size(); // Number of nodes in the ring
107  double areaXY, areaXZ, areaYZ;
108  double x_iatom, y_iatom, z_iatom; // Coordinates of iatom
109  double x_jatom, y_jatom, z_jatom; // Coordinates of jatom
110  // ----------------------------------------
111  // Calculate projected area onto the XY, YZ and XZ planes for basal1
112 
113  // Init the projected area
114  areaXY = 0.0;
115  areaXZ = 0.0;
116  areaYZ = 0.0;
117 
118  jatomIndex = ring[0];
119 
120  // All points except the first pair
121  for (int k = 1; k < ringSize; k++) {
122  iatomIndex = ring[k]; // Current vertex
123 
124  // --------------------------------------------------------------------
125  // SHIFT PARTICLES TEMPORARILY (IN CASE OF UNWRAPPED COORDINATES)
126  gen::unwrappedCoordShift(yCloud, iatomIndex, jatomIndex, &x_iatom, &y_iatom,
127  &z_iatom, &x_jatom, &y_jatom, &z_jatom);
128  // --------------------------------------------------------------------
129 
130  // Add to the polygon area
131  // ------
132  // XY plane
133  areaXY += (x_jatom + x_iatom) * (y_jatom - y_iatom);
134  // ------
135  // XZ plane
136  areaXZ += (x_jatom + x_iatom) * (z_jatom - z_iatom);
137  // ------
138  // YZ plane
139  areaYZ += (y_jatom + y_iatom) * (z_jatom - z_iatom);
140  // ------
141  jatomIndex = iatomIndex;
142  }
143 
144  // Closure point
145  iatomIndex = ring[0];
146  // Unwrapped coordinates needed
147  gen::unwrappedCoordShift(yCloud, iatomIndex, jatomIndex, &x_iatom, &y_iatom,
148  &z_iatom, &x_jatom, &y_jatom, &z_jatom);
149  // ------
150  // XY plane
151  areaXY += (x_jatom + x_iatom) * (y_jatom - y_iatom);
152  // ------
153  // XZ plane
154  areaXZ += (x_jatom + x_iatom) * (z_jatom - z_iatom);
155  // ------
156  // YZ plane
157  areaYZ += (y_jatom + y_iatom) * (z_jatom - z_iatom);
158  // ------
159  // The actual projected area is half of this
160  areaXY *= 0.5;
161  areaXZ *= 0.5;
162  areaYZ *= 0.5;
163 
164  // Absolute area
165  areaXY = fabs(areaXY);
166  areaXZ = fabs(areaXZ);
167  areaYZ = fabs(areaYZ);
168 
169  return {areaXY, areaYZ, areaXZ};
170 }
int unwrappedCoordShift(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatomIndex, int jatomIndex, double *x_i, double *y_i, double *z_i, double *x_j, double *y_j, double *z_j)
Shift particles (unwrapped coordinates)
Definition: generic.cpp:66
Topological network criteria functions.
Definition: ring.hpp:60
std::vector< double > projAreaSingleRing(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > ring)
Calculates the projected area on the XY, YZ and XZ planes.
double normHeightPercent(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int nPrisms, double avgPrismHeight)
std::vector< double > calcCoverageArea(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> rings, double sheetArea)
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