ppitable.h

// -*-c++-*-
#ifndef STORMM_PPI_TABLE_H
#define STORMM_PPI_TABLE_H
 
#include <vector>
#include "copyright.h"
#include "Constants/behavior.h"
#include "Constants/symbol_values.h"
#include "DataTypes/common_types.h"
#include "Math/log_scale_spline.h"
#include "Math/math_enumerators.h"
#include "Parsing/parsing_enumerators.h"
#include "energy_enumerators.h"
#include "pme_util.h"
 
namespace stormm {
namespace energy {
  
using constants::ExceptionResponse;
using constants::PrecisionModel;
using symbols::amber_ancient_bioq;
using card::Hybrid;
using card::HybridKind;
using card::HybridTargetLevel;
using parse::NumberFormat;
using stmath::BasisFunctions;
using stmath::LogScaleSpline;
using stmath::LogSplineForm;
using stmath::LogSplineTable;
using stmath::TableIndexing;
using stmath::doublePrecisionSplineDetailMask;
using stmath::singlePrecisionSplineDetailMask;

template <typename T, typename T4> struct PPIKit {

  PPIKit(NonbondedTheme theme_in, BasisFunctions basis_in, TableIndexing lookup_in,
         int index_bound_in, int excl_offset_in, int index_shift_bits_in, ullint dp_detail_mask_in,
         uint sp_detail_mask_in, T arg_offset_in, const T4* energy_in, const T4* force_in,
         const T4* energy_excl_in, const T4* force_excl_in);

  PPIKit(const PPIKit &original) = default;
  PPIKit(PPIKit &&original) = default;
 
  const NonbondedTheme theme;   
  const BasisFunctions basis;   
  const TableIndexing lookup;   
  const int index_bound;        
  const int excl_offset;        
  const int index_shift_bits;   
  const ullint dp_detail_mask;  
  const uint sp_detail_mask;    
  const T arg_offset;           
  const T4* energy;             
  const T4* force;              
  const T4* energy_excl;        
  const T4* force_excl;         
};

template <typename T> struct PPIeKit {

  PPIeKit(NonbondedTheme theme_in, BasisFunctions basis_in, TableIndexing lookup_in,
          int index_bound_in, int excl_offset_in, int index_shift_bits_in,
          ullint dp_detail_mask_in, uint sp_detail_mask_in, T arg_offset_in, const T* energy_x_in,
          const T* energy_y_in, const T* energy_z_in, const T* energy_w_in, const T* force_x_in,
          const T* force_y_in, const T* force_z_in, const T* force_w_in, const T* energy_excl_x_in,
          const T* energy_excl_y_in, const T* energy_excl_z_in, const T* energy_excl_w_in,
          const T* force_excl_x_in, const T* force_excl_y_in, const T* force_excl_z_in,
          const T* force_excl_w_in);

  PPIeKit(const PPIeKit &original) = default;
  PPIeKit(PPIeKit &&original) = default;
 
  const NonbondedTheme theme;   
  const BasisFunctions basis;   
  const TableIndexing lookup;   
  const int index_bound;        
  const int excl_offset;        
  const int index_shift_bits;   
  const ullint dp_detail_mask;  
  const uint sp_detail_mask;    
  const T arg_offset;           
  const T* energy_x;            
  const T* energy_y;            
  const T* energy_z;            
  const T* energy_w;            
  const T* force_x;             
  const T* force_y;             
  const T* force_z;             
  const T* force_w;             
  const T* energy_excl_x;       
  const T* energy_excl_y;       
  const T* energy_excl_z;       
  const T* energy_excl_w;       
  const T* force_excl_x;        
  const T* force_excl_y;        
  const T* force_excl_z;        
  const T* force_excl_w;        
};
  
class PPITable {
public:

