IMP logo
IMP Reference Guide  2.13.0
The Integrative Modeling Platform
linear_distance_pair_scores.h
Go to the documentation of this file.
1 /**
2  * \file linear_distance_pair_scores.h
3  * \brief A Score on the distance between a pair of particles.
4  *
5  * Copyright 2007-2020 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPNPCTRANSPORT_LINEAR_DISTANCE_PAIR_SCORES_H
9 #define IMPNPCTRANSPORT_LINEAR_DISTANCE_PAIR_SCORES_H
10 
11 #include "npctransport_config.h"
12 #include <IMP/PairScore.h>
13 #include <IMP/pair_macros.h>
14 #include <IMP/core/XYZR.h>
15 #include <IMP/algebra/utility.h>
16 #include "internal/sites.h"
17 
18 #include <boost/array.hpp>
19 
20 IMPNPCTRANSPORT_BEGIN_NAMESPACE
21 
22 #ifndef SWIG
23 /**
24  Evaluates a linear pair potential, such that score = k * (delta_length - x0)
25  is returned.
26  Also, updates the derivative vectors of the particles if da is not null.
27 
28  @param[out] d_xyzr0 pointer to derivative of particle 0
29  @param[out] d_xyzr1 pointer to derivative of particle 1
30  @param[in,out] da accumulator for score derivatives to be updated
31  @param[in] delta a vector that represents the displacement between the
32  two particles
33  @param delta_length the cached length of delta, assumed correct, and required
34  for faster calculation
35  @param x0 resting distance (where score = 0)
36  @param k score linear coefficient. Note that the score is attractive
37  for a positive k, and repulsive for a negative k
38  (assuming the lower the score the better).
39  @return the score
40  */
41 inline double do_evaluate_index(algebra::Sphere3D& d_xyzr0,
42  algebra::Sphere3D& d_xyzr1,
44  const algebra::Vector3D &delta,
45  double delta_length, double x0, double k) {
46  double shifted_length = delta_length - x0;
47  double score = k * shifted_length;
48  static const double MIN_DISTANCE = .00001;
49  if (da && delta_length > MIN_DISTANCE) {
50  algebra::Vector3D deriv = k * delta / delta_length; // deriv magnitude is k
51  Model::add_to_coordinate_derivatives(d_xyzr0, deriv, *da);
52  Model::add_to_coordinate_derivatives(d_xyzr1, -deriv, *da);
53  IMP_LOG(TERSE, "Distance: " << shifted_length << "\nscore: " << score
54  << "\nderiv: " << deriv << std::endl);
55  }
56  IMP_CHECK_CODE ( else {
57  IMP_LOG(TERSE, "Distance: " << shifted_length << "\nscore: " << score
58  << std::endl);
59  } );
60  return score;
61 }
62 #endif
63 
64 /**
65  A soft linear repulsive score between two spheres.
66  If the spheres penetrate each other, the score is
67  ( k_ * [penetration-magnitude] kcal/mol), otherwise the score is 0.
68  */
69 class IMPNPCTRANSPORTEXPORT LinearSoftSpherePairScore : public PairScore {
70  private:
71  // score potential coefficient (larger = stronger repulsion)
72  double k_;
73 
74  protected:
75  //! evaluates the soft linear repulsive score, for spheres s0 and s1,
76  //! outputing derivatives to ds0 and ds1 scaled by da, if da is not null
77  inline double evaluate_index
78  (algebra::Sphere3D const& s0,
79  algebra::Sphere3D const& s1,
80  algebra::Sphere3D& ds0,
81  algebra::Sphere3D& ds1,
82  DerivativeAccumulator *da) const;
83 
84  public:
85 
87  std::string name = "LinearSSPairScore%1%");
88 
89  virtual double evaluate_index(Model *m, const ParticleIndexPair &p,
91 
92  virtual double evaluate_indexes(Model *m,
93  const ParticleIndexPairs &pips,
95  unsigned int lower_bound,
96  unsigned int upper_bound) const IMP_OVERRIDE;
97 
99  ( Model *m,
100  const ParticleIndexPairs &p,
102  double max, unsigned int lower_bound, unsigned int upper_bound) const
103  {
104  double ret = 0;
105  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
106  ret += evaluate_if_good_index(m, p[i], da, max - ret);
107  if (ret > max) return std::numeric_limits<double>::max();
108  }
109  return ret;
110  }
111 
112 
114  const ParticleIndexes &pis) const;
115 
116  //! returns the k for sphere-sphere repulsion
117  double get_k() const { return k_; }
118 
120  ;
121 };
122 
123 #ifndef IMP_DOXYGEN
124 
125 //! evaluates the soft linear repulsive score, for spheres s0 and s1,
126 //! outputing derivatives to ds0 and ds1 scaled by da, if da is not null
127 inline double
128 LinearSoftSpherePairScore::evaluate_index
129 (algebra::Sphere3D const& s0,
130  algebra::Sphere3D const& s1,
131  algebra::Sphere3D& ds0,
132  algebra::Sphere3D& ds1,
133  DerivativeAccumulator *da) const
134 {
136  algebra::Vector3D delta =
137  s0.get_center() - s1.get_center();
138  double delta_length_2 =
139  delta.get_squared_magnitude(); // (work with sqr for efficiency)
140  double x0 =
141  s0.get_radius() + s1.get_radius();
142  double x0_2 = x0 * x0;
143  bool not_penetrating = delta_length_2 > x0_2;
144  if (not_penetrating) return 0;
145  // penetrating spheres ; -k_ in order to get a repulsive score
146  double delta_length = std::sqrt(delta_length_2);
147  return do_evaluate_index(ds0, ds1, da,
148  delta, delta_length, x0, -k_);
149 }
150 
151 // evaluates the soft linear repulsive score, using the internal Model
152 // indexing to access particles data (see Model/include/internal/)
153 inline double
154 LinearSoftSpherePairScore::evaluate_index
155 ( Model *m,
156  const ParticleIndexPair &pp,
157  DerivativeAccumulator *da) const
158 {
160  algebra::Sphere3D const& s0(m->get_sphere(pp[0]));
161  algebra::Sphere3D const& s1(m->get_sphere(pp[1]));
162  algebra::Sphere3D* d_xyzrs=
163  m->access_sphere_derivatives_data();
164  algebra::Sphere3D& ds0( d_xyzrs[pp[0].get_index()] );
165  algebra::Sphere3D& ds1( d_xyzrs[pp[1].get_index()] );
166  return evaluate_index(s0, s1, ds0, ds1, da);
167 }
168 
169 //! evaluate all index pairs between pis[lower_bound] and pis[upper_bound]
170 //! and update their coordinate derivatives scaled by da, if da is not null
171 inline double
172 LinearSoftSpherePairScore::evaluate_indexes
173 ( Model *m,
174  const ParticleIndexPairs &pips,
175  DerivativeAccumulator *da,
176  unsigned int lower_bound,
177  unsigned int upper_bound ) const
178 {
180  algebra::Sphere3D const* xyzrs=
181  m->access_spheres_data();
182  algebra::Sphere3D* d_xyzrs=
183  m->access_sphere_derivatives_data();
184  double ret = 0;
185  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
186  int i0(pips[i][0].get_index());
187  int i1(pips[i][1].get_index());
188  ret += evaluate_index( xyzrs[i0], xyzrs[i1],
189  d_xyzrs[i0], d_xyzrs[i1],
190  da);
191  }
192  return ret;
193 }
194 
195 #endif
196 
197 /**
198  A soft linear attractive / repulsive score between two spheres.
199  The score is 0 if the spheres are beyond the attractive range.
200  Within the attractive range, the score decreases linearly (= attraction)
201  with slope k_attr_ until the spheres touch. Once the spheres begin to
202  penetrate each other, the score rises linearly with slope k_rep_
203  (= repulsion), though it may be negative for small penetration.
204 */
205 class IMPNPCTRANSPORTEXPORT LinearInteractionPairScore : public PairScore {
206 #ifndef SWIG
207  public:
208  /** cache for reusable intermediate calculation */
210  // the squared-distance between particle centers:
211  double particles_delta_squared;
212  // the sum of particle radii:
213  double sum_particles_radii;
214  };
215 #endif
216 
217  private:
218  double range_attr_, // range of attraction between particles
219  k_rep_, // repulsion coeficcient
220  k_attr_; // attraction coefficient
221 
222  // cache for reusable intermediate calcualtions of evaluate_infex
223  mutable EvaluationCache cache_;
224 
225  protected:
226  //! evaluates the linear interaction score, for spheres s0 and s1,
227  //! outputing derivatives to ds0 and ds1 scaled by da, if da is not null
228  inline double evaluate_index
229  (algebra::Sphere3D const& s0,
230  algebra::Sphere3D const& s1,
231  algebra::Sphere3D& ds0,
232  algebra::Sphere3D& ds1,
233  DerivativeAccumulator *da) const;
234 
235  public:
236  /**
237  The score is 0 if the spheres are beyond the attractive range.
238  Within the attractive range, the score decreases linearly (= attraction)
239  with slope k_attr_ until the spheres touch. Once the spheres begin to
240  penetrate each other, the score rises linearly with slope k_rep_
241  (= repulsion), though it may be negative for small penetration.
242 
243  The energy potential difference between the unbound state and when the
244  beads touch is:
245  DELTA-U=0.5*k_attr*range_attr [kcal/mol]
246 
247  @param k_rep repulsion force constant in kcal/mol/A units
248  @param range_attr attraction range between sphere surfaces
249  @param k_attr attraction force constant in kcal/mol/A unit
250  @param name score object names
251  */
252  LinearInteractionPairScore(double k_rep, double range_attr, double k_attr,
253  std::string name = "LinearIDPairScore%1%");
254 
255 #ifndef SWIG
256  // returns cached intermediate computations from last call to
257  // this->evaluate_index()
258  EvaluationCache const &get_evaluation_cache() const { return cache_; }
259 #endif
260 
262  DerivativeAccumulator *da) const
263  IMP_OVERRIDE);
264  IMP_IMPLEMENT_INLINE(double evaluate_if_good_index
265  ( Model *m,
266  const ParticleIndexPair &p,
268  double max ) const,
269  {
270  IMP_UNUSED(max);
271  return evaluate_index(m, p, da);
272  }
273  );
274 
275  //! evaluate all index pairs between pips[lower_bound] and pis[upper_bound]
276  //! and update their coordinate derivatives scaled by da, if da is not null
277  virtual double evaluate_indexes(Model *m,
278  const ParticleIndexPairs &pips,
280  unsigned int lower_bound,
281  unsigned int upper_bound) const IMP_OVERRIDE;
282 
283  double evaluate_if_good_index(Model *m, const ParticleIndexPairs &p,
284  DerivativeAccumulator *da, double max,
285  unsigned int lower_bound,
286  unsigned int upper_bound) const {
287  double ret = 0;
288  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
289  ret += evaluate_if_good_index(m, p[i], da, max - ret);
290  if (ret > max) return std::numeric_limits<double>::max();
291  }
292  return ret;
293  }
294 
295  ModelObjectsTemp do_get_inputs(Model *m, const ParticleIndexes &pis) const
296  IMP_OVERRIDE;
297 
298  //! returns the range for sphere-sphere attraction in A
299  double get_range_attraction() const { return range_attr_; }
300 
301  //! returns the k for sphere-sphere attraction in kcal/mol/A units
302  double get_k_attraction() const { return k_attr_; }
303 
304  //! set k for sphere-sphere attraction in kcal/mol/A units
305  void set_k_attraction(double k_attr) {
306  k_attr_= k_attr;
307  }
308 
309  //! returns the k for sphere-sphere repulsion in kcal/mol/A units
310  double get_k_repulsion() const { return k_rep_; }
311 
312  //! set k for sphere-sphere repulsion in kcal/mol/A units
313  void set_k_repulsion(double k_rep) {
314  k_rep_ = k_rep;
315  }
316 
318 };
319 
320 #ifndef IMP_DOXYGEN
321 
322 //! evaluates the linear interaction score, for spheres s0 and s1,
323 //! outputing derivatives to ds0 and ds1 scaled by da, if da is not null
324 inline double
325 LinearInteractionPairScore::evaluate_index
326 (algebra::Sphere3D const& s0,
327  algebra::Sphere3D const& s1,
328  algebra::Sphere3D& ds0,
329  algebra::Sphere3D& ds1,
330  DerivativeAccumulator *da) const
331 {
332  // Associate intermediate variables with cache_, for further reuse:
333  double &delta_length_2 = cache_.particles_delta_squared;
334  double &x0 = cache_.sum_particles_radii;
335  algebra::Vector3D delta =
336  s0.get_center() - s1.get_center();
337  delta_length_2 = delta.get_squared_magnitude();
338  IMP_LOG(PROGRESS,
339  "LinearInteractionPairScore cached delta2 "
340  << cache_.particles_delta_squared << std::endl);
341  x0 = s0.get_radius() + s1.get_radius();
342  // Terminate immediately if very far, work with squares for speed
343  // equivalent to [delta_length > x0 + attr_range]:
344  if (delta_length_2 > std::pow(x0 + range_attr_, 2)) return 0;
345  double offset = -range_attr_ * k_attr_;
346  double delta_length = std::sqrt(delta_length_2);
347  if (delta_length > x0) { // attractive regime
348  return // decreases with slope k_attr_ as spheres get closer
349  (do_evaluate_index(ds0, ds1, da,
350  delta, delta_length, x0, k_attr_) +
351  offset);
352  } else { // repulsive regime, may be negative for small penetration
353  return (do_evaluate_index(ds0, ds1, da,
354  delta, delta_length, x0, -k_rep_) +
355  offset);
356  }
357 }
358 
359 inline double
360 LinearInteractionPairScore::evaluate_index
361 ( Model *m,
362  const ParticleIndexPair &pp,
363  DerivativeAccumulator *da) const
364 {
366  algebra::Sphere3D const& s0(m->get_sphere(pp[0]));
367  algebra::Sphere3D const& s1(m->get_sphere(pp[1]));
368  algebra::Sphere3D* d_xyzrs=
369  m->access_sphere_derivatives_data();
370  algebra::Sphere3D& ds0( d_xyzrs[pp[0].get_index()] );
371  algebra::Sphere3D& ds1( d_xyzrs[pp[1].get_index()] );
372  return evaluate_index(s0, s1, ds0, ds1, da);
373  }
374 
375 inline double
376 LinearInteractionPairScore::evaluate_indexes
377 ( Model *m,
378  const ParticleIndexPairs &pips,
379  DerivativeAccumulator *da,
380  unsigned int lower_bound,
381  unsigned int upper_bound) const
382 {
384  algebra::Sphere3D const* xyzrs=
385  m->access_spheres_data();
386  algebra::Sphere3D* d_xyzrs=
387  m->access_sphere_derivatives_data();
388  double ret = 0;
389  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
390  int i0(pips[i][0].get_index());
391  int i1(pips[i][1].get_index());
392  ret += evaluate_index( xyzrs[i0], xyzrs[i1],
393  d_xyzrs[i0], d_xyzrs[i1],
394  da);
395  }
396  return ret;
397 }
398 
399 
400 #endif
401 
402 /**
403  A soft linear attractive / repulsive score between two spheres
404  for the backbone of a chain.
405 */
406 class IMPNPCTRANSPORTEXPORT LinearWellPairScore : public PairScore {
407  private:
408  double rest_length_factor_, k_;
409 
410  public:
411  /**
412  a linear well pair potential that keeps two particles around
413  a resting distance relative to their radii
414 
415  @param rest_length_factor the resting distance between particles
416  relative to their sum of radii
417  @param k the force constant (attractive if beyond or repulsive
418  if below rest length) in units kcal/mol/A
419  @param name the name of the score
420  */
421  LinearWellPairScore(double rest_length_factor, double k,
422  std::string name = "LinearIDPairScore%1%");
423 
424  void set_rest_length_factor(double rest_length_factor)
425  { rest_length_factor_ = rest_length_factor; }
426  double get_rest_length_factor() const { return rest_length_factor_; }
427  void set_k(double k)
428  { k_ = k; }
429  double get_k() { return k_; }
430  double evaluate_index(Model *m, const ParticleIndexPair &p,
432  ModelObjectsTemp do_get_inputs(Model *m, const ParticleIndexes &pis) const;
435  ;
436 };
437 
438 #ifndef IMP_DOXYGEN
439 inline double
440 LinearWellPairScore::evaluate_index
441 ( Model *m,
442  const ParticleIndexPair &pp,
443  DerivativeAccumulator *da) const
444 {
446  algebra::Sphere3D const& s0 = m->get_sphere(pp[0]);
447  algebra::Sphere3D const& s1 = m->get_sphere(pp[1]);
448  algebra::Sphere3D* d_xyzrs=
449  m->access_sphere_derivatives_data();
450  algebra::Sphere3D& ds0( d_xyzrs[pp[0].get_index()] );
451  algebra::Sphere3D& ds1( d_xyzrs[pp[1].get_index()] );
452  double x0 = (s0.get_radius() + s1.get_radius()) * rest_length_factor_;
453  algebra::Vector3D delta = s0.get_center() - s1.get_center();
454  double delta_length_2 = delta.get_squared_magnitude();
455  double delta_length = std::sqrt(delta_length_2);
456  if (delta_length > x0) { // attractive regime
457  return // k_ > 0 = get spheres closer
458  do_evaluate_index(ds0, ds1, da,
459  delta, delta_length, x0, k_);
460  } else {
461  return// -k_ < 0 = keep spheres apart
462  do_evaluate_index(ds0, ds1, da,
463  delta, delta_length, x0, -k_);
464  }
465 }
466 #endif
467 
469 
470 
471 IMPNPCTRANSPORT_END_NAMESPACE
472 
473 #endif /* IMPNPCTRANSPORT_LINEAR_DISTANCE_PAIR_SCORES_H */
#define IMP_IMPLEMENT_INLINE(signature, body)
void set_k_attraction(double k_attr)
set k for sphere-sphere attraction in kcal/mol/A units
Abstract class for scoring object(s) of type ParticleIndexPair.
Definition: PairScore.h:37
SphereD< 3 > Sphere3D
Typedef for Python.
Definition: SphereD.h:96
virtual double evaluate_if_good_indexes(Model *m, const ParticleIndexPairs &o, DerivativeAccumulator *da, double max, unsigned int lower_bound, unsigned int upper_bound) const
double get_range_attraction() const
returns the range for sphere-sphere attraction in A
Macros for various classes.
#define IMP_PAIR_SCORE_METHODS(Name)
Definition: pair_macros.h:25
IMP::Vector< ParticleIndexPair > ParticleIndexPairs
Definition: base_types.h:161
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
#define IMP_OBJECT_LOG
Set the log level to the object's log level.
Definition: log_macros.h:284
double get_k_attraction() const
returns the k for sphere-sphere attraction in kcal/mol/A units
#define IMP_IMPLEMENT(signature)
IMP::Vector< IMP::WeakPointer< ModelObject > > ModelObjectsTemp
Definition: base_types.h:82
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Functions to deal with very common math operations.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
#define IMP_UNUSED(variable)
Define PairScore.
virtual ModelObjectsTemp do_get_inputs(Model *m, const ParticleIndexes &pis) const =0
Overload this method to specify the inputs.
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing lists of object pointers.
Definition: object_macros.h:44
#define IMP_CHECK_CODE(expr)
Only compile the code if checks are enabled.
Definition: check_macros.h:117
virtual double evaluate_indexes(Model *m, const ParticleIndexPairs &o, DerivativeAccumulator *da, unsigned int lower_bound, unsigned int upper_bound) const
Compute the score and the derivative if needed over a set.
VectorD< 3 > Vector3D
Definition: VectorD.h:421
void set_k_repulsion(double k_rep)
set k for sphere-sphere repulsion in kcal/mol/A units
virtual double evaluate_if_good_index(Model *m, const ParticleIndexPair &vt, DerivativeAccumulator *da, double max) const
Compute the score and the derivative if needed, only if "good".
Decorator for a sphere-like particle.
double do_evaluate_index(algebra::Sphere3D &d_xyzr0, algebra::Sphere3D &d_xyzr1, DerivativeAccumulator *da, const algebra::Vector3D &delta, double delta_length, double x0, double k)
double get_k_repulsion() const
returns the k for sphere-sphere repulsion in kcal/mol/A units
virtual double evaluate_index(Model *m, const ParticleIndexPair &vt, DerivativeAccumulator *da) const =0
Compute the score and the derivative if needed.
#define IMP_OVERRIDE
Cause a compile error if this method does not override a parent method.
Class for adding derivatives from restraints to the model.
double get_k() const
returns the k for sphere-sphere repulsion