cctbx/geometry_restraints/dihedral.h

#ifndef CCTBX_GEOMETRY_RESTRAINTS_DIHEDRAL_H
#define CCTBX_GEOMETRY_RESTRAINTS_DIHEDRAL_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 dihedral_proxy
  {
    typedef af::tiny<unsigned, 4> i_seqs_type;

    dihedral_proxy() {}

    dihedral_proxy(
      i_seqs_type const& i_seqs_,
      double angle_ideal_,
      double weight_,
      int periodicity_=0)
    :
      angle_ideal(angle_ideal_),
      weight(weight_),
      periodicity(periodicity_),
      i_seqs(i_seqs_)
    {}

    dihedral_proxy(
      i_seqs_type const& i_seqs_,
      dihedral_proxy const& params)
    :
      angle_ideal(params.angle_ideal),
      weight(params.weight),
      periodicity(params.periodicity),
      i_seqs(i_seqs_)
    {}

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

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

    dihedral_proxy
    sort_i_seqs() const
    {
      dihedral_proxy result(*this);
      if (result.i_seqs[0] > result.i_seqs[3]) {
        std::swap(result.i_seqs[0], result.i_seqs[3]);
        if ( sym_ops.get() ) {
          std::swap(result.sym_ops[0], result.sym_ops[3]);
        }
        result.angle_ideal *= -1;
      }
      if (result.i_seqs[1] > result.i_seqs[2]) {
        std::swap(result.i_seqs[1], result.i_seqs[2]);
        if ( sym_ops.get() ) {
          std::swap(result.sym_ops[1], result.sym_ops[2]);
        }
        result.angle_ideal *= -1;
      }
      return result;
    }

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


  class dihedral
  {
    public:
      dihedral() {}

      dihedral(
        af::tiny<scitbx::vec3<double>, 4> const& sites_,
        double angle_ideal_,
        double weight_,
        int periodicity_=0)
      :
        sites(sites_),
        angle_ideal(angle_ideal_),
        weight(weight_),
        periodicity(periodicity_)
      {
        init_angle_model();
      }

      dihedral(
        af::const_ref<scitbx::vec3<double> > const& sites_cart,
        dihedral_proxy const& proxy)
      :
        angle_ideal(proxy.angle_ideal),
        weight(proxy.weight),
        periodicity(proxy.periodicity)
      {
        for(int i=0;i<4;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();
      }

      dihedral(
        uctbx::unit_cell const& unit_cell,
        af::const_ref<scitbx::vec3<double> > const& sites_cart,
        dihedral_proxy const& proxy)
      :
        angle_ideal(proxy.angle_ideal),
        weight(proxy.weight),
        periodicity(proxy.periodicity)
      {
        for(int i=0;i<4;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
      {
        using scitbx::constants::pi_180;
        double term;
        if (periodicity > 0) {
          term = 9600. / (periodicity * periodicity)
               * (1 - std::cos(periodicity * delta * pi_180));
        }
        else {
          term = delta * delta;
        }
        return weight * term;
      }


      af::tiny<scitbx::vec3<double>, 4>
      gradients(double epsilon=1.e-100) const
      {
        af::tiny<scitbx::vec3<double>, 4> result;
        double d_21_norm = d_21.length_sq();
        if (   !have_angle_model
            || d_21_norm < epsilon
            || n_0121_norm < epsilon
            || n_2123_norm < epsilon) {
          result.fill(scitbx::vec3<double>(0,0,0));
        }
        else {
          using scitbx::constants::pi_180;
          double grad_factor;
          if (periodicity > 0) {
            grad_factor = 9600. / periodicity * pi_180
                        * std::sin(periodicity * delta * pi_180);
          }
          else {
            grad_factor = 2 * delta;
          }
          grad_factor *= weight * d_21.length() / pi_180;
          result[0] = -grad_factor/n_0121_norm * n_0121;
          result[3] = grad_factor/n_2123_norm * n_2123;
          double d_01_dot_d_21 = d_01 * d_21;
          double d_21_dot_d_23 = d_21 * d_23;
          result[1] = (d_01_dot_d_21/d_21_norm-1) * result[0]
                    - d_21_dot_d_23/d_21_norm * result[3];
          result[2] = (d_21_dot_d_23/d_21_norm-1) * result[3]
                    - d_01_dot_d_21/d_21_norm * result[0];
        }
        return result;
      }


      void
      add_gradients(
        af::ref<scitbx::vec3<double> > const& gradient_array,
        dihedral_proxy::i_seqs_type const& i_seqs) const
      {
        af::tiny<scitbx::vec3<double>, 4> grads = gradients();
        for(int i=0;i<4;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,
        dihedral_proxy const& proxy) const
      {
        dihedral_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>, 4> grads = gradients();
        for(int i=0;i<4;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>, 4> sites;
      double angle_ideal;
      double weight;
      int periodicity;
      bool have_angle_model;
    protected:
      scitbx::vec3<double> d_01;
      scitbx::vec3<double> d_21;
      scitbx::vec3<double> d_23;
      scitbx::vec3<double> n_0121;
      double n_0121_norm;
      scitbx::vec3<double> n_2123;
      double n_2123_norm;
    public:
      double angle_model;
      double delta;

    protected:
      void
      init_angle_model()
      {
        using scitbx::constants::pi_180;
        have_angle_model = false;
        angle_model = angle_ideal;
        delta = 0;
        d_01 = sites[0] - sites[1];
        d_21 = sites[2] - sites[1];
        d_23 = sites[2] - sites[3];
        n_0121 = d_01.cross(d_21);
        n_0121_norm = n_0121.length_sq();
        if (n_0121_norm == 0) return;
        n_2123 = d_21.cross(d_23);
        n_2123_norm = n_2123.length_sq();
        if (n_2123_norm == 0) return;
        double cos_angle_model = std::max(-1.,std::min(1.,
          n_0121 * n_2123 / std::sqrt(n_0121_norm * n_2123_norm)));
        angle_model = std::acos(cos_angle_model) / pi_180;
        if (d_21 * (n_0121.cross(n_2123)) < 0) {
          angle_model *= -1;
        }
        have_angle_model = true;
        delta = angle_delta_deg(angle_model, angle_ideal, periodicity);
      }
  };

  inline
  std::size_t
  dihedral_count_harmonic(
    af::const_ref<dihedral_proxy> const& proxies)
  {
    std::size_t result = 0;
    for(std::size_t i=0;i<proxies.size();i++) {
      if (proxies[i].periodicity <= 0) result++;
    }
    return result;
  }

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

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

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

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

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

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

}} // namespace cctbx::geometry_restraints

#endif // CCTBX_GEOMETRY_RESTRAINTS_DIHEDRAL_H

---