IMP logo
IMP Reference Guide  2.20.0
The Integrative Modeling Platform
SlabWithToroidalPorePairScore.h
Go to the documentation of this file.
1 /**
2  * \file IMP/npctransport/SlabWithToroidalPorePairScore.h
3  * \brief a score for a slab with a toroidal pore
4  *
5 
6  * Copyright 2007-2022 IMP Inventors. All rights reserved.
7  *
8  */
9 
10 #ifndef IMPNPCTRANSPORT_SLAB_WITH_TOROIDAL_PORE_PAIR_SCORE_H
11 #define IMPNPCTRANSPORT_SLAB_WITH_TOROIDAL_PORE_PAIR_SCORE_H
12 
13 #include "npctransport_config.h"
14 #include "SlabWithToroidalPore.h"
15 #include <IMP/Model.h>
16 #include <IMP/Pointer.h>
17 #include <IMP/check_macros.h>
18 #include <IMP/PairScore.h>
19 #include <IMP/pair_macros.h>
20 #include "IMP/core/XYZR.h"
21 
22 IMPNPCTRANSPORT_BEGIN_NAMESPACE
23 
24 //! Score for a slab with a toroidal pore
25 class IMPNPCTRANSPORTEXPORT SlabWithToroidalPorePairScore
26 : public PairScore
27 {
28  double k_; // repulsion constant in kcal/mol/A
29 
30  // cached variables for currently evaluated slab (therefore, mutable):
31  mutable double top_; // top z - for caching
32  mutable double bottom_; // bottom z - for caching
33  mutable double midZ_; // vertical center of slab
34  mutable double R_; // torus major radius
35  mutable double rh_; // horizontal minor radius (ellipse horizontal semi-axis)
36  mutable double rv_; // vertical minor radius (ellipse vertical semi-axis)
37  mutable bool is_pore_radius_optimized_;
38 
39 public:
40  //! Constructs a horizontal slab with a toroidal pore,
41  //! centered at the z=0 plane
42  /**
43  Constructs a score over a horizontal slab with a ring toroidal pore
44 
45  @param k the slab repulsive force constant in kcal/mol/A
46  */
48  (double k);
49 
50  public:
51  //! evaluate score for particle pi in model m. If da is not null,
52  //! use it to accumulate derivatives in model.
53  virtual double evaluate_index
54  (Model *m,
55  ParticleIndexPair const& pip,
56  DerivativeAccumulator *da) const override;
57 
59  const ParticleIndexes &pis) const
60  override;
61 
62  //! evaluate score for particles pis[lower_bound..upper_bound] in
63  //! model m. If da is not null, use it to accumulate derivatives
64  //! in model.
65  virtual double evaluate_indexes
66  (Model *m,
67  const ParticleIndexPairs &pips,
69  unsigned int lower_bound,
70  unsigned int upper_bound) const override final;
71 
72  /**
73  Evaluate score for particles pis[lower_bound..upper_bound] in
74  model m. If da is not null, use it to accumulate derivatives
75  in model.
76 
77  Abort early and return maximal double value if score is larger
78  than max (note this assumes all score components are
79  non-negative, which is true for slab score)
80  */
82  ( Model *m,
83  const ParticleIndexPairs &p,
85  double max,
86  unsigned int lower_bound,
87  unsigned int upper_bound) const override
88  {
89  double ret = 0;
90  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
91  ret += evaluate_if_good_index(m, p[i], da, max - ret);
92  if (ret > max) return std::numeric_limits<double>::max();
93  }
94  return ret;
95  }
96 
98 
99  private:
100 
101  /**
102  Computes the penetration depth between the specified sphere and an
103  axis aligned ellipsoid defined by the torus vertical minor radius and
104  horizontal minor radius, centered at the specified origin (based on
105  cached slab parameters)
106  Optionally stores a normalized translation vector in *out_translation
107 
108  @param sphere Sphere for which the distance is computed
109  @param origin Ellipsoid origin (a point on the major circle of the torus)
110  @param out_translation A pointer to the normalized translation
111  vector for removing the sphere out of the ellipsoid is
112  stored in *out_translation, if it is not a null pointer.
113  (0,0,0) is stored for 0 overlap.
114 
115  @return The penetration depth or 0 when the sphere does not overlap with
116  the ellipsoid.
117  */
118  inline double
119  get_sphere_ellipsoid_penetration_depth
120  (algebra::Sphere3D const& sphere,
121  algebra::Vector3D const& origin,
122  algebra::Vector3D* out_translation) const;
123 
124  // returns the penetration depth D of a sphere relative to the porus slab surface,
125  // or 0 if the sphere and the slab do not overlap (based on cached slab
126  // parameters)
127  //
128  // @param s the sphere to evaluate
129  // @param out_translation if not null, *out_translation
130  // is used to store the output normalized displacement vector for removing
131  // the sphere out of the slab (or 0,0,0 if no overlap)
132  // In other words, D*out_translation be the shortest displacement vector
133  // needed to remove the sphere out of the slab.
134  inline double get_sphere_penetration_depth
136  algebra::Vector3D* out_translation) const;
137 
138  // update internal variables holding slab params for fast access
139  // based on decorated particle swp
140  void update_cached_slab_params
141  (SlabWithToroidalPore slab) const;
142 };
143 
144 inline void
145 SlabWithToroidalPorePairScore::update_cached_slab_params
146 (SlabWithToroidalPore slab) const
147 {
148  //TODO: support slabs with non-zero x,y,z origin
149  double thickness= slab.get_thickness();
150  top_= 0.5*thickness;
151  bottom_= -0.5*thickness;
152  midZ_= 0;
153  R_= slab.get_pore_radius(); // major radius
154  rh_= slab.get_horizontal_minor_radius();
155  rv_= slab.get_vertical_minor_radius();
156  is_pore_radius_optimized_= slab.get_pore_radius_is_optimized();
157 }
158 
159 
160 //
161 inline double
162 SlabWithToroidalPorePairScore::evaluate_index
163 (Model *m,
164  const ParticleIndexPair& pip,
165  DerivativeAccumulator *da) const
166 {
168  SlabWithToroidalPore slab(m, pip[0]);
169  update_cached_slab_params(slab);
170  IMP::core::XYZR d(m, pip[1]);
172  return false;
173  algebra::Sphere3D d_sphere( d.get_sphere() );
174  algebra::Vector3D displacement;
175  double score=get_sphere_penetration_depth(d.get_sphere(),
176  da ? &displacement : nullptr);
177  IMP_LOG_TERSE("sphere " << d << " score " << score);
178  if(da && score>0.0){
179  algebra::Vector3D derivative_vector = -k_*displacement;
180  IMP_LOG_TERSE(" derivative vector " << derivative_vector);
181  d.add_to_derivatives(derivative_vector, *da);
182  if(is_pore_radius_optimized_){
183  double radial_displacement= displacement[0]*d_sphere[0]+displacement[1]*d_sphere[1]; // TODO: assumes slab origin at 0,0,0 - perhaps in the future extend to general case
184  slab.add_to_pore_radius_derivative(k_*radial_displacement, *da);
185  }
186 
187  }
188  return score;
189 }
190 
191 //
192 inline double
193 SlabWithToroidalPorePairScore::evaluate_indexes
194  (Model *m,
195  const ParticleIndexPairs &pips,
197  unsigned int lower_bound,
198  unsigned int upper_bound) const
199 {
200  IMP_LOG_TERSE("SlabWithToroidalPore singleton - evaluate indexes"
201  << std::endl);
202  if(upper_bound<lower_bound){
203  return 0.0;
204  }
205  double ret(0.0);
206  double radial_displacements(0.0); // sum of pore radius displacemnets
207  algebra::Sphere3D const* spheres_table=
208  m->access_spheres_data();
209  algebra::Sphere3D* sphere_derivatives_table=
210  m->access_sphere_derivatives_data();
211  IMP::internal::BoolAttributeTableTraits::Container const&
212  is_optimizable_table= m->access_optimizeds_data
213  (core::XYZ::get_coordinate_key(0)); // x is indicator
214  ParticleIndex slab_pi(pips[lower_bound][0]);
215  SlabWithToroidalPore slab(m, slab_pi); // TODO: do this only in first round
216  update_cached_slab_params(slab);
217 
218  // Evaluate and sum score and derivative for all particles:
219  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
220  ParticleIndex pi( pips[i][1]);
221  int pi_index=pi.get_index();
222  // Check attributes have valid values:
223  IMP_CHECK_CODE( {
224  IMP_INTERNAL_CHECK(pips[i][0]==slab_pi,
225  "All particles are assumed to be evaluated against"
226  " the same slab");
227  IMP::core::XYZR d(m, pi);
228  algebra::Sphere3D s=spheres_table[pi_index];
229  IMP_INTERNAL_CHECK(d.get_coordinates_are_optimized() == is_optimizable_table[pi],
230  "optimable table inconsistent with d.get_coordinates_are_optimized for particle " << d);
231  IMP_INTERNAL_CHECK((d.get_coordinates() - s.get_center()).get_magnitude()<.001,
232  "Different coords for particle " << d << " *** "
233  << d.get_coordinates() << " vs. " << s.get_center());
234  IMP_INTERNAL_CHECK(d.get_radius() == s.get_radius(),
235  "Different radii for particle " << d << " *** "
236  << d.get_radius() << " vs. " << s.get_radius());
237  } ); // IMP_CHECK_CODE
238  if(!is_optimizable_table[pi]) {
239  continue;
240  }
241  algebra::Vector3D displacement;
242  double cur_score = get_sphere_penetration_depth(spheres_table[pi_index],
243  da ? &displacement : nullptr);
244  IMP_LOG_TERSE("SlabWithToroidalPore singleton score for sphere / displacement "
245  << spheres_table[pi_index] << " is " << cur_score << " / "
246  << displacement << std::endl);
247  ret+= cur_score;
248  if(cur_score>0.0 && da) {
249  algebra::Vector3D derivative_vector = -k_*displacement;
250  // accumulate derivatives directly for speed
251  for(unsigned int j=0; j<3; j++) {
252  sphere_derivatives_table[pi_index][j] += (*da)(derivative_vector[j]);
253  }
254  algebra::Sphere3D const& s=spheres_table[pi_index];
255  radial_displacements+= displacement[0]*s[0] + displacement[1]*s[1]; // TODO: assumes slab origin at 0,0,0 - perhaps in דthe future extend to general case
256 
257  }
258  }
259  if(da && is_pore_radius_optimized_){
260  slab.add_to_pore_radius_derivative(k_ * radial_displacements, *da);
261  }
262  return ret;
263 }
264 
265 // sphere - sphere for which the distance is computed
266 // origin - ellipsoid origin (a point on the major circle of the torus)
267 double
268 SlabWithToroidalPorePairScore::
269 get_sphere_ellipsoid_penetration_depth
270 (algebra::Sphere3D const& sphere,
271  algebra::Vector3D const& origin,
272  algebra::Vector3D* out_translation) const
273 {
274  const double eps= 1e-9;
275  algebra::Vector3D v(sphere.get_center()-origin);
276  double dXY2= v[0]*v[0]+v[1]*v[1];
277  double dZ2= v[2]*v[2];
278  double dv2= dXY2 + dZ2 + eps;
279  // theta = atan(dXY/dZ) = asin(dXY/dv) = acos(dZ/dv)
280  double sinTheta2= dXY2/dv2;
281  double cosTheta2= dZ2/dv2;
282  double cur_r= std::sqrt(rv_*rv_*cosTheta2 + rh_*rh_*sinTheta2);
283  double dv=std::sqrt(dv2);
284  double surface_distance= dv - sphere.get_radius() - cur_r;
285  if(surface_distance>=0 ){
286  if(out_translation){
287  (*out_translation)=IMP::algebra::Vector3D(0,0,0);
288  }
289  return 0;
290  }
291  // Overlap:
292  if(out_translation){
293  (*out_translation)= v.get_unit_vector();
294  }
295  return -surface_distance; // penetration is negative the 'distance'
296 
297 
298 }
299 
300 // return - distance to nearest surface point
301 // out_translation - a normalized vector pointing to the nearest surface point
302 double
303 SlabWithToroidalPorePairScore::get_sphere_penetration_depth
304 (algebra::Sphere3D sphere,
305  algebra::Vector3D* out_translation) const
306 {
307  double const x=sphere[0];
308  double const y=sphere[1];
309  double const z=sphere[2];
310  double const sr=sphere.get_radius();
311  double sphere_dz_top= (z - sr) - top_;
312  double sphere_dz_bottom = (z + sr) - bottom_;
313  bool is_above_top=sphere_dz_top>0;
314  bool is_below_bottom=sphere_dz_bottom<0;
315  if(is_above_top || is_below_bottom) {
316  // early abort I - above or below slab
317  if(out_translation){
318  (*out_translation)=IMP::algebra::Vector3D(0,0,0);
319  }
320  return 0;
321  }
322  double d_xy2 = x*x + y*y; // distance from (0,0,z)
323  bool is_closer_to_top= (sphere_dz_bottom > -sphere_dz_top);
324  double abs_sphere_dz =
325  is_closer_to_top ? -sphere_dz_top : sphere_dz_bottom; // d_z is the mianimal vertical change in z that brings s either fully above or fully below the slab
326  if (d_xy2 > R_*R_) {
327  // early about II - radially outside torus major circle of radius R
328  if(out_translation){
329  (*out_translation)=IMP::algebra::Vector3D(0,0,is_closer_to_top?+1:-1);
330  }
331  return abs_sphere_dz;
332  }
333  double d_xy = std::sqrt(d_xy2);
334  // IV. In slab vertically + radially within torus major circle of radius R
335  const double eps = 1e-9;
336  double ret;
337  if (d_xy>eps) {
338  double scale= R_/d_xy;
339  IMP::algebra::Vector3D origin(x*scale, y*scale, midZ_);
340  ret= get_sphere_ellipsoid_penetration_depth
341  (sphere, origin, out_translation );
342  } else {
343  IMP::algebra::Vector3D origin(R_, 0.0, midZ_);
344  ret= get_sphere_ellipsoid_penetration_depth
345  (sphere, origin, out_translation );
346  }
347  return ret;
348 }
349 
350 
351 
352 IMPNPCTRANSPORT_END_NAMESPACE
353 
354 #endif /* IMPNPCTRANSPORT_SLAB_WITH_TOROIDAL_PORE_PAIR_SCORE_H */
void add_to_derivatives(const algebra::Vector3D &v, DerivativeAccumulator &d)
Add the vector v to the derivative vector of the x,y,z coordinates.
Definition: XYZ.h:82
void add_to_pore_radius_derivative(double v, DerivativeAccumulator &d)
Definition: SlabWithPore.h:102
Abstract class for scoring object(s) of type ParticleIndexPair.
Definition: PairScore.h:44
A decorator for a particle that's a slab with a toroidal pore.
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
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
Storage of a model, its restraints, constraints and particles.
bool get_coordinates_are_optimized() const
Get whether the coordinates are optimized.
Definition: XYZ.h:89
A more IMP-like version of the std::vector.
Definition: Vector.h:42
#define IMP_LOG_TERSE(expr)
Definition: log_macros.h:72
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:86
const algebra::Sphere3D & get_sphere() const
Return a sphere object.
Definition: XYZR.h:67
Float get_pore_radius() const
get cylindrical pore radius
Definition: SlabWithPore.h:91
Define PairScore.
virtual ModelObjectsTemp do_get_inputs(Model *m, const ParticleIndexes &pis) const =0
Overload this method to specify the inputs.
const algebra::Vector3D & get_coordinates() const
Convert it to a vector.
Definition: XYZ.h:109
#define IMP_CHECK_CODE(expr)
Only compile the code if checks are enabled.
Definition: check_macros.h:117
A nullptr-initialized pointer to an IMP Object.
Float get_thickness() const
returns whether the particle last entered the transport moiety from its
Definition: SlabWithPore.h:86
Helper macros for throwing and handling exceptions.
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:425
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.
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.
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27