IMP logo
IMP Reference Guide  develop.36b5e40f34,2024/10/14
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) const override final;
83 
84  /**
85  Evaluate score for particles pis[lower_bound..upper_bound] in
86  model m. If da is not null, use it to accumulate derivatives
87  in model.
88 
89  Abort early and return maximal double value if score is larger
90  than max (note this assumes all score components are
91  non-negative, which is true for slab score)
92  */
94  ( Model *m,
95  const ParticleIndexPairs &pips,
97  double max,
98  unsigned int lower_bound,
99  unsigned int upper_bound) const override
100  {
101  double ret = 0;
102  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
103  ret += evaluate_if_good_index(m, pips[i], da, max - ret);
104  if (ret > max) return std::numeric_limits<double>::max();
105  }
106  return ret;
107  }
108 
109  // IMP_PAIR_SCORE_METHODS(SlabWithCylindricalPorePairScore);
111 
112  private:
113  // evaluate slab for specified sphere, based on most recent cached
114  // slab parameters (/see update_cached_slab_params()).
115  // Return 0 if ball does not penetrate slab.
116  //
117  // @param s the sphere to evaluate
118  // @param out_displacement if not null and the returned score is positive, *out_displacement
119  // is used to store the computed displacement vector from
120  // the surface of the z-axis aligned cylinder to
121  // the center of s. Ignore if score is zero.
122  //
123  inline double evaluate_sphere
125  algebra::Vector3D* out_displacement) const;
126 
127 
128  // computes the displacement from the surface of the z-axis
129  // aligned cylinder to v, based on the most recent cached
130  // slab parameters (/see update_cached_slab_params()).
131  //
132  // @return <distance, a vector pointing out>,
133  // negative distance means v is inside cylinder
134  std::pair<double, algebra::Vector3D> get_displacement_vector(
135  const algebra::Vector3D &v) const;
136 
137  // update internal variables holding slab params for fast access
138  // based on decorated particle slab
139  void update_cached_slab_params
140  (SlabWithCylindricalPore slab) const;
141 };
142 
143 inline void
144 SlabWithCylindricalPorePairScore::update_cached_slab_params
145 (SlabWithCylindricalPore slab) const
146 {
147  //TODO: support slabs with non-zero x,y,z origin
148  thickness_= slab.get_thickness();
149  top_= 0.5*thickness_;
150  bottom_= -0.5*thickness_;
151  midZ_= 0;
152  pore_radius_= slab.get_pore_radius();
153  is_pore_radius_optimized_= slab.get_pore_radius_is_optimized();
154 }
155 
156 //
157 inline double
158 SlabWithCylindricalPorePairScore::evaluate_index
159 (Model *m,
160  const ParticleIndexPair& pip,
161  DerivativeAccumulator *da) const
162 {
164  IMP_USAGE_CHECK(SlabWithCylindricalPore::get_is_setup(m, pip[0]),
165  "pip[0] is not a SlabWithCylindricalPore in evaluate_index()");
166  SlabWithCylindricalPore slab(m, pip[0]);
167  update_cached_slab_params(slab);
168  IMP::core::XYZR d(m, pip[1]);
169  algebra::Sphere3D d_sphere( d.get_sphere() );
171  return false;
172  algebra::Vector3D displacement; // a unit displacement vector - output of evaluate sphere
173  double score=evaluate_sphere(d_sphere,
174  da ? &displacement : nullptr);
175  if(da && score>0.0){
176  algebra::Vector3D derivative_vector = -k_*displacement;
177  IMP_LOG(PROGRESS, "result in " << score << " and " << derivative_vector << std::endl);
178  d.add_to_derivatives(derivative_vector, *da);
179  if(is_pore_radius_optimized_){
180  // TODO: assume that the direction of a positive radial displacement vector is opposite to the sphere x,y vector - is this always true?
181  double radial_displacement_magnitude= // magnitude of the displacement vector projected on the x,y plane
182  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
183  slab.add_to_pore_radius_derivative(-k_*radial_displacement_magnitude, *da);
184  }
185  }
186  return score;
187 }
188 
189 //
190 inline double
191 SlabWithCylindricalPorePairScore::evaluate_indexes
192 (Model *m,
193  const ParticleIndexPairs &pips,
195  unsigned int lower_bound,
196  unsigned int upper_bound) const
197 {
198  if(upper_bound<lower_bound){
199  return 0.0;
200  }
201  double ret(0.0);
202  double radial_displacements_magnitude(0.0); // sum of pore radius displacemnets
203  algebra::Sphere3D const* spheres_table=
204  m->access_spheres_data();
205  algebra::Sphere3D* sphere_derivatives_table=
206  m->access_sphere_derivatives_data();
207  IMP::internal::BoolAttributeTableTraits::Container const& is_optimizable_table=
208  m->access_optimizeds_data(core::XYZ::get_coordinate_key(0)); // use only x coordinate as indicator
209  ParticleIndex slab_pi(pips[lower_bound][0]);
210  SlabWithCylindricalPore slab(m, slab_pi); // TODO: do this only in first round
211  update_cached_slab_params(slab);
212  // Evaluate and sum score and derivative for all particles:
213  for (unsigned int i = lower_bound; i < upper_bound; ++i) {
214  ParticleIndex pi( pips[i][1]);
215  int pi_index=pi.get_index();
216  // Check attributes have valid values:
217  IMP_CHECK_CODE( {
218  IMP_INTERNAL_CHECK(pips[i][0]==slab_pi,
219  "All particles are assumed to be evaluated against"
220  " the same slab");
221  IMP::core::XYZR d(m, pi);
222  algebra::Sphere3D s=spheres_table[pi_index];
223  IMP_INTERNAL_CHECK(d.get_coordinates_are_optimized() == is_optimizable_table[pi],
224  "optimable table inconsistent with d.get_coordinates_are_optimized for particle " << d);
225  IMP_INTERNAL_CHECK((d.get_coordinates() - s.get_center()).get_magnitude()<.001,
226  "Different coords for particle " << d << " *** "
227  << d.get_coordinates() << " vs. " << s.get_center());
228  IMP_INTERNAL_CHECK(d.get_radius() == s.get_radius(),
229  "Different radii for particle " << d << " *** "
230  << d.get_radius() << " vs. " << s.get_radius());
231  } ); // IMP_CHECK_CODE
232  if(!is_optimizable_table[pi]) {
233  continue;
234  }
235  algebra::Vector3D displacement;
236  double cur_score = evaluate_sphere(spheres_table[pi_index],
237  da ? &displacement : nullptr);
238  ret+= cur_score;
239  if(cur_score>0.0 && da) {
240  algebra::Vector3D derivative_vector = -k_*displacement;
241  // accumulate derivatives directly for speed
242  for(unsigned int j=0; j<3; j++) {
243  sphere_derivatives_table[pi_index][j] += (*da)(derivative_vector[j]);
244  }
245  // TODO: assume that the direction of a positive radial displacement vector is opposite to the sphere x,y vector - is this always true?
246  radial_displacements_magnitude+=
247  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
248  }
249  }
250  if(da && is_pore_radius_optimized_){
251  slab.add_to_pore_radius_derivative(-k_ * radial_displacements_magnitude, *da);
252  }
253  return ret;
254 }
255 
256 //
257 double
258 SlabWithCylindricalPorePairScore::evaluate_sphere
260  algebra::Vector3D* out_displacement) const
261 {
263  IMP_LOG(VERBOSE, "evaluate_sphere " << s << std::endl);
264  double const x= s[0];
265  double const y= s[1];
266  double const z= s[2];
267  double const sr= s.get_radius();
268  // early abort if above or below slab
269  if ((z-sr > top_) || (z+sr < bottom_)) {
270  return 0;
271  }
272  double const x2= x*x;
273  double const y2= y*y;
274  double const R= pore_radius_-sr;
275  double const R2= R*R;
276  // early abort if [x,y] within cylinder perimeter
277  if (x2+y2 < R2) {
278  return 0;
279  }
280  std::pair<double, algebra::Vector3D> dp = get_displacement_vector(s.get_center());
281  IMP_LOG(PROGRESS,
282  "At point " << s.get_center() << " have distance " << dp.first
283  << " and direction " << dp.second << std::endl);
284  double const distance = dp.first;
285  if (distance > sr) {
286  return 0;
287  }
288  double const score= k_ * (sr - distance); // must be positive if distance <= sr
289  if(out_displacement){
290  *out_displacement= dp.second;
291  IMP_INTERNAL_CHECK(std::abs(out_displacement->get_magnitude() - 1) < .1,
292  "Not a unit vector");
293  }
294  return score;
295 }
296 
297 // computes the distance and a unit displacement vector of v
298 // from the surface of a z-axis aligned cylinder
299 //
300 // @return <distance, a unit vector pointing outwards>,
301 // negative distance should mean v is inside the cylinder
302 inline std::pair<double, algebra::Vector3D>
303 SlabWithCylindricalPorePairScore::get_displacement_vector(const algebra::Vector3D &v) const {
304  double dXY2 = square(v[0]) + square(v[1]); // r^2 for [x,y] projection
305  double dZ = v[2] - midZ_; // thickness on z-axis from cyl origin
306  IMP_LOG_PROGRESS( dZ << " " << dXY2 << " for " << v << std::endl);
307  if (dXY2 > square(pore_radius_) ||
308  (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
309  double abs_dZ = std::abs(dZ);
310  double abs_dXY = std::sqrt(dXY2);
311  double dR = abs_dXY - pore_radius_; // displacement on [x,y] direction (positive = outside pore perimeter on x,y plane)
312  if (dR + abs_dZ < .5 * thickness_) { // in a cones from (0,0,.5thickness) to (0,.5thicness,0)
313  IMP_LOG_PROGRESS("ring or pore" << std::endl);
314  if (dXY2 < .00001) { // at origin
315  return std::make_pair(pore_radius_, algebra::Vector3D(0, 0, 1));
316  } else {
317  algebra::Vector3D rv(-v[0], -v[1], 0);
318  return std::make_pair(pore_radius_ - abs_dXY, rv.get_unit_vector());
319  }
320  } else { // possibly in the pore but out of the 'double-cone diamond'
321  IMP_LOG_PROGRESS("in or out of slab" << std::endl);
322  if (dZ > 0) {
323  return std::make_pair(v[2] - top_, algebra::Vector3D(0, 0, 1));
324  } else {
325  return std::make_pair(bottom_ - v[2], algebra::Vector3D(0, 0, -1));
326  }
327  }
328  } else { // = outside slab boundaries AND insider pore perimeter on x,y plane
329  IMP_LOG_PROGRESS("channel" << std::endl);
330  if (dXY2 < .00001) { // at central axis
331  IMP_LOG_PROGRESS("in center " << std::endl);
332  if (dZ > 0) {
333  return std::make_pair(v[2] - top_, algebra::Vector3D(0, 0, 1));
334  } else {
335  return std::make_pair(bottom_ - v[2], algebra::Vector3D(0, 0, -1));
336  }
337  }
338  algebra::Vector3D rim =
339  algebra::Vector3D(v[0], v[1], 0).get_unit_vector() * pore_radius_;
340  if (dZ > 0) {
341  rim[2] = top_;
342  } else {
343  rim[2] = bottom_;
344  }
345  IMP_LOG_PROGRESS( "rim is " << rim << std::endl);
346  algebra::Vector3D diff = v - rim;
347  return std::make_pair(diff.get_magnitude(), diff.get_unit_vector());
348  }
349 }
350 
351 
352 
353 
354 IMPNPCTRANSPORT_END_NAMESPACE
355 
356 #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
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
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.
Float get_pore_radius() const
get cylindrical pore radius
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
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: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