cctbx/geometry_restraints/angle.h

#ifndef CCTBX_GEOMETRY_RESTRAINTS_ANGLE_H
#define CCTBX_GEOMETRY_RESTRAINTS_ANGLE_H

#include <cctbx/sgtbx/rt_mx.h>
#include <cctbx/geometry_restraints/utils.h>
#include <scitbx/constants.h>
#include <scitbx/optional_copy.h>

namespace cctbx { namespace geometry_restraints {

  struct angle_proxy
  {
    typedef af::tiny<unsigned, 3> i_seqs_type;

    angle_proxy() {}

    angle_proxy(
      i_seqs_type const& i_seqs_,
      double angle_ideal_,
      double weight_)
    :
      i_seqs(i_seqs_),
      angle_ideal(angle_ideal_),
      weight(weight_)
    {}

    angle_proxy(
      i_seqs_type const& i_seqs_,
      af::shared<sgtbx::rt_mx> const& sym_ops_,
      double angle_ideal_,
      double weight_)
    :
      i_seqs(i_seqs_),
      sym_ops(sym_ops_),
      angle_ideal(angle_ideal_),
      weight(weight_)
    {
      if ( sym_ops.get() ) {
        CCTBX_ASSERT(sym_ops.get()->size() == i_seqs.size());
      }
    }


    angle_proxy(
      i_seqs_type const& i_seqs_,
      angle_proxy const& proxy)
    :
      i_seqs(i_seqs_),
      angle_ideal(proxy.angle_ideal),
      weight(proxy.weight)
    {}


    angle_proxy(
      i_seqs_type const& i_seqs_,
      af::shared<sgtbx::rt_mx> const& sym_ops_,
      angle_proxy const& proxy)
    :
      i_seqs(i_seqs_),
      sym_ops(sym_ops_),
      angle_ideal(proxy.angle_ideal),
      weight(proxy.weight)
    {}

    angle_proxy
    sort_i_seqs() const
    {
      angle_proxy result(*this);
      if (result.i_seqs[0] > result.i_seqs[2]) {
        std::swap(result.i_seqs[0], result.i_seqs[2]);
        if ( sym_ops.get() ) {
          std::swap(result.sym_ops[0], result.sym_ops[2]);
        }
      }
      return result;
    }

    i_seqs_type i_seqs;
    scitbx::optional_copy<af::shared<sgtbx::rt_mx> > sym_ops;
    double angle_ideal;
    double weight;
  };

  class angle
  {
    public:
      angle() {}

      angle(
        af::tiny<scitbx::vec3<double>, 3> const& sites_,
        double angle_ideal_,
        double weight_)
      :
        sites(sites_),
        angle_ideal(angle_ideal_),
        weight(weight_)
      {
        init_angle_model();
      }

      angle(
        af::const_ref<scitbx::vec3<double> > const& sites_cart,
        angle_proxy const& proxy)
      :
        angle_ideal(proxy.angle_ideal),
        weight(proxy.weight)
      {
        for(int i=0;i<3;i++) {
          std::size_t i_seq = proxy.i_seqs[i];
          CCTBX_ASSERT(i_seq < sites_cart.size());
          sites[i] = sites_cart[i_seq];
        }
        init_angle_model();
      }

      angle(
        uctbx::unit_cell const& unit_cell,
        af::const_ref<scitbx::vec3<double> > const& sites_cart,
        angle_proxy const& proxy)
      :
        angle_ideal(proxy.angle_ideal),
        weight(proxy.weight)
      {
        for(int i=0;i<3;i++) {
          std::size_t i_seq = proxy.i_seqs[i];
          CCTBX_ASSERT(i_seq < sites_cart.size());
          sites[i] = sites_cart[i_seq];
          if ( proxy.sym_ops.get() ) {
            sgtbx::rt_mx rt_mx = proxy.sym_ops[i];
            if ( !rt_mx.is_unit_mx() ) {
              sites[i] = unit_cell.orthogonalize(
                rt_mx * unit_cell.fractionalize(sites[i]));
            }
          }
        }
        init_angle_model();
      }


      double
      residual() const { return weight * scitbx::fn::pow2(delta); }


