00001
00005 #ifndef CCTBX_ADPTBX_H
00006 #define CCTBX_ADPTBX_H
00007
00008 #include <cctbx/uctbx.h>
00009 #include <scitbx/matrix/eigensystem.h>
00010 #include <scitbx/matrix/tensor_rank_2.h>
00011 #include <scitbx/array_family/tiny_algebra.h>
00012 #include <scitbx/type_holder.h>
00013 #include <cctbx/error.h>
00014
00015 namespace cctbx {
00016
00018 namespace adptbx {
00019
00020 using scitbx::vec3;
00021 using scitbx::mat3;
00022 using scitbx::sym_mat3;
00023
00024 inline void
00025 throw_not_positive_definite()
00026 {
00027 throw error("anisotropic displacement tensor is not positive definite.");
00028 }
00029
00031 inline double
00032 u_as_b(double u_iso)
00033 {
00034 return u_iso * scitbx::constants::eight_pi_sq;
00035 }
00036
00038 inline double
00039 b_as_u(double b_iso)
00040 {
00041 return b_iso / scitbx::constants::eight_pi_sq;
00042 }
00043
00045 template <typename FloatType>
00046 sym_mat3<FloatType>
00047 u_as_b(sym_mat3<FloatType> const& u_aniso)
00048 {
00049 return
00050 scitbx::constants::eight_pi_sq * sym_mat3<FloatType>(u_aniso);
00051 }
00052
00054 template <typename FloatType>
00055 sym_mat3<FloatType>
00056 b_as_u(sym_mat3<FloatType> const& b_aniso)
00057 {
00058 return FloatType(1. / scitbx::constants::eight_pi_sq) * b_aniso;
00059 }
00060
00062
00076 template <typename FloatType>
00077 sym_mat3<FloatType>
00078 u_cif_as_u_star(
00079 uctbx::unit_cell const& unit_cell,
00080 sym_mat3<FloatType> const& u_cif)
00081 {
00082 af::double6 const& r_params = unit_cell.reciprocal_parameters();
00083 return sym_mat3<FloatType>(
00084 u_cif[0] * (r_params[0] * r_params[0]),
00085 u_cif[1] * (r_params[1] * r_params[1]),
00086 u_cif[2] * (r_params[2] * r_params[2]),
00087 u_cif[3] * (r_params[0] * r_params[1]),
00088 u_cif[4] * (r_params[0] * r_params[2]),
00089 u_cif[5] * (r_params[1] * r_params[2]));
00090 }
00091
00093
00095 template <typename FloatType>
00096 sym_mat3<FloatType>
00097 u_star_as_u_cif(
00098 uctbx::unit_cell const& unit_cell,
00099 sym_mat3<FloatType> const& u_star)
00100 {
00101 af::double6 const& r_params = unit_cell.reciprocal_parameters();
00102 return sym_mat3<FloatType>(
00103 u_star[0] / (r_params[0] * r_params[0]),
00104 u_star[1] / (r_params[1] * r_params[1]),
00105 u_star[2] / (r_params[2] * r_params[2]),
00106 u_star[3] / (r_params[0] * r_params[1]),
00107 u_star[4] / (r_params[0] * r_params[2]),
00108 u_star[5] / (r_params[1] * r_params[2]));
00109 }
00110
00112
00116 template <typename FloatType>
00117 inline sym_mat3<FloatType>
00118 u_cart_as_u_star(
00119 uctbx::unit_cell const& unit_cell,
00120 sym_mat3<FloatType> const& u_cart)
00121 {
00122 return u_cart.tensor_transform(unit_cell.fractionalization_matrix());
00123 }
00124
00126
00128 template <typename FloatType>
00129 inline sym_mat3<FloatType>
00130 u_star_as_u_cart(
00131 uctbx::unit_cell const& unit_cell,
00132 sym_mat3<FloatType> const& u_star)
00133 {
00134 return u_star.tensor_transform(unit_cell.orthogonalization_matrix());
00135 }
00136
00138
00141 template <typename FloatType>
00142 inline sym_mat3<FloatType>
00143 u_cart_as_u_cif(
00144 uctbx::unit_cell const& unit_cell,
00145 sym_mat3<FloatType> const& u_cart)
00146 {
00147 return u_star_as_u_cif(unit_cell, u_cart_as_u_star(unit_cell, u_cart));
00148 }
00149
00151
00154 template <typename FloatType>
00155 inline sym_mat3<FloatType>
00156 u_cif_as_u_cart(
00157 uctbx::unit_cell const& unit_cell,
00158 sym_mat3<FloatType> const& u_cif)
00159 {
00160 return u_star_as_u_cart(unit_cell, u_cif_as_u_star(unit_cell, u_cif));
00161 }
00162
00164
00166 template <typename FloatType>
00167 inline sym_mat3<FloatType>
00168 u_star_as_beta(sym_mat3<FloatType> const& u_star)
00169 {
00170 return FloatType(scitbx::constants::two_pi_sq) * u_star;
00171 }
00172
00174
00176 template <typename FloatType>
00177 inline sym_mat3<FloatType>
00178 beta_as_u_star(sym_mat3<FloatType> const& beta)
00179 {
00180 return beta / FloatType(scitbx::constants::two_pi_sq);
00181 }
00182
00184
00187 template <typename FloatType>
00188 inline sym_mat3<FloatType>
00189 u_cart_as_beta(
00190 uctbx::unit_cell const& unit_cell,
00191 sym_mat3<FloatType> const& u_cart)
00192 {
00193 return u_star_as_beta(u_cart_as_u_star(unit_cell, u_cart));
00194 }
00195
00197
00200 template <typename FloatType>
00201 inline sym_mat3<FloatType>
00202 beta_as_u_cart(
00203 uctbx::unit_cell const& unit_cell,
00204 sym_mat3<FloatType> const& beta)
00205 {
00206 return u_star_as_u_cart(unit_cell, beta_as_u_star(beta));
00207 }
00208
00210
00213 template <typename FloatType>
00214 inline sym_mat3<FloatType>
00215 u_cif_as_beta(
00216 uctbx::unit_cell const& unit_cell,
00217 sym_mat3<FloatType> const& u_cif)
00218 {
00219 return u_star_as_beta(u_cif_as_u_star(unit_cell, u_cif));
00220 }
00221
00223
00226 template <typename FloatType>
00227 inline sym_mat3<FloatType>
00228 beta_as_u_cif(
00229 uctbx::unit_cell const& unit_cell,
00230 sym_mat3<FloatType> const& beta)
00231 {
00232 return u_star_as_u_cif(unit_cell, beta_as_u_star(beta));
00233 }
00234
00236
00239 template <typename FloatType>
00240 inline FloatType
00241 u_cart_as_u_iso(sym_mat3<FloatType> const& u_cart)
00242 {
00243 return (u_cart[0] + u_cart[1] + u_cart[2]) / 3.;
00244 }
00245
00247
00250 template <typename FloatType>
00251 sym_mat3<FloatType>
00252 u_iso_as_u_cart(FloatType const& u_iso)
00253 {
00254 return sym_mat3<FloatType>(u_iso,u_iso,u_iso,0,0,0);
00255 }
00256
00258
00261 template <typename FloatType>
00262 inline FloatType
00263 u_star_as_u_iso(
00264 uctbx::unit_cell const& unit_cell,
00265 sym_mat3<FloatType> const& u_star)
00266 {
00267 return u_cart_as_u_iso(u_star_as_u_cart(unit_cell, u_star));
00268 }
00269
00271
00273 template <typename FloatType>
00274 inline sym_mat3<FloatType>
00275 u_iso_as_u_star(
00276 uctbx::unit_cell const& unit_cell,
00277 FloatType const& u_iso)
00278 {
00279 return u_cart_as_u_star(unit_cell, u_iso_as_u_cart(u_iso));
00280 }
00281
00283
00286 template <typename FloatType>
00287 inline FloatType
00288 u_cif_as_u_iso(
00289 uctbx::unit_cell const& unit_cell,
00290 sym_mat3<FloatType> const& u_cif)
00291 {
00292 return u_cart_as_u_iso(u_cif_as_u_cart(unit_cell, u_cif));
00293 }
00294
00296
00299 template <typename FloatType>
00300 inline sym_mat3<FloatType>
00301 u_iso_as_u_cif(
00302 uctbx::unit_cell const& unit_cell,
00303 FloatType const& u_iso)
00304 {
00305 return u_cart_as_u_cif(unit_cell, u_iso_as_u_cart(u_iso));
00306 }
00307
00309
00312 template <typename FloatType>
00313 inline FloatType
00314 beta_as_u_iso(
00315 uctbx::unit_cell const& unit_cell,
00316 sym_mat3<FloatType> const& beta)
00317 {
00318 return u_cart_as_u_iso(beta_as_u_cart(unit_cell, beta));
00319 }
00320
00322
00325 template <typename FloatType>
00326 inline sym_mat3<FloatType>
00327 u_iso_as_beta(
00328 uctbx::unit_cell const& unit_cell,
00329 FloatType const& u_iso)
00330 {
00331 return u_cart_as_beta(unit_cell, u_iso_as_u_cart(u_iso));
00332 }
00333
00335 template <typename FloatType=double>
00336 struct factor_u_cart_u_iso
00337 {
00338 factor_u_cart_u_iso() {}
00339
00340 factor_u_cart_u_iso(sym_mat3<FloatType> const& u_cart)
00341 :
00342 u_iso(u_cart_as_u_iso(u_cart)),
00343 u_cart_minus_u_iso(u_cart)
00344 {
00345 for(unsigned i=0;i<3;i++) {
00346 u_cart_minus_u_iso[i] -= u_iso;
00347 }
00348 }
00349
00350 FloatType u_iso;
00351 sym_mat3<FloatType> u_cart_minus_u_iso;
00352 };
00353
00355 template <typename FloatType=double>
00356 struct factor_u_star_u_iso
00357 {
00358 factor_u_star_u_iso() {}
00359
00360 factor_u_star_u_iso(
00361 uctbx::unit_cell const& unit_cell,
00362 sym_mat3<FloatType> const& u_star)
00363 {
00364 factor_u_cart_u_iso<FloatType>
00365 f_cart(u_star_as_u_cart(unit_cell, u_star));
00366 u_iso = f_cart.u_iso;
00367 u_star_minus_u_iso = u_cart_as_u_star(
00368 unit_cell, f_cart.u_cart_minus_u_iso);
00369 }
00370
00371 FloatType u_iso;
00372 sym_mat3<FloatType> u_star_minus_u_iso;
00373 };
00374
00376 template <typename FloatType=double>
00377 struct factor_u_cif_u_iso
00378 {
00379 factor_u_cif_u_iso() {}
00380
00381 factor_u_cif_u_iso(
00382 uctbx::unit_cell const& unit_cell,
00383 sym_mat3<FloatType> const& u_cif)
00384 {
00385 factor_u_cart_u_iso<FloatType>
00386 f_cart(u_cif_as_u_cart(unit_cell, u_cif));
00387 u_iso = f_cart.u_iso;
00388 u_cif_minus_u_iso = u_cart_as_u_cif(
00389 unit_cell, f_cart.u_cart_minus_u_iso);
00390 }
00391
00392 FloatType u_iso;
00393 sym_mat3<FloatType> u_cif_minus_u_iso;
00394 };
00395
00397 template <typename FloatType=double>
00398 struct factor_beta_u_iso
00399 {
00400 factor_beta_u_iso() {}
00401
00402 factor_beta_u_iso(
00403 uctbx::unit_cell const& unit_cell,
00404 sym_mat3<FloatType> const& beta)
00405 {
00406 factor_u_cart_u_iso<FloatType>
00407 f_cart(beta_as_u_cart(unit_cell, beta));
00408 u_iso = f_cart.u_iso;
00409 beta_minus_u_iso = u_cart_as_beta(
00410 unit_cell, f_cart.u_cart_minus_u_iso);
00411 }
00412
00413 FloatType u_iso;
00414 sym_mat3<FloatType> beta_minus_u_iso;
00415 };
00416
00418 inline double
00419 debye_waller_factor_exp(const char* u_type, double arg, double max_arg=50)
00420 {
00421 if (arg > max_arg) {
00422 char buf[256];
00423 std::sprintf(buf,
00424 "cctbx::adptbx::debye_waller_factor_exp:"
00425 " max_arg exceeded (%s):"
00426 " arg = %.6g"
00427 " max_arg = %.6g",
00428 u_type, arg, max_arg);
00429 throw std::runtime_error(buf);
00430 }
00431 return std::exp(arg);
00432 }
00433
00435 inline double
00436 debye_waller_factor_b_iso(
00437 double stol_sq,
00438 double b_iso)
00439 {
00440 return debye_waller_factor_exp("isotropic", -b_iso * stol_sq);
00441 }
00442
00444 inline double
00445 debye_waller_factor_u_iso(
00446 double stol_sq,
00447 double u_iso)
00448 {
00449 return debye_waller_factor_b_iso(stol_sq, u_as_b(u_iso));
00450 }
00451
00453 inline double
00454 debye_waller_factor_b_iso(
00455 uctbx::unit_cell const& unit_cell,
00456 miller::index<> const& h,
00457 double b_iso)
00458 {
00459 return debye_waller_factor_b_iso(unit_cell.d_star_sq(h) / 4., b_iso);
00460 }
00461
00463 inline double
00464 debye_waller_factor_u_iso(
00465 uctbx::unit_cell const& unit_cell,
00466 miller::index<> const& h,
00467 double u_iso)
00468 {
00469 return debye_waller_factor_b_iso(unit_cell, h, u_as_b(u_iso));
00470 }
00471
00473 template <typename FloatType>
00474 inline FloatType
00475 debye_waller_factor_u_star(
00476 miller::index<> const& h,
00477 sym_mat3<FloatType> const& u_star)
00478 {
00479 return debye_waller_factor_exp(
00480 "anisotropic",
00481 -scitbx::constants::two_pi_sq * (
00482 (h[0] * h[0]) * u_star[0]
00483 + (h[1] * h[1]) * u_star[1]
00484 + (h[2] * h[2]) * u_star[2]
00485 + (2 * h[0] * h[1]) * u_star[3]
00486 + (2 * h[0] * h[2]) * u_star[4]
00487 + (2 * h[1] * h[2]) * u_star[5]));
00488 }
00489
00491
00508 template <typename NumType>
00509 inline sym_mat3<NumType>
00510 debye_waller_factor_u_star_gradient_coefficients(
00511 miller::index<> const& h,
00512 scitbx::type_holder<NumType> holder=scitbx::type_holder<NumType>())
00513 {
00514 return sym_mat3<NumType>(
00515 (h[0] * h[0]),
00516 (h[1] * h[1]),
00517 (h[2] * h[2]),
00518 (2 * h[0] * h[1]),
00519 (2 * h[0] * h[2]),
00520 (2 * h[1] * h[2]));
00521 }
00522
00524
00530 template <typename NumType>
00531 af::shared<NumType>
00532 debye_waller_factor_u_star_curvature_coefficients(
00533 miller::index<> const& h,
00534 scitbx::type_holder<NumType>)
00535 {
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 af::shared<NumType> result(6*(6+1)/2, 0);
00548 NumType h0h0 = h[0] * h[0];
00549 NumType h1h1 = h[1] * h[1];
00550 NumType h2h2 = h[2] * h[2];
00551 NumType h0h1 = 2 * h[0] * h[1];
00552 NumType h0h2 = 2 * h[0] * h[2];
00553 NumType h1h2 = 2 * h[1] * h[2];
00554 result[0] = h0h0 * h0h0;
00555 result[1] = h0h0 * h1h1;
00556 result[2] = h0h0 * h2h2;
00557 result[3] = h0h0 * h0h1;
00558 result[4] = h0h0 * h0h2;
00559 result[5] = h0h0 * h1h2;
00560 result[6] = h1h1 * h1h1;
00561 result[7] = h1h1 * h2h2;
00562 result[8] = h1h1 * h0h1;
00563 result[9] = h1h1 * h0h2;
00564 result[10] = h1h1 * h1h2;
00565 result[11] = h2h2 * h2h2;
00566 result[12] = h2h2 * h0h1;
00567 result[13] = h2h2 * h0h2;
00568 result[14] = h2h2 * h1h2;
00569 result[15] = h0h1 * h0h1;
00570 result[16] = h0h1 * h0h2;
00571 result[17] = h0h1 * h1h2;
00572 result[18] = h0h2 * h0h2;
00573 result[19] = h0h2 * h1h2;
00574 result[20] = h1h2 * h1h2;
00575 return result;
00576 }
00577
00579 template <typename FloatType>
00580 inline FloatType
00581 debye_waller_factor_beta(
00582 miller::index<> const& h,
00583 sym_mat3<FloatType> const& beta)
00584 {
00585 return debye_waller_factor_u_star(h, beta_as_u_star(beta));
00586 }
00587
00589 template <typename FloatType>
00590 inline FloatType
00591 debye_waller_factor_u_cif(
00592 uctbx::unit_cell const& unit_cell,
00593 miller::index<> const& h,
00594 sym_mat3<FloatType> const& u_cif)
00595 {
00596 return debye_waller_factor_u_star(h, u_cif_as_u_star(unit_cell, u_cif));
00597 }
00598
00600 template <typename FloatType>
00601 inline FloatType
00602 debye_waller_factor_u_cart(
00603 uctbx::unit_cell const& unit_cell,
00604 miller::index<> const& h,
00605 sym_mat3<FloatType> const& u_cart)
00606 {
00607 return debye_waller_factor_u_star(h, u_cart_as_u_star(unit_cell, u_cart));
00608 }
00609
00611
00613 template <typename FloatType>
00614 inline sym_mat3<FloatType>
00615 grad_u_star_as_u_cart(
00616 uctbx::unit_cell const& unit_cell,
00617 sym_mat3<FloatType> const& grad_u_star)
00618 {
00619 return scitbx::matrix::tensor_rank_2::gradient_transform(
00620 unit_cell.fractionalization_matrix(), grad_u_star);
00621 }
00622
00624
00626 template <typename FloatType>
00627 af::shared<sym_mat3<FloatType> >
00628 grad_u_star_as_u_cart(
00629 uctbx::unit_cell const& unit_cell,
00630 af::const_ref<sym_mat3<FloatType> > const& grad_u_star)
00631 {
00632 af::shared<sym_mat3<FloatType> > result((af::reserve(grad_u_star.size())));
00633 for(std::size_t i=0;i<grad_u_star.size();i++) {
00634 result.push_back(grad_u_star_as_u_cart(unit_cell, grad_u_star[i]));
00635 }
00636 return result;
00637 }
00638
00640
00642 template <typename FloatType>
00643 inline sym_mat3<FloatType>
00644 grad_u_cart_as_u_star(
00645 uctbx::unit_cell const& unit_cell,
00646 sym_mat3<FloatType> const& grad_u_cart)
00647 {
00648 return scitbx::matrix::tensor_rank_2::gradient_transform(
00649 unit_cell.orthogonalization_matrix(), grad_u_cart);
00650 }
00651
00653
00655 template <typename FloatType>
00656 af::shared<sym_mat3<FloatType> >
00657 grad_u_cart_as_u_star(
00658 uctbx::unit_cell const& unit_cell,
00659 af::const_ref<sym_mat3<FloatType> > const& grad_u_cart)
00660 {
00661 af::shared<sym_mat3<FloatType> > result((af::reserve(grad_u_cart.size())));
00662 for(std::size_t i=0;i<grad_u_cart.size();i++) {
00663 result.push_back(grad_u_cart_as_u_star(unit_cell, grad_u_cart[i]));
00664 }
00665 return result;
00666 }
00667
00669 template <typename FloatType>
00670 class eigensystem
00671 {
00672 public:
00674 eigensystem() {}
00675
00686 eigensystem(sym_mat3<FloatType> const& adp)
00687 {
00688 scitbx::matrix::eigensystem::real_symmetric<FloatType> es(adp);
00689 for(std::size_t i=0;i<3;i++) {
00690 vectors_[i] = vec3<FloatType>(&es.vectors()[i*3]);
00691 }
00692 values_ = vec3<FloatType>(es.values().begin());
00693 }
00694
00696
00698 vec3<FloatType> const&
00699 vectors(std::size_t i) const
00700 {
00701 if (i >= vectors_.size()) throw error_index();
00702 return vectors_[i];
00703 }
00704
00706 vec3<FloatType> const&
00707 values() const { return values_; }
00708
00709 private:
00710 af::tiny<vec3<FloatType>, 3> vectors_;
00711 vec3<FloatType> values_;
00712 };
00713
00715
00717 template <typename FloatType>
00718 vec3<FloatType>
00719 eigenvalues(sym_mat3<FloatType> const& adp)
00720 {
00721 return eigensystem<FloatType>(adp).values();
00722 }
00723
00731 template <typename FloatType>
00732 bool
00733 is_positive_definite(
00734 vec3<FloatType> const& adp_eigenvalues)
00735 {
00736 return scitbx::af::min(adp_eigenvalues.const_ref()) > 0;
00737 }
00738
00746 template <typename FloatType>
00747 bool
00748 is_positive_definite(
00749 vec3<FloatType> const& adp_eigenvalues,
00750 FloatType const& tolerance)
00751 {
00752 return scitbx::af::min(adp_eigenvalues.const_ref()) >= -tolerance;
00753 }
00754
00760 template <typename FloatType>
00761 bool
00762 is_positive_definite(
00763 sym_mat3<FloatType> const& adp)
00764 {
00765 return is_positive_definite(eigenvalues(adp));
00766 }
00767
00773 template <typename FloatType>
00774 bool
00775 is_positive_definite(
00776 sym_mat3<FloatType> const& adp,
00777 FloatType const& tolerance)
00778 {
00779 return is_positive_definite(eigenvalues(adp), tolerance);
00780 }
00781
00786 template <typename FloatType>
00787 af::shared<bool>
00788 is_positive_definite(
00789 af::const_ref<sym_mat3<FloatType> > const& adp,
00790 FloatType const& tolerance)
00791 {
00792 af::shared<bool> result((af::reserve(adp.size())));
00793 for(std::size_t i=0;i<adp.size();i++) {
00794 result.push_back(is_positive_definite(eigenvalues(adp[i]), tolerance));
00795 }
00796 return result;
00797 }
00798
00800
00802 template <typename FloatType>
00803 sym_mat3<FloatType>
00804 eigenvalue_filtering(
00805 sym_mat3<FloatType> const& u_cart,
00806 FloatType const& u_min=0,
00807 FloatType const& u_max=0)
00808 {
00809 scitbx::matrix::eigensystem::real_symmetric<FloatType> es(u_cart);
00810 scitbx::vec3<FloatType> es_val(es.values().begin());
00811 for(std::size_t i=0;i<3;i++) if (es_val[i] < u_min) es_val[i] = u_min;
00812 if (u_max > 0) {
00813 for(std::size_t i=0;i<3;i++) if (es_val[i] > u_max) es_val[i] = u_max;
00814 }
00815 scitbx::mat3<FloatType> es_vec(es.vectors().begin());
00816 scitbx::mat3<FloatType> es_vec_inv = es_vec.inverse();
00817 return sym_mat3<FloatType>(es_val).tensor_transform(es_vec_inv);
00818 }
00819
00821
00823 template <typename FloatType>
00824 inline sym_mat3<FloatType>
00825 c_u_c_transpose(mat3<FloatType> const& c, sym_mat3<FloatType> const& u)
00826 {
00827 return u.tensor_transform(c);
00828 }
00829
00830 }}
00831
00832 #endif // CCTBX_ADPTBX_H
