IMP logo
IMP Reference Guide  develop.562c2f4da8,2025/03/10
The Integrative Modeling Platform
SlabWithCylindricalPorePairScore.h
Go to the documentation of this file.
1 /**
2  * \file SlabWithCylindricalPorePairScore.h
3  * \brief XXXXXXXXXXXXXX
4  *
5  * Copyright 2007-2022 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPNPCTRANSPORT_SLAB_PAIR_SCORE_H
9 #define IMPNPCTRANSPORT_SLAB_PAIR_SCORE_H
10 
11 #include "npctransport_config.h"
13 #include <IMP/check_macros.h>
14 #include <IMP/PairScore.h>
15 #include <IMP/pair_macros.h>
16 #include "IMP/core/XYZR.h"
17 
18 IMPNPCTRANSPORT_BEGIN_NAMESPACE
19 
20 //! XXXX
21 /** An origin centered slab with a pore in the vertical direction,
22  for z = [-0.5*thickness_...0.5*thickness_]
23 
24  The score evaluates to 0 for all particles fully beyond z range
25  or fully within slab radius from the origin in the [X,Y] plane
26  For particles that penetrate the slab, the score increases linearly
27  with the magnitude of penetration with a slope k in units of kcal/mol/A.
28  Conseuqently, the gradient is constant and it is oriented so as to repulse
29  the particle towards the nearest point on the slab surface.
30  */
31 class IMPNPCTRANSPORTEXPORT
33  private:
34  double k_; // coefficient for violation of a slab restraint in kcal/mol/A
35 
36  // cache variables (therefore, all are mutable, as they are only used for performance purposes)
37  mutable double thickness_; // thickness of slab
38  mutable double pore_radius_; // radius of slab cylinder
39  mutable double top_; // top of slab on z-axis
40  mutable double bottom_; // bottom of slab on x-axis
41  mutable double midZ_; // (top + bottom) / 2, for caching some calculations
42  mutable bool is_pore_radius_optimized_;
43 
44  public:
45  //! Constructs a slab pair score that acts on a
46  //! SlabWithCylindricalPore particle and a diffusing particle.
47  //! The slab applies a repulsive force constant k specified in units
48  //! of kcal/mol/A (linear potential)
50 
51  //! returns the direction vector for the displacement of point v relative to the slab surface
52  algebra::Vector3D get_displacement_direction
53  (SlabWithCylindricalPore const& slab, const algebra::Vector3D &v) const;
54 
55  //! returns the displacement magnitude of point v relative to the slab surface
56  double get_displacement_magnitude
57  (SlabWithCylindricalPore const&slab, const algebra::Vector3D &v) const;
58 
59  public:
60  //! evaluate score for particle pair pip in model m
61  /** evaluate score for particle pair pip in model m, where the first particle
62  is assumed to be a cylindrical slab. If da is not null,
63  use it to accumulate derivatives in the model.
64  */
65  virtual double evaluate_index
66  (Model *m,
67  const ParticleIndexPair& pip,
68  DerivativeAccumulator *da) const override;
69 
71  const ParticleIndexes &pis) const
72  override;
73 
74  //! evaluate score for particles pis[lower_bound..upper_bound] in
75  //! model m. If da is not null, use it to accumulate derivatives
76  //! in model.
77  virtual double evaluate_indexes
78  (Model *m,
79  const ParticleIndexPairs &pips,
81  unsigned int lower_bound,
82  unsigned int upper_bound,
83  bool all_indexes_checked=false) const override final;
84 
85  /**
86  Evaluate score for particles pis[lower_bound..upper_bound] in
87  model m. If da is not null, use it to accumulate derivatives
88  in model.
89 
90  Abort early and return maximal double value if score is larger
91  than max (note this assumes all score components are
92  non-negative, which is true for slab score)
93  */
95  ( Model *m,
96  const ParticleIndexPairs &pips,
98  double max,
99  unsigned int lower_bound,
100  unsigned int upper_bound,
101  bool all_indexes_checked=false) const override
102  {
103  IMP_UNUSED(all_indexes_checked);
104  double ret = 0;
105  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
106  ret += evaluate_if_good_index(m, pips[i], da, max - ret);
107  if (ret > max) return std::numeric_limits<double>::max();
108  }
109  return ret;
110  }
111 
112  // IMP_PAIR_SCORE_METHODS(SlabWithCylindricalPorePairScore);
114 
115  private:
116  // evaluate slab for specified sphere, based on most recent cached
117  // slab parameters (/see update_cached_slab_params()).
118  // Return 0 if ball does not penetrate slab.
119  //
120  // @param s the sphere to evaluate
121  // @param out_displacement if not null and the returned score is positive, *out_displacement
122  // is used to store the computed displacement vector from
123  // the surface of the z-axis aligned cylinder to
124  // the center of s. Ignore if score is zero.
125  //
126  inline double evaluate_sphere
128  algebra::Vector3D* out_displacement) const;
129 
130 
131  // computes the displacement from the surface of the z-axis
132  // aligned cylinder to v, based on the most recent cached
133  // slab parameters (/see update_cached_slab_params()).
134  //
135  // @return <distance, a vector pointing out>,
136  // negative distance means v is inside cylinder
137  std::pair<double, algebra::Vector3D> get_displacement_vector(
138  const algebra::Vector3D &v) const;
139 
140  // update internal variables holding slab params for fast access
141  // based on decorated particle slab
142  void update_cached_slab_params
143  (SlabWithCylindricalPore slab) const;
144 };
145 
146 inline void
147 SlabWithCylindricalPorePairScore::update_cached_slab_params
148 (SlabWithCylindricalPore slab) const
149 {
150  //TODO: support slabs with non-zero x,y,z origin
151  thickness_= slab.get_thickness();
152  top_= 0.5*thickness_;
153  bottom_= -0.5*thickness_;
154  midZ_= 0;
155  pore_radius_= slab.get_pore_radius();
156  is_pore_radius_optimized_= slab.get_pore_radius_is_optimized();
157 }
158 
159 //
160 inline double
161 SlabWithCylindricalPorePairScore::evaluate_index
162 (Model *m,
163  const ParticleIndexPair& pip,
164  DerivativeAccumulator *da) const
165 {
167  IMP_USAGE_CHECK(SlabWithCylindricalPore::get_is_setup(m, pip[0]),
168  "pip[0] is not a SlabWithCylindricalPore in evaluate_index()");
169  SlabWithCylindricalPore slab(m, pip[0]);
170  update_cached_slab_params(slab);
171  IMP::core::XYZR d(m, pip[1]);
172  algebra::Sphere3D d_sphere( d.get_sphere() );
174  return false;
175  algebra::Vector3D displacement; // a unit displacement vector - output of evaluate sphere
176  double score=evaluate_sphere(d_sphere,
177  da ? &displacement : nullptr);
178  if(da && score>0.0){
179  algebra::Vector3D derivative_vector = -k_*displacement;
180  IMP_LOG(PROGRESS, "result in " << score << " and " << derivative_vector << std::endl);
181  d.add_to_derivatives(derivative_vector, *da);
182  if(is_pore_radius_optimized_){
183  // TODO: assume that the direction of a positive radial displacement vector is opposite to the sphere x,y vector - is this always true?
184  double radial_displacement_magnitude= // magnitude of the displacement vector projected on the x,y plane
185  std::sqrt(displacement[0]*displacement[0]+displacement[1]*displacement[1]); // TODO: currently we assume slab origin at 0,0,0 - perhaps in the future extend to general case
186  slab.add_to_pore_radius_derivative(-k_*radial_displacement_magnitude, *da);
187  }
188  }
189  return score;
190 }
191 
192 //
193 inline double
194 SlabWithCylindricalPorePairScore::evaluate_indexes
195 (Model *m,
196  const ParticleIndexPairs &pips,
198  unsigned int lower_bound,
199  unsigned int upper_bound, bool) const
200 {
201  if(upper_bound<lower_bound){
202  return 0.0;
203  }
204  double ret(0.0);
205  double radial_displacements_magnitude(0.0); // sum of pore radius displacemnets
206  algebra::Sphere3D const* spheres_table=
207  m->access_spheres_data();
208  algebra::Sphere3D* sphere_derivatives_table=
209  m->access_sphere_derivatives_data();
210  IMP::internal::BoolAttributeTableTraits::Container const& is_optimizable_table=
211  m->access_optimizeds_data(core::XYZ::get_coordinate_key(0)); // use only x coordinate as indicator
212  ParticleIndex slab_pi(pips[lower_bound][0]);
213  SlabWithCylindricalPore slab(m, slab_pi); // TODO: do this only in first round
214  update_cached_slab_params(slab);
215  // Evaluate and sum score and derivative for all particles:
216  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
217  ParticleIndex pi( pips[i][1]);
218  int pi_index=pi.get_index();
219  // Check attributes have valid values:
220  IMP_CHECK_CODE( {
221  IMP_INTERNAL_CHECK(pips[i][0]==slab_pi,
222  "All particles are assumed to be evaluated against"
223  " the same slab");
224  IMP::core::XYZR d(m, pi);
225  algebra::Sphere3D s=spheres_table[pi_index];
226  IMP_INTERNAL_CHECK(d.get_coordinates_are_optimized() == is_optimizable_table[pi],
227  "optimable table inconsistent with d.get_coordinates_are_optimized for particle " << d);
228  IMP_INTERNAL_CHECK((d.get_coordinates() - s.get_center()).get_magnitude()<.001,
229  "Different coords for particle " << d << " *** "
230  << d.get_coordinates() << " vs. " << s.get_center());
231  IMP_INTERNAL_CHECK(d.get_radius() == s.get_radius(),
232  "Different radii for particle " << d << " *** "
233  << d.get_radius() << " vs. " << s.get_radius());
234  } ); // IMP_CHECK_CODE
235  if(!is_optimizable_table[pi]) {
236  continue;
237  }
238  algebra::Vector3D displacement;
239  double cur_score = evaluate_sphere(spheres_table[pi_index],
240  da ? &displacement : nullptr);
241  ret+= cur_score;
242  if(cur_score>0.0 && da) {
243  algebra::Vector3D derivative_vector = -k_*displacement;
244  // accumulate derivatives directly for speed
245  for(unsigned int j=0; j<3; j++) {
246  sphere_derivatives_table[pi_index][j] += (*da)(derivative_vector[j]);
247  }
248  // TODO: assume that the direction of a positive radial displacement vector is opposite to the sphere x,y vector - is this always true?
249  radial_displacements_magnitude+=
250  std::sqrt(displacement[0]*displacement[0] + displacement[1]*displacement[1]); // TODO: assumes slab origin at 0,0,0 - perhaps in דthe future extend to general case
251  }
252  }
253  if(da && is_pore_radius_optimized_){
254  slab.add_to_pore_radius_derivative(-k_ * radial_displacements_magnitude, *da);
255  }
256  return ret;
257 }
258 
259 //
260 double
261 SlabWithCylindricalPorePairScore::evaluate_sphere
263  algebra::Vector3D* out_displacement) const
264 {
266  IMP_LOG(VERBOSE, "evaluate_sphere " << s << std::endl);
267  double const x= s[0];
268  double const y= s[1];
269  double const z= s[2];
270  double const sr= s.get_radius();
271  // early abort if above or below slab
272  if ((z-sr > top_) || (z+sr < bottom_)) {
273  return 0;
274  }
275  double const x2= x*x;
276  double const y2= y*y;
277  double const R= pore_radius_-sr;
278  double const R2= R*R;
279  // early abort if [x,y] within cylinder perimeter
280  if (x2+y2 < R2) {
281  return 0;
282  }
283  std::pair<double, algebra::Vector3D> dp = get_displacement_vector(s.get_center());
284  IMP_LOG(PROGRESS,
285  "At point " << s.get_center() << " have distance " << dp.first
286  << " and direction " << dp.second << std::endl);
287  double const distance = dp.first;
288  if (distance > sr) {
289  return 0;
290  }
291  double const score= k_ * (sr - distance); // must be positive if distance <= sr
292  if(out_displacement){
293  *out_displacement= dp.second;
294  IMP_INTERNAL_CHECK(std::abs(out_displacement->get_magnitude() - 1) < .1,
295  "Not a unit vector");
296  }
297  return score;
298 }
299 
300 // computes the distance and a unit displacement vector of v
301 // from the surface of a z-axis aligned cylinder
302 //
303 // @return <distance, a unit vector pointing outwards>,
304 // negative distance should mean v is inside the cylinder
305 inline std::pair<double, algebra::Vector3D>
306 SlabWithCylindricalPorePairScore::get_displacement_vector(const algebra::Vector3D &v) const {
307  double dXY2 = square(v[0]) + square(v[1]); // r^2 for [x,y] projection
308  double dZ = v[2] - midZ_; // thickness on z-axis from cyl origin
309  IMP_LOG_PROGRESS( dZ << " " << dXY2 << " for " << v << std::endl);
310  if (dXY2 > square(pore_radius_) ||
311  (v[2] <= top_ && v[2] >= bottom_)) { // inside vertical slab boundaries and pore perimeter on x,y plane OR outside pore perimeter in any vertical position
312  double abs_dZ = std::abs(dZ);
313  double abs_dXY = std::sqrt(dXY2);
314  double dR = abs_dXY - pore_radius_; // displacement on [x,y] direction (positive = outside pore perimeter on x,y plane)
315  if (dR + abs_dZ < .5 * thickness_) { // in a cones from (0,0,.5thickness) to (0,.5thicness,0)
316  IMP_LOG_PROGRESS("ring or pore" << std::endl);
317  if (dXY2 < .00001) { // at origin
318  return std::make_pair(pore_radius_, algebra::Vector3D(0, 0, 1));
319  } else {
320  algebra::Vector3D rv(-v[0], -v[1], 0);
321  return std::make_pair(pore_radius_ - abs_dXY, rv.get_unit_vector());
322  }
323  } else { // possibly in the pore but out of the 'double-cone diamond'
324  IMP_LOG_PROGRESS("in or out of slab" << std::endl);
325  if (dZ > 0) {
326  return std::make_pair(v[2] - top_, algebra::Vector3D(0, 0, 1));
327  } else {
328  return std::make_pair(bottom_ - v[2], algebra::Vector3D(0, 0, -1));
329  }
330  }
331  } else { // = outside slab boundaries AND insider pore perimeter on x,y plane
332  IMP_LOG_PROGRESS("channel" << std::endl);
333  if (dXY2 < .00001) { // at central axis
334  IMP_LOG_PROGRESS("in center " << std::endl);
335  if (dZ > 0) {
336  return std::make_pair(v[2] - top_, algebra::Vector3D(0, 0, 1));
337  } else {
338  return std::make_pair(bottom_ - v[2], algebra::Vector3D(0, 0, -1));
339  }
340  }
341  algebra::Vector3D rim =
342  algebra::Vector3D(v[0], v[1], 0).get_unit_vector() * pore_radius_;
343  if (dZ > 0) {
344  rim[2] = top_;
345  } else {
346  rim[2] = bottom_;
347  }
348  IMP_LOG_PROGRESS( "rim is " << rim << std::endl);
349  algebra::Vector3D diff = v - rim;
350  return std::make_pair(diff.get_magnitude(), diff.get_unit_vector());
351  }
352 }
353 
354 
355 
356 
357 IMPNPCTRANSPORT_END_NAMESPACE
358 
359 #endif /* IMPNPCTRANSPORT_SLAB_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)
Abstract class for scoring object(s) of type ParticleIndexPair.
Definition: PairScore.h:44
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
bool get_coordinates_are_optimized() const
Get whether the coordinates are optimized.
Definition: XYZ.h:89
#define IMP_LOG_PROGRESS(expr)
Definition: log_macros.h:94
A more IMP-like version of the std::vector.
Definition: Vector.h:50
#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
A decorator for a particle that's a slab with a cylindrical pore.
virtual double evaluate_if_good_indexes(Model *m, const ParticleIndexPairs &o, DerivativeAccumulator *da, double max, unsigned int lower_bound, unsigned int upper_bound, bool all_indexes_checked=false) const
Float get_pore_radius() const
get cylindrical pore radius
#define IMP_UNUSED(variable)
Produce copious output to allow someone to trace through the computation.
Definition: enums.h:33
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
Float get_thickness() const
returns whether the particle last entered the transport moiety from its
virtual double evaluate_indexes(Model *m, const ParticleIndexPairs &o, DerivativeAccumulator *da, unsigned int lower_bound, unsigned int upper_bound, bool all_indexes_checked=false) const
Compute the score and the derivative if needed over a set.
Helper macros for throwing and handling exceptions.
VectorD< 3 > Vector3D
Definition: VectorD.h:408
#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_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