atom_algorithms.hpp

Go to the documentation of this file.

// -*- mode: c++; -*-
//
//The Biomolecule Toolkit (BTK) is a C++ library for use in the
//modeling, analysis, and design of biological macromolecules.
//Copyright (C) 2001, 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.



#include <cmath>

#include <iterator>

namespace BTK{
  namespace Internal{

    template <typename AtomIterator1,typename AtomIterator2>
    void
    setup_rmsd(AtomIterator1 first1,
               AtomIterator1 last1,
               AtomIterator2 first2,
               AtomIterator2 last2,
               BTKVector& T1,
               BTKVector& T2,
               BTKMatrix& R,
               double& Eo,
               unsigned& N);
  }
}


template <typename AtomIterator>
BTK::BTKVector 
BTK::
center_of_mass(AtomIterator first,
               AtomIterator last)
{
  typedef typename std::iterator_traits<AtomIterator>::value_type InputT;
  boost::function_requires<boost::ConvertibleConcept<InputT,Atom> >();

  BTKVector COM = new_vector();
  for(AtomIterator i=first;i!=last;++i){ COM += i->position(); }
  
  return COM/std::distance(first,last);
}


template <typename AtomIterator>
void 
BTK::
center(AtomIterator first,
       AtomIterator last) 
{
  BTKVector COM(center_of_mass(first,last));

  translate(first,last,-COM);
}


template <typename AtomIterator1,typename AtomIterator2>
void 
BTK::
align_centers(AtomIterator1 fixed_first,
              AtomIterator1 fixed_last,
              AtomIterator2 align_first,
              AtomIterator2 align_last)
{
  BTKVector fixedCOM(center_of_mass(fixed_first,fixed_last));
  BTKVector COM(center_of_mass(align_first,align_last));
  translate(align_first,align_last,fixedCOM-COM);
}



template <typename AtomIterator>
void 
BTK::
translate(AtomIterator first,
          AtomIterator last,
          const BTKVector& T) 
{
  typedef typename std::iterator_traits<AtomIterator>::value_type InputT;
  boost::function_requires<boost::ConvertibleConcept<InputT,Atom> >();

  for (;first!=last;++first) {
    first->set_position(first->position()+T);
  }
}


template <typename AtomIterator>
void 
BTK::
rotate(AtomIterator first,
       AtomIterator last,
       double phi,
       double theta,
       double psi,
       const BTKVector& fixed_point)
{
  BTKMatrix R(create_rotation_matrix(phi,theta,psi));
  BTKVector T(prec_prod(R,-fixed_point)+fixed_point);  

  for(; first != last; ++first) {
    first->set_position(prec_prod(R,first->position())+T);
  }
}


template <typename AtomIterator>
void 
BTK::
rotate(AtomIterator first,
       AtomIterator last,
       const BTKVector& axis,
       double theta,
       const BTKVector& fixed_point)
{
  BTKMatrix R(create_rotation_matrix(axis,theta));
  BTKVector T(prec_prod(R,-fixed_point)+fixed_point);  

  for(; first != last; ++first) {
    first->set_position(prec_prod(R,first->position())+T);
  }
}


template <typename AtomIterator>
void 
BTK::
rotate(AtomIterator first,
       AtomIterator last,
       const BTKMatrix& R)
{
  for (; first != last; ++first) {
    first->set_position(prec_prod(R,first->position()));
  }
}


template <typename AtomIterator>
void
BTK::
transform(AtomIterator first,
          AtomIterator last,
          const BTKMatrix& R,
          const BTKVector& T)
{
  for(;first!=last;++first){
    first->set_position(prec_prod(R,first->position())+T);
  }
}



template <typename AtomIterator>
void 
BTK::
principal_axes(AtomIterator first,
               AtomIterator last,
               std::vector<EigenState>& axes,
               BTKVector& COM)
{
  COM = center_of_mass(first,last);

  BTKMatrix covar(3,3);
  covar.clear();

  for(AtomIterator i = first; i != last; ++i){
    BTKVector tmp(i->position() - COM);
    covar += outer_prod(tmp,tmp);
  }

  axes = diagonalize_symmetric(covar);
}


