cctbx/uctbx.h

#ifndef CCTBX_UCTBX_H
#define CCTBX_UCTBX_H

#include <cmath>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/optional.hpp>
#include <scitbx/constants.h>
#include <scitbx/sym_mat3.h>
#include <scitbx/array_family/tiny_types.h>
#include <scitbx/array_family/small.h>
#include <scitbx/array_family/shared.h>
#include <scitbx/math/utils.h>
#include <cctbx/coordinates.h>
#include <cctbx/miller.h>

namespace cctbx {

  // forward declaration
  namespace sgtbx
  {
    class rot_mx;
    class change_of_basis_op;
  }

  typedef scitbx::vec3<double> uc_vec3;
  typedef scitbx::mat3<double> uc_mat3;
  typedef scitbx::sym_mat3<double> uc_sym_mat3;

  namespace uctbx {

  inline double d_star_sq_as_stol_sq(double d_star_sq)
  {
    return d_star_sq * .25;
  }

  inline double d_star_sq_as_two_stol(double d_star_sq)
  {
    return std::sqrt(d_star_sq);
  }

  inline double d_star_sq_as_stol(double d_star_sq)
  {
    return std::sqrt(d_star_sq) * .5;
  }

  inline double d_star_sq_as_d(double d_star_sq)
  {
    if (d_star_sq == 0.) return -1.;
    return 1. / std::sqrt(d_star_sq);
  }

  inline double d_star_sq_as_two_theta(double d_star_sq, double wavelength,
                                       bool deg=false)
  {
    double sin_theta = d_star_sq_as_stol(d_star_sq) * wavelength;
    CCTBX_ASSERT(sin_theta <= 1.0);
    double result = 2. * std::asin(sin_theta);
    if (deg) return scitbx::rad_as_deg(result);
    return result;
  }

  template <typename FloatType>
  FloatType
  mean_square_difference(
    scitbx::mat3<FloatType> const& a,
    scitbx::mat3<FloatType> const& b)
  {
    FloatType result = 0;
    for(std::size_t i=0;i<9;i++) {
      result += scitbx::fn::pow2(a[i] - b[i]);
    }
    return result;
  }


  class unit_cell
  {
    public:

      unit_cell() : volume_(0) {}


      explicit
      unit_cell(af::small<double, 6> const& parameters);


      explicit
      unit_cell(af::double6 const& parameters);


      explicit
      unit_cell(uc_sym_mat3 const& metrical_matrix);


      explicit
      unit_cell(uc_mat3 const& orthogonalization_matrix);

      af::double6 const&
      parameters() const { return params_; }

      af::double6 const&
      reciprocal_parameters() const { return r_params_; }

      uc_sym_mat3 const&
      metrical_matrix() const { return metr_mx_; }

      uc_sym_mat3 const&
      reciprocal_metrical_matrix() const { return r_metr_mx_; }

      double
      volume() const { return volume_; }

      unit_cell
      reciprocal() const;

      double
      longest_vector_sq() const;


      double
      shortest_vector_sq() const;

      bool
      is_degenerate(double min_min_length_over_max_length=1.e-10,
                    double min_volume_over_min_length=1.e-5);

      bool
      is_similar_to(unit_cell const& other,
                    double relative_length_tolerance=0.01,
                    double absolute_angle_tolerance=1.) const;

      af::shared<scitbx::mat3<int> >
      similarity_transformations(
        unit_cell const& other,
        double relative_length_tolerance=0.02,
        double absolute_angle_tolerance=2.,
        int unimodular_generator_range=1) const;


      uc_mat3 const& fractionalization_matrix() const { return frac_; }


      uc_mat3 const& orthogonalization_matrix() const { return orth_; }


      template <class IntType>
      uc_mat3
      grid_index_as_site_cart_matrix(
        scitbx::vec3<IntType> const& gridding) const
      {
        uc_mat3 result = orth_;
        for(unsigned i=0;i<3;i++) {
          CCTBX_ASSERT(gridding[i] > 0);
          double f = 1. / boost::numeric_cast<double>(gridding[i]);
          for(unsigned j=0;j<9;j+=3) result[i+j] *= f;
        }
        return result;
      }

