atomgraph_refinement.h

// -*-c++-*-
#ifndef STORMM_AG_REFINEMENT_H
#define STORMM_AG_REFINEMENT_H
 
#include <string>
#include <vector>
#include "copyright.h"
#include "Constants/behavior.h"
#include "DataTypes/stormm_vector_types.h"
#include "Math/series_ops.h"
#include "Math/vector_ops.h"
#include "Parsing/parse.h"
#include "Structure/structure_enumerators.h"
#include "atomgraph_enumerators.h"
 
namespace stormm {
namespace topology {
 
using constants::ExceptionResponse;
using parse::WildCardKind;
using stmath::minValue;
using stmath::incrementingSeries;
using structure::ApplyConstraints;

template <typename T> class ParameterUnion {
public:
  
  explicit ParameterUnion(const T* comp_xa, const T* comp_ya, const T* comp_za, const T* comp_wa,
                          const T* comp_va, int aparam_count, const T* comp_xb, const T* comp_yb,
                          const T* comp_zb, const T* comp_wb, const T* comp_vb, int bparam_count,
                          double match_tol = constants::small);
 
  explicit ParameterUnion(const T* comp_xa, const T* comp_ya, const T* comp_za,
                          int aparam_count_in, const T* comp_xb, const T* comp_yb,
                          const T* comp_zb, int bparam_count_in,
                          double match_tol = constants::small);
 
  explicit ParameterUnion(const T* comp_xa, const T* comp_ya, int aparam_count_in,
                          const T* comp_xb, const T* comp_yb, int bparam_count_in,
                          double match_tol = constants::small);
 
  explicit ParameterUnion(const T* comp_xa, int aparam_count_in, const T* comp_xb,
                          int bparam_count_in, double match_tol = constants::small);
 
  explicit ParameterUnion(const T* comp_xa, const T* comp_ya, const T* comp_za, const T* comp_wa,
                          const T* comp_va, int aparam_count_in);
 
  explicit ParameterUnion(const T* comp_xa, const T* comp_ya, const T* comp_za,
                          int aparam_count_in);
 
  explicit ParameterUnion(const T* comp_xa, const T* comp_ya, int aparam_count_in);
 
  explicit ParameterUnion(const T* comp_xa, int aparam_count_in);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya,
                          const std::vector<T> &comp_za, const std::vector<T> &comp_wa,
                          const std::vector<T> &comp_va, const std::vector<T> &comp_xb,
                          const std::vector<T> &comp_yb, const std::vector<T> &comp_zb,
                          const std::vector<T> &comp_wb, const std::vector<T> &comp_vb,
                          double match_tol = constants::small);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya,
                          const std::vector<T> &comp_za, const std::vector<T> &comp_xb,
                          const std::vector<T> &comp_yb, const std::vector<T> &comp_zb,
                          double match_tol = constants::small);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya,
                          const std::vector<T> &comp_xb, const std::vector<T> &comp_yb,
                          double match_tol = constants::small);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya,
                          double match_tol);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya,
                          const std::vector<T> &comp_za, const std::vector<T> &comp_wa,
                          const std::vector<T> &comp_va);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya,
                          const std::vector<T> &comp_za);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa, const std::vector<T> &comp_ya);
 
  explicit ParameterUnion(const std::vector<T> &comp_xa);

  ParameterUnion(const ParameterUnion &original) = default;
  ParameterUnion(ParameterUnion &&original) = default;
  ParameterUnion& operator=(const ParameterUnion &original) = default;
  ParameterUnion& operator=(ParameterUnion &&original) = default;

  int getInputSetParameterCount(int set_index) const;
  
  int getUniqueParameterCount() const;

  std::vector<T> getUnion(int comp_idx) const;

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

  void addSet(const T* comp_x, const T* comp_y, const T* comp_z, const T* comp_w, const T* comp_v,
              int nparm, double tol = constants::small);
 
  void addSet(const T* comp_x, const T* comp_y, const T* comp_z, int nparm,
              double tol = constants::small);
 
  void addSet(const T* comp_x, const T* comp_y, int nparm, double tol = constants::small);
 
  void addSet(const T* comp_x, int nparm, double tol = constants::small);
 
  void addSet(const std::vector<T> &comp_x, const std::vector<T> &comp_y,
              const std::vector<T> &comp_z, const std::vector<T> &comp_w,
              const std::vector<T> &comp_v, double tol = constants::small);
 
  void addSet(const std::vector<T> &comp_x, const std::vector<T> &comp_y,
              const std::vector<T> &comp_z, double tol = constants::small);
 
  void addSet(const std::vector<T> &comp_x, const std::vector<T> &comp_y,
              double tol = constants::small);
  
