src/math/RMSD.cpp

#include "math/RMSD.hpp"
#include "math/SVD.hpp"

using namespace oopse;
using namespace JAMA;

RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov,
                              Vector3d mov_com,
                              Vector3d mov_to_ref,
                              RotMat3x3d U) {

  assert(mov.size() == ref_.size());
  int n;
  int n_vec = ref_.size();

  /* calculate the centre of mass */
  mov_com = V3Zero;
  
  for (n=0; n < n_vec; n++) {
    mov_com += mov[n];
  }
  
  mov_com /= (RealType)n_vec;
  
  mov_to_ref = ref_com - mov_com;
  
  /* shift mov to center of mass */

  for (n=0; n < n_vec; n++) {
    mov[n] -= mov_com;
  }
  
  /* initialize */
  Mat3x3d R = Mat3x3d(0.0);
  RealType E0 = 0.0;
  
  for (n=0; n < n_vec; n++) {
    
    /* 
     * E0 = 1/2 * sum(over n): y(n)*y(n) + x(n)*x(n) 
     */
    E0 += dot(mov[n], mov[n]) + dot(ref_[n], ref_[n]);
    
    /*
     * correlation matrix R:   
     *   R(i,j) = sum(over n): y(n,i) * x(n,j)  
     *   where x(n) and y(n) are two vector sets   
     */

    R += outProduct(mov[n], ref_[n]);

  }
  E0 *= 0.5;

  RectMatrix<RealType, n_vec, 3> v;
  Vector3d s;
  Mat3x3d w;

  SVD<RealType, n_vec, 3> svd = SVD<RealType, n_vec, 3>(R);
  svd.getU(v);
  svd.getSingularValues(s);
  svd.getV(w);
    
  int is_reflection = (v.determinant() * w.determinant()) < 0.0;
  if (is_reflection)
    s(2) = -s(2);

  RealType rmsd_sq = (E0 - 2.0 * s.sum() )/ (RealType)n_vec;
  rmsd_sq = max(rmsd_sq,0.0);
  RealType rmsd = sqrt(rmsd_sq);
  return rmsd;
}
                            

RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov) {
}

Revision:	1335
Committed:	Fri Apr 3 17:44:57 2009 UTC (16 years, 1 month ago) by gezelter
File size:	1636 byte(s)
Log Message:	Beginning import of SVD and RMSD code to handle orientational restraints of flexible molecules. Work in progress, doesn't compile yet.
#	User	Rev	Content
1	gezelter	1335	#include "math/RMSD.hpp"
2			#include "math/SVD.hpp"
3
4			using namespace oopse;
5			using namespace JAMA;
6
7			RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov,
8			Vector3d mov_com,
9			Vector3d mov_to_ref,
10			RotMat3x3d U) {
11
12			assert(mov.size() == ref_.size());
13			int n;
14			int n_vec = ref_.size();
15
16			/* calculate the centre of mass */
17			mov_com = V3Zero;
18
19			for (n=0; n < n_vec; n++) {
20			mov_com += mov[n];
21			}
22
23			mov_com /= (RealType)n_vec;
24
25			mov_to_ref = ref_com - mov_com;
26
27			/* shift mov to center of mass */
28
29			for (n=0; n < n_vec; n++) {
30			mov[n] -= mov_com;
31			}
32
33			/* initialize */
34			Mat3x3d R = Mat3x3d(0.0);
35			RealType E0 = 0.0;
36
37			for (n=0; n < n_vec; n++) {
38
39			/*
40			* E0 = 1/2 * sum(over n): y(n)y(n) + x(n)x(n)
41			*/
42			E0 += dot(mov[n], mov[n]) + dot(ref_[n], ref_[n]);
43
44			/*
45			* correlation matrix R:
46			* R(i,j) = sum(over n): y(n,i) * x(n,j)
47			* where x(n) and y(n) are two vector sets
48			*/
49
50			R += outProduct(mov[n], ref_[n]);
51
52			}
53			E0 *= 0.5;
54
55			RectMatrix<RealType, n_vec, 3> v;
56			Vector3d s;
57			Mat3x3d w;
58
59			SVD<RealType, n_vec, 3> svd = SVD<RealType, n_vec, 3>(R);
60			svd.getU(v);
61			svd.getSingularValues(s);
62			svd.getV(w);
63
64			int is_reflection = (v.determinant() * w.determinant()) < 0.0;
65			if (is_reflection)
66			s(2) = -s(2);
67
68			RealType rmsd_sq = (E0 - 2.0 * s.sum() )/ (RealType)n_vec;
69			rmsd_sq = max(rmsd_sq,0.0);
70			RealType rmsd = sqrt(rmsd_sq);
71			return rmsd;
72			}
73
74
75			RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov) {
76			}
77