      template <class FloatType>
      scitbx::vec3<FloatType>
      fractionalize(scitbx::vec3<FloatType> const& site_cart) const
      {
        // take advantage of the fact that frac_ is upper-triangular.
        return fractional<FloatType>(
            frac_[0] * site_cart[0]
          + frac_[1] * site_cart[1]
          + frac_[2] * site_cart[2],
            frac_[4] * site_cart[1]
          + frac_[5] * site_cart[2],
            frac_[8] * site_cart[2]);
      }

      template <class FloatType>
      scitbx::vec3<FloatType>
      orthogonalize(scitbx::vec3<FloatType> const& site_frac) const
      {
        // take advantage of the fact that orth_ is upper-triangular.
        return cartesian<FloatType>(
            orth_[0] * site_frac[0]
          + orth_[1] * site_frac[1]
          + orth_[2] * site_frac[2],
            orth_[4] * site_frac[1]
          + orth_[5] * site_frac[2],
            orth_[8] * site_frac[2]);
      }

      template <typename FloatType>
      af::shared<scitbx::vec3<FloatType> >
      fractionalize(
        af::const_ref<scitbx::vec3<FloatType> > const& sites_cart) const
      {
        af::shared<scitbx::vec3<FloatType> > result(
          sites_cart.size(),
          af::init_functor_null<scitbx::vec3<FloatType> >());
        const scitbx::vec3<FloatType>* si = sites_cart.begin();
        scitbx::vec3<FloatType>* ri = result.begin();
        for(std::size_t i=0;i<sites_cart.size();i++,ri++,si++) {
          // take advantage of the fact that frac_ is upper-triangular.
          (*ri)[0] = frac_[0] * (*si)[0]
                   + frac_[1] * (*si)[1]
                   + frac_[2] * (*si)[2];
          (*ri)[1] = frac_[4] * (*si)[1]
                   + frac_[5] * (*si)[2];
          (*ri)[2] = frac_[8] * (*si)[2];
        }
        return result;
      }

      template <typename FloatType>
      af::shared<scitbx::vec3<FloatType> >
      orthogonalize(
        af::const_ref<scitbx::vec3<FloatType> > const& sites_frac) const
      {
        af::shared<scitbx::vec3<FloatType> > result(
          sites_frac.size(),
          af::init_functor_null<scitbx::vec3<FloatType> >());
        const scitbx::vec3<FloatType>* si = sites_frac.begin();
        scitbx::vec3<FloatType>* ri = result.begin();
        for(std::size_t i=0;i<sites_frac.size();i++,ri++,si++) {
          // take advantage of the fact that orth_ is upper-triangular.
          (*ri)[0] = orth_[0] * (*si)[0]
                   + orth_[1] * (*si)[1]
                   + orth_[2] * (*si)[2];
          (*ri)[1] = orth_[4] * (*si)[1]
                   + orth_[5] * (*si)[2];
          (*ri)[2] = orth_[8] * (*si)[2];
        }
        return result;
      }


      template <class FloatType>
      scitbx::vec3<FloatType>
      v_times_fractionalization_matrix_transpose(
        scitbx::vec3<FloatType> const& v) const
      {
        // take advantage of the fact that frac_ is upper-triangular.
        return fractional<FloatType>(
            frac_[0] * v[0],
            frac_[1] * v[0]
          + frac_[4] * v[1],
            frac_[2] * v[0]
          + frac_[5] * v[1]
          + frac_[8] * v[2]);
      }

      // ! gradient wrt Cartesian coordinates from gradient wrt fractional ones
      template <class FloatType>
      scitbx::vec3<FloatType>
      orthogonalize_gradient (scitbx::vec3<FloatType> const& v) const
      {
        return v_times_fractionalization_matrix_transpose(v);
      }

      // ! gradient wrt fractional coordinates from gradient wrt Cartesian ones
      template <class FloatType>
      scitbx::vec3<FloatType>
      fractionalize_gradient (const scitbx::vec3<FloatType>& g) const
      {
        const uc_mat3& O = orth_;
        return fractional<FloatType>(O[0]*g[0],
                                     O[1]*g[0] + O[4]*g[1],
                                     O[2]*g[0] + O[5]*g[1] + O[8]*g[2]);
      }