template <typename AtomIterator1, typename AtomIterator2>
void 
BTK::
Internal::
setup_rmsd(AtomIterator1 first1,
           AtomIterator1 last1,
           AtomIterator2 first2,
           AtomIterator2 last2,
           BTKVector& T1,
           BTKVector& T2,
           BTKMatrix& R,
           double& Eo,
           unsigned& N) 
{
  // Get T1, the transform from the moved set COM to the origin
  //   & T2, the transform from the origin to the fixed set COM
  T1 *= 0;
  T2 *= 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;
  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.;
}



// Fast calculation of rmsd w/o calculating a rotation matrix.
//
//   Chris Saunders 11/2002 - Fast rmsd calculation by the method of 
// Kabsch 1978, where the required eigenvalues are found by an 
// analytical, rather than iterative, method to save time. 
// The cubic factorization used to accomplish this only produces 
// stable eigenvalues for the transpose(R)*R matrix of a typical 
// protein after the whole matrix has been normalized. Note that 
// the normalization process used here is completely empirical 
// and that, at the present time, there are **no checks** or 
// warnings on the quality of the (potentially unstable) cubic 
// factorization. 
//
template <typename AtomIterator1, typename AtomIterator2>
double 
BTK::
get_rmsd(AtomIterator1 a1_first,
         AtomIterator1 a1_last,
         AtomIterator2 a2_first,
         AtomIterator2 a2_last)
{
  double Eo;
  BTKMatrix R(3,3);
  unsigned N;
  BTKVector T1(3),T2(3);
  Internal::setup_rmsd(a1_first,a1_last,a2_first,a2_last,T1,T2,R,Eo,N);

  // check if the determinant is greater than 0
  double omega = (inner_prod(row(R,0),cross(row(R,1),row(R,2))) > 0 ? 1 : -1);

  //  get elements we need from trans(R)*R 
  //   (funky matrix naming relic of first attempt using pivots)
  //          matrix = d0 e0 f0
  //                      d1 e1
  //                         d2
  double d0,d1,d2,e0,e1,f0;

  // divide matrix by d0, so that cubic root algorithm can handle it
  d0 = R(0,0)*R(0,0) + R(1,0)*R(1,0) + R(2,0)*R(2,0);
  d1 = (R(0,1)*R(0,1) + R(1,1)*R(1,1) + R(2,1)*R(2,1))/d0;
  d2 = (R(0,2)*R(0,2) + R(1,2)*R(1,2) + R(2,2)*R(2,2))/d0;

  e0 = (R(0,0)*R(0,1) + R(1,0)*R(1,1) + R(2,0)*R(2,1))/d0;
  e1 = (R(0,1)*R(0,2) + R(1,1)*R(1,2) + R(2,1)*R(2,2))/d0;

  f0 = (R(0,0)*R(0,2) + R(1,0)*R(1,2) + R(2,0)*R(2,2))/d0;

  // cubic roots
  double r1,r2,r3;
  {
    // solving for eigenvalues as det(A-I*lambda) = 0
    // yeilds the values below corresponding to:
    // lambda**3 + B*lambda**2 + C*lambda + D = 0
    //   (given that d0=1.)
    //
    double B = -1-d1-d2;
    double C = d1+d2+d1*d2-e0*e0-f0*f0-e1*e1;
    double D = e0*e0*d2+e1*e1+f0*f0*d1-d1*d2-2*e0*f0*e1;

    // cubic root method of Viete w/ all the safety belts off
    double q = (B*B-3.0*C)/9.0;
    double q3 = q*q*q;
    double r = (2.0*B*B*B-9.0*B*C+27.0*D)/54.0;
    double theta = acos(r/sqrt(q3));
    r1 = r2 = r3 = -2.0*sqrt(q);
    r1 *= cos(theta/3.0);
    r2 *= cos((theta+2*PI)/3.0);
    r3 *= cos((theta-2*PI)/3.0);
    r1 -= B/3.0;
    r2 -= B/3.0;
    r3 -= B/3.0;
  }

  // undo the d0 norm to get eigenvalues
  r1 = r1*d0;
  r2 = r2*d0;
  r3 = r3*d0;

  // set rlow to lowest eigenval; set other two to r1,r2 
  double rlow;
  if(r3<r1 && r3<r2){
    rlow = r3;
  } else if(r2<r1 && r2<r3){
    rlow = r2; r2=r3;
  } else { 
    rlow = r1; r1=r3;
  }
  double rmsd = (2./static_cast<double>(N))*
                (Eo - std::sqrt(r1) - std::sqrt(r2) - omega*std::sqrt(rlow));

  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 
BTK::
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){
    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;
}