private:
  int total_parameters;              
  int component_count;               
  int set_count;                     
  std::vector<T> union_parameter_x;  
  std::vector<T> union_parameter_y;  
  std::vector<T> union_parameter_z;  
  std::vector<T> union_parameter_w;  
  std::vector<T> union_parameter_v;  
  
  std::vector<std::vector<int>> set_to_consensus_map;

  void validateSetIndex(int set_index, const char* caller = nullptr) const;
};

class CmapSurfaceUnion {
public:
  
  CmapSurfaceUnion();
 
  CmapSurfaceUnion(const double* surf_a, const int* dim_a, int nmap_a, const double* surf_b,
                   const int* dim_b, int nmap_b, double match_tol = constants::small);
 
  CmapSurfaceUnion(const double* surf_a, const int* dim_a, int nmap_a);
 
  CmapSurfaceUnion(const std::vector<double> &surf_a, const std::vector<int> &dim_a,
                   const std::vector<double> &surf_b, const std::vector<int> &dim_b,
                   double match_tol = constants::small);
 
  CmapSurfaceUnion(const std::vector<double> &surf_a, const std::vector<int> &dim_a);

  CmapSurfaceUnion(const CmapSurfaceUnion &orig) = default;
  CmapSurfaceUnion(CmapSurfaceUnion &&orig) = default;
  CmapSurfaceUnion& operator=(const CmapSurfaceUnion &orig) = default;
  CmapSurfaceUnion& operator=(CmapSurfaceUnion &&orig) = default;

  int getUniqueSurfaceCount() const;

  int getContributingTopologyCount() const;
  
  std::vector<double> getCmapSurface(int surf_index) const;

  const std::vector<double>& getAllSurfaces() const;

  const std::vector<int>& getSurfaceDimensions() const;
  int getSurfaceDimensions(int surf_index) const;

  const std::vector<int>& getCorrespondence(int set_index) const;
  int getCorrespondence(int set_index, int surf_index) const;
  
  void addSet(const double* surf_v, const int* dim_v, int nmap, double tol = constants::small);
  
private:
 
  int surface_count;                   
  int set_count;                       
  std::vector<double> surface_values;  
  std::vector<int> surface_dims;       
  std::vector<int> surface_bounds;     

  std::vector<std::vector<int>> set_to_consensus_map;

  void populateByOneSet(const double* surf_v, const int* dim_v, int nmap);

  void validateParameterIndex(int surf_index, const char* caller) const;

  void validateSetIndex(int set_index, const char* caller) const;
};
  
class BasicValenceTable {
public:

  BasicValenceTable();
  BasicValenceTable(int natom_in, int nbond_in, int nangl_in, int ndihe_in,
                    const std::vector<int> &bond_i_atoms_in = {},
                    const std::vector<int> &bond_j_atoms_in = {},
                    const std::vector<int> &bond_param_idx_in = {},
                    const std::vector<int> &angl_i_atoms_in = {},
                    const std::vector<int> &angl_j_atoms_in = {},
                    const std::vector<int> &angl_k_atoms_in = {},
                    const std::vector<int> &angl_param_idx_in = {},
                    const std::vector<int> &dihe_i_atoms_in = {},
                    const std::vector<int> &dihe_j_atoms_in = {},
                    const std::vector<int> &dihe_k_atoms_in = {},
                    const std::vector<int> &dihe_l_atoms_in = {},
                    const std::vector<int> &dihe_param_idx_in = {},
                    const std::vector<char4> &bond_mods_in = {},
                    const std::vector<char4> &angl_mods_in = {},
                    const std::vector<char4> &dihe_mods_in = {});