      template <class FloatType>
      FloatType
      length_sq(fractional<FloatType> const& site_frac) const
      {
        return orthogonalize(site_frac).length_sq();
      }

      template <class FloatType>
      FloatType
      length(fractional<FloatType> const& site_frac) const
      {
        return std::sqrt(length_sq(site_frac));
      }


      template <class FloatType>
      FloatType
      distance_sq(fractional<FloatType> const& site_frac_1,
                  fractional<FloatType> const& site_frac_2) const
      {
        return length_sq(fractional<FloatType>(site_frac_1 - site_frac_2));
      }

      template <class FloatType>
      FloatType
      distance(fractional<FloatType> const& site_frac_1,
               fractional<FloatType> const& site_frac_2) const
      {
        return length(fractional<FloatType>(site_frac_1 - site_frac_2));
      }

      template <class FloatType>
      boost::optional<FloatType>
      angle(fractional<FloatType> const& site_frac_1,
            fractional<FloatType> const& site_frac_2,
            fractional<FloatType> const& site_frac_3) const
      {
        cartesian<FloatType> vec_12 = orthogonalize(site_frac_1 - site_frac_2);
        cartesian<FloatType> vec_32 = orthogonalize(site_frac_3 - site_frac_2);
        FloatType length_12 = vec_12.length();
        FloatType length_32 = vec_32.length();
        if (length_12 == 0 || length_32 == 0) {
          return boost::optional<FloatType>();
        }
        const FloatType cos_angle = std::max(-1.,std::min(1.,
          (vec_12 * vec_32)/(length_12 * length_32)));
        return boost::optional<FloatType>(
          std::acos(cos_angle) / scitbx::constants::pi_180);
      }


      template <class FloatType>
      boost::optional<FloatType>
      dihedral(fractional<FloatType> const& site_frac_1,
               fractional<FloatType> const& site_frac_2,
               fractional<FloatType> const& site_frac_3,
               fractional<FloatType> const& site_frac_4) const
      {
        cartesian<FloatType> u = orthogonalize(site_frac_1 - site_frac_2);
        cartesian<FloatType> v = orthogonalize(site_frac_3 - site_frac_2);
        cartesian<FloatType> w = orthogonalize(site_frac_3 - site_frac_4);
        cartesian<FloatType> u_cross_v = u.cross(v);
        cartesian<FloatType> v_cross_w = v.cross(w);
        FloatType norm_u_cross_v = u_cross_v.length_sq();
        if (norm_u_cross_v == 0) return boost::optional<FloatType>();
        FloatType norm_v_cross_w = v_cross_w.length_sq();
        if (norm_v_cross_w == 0) return boost::optional<FloatType>();
        FloatType cos_angle = std::max(-1.,std::min(1.,
          u_cross_v * v_cross_w / std::sqrt(norm_u_cross_v * norm_v_cross_w)));
        FloatType angle = std::acos(cos_angle) / scitbx::constants::pi_180;
        if (u_cross_v * w < 0) {
          angle *= -1;
        }
        return boost::optional<FloatType>(angle);
      }

      template <class FloatType>
      FloatType
      mod_short_length_sq(fractional<FloatType> const& site_frac) const
      {
        return length_sq(site_frac.mod_short());
      }

      template <class FloatType>
      FloatType
      mod_short_length(fractional<FloatType> const& site_frac) const
      {
        return std::sqrt(mod_short_length_sq(site_frac));
      }

      template <class FloatType>
      FloatType
      mod_short_distance_sq(fractional<FloatType> const& site_frac_1,
                            fractional<FloatType> const& site_frac_2) const
      {
        return mod_short_length_sq(
          fractional<FloatType>(site_frac_1 - site_frac_2));
      }

      template <class FloatType>
      FloatType
      mod_short_distance(fractional<FloatType> const& site_frac_1,
                         fractional<FloatType> const& site_frac_2) const
      {
        return std::sqrt(mod_short_distance_sq(site_frac_1, site_frac_2));
      }

