IMP  2.0.1
The Integrative Modeling Platform
WormLikeChain.h
Go to the documentation of this file.
1 /**
2  * \file IMP/misc/WormLikeChain.h
3  * \brief Worm-like-chain score for polymer chains.
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPMISC_WORM_LIKE_CHAIN_H
9 #define IMPMISC_WORM_LIKE_CHAIN_H
10 
11 #include <IMP/misc/misc_config.h>
12 #include <IMP/UnaryFunction.h>
13 #include <IMP/internal/constants.h>
14 #include <IMP/internal/units.h>
15 
16 IMPMISC_BEGIN_NAMESPACE
17 
18 //! Worm-like-chain energy for polymer chains
19 /** This function implements one polymer force/extension curve. It
20  is a physical energy and all the inputs are in angstroms
21  and the outputs in kcal/mol (/angstrom).
22 
23  \note The actual worm like chain force blows up as the extension approaches
24  the maximum. Since that makes optimization problematic, we
25  just linearly extend the force after 99% of maximum.
26  */
28 {
29 public:
30  //! Define the energy term
31  /** \param[in] l_max maximum length of the chain in angstroms
32  \param[in] lp persistence length in angstroms
33  */
34  WormLikeChain(Float l_max, Float lp) : lmax_(l_max), lp_(lp) {
35  IMP_USAGE_CHECK(l_max > lp, "The persistence length should be less "
36  << "than the total length for this model");
37  }
38 
39  virtual DerivativePair evaluate_with_derivative(double feature) const;
40 
41  virtual double evaluate(double feature) const;
42 
44 
45  void do_show(std::ostream &out) const;
46 
47 private:
48  unit::Piconewton cderiv(unit::Angstrom l) const {
49  unit::Piconewton pn= IMP::internal::KB*IMP::internal::DEFAULT_TEMPERATURE
50  /lp_*(.25/ square(1.0-(l/lmax_))
51  -.25+(l/lmax_));
52  return pn;
53  }
54 
55  unit::Picojoule eval(unit::Angstrom m) const {
56  unit::Picojoule J
57  = IMP::internal::KB *
58  IMP::internal::DEFAULT_TEMPERATURE/lp_*(.25*square(lmax_)
59  /(lmax_-m)
60  -m*.25
61  +.5*square(m)
62  /lmax_);
63  return J;
64  }
65 
66  unit::Angstrom cutoff() const {
67  return .99*lmax_;
68  }
69 
70  unit::Angstrom lmax_, lp_;
71 };
72 
73 #ifndef IMP_DOXYGEN
74 inline double WormLikeChain::evaluate(double v) const {
75  return evaluate_with_derivative(v).first;
76 }
77 
78 inline DerivativePair WormLikeChain::evaluate_with_derivative(double v) const {
79  static const unit::Picojoule zero=eval(unit::Angstrom(0));
80  unit::Angstrom l(v);
81  if (l < unit::Angstrom(0)) l=unit::Angstrom(0);
82  unit::Picojoule ret;
83 
84  unit::Piconewton doubled;
85  if (l < cutoff()) {
86  ret= (eval(l) - zero);
87  doubled= cderiv(l);
88  } else {
89  unit::Picojoule springterm=(l-cutoff())*cderiv(cutoff());
90  ret= (eval(cutoff())+ springterm -zero);
91  doubled= cderiv(cutoff());
92  IMP_LOG_VERBOSE( "Overstretched " << cderiv(cutoff()) << " " << doubled
93  << " " << l << " " << lmax_ << " " << cutoff()
94  << std::endl);
95  }
96  unit::YoctoKilocalorie zc= convert_J_to_Cal(ret);
97  double value=(zc*unit::ATOMS_PER_MOL).get_value();
98  //std::cout << "Force is " << doubled << " for length " << l << std::endl;
99  // convert from picoNewton
100  unit::YoctoKilocaloriePerAngstrom du= unit::convert_J_to_Cal(doubled);
101 
102  double deriv = (du*unit::ATOMS_PER_MOL).get_value();
103  //std::cout << "Which converts to " << d << std::endl;
104  return std::make_pair(value, deriv);
105 }
106 
107 inline void WormLikeChain::do_show(std::ostream &out) const {
108  out << "params " << lmax_ << " " << lp_ << std::endl;
109 }
110 #endif
111 
112 IMPMISC_END_NAMESPACE
113 
114 #endif /* IMPMISC_WORM_LIKE_CHAIN_H */