rmsd.hpp

Go to the documentation of this file.

// -*- mode: c++; indent-tabs-mode: nil; -*-
//
//The Biomolecule Toolkit (BTK) is a C++ library for use in the
//modeling, analysis, and design of biological macromolecules.
//Copyright (C) 2001-2006, Chris Saunders <ctsa@users.sourceforge.net>,
//                         Tim Robertson <kid50@users.sourceforge.net>
//
//This program is free software; you can redistribute it and/or modify
//it under the terms of the GNU Lesser General Public License as published
//by the Free Software Foundation; either version 2.1 of the License, or (at
//your option) any later version.
//
//This program is distributed in the hope that it will be useful,  but
//WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//Lesser General Public License for more details.
//
//You should have received a copy of the GNU Lesser General Public License
//along with this program; if not, write to the Free Software Foundation,
//Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


#ifndef BTK_ALGORITHMS_RMSD_HPP
#define BTK_ALGORITHMS_RMSD_HPP

#include <cmath>
#include <functional>
#include <algorithm>

#include <btk/core/math/btk_vector.hpp>
#include <btk/core/math/btk_matrix.hpp>
#include <btk/core/math/linear_algebra.hpp>

namespace BTK {
namespace ALGORITHMS {


#ifndef DOXYGEN_SHOULD_SKIP_THIS
namespace internal {

// This routine will calculate the eigenvectors corresponding
// to pre-calculated eigenvalues for use in RMSD calculations.
// It is not intended to be a solution to the full eigenvector
// calculation problem for 3D matrices.
//
// I think this will work for any 3D matrix that produces three
// real eigenvalues, but I make no guarantees. For RMSD calculations,
// it will work b/c we know that we are dealing with three orthogonal
// real vectors.  If you can't make that assumption
// in your application, don't use this code!
//
// The code checks for eigenvalues with multiplicity > 1, but the
// method used in this case may not be correct, so be careful (as I've
// already said, this won't be the case for 3D structural data).
//
void fast_solve_3D_eigenvectors(BTK::MATH::BTKMatrix & M,
                                double const evalues[3],
                                BTK::MATH::BTKVector evecs[3]);

// Solve eigenvalues (ev) for a 3x3 matrix M quickly by analytically
// finding roots of the characteristic equation: det(M - evI) = 0.
//
// It is possible that the factorization used here
// is numerically unstable (but it likely doesn't matter for RMSD
// calculations).  At the present time, there are no checks or
// warnings on the quality of the result.
void fast_solve_3D_eigenvalues(BTK::MATH::BTKMatrix const & M,
                               double e_vals[3]);


// Solves the determinant for a 3x3 matrix -- not general!
inline double determinant(BTK::MATH::BTKMatrix const & M)
{
  using boost::numeric::ublas::inner_prod;
  using BTK::MATH::cross;
  return inner_prod(row(M,0),cross(row(M,1),row(M,2)));
}

template <typename AtomIterator1, typename AtomIterator2>
 void
 setup_rmsd(AtomIterator1 first1,
            AtomIterator1 last1,
            AtomIterator2 first2,
            AtomIterator2 last2,
            BTK::MATH::BTKVector & T1,
            BTK::MATH::BTKVector & T2,
            BTK::MATH::BTKMatrix & R,
            BTK::MATH::BTKMatrix & RtR,
            double & Eo,
            unsigned & N,
            double & omega,
            double e_vals[3])
{
  // Get T1, the transform from the moved set COM to the origin
  //   & T2, the transform from the origin to the fixed set COM
  T1 = BTK::MATH::BTKVector(0.0);
  T2 = BTK::MATH::BTKVector(0.0);
  AtomIterator1 i1(first1);
  AtomIterator2 i2(first2);

  for (N = 0;
       i1 != last1 && i2 != last2;
       ++i1,++i2) {
    T1 += i2->position();
    T2 += i1->position();
    N++;
  }

  if (i2 != last2 || i1 != last1) {
    std::cerr << "Warning: attempting to compute RMSD between "
              << "structures of different lengths! " << std::endl
              << "Assuming atom length of smaller structure (" << N << ")...."
              << std::endl;
  }

  T1 /= -1.*static_cast<double>(N);
  T2 /= static_cast<double>(N);

  //
  // Get Eo and R
  //
  int i = 0;
  BTK::MATH::BTKMatrix tmp1(3,N),tmp2(3,N);
  for (i1 = first1, i2 = first2;
       i1 != last1 && i2 != last2;
       ++i1,++i2) {
    column(tmp1,i) = i1->position() - T2;
    column(tmp2,i) = i2->position() + T1;
    i++;
  }

  Eo = 0;
  for (unsigned j = 0; j < 3; ++j) {
    Eo += inner_prod(row(tmp1,j),row(tmp1,j)) +
      inner_prod(row(tmp2,j),row(tmp2,j));
    for (unsigned k = 0; k < 3; ++k) {
      R(j,k) = inner_prod(row(tmp1,j),row(tmp2,k));
    }
  }
  Eo /= 2.0;

  // if det(R) < 0, the rotation is degenerate (a flip).
  if (determinant(R) > 0) omega = 1;
  else omega = -1;

  RtR = prec_prod(trans(R),R);

  internal::fast_solve_3D_eigenvalues(RtR,e_vals);

  // The eigenvalues of R (what we really want) are equal to the square root
  // of the eigenvalues of RtR. The fast eigenvalue routine is somewhat 
  // unstable, and the if statements correct for small-magnitude negative
  // values that will occasionally come out of the routine when the true
  // eigenvalue is zero (b/c RtR is a positive definite symmetric real matrix
  // we know in advance that its eigenvalues will always be real and posiive, 
  // so this correction is always mathematically valid.)
  for (unsigned i = 0; i < 3; ++i) {
    if (e_vals[i] > 0) e_vals[i] = std::sqrt(e_vals[i]);
    else e_vals[i] = 0.0;
  }

  std::sort(e_vals,e_vals + 3,
            std::greater<double>());
}

} // namespace internal
#endif // DOXYGEN_SHOULD_SKIP_THIS

template <typename AtomIterator1, typename AtomIterator2>
double
leastsquares_superposition(AtomIterator1 fixed_first,
                           AtomIterator1 fixed_last,
                           AtomIterator2 superimpose_first,
                           AtomIterator2 superimpose_last,
                           BTK::MATH::BTKMatrix & U,
                           BTK::MATH::BTKVector & T)
 {
   double Eo, omega;
   double e_vals[3];
   BTK::MATH::BTKMatrix R,RtR;
   BTK::MATH::BTKVector T1;
   BTK::MATH::BTKVector T2;
   unsigned N;

   internal::setup_rmsd(fixed_first,
                        fixed_last,
                        superimpose_first,
                        superimpose_last,
                        T1,T2,R,RtR,Eo,N,
                        omega,e_vals);

   BTK::MATH::BTKVector right_evecs[3];

   // In order to create the rotation matrix, we need the
   // eigenvectors of the RtR matrix. To do this, we need
   // the eigenvalues of this matrix.  These are the
   // squares of the eigenvalues of R.
   e_vals[0] *= e_vals[0];
   e_vals[1] *= e_vals[1];
   e_vals[2] *= e_vals[2];

   internal::fast_solve_3D_eigenvectors(RtR,e_vals,right_evecs);

   // Useful algorithm information:
   // The eigenvectors of RtR are equivalent to the right eigenvectors of
   // matrix R. We need the left eigenvectors of R, which are determined by
   // Left = R * Right.
   BTK::MATH::BTKVector left_evecs[3];

   left_evecs[0] = prec_prod(R,right_evecs[0]);
   left_evecs[1] = prec_prod(R,right_evecs[1]);

   // normalize the vectors
   left_evecs[0] /= norm_2(left_evecs[0]);
   left_evecs[1] /= norm_2(left_evecs[1]);

   // in order to guarantee a right-handed system, we'll
   // determine the Z vector as the cross of X and Y.
   left_evecs[2] = cross(left_evecs[0],left_evecs[1]);

   // Construct the optimal rotation matrix.
   //
   // More algorithm details:
   // The optimal rotation matrix is the product of matrices
   // containing the left and right eigenvectors of RtR.
   //
   // Said another way: if  RtR = BSV
   //  (this is the result of singular value decomposition on RtR)
   // the optimal rotation is U = B * t(V)
   //
   U.clear();
   for (unsigned i = 0; i < 3; ++i) {
     for (unsigned j = 0; j < 3; ++j) {
       U(i,j) = (left_evecs[0][i]*right_evecs[0][j]) +
                (left_evecs[1][i]*right_evecs[1][j]) +
                (left_evecs[2][i]*right_evecs[2][j]);
     }
   }

   // construct the post-rotation translation vector
   T = prec_prod(U,T1)+T2;


   // compute the rmsd between the aligned structures.
   double rmsd = (2.0/static_cast<double>(N))*
     (Eo - e_vals[0] - e_vals[1] - omega*e_vals[2]);

   double const sqrt_min = 1e-8;
   if(rmsd > sqrt_min) { rmsd = std::sqrt(rmsd); }
   else                { rmsd = 0; }

   return rmsd;
}

template <typename AtomIterator1, typename AtomIterator2, typename AtomIterator3>
double
leastsquares_superposition(AtomIterator1 fixed_first,
                           AtomIterator1 fixed_last,
                           AtomIterator2 superimpose_first,
                           AtomIterator2 superimpose_last,
                           AtomIterator3 transform_first,
                           AtomIterator3 transform_last)
{
  BTK::MATH::BTKVector T;
  BTK::MATH::BTKMatrix U;
  double rmsd =
    leastsquares_superposition(fixed_first,fixed_last,
                               superimpose_first,superimpose_last,
                               U,T);

  for(AtomIterator3 ai=transform_first; ai != transform_last; ++ai)
    ai->set_position(prec_prod(U,ai->position())+T);

  return rmsd;
}

template <typename AtomIterator1, typename AtomIterator2>
double
get_rmsd(AtomIterator1 a1_first,
         AtomIterator1 a1_last,
         AtomIterator2 a2_first,
         AtomIterator2 a2_last)
{
  double Eo, omega;
  double e_vals[3];
  BTK::MATH::BTKMatrix R,RtR;
  unsigned N;
  BTK::MATH::BTKVector T1,T2;

  internal::setup_rmsd(a1_first,a1_last,a2_first,a2_last,
                       T1,T2,R,RtR,Eo,N,omega,e_vals);

  double rmsd = (2.0/static_cast<double>(N))*
    (Eo - e_vals[0] - e_vals[1] - omega*e_vals[2]);

  double const sqrt_min = 1e-8;
  if(rmsd > sqrt_min) { rmsd = std::sqrt(rmsd); }
  else                { rmsd = 0; }

  return rmsd;
}

template <typename AtomIterator1, typename AtomIterator2>
double
get_trivial_rmsd(AtomIterator1 a1_first,
                 AtomIterator1 a1_last,
                 AtomIterator2 a2_first,
                 AtomIterator2 a2_last)
{
  double rmsd2 = 0.;
  unsigned N = 0;

  AtomIterator1 a1 = a1_first;
  AtomIterator2 a2 = a2_first;

  for(;a1!=a1_last && a2!=a2_last;++a1,++a2){
    BTK::MATH::BTKVector v(a1->position()-a2->position());
    rmsd2 += inner_prod(v,v);
    N++;
  }

  if(a1!=a1_last || a2!=a2_last) {
    std::cerr << "Warning: attempting to compute RMSD between "
              << "structures of different lengths! " << std::endl
              << "Assuming atom length of smaller structure (" << N << ")...."
              << std::endl;
  }

  if(N) rmsd2 /= static_cast<double>(N);

  double const sqrt_min = 1e-8;
  if(rmsd2 > sqrt_min) { rmsd2 = std::sqrt(rmsd2); }
  else                 { rmsd2 = 0; }

  return rmsd2;
}

} // ALGORITHMS
} // BTK

#endif

---