restraint_builder.h

// -*-c++-*-
#ifndef STORMM_RESTRAINT_BUILDER_H
#define STORMM_RESTRAINT_BUILDER_H
 
#include <vector>
#include "copyright.h"
#include "Chemistry/atommask.h"
#include "bounded_restraint.h"
 
namespace stormm {
namespace restraints {
 
using chemistry::AtomMask;
using topology::AtomGraph;
using trajectory::CoordinateFrameReader;

constexpr double default_hb_prevention_limit = 3.1;
  
struct FlatBottomPlan {

  FlatBottomPlan();
  FlatBottomPlan(double k_in, double r_in);
  FlatBottomPlan(double k_in, double r2_in, double r3_in);
  FlatBottomPlan(double k2_in, double k3_in, double r2_in, double r3_in);
  FlatBottomPlan(double k2_in, double k3_in, double r1_in, double r2_in, double r3_in,
                 double r4_in, int step_in);
  
  int activation_step;  // Step number at which this state of the restraint will take (full) effect
  double k2;            // Left-hand harmonic stiffness
  double k3;            // Right-hand harmonic stiffness
  double r1;            // Leftmost boundary between linear and harmonic behavior
  double r2;            // Leftmost boundary of the flat-bottom potential
  double r3;            // Rightmost boundary of the flat-bottom potential
  double r4;            // Rightmost boundard between harmonic and linear behavior
};
  
BoundedRestraint applyDistanceRestraint(const AtomGraph *ag, int atom_i, int atom_j,
                                        const FlatBottomPlan fb_init);
BoundedRestraint applyDistanceRestraint(const AtomGraph *ag, int atom_i, int atom_j,
                                        const FlatBottomPlan fb_init,
                                        const FlatBottomPlan fb_final);

BoundedRestraint applyAngleRestraint(const AtomGraph *ag, int atom_i, int atom_j, int atom_k,
                                     const FlatBottomPlan fb_init);
BoundedRestraint applyAngleRestraint(const AtomGraph *ag, int atom_i, int atom_j, int atom_k,
                                     const FlatBottomPlan fb_init, const FlatBottomPlan fb_final);

BoundedRestraint applyDihedralRestraint(const AtomGraph *ag, int atom_i, int atom_j, int atom_k,
                                        int atom_l, const FlatBottomPlan fb_init);
BoundedRestraint applyDihedralRestraint(const AtomGraph *ag, int atom_i, int atom_j, int atom_k,
                                        int atom_l, const FlatBottomPlan fb_init,
                                        const FlatBottomPlan fb_final);

void restraintTopologyChecks(const AtomGraph *ag, const CoordinateFrameReader &cframe,
                             const AtomMask &mask);
  
std::vector<BoundedRestraint>
applyPositionalRestraints(const AtomGraph *ag, const CoordinateFrameReader &ref_cfr,
                          const AtomMask &mask, double displacement_penalty,
                          double displacement_onset = 0.0, double displacement_plateau = 16.0,
                          double proximity_penalty = 0.0, double proximity_onset = 0.0,
                          double proximity_plateau = 0.0);
 
std::vector<BoundedRestraint>
applyPositionalRestraints(const AtomGraph *ag, const CoordinateFrameReader &ref_cfr,
                          const std::string &mask, double displacement_penalty,
                          double displacement_onset = 0.0, double displacement_plateau = 16.0,
                          double proximity_penalty = 0.0, double proximity_onset = 0.0,
                          double proximity_plateau = 0.0);
 
std::vector<BoundedRestraint>
applyPositionalRestraints(const AtomGraph *ag, const CoordinateFrameReader &ref_cfr,
                          const std::vector<int> &masked_atoms, double displacement_penalty,
                          double displacement_onset = 0.0, double displacement_plateau = 16.0,
                          double proximity_penalty = 0.0, double proximity_onset = 0.0,
                          double proximity_plateau = 0.0);


std::vector<BoundedRestraint>
applyDistanceRestraints(const AtomGraph *ag, const CoordinateFrameReader &ref_cfr,
                        const AtomMask &mask, double penalty, double flat_bottom_half_width,
                        double cutoff);
 
std::vector<BoundedRestraint>
applyDistanceRestraints(const AtomGraph *ag, const CoordinateFrameReader &ref_cfr,
                        const std::string &mask, double penalty, double flat_bottom_half_width,
                        double cutoff);
 
std::vector<BoundedRestraint>
applyDistanceRestraints(const AtomGraph *ag, const CoordinateFrameReader &ref_cfr,
                        const std::vector<int> &masked_atoms, double penalty,
                        double flat_bottom_half_width, double cutoff);
  
std::vector<BoundedRestraint>
applyHoldingRestraints(const AtomGraph *ag, const CoordinateFrameReader &cfr,
                       const AtomMask &mask, double penalty, double flat_bottom_half_width = 0.0,
                       double harmonic_maximum = 16.0);
 
std::vector<BoundedRestraint>
applyHoldingRestraints(const AtomGraph *ag, const CoordinateFrameReader &cfr,
                       const std::string &mask, double penalty,
                       double flat_bottom_half_width = 0.0, double harmonic_maximum = 16.0);

std::vector<BoundedRestraint>
applyHydrogenBondPreventors(const AtomGraph *ag, const ChemicalFeatures &chemfe,
                            const double penalty,
                            const double approach_point = default_hb_prevention_limit);
 
std::vector<BoundedRestraint>
applyHydrogenBondPreventors(const AtomGraph *ag, const CoordinateFrameReader &cfr,
                            const double penalty,
                            const double approach_point = default_hb_prevention_limit);
  
} // namespace restraints
} // namespace stormm
 
#endif