  void checkBondAngleRelevance(ApplyConstraints use_shake, ApplyConstraints use_settle,
                               const std::vector<int> &atomic_numbers);
  
  void makeAtomAssignments();
 
  // The member variables of this class are public.
  int total_bonds;                        
  int total_angls;                        
  int total_dihes;                        
  std::vector<int> bond_i_atoms;          
  std::vector<int> bond_j_atoms;          
  std::vector<int> bond_param_idx;        
  std::vector<int> angl_i_atoms;          
  std::vector<int> angl_j_atoms;          
  std::vector<int> angl_k_atoms;          
  std::vector<int> angl_param_idx;        
  std::vector<int> dihe_i_atoms;          
  std::vector<int> dihe_j_atoms;          
  std::vector<int> dihe_k_atoms;          
  std::vector<int> dihe_l_atoms;          
  std::vector<int> dihe_param_idx;        
  std::vector<char4> bond_mods;           
  std::vector<char4> angl_mods;           
  std::vector<char4> dihe_mods;           
  std::vector<int> bond_assigned_atoms;   
  std::vector<int> bond_assigned_index;   
  std::vector<int> bond_assigned_terms;   
  std::vector<int> bond_assigned_bounds;  
  std::vector<int> angl_assigned_atoms;   
  std::vector<int> angl_assigned_index;   
  std::vector<int> angl_assigned_terms;   
  std::vector<int> angl_assigned_bounds;  
  std::vector<int> dihe_assigned_atoms;   
  std::vector<int> dihe_assigned_index;   
  std::vector<int> dihe_assigned_terms;   
  std::vector<int> dihe_assigned_bounds;  
  std::vector<char4> bond_assigned_mods;  
  std::vector<char4> angl_assigned_mods;  
  std::vector<char4> dihe_assigned_mods;  
  std::vector<uint> bond_relevance;       
  std::vector<uint> angl_relevance;       
};

class CharmmValenceTable {
public:
  
  CharmmValenceTable();
  CharmmValenceTable(int natom_in, int nubrd_in, int nimpr_in, int ncmap_in,
                     const std::vector<int> &ubrd_i_atoms_in = {},
                     const std::vector<int> &ubrd_k_atoms_in = {},
                     const std::vector<int> &ubrd_param_idx_in = {},
                     const std::vector<int> &impr_i_atoms_in = {},
                     const std::vector<int> &impr_j_atoms_in = {},
                     const std::vector<int> &impr_k_atoms_in = {},
                     const std::vector<int> &impr_l_atoms_in = {},
                     const std::vector<int> &impr_param_idx_in = {},
                     const std::vector<int> &cmap_i_atoms_in = {},
                     const std::vector<int> &cmap_j_atoms_in = {},
                     const std::vector<int> &cmap_k_atoms_in = {},
                     const std::vector<int> &cmap_l_atoms_in = {},
                     const std::vector<int> &cmap_m_atoms_in = {},
                     const std::vector<int> &cmap_param_idx_in = {});
  
  void makeAtomAssignments();
 
