indigo.h

// -*-c++-*-
#ifndef STORMM_INDIGO_H
#define STORMM_INDIGO_H
 
#include <map>
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "DataTypes/common_types.h"
#include "DataTypes/mixed_types.h"
#include "Topology/atomgraph.h"
 
namespace stormm {
namespace chemistry {
 
using card::Hybrid;
using card::HybridTargetLevel;
using topology::AtomGraph;

constexpr int minimum_formal_charge = -8;
constexpr int maximum_formal_charge =  7;
constexpr int formal_charge_range   = maximum_formal_charge - minimum_formal_charge + 1;
constexpr int formal_charge_bits    = 4;

constexpr int default_indigo_energy_gap = 3187;

constexpr int maximum_negative_carbons = 1;
  
constexpr int maximum_bond_order = 3;
constexpr int bond_order_range   = maximum_bond_order + 1;
  
constexpr int maximum_indigo_atomic_number = 35;
constexpr int indigo_atomic_number_range   = maximum_indigo_atomic_number + 1;
  
int IndigoFormalChargeKey(int atomic_number, int formal_charge);

int IndigoBondOrderKey(int atomic_number_i, int atomic_number_j, int bond_order);

std::vector<int> indigoFormalChargeScores();
  
std::vector<int> indigoBondOrderScores();

class IndigoAtomCenter {
public:

  IndigoAtomCenter(int table_index_in, int z_number_in, int bond_count_in,
                   int valence_electrons_in, const std::vector<int> &relevant_bonds_in,
                   const std::vector<int> &partner_atoms_in, const std::vector<int> &atom_scores,
                   const std::vector<int> &bond_scores);

  int getIndigoTableIndex() const;
  
  int getStateCount() const;

  uint getState(int state_index) const;

  int getCharge(int state_index) const;

  int getAtomicNumber() const;
  
  int getScore(int state_index) const;

  void sortStates();
  
  int getPartnerCount() const;

  int getRelevantBond(int index) const;

  int getPartner(int index) const;

  int findPartnerIndex(int partner_atom_number) const;

  int getBondOrderOfState(int bond_index, int state_index) const;  
  
  bool queryBondToPartner(int partner_atom_number, int partner_bo) const;

  int cullByPartners(const std::vector<IndigoAtomCenter> &acen);
  
private:
  int table_index;                  
  int z_number;                     
  int bond_count;                   
  int valence_electrons;            
  int possible_states;              
  std::vector<int> relevant_bonds;  
  std::vector<int> partner_atoms;   
  std::vector<uint> states;         
  std::vector<int> scores;          

  uint computeState(const int fc_value, const std::vector<uint> &bond_orders);

  int computeScore(const uint state, const std::vector<int> &atom_scores,
                   const std::vector<int> &bond_scores);
};

class IndigoFragment {
public:

  IndigoFragment(const std::vector<int> &centers_list_in,
                 const std::vector<IndigoAtomCenter> &all_centers,
                 int score_delta = default_indigo_energy_gap);

  int getCenterCount() const;
  
  int getStateCount() const;
  
  std::vector<int2> getState(int state_index) const;

  int getGlobalCenter(int center_index) const;
  
  int getCharge(int state_index) const;

  int getScore(int state_index) const;

  std::vector<int> getChargeStates() const;

  std::vector<int3> getChargesAndBestEnergies() const;

  std::vector<int2> getStatesBearingCharge(int charge_value) const;

  bool testEquivalence(const IndigoFragment &other, const std::vector<int> &real_atom_map,
                       const AtomGraph *ag_pointer,
                       const std::vector<IndigoAtomCenter> &atom_centers) const;

  int cullStatesBearingCharge(int charge_value);  

  int cullHigherEnergyStatesByCharge(int charge_value,
                                     int score_delta = default_indigo_energy_gap);
  
private:
  int center_count;               
  int possible_states;            
  std::vector<int> centers_list;  
  std::vector<int*> states;       
  std::vector<int> states_data;   
  std::vector<int> net_charges;   
  std::vector<int> scores;        
};
  
class IndigoTable {
public:

  IndigoTable(const AtomGraph *ag_in, int molecule_index = 0, double temperature_in = 300.0);

  int getAtomCount() const;

  int getBondCount() const;

  int getNetCharge() const;

  std::vector<CombineIDp> getGroundStateFormalCharges() const;

  std::vector<CombineIDp> getGroundStateBondOrders() const;

  std::vector<CombineIDp> getGroundStateFreeElectrons() const;

  std::vector<int2> getZeroKelvinFormalCharges() const;

  std::vector<int2> getZeroKelvinBondOrders() const;

  std::vector<int2> getZeroKelvinFreeElectrons() const;
  
  double getGroundStateEnergy() const;  
  
private:
  int atom_count;                       
  int bond_count;                       
  int net_charge;                       
  int fragment_count;                   
  double temperature;                   
  double ground_state_energy;           
  std::vector<int> real_atom_map;       
  std::vector<int> real_bond_map;       
  std::vector<int> bond_i_atoms;        
  std::vector<int> bond_j_atoms;        
  std::vector<int> atom_connect_bounds; 
  std::vector<int> atom_relevant_bonds; 
  std::vector<int> atom_bond_partners;  
  std::vector<int> atom_scores;         
  std::vector<int> bond_scores;         
  std::vector<int> valence_electrons;   

  std::vector<IndigoAtomCenter> atom_centers;

  std::vector<IndigoFragment> mutable_fragments;

  std::vector<double> ground_state_formal_charges;

  std::vector<double> ground_state_bond_orders;

  std::vector<double> ground_state_free_electrons;

  std::vector<double> zerok_formal_charges;

  std::vector<double> zerok_bond_orders;

  std::vector<double> zerok_free_electrons;

  const AtomGraph *ag_pointer;

  double addToGroundState(const IndigoAtomCenter &ac, int state_index,
                          std::vector<double> *acc_formal_charges,
                          std::vector<double> *acc_bond_orders,
                          std::vector<double> *acc_free_electrons, int accumulation = 0,
                          double probability = 1.0);
};
 
} // namespace chemistry
} // namespace stormm
 
#endif