layered_potential_metrics.h

// -*-c++-*-
#ifndef STORMM_LAYERED_POTENTIAL_METRICS_H
#define STORMM_LAYERED_POTENTIAL_METRICS_H
 
#include <vector>
#include "copyright.h"
#include "Constants/scaling.h"
#include "Constants/symbol_values.h"
#include "Math/log_scale_spline.h"
#include "Structure/structure_enumerators.h"
#include "energy_enumerators.h"
 
namespace stormm {
namespace energy {
 
using structure::BoundaryCondition;
using symbols::amber_ancient_bioq;
  
constexpr double hail_minimum_range_limit = 4.0;
constexpr double hail_minimum_range_multiplier = 1.5;
constexpr double hail_maximum_range_multiplier = 3.0;
constexpr double hail_default_vdw_xtn_midpoint = 7.0;
constexpr double hail_default_vdw_xtn_intensity = 3.0;
constexpr double hail_default_ewald_coefficient = 0.3;
  
class LayeredPotentialMetrics {
public:

  LayeredPotentialMetrics(DecomposablePotential form_in = DecomposablePotential::ELECTROSTATIC,
                          BoundaryCondition boundaries_in = BoundaryCondition::ISOLATED,
                          double coulomb_in = amber_ancient_bioq,
                          VdwCombiningRule mixing_rule_in = VdwCombiningRule::GEOMETRIC,
                          const std::vector<double> &exponent_c_in = { 1.0, 1.5, 2.0 },
                          double particle_particle_cutoff = 5.0,
                          double range_multiplier_in = 2.0,
                          double range_limit_in = 1000.0,
                          double vdw_transition_midpoint_in = hail_default_vdw_xtn_midpoint,
                          double vdw_transition_intensity_in = hail_default_vdw_xtn_intensity,
                          double ewald_coefficient_in = hail_default_ewald_coefficient,
                          const std::vector<double> &box_vectors_in = { 50.0,  0.0,  0.0,
                                                                         0.0, 50.0,  0.0,
                                                                         0.0,  0.0, 50.0 });

  LayeredPotentialMetrics(const LayeredPotentialMetrics &original) = default;
  LayeredPotentialMetrics(LayeredPotentialMetrics &&original) = default;
  LayeredPotentialMetrics& operator=(const LayeredPotentialMetrics &other) = default;
  LayeredPotentialMetrics& operator=(LayeredPotentialMetrics &&other) = default;
  
  DecomposablePotential getForm() const;

  BoundaryCondition getBoundaryCondition() const;
  
  int getLayerCount() const;
  
  double getCutoff(int layer_index) const;

  double getCoulombConstant() const;

  VdwCombiningRule getMixingRule() const;

  double getExponentFactor(int factor_index, int layer_index = 1) const;

  double getRangeCompounding() const;

  double getMaximumRange() const;

  double getVdwTransitionMidpoint() const;

  double getVdwTransitionIntensity() const;

  double getEwaldCoefficient() const;
  
  const std::vector<double>& getUnitCellVectors() const;

  void setForm(DecomposablePotential form_in);

  void setBoundaryCondition(BoundaryCondition boundary_in);

  void setCutoff(double cutoff_in);

  void setCoulombConstant(double coulomb_in);

  void setMixingRule(VdwCombiningRule mixing_rule_in);

  void setExponentFactor(int factor_index, double exp_c_in);
  void setExponentFactor(const std::vector<double> &exp_c_in);
  void setExponentFactor(const double3 exp_c_in);

  void setRangeCompounding(double range_multiplier_in);

  void setMaximumRange(double range_limit_in);
  
  void setVdwTransitionMidpoint(double midpoint_in);

  void setVdwTransitionIntensity(double intensity_in);

  void setEwaldCoefficient(double ewald_coefficient_in);
  
  void setUnitCellVectors(const std::vector<double> &box_vectors_in);
 
private:
 
  DecomposablePotential form;       
  BoundaryCondition boundaries;     
  int layer_count;                  
  std::vector<double> cutoffs;      
  double coulomb_constant;          
  VdwCombiningRule mixing_rule;     
  std::vector<double3> exponent_c;  
  double range_multiplier;          
  double range_limit;               
  double vdw_transition_midpoint;   
  double vdw_transition_intensity;  
  double ewald_coefficient;         
  std::vector<double> box_vectors;  

  void validateLayerIndex(int layer_index, const char* caller = nullptr) const;

  void computeLayers();
};
 
} // namespace energy
} // namespace stormm
 
#endif