  // The member variables of this class are public.
  int total_ub_angles;
  int total_impropers;
  int total_cmaps;
  std::vector<int> ubrd_i_atoms;
  std::vector<int> ubrd_k_atoms;
  std::vector<int> ubrd_param_idx;
  std::vector<int> impr_i_atoms;
  std::vector<int> impr_j_atoms;
  std::vector<int> impr_k_atoms;
  std::vector<int> impr_l_atoms;
  std::vector<int> impr_param_idx;
  std::vector<int> cmap_i_atoms;
  std::vector<int> cmap_j_atoms;
  std::vector<int> cmap_k_atoms;
  std::vector<int> cmap_l_atoms;
  std::vector<int> cmap_m_atoms;
  std::vector<int> cmap_param_idx;
  std::vector<int> ubrd_assigned_atoms;
  std::vector<int> ubrd_assigned_index;
  std::vector<int> ubrd_assigned_terms;
  std::vector<int> ubrd_assigned_bounds;
  std::vector<int> impr_assigned_atoms;
  std::vector<int> impr_assigned_index;
  std::vector<int> impr_assigned_terms;
  std::vector<int> impr_assigned_bounds;
  std::vector<int> cmap_assigned_atoms;
  std::vector<int> cmap_assigned_index;
  std::vector<int> cmap_assigned_terms;
  std::vector<int> cmap_assigned_bounds;
};

struct AttenuationParameterSet {

  AttenuationParameterSet(int set_count_in = 0);
  
  int total_14_sets;                          
  std::vector<int> dihe14_parameter_indices;  
  std::vector<double> elec_screening_factors; 
  std::vector<double> vdw_screening_factors;  
};

struct CondensedExclusions {

  CondensedExclusions();
  CondensedExclusions(int natom_in, int total_exclusions_in);
 
  int total_exclusions;
  std::vector<int> atom_excl_bounds;
  std::vector<int> atom_excl_list;
};

struct VirtualSiteTable {

  VirtualSiteTable();
 
  VirtualSiteTable(int natom_in, int vs_count_in);
 
  VirtualSiteTable(int natom_in, const std::vector<int> &vs_atoms_in,
                   const std::vector<int> &frame_types_in, const std::vector<int> &frame1_atoms_in,
                   const std::vector<int> &frame2_atoms_in,
                   const std::vector<int> &frame3_atoms_in,
                   const std::vector<int> &frame4_atoms_in,const std::vector<int> &param_idx_in,
                   const std::vector<double> &frame_dim1_in,
                   const std::vector<double> &frame_dim2_in,
                   const std::vector<double> &frame_dim3_in);
 
  int vs_count;                     
  std::vector<int> vs_numbers;      
  std::vector<int> vs_atoms;        
  std::vector<int> frame_types;     
  std::vector<int> frame1_atoms;    
  std::vector<int> frame2_atoms;    
  std::vector<int> frame3_atoms;    
  std::vector<int> frame4_atoms;    
  std::vector<int> param_idx;       
  std::vector<double> frame_dim1;   
  std::vector<double> frame_dim2;   
  std::vector<double> frame_dim3;   
};

class Map1234 {
public:
  
  Map1234();
  Map1234(int natom_in, int nb11_count_in, int nb12_count_in, int nb13_count_in,
          int nb14_count_in);

  void addExclusions(const std::vector<int3> &plus_excl,
                     ExceptionResponse policy = ExceptionResponse::DIE);
 
  void addExclusions(const std::vector<int2> &plus_excl,
                     ExceptionResponse policy = ExceptionResponse::DIE);

  void removeExclusions(const std::vector<int3> &minus_excl,
                        ExceptionResponse policy = ExceptionResponse::DIE);
  
  void removeExclusions(const std::vector<int2> &minus_excl,
                        ExceptionResponse policy = ExceptionResponse::DIE);
 
  std::vector<int> nb11_excl_bounds;  
  std::vector<int> nb11_excl_list;    
  std::vector<int> nb12_excl_bounds;  
  std::vector<int> nb12_excl_list;    
  std::vector<int> nb13_excl_bounds;  
  std::vector<int> nb13_excl_list;    
  std::vector<int> nb14_excl_bounds;  
  std::vector<int> nb14_excl_list;    
};

struct CmapAccessories {

  CmapAccessories();
  CmapAccessories(const std::vector<int> &cmap_dimensions_in);
  