  PPITable(NonbondedTheme theme_in = NonbondedTheme::ELECTROSTATIC,
           BasisFunctions basis_set_in = BasisFunctions::POLYNOMIAL,
           TableIndexing indexing_method_in = TableIndexing::SQUARED_ARG,
           double cutoff_in = default_pme_cutoff, double argument_offset_in = 0.0,
           double dsum_tol_in = default_dsum_tol, int mantissa_bits_in = 5,
           double coulomb_in = amber_ancient_bioq, double min_range_in = 0.015625);
  
  template <typename T4>
  PPITable(const LogScaleSpline<T4> &spl_a, const LogScaleSpline<T4> &spl_b,
           const LogScaleSpline<T4> &spl_c, const LogScaleSpline<T4> &spl_d,
           double cutoff_in = default_pme_cutoff);
 
  template <typename T4>
  PPITable(const LogScaleSpline<T4> &spl_a, const LogScaleSpline<T4> &spl_b,
           const LogScaleSpline<T4> &spl_c);
 
  template <typename T4>
  PPITable(const LogScaleSpline<T4> &spl_a, const LogScaleSpline<T4> &spl_b);
 
  template <typename T4>
  PPITable(const LogScaleSpline<T4> &spl_a);
  
  PPITable(const PPITable &original);
  PPITable(PPITable &&original);
  PPITable& operator=(const PPITable &original);
  PPITable& operator=(PPITable &&original);

  NonbondedTheme getTheme() const;

  double getCutoff() const;

  double getMaximumRange() const;

  double getIndexingOffset() const;
  
  double getDirectSumTolerance() const;

  int getBitStride() const;

  double getEwaldCoefficient() const;

  double getGaussianWidth() const;

  int getTableIndex(double arg, PrecisionModel prec = PrecisionModel::SINGLE) const;

  int getTableIndexByRealArg(double arg, PrecisionModel prec = PrecisionModel::SINGLE) const;

  double evaluate(double arg, LogSplineForm kind, PrecisionModel prec = PrecisionModel::SINGLE,
                  bool use_elemental_tables = false) const;

  double evaluateByRealArg(double arg, LogSplineForm kind,
                           PrecisionModel prec = PrecisionModel::SINGLE,
                           bool use_elemental_tables = false) const;
  
  const PPIKit<double, double4> dpData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const PPIKit<float, float4> spData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const PPIeKit<double> dpeData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const PPIeKit<float> speData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;
 
#ifdef STORMM_USE_HPC
  void upload();
  
  void download();
#endif
 
private:
 
  NonbondedTheme theme;           
  BasisFunctions basis_set;       
  TableIndexing indexing_method;  
  double cutoff;                  
  double max_range;               
  double min_range;               
  double argument_offset;         
  double dsum_tol;                
  double ew_coeff;                
  int mantissa_bits;              
  double coulomb;                 
  int dp_exclusion_offset;        
  int sp_exclusion_offset;        
  ullint dp_detail_bitmask;       
  uint sp_detail_bitmask;         
 
  // The main arrays, in full (double) precision
  Hybrid<double4> energy;  
  Hybrid<double4> force;   

  Hybrid<double4> energy_with_exclusions;

  Hybrid<double4> force_with_exclusions;
 
  // The main arrays, in single precision
  Hybrid<float4> sp_energy;  
  Hybrid<float4> sp_force;   
  
  Hybrid<float4> sp_energy_with_exclusions;

  Hybrid<float4> sp_force_with_exclusions;

  Hybrid<double4> coeffs;

  Hybrid<float4> sp_coeffs;

  Hybrid<double> energy_x;
  Hybrid<double> energy_y;
  Hybrid<double> energy_z;
  Hybrid<double> energy_w;

  Hybrid<double> force_x;
  Hybrid<double> force_y;
  Hybrid<double> force_z;
  Hybrid<double> force_w;

  Hybrid<double> energy_with_excl_x;
  Hybrid<double> energy_with_excl_y;
  Hybrid<double> energy_with_excl_z;
  Hybrid<double> energy_with_excl_w;

