IMP logo
IMP Reference Guide  2.14.0
The Integrative Modeling Platform
HarmonicSpringSingletonScore.h
1 /**
2  * \file HarmonicSpringSingletonScoren.h
3  * \brief Harmonic scores on the spring between a pair of particles.
4  * The spring resting length is dynamic
5  * Copyright 2007-2020 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPNPCTRANSPORT_HARMONIC_SPRING_SINGLETON_SCORE_H
9 #define IMPNPCTRANSPORT_HARMONIC_SPRING_SINGLETON_SCORE_H
10 
11 #include "npctransport_config.h"
12 #include "RelaxingSpring.h"
13 #include <IMP/compiler_macros.h>
14 #include <IMP/SingletonScore.h>
15 #include <IMP/pair_macros.h>
16 #include <IMP/core/XYZR.h>
17 #include <IMP/algebra/utility.h>
18 
19 #include <boost/array.hpp>
20 
21 IMPNPCTRANSPORT_BEGIN_NAMESPACE
22 
23 
24 /**
25  A harmonic score between two bonded spheres
26  connected by a spring whose resting length is dynamic
27 */
28 class IMPNPCTRANSPORTEXPORT HarmonicSpringSingletonScore
29 : public SingletonScore
30 {
31  private:
32  double k1_; // see ctr
33  double k2_; // see ctr
34 
35  public:
36  /**
37  a harmonic well pair potential that keeps two particles connected
38  by a spring particle whose resting length itself is dynamic
39 
40  @param k1 the force constant applied by the string on tethered
41  particles in units of kcal/mol/A^2 - attractive if spring
42  is longer than rest length, repulsive if shorter than rest
43  length. The spring applies a force of k1*dX on each particle,
44  and a counter force of 2*k1*dX acts on the spring, where
45  dX=|distance(p1,p2)-(rest_length)|
46  @param k2 the force constant for relaxation of the spring
47  rest length to its equilibrium value in kcal/mol/A^2.
48  The force equals k2*dX where dX=|(rest_length)-(equilibrium length)|
49  @param name the name of the score
50  */
52  ( double k1,
53  double k2,
54  std::string name = "HarmonicSpringSingletonScore%1%");
55 
56  //! sets the force constant to bring particles to the current rest length
57  //! in kcal/mol/A^2
58  void set_k1(double k1)
59  { k1_ = k1; }
60 
61  //! returns the force constant to bring particles to the current rest length
62  //! in kcal/mol/A^2
63  double get_k1() const
64  { return k1_; }
65 
66  //! sets the force constant to bring the current rest length to equilibrium
67  //! in kcal/mol/A^2
68  void set_k2(double k2)
69  { k2_ = k2; }
70 
71  //! returns the force constant to bring the current rest length to equilibrium
72  //! in kcal/mol/A^2
73  double get_k2() const
74  { return k2_; }
75 
76  virtual double
78  (Model *m,
79  ParticleIndex pi,
81 
82  virtual ModelObjectsTemp
83  do_get_inputs(Model *m, const ParticleIndexes &pis) const IMP_OVERRIDE;
84 
86 
88  ;
89 };
90 
91 #ifndef IMP_DOXYGEN
92 
93 inline double
94 HarmonicSpringSingletonScore
95 ::evaluate_index
96 ( Model *m, ParticleIndex pi, DerivativeAccumulator *da) const
97 {
99 
100  IMP_USAGE_CHECK(RelaxingSpring::get_is_setup(m,pi),
101  "particle 0 is expected to be string in HarmonicSpringSingletonScore");
102  RelaxingSpring s(m, pi);
103  ParticleIndex pi0(s.get_bonded_particle_index_0());
104  ParticleIndex pi1(s.get_bonded_particle_index_1());
105  algebra::Sphere3D s0 = m->get_sphere(pi0);
106  algebra::Sphere3D s1 = m->get_sphere(pi1);
107  double rest_delta_length= s.get_rest_length();
108  algebra::Vector3D delta_1_to_0 = s0.get_center() - s1.get_center();
109  double delta_length_2 = delta_1_to_0.get_squared_magnitude();
110  double delta_length = std::sqrt(delta_length_2);
111  double dDelta = delta_length - rest_delta_length; // positive if spring is extended, negative if compressed
112  double scoreDelta = 2 * 0.5 * k1_ * dDelta * dDelta; // x2 because each particle applies force on the spring independently in this model of a relaxing spring
113  double eq_rest_length= s.get_equilibrium_rest_length_factor()
114  * (s0.get_radius() + s1.get_radius());
115  double dEq= rest_delta_length - eq_rest_length; // positive is rest length is stretched relative to equilibrium rest length
116  double scoreEq = 0.5 * k2_ * dEq * dEq;
117  bool is_tiny_rest_length= (rest_delta_length<0.1*eq_rest_length && rest_delta_length<1.0);
118  if(IMP_UNLIKELY(is_tiny_rest_length)) {
119  double threshold=std::min(0.1*eq_rest_length, 1.0);
120  double dThreshold= threshold-rest_delta_length;
121  dEq+= std::pow(10.0 * k2_ * dThreshold / threshold, 4);
122  }
123  double score= scoreDelta + scoreEq;
124  IMP_LOG(TERSE, "dDelta: " << dDelta << " scoreDelta: " << scoreDelta
125  << " dEq: " << dEq << " scoreEq: " << scoreEq
126  << " total: " << score);
127 
128  // Compute derivatives:
129  static const double MIN_DISTANCE = .00001;
130  if (IMP_LIKELY( da && delta_length > MIN_DISTANCE )) { // Profiling note on use of likely(): in BD simulations, the simulation bottleneck is when da is true, and the spring is likely out of equilibrium
131  double fParticles= k1_*dDelta; // force pulling particles closer together for a positive force (or apart if negative)
132  double fSpring= k2_*dEq - 2*fParticles; // force pulling rest length down for a positive force (or up for negtive); the -2*f1 is the coutnerforce exerted on the spring by the tnwo tethered particles
133  if(IMP_UNLIKELY(is_tiny_rest_length)) {
134  double threshold=std::min(0.1*eq_rest_length, 1.0);
135  double dThreshold= threshold-rest_delta_length;
136  dEq+= 40 * k2_ * std::pow(10.0 * k2_ * dThreshold / threshold, 3); // f(x)=(10*k*x)^4; f'(x)= 40*k*(10*k*x)^3
137  }
138  s.add_to_rest_length_derivative(fSpring, *da);
139  algebra::Vector3D deriv0( delta_1_to_0 * (fParticles / delta_length) );
140  m->add_to_coordinate_derivatives(pi0, deriv0, *da);
141  m->add_to_coordinate_derivatives(pi1, -deriv0, *da);
142  IMP_LOG(TERSE, "\nderiv on pi0: " << deriv0);
143  IMP_LOG(TERSE, "\nderiv on spring: " << -fSpring);
144  }
145  IMP_LOG(TERSE, std::endl);
146 
147  return score;
148 }
149 
150 #endif
151 
153 
154 
155 IMPNPCTRANSPORT_END_NAMESPACE
156 
157 #endif /* IMPNPCTRANSPORT_HARMONIC_SPRING_SINGLETON_SCORE_H */
a spring between two diffusing particles whose resting length is relaxing towards some equilibrium va...
Macros for various classes.
#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
A more IMP-like version of the std::vector.
Definition: Vector.h:39
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Functions to deal with very common math operations.
Abstract class for scoring object(s) of type ParticleIndex.
#define IMP_SINGLETON_SCORE_METHODS(Name)
virtual ModelObjectsTemp do_get_inputs(Model *m, const ParticleIndexes &pis) const =0
Overload this method to specify the inputs.
Define SingletonScore.
A decorator for a spring particle connecting two diffusing particles.
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing lists of object pointers.
Definition: object_macros.h:44
VectorD< 3 > Vector3D
Definition: VectorD.h:421
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
virtual double evaluate_index(Model *m, ParticleIndex vt, DerivativeAccumulator *da) const =0
Compute the score and the derivative if needed.
Various compiler workarounds.
Decorator for a sphere-like particle.
#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.