lennard_jones_analysis.h

// -*-c++-*-
#ifndef LENNARD_JONES_ANALYSIS_H
#define LENNARD_JONES_ANALYSIS_H
 
#include <map>
#include <vector>
#include "copyright.h"
#include "Constants/behavior.h"
#include "DataTypes/stormm_vector_types.h"
#include "Potential/energy_enumerators.h"
#include "Synthesis/atomgraph_synthesis.h"
 
namespace stormm {
namespace topology {
 
using constants::ExceptionResponse;
using energy::VdwCombiningRule;
using synthesis::AtomGraphSynthesis;
  
struct PairLJInteraction {
public:

  PairLJInteraction(const char4 type_a_in = { ' ', ' ', ' ', ' ' },
                    const char4 type_b_in = { ' ', ' ', ' ', ' ' }, double lja_in = 0.0,
                    double ljb_in = 0.0, double lja_14_in = 0.0, double ljb_14_in = 0.0);

  PairLJInteraction(const PairLJInteraction &original) = default;
  PairLJInteraction(PairLJInteraction &&original) = default;
  PairLJInteraction& operator=(const PairLJInteraction &original) = default;
  PairLJInteraction& operator=(PairLJInteraction &&original) = default;
 
  // All member variables are public.
  char4 type_a;   
  char4 type_b;   
  double lja;     
  double ljb;     
  double lja_14;  
  double ljb_14;  
};
  
class LennardJonesAnalysis {
public:

  LennardJonesAnalysis(const NonbondedKit<double> &nbk,
                       const std::vector<std::vector<char4>> &atom_type_aliases_in);

  LennardJonesAnalysis(const LennardJonesAnalysis &original) = default;
  LennardJonesAnalysis(LennardJonesAnalysis &&original) = default;
  LennardJonesAnalysis& operator=(const LennardJonesAnalysis &original) = default;
  LennardJonesAnalysis& operator=(LennardJonesAnalysis &&original) = default;
  
  int getLJTypeCount() const;

  int getAtomTypeCount() const;

  int getSetCount() const;
  
  VdwCombiningRule getMostPrevalentCombiningRule() const;

  const std::vector<char4>& getLJAliases(int consensus_index) const;
 
  const std::vector<char4>& getLJAliases(double sigma_query, double epsilon_query,
                                         double tolerance = 1.0e-4) const;
 
  const std::vector<char4>& getLJAliases(const char4 atom_type_query,
                                         ExceptionResponse policy = ExceptionResponse::DIE) const;

  double getLJSigma(int consensus_index) const;

  double getLJ14Sigma(int consensus_index) const;

  double getLJEpsilon(int consensus_index) const;

  double getLJ14Epsilon(int consensus_index) const;

  double2 getLJParameters(int consensus_index) const;
  double2 getLJParameters(const char4 atom_type_query,
                          ExceptionResponse policy = ExceptionResponse::DIE) const;

  double2 getLJ14Parameters(int consensus_index) const;
  double2 getLJ14Parameters(const char4 atom_type_query,
                            ExceptionResponse policy = ExceptionResponse::DIE) const;

  double3 getLJCoefficients(int index_i, int index_j) const;
  double3 getLJCoefficients(const char4 atype_i, const char4 atype_j,
                            ExceptionResponse policy = ExceptionResponse::DIE) const;

  double3 getLJ14Coefficients(int index_i, int index_j) const;
  double3 getLJ14Coefficients(const char4 atype_i, const char4 atype_j,
                              ExceptionResponse policy = ExceptionResponse::DIE) const;

  int getCorrespondence(int set_index, int type_index) const;
  
  const std::vector<int>& getSetCorrespondence(int set_index) const;

  const std::vector<int2>& getInputInstances(int consensus_index) const;
  
  void addSet(const NonbondedKit<double> &nbk,
              const std::vector<std::vector<char4>> &othr_type_aliases);
 
private:

  int lj_type_count;

  int atom_type_count;

  int set_count;
  
  VdwCombiningRule prevalent_rule;

  VdwCombiningRule absolute_rule;
  
  std::vector<double> sigma;

  std::vector<double> sigma_14;

  std::vector<double> epsilon;

  std::vector<double> epsilon_14;

  std::vector<double> lja_coeff;
  std::vector<double> ljb_coeff;
  std::vector<double> ljc_coeff;
  std::vector<double> lja_14_coeff;
  std::vector<double> ljb_14_coeff;
  std::vector<double> ljc_14_coeff;

  std::vector<PairLJInteraction> edits;

  std::vector<int> ag_index_origins_bounds;

  std::vector<std::vector<char4>> atom_type_aliases;

  std::map<uint, std::vector<int>> atom_type_map;
  
  std::vector<std::vector<int>> set_to_consensus_map;

  std::vector<std::vector<int2>> consensus_to_set_map;

  void validateConsensusIndexQuery(int ljt_query) const;

  void confirmDistinctType(const char4 atom_type_query,
                           const ExceptionResponse policy = ExceptionResponse::DIE,
                           const char* caller = nullptr) const;
};

std::string listAtomTypesAsString(const std::vector<char4> &atyp_list);

int inferLennardJonesTypeCount(const int length_a, const int length_b,
                               const char* caller = nullptr);
  
template <typename T>
VdwCombiningRule inferCombiningRule(const T* lja, const T* ljb, int lj_type_count,
                                    ExceptionResponse policy = ExceptionResponse::WARN,
                                    bool seek_prevalent = false);
 
template <typename T>
VdwCombiningRule inferCombiningRule(const std::vector<T> &lja, const std::vector<T> &ljb,
                                    ExceptionResponse policy = ExceptionResponse::WARN,
                                    bool seek_prevalent = false);
 
template <typename T>
VdwCombiningRule inferCombiningRule(const Hybrid<T> &lja, const Hybrid<T> &ljb,
                                    ExceptionResponse policy = ExceptionResponse::WARN,
                                    bool seek_prevalent = false);

template <typename T>
std::vector<bool> findPairSpecificParticipation(const T* lja, const T* ljb, int lj_type_count);
 
template <typename T>
VdwCombiningRule findPairSpecificParticipation(const std::vector<T> &lja,
                                               const std::vector<T> &ljb);
 
template <typename T>
VdwCombiningRule findPairSpecificParticipation(const Hybrid<T> &lja, const Hybrid<T> &ljb);
  
} // namespace topology
} // namespace stormm
 
#include "lennard_jones_analysis.tpp"
 
#endif