atom_equivalence.h

// -*-c++-*-
#ifndef STORMM_ATOM_EQUIVALENCE_H
#define STORMM_ATOM_EQUIVALENCE_H
 
#include <string>
#include <vector>
#include "copyright.h"
#include "Topology/atomgraph.h"
#include "Trajectory/coordinateframe.h"
#include "UnitTesting/stopwatch.h"
#include "chemical_features.h"
#include "chemistry_enumerators.h"
 
namespace stormm {
namespace chemistry {
 
using topology::AtomGraph;
using testing::StopWatch;

class AtomRank {
public:

  AtomRank(const AtomGraph *ag_in = nullptr);
  
  AtomRank(const AtomGraph *ag_in, const std::vector<double> &formal_charges,
           const std::vector<double> &free_electrons, const std::vector<ullint> &ring_inclusion,
           const std::vector<ChiralOrientation> &chiralities, std::vector<int> *a_idx_tree,
           std::vector<int> *b_idx_tree, std::vector<int> *a_zn_tree, std::vector<int> *b_zn_tree,
           std::vector<double> *a_fc_tree, std::vector<double> *b_fc_tree,
           std::vector<double> *a_fe_tree, std::vector<double> *b_fe_tree,
           std::vector<ullint> *a_ri_tree, std::vector<ullint> *b_ri_tree,
           std::vector<ChiralOrientation> *a_ch_tree, std::vector<ChiralOrientation> *b_ch_tree,
           std::vector<int> *a_coverage, std::vector<int> *b_coverage, int low_molecule_index = 0,
           int high_molecule_index = -1);
 
  AtomRank(const AtomGraph &ag_in, const std::vector<double> &formal_charges,
           const std::vector<double> &free_electrons, const std::vector<ullint> &ring_inclusion,
           const std::vector<ChiralOrientation> &chiralities, int low_molecule_index = 0,
           int high_molecule_index = -1);
 
  AtomRank(const AtomGraph *ag_in, const ChemicalFeatures &chemfe, int low_molecule_index = 0,
           int high_molecule_index = -1);
 
  AtomRank(const AtomGraph &ag_in, const ChemicalFeatures &chemfe, int low_molecule_index = 0,
           int high_molecule_index = -1);

  const std::vector<int>& getRanks() const;
  std::vector<int> getRanks(int low_index, int high_index) const;
  int getRank(int atom_index) const;

  std::vector<int> getAtomsWithRank(const int rank_value) const;

  std::vector<int> getRankPartners(const int atom_index) const;
  
private:
  int maximum_rank;
  std::vector<int> ranks;
  std::vector<int> rank_degeneracy_bounds;
  std::vector<int> rank_instances;
  AtomGraph *ag_pointer;
};
  
class AtomEquivalence {
public:

  AtomEquivalence();
  
  AtomEquivalence(const AtomGraph *ag_in, const std::vector<double> &formal_charges,
                  const std::vector<double> &free_electrons,
                  const std::vector<ullint> &ring_inclusion,
                  const std::vector<ChiralOrientation> &chiralities, StopWatch *timer = nullptr,
                  int low_molecule_index = 0, int high_molecule_index = -1);
 
  AtomEquivalence(const AtomGraph &ag_in, const std::vector<double> &formal_charges,
                  const std::vector<double> &free_electrons,
                  const std::vector<ullint> &ring_inclusion,
                  const std::vector<ChiralOrientation> &chiralities, StopWatch *timer = nullptr,
                  int low_molecule_index = 0, int high_molecule_index = -1);
  
  AtomEquivalence(const ChemicalFeatures &chemfe_in, StopWatch *timer = nullptr,
                  int low_molecule_index = 0, int high_molecule_index = -1);
 
  AtomEquivalence(const AtomGraph &ag_in, const CoordinateFrame &cf, StopWatch *timer = nullptr,
                  int low_molecule_index = 0, int high_molecule_index = -1);

  AtomEquivalence(const AtomEquivalence &original) = default;
  AtomEquivalence(AtomEquivalence &&original) = default;
  AtomEquivalence& operator=(const AtomEquivalence &original) = default;
  AtomEquivalence& operator=(AtomEquivalence &&original) = default;
  
  int getGroupCount() const;

  int getSymmetryDepth() const;
  
  std::vector<int> getGroup(int group_index) const;

  int getSymmetryRelatedAtom(int group_index, int domain_index, int atom_index) const;

  const int* getGroupPointer(int group_index) const;

  int getGroupSize(int group_index) const;

