spring_functions.hpp

/*
  This file is part of LOOS.
 
  LOOS (Lightweight Object-Oriented Structure library)
  Copyright (c) 2010 Tod D. Romo
  Department of Biochemistry and Biophysics
  School of Medicine & Dentistry, University of Rochester
 
  This package (LOOS) is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation under version 3 of the License.
 
  This package is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
*/
 
 
#if !defined(LOOS_SPRING_FUNCTIONS_HPP)
#define LOOS_SPRING_FUNCTIONS_HPP
 
 
 
 
#include <loos.hpp>
 
 
namespace ENM {
 
  // Exceptions classes (primarily for springFactory())
  class BadSpringFunction : public std::runtime_error
  {
  public:
    BadSpringFunction(const std::string& s) : runtime_error(s) { }
  };
 
  class BadSpringParameter : public std::runtime_error
  {
  public:
    BadSpringParameter(const std::string& s) : runtime_error(s) { }
  };
 
 
 
 
  class SpringFunction {
  public:
    typedef std::vector<double>      Params;
  public:
    SpringFunction() : warned(false) { }
    virtual ~SpringFunction() { }
 
    virtual std::string name() const =0;
 
    virtual Params setParams(const Params& konst) =0;
 
    virtual bool validParams() const =0;
 
    virtual uint paramSize() const =0;
 
 
  
    virtual loos::DoubleMatrix constant(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d)  =0;
 
  protected:
 
 
    double checkConstant(double d) {
      if (d < 0.0) {
        if (!warned) {
          warned = true;
          std::cerr << "Warning- negative spring constants found in " << name() << ".  Setting to 0.\n";
        }
        d = 0.0;
      }
 
      return(d);
    }
 
  private:
    bool warned;
  };
 
 
 
 
  class UniformSpringFunction : public SpringFunction {
  public:
 
    loos::DoubleMatrix constant(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) {
      double k = checkConstant(constantImpl(u, v, d));
      loos::DoubleMatrix B(3, 3);
      for (uint i=0; i<9; ++i)
        B[i] = k;
 
      return(B);
    }
 
  private:
 
    virtual double constantImpl(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) =0;
  };
 
 
 
  // The following are all uniform spring constants...
 
 
 
  class DistanceCutoff : public UniformSpringFunction {
  public:
    DistanceCutoff(const double& r) : radius(r*r) { }
    DistanceCutoff() : radius(15.0*15.0) { }
 
    std::string name() const { return("DistanceCutoff"); }
 
    Params setParams(const Params& p) {
      if (p.size() < 1)
        throw(BadSpringParameter("Insufficient number of spring parameters"));
 
      Params q(p);
      radius = q.back();
      radius *= radius;
      q.pop_back();
      return(q);
    }
 
    bool validParams() const { return(radius > 0.0); }
 
    uint paramSize() const { return(1); }
 
    double constantImpl(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) {
      double s = d.length2();
      if (s <= radius)
        return(1./s);
 
      return(0.0);
    }
 
  private:
    double radius;
  };
 
 
  class DistanceWeight : public UniformSpringFunction {
  public:
    DistanceWeight(const double p) : power(p) { }
    DistanceWeight() : power(-2.0) { }
 
    std::string name() const { return("DistanceWeight"); }
 
 
    Params setParams(const Params& p) {
      if (p.size() < 1)
        throw(BadSpringParameter("Insufficient number of spring parameters"));
 
      Params q(p);
      power = q.back();
      q.pop_back();
      return(q);
    }
 
    bool validParams() const { return(power < 0.0); }
 
    uint paramSize() const { return(1); }
 
 
    double constantImpl(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) {
      double s = d.length();
      return(pow(s, power));
    }
 
  private:
    double power;
  };
 
 
 
 
  class ExponentialDistance : public UniformSpringFunction {
  public:
    ExponentialDistance(const double s) : scale(s) { }
    ExponentialDistance() : scale(-1.5) { }
 
    std::string name() const { return("ExponentialDistance"); }
 
    Params setParams(const Params& p) {
      if (p.size() < 1)
        throw(BadSpringParameter("Insufficient number of spring parameters"));
 
      Params q(p);
      scale = q.back();
      q.pop_back();
      return(q);
    }
 
    bool validParams() const { return(scale != 0.0); }
 
    uint paramSize() const { return(1); }
 
 
    double constantImpl(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) {
      double s = d.length();
      return(exp(scale * s));
    }
 
  private:
    double scale;
  };
 
 
 
 
 
 
  class HCA : public UniformSpringFunction {
  public:
    HCA(const double rc, const double a, const double b, const double c, const double d) :
      rcut(rc), k1(a), k2(b), k3(c), k4(d) { }
 
    HCA() :
      rcut(4.0), k1(205.5), k2(571.2), k3(305.9e3), k4(6.0) { }
 
    std::string name() const { return("HCA"); }
 
    Params setParams(const Params& p) {
      if (p.size() < 5)
        throw(BadSpringParameter("Insufficient number of spring parameters"));
 
      Params q(p);
    
      k4 = q.back();
      q.pop_back();
 
      k3 = q.back();
      q.pop_back();
 
      k2 = q.back();
      q.pop_back();
 
      k1 = q.back();
      q.pop_back();
 
      rcut = q.back();
      q.pop_back();
 
      return(q);
    }
 
    bool validParams() const { return(rcut >= 0.0 && k4 >= 0.0); }
 
    uint paramSize() const { return(5); }
 
 
    double constantImpl(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) {
      double s = d.length();
      double k;
 
      if (s <= rcut)
        k = k1 * s - k2;
      else
        k = k3 * pow(s, -k4);
 
      return(k);
    }
 
  private:
    double rcut, k1, k2, k3, k4;
  };
 
class ConstBonded : public UniformSpringFunction {
public:
  ConstBonded(const double& s) : scale(s) { }
  ConstBonded() : scale(1) { }
  
  std::string name() const { return("ConstBonded"); }
 
  Params setParams(const Params& p) {
    if (p.size() < 1)
      throw(BadSpringParameter("Insufficient number of spring parameters"));
 
    Params q(p);
    scale = q.back();
    q.pop_back();
    return(q);
  }
 
  bool validParams() const { return(scale > 0.0); }
 
  uint paramSize() const { return(1); }
 
  double constantImpl(const loos::GCoord& u, const loos::GCoord& v, const loos::GCoord& d) {
    //std::cerr << "In impl in constbonded :)\n";
    return(scale);
  }
 
private:
  double scale;
};
 
  SpringFunction* springFactory(const std::string& spring_desc);
 
  std::vector<std::string> springNames();
 
};
 
#endif