cctbx/maptbx/grid_indices_around_sites.h

#ifndef CCTBX_MAPTBX_GRID_INDICES_AROUND_SITES_H
#define CCTBX_MAPTBX_GRID_INDICES_AROUND_SITES_H

#include <cctbx/uctbx.h>
#include <scitbx/math/modulo.h>
#include <boost/scoped_array.hpp>

namespace cctbx { namespace maptbx {

  template <typename SetType>
  std::auto_ptr<SetType>
  grid_indices_around_sites(
    uctbx::unit_cell const& unit_cell,
    af::tiny<int, 3> const& fft_n_real,
    af::tiny<int, 3> const& fft_m_real,
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<double> const& site_radii)
  {
    CCTBX_ASSERT(unit_cell.volume() > 0);
    CCTBX_ASSERT(fft_n_real.const_ref().all_gt(0));
    CCTBX_ASSERT(fft_m_real.const_ref().all_ge(fft_n_real.const_ref()));
    CCTBX_ASSERT(site_radii.size() == sites_cart.size());
    typedef typename SetType::value_type svt;
    {
      svt m[3];
      for(unsigned i=0;i<3;i++) {
        m[i] = boost::numeric_cast<svt>(fft_m_real[i]);
      }
      if (scitbx::math::unsigned_product_leads_to_overflow(m, 3U)) {
        throw std::runtime_error(
          "product of fft_m_real grid dimensions leads to SetType::value_type"
          " integer overflow (i.e. the grid is too large).");
      }
    }
    std::auto_ptr<SetType> result_ap(new SetType);
    SetType* result = result_ap.get();
    scitbx::vec3<double> fft_n_real_f;
    scitbx::vec3<double> one_over_grid_spacing;
    for(unsigned i=0;i<3;i++) {
      fft_n_real_f[i] = boost::numeric_cast<double>(fft_n_real[i]);
      double unit_cell_parameter = unit_cell.parameters()[i];
      CCTBX_ASSERT(unit_cell_parameter > 0);
      one_over_grid_spacing[i] = fft_n_real_f[i] / unit_cell_parameter;
    }
    svt n0 = boost::numeric_cast<svt>(fft_n_real[0]);
    svt n1 = boost::numeric_cast<svt>(fft_n_real[1]);
    svt n2 = boost::numeric_cast<svt>(fft_n_real[2]);
    svt m1 = boost::numeric_cast<svt>(fft_m_real[1]);
    svt m2 = boost::numeric_cast<svt>(fft_m_real[2]);
    for(std::size_t i_site=0;i_site<sites_cart.size();i_site++) {
      scitbx::vec3<double> uvw = unit_cell.fractionalize(sites_cart[i_site]);
      for(unsigned i=0;i<3;i++) uvw[i] *= fft_n_real_f[i];
      double site_radius = site_radii[i_site];
      af::tiny<svt, 3> box_min;
      af::tiny<svt, 3> box_max;
      for(unsigned i=0;i<3;i++) {
        double grid_radius = std::min(
          site_radius * one_over_grid_spacing[i],
          fft_n_real_f[i] * (0.5+1.e-5));
        typedef scitbx::math::float_int_conversions<double, int> fic;
        int box_min_i = fic::ifloor(uvw[i] - grid_radius);
        int box_max_i = fic::iceil( uvw[i] + grid_radius);
        int box_min_mod = scitbx::math::mod_positive(box_min_i, fft_n_real[i]);
        int shift = box_min_mod - box_min_i;
        uvw[i] += boost::numeric_cast<double>(shift);
        box_min[i] = boost::numeric_cast<svt>(box_min_i + shift);
        box_max[i] = boost::numeric_cast<svt>(box_max_i + shift);
      }
      double site_radius_sq = site_radius * site_radius;
      scitbx::mat3<double> orth = unit_cell.orthogonalization_matrix();
      CCTBX_ASSERT(orth[3] == 0);
      CCTBX_ASSERT(orth[6] == 0);
      CCTBX_ASSERT(orth[7] == 0);
      orth[0] /= fft_n_real_f[0];
      orth[1] /= fft_n_real_f[1];
      orth[4] /= fft_n_real_f[1];
      orth[2] /= fft_n_real_f[2];
      orth[5] /= fft_n_real_f[2];
      orth[8] /= fft_n_real_f[2];
      svt box_n1_2 = (box_max[1] + 1 - box_min[1]) * 2;
      svt box_n2_3 = (box_max[2] + 1 - box_min[2]) * 3;
      boost::scoped_array<double> buffer(new double[box_n1_2 + box_n2_3]);
      double* o14d1_beg = buffer.get();
      double* o258d2_beg = o14d1_beg + box_n1_2;
      const double* o258d2_end = o258d2_beg + box_n2_3;
      {
        double* o14d1 = o14d1_beg;
        for(svt b1=box_min[1];b1<=box_max[1];b1++) {
          double diff1 = boost::numeric_cast<double>(b1) - uvw[1];
          *o14d1++ = orth[1] * diff1;
          *o14d1++ = orth[4] * diff1;
        }
      }
      {
        double* o258d2 = o258d2_beg;
        for(svt b2=box_min[2];b2<=box_max[2];b2++) {
          double diff2 = boost::numeric_cast<double>(b2) - uvw[2];
          *o258d2++ =                  orth[2] * diff2;
          *o258d2++ =                  orth[5] * diff2;
          *o258d2++ = scitbx::fn::pow2(orth[8] * diff2);
        }
      }
      svt g0 = box_min[0];
      svt g1b = box_min[1];
      svt g2b = box_min[2];
      for(svt b0=box_min[0];b0<=box_max[0];b0++,g0++) {
        double diff0 = boost::numeric_cast<double>(b0) - uvw[0];
        double orth0_diff0 = orth[0] * diff0;
        if (g0 == n0) g0 = 0;
        svt gm0 = g0 * m1;
        svt g1 = g1b;
        for(const double* o14d1=o14d1_beg;o14d1!=o258d2_beg;g1++) {
          double orth0_diff0_orth1_diff_1 = orth0_diff0 + *o14d1++;
          double orth4_diff1 = *o14d1++;
          if (g1 == n1) g1 = 0;
          svt gm01 = (gm0 + g1) * m2;
          svt g2 = g2b;
          for(const double* o258d2=o258d2_beg;o258d2!=o258d2_end;g2++) {
            if (g2 == n2) g2 = 0;
            double dc_sq = orth0_diff0_orth1_diff_1 + *o258d2++;
            dc_sq *= dc_sq;
            double dc1 = orth4_diff1 + *o258d2++;
            dc_sq += dc1 * dc1;
            dc_sq += *o258d2++;
            if (dc_sq > site_radius_sq) continue;
            result->insert(gm01 + g2);
          }
        }
      }
    }
    return result_ap;
  }

}} // namespace cctbx::maptbx

#endif // GUARD

---