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