  int getGroupOrder(int group_index) const;

  int getGroupLevel(int group_index) const;
  
  EquivalenceSwap getGroupRule(int group_index) const;
  
  std::vector<int> getGroupDependencies(int group_index) const;

  int getAsymmetricAtomCount() const;
  
  const std::vector<int>& getAsymmetricAtoms() const;

  const AtomGraph* getTopologyPointer() const;
 
private:
  int group_count;                     
  int symmetry_depth;                  
  std::vector<int> group_atoms;        
  std::vector<int> group_sizes;        
  std::vector<int> group_orders;       
  std::vector<int> group_bounds;       
  std::vector<int> group_levels;       
  std::vector<int> dependencies;       
  std::vector<int> dependency_bounds;  
  std::vector<int> asymmetric_atoms;   

  std::vector<EquivalenceSwap> group_rules;

  AtomGraph *ag_pointer;

  void findEquivalentAtoms(const std::vector<int> &subset_atoms, std::vector<int> *map_to_subset,
                           const std::vector<double> &formal_charges,
                           const std::vector<double> &free_electrons,
                           const std::vector<ullint> &ring_inclusion,
                           const std::vector<ChiralOrientation> &chiralities,
                           const AtomRank &arnks, std::vector<int> *domain_coverage,
                           std::vector<int> *allowed_atoms, std::vector<int> *candidate_hopper,
                           std::vector<int> *domain_assignment, std::vector<int> *jumbled_groups,
                           std::vector<int> *aligned_groups, std::vector<int> *layer_bounds);

  void drawEquivalentGroups(const std::vector<int> &partners, const std::vector<int> &subset_atoms,
                            const std::vector<double> &formal_charges,
                            const std::vector<double> &free_electrons,
                            const std::vector<ullint> &ring_inclusion,
                            const std::vector<ChiralOrientation> &chiralities,
                            const AtomRank &arnks, std::vector<int> *domain_coverage,
                            std::vector<int> *allowed_atoms, std::vector<int> *candidate_hopper,
                            std::vector<int> *domain_assignment, std::vector<int> *jumbled_groups,
                            std::vector<int> *aligned_groups, std::vector<int> *layer_bounds);

  void addSymmetryGroup(const std::vector<int> &all_domains, int domain_size,
                        const std::vector<int> &selected_domains, EquivalenceSwap plan);
  
  void freeForAllSymmetryGroup(const std::vector<int> &all_domains, int domain_size,
                               const std::vector<int> &selected_domains);

  void rotarySymmetryGroup(const std::vector<int> &all_domains, int domain_size,
                           const std::vector<int> &selected_domains,
                           const std::vector<bool> &touch_table);

  void cullDuplicateGroups();
  
  void findDependencies();

  void orderAllGroups();
};

std::vector<int> rankAtoms(const AtomGraph *ag, const std::vector<double> &formal_charges,
                           const std::vector<double> &free_electrons,
                           const std::vector<ullint> &ring_inclusion,
                           const std::vector<ChiralOrientation> &chiralities,
                           std::vector<int> *a_idx_tree, std::vector<int> *b_idx_tree,
                           std::vector<int> *a_zn_tree, std::vector<int> *b_zn_tree,
                           std::vector<double> *a_fc_tree, std::vector<double> *b_fc_tree,
                           std::vector<double> *a_fe_tree, std::vector<double> *b_fe_tree,
                           std::vector<ullint> *a_ri_tree, std::vector<ullint> *b_ri_tree,
                           std::vector<ChiralOrientation> *a_ch_tree,
                           std::vector<ChiralOrientation> *b_ch_tree,
                           std::vector<int> *a_coverage, std::vector<int> *b_coverage,
                           int low_molecule_index = 0, int high_molecule_index = -1);
 
std::vector<int> rankAtoms(const AtomGraph &ag, const std::vector<double> &formal_charges,
                           const std::vector<double> &free_electrons,
                           const std::vector<ullint> &ring_inclusion,
                           const std::vector<ChiralOrientation> &chiralities,
                           int low_molecule_index = 0, int high_molecule_index = -1);
 
std::vector<int> rankAtoms(const AtomGraph *ag, const ChemicalFeatures &chemfe,
                           int low_molecule_index = 0, int high_molecule_index = -1);
 
std::vector<int> rankAtoms(const AtomGraph &ag, const ChemicalFeatures &chemfe,
                           int low_molecule_index = 0, int high_molecule_index = -1);
 
} // namespace chemistry
} // namespace stormm
 
#endif