      template <class FloatType>
      FloatType
      min_mod_short_distance_sq(
        af::const_ref<scitbx::vec3<FloatType> > const& sites_frac_1,
        fractional<FloatType> const& site_frac_2) const
      {
        FloatType
          result = mod_short_distance_sq(
            fractional<FloatType>(sites_frac_1[0]), site_frac_2);
        for(std::size_t i=1;i<sites_frac_1.size();i++) {
          scitbx::math::update_min(
            result, mod_short_distance_sq(
              fractional<FloatType>(sites_frac_1[i]), site_frac_2));
        }
        return result;
      }

      template <class FloatType>
      FloatType
      min_mod_short_distance(
        af::const_ref<scitbx::vec3<FloatType> > const& sites_frac_1,
        fractional<FloatType> const& site_frac_2) const
      {
        return std::sqrt(min_mod_short_distance_sq(sites_frac_1, site_frac_2));
      }


      unit_cell
      change_basis(uc_mat3 const& c_inv_r, double r_den=1.) const;


      unit_cell
      change_basis(sgtbx::rot_mx const& c_inv_r) const;

      unit_cell
      change_basis(sgtbx::change_of_basis_op const& cb_op) const;

      miller::index<>
      max_miller_indices(double d_min, double tolerance=1.e-4) const;

      template <typename NumType>
      double
      d_star_sq(miller::index<NumType> const& miller_index) const
      {
        return
            (miller_index[0] * miller_index[0]) * r_metr_mx_[0]
          + (miller_index[1] * miller_index[1]) * r_metr_mx_[1]
          + (miller_index[2] * miller_index[2]) * r_metr_mx_[2]
          + (2 * miller_index[0] * miller_index[1]) * r_metr_mx_[3]
          + (2 * miller_index[0] * miller_index[2]) * r_metr_mx_[4]
          + (2 * miller_index[1] * miller_index[2]) * r_metr_mx_[5];
      }

      template <typename NumType>
      af::shared<double>
      d_star_sq(
        af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        af::shared<double> result(
          miller_indices.size(), af::init_functor_null<double>());
        for(std::size_t i=0;i<miller_indices.size();i++) {
          result[i] = d_star_sq(miller_indices[i]);
        }
        return result;
      }

      template <typename NumType>
      double
      max_d_star_sq(
        af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        double result = 0;
        for(std::size_t i=0;i<miller_indices.size();i++) {
          scitbx::math::update_max(result, d_star_sq(miller_indices[i]));
        }
        return result;
      }

      template <typename NumType>
      af::double2
      min_max_d_star_sq(
        af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        af::double2 result(0, 0);
        if (miller_indices.size()) {
          result.fill(d_star_sq(miller_indices[0]));
          for(std::size_t i=1;i<miller_indices.size();i++) {
            double q = d_star_sq(miller_indices[i]);
            scitbx::math::update_min(result[0], q);
            scitbx::math::update_max(result[1], q);
          }
        }
        return result;
      }

      template <typename NumType>
      double
      stol_sq(miller::index<NumType> const& miller_index) const
      {
        return d_star_sq_as_stol_sq(d_star_sq(miller_index));
      }

      template <typename NumType>
      af::shared<double>
      stol_sq(
        af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        af::shared<double> result(
          miller_indices.size(), af::init_functor_null<double>());
        for(std::size_t i=0;i<miller_indices.size();i++) {
          result[i] = stol_sq(miller_indices[i]);
        }
        return result;
      }

      template <typename NumType>
      double
      two_stol(miller::index<NumType> const& miller_index) const
      {
        return d_star_sq_as_two_stol(d_star_sq(miller_index));
      }

      template <typename NumType>
      af::shared<double>
      two_stol(
        af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        af::shared<double> result(
          miller_indices.size(), af::init_functor_null<double>());
        for(std::size_t i=0;i<miller_indices.size();i++) {
          result[i] = two_stol(miller_indices[i]);
        }
        return result;
      }

      template <typename NumType>
      double
      stol(miller::index<NumType> const& miller_index) const
      {
        return d_star_sq_as_stol(d_star_sq(miller_index));
      }

      template <typename NumType>
      af::shared<double>
      stol(af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        af::shared<double> result(
          miller_indices.size(), af::init_functor_null<double>());
        for(std::size_t i=0;i<miller_indices.size();i++) {
          result[i] = stol(miller_indices[i]);
        }
        return result;
      }