template <typename AtomIterator1, typename AtomIterator2>
void 
BTK::
leastsquares_superposition(AtomIterator1 fixed_first,
                           AtomIterator1 fixed_last, 
                           AtomIterator2 superimpose_first,
                           AtomIterator2 superimpose_last,
                           BTKMatrix& U,
                           BTKVector& T,
                           double& rmsd)
{
  double Eo;
  BTKMatrix R(3,3);
  BTKVector T1(3);
  BTKVector T2(3);
  unsigned N;
  Internal::setup_rmsd(fixed_first,fixed_last,superimpose_first,superimpose_last,T1,T2,R,Eo,N);
  
  // 1) construct the rotation matrix
  
  // Diagonalize transpose(R)*R
  std::vector<EigenState> ES = diagonalize_symmetric(prec_prod(trans(R),R));

  // sort from highest to lowest eigenvalue
  std::sort(ES.rbegin(),ES.rend());
  

  // Set eigenvector 2 = 0 x 1 to guarantee a r.h. system
  ES[2].eigenvector = cross(ES[0].eigenvector,ES[1].eigenvector);

  double mu0,mu1,mu2;

  mu0 = std::sqrt(ES[0].eigenvalue);
  mu1 = std::sqrt(ES[1].eigenvalue);
  mu2 = std::sqrt(ES[2].eigenvalue);

  BTKVector b0(3),b1(3),b2(3);

  b0 = prec_prod(R,ES[0].eigenvector/mu0);
  b1 = prec_prod(R,ES[1].eigenvector/mu1);
  b2 = cross(b0,b1);

  // store omega to distinguish reflection from rotation for rmsd calc
  double omega = (inner_prod(b2,(prec_prod(R,ES[2].eigenvector))) < 0 ? -1 : 1);

  U.clear();
  for (unsigned i = 0; i < 3; ++i) {
    for (unsigned j = 0; j < 3; ++j) {
      U(i,j) = (b0[i]*ES[0].eigenvector[j]) + 
               (b1[i]*ES[1].eigenvector[j]) +
               (b2[i]*ES[2].eigenvector[j]);
    }
  }

  T = prec_prod(U,T1)+T2;
    
  // Eo-sum(muX) is 1/2 the sum of squared differences between atoms. 
  // To get rmsd, multiply this by 2 and divide by the number of points, 
  // then take the square root.
  rmsd = (2./static_cast<double>(N))*(Eo - mu0 - mu1 - omega*mu2);
  // added a check for the smallest value that can safely be sqrt'ed
  double const sqrt_min = 1.e-8;
  if(rmsd > sqrt_min) { rmsd = std::sqrt(rmsd); }
  else                { rmsd = 0; }


#ifndef NDEBUG

  // crap code for verification below here
  
  // check rot and trans
  //
  double rmsd2 = 0.;
  AtomIterator1 fi;
  AtomIterator2 mi;
  for(fi = fixed_first,mi = superimpose_first;
      fi != fixed_last && mi != superimpose_last;
      fi++,mi++){
    BTKVector foo(3);
    foo = fi->position()-(prec_prod(U,mi->position())+T);
    rmsd2 += inner_prod(foo,foo);  // = length of vector squared
  }
  rmsd2 = std::sqrt(rmsd2/static_cast<double>(N));

  std::cerr << "fastrmsd: " << rmsd 
            << " superimposed directrmsd: " << rmsd2 << std::endl;
#endif

}


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

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

---