ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/math/RMSD.cpp
Revision: 1349
Committed: Wed May 20 19:35:05 2009 UTC (16 years, 5 months ago) by cli2
File size: 3196 byte(s)
Log Message: 
Modifications to support RMSD calculations, removing classes we aren't using.

File Contents

# 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 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