generic.hpp

Go to the documentation of this file.

//-----------------------------------------------------------------------------------
// d-SEAMS - Deferred Structural Elucidation Analysis for Molecular Simulations
//
// Copyright (c) 2018--present d-SEAMS core team
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the MIT License as published by
// the Open Source Initiative.
//
// A copy of the MIT License is included in the LICENSE file of this repository.
// You should have received a copy of the MIT License along with this program.
// If not, see <https://opensource.org/licenses/MIT>.
//-----------------------------------------------------------------------------------
 
#ifndef SEAMS_GENERIC_H_
#define SEAMS_GENERIC_H_
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <filesystem>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <mol_sys.hpp>
 
#ifdef SEAMS_HAS_MINIMAGE
#include <minimage.h>
#endif
 
// C++20
#include <numbers>
// Eigen
#include <Eigen/Core>
#include <Eigen/Dense>



 
namespace gen {

constexpr double pi = std::numbers::pi;

inline double radDeg(double angle) { return (angle * 180) / gen::pi; }

double eigenVecAngle(std::vector<double> OO, std::vector<double> OH);

double getAverageWithoutOutliers(std::vector<double> inpVec);

inline double calcMedian(std::vector<double> *input) {
  int n = (*input).size(); // Number of elements
  double median;           // Output median value
 
  // Sort a copy (avoid mutating input)
  std::vector<double> sorted = *input;
  std::sort(sorted.begin(), sorted.end());
 
  // Calculate the median
  // For even values, the median is the average of the two middle values
  if (n % 2 == 0) {
    median = 0.5 * (sorted[n / 2] + sorted[n / 2 - 1]);
  } else {
    median = sorted[(n + 1) / 2 - 1];
  }
 
  return median;
}
 
#ifdef SEAMS_HAS_MINIMAGE
// Dump bound spans plus optional tilt onto an mi_cell.
inline mi_cell pointCloudCell(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud) {
  const auto &box = yCloud.box;
  const auto &boxLow = yCloud.boxLow;
  const double xlo_b = boxLow.size() > 0 ? boxLow[0] : 0.0;
  const double ylo_b = boxLow.size() > 1 ? boxLow[1] : 0.0;
  const double zlo_b = boxLow.size() > 2 ? boxLow[2] : 0.0;
  mi_cell c = mi_cell_ortho(0.0, 0.0, 0.0);
  if (box.size() >= 6) {
    mi_cell_from_lammps_bounds(box[0], box[1], box[2], box[3], box[4], box[5],
                               xlo_b, ylo_b, zlo_b, &c);
  } else if (box.size() >= 3) {
    c = mi_cell_ortho(box[0], box[1], box[2]);
    c.ox = xlo_b;
    c.oy = ylo_b;
    c.oz = zlo_b;
  }
  return c;
}
#endif
 
// Recover H (columns a, b, c) and origin from a LAMMPS dump box with the
// same mapping as nneigh::lammpsBoxToLcCell: box[0..3] are bound spans,
// box[3..6] are tilt factors xy, xz, yz, boxLow is the bound lo.
// Convert each cartesian point to fractional coordinates relative to
// origin, wrap the fractional difference, and map back through H.
inline std::array<double, 3> triclinicMinImage(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, double xi,
    double yi, double zi, double xj, double yj, double zj) {
#ifdef SEAMS_HAS_MINIMAGE
  const mi_cell cell = pointCloudCell(yCloud);
  const double p[3] = {xj, yj, zj};
  const double q[3] = {xi, yi, zi};
  double dr[3] = {0.0, 0.0, 0.0};
  mi_displacement(&cell, p, q, dr);
  return {dr[0], dr[1], dr[2]};
#else
  const auto &box = yCloud.box;
  const auto &boxLow = yCloud.boxLow;
  const double xspan = box[0];
  const double yspan = box[1];
  const double zspan = box[2];
  const double xlo_b = boxLow.size() > 0 ? boxLow[0] : 0.0;
  const double ylo_b = boxLow.size() > 1 ? boxLow[1] : 0.0;
  const double zlo_b = boxLow.size() > 2 ? boxLow[2] : 0.0;
  const double xy = box[3];
  const double xz = box[4];
  const double yz = box[5];
  const double xmin = std::min(std::min(0.0, xy), std::min(xz, xy + xz));
  const double xmax = std::max(std::max(0.0, xy), std::max(xz, xy + xz));
  const double ymin = std::min(0.0, yz);
  const double ymax = std::max(0.0, yz);
  const double lx = xspan - xmax + xmin;
  const double ly = yspan - ymax + ymin;
  const double lz = zspan;
  const double ox = xlo_b - xmin;
  const double oy = ylo_b - ymin;
  const double oz = zlo_b;
 
  auto toFrac = [&](double x, double y, double z) {
    const double sz = (z - oz) / lz;
    const double sy = (y - oy - yz * sz) / ly;
    const double sx = (x - ox - xy * sy - xz * sz) / lx;
    return std::array<double, 3>{sx, sy, sz};
  };
  const auto si = toFrac(xi, yi, zi);
  const auto sj = toFrac(xj, yj, zj);
  double dsx = si[0] - sj[0];
  double dsy = si[1] - sj[1];
  double dsz = si[2] - sj[2];
  dsx -= std::round(dsx);
  dsy -= std::round(dsy);
  dsz -= std::round(dsz);
  return {lx * dsx + xy * dsy + xz * dsz, ly * dsy + yz * dsz, lz * dsz};
#endif
}
 
// Generic function for getting the unwrapped distance
inline double
periodicDistSq(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
               int iatom, int jatom) {
#ifdef SEAMS_HAS_MINIMAGE
  const mi_cell cell = pointCloudCell(yCloud);
  const double p[3] = {yCloud.pts[iatom].x, yCloud.pts[iatom].y,
                       yCloud.pts[iatom].z};
  const double q[3] = {yCloud.pts[jatom].x, yCloud.pts[jatom].y,
                       yCloud.pts[jatom].z};
  double out = 0.0;
  mi_dist2(&cell, p, q, &out);
  return out;
#else
  if (yCloud.box.size() >= 6) {
    const auto dr = triclinicMinImage(
        yCloud, yCloud.pts[iatom].x, yCloud.pts[iatom].y, yCloud.pts[iatom].z,
        yCloud.pts[jatom].x, yCloud.pts[jatom].y, yCloud.pts[jatom].z);
    return dr[0] * dr[0] + dr[1] * dr[1] + dr[2] * dr[2];
  }
 
  std::array<double, 3> dr;
  double r2 = 0.0; // Squared absolute distance
 
  // Get x1-x2 etc
  dr[0] = std::fabs(yCloud.pts[iatom].x - yCloud.pts[jatom].x);
  dr[1] = std::fabs(yCloud.pts[iatom].y - yCloud.pts[jatom].y);
  dr[2] = std::fabs(yCloud.pts[iatom].z - yCloud.pts[jatom].z);
 
  // Three-axis wrap for an orthorhombic length-3 box
  for (int k = 0; k < 3; k++) {
    dr[k] -= yCloud.box[k] * std::round(dr[k] / yCloud.box[k]);
    r2 += dr[k] * dr[k];
  }
 
  return r2;
#endif
}

inline double
periodicDist(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
             int iatom, int jatom) {
  return std::sqrt(periodicDistSq(yCloud, iatom, jatom));
}
 
// Scalar pair batch. Highway BatchPeriodicDistSq is ortho-only; a
// tilt dump (box.size() >= 6) must go through periodicDistSq.
inline void batchPeriodicDistSq(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int iatom,
    const int *jatom, std::size_t n, double *distSq) {
  for (std::size_t k = 0; k < n; k++) {
    distSq[k] = periodicDistSq(yCloud, iatom, jatom[k]);
  }
}
 
// Bound spans, then tilt when box.size() >= 6.
inline std::string formatDumpBox(const std::vector<double> &box) {
  std::ostringstream oss;
  if (box.size() >= 3) {
    oss << box[0] << ' ' << box[1] << ' ' << box[2];
  }
  if (box.size() >= 6) {
    oss << " xy " << box[3] << " xz " << box[4] << " yz " << box[5];
  }
  return oss.str();
}
 
// ITEM line plus three bound lines. Tilt is a third field per line.
inline void writeDumpBoxBounds(
    std::ostream &os,
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud) {
  const bool tilt = yCloud.box.size() >= 6;
  if (tilt) {
    os << "ITEM: BOX BOUNDS xy xz yz pp pp pp\n";
  } else {
    os << "ITEM: BOX BOUNDS pp pp pp\n";
  }
  for (int k = 0; k < 3; k++) {
    const double lo = (static_cast<std::size_t>(k) < yCloud.boxLow.size())
                          ? yCloud.boxLow[k]
                          : 0.0;
    const double len = (static_cast<std::size_t>(k) < yCloud.box.size())
                           ? yCloud.box[k]
                           : 0.0;
    os << lo << ' ' << lo + len;
    if (tilt) {
      os << ' ' << yCloud.box[static_cast<std::size_t>(k + 3)];
    }
    os << '\n';
  }
}

inline std::array<double, 3> relDistFromPoint(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int iatom,
    double xj, double yj, double zj) {
#ifdef SEAMS_HAS_MINIMAGE
  const mi_cell cell = pointCloudCell(yCloud);
  const double p[3] = {xj, yj, zj};
  const double q[3] = {yCloud.pts[iatom].x, yCloud.pts[iatom].y,
                       yCloud.pts[iatom].z};
  double dr[3] = {0.0, 0.0, 0.0};
  mi_displacement(&cell, p, q, dr);
  return {dr[0], dr[1], dr[2]};
#else
  if (yCloud.box.size() >= 6) {
    return triclinicMinImage(yCloud, yCloud.pts[iatom].x, yCloud.pts[iatom].y,
                             yCloud.pts[iatom].z, xj, yj, zj);
  }
 
  std::array<double, 3> dr = {yCloud.pts[iatom].x - xj, yCloud.pts[iatom].y - yj,
                              yCloud.pts[iatom].z - zj};
  for (int k = 0; k < 3; k++) {
    if (dr[k] < -yCloud.box[k] * 0.5) {
      dr[k] += yCloud.box[k];
    }
    if (dr[k] >= yCloud.box[k] * 0.5) {
      dr[k] -= yCloud.box[k];
    }
  }
  return dr;
#endif
}
 
inline double unWrappedDistFromPoint(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int iatom,
    const std::vector<double> &singlePoint) {
  const auto dr =
      relDistFromPoint(yCloud, iatom, singlePoint[0], singlePoint[1],
                       singlePoint[2]);
  return std::sqrt(dr[0] * dr[0] + dr[1] * dr[1] + dr[2] * dr[2]);
}
 
// Generic function for getting the distance (no PBCs applied)
inline double
distance(const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int iatom,
         int jatom) {
  std::array<double, 3> dr;
  double r2 = 0.0; // Squared absolute distance
 
  // Get x1-x2 etc
  dr[0] = fabs(yCloud.pts[iatom].x - yCloud.pts[jatom].x);
  dr[1] = fabs(yCloud.pts[iatom].y - yCloud.pts[jatom].y);
  dr[2] = fabs(yCloud.pts[iatom].z - yCloud.pts[jatom].z);
 
  // Get the squared absolute distance
  for (int k = 0; k < 3; k++) {
    r2 += pow(dr[k], 2.0);
  }
 
  return sqrt(r2);
}
 
// Generic function for getting the relative coordinates
inline std::array<double, 3>
relDist(const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int iatom,
        int jatom) {
  return relDistFromPoint(yCloud, iatom, yCloud.pts[jatom].x,
                          yCloud.pts[jatom].y, yCloud.pts[jatom].z);
}
 
// Function for sorting according to atom ID
// Comparator for std::sort
inline bool compareByAtomID(const molSys::Point<double> &a,
                            const molSys::Point<double> &b) {
  return a.atomID < b.atomID;
}

[[nodiscard]] int prettyPrintYoda(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                    std::string outFile);

[[nodiscard]] int unwrappedCoordShift(
    const 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);

double angDistDegQuaternions(std::vector<double> quat1,
                             std::vector<double> quat2);

inline std::vector<std::string> tokenizer(std::string line) {
  std::istringstream iss(line);
  std::vector<std::string> tokens{std::istream_iterator<std::string>{iss},
                                  std::istream_iterator<std::string>{}};
  return tokens;
}

inline std::vector<double> tokenizerDouble(std::string line) {
  std::istringstream iss(line);
  std::vector<double> tokens;
  double number; // Each number being read in from the line
  while (iss >> number) {
    tokens.push_back(number);
  }
  return tokens;
}

inline std::vector<int> tokenizerInt(std::string line) {
  std::istringstream iss(line);
  std::vector<int> tokens;
  int number; // Each number being read in from the line
  while (iss >> number) {
    tokens.push_back(number);
  }
  return tokens;
}

inline bool file_exists(const std::string &name) {
  return std::filesystem::exists(name);
}

inline std::vector<std::complex<double>>
avgVector(std::vector<std::complex<double>> v, int l, int neigh) {
  if (neigh == 0) {
    return v;
  }
  for (int m = 0; m < 2 * l + 1; m++) {
    v[m] = (1.0 / static_cast<double>(neigh)) * v[m];
  }
 
  return v;
}
 
} // namespace gen
 
#endif // SEAMS_GENERIC_H_