      void
      grads_and_curvs_impl(
        scitbx::vec3<double>* grads,
        scitbx::vec3<double>* curvs,
        double epsilon=1.e-100) const
      {
        if (!have_angle_model) {
          std::fill_n(grads, 3U, scitbx::vec3<double>(0,0,0));
          if (curvs != 0) {
            std::fill_n(curvs, 3U, scitbx::vec3<double>(0,0,0));
          }
        }
        else {
          double
          sin_angle_model = std::sqrt(1-scitbx::fn::pow2(cos_angle_model));
          if (sin_angle_model < epsilon) {
            std::fill_n(grads, 3U, scitbx::vec3<double>(0,0,0));
            if (curvs != 0) {
              std::fill_n(curvs, 3U, scitbx::vec3<double>(0,0,0));
            }
          }
          else {
            using scitbx::constants::pi_180;
            scitbx::vec3<double>
                d_angle_d_site0, d_angle_d_site1, d_angle_d_site2;
            double grad_factor = -2 * weight * delta / pi_180;
            d_angle_d_site0 = (d_01_unit * cos_angle_model - d_21_unit) /
                              (sin_angle_model * d_01_abs);
            grads[0] = grad_factor * d_angle_d_site0;
            d_angle_d_site2 = (d_21_unit * cos_angle_model - d_01_unit) /
                              (sin_angle_model * d_21_abs);
            grads[2] = grad_factor * d_angle_d_site2;
            grads[1] = -(grads[0] + grads[2]);
            if (curvs != 0)
            {
              d_angle_d_site1 = grads[1] / grad_factor;
              scitbx::vec3<double> v111(1,1,1);
              double sinsqr = scitbx::fn::pow2(sin_angle_model);

              scitbx::vec3<double> d2_angle_d_site00 =
                (2. * d_21_unit.each_mul(d_01_unit) +
                cos_angle_model * (v111 * sinsqr - d_21_unit.each_mul(d_21_unit)
                - (1. + 2. * sinsqr) * d_01_unit.each_mul(d_01_unit))) /
                (scitbx::fn::pow2(d_01_abs)*sin_angle_model*sinsqr);

              scitbx::vec3<double> d2_angle_d_site22 =
                (2. * d_01_unit.each_mul(d_21_unit) +
                cos_angle_model * (v111 * sinsqr - d_01_unit.each_mul(d_01_unit)
                - (1. + 2. * sinsqr) * d_21_unit.each_mul(d_21_unit))) /
                (scitbx::fn::pow2(d_21_abs)*sin_angle_model*sinsqr);

               double d01sqr = scitbx::fn::pow2(d_01_abs);
               double d21sqr = scitbx::fn::pow2(d_21_abs);
               scitbx::vec3<double> term1 = d_angle_d_site1 / sin_angle_model;
               scitbx::vec3<double> tvec1 = d_01_abs * d_21_unit;
               scitbx::vec3<double> tvec2 = d_21_abs * d_01_unit;
               scitbx::vec3<double> sumvec = d_01 + d_21;
               scitbx::vec3<double> d2_angle_d_site11 =
                 (-2. * cos_angle_model *
                 (tvec1.each_mul(tvec1) + tvec2.each_mul(tvec2)) +
                 (tvec1+tvec2).each_mul(sumvec) +
                 (d01sqr + d21sqr) * v111 * cos_angle_model +
                 (d01sqr * d_21 + d21sqr * d_01).each_mul(term1) -
                 d_01_abs * d_21_abs *
                 (2. * v111 + cos_angle_model * term1.each_mul(sumvec))) /
                 (sin_angle_model * d01sqr * d21sqr);

               double curvfac = 2. * weight / pi_180;
               curvs[0] = curvfac *
                 (d_angle_d_site0.each_mul(d_angle_d_site0) / pi_180 -
                 delta * d2_angle_d_site00);
               curvs[1] = curvfac *
                 (d_angle_d_site1.each_mul(d_angle_d_site1) / pi_180 -
                 delta * d2_angle_d_site11);
               curvs[2] = curvfac *
                 (d_angle_d_site2.each_mul(d_angle_d_site2) / pi_180 -
                 delta * d2_angle_d_site22);
            }
          }
        }
      }


      af::shared<scitbx::vec3<double> >
      grads_and_curvs(double epsilon=1.e-100) const
      {
        af::shared<scitbx::vec3<double> > result(6);
        grads_and_curvs_impl(&result[0], &result[3], epsilon);
        return result;
      }


      af::tiny<scitbx::vec3<double>, 3>
      gradients(double epsilon=1.e-100) const
      {
        af::tiny<scitbx::vec3<double>, 3> result;
        grads_and_curvs_impl(result.begin(), 0, epsilon);
        return result;
      }


