IMP logo
IMP Reference Guide  develop.4390929e3e,2024/11/12
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-2022 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  //! outputting 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,
90  DerivativeAccumulator *da) const override;
91 
92  virtual double evaluate_indexes(Model *m,
93  const ParticleIndexPairs &pips,
95  unsigned int lower_bound,
96  unsigned int upper_bound) const override;
97 
99  ( Model *m,
100  const ParticleIndexPairs &p,
102  double max, unsigned int lower_bound,
103  unsigned int upper_bound) const override
104  {
105  double ret = 0;
106  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
107  ret += evaluate_if_good_index(m, p[i], da, max - ret);
108  if (ret > max) return std::numeric_limits<double>::max();
109  }
110  return ret;
111  }
112 
113 
115  const ParticleIndexes &pis) const override;
116 
117  //! returns the k for sphere-sphere repulsion
118  double get_k() const { return k_; }
119 
121  ;
122 };
123 
124 #ifndef IMP_DOXYGEN
125 
126 //! evaluates the soft linear repulsive score, for spheres s0 and s1,
127 //! outputting derivatives to ds0 and ds1 scaled by da, if da is not null
128 inline double
129 LinearSoftSpherePairScore::evaluate_index
130 (algebra::Sphere3D const& s0,
131  algebra::Sphere3D const& s1,
132  algebra::Sphere3D& ds0,
133  algebra::Sphere3D& ds1,
134  DerivativeAccumulator *da) const
135 {
137  algebra::Vector3D delta =
138  s0.get_center() - s1.get_center();
139  double delta_length_2 =
140  delta.get_squared_magnitude(); // (work with sqr for efficiency)
141  double x0 =
142  s0.get_radius() + s1.get_radius();
143  double x0_2 = x0 * x0;
144  bool not_penetrating = delta_length_2 > x0_2;
145  if (not_penetrating) return 0;
146  // penetrating spheres ; -k_ in order to get a repulsive score
147  double delta_length = std::sqrt(delta_length_2);
148  return do_evaluate_index(ds0, ds1, da,
149  delta, delta_length, x0, -k_);
150 }
151 
152 // evaluates the soft linear repulsive score, using the internal Model
153 // indexing to access particles data (see Model/include/internal/)
154 inline double
155 LinearSoftSpherePairScore::evaluate_index
156 ( Model *m,
157  const ParticleIndexPair &pp,
158  DerivativeAccumulator *da) const
159 {
161  algebra::Sphere3D const& s0(m->get_sphere(pp[0]));
162  algebra::Sphere3D const& s1(m->get_sphere(pp[1]));
163  algebra::Sphere3D* d_xyzrs=
164  m->access_sphere_derivatives_data();
165  algebra::Sphere3D& ds0( d_xyzrs[pp[0].get_index()] );
166  algebra::Sphere3D& ds1( d_xyzrs[pp[1].get_index()] );
167  return evaluate_index(s0, s1, ds0, ds1, da);
168 }
169 
170 //! evaluate all index pairs between pis[lower_bound] and pis[upper_bound]
171 //! and update their coordinate derivatives scaled by da, if da is not null
172 inline double
173 LinearSoftSpherePairScore::evaluate_indexes
174 ( Model *m,
175  const ParticleIndexPairs &pips,
176  DerivativeAccumulator *da,
177  unsigned int lower_bound,
178  unsigned int upper_bound ) const
179 {
181  algebra::Sphere3D const* xyzrs=
182  m->access_spheres_data();
183  algebra::Sphere3D* d_xyzrs=
184  m->access_sphere_derivatives_data();
185  double ret = 0;
186  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
187  int i0(pips[i][0].get_index());
188  int i1(pips[i][1].get_index());
189  ret += evaluate_index( xyzrs[i0], xyzrs[i1],
190  d_xyzrs[i0], d_xyzrs[i1],
191  da);
192  }
193  return ret;
194 }
195 
196 #endif
197 
198 /**
199  A soft linear attractive / repulsive score between two spheres.
200  The score is 0 if the spheres are beyond the attractive range.
201  Within the attractive range, the score decreases linearly (= attraction)
202  with slope k_attr_ until the spheres touch. Once the spheres begin to
203  penetrate each other, the score rises linearly with slope k_rep_
204  (= repulsion), though it may be negative for small penetration.
205 */
206 class IMPNPCTRANSPORTEXPORT LinearInteractionPairScore : public PairScore {
207 #ifndef SWIG
208  public:
209  /** cache for reusable intermediate calculation */
211  // the squared-distance between particle centers:
212  double particles_delta_squared;
213  // the sum of particle radii:
214  double sum_particles_radii;
215  };
216 #endif
217 
218  private:
219  double range_attr_, // range of attraction between particles
220  k_rep_, // repulsion coeficcient
221  k_attr_; // attraction coefficient
222 
223  // cache for reusable intermediate calculations of evaluate_infex
224  mutable EvaluationCache cache_;
225 
226  protected:
227  //! evaluates the linear interaction score, for spheres s0 and s1,
228  //! outputting derivatives to ds0 and ds1 scaled by da, if da is not null
229  inline double evaluate_index
230  (algebra::Sphere3D const& s0,
231  algebra::Sphere3D const& s1,
232  algebra::Sphere3D& ds0,
233  algebra::Sphere3D& ds1,
234  DerivativeAccumulator *da) const;
235 
236  public:
237  /**
238  The score is 0 if the spheres are beyond the attractive range.
239  Within the attractive range, the score decreases linearly (= attraction)
240  with slope k_attr_ until the spheres touch. Once the spheres begin to
241  penetrate each other, the score rises linearly with slope k_rep_
242  (= repulsion), though it may be negative for small penetration.
243 
244  The energy potential difference between the unbound state and when the
245  beads touch is:
246  DELTA-U=0.5*k_attr*range_attr [kcal/mol]
247 
248  @param k_rep repulsion force constant in kcal/mol/A units
249  @param range_attr attraction range between sphere surfaces
250  @param k_attr attraction force constant in kcal/mol/A unit
251  @param name score object names
252  */
253  LinearInteractionPairScore(double k_rep, double range_attr, double k_attr,
254  std::string name = "LinearIDPairScore%1%");
255 
256 #ifndef SWIG
257  // returns cached intermediate computations from last call to
258  // this->evaluate_index()
259  EvaluationCache const &get_evaluation_cache() const { return cache_; }
260 #endif
261 
263  DerivativeAccumulator *da) const
264  override);
265  IMP_IMPLEMENT_INLINE(double evaluate_if_good_index
266  ( Model *m,
267  const ParticleIndexPair &p,
269  double max ) const,
270  {
271  IMP_UNUSED(max);
272  return evaluate_index(m, p, da);
273  }
274  );
275 
276  //! evaluate all index pairs between pips[lower_bound] and pis[upper_bound]
277  //! and update their coordinate derivatives scaled by da, if da is not null
278  virtual double evaluate_indexes(Model *m,
279  const ParticleIndexPairs &pips,
281  unsigned int lower_bound,
282  unsigned int upper_bound) const override;
283 
284  double evaluate_if_good_index(Model *m, const ParticleIndexPairs &p,
285  DerivativeAccumulator *da, double max,
286  unsigned int lower_bound,
287  unsigned int upper_bound) const {
288  double ret = 0;
289  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
290  ret += evaluate_if_good_index(m, p[i], da, max - ret);
291  if (ret > max) return std::numeric_limits<double>::max();
292  }
293  return ret;
294  }
295 
296  ModelObjectsTemp do_get_inputs(Model *m, const ParticleIndexes &pis) const
297  override;
298 
299  //! returns the range for sphere-sphere attraction in A
300  double get_range_attraction() const { return range_attr_; }
301 
302  //! returns the k for sphere-sphere attraction in kcal/mol/A units
303  double get_k_attraction() const { return k_attr_; }
304 
305  //! set k for sphere-sphere attraction in kcal/mol/A units
306  void set_k_attraction(double k_attr) {
307  k_attr_= k_attr;
308  }
309 
310  //! returns the k for sphere-sphere repulsion in kcal/mol/A units
311  double get_k_repulsion() const { return k_rep_; }
312 
313  //! set k for sphere-sphere repulsion in kcal/mol/A units
314  void set_k_repulsion(double k_rep) {
315  k_rep_ = k_rep;
316  }
317 
319 };
320 
321 #ifndef IMP_DOXYGEN
322 
323 //! evaluates the linear interaction score, for spheres s0 and s1,
324 //! outputting derivatives to ds0 and ds1 scaled by da, if da is not null
325 inline double
326 LinearInteractionPairScore::evaluate_index
327 (algebra::Sphere3D const& s0,
328  algebra::Sphere3D const& s1,
329  algebra::Sphere3D& ds0,
330  algebra::Sphere3D& ds1,
331  DerivativeAccumulator *da) const
332 {
333  // Associate intermediate variables with cache_, for further reuse:
334  double &delta_length_2 = cache_.particles_delta_squared;
335  double &x0 = cache_.sum_particles_radii;
336  algebra::Vector3D delta =
337  s0.get_center() - s1.get_center();
338  delta_length_2 = delta.get_squared_magnitude();
339  IMP_LOG(PROGRESS,
340  "LinearInteractionPairScore cached delta2 "
341  << cache_.particles_delta_squared << std::endl);
342  x0 = s0.get_radius() + s1.get_radius();
343  // Terminate immediately if very far, work with squares for speed
344  // equivalent to [delta_length > x0 + attr_range]:
345  if (delta_length_2 > std::pow(x0 + range_attr_, 2)) return 0;
346  double offset = -range_attr_ * k_attr_;
347  double delta_length = std::sqrt(delta_length_2);
348  if (delta_length > x0) { // attractive regime
349  return // decreases with slope k_attr_ as spheres get closer
350  (do_evaluate_index(ds0, ds1, da,
351  delta, delta_length, x0, k_attr_) +
352  offset);
353  } else { // repulsive regime, may be negative for small penetration
354  return (do_evaluate_index(ds0, ds1, da,
355  delta, delta_length, x0, -k_rep_) +
356  offset);
357  }
358 }
359 
360 inline double
361 LinearInteractionPairScore::evaluate_index
362 ( Model *m,
363  const ParticleIndexPair &pp,
364  DerivativeAccumulator *da) const
365 {
367  algebra::Sphere3D const& s0(m->get_sphere(pp[0]));
368  algebra::Sphere3D const& s1(m->get_sphere(pp[1]));
369  algebra::Sphere3D* d_xyzrs=
370  m->access_sphere_derivatives_data();
371  algebra::Sphere3D& ds0( d_xyzrs[pp[0].get_index()] );
372  algebra::Sphere3D& ds1( d_xyzrs[pp[1].get_index()] );
373  return evaluate_index(s0, s1, ds0, ds1, da);
374  }
375 
376 inline double
377 LinearInteractionPairScore::evaluate_indexes
378 ( Model *m,
379  const ParticleIndexPairs &pips,
380  DerivativeAccumulator *da,
381  unsigned int lower_bound,
382  unsigned int upper_bound) const
383 {
385  algebra::Sphere3D const* xyzrs=
386  m->access_spheres_data();
387  algebra::Sphere3D* d_xyzrs=
388  m->access_sphere_derivatives_data();
389  double ret = 0;
390  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
391  int i0(pips[i][0].get_index());
392  int i1(pips[i][1].get_index());
393  ret += evaluate_index( xyzrs[i0], xyzrs[i1],
394  d_xyzrs[i0], d_xyzrs[i1],
395  da);
396  }
397  return ret;
398 }
399 
400 
401 #endif
402 
403 /**
404  A soft linear attractive / repulsive score between two spheres
405  for the backbone of a chain.
406 */
407 class IMPNPCTRANSPORTEXPORT LinearWellPairScore : public PairScore {
408  private:
409  double rest_length_factor_, k_;
410 
411  public:
412  /**
413  a linear well pair potential that keeps two particles around
414  a resting distance relative to their radii
415 
416  @param rest_length_factor the resting distance between particles
417  relative to their sum of radii
418  @param k the force constant (attractive if beyond or repulsive
419  if below rest length) in units kcal/mol/A
420  @param name the name of the score
421  */
422  LinearWellPairScore(double rest_length_factor, double k,
423  std::string name = "LinearIDPairScore%1%");
424 
425  void set_rest_length_factor(double rest_length_factor)
426  { rest_length_factor_ = rest_length_factor; }
427  double get_rest_length_factor() const { return rest_length_factor_; }
428  void set_k(double k)
429  { k_ = k; }
430  double get_k() { return k_; }
431  double evaluate_index(Model *m, const ParticleIndexPair &p,
432  DerivativeAccumulator *da) const override;
434  const ParticleIndexes &pis) const override;
437  ;
438 };
439 
440 #ifndef IMP_DOXYGEN
441 inline double
442 LinearWellPairScore::evaluate_index
443 ( Model *m,
444  const ParticleIndexPair &pp,
445  DerivativeAccumulator *da) const
446 {
448  algebra::Sphere3D const& s0 = m->get_sphere(pp[0]);
449  algebra::Sphere3D const& s1 = m->get_sphere(pp[1]);
450  algebra::Sphere3D* d_xyzrs=
451  m->access_sphere_derivatives_data();
452  algebra::Sphere3D& ds0( d_xyzrs[pp[0].get_index()] );
453  algebra::Sphere3D& ds1( d_xyzrs[pp[1].get_index()] );
454  double x0 = (s0.get_radius() + s1.get_radius()) * rest_length_factor_;
455  algebra::Vector3D delta = s0.get_center() - s1.get_center();
456  double delta_length_2 = delta.get_squared_magnitude();
457  double delta_length = std::sqrt(delta_length_2);
458  if (delta_length > x0) { // attractive regime
459  return // k_ > 0 = get spheres closer
460  do_evaluate_index(ds0, ds1, da,
461  delta, delta_length, x0, k_);
462  } else {
463  return// -k_ < 0 = keep spheres apart
464  do_evaluate_index(ds0, ds1, da,
465  delta, delta_length, x0, -k_);
466  }
467 }
468 #endif
469 
471 
472 
473 IMPNPCTRANSPORT_END_NAMESPACE
474 
475 #endif /* IMPNPCTRANSPORT_LINEAR_DISTANCE_PAIR_SCORES_H */
IMP::Vector< ParticleIndexPair, std::allocator< ParticleIndexPair > > ParticleIndexPairs
Definition: base_types.h:186
#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:44
SphereD< 3 > Sphere3D
Typedef for Python.
Definition: SphereD.h:104
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
#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:106
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:86
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:408
void set_k_repulsion(double k_rep)
set k for sphere-sphere repulsion in kcal/mol/A units
Output a line or two per evaluation call.
Definition: enums.h:31
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.
Class for adding derivatives from restraints to the model.
double get_k() const
returns the k for sphere-sphere repulsion