layered_potential.h

// -*-c++-*-
#ifndef STORMM_LAYERED_POTENTIAL_H
#define STORMM_LAYERED_POTENTIAL_H
 
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "Math/formulas.h"
#include "Math/log_scale_spline.h"
#include "Math/matrix_ops.h"
#include "Potential/pme_util.h"
#include "Structure/structure_enumerators.h"
#include "energy_enumerators.h"
#include "layered_potential_metrics.h"
 
namespace stormm {
namespace energy {
 
using card::Hybrid;
using card::HybridTargetLevel;
using energy::elecPMEDirectSpace;
using stmath::sigmoid;
using stmath::sigmoidf;
using stmath::LogScaleSpline;
using stmath::LogSplineTable;
using stmath::qrSolver;

template <typename T4>
class NrgLayerKit {

  NrgLayerKit(const LogSplineTable<T4> &u_in, const LogSplineTable<T4> &du_in,
              const LogSplineTable<T4> &d2u_in, const LogSplineTable<T4> &d3u_in,
              double cutoff_in);

  NrgLayerKit(const NrgLayerKit &original) = default;
  NrgLayerKit(NrgLayerKit &&original) = default;
  NrgLayerKit& operator=(const NrgLayerKit &original) = default;
  NrgLayerKit& operator=(NrgLayerKit &&original) = default;
 
  // A total of four logarithmic spline tables are contained within the abstract.
  const LogSplineTable<T4> u;    
  const LogSplineTable<T4> du;   
  const LogSplineTable<T4> d2u;  
  const LogSplineTable<T4> d3u;  
  const double cutoff;           
};
  
template <typename T, typename T4> class LayeredPotential {
public:

  LayeredPotential(const LayeredPotentialMetrics &parameters_in = LayeredPotentialMetrics());

  LayeredPotential(const LayeredPotential &original) = default;
  LayeredPotential(LayeredPotential &&original) = default;
  LayeredPotential& operator=(const LayeredPotential &other) = default;
  LayeredPotential& operator=(LayeredPotential &&other) = default;

  const LayeredPotentialMetrics& getParameters() const;

  double4 getSmoothingCoefficients(int layer) const;
  
  T getSplinedValue(int layer, T r, T r2 = -1.0) const;

  T getAnalyticValue(int layer, T r, T r2 = -1.0) const;

  T getSplinedDerivative(int layer, T r, T r2 = -1.0, int order = 1) const;

  T getAnalyticDerivative(int layer, T r, T r2 = -1.0, int order = 1) const;

  const NrgLayerKit<T4> data(int layer_index,
                             HybridTargetLevel tier = HybridTargetLevel::HOST) const;
  
private:
  LayeredPotentialMetrics parameters;
  std::vector<LogScaleSpline<T4>> spline_tables;
  std::vector<double4> smoothing_coefficients;

  void validateLayerIndex(int layer_index, const char* caller = nullptr) const;
  
  void solveLayerCoefficients();
};
 
} // namespace energy
} // namespace stormm
 
#include "layered_potential.tpp"
 
#endif