  Hybrid<double> force_with_excl_x;
  Hybrid<double> force_with_excl_y;
  Hybrid<double> force_with_excl_z;
  Hybrid<double> force_with_excl_w;

  Hybrid<float> sp_energy_x;
  Hybrid<float> sp_energy_y;
  Hybrid<float> sp_energy_z;
  Hybrid<float> sp_energy_w;

  Hybrid<float> sp_force_x;
  Hybrid<float> sp_force_y;
  Hybrid<float> sp_force_z;
  Hybrid<float> sp_force_w;

  Hybrid<float> sp_energy_with_excl_x;
  Hybrid<float> sp_energy_with_excl_y;
  Hybrid<float> sp_energy_with_excl_z;
  Hybrid<float> sp_energy_with_excl_w;

  Hybrid<float> sp_force_with_excl_x;
  Hybrid<float> sp_force_with_excl_y;
  Hybrid<float> sp_force_with_excl_z;
  Hybrid<float> sp_force_with_excl_w;

  Hybrid<double> elemental_coeffs;
  
  Hybrid<float> sp_elemental_coeffs;
  
  template <typename T4> NonbondedTheme findTheme(const LogScaleSpline<T4> &spl) const;

  template <typename T4> std::vector<LogScaleSpline<T4>> buildAllSplineTables() const;
  
  template <typename T4>
  const LogScaleSpline<T4>& findNonExclPotential(const LogScaleSpline<T4> &spl_a,
                                                 const LogScaleSpline<T4> &spl_b,
                                                 const LogScaleSpline<T4> &spl_c,
                                                 const LogScaleSpline<T4> &spl_d) const;

  template <typename T4>
  const LogScaleSpline<T4>& findExclPotential(const LogScaleSpline<T4> &spl_a,
                                              const LogScaleSpline<T4> &spl_b,
                                              const LogScaleSpline<T4> &spl_c,
                                              const LogScaleSpline<T4> &spl_d) const;

  template <typename T4>
  const LogScaleSpline<T4>& findNonExclForce(const LogScaleSpline<T4> &spl_a,
                                             const LogScaleSpline<T4> &spl_b,
                                             const LogScaleSpline<T4> &spl_c,
                                             const LogScaleSpline<T4> &spl_d) const;

  template <typename T4>
  const LogScaleSpline<T4>& findExclForce(const LogScaleSpline<T4> &spl_a,
                                          const LogScaleSpline<T4> &spl_b,
                                          const LogScaleSpline<T4> &spl_c,
                                          const LogScaleSpline<T4> &spl_d) const;

  template <typename T4> uint checkPriority(const LogScaleSpline<T4> &spl_x) const;

  template <typename T4> void checkSplineCompatibility(const LogScaleSpline<T4> &spl_a,
                                                       const LogScaleSpline<T4> &spl_b) const;
  
  LogSplineForm findMissingForm(const uint holdings) const;
  
  template <typename T4>
  LogScaleSpline<T4> getTablePriority(const LogScaleSpline<T4> &spl_a,
                                      const LogScaleSpline<T4> &spl_b,
                                      const LogScaleSpline<T4> &spl_c) const;
 
  template <typename T4>
  LogScaleSpline<T4> getTablePriority(const LogScaleSpline<T4> &spl_a,
                                      const LogScaleSpline<T4> &spl_b) const;
 
  template <typename T4>
  LogScaleSpline<T4> getTablePriority(const LogScaleSpline<T4> &spl_a) const;

  template <typename T4>
  void populateCoefficients(const LogScaleSpline<T4> &u, const LogScaleSpline<T4> &du,
                            const LogScaleSpline<T4> &ux, const LogScaleSpline<T4> &dux);
};
 
} // namespace energy
} // namespace stormm
 
#include "ppitable.tpp"
 
#endif