      void
      add_gradients(
        af::ref<scitbx::vec3<double> > const& gradient_array,
        angle_proxy::i_seqs_type const& i_seqs) const
      {
        af::tiny<scitbx::vec3<double>, 3> grads;
        grads_and_curvs_impl(grads.begin(), 0);
        for(int i=0;i<3;i++) {
          gradient_array[i_seqs[i]] += grads[i];
        }
      }


      void
      add_gradients(
        uctbx::unit_cell const& unit_cell,
        af::ref<scitbx::vec3<double> > const& gradient_array,
        angle_proxy const& proxy) const
      {
        angle_proxy::i_seqs_type const& i_seqs = proxy.i_seqs;
        scitbx::optional_copy<af::shared<sgtbx::rt_mx> > const&
          sym_ops = proxy.sym_ops;
        af::tiny<scitbx::vec3<double>, 3> grads;
        grads_and_curvs_impl(grads.begin(), 0);
        for(int i=0;i<3;i++) {
          if ( sym_ops.get() && !sym_ops[i].is_unit_mx() ) {
            scitbx::mat3<double> r_inv_cart_ = r_inv_cart(
              unit_cell, sym_ops[i]);
            gradient_array[i_seqs[i]] += grads[i] * r_inv_cart_;
          }
          else { gradient_array[i_seqs[i]] += grads[i]; }
        }
      }

      af::tiny<scitbx::vec3<double>, 3> sites;
      double angle_ideal;
      double weight;
      bool have_angle_model;
    protected:
      double d_01_abs;
      double d_21_abs;
      scitbx::vec3<double> d_01;
      scitbx::vec3<double> d_21;
      scitbx::vec3<double> d_01_unit;
      scitbx::vec3<double> d_21_unit;
      double cos_angle_model;
    public:
      double angle_model;
      double delta;

    protected:
      void
      init_angle_model()
      {
        have_angle_model = false;
        d_01_abs = 0;
        d_21_abs = 0;
        d_01.fill(0);
        d_21.fill(0);
        d_01_unit.fill(0);
        d_21_unit.fill(0);
        cos_angle_model = -9;
        angle_model = angle_ideal;
        delta = 0;
        d_01 = sites[0] - sites[1];
        d_01_abs = d_01.length();
        if (d_01_abs > 0) {
          d_21 = sites[2] - sites[1];
          d_21_abs = d_21.length();
          if (d_21_abs > 0) {
            d_01_unit = d_01 / d_01_abs;
            d_21_unit = d_21 / d_21_abs;
            cos_angle_model = std::max(-1.,std::min(1.,d_01_unit*d_21_unit));
            angle_model = std::acos(cos_angle_model)
                        / scitbx::constants::pi_180;
            have_angle_model = true;
            delta = angle_delta_deg(angle_model, angle_ideal);
          }
        }
      }
  };

  inline
  af::shared<double>
  angle_deltas(
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<angle_proxy> const& proxies)
  {
    return detail::generic_deltas<angle_proxy, angle>::get(
      sites_cart, proxies);
  }

  inline
  af::shared<double>
  angle_deltas(
    uctbx::unit_cell const& unit_cell,
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<angle_proxy> const& proxies)
  {
    return detail::generic_deltas<angle_proxy, angle>::get(
      unit_cell, sites_cart, proxies);
  }

  inline
  af::shared<double>
  angle_residuals(
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<angle_proxy> const& proxies)
  {
    return detail::generic_residuals<angle_proxy, angle>::get(
      sites_cart, proxies);
  }

  inline
  af::shared<double>
  angle_residuals(
    uctbx::unit_cell const& unit_cell,
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<angle_proxy> const& proxies)
  {
    return detail::generic_residuals<angle_proxy, angle>::get(
      unit_cell, sites_cart, proxies);
  }

  inline
  double
  angle_residual_sum(
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<angle_proxy> const& proxies,
    af::ref<scitbx::vec3<double> > const& gradient_array)
  {
    return detail::generic_residual_sum<angle_proxy, angle>::get(
      sites_cart, proxies, gradient_array);
  }

  inline
  double
  angle_residual_sum(
    uctbx::unit_cell const& unit_cell,
    af::const_ref<scitbx::vec3<double> > const& sites_cart,
    af::const_ref<angle_proxy> const& proxies,
    af::ref<scitbx::vec3<double> > const& gradient_array)
  {
    return detail::generic_residual_sum<angle_proxy, angle>::get(
      unit_cell, sites_cart, proxies, gradient_array);
  }

}} // namespace cctbx::geometry_restraints

#endif // CCTBX_GEOMETRY_RESTRAINTS_ANGLE_H

---