26 outputFile.
open(outFile);
30 outputFile <<
"# Frame\tAtomID\tx\ty\tz\tcij\ticeType\n";
32 for (
int i = 0; i < yCloud->nop; i++) {
33 outputFile << yCloud->currentFrame <<
"\t" << yCloud->pts[i].atomID
34 <<
"\t" << yCloud->pts[i].x <<
"\t" << yCloud->pts[i].y <<
"\t"
35 << yCloud->pts[i].z <<
"\t";
39 outputFile << yCloud->pts[i].iceType <<
"\n";
68 int jatomIndex,
double *x_i,
double *y_i,
double *z_i,
double *x_j,
69 double *y_j,
double *z_j) {
71 double x_iatom, y_iatom, z_iatom;
72 double x_jatom, y_jatom, z_jatom;
73 double x_ij, y_ij, z_ij;
75 double xPBC, yPBC, zPBC;
80 x_iatom = yCloud->pts[iatomIndex].x;
81 y_iatom = yCloud->pts[iatomIndex].y;
82 z_iatom = yCloud->pts[iatomIndex].z;
84 x_jatom = yCloud->pts[jatomIndex].x;
85 y_jatom = yCloud->pts[jatomIndex].y;
86 z_jatom = yCloud->pts[jatomIndex].z;
89 x_ij = x_iatom - x_jatom;
90 y_ij = y_iatom - y_jatom;
91 z_ij = z_iatom - z_jatom;
95 if (fabs(x_ij) > 0.5 * box[0]) {
97 xPBC = box[0] - fabs(x_ij);
99 x_jatom = x_iatom - xPBC;
102 x_jatom = x_iatom + xPBC;
107 if (fabs(y_ij) > 0.5 * box[1]) {
109 yPBC = box[1] - fabs(y_ij);
111 y_jatom = y_iatom - yPBC;
114 y_jatom = y_iatom + yPBC;
119 if (fabs(z_ij) > 0.5 * box[2]) {
121 zPBC = box[2] - fabs(z_ij);
123 z_jatom = z_iatom - zPBC;
126 z_jatom = z_iatom + zPBC;
151 Eigen::Vector3d eigOO = Eigen::Map<Eigen::Vector3d>(OO.
data(), OO.
size());
152 Eigen::Vector3d eigOH = Eigen::Map<Eigen::Vector3d>(OH.
data(), OH.
size());
154 angle = acos(eigOO.dot(eigOH) / (eigOO.norm() * eigOH.norm()));
168 int n = inpVec.
size();
171 double firstQuartile, thirdQuartile;
173 double outlierLimLow, outlierLimHigh;
174 int numOfObservations = 0;
182 for (
int i = 0; i < n / 2; i++) {
191 int halfN = (n + 1) / 2;
193 for (
int i = 0; i < halfN; i++) {
208 iqr = thirdQuartile - firstQuartile;
216 outlierLimLow = firstQuartile - 1.5 * iqr;
217 outlierLimHigh = thirdQuartile + 1.5 * iqr;
220 for (
int i = 0; i < n; i++) {
222 if (inpVec[i] < outlierLimLow) {
225 else if (inpVec[i] > outlierLimHigh) {
237 avgVal /= numOfObservations;
240 if (numOfObservations == 0) {
243 for (
int i = 0; i <= n; i++) {
268 prod = quat1[0] * quat2[0] - quat1[1] * quat2[1] - quat1[2] * quat2[2] -
271 double angDist = 2 * acos(prod) * 180.0 / (
gen::pi);
File for containing generic or common functions.
double eigenVecAngle(std::vector< double > OO, std::vector< double > OH)
Eigen function for getting the angle (in radians) between the O–O and O-H vectors.
double calcMedian(std::vector< double > *input)
Inline generic function for calculating the median given a vector of the values.
double angDistDegQuaternions(std::vector< double > quat1, std::vector< double > quat2)
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)
int prettyPrintYoda(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string outFile)
Generic function for printing all the struct information.
double getAverageWithoutOutliers(std::vector< double > inpVec)
Get the average, after excluding the outliers, using quartiles.
This contains a collection of points; contains information for a particular frame.
This contains per-particle information.