00001 #ifndef CCTBX_CRYSTAL_INCREMENTAL_PAIRS_H
00002 #define CCTBX_CRYSTAL_INCREMENTAL_PAIRS_H
00003
00004 #include <cctbx/crystal/pair_tables.h>
00005
00006 namespace cctbx { namespace crystal {
00007
00008 template <typename FloatType=double, typename IntShiftType=int>
00009 class incremental_pairs
00010 {
00011 protected:
00012 FloatType distance_cutoff_;
00013 FloatType distance_cutoff_sq_;
00014 FloatType cubicle_epsilon_;
00015 typedef
00016 direct_space_asu::asu_mappings<FloatType, IntShiftType>
00017 asu_mappings_t;
00018 boost::shared_ptr<asu_mappings_t> asu_mappings_owner_;
00019 asu_mappings_t* asu_mappings_;
00020 typedef std::vector<direct_space_asu::asu_mapping_index>
00021 cubicle_content_t;
00022 scitbx::cubicles<cubicle_content_t, FloatType> cubicles_;
00023 typedef crystal::pair_asu_table<FloatType, IntShiftType>
00024 pair_asu_table_t;
00025 boost::shared_ptr<pair_asu_table_t> pair_asu_table_owner_;
00026 pair_asu_table_t* pair_asu_table_;
00027
00028 public:
00030 incremental_pairs() {}
00031
00032 incremental_pairs(
00033 sgtbx::space_group const& space_group,
00034 direct_space_asu::float_asu<FloatType> const& asu,
00035 FloatType const& distance_cutoff,
00036 FloatType const& asu_mappings_buffer_thickness=-1,
00037 FloatType const& cubicle_epsilon=-1)
00038 :
00039 distance_cutoff_(distance_cutoff),
00040 distance_cutoff_sq_(distance_cutoff*distance_cutoff),
00041 cubicle_epsilon_(cubicle_epsilon >= 0
00042 ? cubicle_epsilon
00043 : asu.is_inside_epsilon()),
00044 asu_mappings_owner_(new asu_mappings_t(
00045 space_group,
00046 asu,
00047 asu_mappings_buffer_thickness >= 0
00048 ? asu_mappings_buffer_thickness
00049 : distance_cutoff)),
00050 asu_mappings_(asu_mappings_owner_.get()),
00051 cubicles_(
00052 asu_mappings_->asu_buffer().box_min( true),
00053 asu_mappings_->asu_buffer().box_span( true),
00054 distance_cutoff_,
00055 cubicle_epsilon_),
00056 pair_asu_table_owner_(new pair_asu_table_t(asu_mappings_owner_)),
00057 pair_asu_table_(pair_asu_table_owner_.get()),
00058 min_distance_sym_equiv(0.5),
00059 assert_min_distance_sym_equiv(true)
00060 {
00061 CCTBX_ASSERT(distance_cutoff_ > 0);
00062 CCTBX_ASSERT(asu_mappings_->buffer_thickness() >= distance_cutoff_);
00063 }
00064
00065 boost::shared_ptr<
00066 direct_space_asu::asu_mappings<FloatType, IntShiftType> >
00067 asu_mappings() const { return asu_mappings_owner_; }
00068
00069 boost::shared_ptr<
00070 crystal::pair_asu_table<FloatType, IntShiftType> >
00071 pair_asu_table() const { return pair_asu_table_owner_; }
00072
00073 FloatType min_distance_sym_equiv;
00074 bool assert_min_distance_sym_equiv;
00075
00076 void
00077 process_site_frac(
00078 fractional<FloatType> const& original_site,
00079 sgtbx::site_symmetry_ops const& site_symmetry_ops)
00080 {
00081 direct_space_asu::asu_mapping_index mi;
00082 mi.i_seq = asu_mappings_->mappings().size();
00083 asu_mappings_->process(original_site, site_symmetry_ops);
00084 pair_asu_table_->grow(1);
00085 sgtbx::rt_mx
00086 rt_mx_i_inverse = asu_mappings_->get_rt_mx(mi.i_seq, 0).inverse();
00087 af::const_ref<
00088 typename asu_mappings_t::array_of_mappings_for_one_site> const&
00089 mappings = asu_mappings_->mappings_const_ref();
00090 typename asu_mappings_t::array_of_mappings_for_one_site const&
00091 mappings_i = mappings.back();
00092 scitbx::vec3<unsigned> n_cub = cubicles_.ref.accessor();
00093 for(mi.i_sym=0; mi.i_sym<mappings_i.size(); mi.i_sym++) {
00094 scitbx::vec3<unsigned>
00095 i_cub = cubicles_.i_cubicle(mappings_i[mi.i_sym].mapped_site());
00096 scitbx::vec3<unsigned> j_cub_min;
00097 scitbx::vec3<unsigned> j_cub_max;
00098 j_cub_min[0] = (i_cub[0] == 0 ? 0 : i_cub[0]-1);
00099 j_cub_max[0] = (i_cub[0] == n_cub[0]-1 ? i_cub[0] : i_cub[0]+1);
00100 j_cub_min[1] = (i_cub[1] == 0 ? 0 : i_cub[1]-1);
00101 j_cub_max[1] = (i_cub[1] == n_cub[1]-1 ? i_cub[1] : i_cub[1]+1);
00102 j_cub_min[2] = (i_cub[2] == 0 ? 0 : i_cub[2]-1);
00103 j_cub_max[2] = (i_cub[2] == n_cub[2]-1 ? i_cub[2] : i_cub[2]+1);
00104 scitbx::vec3<unsigned> j_cub;
00105 for(j_cub[0]=j_cub_min[0];j_cub[0]<=j_cub_max[0];j_cub[0]++) {
00106 for(j_cub[1]=j_cub_min[1];j_cub[1]<=j_cub_max[1];j_cub[1]++) {
00107 for(j_cub[2]=j_cub_min[2];j_cub[2]<=j_cub_max[2];j_cub[2]++) {
00108 const cubicle_content_t* cub_j = &cubicles_.ref(j_cub);
00109 typename cubicle_content_t::const_iterator cub_ji;
00110 if (mi.i_sym == 0) {
00111 for(cub_ji=cub_j->begin();
00112 cub_ji!=cub_j->end();
00113 cub_ji++) {
00114 scitbx::vec3<FloatType> diff_vec =
00115 mappings[cub_ji->i_seq][cub_ji->i_sym].mapped_site()
00116 - mappings_i[0].mapped_site();
00117 if (diff_vec.length_sq() > distance_cutoff_sq_) continue;
00118 pair_asu_table_->add_pair(
00119 mi.i_seq,
00120 cub_ji->i_seq,
00121 rt_mx_i_inverse.multiply(
00122 asu_mappings_->get_rt_mx(cub_ji->i_seq, cub_ji->i_sym)));
00123 }
00124 }
00125 else {
00126 for(cub_ji=cub_j->begin();
00127 cub_ji!=cub_j->end();
00128 cub_ji++) {
00129 if (cub_ji->i_sym != 0) continue;
00130 scitbx::vec3<FloatType> diff_vec =
00131 mappings_i[mi.i_sym].mapped_site()
00132 - mappings[cub_ji->i_seq][0].mapped_site();
00133 if (diff_vec.length_sq() > distance_cutoff_sq_) continue;
00134 pair_asu_table_->add_pair(
00135 cub_ji->i_seq,
00136 mi.i_seq,
00137 asu_mappings_->get_rt_mx(cub_ji->i_seq, 0)
00138 .inverse().multiply(
00139 asu_mappings_->get_rt_mx(mi.i_seq, mi.i_sym)));
00140 }
00141 }
00142 }}}
00143 std::size_t i1d_cub = cubicles_.ref.accessor()(i_cub);
00144 cubicles_.ref[i1d_cub].push_back(mi);
00145 }
00146 }
00147
00148 void
00149 process_site_frac(
00150 fractional<FloatType> const& original_site)
00151 {
00152 process_site_frac(
00153 original_site,
00154 sgtbx::site_symmetry(
00155 asu_mappings_->asu().unit_cell(),
00156 asu_mappings_->space_group(),
00157 original_site,
00158 min_distance_sym_equiv,
00159 assert_min_distance_sym_equiv));
00160 }
00161
00163 void
00164 reserve_additional(std::size_t n)
00165 {
00166 n += asu_mappings_->mappings_const_ref().size();
00167 asu_mappings_->reserve(n);
00168 pair_asu_table_->reserve(n);
00169 }
00170
00171 void
00172 process_sites_frac(
00173 af::const_ref<scitbx::vec3<FloatType> > const& original_sites,
00174 sgtbx::site_symmetry_table const& site_symmetry_table)
00175 {
00176 CCTBX_ASSERT(site_symmetry_table.indices_const_ref().size()
00177 == original_sites.size());
00178 reserve_additional(original_sites.size());
00179 for(std::size_t i=0;i<original_sites.size();i++) {
00180 process_site_frac(original_sites[i], site_symmetry_table.get(i));
00181 }
00182 }
00183
00184 void
00185 process_sites_frac(
00186 af::const_ref<scitbx::vec3<FloatType> > const& original_sites)
00187 {
00188 reserve_additional(original_sites.size());
00189 for(std::size_t i=0;i<original_sites.size();i++) {
00190 process_site_frac(original_sites[i]);
00191 }
00192 }
00193
00194 void
00195 process_sites_cart(
00196 af::const_ref<scitbx::vec3<FloatType> > const& original_sites,
00197 sgtbx::site_symmetry_table const& site_symmetry_table)
00198 {
00199 CCTBX_ASSERT(site_symmetry_table.indices_const_ref().size()
00200 == original_sites.size());
00201 reserve_additional(original_sites.size());
00202 uctbx::unit_cell const& uc = asu_mappings_->asu().unit_cell();
00203 for(std::size_t i=0;i<original_sites.size();i++) {
00204 process_site_frac(
00205 uc.fractionalize(original_sites[i]), site_symmetry_table.get(i));
00206 }
00207 }
00208
00209 void
00210 process_sites_cart(
00211 af::const_ref<scitbx::vec3<FloatType> > const& original_sites)
00212 {
00213 reserve_additional(original_sites.size());
00214 uctbx::unit_cell const& uc = asu_mappings_->asu().unit_cell();
00215 for(std::size_t i=0;i<original_sites.size();i++) {
00216 process_site_frac(uc.fractionalize(original_sites[i]));
00217 }
00218 }
00219
00220 std::map<long, long>
00221 cubicle_size_counts() const { return cubicles_.cubicle_size_counts(); }
00222 };
00223
00224 }}
00225
00226 #endif // CCTBX_CRYSTAL_INCREMENTAL_PAIRS_H