      template <typename NumType>
      double
      d(miller::index<NumType> const& miller_index) const
      {
        return d_star_sq_as_d(d_star_sq(miller_index));
      }

      template <typename NumType>
      af::shared<double>
      d(af::const_ref<miller::index<NumType> > const& miller_indices) const
      {
        af::shared<double> result(
          miller_indices.size(), af::init_functor_null<double>());
        for(std::size_t i=0;i<miller_indices.size();i++) {
          result[i] = d(miller_indices[i]);
        }
        return result;
      }

      template <typename NumType>
      double
      two_theta(
        miller::index<NumType> const& miller_index,
        double wavelength,
        bool deg=false) const
      {
        return d_star_sq_as_two_theta(
          d_star_sq(miller_index), wavelength, deg);
      }

      template <typename NumType>
      af::shared<double>
      two_theta(
        af::const_ref<miller::index<NumType> > const& miller_indices,
        double wavelength,
        bool deg=false) const
      {
        af::shared<double> result(
          miller_indices.size(), af::init_functor_null<double>());
        for(std::size_t i=0;i<miller_indices.size();i++) {
          result[i] = two_theta(miller_indices[i], wavelength, deg);
        }
        return result;
      }


      double
      bases_mean_square_difference(unit_cell const& other) const
      {
        return mean_square_difference(orth_, other.orth_) / 3;
      }

      int
      compare_orthorhombic(
        const unit_cell& other) const;

      int
      compare_monoclinic(
        const unit_cell& other,
        unsigned unique_axis,
        double angular_tolerance) const;

      sgtbx::change_of_basis_op const&
      change_of_basis_op_for_best_monoclinic_beta() const;

    protected:
      void init_volume();
      void init_reciprocal();
      void init_metrical_matrices();
      void init_orth_and_frac_matrices();
      void initialize();

      af::double6 params_;
      af::double3 sin_ang_;
      af::double3 cos_ang_;
      double volume_;
      uc_sym_mat3 metr_mx_;
      af::double6 r_params_;
      af::double3 r_sin_ang_;
      af::double3 r_cos_ang_;
      uc_sym_mat3 r_metr_mx_;
      uc_mat3 frac_;
      uc_mat3 orth_;
      mutable double longest_vector_sq_;
      mutable double shortest_vector_sq_;

      // used by reciprocal()
      unit_cell(
        af::double6 const& params,
        af::double3 const& sin_ang,
        af::double3 const& cos_ang,
        double volume,
        uc_sym_mat3 const& metr_mx,
        af::double6 const& r_params,
        af::double3 const& r_sin_ang,
        af::double3 const& r_cos_ang,
        uc_sym_mat3 const& r_metr_mx);
  };

  template <typename FloatType>
  class incremental_d_star_sq
  {
    public:
      incremental_d_star_sq() {}


      incremental_d_star_sq(unit_cell const& ucell)
      {
        initialize(ucell.reciprocal_metrical_matrix());
      }

      void update0(int h0)
      {
        h0_ = h0;
        im0_ = (h0_ * h0_) * r_g00_;
      }

      void update1(int h1)
      {
        h1_ = h1;
        im1_ = im0_ + (h1_ * h1_) * r_g11_
                    + (2 * h0_ * h1_) * r_g01_;
      }

      FloatType get(int h2)
      {
        return im1_ + (h2 * h2) * r_g22_
                    + (2 * h0_ * h2) * r_g02_
                    + (2 * h1_ * h2) * r_g12_;
      }

    protected:
      FloatType r_g00_, r_g11_, r_g22_, r_g01_, r_g02_, r_g12_;
      int h0_, h1_;
      FloatType im0_, im1_;

      void initialize(uc_sym_mat3 const& r_g)
      {
        r_g00_ = r_g[0];
        r_g11_ = r_g[1];
        r_g22_ = r_g[2];
        r_g01_ = r_g[3];
        r_g02_ = r_g[4];
        r_g12_ = r_g[5];
      }
  };

}} // namespace cctbx::uctbx

#endif // CCTBX_UCTBX_H

---