  std::vector<double> phi_derivatives;      
  std::vector<double> psi_derivatives;      
  std::vector<double> phi_psi_derivatives;  
  std::vector<int> patch_matrix_bounds;     
  std::vector<double> patch_matrix_form;    
};

struct SettleParm {
  double mo;     
  double mh;     
  double moh;    
  double mormt;  
  double mhrmt;  
  double ra;     
  double rb;     
  double rc;     
};

struct ConstraintTable {

  ConstraintTable(const std::vector<int> &atomic_numbers, const std::vector<double> &atomic_masses,
                  const std::vector<int> &mol_limits, const std::vector<int> &mol_contents,
                  const std::vector<int> &mol_home, const BasicValenceTable &bvt,
                  const Map1234 &all_nb_excl, const std::vector<double> &bond_equilibria,
                  const std::vector<double> &angl_equilibria);
 
  int cnst_group_count;                   
  int settle_group_count;                 
  int cnst_parameter_count;               
  int settle_parameter_count;             
  std::vector<int> cnst_group_list;       
  std::vector<int> cnst_group_bounds;     
  std::vector<int> cnst_group_param_idx;  
  std::vector<int> settle_ox_atoms;       
  std::vector<int> settle_h1_atoms;       
  std::vector<int> settle_h2_atoms;       
  std::vector<int> settle_param_idx;      

  std::vector<double2> cnst_parameter_list;

  std::vector<int> cnst_parameter_bounds;
 
  // Settle parameter sets
  std::vector<SettleParm> settle_measurements;
};

struct AttenuatedPair {
  int atom_i;           
  int atom_j;           
  double lj_scaling;    
  double elec_scaling;  
};
  
void smoothCharges(std::vector<double> *q, std::vector<double> *tmp_charge_parameters,
                   std::vector<int> *tmp_charge_indices, int *charge_parameter_count,
                   double rounding_tol, double precision,
                   const std::string &filename = std::string(""));
  
void expandLennardJonesTables(std::vector<double> *lj_a_values, std::vector<double> *lj_b_values,
                              std::vector<double> *lj_c_values,
                              std::vector<double> *lj_14_a_values,
                              std::vector<double> *lj_14_b_values,
                              std::vector<double> *lj_14_c_values,
                              std::vector<double> *hb_a_values, std::vector<double> *hb_b_values,
                              std::vector<double> *hb_c_values, int n_lj_types,
                              const std::vector<int> &nb_param_index);
  
CondensedExclusions processExclusions(const std::vector<int> &raw_counts,
                                      const std::vector<int> &raw_exclusions,
                                      const std::string &file_name);

CondensedExclusions calculatePrmtopExclusions(const Map1234 nb_excl);

int countPrmtopZeroExclusions(const Map1234 nb_excl);

BasicValenceTable basicValenceIndexing(int atom_count, const std::vector<int> &tmp_bond_atoms_h,
                                       const std::vector<int> &tmp_bond_atoms_noh,
                                       const std::vector<int> &tmp_angl_atoms_h,
                                       const std::vector<int> &tmp_angl_atoms_noh,
                                       const std::vector<int> &tmp_dihe_atoms_h,
                                       const std::vector<int> &tmp_dihe_atoms_noh);

CharmmValenceTable charmmValenceIndexing(int atom_count, const std::vector<int> &tmp_ub_atoms,
                                         const std::vector<int> &tmp_charmm_impr_atoms,
                                         const std::vector<int> &tmp_cmap_atoms);

AttenuationParameterSet condenseScreeningFactors(const BasicValenceTable &bvt,
                                                 const std::vector<double> &dihe_elec_screenings,
                                                 const std::vector<double> &dihe_vdw_screenings,
                                                 const double default_elec14_screening,
                                                 const double default_vdw14_screening);

