mesh_forcefield.h

// -*-c++-*-
#ifndef STORMM_MESH_FORCEFIELD_H
#define STORMM_MESH_FORCEFIELD_H
 
#ifdef STORMM_USE_HPC
#  ifdef STORMM_USE_CUDA
#    include <cuda.h>
#    include <cuda_runtime.h>
#  endif
#endif
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "Constants/symbol_values.h"
#include "Math/math_enumerators.h"
#include "Potential/energy_enumerators.h"
#include "Potential/soft_core_potentials.h"
#include "Reporting/error_format.h"
#include "Topology/atomgraph.h"
 
namespace stormm {
namespace structure {
 
using card::Hybrid;
using card::HybridKind;
using card::HybridTargetLevel;
using energy::cubicSoftCore;
using energy::quinticSoftCore;
using energy::VdwCombiningRule;
using stmath::Interpolant;
using symbols::amber_ancient_bioq;
using topology::AtomGraph;

template <typename T> struct MeshFFKit {

  MeshFFKit();
  MeshFFKit(VdwCombiningRule lj_rule_in, double coulomb_in, double clash_ratio_in,
            double clash_distance_in, const T* probe_lja_in, const T* probe_ljb_in,
            const T* probe_ljsig_in, const T* softcore_lja_in, const T* softcore_ljb_in,
            const T* softcore_ljc_in, const T* softcore_ljd_in, const T* softcore_lje_in,
            const T* softcore_ljf_in, const T* softcore_qq_in);

  MeshFFKit(const MeshFFKit<T> &original) = default;
  MeshFFKit(MeshFFKit<T> &&original) = default;
 
  // None of the member variables are mutable.
  const VdwCombiningRule ljrule;  
  const double coulomb;           
  const float coulomb_f;          
  const double clash_ratio;       
  const float clash_ratio_f;      
  const double clash_distance;    
  const float clash_distance_f;   
  const T* probe_lja;             
  const T* probe_ljb;             
  const T* probe_ljsig;           
  const T* softcore_lja;          
  const T* softcore_ljb;          
  const T* softcore_ljc;          
  const T* softcore_ljd;          
  const T* softcore_lje;          
  const T* softcore_ljf;          
  const T* softcore_qq;           
};

template <typename T> class MeshForceField {
public:

  MeshForceField(VdwCombiningRule mixing_protocol_in = VdwCombiningRule::LORENTZ_BERTHELOT,
                 double coulomb_constant_in = amber_ancient_bioq,
                 double clash_ratio_in = 0.0, double clash_distance_in = 0.0,
                 int lj_type_count_in = 0);
 
  MeshForceField(VdwCombiningRule mixing_protocol_in, double clash_ratio_in,
                 double clash_distance_in, const AtomGraph *ag);
 
  MeshForceField(VdwCombiningRule mixing_protocol_in, double clash_ratio_in,
                 double clash_distance_in, const AtomGraph &ag);
 
  MeshForceField(VdwCombiningRule mixing_protocol_in, double coulomb_constant_in,
                 double clash_ratio_in, double clash_distance_in,
                 const std::vector<double> &probe_sigma,
                 const std::vector<double> &probe_epsilon);

  MeshForceField(const MeshForceField<T> &original);
  MeshForceField(MeshForceField<T> &&original);
  MeshForceField& operator=(const MeshForceField<T> &other);
  MeshForceField& operator=(MeshForceField<T> &&other);

  VdwCombiningRule getCombiningRule() const;

  double getClashRatio() const;

  double getClashDistance() const;

  const MeshFFKit<T> data(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const MeshFFKit<double> referenceData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const MeshFFKit<void> templateFreeData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;
 
#ifdef STORMM_USE_HPC
#  ifdef STORMM_USE_CUDA
  const MeshFFKit<T> deviceViewToHostData() const;

  const MeshFFKit<double> deviceViewToReferenceHostData() const;

  const MeshFFKit<void> deviceViewToTemplateFreeHostData() const;
#  endif

  void upload();

  void download();
#endif

  void setCombiningRule(VdwCombiningRule mixing_protocol_in);

  void setClashDistance(double clash_distance_in);

  void setClashRatio(double clash_ratio_in);

  void setElecSoftcoreParameter(double coef, int pos);

  void setElecSoftcoreParameters(Interpolant stencil_kind);

  void setLJCoefficients(int index, double lja_in, double ljb_in);

  void setLJCoefficients(const AtomGraph *ag, const std::vector<double> &probe_radius,
                         const std::vector<double> &probe_well_depth);
 
  void setLJCoefficients(const AtomGraph &ag, const std::vector<double> &probe_radius,
                         const std::vector<double> &probe_well_depth);
 
  void setLJCoefficients(const AtomGraph *ag, double probe_radius, double probe_well_depth);
 
  void setLJCoefficients(const AtomGraph &ag, double probe_radius, double probe_well_depth);

  void setLJSoftcoreParameter(int index, double coef, int pos);

  void setLJSoftcoreParameters(Interpolant stencil_kind);
  
  void reallocate(const int lj_type_count_in);
  
private:

  int lj_type_count;
  
  VdwCombiningRule mixing_protocol;

  double coulomb_constant;

  double clash_ratio;

  double clash_distance;
 
  // The following arrays store atom type-specific Lennard-Jones parameters for meshes detailing
  // steric clashes or Lennard-Jones potentials.  If the combining rule is NBFIX, the pairwise
  // parameters for interactions of the probe with each type of atom in the underlying topology
  // must be provided.  Otherwise, these parameters will be inferred based on the topology and a
  // single set of parameters for the solvent probe using one of the other combining rules.
  Hybrid<T> probe_lja;           
  Hybrid<T> probe_ljb;           
  Hybrid<T> probe_lj_sigma;      
  Hybrid<T> probe_softcore_lja;  
  Hybrid<T> probe_softcore_ljb;  
  Hybrid<T> probe_softcore_ljc;  
  Hybrid<T> probe_softcore_ljd;  
  Hybrid<T> probe_softcore_lje;  
  Hybrid<T> probe_softcore_ljf;  
  Hybrid<T> elec_softcore;       

  Hybrid<T> softcore_data;
  
  Hybrid<double> ref_probe_lja;
  Hybrid<double> ref_probe_ljb;
  Hybrid<double> ref_probe_lj_sigma;
  Hybrid<double> ref_probe_softcore_lja;
  Hybrid<double> ref_probe_softcore_ljb;
  Hybrid<double> ref_probe_softcore_ljc;
  Hybrid<double> ref_probe_softcore_ljd;
  Hybrid<double> ref_probe_softcore_lje;
  Hybrid<double> ref_probe_softcore_ljf;
  Hybrid<double> ref_elec_softcore;
  Hybrid<double> ref_softcore_data;

  void allocate(int lj_type_count_in);
  
  void rebasePointers();

  void validateTypeIndex(int index);
};

template <typename T> MeshFFKit<T> restoreType(const MeshFFKit<void> *rasa);
template <typename T> MeshFFKit<T> restoreType(const MeshFFKit<void> &rasa);
 
} // namespace structure
} // namespace stormm
 
#include "mesh_forcefield.tpp"
 
#endif