site.hpp

Go to the documentation of this file.

//-----------------------------------------------------------------------------------
// d-SEAMS - Deferred Structural Elucidation Analysis for Molecular Simulations
// SPDX-License-Identifier: MIT
//-----------------------------------------------------------------------------------
 
#ifndef SEAMS_SITE_H_
#define SEAMS_SITE_H_
 
#include <mol_sys.hpp>
 
#include <array>
#include <string_view>
#include <unordered_map>
#include <vector>

 
namespace site {
 
enum class Kind {
  unspecified,
  cationHead,
  anion,
  tail,
  donorH,
  acceptor,
  polar,
  apolar,
  waterO,
  waterH,
  solvent
};
 
enum class Family {
  waterIce,      // default; CHILL/TUM allowed
  ionicLiquid,   // CHILL/TUM refused
  moltenSalt,    // CHILL/TUM refused
  des,           // CHILL/TUM refused
  electrolyte,   // CHILL/TUM only if the caller also names a waterIce subset
  confinedIL,    // CHILL/TUM refused
  confinedWater, // 2D RDF / monolayer rings; bulk CHILL refused
  networkFormer  // silica / BeF2; Franzblau yes, CHILL no
};
 
struct Table {
  Family family = Family::waterIce;
  std::unordered_map<int, Kind> typeToKind;   // LAMMPS type ID
  std::unordered_map<int, Kind> atomOverride; // atom ID wins over type
  Kind of(const molSys::Point<double> &p) const;
  Kind ofType(int typeId) const;
};
 
std::vector<int>
indicesOf(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
          const Table &table, Kind kind);
 
int lammpsTypeOfKind(const Table &table, Kind kind); // error if not unique
 
// One vertex per molID that carries ionKind (cationHead or anion).
// Coordinates: unweighted geometric COM of atoms of that molID whose
// kind is ionKind, unwrapped with gen::relDist to the first atom of
// the group (createMolIDAtomIDMultiMap, mol_sys.hpp:187-190). Copies
// box / boxLow from src. Point::type is 1 for cationHead molecules
// and 2 for anion molecules so neighList(merged, 1, 2) is legal.
// A molecule that is already one tagged site (UA, or one designated
// type per ion) is a no-op: the COM is that site.
molSys::PointCloud<molSys::Point<double>, double>
ionCloud(const molSys::PointCloud<molSys::Point<double>, double> &src,
         const Table &table);
 
Table parseSiteSpec(std::string_view spec);
 

enum class IonState { liquid = 0, front = 1, ice = 2 };
struct IonEnvironment {
  std::vector<int> ion;            
  std::vector<int> shell;          
  std::vector<double> iceFraction; 
  std::vector<IonState> state;
  std::vector<std::vector<int>> members;  
  int nIce = 0;
  int nFront = 0;
  int nLiquid = 0;
};

IonEnvironment
ionEnvironment(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
               const std::vector<bool> &iceFlag, const std::vector<int> &ionIndices,
               int waterType, double cutoff);

std::vector<int> shellRingCensus(const std::vector<std::vector<int>> &rings,
                                 const std::vector<int> &shell, int maxRingSize);

struct GuestOccupancy {
  std::vector<int> guestsPerCage;        
  std::vector<int> cageOfGuest;          
  std::vector<double> centreDistance;    
  int occupied = 0;                      
  int multiply = 0;                      
  int free = 0;                          
  std::vector<int> occupancyHistogram;   
};

GuestOccupancy
guestOccupancy(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
               const std::vector<std::vector<int>> &cages,
               const std::vector<int> &guestIndices, double radius);

std::array<double, 3>
periodicCentroid(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                 const std::vector<int> &atoms);

GuestOccupancy
guestOccupancyInside(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                     const std::vector<std::vector<int>> &rings,
                     const std::vector<std::vector<int>> &cageFaces,
                     const std::vector<int> &guestIndices);

struct IceClusterIons {
  std::vector<int> clusterOf;      
  std::vector<int> ionsInCluster;  
  std::vector<int> clusterOfIon;   
  int nClusters = 0;
};
 
IceClusterIons
iceClusterIonCensus(const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
                    const std::vector<bool> &ice,
                    const std::vector<std::vector<int>> &nListByIndex,
                    const std::vector<int> &ionIndices, double cutoff);
} // namespace site
 
#endif // SEAMS_SITE_H_