VirtualSiteTable listVirtualSites(int expected_vsite_count, const std::string &file_name,
                                  const std::vector<double> &tmp_masses,
                                  const BasicValenceTable &bvt,
                                  const std::vector<double> &tmp_bond_equilibria,
                                  const std::vector<double> &tmp_angl_equilibria,
                                  const std::vector<char4> &tmp_atom_types,
                                  const std::vector<char4> &tmp_atom_names,
                                  const std::vector<int> &vsite_custom_frames,
                                  const std::vector<double> &vsite_custom_details);

void accumulateExclusionBounds(std::vector<int> *current_bounds, const int atom_a,
                               const int atom_b, const std::vector<int> &vsite_child_bounds,
                               const std::vector<int> &vsite_child_list);

void markPairExclusion(std::vector<int> *excl_list, std::vector<int> *counters,
                       const std::vector<int> &bounds, const int atom_i, const int atom_j);

void markExclusions(std::vector<int> *excl_list, std::vector<int> *excl_counters,
                    const std::vector<int> &excl_bounds, const int atom_a, const int atom_b,
                    const std::vector<int> &vsite_child_bounds,
                    const std::vector<int> &vsite_child_list);

void cullDuplicateExclusions(std::vector<int> *excl_list, std::vector<int> *excl_bounds);

void cullPriorExclusions(std::vector<int> *excl_list, std::vector<int> *excl_bounds,
                         const std::vector<int> &prior_list,
                         const std::vector<int> &prior_bounds);

Map1234 mapExclusions(int atom_count, const std::vector<int> &vs_particles,
                      const std::vector<int> &vs_parent_atoms,
                      const std::vector<int> &bond_i_atoms, const std::vector<int> &bond_j_atoms);
  
Map1234 mapExclusions(int atom_count, const BasicValenceTable &bvt, const VirtualSiteTable &vst);

void exclusionSearchLoop(std::vector<int> &coverage, const std::vector<int> &bounds,
                         const std::vector<int> &excl_list, const CondensedExclusions &ce);

void checkExclusions(const CondensedExclusions &ce, const Map1234 &all_nb_excl,
                     const std::string &file_name);

std::vector<int> atomicNumbersFromMasses(const std::vector<double> &masses,
                                         const std::vector<char4> &atom_names,
                                         const std::vector<int> &nb12_list,
                                         const std::vector<int> &nb12_bounds,
                                         const std::string &source,
                                         const ExceptionResponse policy);

std::vector<int> traceBondedPatterns(const Map1234 &all_nb_excl);

void mapMolecules(const int atom_count, int *molecule_count, int *residue_count,
                  const Map1234 &all_nb_excl, std::vector<int> *molecule_membership,
                  std::vector<int> *molecule_limits, std::vector<int> *molecule_contents,
                  std::vector<int> *residue_limits, std::vector<char4> *residue_names);

std::vector<double> cubicSplineDerivativeStencil(const int npts);

CmapAccessories computeCmapDerivatives(int cmap_surf_count,
                                       const std::vector<int> tmp_cmap_surf_dims,
                                       const std::vector<int> tmp_cmap_surf_bounds,
                                       const std::vector<double> tmp_cmap_surfaces);

std::vector<int3> checkDihedral14Coverage(int atom_count, const std::vector<int> &atomic_numbers,
                                          const BasicValenceTable &bvt, const Map1234 &all_nb_excl,
                                          const VirtualSiteTable &vs_tbl,
                                          const AttenuationParameterSet &attn_parm,
                                          ExceptionResponse policy);

int reviewLargestResidue(const std::vector<int> residue_limits, int current_lr_size,
                         ExceptionResponse policy);

std::vector<int> matchExtendedName(const char4* overflow_names, int n_overflow,
                                   const std::string &query,
                                   const std::vector<WildCardKind> &wildcards = {});

SettleParm getSettleParameters(int ox_idx, int h1_idx, int h2_idx,
                               const std::vector<double> &atomic_masses,
                               const BasicValenceTable &bvt,
                               const std::vector<double> &bond_equilibria,
                               const std::vector<double> &angl_equilibria);
 
} // namespace topology
} // namespace stormm
 
#include "atomgraph_refinement.tpp"
 
#endif