bop.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_BOP_H_
#define SEAMS_BOP_H_
 
#include <array>
#include <cmath>
#include <complex>
#include <string>
#include <generic.hpp>
#include <mol_sys.hpp>
#include <neighbours.hpp>
#include <seams_output.hpp>



 
namespace chill {
 
// 2*l+1 length complex vector
struct YlmAtom {
  std::vector<std::complex<double>> ylm;
};
 
// Vector of 2*l+1 averaged over 4 nearest neighbours
struct QlmAtom {
  std::vector<YlmAtom> ptq; // Averaged over neighbours
};

struct BondClassifier {
  double staggeredMax;    
  double eclipsedMin;     
  double eclipsedMax;     
  int coordinationNumber; 
};

[[nodiscard]] BondClassifier chillRule();

[[nodiscard]] BondClassifier chillPlusRule();

[[nodiscard]] BondClassifier bondClassifier(const std::string &name);

void registerBondClassifier(const std::string &name,
                            const BondClassifier &rule);

[[nodiscard]] std::vector<std::string> bondClassifierNames();

void classifyBonds(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                   const std::vector<std::vector<int>> &nList,
                   const BondClassifier &rule, bool isSlice = false);

void
getCorrel(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
          const std::vector<std::vector<int>> &nList, bool isSlice = false,
          int coordinationNumber = 4);

void
getIceTypeNoPrint(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                  const std::vector<std::vector<int>> &nList, bool isSlice = false);
 
// Classifies each atom according to the CHILL algorithm
void
getIceType(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
           const std::vector<std::vector<int>> &nList, std::string path,
           int firstFrame, bool isSlice = false,
           std::string outputFileName = "chill.txt");

void
getCorrelPlus(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
              const std::vector<std::vector<int>> &nList, bool isSlice = false,
              int coordinationNumber = 4);

void
getIceTypePlus(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
               const std::vector<std::vector<int>> &nList, std::string path,
               int firstFrame, bool isSlice = false,
               std::string outputFileName = "chillPlus.txt");

void
getIceTypePlusNoPrint(
    molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    const std::vector<std::vector<int>> &nList, bool isSlice = false);

std::vector<double>
getq6(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
      const std::vector<std::vector<int>> &nList, bool isSlice = false);

void
reclassifyWater(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                std::vector<double> &q6);

[[nodiscard]] int printIceType(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                 std::string path, int firstFrame, bool isSlice = false,
                 std::string outputFileName = "superChill.txt");

bool isInterfacial(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                   const std::vector<std::vector<int>> &nList, int iatom,
                   int num_staggrd, int num_eclipsd, bool chillPlus = false);

[[nodiscard]] int numStaggered(molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                 const std::vector<std::vector<int>> &nList, int jatom);

struct SteinhardtQl {
  std::vector<double> ql;    
  std::vector<double> qlBar; 
};

[[nodiscard]] SteinhardtQl
steinhardtQl(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
             const std::vector<std::vector<int>> &nList, int orderL);
 
} // namespace chill

namespace sph {
 
// 7 is for Q3, orderL=3
 
std::vector<std::complex<double>>
spheriHarmo(int orderL, std::array<double, 2> radialCoord);
 
std::array<double, 2> radialCoord(std::array<double, 3> cartCoord);

std::vector<std::complex<double>>
lookupTableQ3Vec(std::array<double, 2> angles);

std::complex<double> lookupTableQ3(int m, std::array<double, 2> angles);

std::vector<std::complex<double>>
lookupTableQ4Vec(std::array<double, 2> angles);

std::complex<double> lookupTableQ4(int m, std::array<double, 2> angles);

std::vector<std::complex<double>>
lookupTableQ6Vec(std::array<double, 2> angles);

std::vector<std::complex<double>>
lookupTableQ8Vec(std::array<double, 2> angles);

std::complex<double> lookupTableQ8(int m, std::array<double, 2> angles);

std::complex<double> lookupTableQ6(int m, std::array<double, 2> angles);
 
} // namespace sph
 
#endif // SEAMS_BOP_H_