src/math/RMSD.cpp

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

using namespace OpenMD;
using namespace JAMA;

RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov,
                              Vector3d mov_com,
                              Vector3d mov_to_ref) {
  
  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(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;

  // Pack everything into Dynamic arrays for the SVD:
  DynamicRectMatrix<RealType> Rtmp(3,3,0.0);
  DynamicRectMatrix<RealType> vtmp(3, 3);
  DynamicVector<RealType> stmp(3);
  DynamicRectMatrix<RealType> wtmp(3, 3);

  Rtmp.setSubMatrix(0, 0, R);
  SVD<RealType> svd(Rtmp);

  svd.getU(vtmp);
  svd.getSingularValues(stmp);
  svd.getV(wtmp);

  Mat3x3d v;
  Vector3d s;
  Mat3x3d w;

  vtmp.getSubMatrix(0,0,v);
  stmp.getSubVector(0,s);
  wtmp.getSubMatrix(0,0,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,
                                       Vector3d mov_com,
                                       Vector3d mov_to_ref) {
    
  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(0.0);
  RealType E0 = 0.0;
  
  for (n=0; n < n_vec; 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]);
  }

  // Pack everything into Dynamic arrays for the SVD:
  DynamicRectMatrix<RealType> Rtmp(3,3,0.0);
  DynamicRectMatrix<RealType> vtmp(3, 3);
  DynamicVector<RealType> stmp(3);
  DynamicRectMatrix<RealType> wtmp(3, 3);

  Rtmp.setSubMatrix(0, 0, R);
  SVD<RealType> svd(Rtmp);

  svd.getU(vtmp);
  svd.getSingularValues(stmp);
  svd.getV(wtmp);

  Mat3x3d v;
  Vector3d s;
  Mat3x3d w;

  vtmp.getSubMatrix(0,0,v);
  stmp.getSubVector(0,s);
  wtmp.getSubMatrix(0,0,w);
    
  int is_reflection = (v.determinant() * w.determinant()) < 0.0;
  if (is_reflection)
    s(2) = -s(2);

  RotMat3x3d U = v*w;
  return U;
}

Revision:	1442
Committed:	Mon May 10 17:28:26 2010 UTC (14 years, 11 months ago) by gezelter
File size:	3197 byte(s)
Log Message:	Adding property set to svn entries
#	User	Rev	Content
1	gezelter	1335	#include "math/RMSD.hpp"
2			#include "math/SVD.hpp"
3
4	gezelter	1390	using namespace OpenMD;
5	gezelter	1335	using namespace JAMA;
6
7			RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov,
8			Vector3d mov_com,
9	cli2	1349	Vector3d mov_to_ref) {
10
11	gezelter	1335	assert(mov.size() == ref_.size());
12			int n;
13			int n_vec = ref_.size();
14
15			/* calculate the centre of mass */
16			mov_com = V3Zero;
17
18			for (n=0; n < n_vec; n++) {
19			mov_com += mov[n];
20			}
21
22			mov_com /= (RealType)n_vec;
23
24			mov_to_ref = ref_com - mov_com;
25
26			/* shift mov to center of mass */
27
28			for (n=0; n < n_vec; n++) {
29			mov[n] -= mov_com;
30			}
31
32			/* initialize */
33	cli2	1349	Mat3x3d R(0.0);
34	gezelter	1335	RealType E0 = 0.0;
35
36			for (n=0; n < n_vec; n++) {
37
38			/*
39			* E0 = 1/2 * sum(over n): y(n)y(n) + x(n)x(n)
40			*/
41			E0 += dot(mov[n], mov[n]) + dot(ref_[n], ref_[n]);
42
43			/*
44			* correlation matrix R:
45			* R(i,j) = sum(over n): y(n,i) * x(n,j)
46			* where x(n) and y(n) are two vector sets
47			*/
48
49			R += outProduct(mov[n], ref_[n]);
50
51			}
52			E0 *= 0.5;
53
54	cli2	1349	// Pack everything into Dynamic arrays for the SVD:
55			DynamicRectMatrix<RealType> Rtmp(3,3,0.0);
56			DynamicRectMatrix<RealType> vtmp(3, 3);
57			DynamicVector<RealType> stmp(3);
58			DynamicRectMatrix<RealType> wtmp(3, 3);
59
60			Rtmp.setSubMatrix(0, 0, R);
61			SVD<RealType> svd(Rtmp);
62
63			svd.getU(vtmp);
64			svd.getSingularValues(stmp);
65			svd.getV(wtmp);
66
67			Mat3x3d v;
68	gezelter	1335	Vector3d s;
69			Mat3x3d w;
70
71	cli2	1349	vtmp.getSubMatrix(0,0,v);
72			stmp.getSubVector(0,s);
73			wtmp.getSubMatrix(0,0,w);
74	gezelter	1335
75			int is_reflection = (v.determinant() * w.determinant()) < 0.0;
76			if (is_reflection)
77			s(2) = -s(2);
78
79			RealType rmsd_sq = (E0 - 2.0 * s.sum() )/ (RealType)n_vec;
80			rmsd_sq = max(rmsd_sq,0.0);
81			RealType rmsd = sqrt(rmsd_sq);
82			return rmsd;
83			}
84
85	cli2	1349	RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov,
86			Vector3d mov_com,
87			Vector3d mov_to_ref) {
88
89			assert(mov.size() == ref_.size());
90			int n;
91			int n_vec = ref_.size();
92	gezelter	1335
93	cli2	1349	/* calculate the centre of mass */
94			mov_com = V3Zero;
95
96			for (n=0; n < n_vec; n++) {
97			mov_com += mov[n];
98			}
99
100			mov_com /= (RealType)n_vec;
101
102			mov_to_ref = ref_com - mov_com;
103
104			/* shift mov to center of mass */
105
106			for (n=0; n < n_vec; n++) {
107			mov[n] -= mov_com;
108			}
109
110			/* initialize */
111			Mat3x3d R(0.0);
112			RealType E0 = 0.0;
113
114			for (n=0; n < n_vec; n++) {
115
116			/*
117			* correlation matrix R:
118			* R(i,j) = sum(over n): y(n,i) * x(n,j)
119			* where x(n) and y(n) are two vector sets
120			*/
121
122			R += outProduct(mov[n], ref_[n]);
123			}
124
125			// Pack everything into Dynamic arrays for the SVD:
126			DynamicRectMatrix<RealType> Rtmp(3,3,0.0);
127			DynamicRectMatrix<RealType> vtmp(3, 3);
128			DynamicVector<RealType> stmp(3);
129			DynamicRectMatrix<RealType> wtmp(3, 3);
130
131			Rtmp.setSubMatrix(0, 0, R);
132			SVD<RealType> svd(Rtmp);
133
134			svd.getU(vtmp);
135			svd.getSingularValues(stmp);
136			svd.getV(wtmp);
137
138			Mat3x3d v;
139			Vector3d s;
140			Mat3x3d w;
141
142			vtmp.getSubMatrix(0,0,v);
143			stmp.getSubVector(0,s);
144			wtmp.getSubMatrix(0,0,w);
145
146			int is_reflection = (v.determinant() * w.determinant()) < 0.0;
147			if (is_reflection)
148			s(2) = -s(2);
149
150			RotMat3x3d U = v*w;
151			return U;
152	gezelter	1335	}
153