00001 #ifndef CCTBX_XRAY_GRADIENTS_DIRECT_H
00002 #define CCTBX_XRAY_GRADIENTS_DIRECT_H
00003
00004 #include <cctbx/xray/scattering_type_registry.h>
00005 #include <cctbx/xray/hr_ht_cache.h>
00006 #include <cctbx/xray/packing_order.h>
00007 #include <cctbx/sgtbx/site_symmetry_table.h>
00008 #include <cctbx/math/cos_sin_table.h>
00009
00010 namespace cctbx { namespace xray { namespace structure_factors {
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 template <typename CosSinType, typename ScattererType>
00021 struct gradients_direct_one_h_one_scatterer
00022 {
00023 typedef typename ScattererType::float_type float_type;
00024
00025 gradients_direct_one_h_one_scatterer(
00026 CosSinType const& cos_sin,
00027 sgtbx::space_group const& space_group,
00028 hr_ht_cache<float_type> const& hr_ht,
00029 float_type d_star_sq,
00030 float_type f0,
00031 ScattererType const& scatterer,
00032 std::complex<float_type> const& d_target_d_f_calc)
00033 :
00034 const_h_sum(0,0)
00035 {
00036 typedef float_type f_t;
00037 typedef std::complex<f_t> c_t;
00038 if (scatterer.flags.grad_site()) d_target_d_site_frac.fill(0);
00039 if (scatterer.flags.grad_u_aniso() &&
00040 scatterer.flags.use_u_aniso()) d_target_d_u_star.fill(0);
00041 fractional<float_type> dtds_term;
00042 scitbx::sym_mat3<f_t> dw_coeff;
00043 f0_fp_fdp = f0 + c_t(scatterer.fp, scatterer.fdp);
00044 c_t f0_fp_fdp_w = f0_fp_fdp * scatterer.weight();
00045 f_t dw_iso;
00046 if (scatterer.flags.use_u_iso()) {
00047 dw_iso = adptbx::debye_waller_factor_u_iso(
00048 d_star_sq/4, scatterer.u_iso);
00049 }
00050 for(std::size_t i=0;i<hr_ht.groups.size();i++) {
00051 hr_ht_group<f_t> const& g = hr_ht.groups[i];
00052 f_t hrx = g.hr * scatterer.site;
00053 f_t ht = g.ht;
00054 if (scatterer.flags.use_u_aniso()) {
00055
00056 dw_coeff = adptbx::debye_waller_factor_u_star_gradient_coefficients(
00057 g.hr, scitbx::type_holder<f_t>());
00058 }
00059 c_t sum_inv(0,0);
00060 if (scatterer.flags.grad_site()) dtds_term.fill(0);
00061 for(std::size_t i=0;i<space_group.f_inv();i++) {
00062 if (i) {
00063 hrx *= -1;
00064 ht = hr_ht.h_inv_t - ht;
00065 }
00066 c_t term = cos_sin.get(hrx + ht);
00067 if (scatterer.flags.grad_site()) {
00068 c_t f = f0_fp_fdp_w * term;
00069 f_t c = d_target_d_f_calc.imag() * f.real()
00070 - d_target_d_f_calc.real() * f.imag();
00071 if (i) c *= -1;
00072 for(std::size_t i=0;i<3;i++) {
00073 dtds_term[i] += g.hr[i] * c;
00074 }
00075 }
00076 sum_inv += term;
00077 }
00078 f_t dw = (scatterer.flags.use_u_iso() ? dw_iso : 1);
00079 if (scatterer.flags.use_u_aniso()) {
00080 dw *= adptbx::debye_waller_factor_u_star(g.hr, scatterer.u_star);
00081 }
00082 sum_inv *= dw;
00083 if (scatterer.flags.grad_site()) dtds_term *= dw;
00084 if (scatterer.flags.grad_site()) d_target_d_site_frac += dtds_term;
00085 if (scatterer.flags.grad_u_aniso() && scatterer.flags.use_u_aniso()) {
00086 c_t f = f0_fp_fdp_w * sum_inv;
00087 f_t c = d_target_d_f_calc.real() * f.real()
00088 + d_target_d_f_calc.imag() * f.imag();
00089 d_target_d_u_star += dw_coeff * c;
00090 }
00091 const_h_sum += sum_inv;
00092 }
00093 }
00094
00095 std::complex<float_type> f0_fp_fdp;
00096 std::complex<float_type> const_h_sum;
00097 fractional<float_type> d_target_d_site_frac;
00098 scitbx::sym_mat3<float_type> d_target_d_u_star;
00099 };
00100
00101 namespace detail {
00102
00103 template <typename FloatType>
00104 struct gradient_refs
00105 {
00106 gradient_refs(
00107 af::ref<scitbx::vec3<FloatType> > site_,
00108 af::ref<FloatType> u_iso_,
00109 af::ref<scitbx::sym_mat3<FloatType> > u_star_,
00110 af::ref<FloatType> occupancy_,
00111 af::ref<FloatType> fp_,
00112 af::ref<FloatType> fdp_)
00113 :
00114 site(site_),
00115 u_iso(u_iso_),
00116 u_star(u_star_),
00117 occupancy(occupancy_),
00118 fp(fp_),
00119 fdp(fdp_)
00120 {}
00121
00122 af::ref<scitbx::vec3<FloatType> > site;
00123 af::ref<FloatType> u_iso;
00124 af::ref<scitbx::sym_mat3<FloatType> > u_star;
00125 af::ref<FloatType> occupancy;
00126 af::ref<FloatType> fp;
00127 af::ref<FloatType> fdp;
00128 };
00129
00130 template <typename ElementType, typename FactorType>
00131 inline
00132 void
00133 in_place_multiply(af::ref<ElementType> const& array, FactorType factor)
00134 {
00135 ElementType* p = array.begin();
00136 ElementType* e = array.end();
00137 while (p != e) (*p++) *= factor;
00138 }
00139
00140 }
00141
00142 template <typename CosSinType, typename ScattererType>
00143 struct gradients_direct_one_h
00144 {
00145 typedef typename ScattererType::float_type float_type;
00146
00147 gradients_direct_one_h(
00148 CosSinType const& cos_sin,
00149 sgtbx::space_group const& space_group,
00150 af::const_ref<ScattererType> const& scatterers,
00151 af::const_ref<std::size_t> scattering_type_indices,
00152 miller::index<> h,
00153 float_type d_star_sq,
00154 af::const_ref<double> const& form_factors,
00155 std::complex<float_type> const& d_target_d_f_calc,
00156 detail::gradient_refs<float_type> gr_refs)
00157 {
00158 typedef float_type f_t;
00159 typedef std::complex<float_type> c_t;
00160 hr_ht_cache<f_t> hr_ht(space_group, h);
00161 for(std::size_t i_sc=0;i_sc<scatterers.size();i_sc++) {
00162 ScattererType const& scatterer = scatterers[i_sc];
00163 f_t f0 = form_factors[scattering_type_indices[i_sc]];
00164 bool gf_u_iso = scatterer.flags.grad_u_iso() &&
00165 scatterer.flags.use_u_iso();
00166 bool gf_u_aniso = scatterer.flags.grad_u_aniso() &&
00167 scatterer.flags.use_u_aniso();
00168 gradients_direct_one_h_one_scatterer<CosSinType, ScattererType> sf(
00169 cos_sin,
00170 space_group,
00171 hr_ht,
00172 d_star_sq,
00173 f0,
00174 scatterer,
00175 d_target_d_f_calc);
00176 if (scatterer.flags.grad_site()) {
00177 gr_refs.site[i_sc] += sf.d_target_d_site_frac;
00178 }
00179 if (gf_u_aniso) {
00180 gr_refs.u_star[i_sc] += sf.d_target_d_u_star;
00181 }
00182 if (gf_u_iso || scatterer.flags.grad_occupancy()) {
00183 c_t t = sf.const_h_sum * sf.f0_fp_fdp;
00184 f_t d = d_target_d_f_calc.real() * t.real()
00185 + d_target_d_f_calc.imag() * t.imag();
00186 d *= scatterer.weight_without_occupancy();
00187 if (gf_u_iso) {
00188 gr_refs.u_iso[i_sc] += d * scatterer.occupancy * d_star_sq;
00189 }
00190 if (scatterer.flags.grad_occupancy()) {
00191 gr_refs.occupancy[i_sc] += d;
00192 }
00193 }
00194 if (scatterer.flags.grad_fp() || scatterer.flags.grad_fdp()) {
00195 c_t f = sf.const_h_sum * scatterer.weight();
00196 if (scatterer.flags.grad_fp()) {
00197 gr_refs.fp[i_sc] += d_target_d_f_calc.real() * f.real()
00198 + d_target_d_f_calc.imag() * f.imag();
00199 }
00200 if (scatterer.flags.grad_fdp()) {
00201 gr_refs.fdp[i_sc] += d_target_d_f_calc.imag() * f.real()
00202 - d_target_d_f_calc.real() * f.imag();
00203 }
00204 }
00205 }
00206 }
00207 };
00208
00209 template <typename ScattererType=scatterer<> >
00210 class gradients_direct
00211 {
00212 public:
00213 typedef ScattererType scatterer_type;
00214 typedef typename ScattererType::float_type float_type;
00215
00216 gradients_direct() {}
00217
00218 gradients_direct(
00219 uctbx::unit_cell const& unit_cell,
00220 sgtbx::space_group const& space_group,
00221 af::const_ref<miller::index<> > const& miller_indices,
00222 af::const_ref<ScattererType> const& scatterers,
00223 af::const_ref<float_type> const& u_iso_refinable_params,
00224 xray::scattering_type_registry const& scattering_type_registry,
00225 sgtbx::site_symmetry_table const& site_symmetry_table,
00226 af::const_ref<std::complex<float_type> > const& d_target_d_f_calc,
00227 std::size_t n_parameters=0)
00228 {
00229 math::cos_sin_exact<float_type> cos_sin;
00230 compute(cos_sin, unit_cell, space_group, miller_indices,
00231 scatterers, u_iso_refinable_params,
00232 scattering_type_registry, site_symmetry_table,
00233 d_target_d_f_calc, n_parameters);
00234 }
00235
00236 gradients_direct(
00237 math::cos_sin_table<float_type> const& cos_sin,
00238 uctbx::unit_cell const& unit_cell,
00239 sgtbx::space_group const& space_group,
00240 af::const_ref<miller::index<> > const& miller_indices,
00241 af::const_ref<ScattererType> const& scatterers,
00242 af::const_ref<float_type> const& u_iso_refinable_params,
00243 xray::scattering_type_registry const& scattering_type_registry,
00244 sgtbx::site_symmetry_table const& site_symmetry_table,
00245 af::const_ref<std::complex<float_type> > const& d_target_d_f_calc,
00246 std::size_t n_parameters=0)
00247 {
00248 compute(cos_sin, unit_cell, space_group, miller_indices,
00249 scatterers, u_iso_refinable_params,
00250 scattering_type_registry, site_symmetry_table,
00251 d_target_d_f_calc, n_parameters);
00252 }
00253
00254 af::shared<float_type>
00255 packed() const { return packed_; }
00256
00257 af::shared<scitbx::vec3<float_type> >
00258 d_target_d_site_frac() const { return d_target_d_site_frac_; }
00259
00260 af::shared<float_type>
00261 d_target_d_u_iso() const { return d_target_d_u_iso_; }
00262
00263 af::shared<scitbx::sym_mat3<float_type> >
00264 d_target_d_u_star() const { return d_target_d_u_star_; }
00265
00266 af::shared<float_type>
00267 d_target_d_occupancy() const { return d_target_d_occupancy_; }
00268
00269 af::shared<float_type>
00270 d_target_d_fp() const { return d_target_d_fp_; }
00271
00272 af::shared<float_type>
00273 d_target_d_fdp() const { return d_target_d_fdp_; }
00274
00275 protected:
00276 af::shared<float_type> packed_;
00277 af::shared<scitbx::vec3<float_type> > d_target_d_site_frac_;
00278 af::shared<float_type> d_target_d_u_iso_;
00279 af::shared<scitbx::sym_mat3<float_type> > d_target_d_u_star_;
00280 af::shared<float_type> d_target_d_occupancy_;
00281 af::shared<float_type> d_target_d_fp_;
00282 af::shared<float_type> d_target_d_fdp_;
00283
00284
00285 static
00286 void
00287 average_special_position_site_gradients(
00288 sgtbx::site_symmetry_table const& site_symmetry_table,
00289 af::ref<scitbx::vec3<float_type> > gradients)
00290 {
00291 CCTBX_ASSERT(gradients.size()
00292 == site_symmetry_table.indices_const_ref().size());
00293 af::const_ref<std::size_t> sp_indices = site_symmetry_table
00294 .special_position_indices().const_ref();
00295 for(std::size_t i_sp=0;i_sp<sp_indices.size();i_sp++) {
00296 std::size_t i_seq = sp_indices[i_sp];
00297 gradients[i_seq] = gradients[i_seq]
00298 * site_symmetry_table.get(i_seq).special_op().r();
00299 }
00300 }
00301
00302 template <typename CosSinType>
00303 void
00304 compute(
00305 CosSinType const& cos_sin,
00306 uctbx::unit_cell const& unit_cell,
00307 sgtbx::space_group const& space_group,
00308 af::const_ref<miller::index<> > const& miller_indices,
00309 af::const_ref<ScattererType> const& scatterers,
00310 af::const_ref<float_type> const& u_iso_refinable_params,
00311 xray::scattering_type_registry const& scattering_type_registry,
00312 sgtbx::site_symmetry_table const& site_symmetry_table,
00313 af::const_ref<std::complex<float_type> > const& d_target_d_f_calc,
00314 std::size_t n_parameters)
00315 {
00316 CCTBX_ASSERT(d_target_d_f_calc.size() == miller_indices.size());
00317 cctbx::xray::scatterer_grad_flags_counts grad_flags_counts(scatterers);
00318 if(grad_flags_counts.tan_u_iso != 0 && grad_flags_counts.u_iso != 0) {
00319 CCTBX_ASSERT(u_iso_refinable_params.size() == scatterers.size());
00320 }
00321 CCTBX_ASSERT(grad_flags_counts.n_parameters() != 0);
00322 typedef float_type f_t;
00323 if (n_parameters != 0) {
00324 packed_.reserve(n_parameters);
00325 }
00326 if (grad_flags_counts.site != 0) {
00327 d_target_d_site_frac_.resize(
00328 scatterers.size(), scitbx::vec3<f_t>(0,0,0));
00329 }
00330 if (grad_flags_counts.u_iso != 0) {
00331 d_target_d_u_iso_.resize(scatterers.size(), 0);
00332 }
00333 if (grad_flags_counts.u_aniso != 0) {
00334 d_target_d_u_star_.resize(
00335 scatterers.size(), scitbx::sym_mat3<f_t>(0,0,0,0,0,0));
00336 }
00337 if (grad_flags_counts.occupancy != 0) {
00338 d_target_d_occupancy_.resize(scatterers.size(), 0);
00339 }
00340 if (grad_flags_counts.fp != 0) {
00341 d_target_d_fp_.resize(scatterers.size(), 0);
00342 }
00343 if (grad_flags_counts.fdp != 0) {
00344 d_target_d_fdp_.resize(scatterers.size(), 0);
00345 }
00346 detail::gradient_refs<f_t> gr_refs(
00347 d_target_d_site_frac_.ref(),
00348 d_target_d_u_iso_.ref(),
00349 d_target_d_u_star_.ref(),
00350 d_target_d_occupancy_.ref(),
00351 d_target_d_fp_.ref(),
00352 d_target_d_fdp_.ref());
00353 af::shared<std::size_t>
00354 scattering_type_indices_memory
00355 = scattering_type_registry.unique_indices(scatterers);
00356 af::const_ref<std::size_t>
00357 scattering_type_indices = scattering_type_indices_memory.const_ref();
00358 for(std::size_t i=0;i<miller_indices.size();i++) {
00359 miller::index<> h = miller_indices[i];
00360 f_t d_star_sq = unit_cell.d_star_sq(h);
00361 af::shared<double>
00362 form_factors_memory
00363 = scattering_type_registry.unique_form_factors_at_d_star_sq(
00364 d_star_sq);
00365 af::const_ref<double> form_factors = form_factors_memory.const_ref();
00366 gradients_direct_one_h<CosSinType, ScattererType>(
00367 cos_sin, space_group, scatterers, scattering_type_indices,
00368 h, d_star_sq, form_factors,
00369 d_target_d_f_calc[i], gr_refs);
00370 }
00371 f_t n_ltr = static_cast<f_t>(space_group.n_ltr());
00372 if (grad_flags_counts.site != 0) {
00373 detail::in_place_multiply(gr_refs.site,
00374 n_ltr * static_cast<f_t>(scitbx::constants::two_pi));
00375 }
00376 if (grad_flags_counts.u_iso != 0) {
00377 detail::in_place_multiply(gr_refs.u_iso,
00378 n_ltr * static_cast<f_t>(-scitbx::constants::two_pi_sq));
00379 }
00380 if (grad_flags_counts.u_aniso != 0) {
00381 detail::in_place_multiply(gr_refs.u_star,
00382 n_ltr * static_cast<f_t>(-scitbx::constants::two_pi_sq));
00383 }
00384 if (grad_flags_counts.occupancy != 0) {
00385 detail::in_place_multiply(gr_refs.occupancy, n_ltr);
00386 }
00387 if (grad_flags_counts.fp != 0) {
00388 detail::in_place_multiply(gr_refs.fp, n_ltr);
00389 }
00390 if (grad_flags_counts.fdp != 0) {
00391 detail::in_place_multiply(gr_refs.fdp, n_ltr);
00392 }
00393 if (grad_flags_counts.u_iso != 0 && grad_flags_counts.tan_u_iso != 0) {
00394 float_type* d_t_d_u = d_target_d_u_iso_.begin();
00395 for(std::size_t i=0;i<scatterers.size();i++,d_t_d_u++) {
00396 ScattererType const& scatterer = scatterers[i];
00397 if (scatterer.flags.tan_u_iso() && scatterer.flags.use_u_iso() &&
00398 scatterer.flags.param > 0 && scatterer.flags.grad_u_iso()) {
00399 f_t pi = scitbx::constants::pi;
00400 f_t u_iso_max = adptbx::b_as_u(
00401 static_cast<f_t>(scatterer.flags.param));
00402 (*d_t_d_u) *= u_iso_max/pi/(1.+u_iso_refinable_params[i]*
00403 u_iso_refinable_params[i]);
00404 }
00405 }
00406 }
00407 if (grad_flags_counts.site != 0) {
00408 average_special_position_site_gradients(
00409 site_symmetry_table,
00410 d_target_d_site_frac_.ref());
00411 }
00412 if (n_parameters != 0) {
00413 BOOST_STATIC_ASSERT(packing_order_convention == 2);
00414 for(std::size_t i=0;i<scatterers.size();i++) {
00415 ScattererType const& scatterer = scatterers[i];
00416 if (scatterer.flags.grad_site()) {
00417 scitbx::vec3<f_t> d_target_d_site_cart =
00418 unit_cell.v_times_fractionalization_matrix_transpose(
00419 gr_refs.site[i]);
00420 for(std::size_t j=0;j<3;j++) {
00421 packed_.push_back(d_target_d_site_cart[j]);
00422 }
00423 }
00424 if(scatterer.flags.use_u_iso()) {
00425
00426 if (scatterer.flags.grad_u_iso()) {
00427 packed_.push_back(gr_refs.u_iso[i]);
00428 }
00429 }
00430 if(scatterer.flags.use_u_aniso()) {
00431 if(scatterer.flags.grad_u_aniso()) {
00432
00433 scitbx::sym_mat3<double> d_target_d_u_cart =
00434 adptbx::grad_u_star_as_u_cart(
00435 unit_cell, gr_refs.u_star[i]);
00436 for(std::size_t j=0;j<6;j++) {
00437 packed_.push_back(d_target_d_u_cart[j]);
00438 }
00439 }
00440 }
00441 if (scatterer.flags.grad_occupancy()) {
00442 packed_.push_back(gr_refs.occupancy[i]);
00443 }
00444
00445
00446 if (grad_flags_counts.fp != 0) {
00447 packed_.push_back(gr_refs.fp[i]);
00448 }
00449 if (grad_flags_counts.fdp != 0) {
00450 packed_.push_back(gr_refs.fdp[i]);
00451 }
00452 }
00453 CCTBX_ASSERT(packed_.size() == n_parameters);
00454 d_target_d_site_frac_ = af::shared<scitbx::vec3<f_t> >();
00455 d_target_d_u_iso_ = af::shared<f_t>();
00456 d_target_d_u_star_ = af::shared<scitbx::sym_mat3<f_t> >();
00457 d_target_d_occupancy_ = af::shared<f_t>();
00458 d_target_d_fp_ = af::shared<f_t>();
00459 d_target_d_fdp_ = af::shared<f_t>();
00460 }
00461 }
00462 };
00463
00464 }}}
00465
00466 #endif // CCTBX_XRAY_GRADIENTS_DIRECT_H