Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "RMSD.h"
23 : #include "PDB.h"
24 : #include "Log.h"
25 : #include "Exception.h"
26 : #include <cmath>
27 : #include <iostream>
28 : #include "Tools.h"
29 : using namespace std;
30 : namespace PLMD {
31 :
32 1923 : RMSD::RMSD() : alignmentMethod(SIMPLE),reference_center_is_calculated(false),reference_center_is_removed(false),positions_center_is_calculated(false),positions_center_is_removed(false) {}
33 :
34 : ///
35 : /// general method to set all the rmsd property at once by using a pdb where occupancy column sets the weights for the atoms involved in the
36 : /// alignment and beta sets the weight that are used for calculating the displacement.
37 : ///
38 949 : void RMSD::set(const PDB&pdb, const string & mytype, bool remove_center, bool normalize_weights ) {
39 :
40 949 : set(pdb.getOccupancy(),pdb.getBeta(),pdb.getPositions(),mytype,remove_center,normalize_weights);
41 :
42 949 : }
43 1877 : void RMSD::set(const std::vector<double> & align, const std::vector<double> & displace, const std::vector<Vector> & reference, const string & mytype, bool remove_center, bool normalize_weights ) {
44 :
45 1877 : setReference(reference); // this by default remove the com and assumes uniform weights
46 1877 : setAlign(align, normalize_weights, remove_center); // this recalculates the com with weights. If remove_center=false then it restore the center back
47 1877 : setDisplace(displace, normalize_weights); // this is does not affect any calculation of the weights
48 1877 : setType(mytype);
49 :
50 1877 : }
51 :
52 1877 : void RMSD::setType(const string & mytype) {
53 :
54 1877 : alignmentMethod=SIMPLE; // initialize with the simplest case: no rotation
55 1877 : if (mytype=="SIMPLE") {
56 0 : alignmentMethod=SIMPLE;
57 : }
58 1877 : else if (mytype=="OPTIMAL") {
59 1877 : alignmentMethod=OPTIMAL;
60 : }
61 0 : else if (mytype=="OPTIMAL-FAST") {
62 0 : alignmentMethod=OPTIMAL_FAST;
63 : }
64 0 : else plumed_merror("unknown RMSD type" + mytype);
65 :
66 1877 : }
67 :
68 923 : void RMSD::clear() {
69 923 : reference.clear();
70 923 : reference_center_is_calculated=false;
71 923 : reference_center_is_removed=false;
72 923 : align.clear();
73 923 : displace.clear();
74 923 : positions_center_is_calculated=false;
75 923 : positions_center_is_removed=false;
76 923 : }
77 :
78 4 : string RMSD::getMethod() {
79 4 : string mystring;
80 4 : switch(alignmentMethod) {
81 0 : case SIMPLE: mystring.assign("SIMPLE"); break;
82 4 : case OPTIMAL: mystring.assign("OPTIMAL"); break;
83 0 : case OPTIMAL_FAST: mystring.assign("OPTIMAL-FAST"); break;
84 : }
85 4 : return mystring;
86 : }
87 : ///
88 : /// this calculates the center of mass for the reference and removes it from the reference itself
89 : /// considering uniform weights for alignment
90 : ///
91 1877 : void RMSD::setReference(const vector<Vector> & reference) {
92 1877 : unsigned n=reference.size();
93 1877 : this->reference=reference;
94 1877 : plumed_massert(align.empty(),"you should first clear() an RMSD object, then set a new reference");
95 1877 : plumed_massert(displace.empty(),"you should first clear() an RMSD object, then set a new reference");
96 1877 : align.resize(n,1.0/n);
97 1877 : displace.resize(n,1.0/n);
98 1877 : for(unsigned i=0; i<n; i++) reference_center+=this->reference[i]*align[i];
99 1877 : for(unsigned i=0; i<n; i++) this->reference[i]-=reference_center;
100 1877 : reference_center_is_calculated=true;
101 1877 : reference_center_is_removed=true;
102 1877 : }
103 : ///
104 : /// the alignment weights are here normalized to 1 and the center of the reference is removed accordingly
105 : ///
106 1877 : void RMSD::setAlign(const vector<double> & align, bool normalize_weights, bool remove_center) {
107 1877 : unsigned n=reference.size();
108 1877 : plumed_massert(this->align.size()==align.size(),"mismatch in dimension of align/displace arrays");
109 1877 : this->align=align;
110 1877 : if(normalize_weights) {
111 1873 : double w=0.0;
112 1873 : for(unsigned i=0; i<n; i++) w+=this->align[i];
113 1873 : double inv=1.0/w;
114 1873 : for(unsigned i=0; i<n; i++) this->align[i]*=inv;
115 : }
116 : // recalculate the center anyway
117 : // just remove the center if that is asked
118 : // if the center was removed before, then add it and store the new one
119 1877 : if(reference_center_is_removed) {
120 1877 : plumed_massert(reference_center_is_calculated," seems that the reference center has been removed but not calculated and stored!");
121 1877 : addCenter(reference,reference_center);
122 : }
123 1877 : reference_center=calculateCenter(reference,this->align);
124 1877 : reference_center_is_calculated=true;
125 1877 : if(remove_center) {
126 1873 : removeCenter(reference,reference_center);
127 1873 : reference_center_is_removed=true;
128 : } else {
129 4 : reference_center_is_removed=false;
130 : }
131 1877 : }
132 : ///
133 : /// here the weigth for normalized weighths are normalized and set
134 : ///
135 1877 : void RMSD::setDisplace(const vector<double> & displace, bool normalize_weights) {
136 1877 : unsigned n=reference.size();
137 1877 : plumed_massert(this->displace.size()==displace.size(),"mismatch in dimension of align/displace arrays");
138 1877 : this->displace=displace;
139 1877 : double w=0.0;
140 1877 : for(unsigned i=0; i<n; i++) w+=this->displace[i];
141 1877 : double inv=1.0/w;
142 1877 : if(normalize_weights) {for(unsigned i=0; i<n; i++) this->displace[i]*=inv;}
143 1877 : }
144 : ///
145 : /// This is the main workhorse for rmsd that decides to use specific optimal alignment versions
146 : ///
147 285145 : double RMSD::calculate(const std::vector<Vector> & positions,std::vector<Vector> &derivatives, bool squared)const {
148 :
149 285145 : double ret=0.;
150 :
151 285145 : switch(alignmentMethod) {
152 : case SIMPLE : {
153 : // do a simple alignment without rotation
154 0 : std::vector<Vector> displacement( derivatives.size() );
155 0 : ret=simpleAlignment(align,displace,positions,reference,derivatives,displacement,squared);
156 0 : break;
157 : } case OPTIMAL_FAST : {
158 : // this is calling the fastest option:
159 0 : if(align==displace) ret=optimalAlignment<false,true>(align,displace,positions,reference,derivatives,squared);
160 0 : else ret=optimalAlignment<false,false>(align,displace,positions,reference,derivatives,squared);
161 0 : break;
162 :
163 : } case OPTIMAL : {
164 : // this is the fast routine but in the "safe" mode, which gives less numerical error:
165 285145 : if(align==displace) ret=optimalAlignment<true,true>(align,displace,positions,reference,derivatives,squared);
166 1 : else ret=optimalAlignment<true,false>(align,displace,positions,reference,derivatives,squared);
167 285145 : break;
168 : }
169 : }
170 :
171 285145 : return ret;
172 :
173 : }
174 :
175 :
176 : /// convenience method for calculating the standard derivatives and the derivative of the rmsd respect to the reference position
177 1 : double RMSD::calc_DDistDRef( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, std::vector<Vector>& DDistDRef, const bool squared ) {
178 1 : double ret=0.;
179 1 : switch(alignmentMethod) {
180 : case SIMPLE:
181 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
182 : break;
183 : case OPTIMAL_FAST:
184 0 : if(align==displace) ret=optimalAlignment_DDistDRef<false,true>(align,displace,positions,reference,derivatives,DDistDRef, squared);
185 0 : else ret=optimalAlignment_DDistDRef<false,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
186 0 : break;
187 : case OPTIMAL:
188 1 : if(align==displace) ret=optimalAlignment_DDistDRef<true,true>(align,displace,positions,reference,derivatives,DDistDRef,squared);
189 1 : else ret=optimalAlignment_DDistDRef<true,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
190 1 : break;
191 : }
192 1 : return ret;
193 :
194 : }
195 :
196 : /// convenience method for calculating the standard derivatives and the derivative of the rmsd respect to the reference position without the matrix contribution
197 : /// as required by SOMA
198 0 : double RMSD::calc_SOMA( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, std::vector<Vector>& DDistDRef, const bool squared ) {
199 0 : double ret=0.;
200 0 : switch(alignmentMethod) {
201 : case SIMPLE:
202 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
203 : break;
204 : case OPTIMAL_FAST:
205 0 : if(align==displace) ret=optimalAlignment_SOMA<false,true>(align,displace,positions,reference,derivatives,DDistDRef, squared);
206 0 : else ret=optimalAlignment_SOMA<false,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
207 0 : break;
208 : case OPTIMAL:
209 0 : if(align==displace) ret=optimalAlignment_SOMA<true,true>(align,displace,positions,reference,derivatives,DDistDRef,squared);
210 0 : else ret=optimalAlignment_SOMA<true,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
211 0 : break;
212 : }
213 0 : return ret;
214 :
215 : }
216 :
217 0 : double RMSD::calc_DDistDRef_Rot_DRotDPos( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, std::vector<Vector>& DDistDRef, Tensor & Rot, Matrix<std::vector<Vector> > &DRotDPos, const bool squared ) {
218 0 : double ret=0.;
219 0 : switch(alignmentMethod) {
220 : case SIMPLE:
221 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
222 : break;
223 : case OPTIMAL_FAST:
224 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos<false,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
225 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos<false,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
226 0 : break;
227 : case OPTIMAL:
228 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos<true,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
229 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos<true,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
230 0 : break;
231 : }
232 0 : return ret;
233 : }
234 :
235 0 : double RMSD::calc_DDistDRef_Rot_DRotDPos_DRotDRef( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, std::vector<Vector>& DDistDRef, Tensor & Rot, Matrix<std::vector<Vector> > &DRotDPos, Matrix<std::vector<Vector> > &DRotDRef, const bool squared ) {
236 0 : double ret=0.;
237 0 : switch(alignmentMethod) {
238 : case SIMPLE:
239 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
240 : break;
241 : case OPTIMAL_FAST:
242 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<false,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
243 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<false,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
244 0 : break;
245 : case OPTIMAL:
246 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<true,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
247 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<true,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
248 0 : break;
249 : }
250 0 : return ret;
251 : }
252 :
253 1670 : double RMSD::calc_PCAelements( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, Tensor & Rotation, Matrix<std::vector<Vector> > & DRotDPos,std::vector<Vector> & alignedpositions, std::vector<Vector> & centeredpositions, std::vector<Vector> ¢eredreference, const bool& squared ) const {
254 1670 : double ret=0.;
255 1670 : switch(alignmentMethod) {
256 : case SIMPLE:
257 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
258 : break;
259 : case OPTIMAL_FAST:
260 0 : if(align==displace) ret=optimalAlignment_PCA<false,true>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
261 0 : else ret=optimalAlignment_PCA<false,false>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
262 0 : break;
263 : case OPTIMAL:
264 1670 : if(align==displace) ret=optimalAlignment_PCA<true,true>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
265 0 : else ret=optimalAlignment_PCA<true,false>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
266 1670 : break;
267 : }
268 1670 : return ret;
269 : }
270 :
271 :
272 48 : double RMSD::calc_FitElements( const std::vector<Vector>& positions, Tensor & Rotation, Matrix<std::vector<Vector> > & DRotDPos, std::vector<Vector> & centeredpositions, Vector ¢er_positions, const bool& squared ) {
273 48 : double ret=0.;
274 48 : switch(alignmentMethod) {
275 : case SIMPLE:
276 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
277 : break;
278 : case OPTIMAL_FAST:
279 0 : if(align==displace)ret=optimalAlignment_Fit<false,true>(align,displace,positions,reference, Rotation,DRotDPos,centeredpositions,center_positions,squared);
280 0 : else ret=optimalAlignment_Fit<false,false>(align,displace,positions,reference, Rotation,DRotDPos,centeredpositions,center_positions,squared);
281 0 : break;
282 : case OPTIMAL:
283 48 : if(align==displace)ret=optimalAlignment_Fit<true,true>(align,displace,positions,reference,Rotation,DRotDPos,centeredpositions,center_positions,squared);
284 48 : else ret=optimalAlignment_Fit<true,false>(align,displace,positions,reference,Rotation,DRotDPos,centeredpositions,center_positions,squared);
285 48 : break;
286 : }
287 48 : return ret;
288 : }
289 :
290 :
291 :
292 :
293 :
294 :
295 148 : double RMSD::simpleAlignment(const std::vector<double> & align,
296 : const std::vector<double> & displace,
297 : const std::vector<Vector> & positions,
298 : const std::vector<Vector> & reference,
299 : std::vector<Vector> & derivatives,
300 : std::vector<Vector> & displacement,
301 : bool squared)const {
302 :
303 148 : double dist(0);
304 148 : unsigned n=reference.size();
305 :
306 148 : Vector apositions;
307 148 : Vector areference;
308 148 : Vector dpositions;
309 148 : Vector dreference;
310 :
311 2417 : for(unsigned i=0; i<n; i++) {
312 2269 : double aw=align[i];
313 2269 : double dw=displace[i];
314 2269 : apositions+=positions[i]*aw;
315 2269 : areference+=reference[i]*aw;
316 2269 : dpositions+=positions[i]*dw;
317 2269 : dreference+=reference[i]*dw;
318 : }
319 :
320 148 : Vector shift=((apositions-areference)-(dpositions-dreference));
321 2417 : for(unsigned i=0; i<n; i++) {
322 2269 : displacement[i]=(positions[i]-apositions)-(reference[i]-areference);
323 2269 : dist+=displace[i]*displacement[i].modulo2();
324 2269 : derivatives[i]=2*(displace[i]*displacement[i]+align[i]*shift);
325 : }
326 :
327 148 : if(!squared) {
328 : // sqrt
329 105 : dist=sqrt(dist);
330 : ///// sqrt on derivatives
331 105 : for(unsigned i=0; i<n; i++) {derivatives[i]*=(0.5/dist);}
332 : }
333 148 : return dist;
334 : }
335 :
336 : // this below enable the standard case for rmsd where the rmsd is calculated and the derivative of rmsd respect to positions is retrieved
337 : // additionally this assumes that the com of the reference is already subtracted.
338 : #define OLDRMSD
339 : #ifdef OLDRMSD
340 : // notice that in the current implementation the safe argument only makes sense for
341 : // align==displace
342 : template <bool safe,bool alEqDis>
343 476203 : double RMSD::optimalAlignment(const std::vector<double> & align,
344 : const std::vector<double> & displace,
345 : const std::vector<Vector> & positions,
346 : const std::vector<Vector> & reference,
347 : std::vector<Vector> & derivatives, bool squared)const {
348 476203 : double dist(0);
349 476203 : const unsigned n=reference.size();
350 : // This is the trace of positions*positions + reference*reference
351 476665 : double rr00(0);
352 476665 : double rr11(0);
353 : // This is positions*reference
354 476665 : Tensor rr01;
355 :
356 477309 : derivatives.resize(n);
357 :
358 477506 : Vector cpositions;
359 :
360 : // first expensive loop: compute centers
361 21363800 : for(unsigned iat=0; iat<n; iat++) {
362 20887338 : double w=align[iat];
363 20890490 : cpositions+=positions[iat]*w;
364 : }
365 :
366 : // second expensive loop: compute second moments wrt centers
367 21368884 : for(unsigned iat=0; iat<n; iat++) {
368 20891098 : double w=align[iat];
369 20891399 : rr00+=dotProduct(positions[iat]-cpositions,positions[iat]-cpositions)*w;
370 20888944 : rr11+=dotProduct(reference[iat],reference[iat])*w;
371 20890731 : rr01+=Tensor(positions[iat]-cpositions,reference[iat])*w;
372 : }
373 :
374 477786 : Matrix<double> m=Matrix<double>(4,4);
375 477797 : m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
376 477697 : m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
377 477716 : m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
378 477731 : m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
379 477716 : m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
380 477758 : m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
381 477731 : m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
382 477710 : m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
383 477747 : m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
384 477751 : m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
385 477753 : m[1][0] = m[0][1];
386 477757 : m[2][0] = m[0][2];
387 477767 : m[2][1] = m[1][2];
388 477752 : m[3][0] = m[0][3];
389 477762 : m[3][1] = m[1][3];
390 477752 : m[3][2] = m[2][3];
391 :
392 477734 : Tensor dm_drr01[4][4];
393 : if(!alEqDis) {
394 114 : dm_drr01[0][0] = 2.0*Tensor(-1.0, 0.0, 0.0, 0.0,-1.0, 0.0, 0.0, 0.0,-1.0);
395 114 : dm_drr01[1][1] = 2.0*Tensor(-1.0, 0.0, 0.0, 0.0,+1.0, 0.0, 0.0, 0.0,+1.0);
396 114 : dm_drr01[2][2] = 2.0*Tensor(+1.0, 0.0, 0.0, 0.0,-1.0, 0.0, 0.0, 0.0,+1.0);
397 114 : dm_drr01[3][3] = 2.0*Tensor(+1.0, 0.0, 0.0, 0.0,+1.0, 0.0, 0.0, 0.0,-1.0);
398 114 : dm_drr01[0][1] = 2.0*Tensor( 0.0, 0.0, 0.0, 0.0, 0.0,-1.0, 0.0,+1.0, 0.0);
399 114 : dm_drr01[0][2] = 2.0*Tensor( 0.0, 0.0,+1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0);
400 114 : dm_drr01[0][3] = 2.0*Tensor( 0.0,-1.0, 0.0, +1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
401 114 : dm_drr01[1][2] = 2.0*Tensor( 0.0,-1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
402 114 : dm_drr01[1][3] = 2.0*Tensor( 0.0, 0.0,-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0);
403 114 : dm_drr01[2][3] = 2.0*Tensor( 0.0, 0.0, 0.0, 0.0, 0.0,-1.0, 0.0,-1.0, 0.0);
404 114 : dm_drr01[1][0] = dm_drr01[0][1];
405 114 : dm_drr01[2][0] = dm_drr01[0][2];
406 114 : dm_drr01[2][1] = dm_drr01[1][2];
407 114 : dm_drr01[3][0] = dm_drr01[0][3];
408 114 : dm_drr01[3][1] = dm_drr01[1][3];
409 114 : dm_drr01[3][2] = dm_drr01[2][3];
410 : }
411 :
412 955557 : vector<double> eigenvals;
413 955087 : Matrix<double> eigenvecs;
414 477793 : int diagerror=diagMat(m, eigenvals, eigenvecs );
415 :
416 476878 : if (diagerror!=0) {
417 0 : string sdiagerror;
418 0 : Tools::convert(diagerror,sdiagerror);
419 0 : string msg="DIAGONALIZATION FAILED WITH ERROR CODE "+sdiagerror;
420 0 : plumed_merror(msg);
421 : }
422 :
423 476878 : dist=eigenvals[0]+rr00+rr11;
424 :
425 954751 : Matrix<double> ddist_dm(4,4);
426 :
427 477790 : Vector4d q(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
428 :
429 477427 : Tensor dq_drr01[4];
430 : if(!alEqDis) {
431 : double dq_dm[4][4][4];
432 7410 : for(unsigned i=0; i<4; i++) for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
433 7296 : double tmp=0.0;
434 : // perturbation theory for matrix m
435 7296 : for(unsigned l=1; l<4; l++) tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
436 7296 : dq_dm[i][j][k]=tmp;
437 : }
438 : // propagation to _drr01
439 570 : for(unsigned i=0; i<4; i++) {
440 456 : Tensor tmp;
441 7752 : for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
442 7296 : tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
443 : }
444 456 : dq_drr01[i]=tmp;
445 : }
446 : }
447 :
448 : // This is the rotation matrix that brings reference to positions
449 : // i.e. matmul(rotation,reference[iat])+shift is fitted to positions[iat]
450 :
451 477772 : Tensor rotation;
452 477664 : rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
453 477606 : rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
454 477562 : rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
455 477610 : rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
456 477679 : rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
457 477656 : rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
458 477647 : rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
459 477203 : rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
460 477433 : rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
461 :
462 :
463 477302 : Tensor drotation_drr01[3][3];
464 : if(!alEqDis) {
465 114 : drotation_drr01[0][0]=2*q[0]*dq_drr01[0]+2*q[1]*dq_drr01[1]-2*q[2]*dq_drr01[2]-2*q[3]*dq_drr01[3];
466 114 : drotation_drr01[1][1]=2*q[0]*dq_drr01[0]-2*q[1]*dq_drr01[1]+2*q[2]*dq_drr01[2]-2*q[3]*dq_drr01[3];
467 114 : drotation_drr01[2][2]=2*q[0]*dq_drr01[0]-2*q[1]*dq_drr01[1]-2*q[2]*dq_drr01[2]+2*q[3]*dq_drr01[3];
468 114 : drotation_drr01[0][1]=2*(+(q[0]*dq_drr01[3]+dq_drr01[0]*q[3])+(q[1]*dq_drr01[2]+dq_drr01[1]*q[2]));
469 114 : drotation_drr01[0][2]=2*(-(q[0]*dq_drr01[2]+dq_drr01[0]*q[2])+(q[1]*dq_drr01[3]+dq_drr01[1]*q[3]));
470 114 : drotation_drr01[1][2]=2*(+(q[0]*dq_drr01[1]+dq_drr01[0]*q[1])+(q[2]*dq_drr01[3]+dq_drr01[2]*q[3]));
471 114 : drotation_drr01[1][0]=2*(-(q[0]*dq_drr01[3]+dq_drr01[0]*q[3])+(q[1]*dq_drr01[2]+dq_drr01[1]*q[2]));
472 114 : drotation_drr01[2][0]=2*(+(q[0]*dq_drr01[2]+dq_drr01[0]*q[2])+(q[1]*dq_drr01[3]+dq_drr01[1]*q[3]));
473 114 : drotation_drr01[2][1]=2*(-(q[0]*dq_drr01[1]+dq_drr01[0]*q[1])+(q[2]*dq_drr01[3]+dq_drr01[2]*q[3]));
474 : }
475 :
476 477753 : double prefactor=2.0;
477 :
478 477639 : if(!squared && alEqDis) prefactor*=0.5/sqrt(dist);
479 :
480 : // if "safe", recompute dist here to a better accuracy
481 308634 : if(safe || !alEqDis) dist=0.0;
482 :
483 : // If safe is set to "false", MSD is taken from the eigenvalue of the M matrix
484 : // If safe is set to "true", MSD is recomputed from the rotational matrix
485 : // For some reason, this last approach leads to less numerical noise but adds an overhead
486 :
487 477753 : Tensor ddist_drotation;
488 477778 : Vector ddist_dcpositions;
489 :
490 : // third expensive loop: derivatives
491 21369757 : for(unsigned iat=0; iat<n; iat++) {
492 20891963 : Vector d(positions[iat]-cpositions - matmul(rotation,reference[iat]));
493 : if(alEqDis) {
494 : // there is no need for derivatives of rotation and shift here as it is by construction zero
495 : // (similar to Hellman-Feynman forces)
496 20889971 : derivatives[iat]= prefactor*align[iat]*d;
497 4333883 : if(safe) dist+=align[iat]*modulo2(d);
498 : } else {
499 : // the case for align != displace is different, sob:
500 2176 : dist+=displace[iat]*modulo2(d);
501 : // these are the derivatives assuming the roto-translation as frozen
502 2176 : derivatives[iat]=2*displace[iat]*d;
503 : // here I accumulate derivatives wrt rotation matrix ..
504 2176 : ddist_drotation+=-2*displace[iat]*extProduct(d,reference[iat]);
505 : // .. and cpositions
506 2176 : ddist_dcpositions+=-2*displace[iat]*d;
507 : }
508 : }
509 :
510 : if(!alEqDis) {
511 114 : Tensor ddist_drr01;
512 114 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
513 2290 : for(unsigned iat=0; iat<n; iat++) {
514 : // this is propagating to positions.
515 : // I am implicitly using the derivative of rr01 wrt positions here
516 2176 : derivatives[iat]+=matmul(ddist_drr01,reference[iat])*align[iat];
517 2176 : derivatives[iat]+=ddist_dcpositions*align[iat];
518 : }
519 : }
520 477794 : if(!squared) {
521 74160 : dist=sqrt(dist);
522 : if(!alEqDis) {
523 114 : double xx=0.5/dist;
524 114 : for(unsigned iat=0; iat<n; iat++) derivatives[iat]*=xx;
525 : }
526 : }
527 :
528 955589 : return dist;
529 : }
530 : #else
531 : /// note that this method is intended to be repeatedly invoked
532 : /// when the reference does already have the center subtracted
533 : /// but the position has not calculated center and not subtracted
534 : template <bool safe,bool alEqDis>
535 : double RMSD::optimalAlignment(const std::vector<double> & align,
536 : const std::vector<double> & displace,
537 : const std::vector<Vector> & positions,
538 : const std::vector<Vector> & reference,
539 : std::vector<Vector> & derivatives,
540 : bool squared) const {
541 : //std::cerr<<"setting up the core data \n";
542 : RMSDCoreData cd(align,displace,positions,reference);
543 :
544 : // transfer the settings for the center to let the CoreCalc deal with it
545 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
546 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
547 : else {cd.calcPositionsCenter();};
548 :
549 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
550 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
551 : else {cd.setReferenceCenter(reference_center);}
552 :
553 : // Perform the diagonalization and all the needed stuff
554 : cd.doCoreCalc(safe,alEqDis);
555 : // make the core calc distance
556 : double dist=cd.getDistance(squared);
557 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
558 : derivatives=cd.getDDistanceDPositions();
559 : return dist;
560 : }
561 : #endif
562 : template <bool safe,bool alEqDis>
563 1 : double RMSD::optimalAlignment_DDistDRef(const std::vector<double> & align,
564 : const std::vector<double> & displace,
565 : const std::vector<Vector> & positions,
566 : const std::vector<Vector> & reference,
567 : std::vector<Vector> & derivatives,
568 : std::vector<Vector> & ddistdref,
569 : bool squared) const {
570 : //initialize the data into the structure
571 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
572 1 : RMSDCoreData cd(align,displace,positions,reference);
573 : // transfer the settings for the center to let the CoreCalc deal with it
574 : // transfer the settings for the center to let the CoreCalc deal with it
575 1 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
576 1 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
577 1 : else {cd.calcPositionsCenter();};
578 :
579 1 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
580 1 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
581 1 : else {cd.setReferenceCenter(reference_center);}
582 :
583 : // Perform the diagonalization and all the needed stuff
584 1 : cd.doCoreCalc(safe,alEqDis);
585 : // make the core calc distance
586 1 : double dist=cd.getDistance(squared);
587 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
588 1 : derivatives=cd.getDDistanceDPositions();
589 1 : ddistdref=cd.getDDistanceDReference();
590 1 : return dist;
591 : }
592 :
593 : template <bool safe,bool alEqDis>
594 0 : double RMSD::optimalAlignment_SOMA(const std::vector<double> & align,
595 : const std::vector<double> & displace,
596 : const std::vector<Vector> & positions,
597 : const std::vector<Vector> & reference,
598 : std::vector<Vector> & derivatives,
599 : std::vector<Vector> & ddistdref,
600 : bool squared) const {
601 : //initialize the data into the structure
602 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
603 0 : RMSDCoreData cd(align,displace,positions,reference);
604 : // transfer the settings for the center to let the CoreCalc deal with it
605 : // transfer the settings for the center to let the CoreCalc deal with it
606 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
607 0 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
608 0 : else {cd.calcPositionsCenter();};
609 :
610 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
611 0 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
612 0 : else {cd.setReferenceCenter(reference_center);}
613 :
614 : // Perform the diagonalization and all the needed stuff
615 0 : cd.doCoreCalc(safe,alEqDis);
616 : // make the core calc distance
617 0 : double dist=cd.getDistance(squared);
618 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
619 0 : derivatives=cd.getDDistanceDPositions();
620 0 : ddistdref=cd.getDDistanceDReferenceSOMA();
621 0 : return dist;
622 : }
623 :
624 :
625 : template <bool safe,bool alEqDis>
626 0 : double RMSD::optimalAlignment_DDistDRef_Rot_DRotDPos(const std::vector<double> & align,
627 : const std::vector<double> & displace,
628 : const std::vector<Vector> & positions,
629 : const std::vector<Vector> & reference,
630 : std::vector<Vector> & derivatives,
631 : std::vector<Vector> & ddistdref,
632 : Tensor & Rotation,
633 : Matrix<std::vector<Vector> > &DRotDPos,
634 : bool squared) const {
635 : //initialize the data into the structure
636 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
637 0 : RMSDCoreData cd(align,displace,positions,reference);
638 : // transfer the settings for the center to let the CoreCalc deal with it
639 : // transfer the settings for the center to let the CoreCalc deal with it
640 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
641 0 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
642 0 : else {cd.calcPositionsCenter();};
643 :
644 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
645 0 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
646 0 : else {cd.setReferenceCenter(reference_center);}
647 :
648 : // Perform the diagonalization and all the needed stuff
649 0 : cd.doCoreCalc(safe,alEqDis);
650 : // make the core calc distance
651 0 : double dist=cd.getDistance(squared);
652 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
653 0 : derivatives=cd.getDDistanceDPositions();
654 0 : ddistdref=cd.getDDistanceDReference();
655 : // get the rotation matrix
656 0 : Rotation=cd.getRotationMatrixReferenceToPositions();
657 : // get its derivative
658 0 : DRotDPos=cd.getDRotationDPositions();
659 0 : return dist;
660 : }
661 :
662 : template <bool safe,bool alEqDis>
663 0 : double RMSD::optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef(const std::vector<double> & align,
664 : const std::vector<double> & displace,
665 : const std::vector<Vector> & positions,
666 : const std::vector<Vector> & reference,
667 : std::vector<Vector> & derivatives,
668 : std::vector<Vector> & ddistdref,
669 : Tensor & Rotation,
670 : Matrix<std::vector<Vector> > &DRotDPos,
671 : Matrix<std::vector<Vector> > &DRotDRef,
672 : bool squared) const {
673 : //initialize the data into the structure
674 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
675 0 : RMSDCoreData cd(align,displace,positions,reference);
676 : // transfer the settings for the center to let the CoreCalc deal with it
677 : // transfer the settings for the center to let the CoreCalc deal with it
678 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
679 0 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
680 0 : else {cd.calcPositionsCenter();};
681 :
682 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
683 0 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
684 0 : else {cd.setReferenceCenter(reference_center);}
685 :
686 : // Perform the diagonalization and all the needed stuff
687 0 : cd.doCoreCalc(safe,alEqDis);
688 : // make the core calc distance
689 0 : double dist=cd.getDistance(squared);
690 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
691 0 : derivatives=cd.getDDistanceDPositions();
692 0 : ddistdref=cd.getDDistanceDReference();
693 : // get the rotation matrix
694 0 : Rotation=cd.getRotationMatrixReferenceToPositions();
695 : // get its derivative
696 0 : DRotDPos=cd.getDRotationDPositions();
697 0 : DRotDRef=cd.getDRotationDReference();
698 0 : return dist;
699 : }
700 :
701 : template <bool safe,bool alEqDis>
702 1670 : double RMSD::optimalAlignment_PCA(const std::vector<double> & align,
703 : const std::vector<double> & displace,
704 : const std::vector<Vector> & positions,
705 : const std::vector<Vector> & reference,
706 : std::vector<Vector> & alignedpositions,
707 : std::vector<Vector> & centeredpositions,
708 : std::vector<Vector> & centeredreference,
709 : Tensor & Rotation,
710 : std::vector<Vector> & DDistDPos,
711 : Matrix<std::vector<Vector> > & DRotDPos,
712 : bool squared) const {
713 : //initialize the data into the structure
714 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
715 1670 : RMSDCoreData cd(align,displace,positions,reference);
716 : // transfer the settings for the center to let the CoreCalc deal with it
717 1670 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
718 1670 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
719 1670 : else {cd.calcPositionsCenter();};
720 :
721 1670 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
722 1670 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
723 1670 : else {cd.setReferenceCenter(reference_center);}
724 :
725 : // Perform the diagonalization and all the needed stuff
726 1670 : cd.doCoreCalc(safe,alEqDis);
727 : // make the core calc distance
728 1670 : double dist=cd.getDistance(squared);
729 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
730 1670 : DDistDPos=cd.getDDistanceDPositions();
731 : // get the rotation matrix
732 1670 : Rotation=cd.getRotationMatrixPositionsToReference();
733 : // get its derivative
734 1670 : DRotDPos=cd.getDRotationDPositions(true); // this gives back the inverse
735 : // get aligned positions
736 1670 : alignedpositions=cd.getAlignedPositionsToReference();
737 : // get centered positions
738 1670 : centeredpositions=cd.getCenteredPositions();
739 : // get centered reference
740 1670 : centeredreference=cd.getCenteredReference();
741 1670 : return dist;
742 : }
743 :
744 :
745 : template <bool safe,bool alEqDis>
746 48 : double RMSD::optimalAlignment_Fit(const std::vector<double> & align,
747 : const std::vector<double> & displace,
748 : const std::vector<Vector> & positions,
749 : const std::vector<Vector> & reference,
750 : Tensor & Rotation,
751 : Matrix<std::vector<Vector> > & DRotDPos,
752 : std::vector<Vector> & centeredpositions,
753 : Vector & center_positions,
754 : bool squared) {
755 : //initialize the data into the structure
756 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
757 48 : RMSDCoreData cd(align,displace,positions,reference);
758 : // transfer the settings for the center to let the CoreCalc deal with it
759 48 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
760 48 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
761 48 : else {cd.calcPositionsCenter();};
762 :
763 48 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
764 48 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
765 48 : else {cd.setReferenceCenter(reference_center);}
766 :
767 : // Perform the diagonalization and all the needed stuff
768 48 : cd.doCoreCalc(safe,alEqDis);
769 : // make the core calc distance
770 48 : double dist=cd.getDistance(squared);
771 : // get the rotation matrix
772 48 : Rotation=cd.getRotationMatrixPositionsToReference();
773 : // get its derivative
774 48 : DRotDPos=cd.getDRotationDPositions(true); // this gives back the inverse
775 : // get centered positions
776 48 : centeredpositions=cd.getCenteredPositions();
777 : // get center
778 48 : center_positions=cd.getPositionsCenter();
779 48 : return dist;
780 : }
781 :
782 :
783 :
784 :
785 :
786 :
787 : /// This calculates the elements needed by the quaternion to calculate everything that is needed
788 : /// additional calls retrieve different components
789 : /// note that this considers that the centers of both reference and positions are already setted
790 : /// but automatically should properly account for non removed components: if not removed then it
791 : /// removes prior to calculation of the alignment
792 1719 : void RMSDCoreData::doCoreCalc(bool safe,bool alEqDis, bool only_rotation) {
793 :
794 1719 : retrieve_only_rotation=only_rotation;
795 1719 : const unsigned n=static_cast<unsigned int>(reference.size());
796 :
797 1719 : plumed_massert(creference_is_calculated,"the center of the reference frame must be already provided at this stage");
798 1719 : plumed_massert(cpositions_is_calculated,"the center of the positions frame must be already provided at this stage");
799 :
800 : // This is the trace of positions*positions + reference*reference
801 1719 : rr00=0.;
802 1719 : rr11=0.;
803 : // This is positions*reference
804 1719 : Tensor rr01;
805 : // center of mass managing: must subtract the center from the position or not?
806 1719 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
807 1719 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
808 : // second expensive loop: compute second moments wrt centers
809 32976 : for(unsigned iat=0; iat<n; iat++) {
810 31257 : double w=align[iat];
811 31257 : rr00+=dotProduct(positions[iat]-cp,positions[iat]-cp)*w;
812 31257 : rr11+=dotProduct(reference[iat]-cr,reference[iat]-cr)*w;
813 31257 : rr01+=Tensor(positions[iat]-cp,reference[iat]-cr)*w;
814 : }
815 :
816 : // the quaternion matrix: this is internal
817 1719 : Matrix<double> m=Matrix<double>(4,4);
818 :
819 1719 : m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
820 1719 : m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
821 1719 : m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
822 1719 : m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
823 1719 : m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
824 1719 : m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
825 1719 : m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
826 1719 : m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
827 1719 : m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
828 1719 : m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
829 1719 : m[1][0] = m[0][1];
830 1719 : m[2][0] = m[0][2];
831 1719 : m[2][1] = m[1][2];
832 1719 : m[3][0] = m[0][3];
833 1719 : m[3][1] = m[1][3];
834 1719 : m[3][2] = m[2][3];
835 :
836 :
837 1719 : Tensor dm_drr01[4][4];
838 1719 : if(!alEqDis or !retrieve_only_rotation) {
839 1719 : dm_drr01[0][0] = 2.0*Tensor(-1.0, 0.0, 0.0, 0.0,-1.0, 0.0, 0.0, 0.0,-1.0);
840 1719 : dm_drr01[1][1] = 2.0*Tensor(-1.0, 0.0, 0.0, 0.0,+1.0, 0.0, 0.0, 0.0,+1.0);
841 1719 : dm_drr01[2][2] = 2.0*Tensor(+1.0, 0.0, 0.0, 0.0,-1.0, 0.0, 0.0, 0.0,+1.0);
842 1719 : dm_drr01[3][3] = 2.0*Tensor(+1.0, 0.0, 0.0, 0.0,+1.0, 0.0, 0.0, 0.0,-1.0);
843 1719 : dm_drr01[0][1] = 2.0*Tensor( 0.0, 0.0, 0.0, 0.0, 0.0,-1.0, 0.0,+1.0, 0.0);
844 1719 : dm_drr01[0][2] = 2.0*Tensor( 0.0, 0.0,+1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0);
845 1719 : dm_drr01[0][3] = 2.0*Tensor( 0.0,-1.0, 0.0, +1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
846 1719 : dm_drr01[1][2] = 2.0*Tensor( 0.0,-1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
847 1719 : dm_drr01[1][3] = 2.0*Tensor( 0.0, 0.0,-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0);
848 1719 : dm_drr01[2][3] = 2.0*Tensor( 0.0, 0.0, 0.0, 0.0, 0.0,-1.0, 0.0,-1.0, 0.0);
849 1719 : dm_drr01[1][0] = dm_drr01[0][1];
850 1719 : dm_drr01[2][0] = dm_drr01[0][2];
851 1719 : dm_drr01[2][1] = dm_drr01[1][2];
852 1719 : dm_drr01[3][0] = dm_drr01[0][3];
853 1719 : dm_drr01[3][1] = dm_drr01[1][3];
854 1719 : dm_drr01[3][2] = dm_drr01[2][3];
855 : }
856 :
857 :
858 1719 : int diagerror=diagMat(m, eigenvals, eigenvecs );
859 :
860 1719 : if (diagerror!=0) {
861 0 : std::string sdiagerror;
862 0 : Tools::convert(diagerror,sdiagerror);
863 0 : std::string msg="DIAGONALIZATION FAILED WITH ERROR CODE "+sdiagerror;
864 0 : plumed_merror(msg);
865 : }
866 :
867 1719 : Vector4d q(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
868 : // cerr<<"EIGENVAL "<<eigenvals[0]<<" "<<eigenvals[1]<<" "<<eigenvals[2]<<" "<<eigenvals[3]<<"\n";
869 :
870 1719 : Tensor dq_drr01[4];
871 1719 : if(!alEqDis or !only_rotation) {
872 : double dq_dm[4][4][4];
873 111735 : for(unsigned i=0; i<4; i++) for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
874 110016 : double tmp=0.0;
875 : // perturbation theory for matrix m
876 110016 : for(unsigned l=1; l<4; l++) tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
877 110016 : dq_dm[i][j][k]=tmp;
878 : }
879 : // propagation to _drr01
880 8595 : for(unsigned i=0; i<4; i++) {
881 6876 : Tensor tmp;
882 116892 : for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
883 110016 : tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
884 : }
885 6876 : dq_drr01[i]=tmp;
886 : }
887 : }
888 :
889 : // This is the rotation matrix that brings reference to positions
890 : // i.e. matmul(rotation,reference[iat])+shift is fitted to positions[iat]
891 :
892 1719 : rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
893 1719 : rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
894 1719 : rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
895 1719 : rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
896 1719 : rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
897 1719 : rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
898 1719 : rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
899 1719 : rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
900 1719 : rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
901 :
902 :
903 1719 : if(!alEqDis or !only_rotation) {
904 1719 : drotation_drr01[0][0]=2*q[0]*dq_drr01[0]+2*q[1]*dq_drr01[1]-2*q[2]*dq_drr01[2]-2*q[3]*dq_drr01[3];
905 1719 : drotation_drr01[1][1]=2*q[0]*dq_drr01[0]-2*q[1]*dq_drr01[1]+2*q[2]*dq_drr01[2]-2*q[3]*dq_drr01[3];
906 1719 : drotation_drr01[2][2]=2*q[0]*dq_drr01[0]-2*q[1]*dq_drr01[1]-2*q[2]*dq_drr01[2]+2*q[3]*dq_drr01[3];
907 1719 : drotation_drr01[0][1]=2*(+(q[0]*dq_drr01[3]+dq_drr01[0]*q[3])+(q[1]*dq_drr01[2]+dq_drr01[1]*q[2]));
908 1719 : drotation_drr01[0][2]=2*(-(q[0]*dq_drr01[2]+dq_drr01[0]*q[2])+(q[1]*dq_drr01[3]+dq_drr01[1]*q[3]));
909 1719 : drotation_drr01[1][2]=2*(+(q[0]*dq_drr01[1]+dq_drr01[0]*q[1])+(q[2]*dq_drr01[3]+dq_drr01[2]*q[3]));
910 1719 : drotation_drr01[1][0]=2*(-(q[0]*dq_drr01[3]+dq_drr01[0]*q[3])+(q[1]*dq_drr01[2]+dq_drr01[1]*q[2]));
911 1719 : drotation_drr01[2][0]=2*(+(q[0]*dq_drr01[2]+dq_drr01[0]*q[2])+(q[1]*dq_drr01[3]+dq_drr01[1]*q[3]));
912 1719 : drotation_drr01[2][1]=2*(-(q[0]*dq_drr01[1]+dq_drr01[0]*q[1])+(q[2]*dq_drr01[3]+dq_drr01[2]*q[3]));
913 : }
914 :
915 1719 : d.resize(n);
916 :
917 : // calculate rotation matrix derivatives and components distances needed for components only when align!=displacement
918 1719 : if(!alEqDis)ddist_drotation.zero();
919 32976 : for(unsigned iat=0; iat<n; iat++) {
920 : // components differences: this is useful externally
921 31257 : d[iat]=positions[iat]-cp - matmul(rotation,reference[iat]-cr);
922 : //cerr<<"D "<<iat<<" "<<d[iat][0]<<" "<<d[iat][1]<<" "<<d[iat][2]<<"\n";
923 :
924 : // ddist_drotation if needed
925 31257 : if(!alEqDis or !only_rotation) ddist_drotation+=-2*displace[iat]*extProduct(d[iat],reference[iat]-cr);
926 : }
927 :
928 1719 : if(!alEqDis or !only_rotation) {
929 1719 : ddist_drr01.zero();
930 1719 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
931 : }
932 : // transfer this bools to the cd so that this settings will be reflected in the other calls
933 1719 : this->alEqDis=alEqDis;
934 1719 : this->safe=safe;
935 1719 : isInitialized=true;
936 :
937 1719 : }
938 : /// just retrieve the distance already calculated
939 1719 : double RMSDCoreData::getDistance( bool squared) {
940 :
941 1719 : if(!isInitialized)plumed_merror("getDistance cannot calculate the distance without being initialized first by doCoreCalc ");
942 1719 : dist=eigenvals[0]+rr00+rr11;
943 :
944 1719 : if(safe || !alEqDis) dist=0.0;
945 1719 : const unsigned n=static_cast<unsigned int>(reference.size());
946 32976 : for(unsigned iat=0; iat<n; iat++) {
947 31257 : if(alEqDis) {
948 30162 : if(safe) dist+=align[iat]*modulo2(d[iat]);
949 : } else {
950 1095 : dist+=displace[iat]*modulo2(d[iat]);
951 : }
952 : }
953 1719 : if(!squared) {
954 49 : dist=sqrt(dist);
955 49 : distanceIsMSD=false;
956 : } else {
957 1670 : distanceIsMSD=true;
958 : }
959 1719 : hasDistance=true;
960 1719 : return dist;
961 : }
962 :
963 1671 : std::vector<Vector> RMSDCoreData::getDDistanceDPositions() {
964 1671 : std::vector<Vector> derivatives;
965 1671 : const unsigned n=static_cast<unsigned int>(reference.size());
966 1671 : Vector ddist_dcpositions;
967 1671 : derivatives.resize(n);
968 1671 : double prefactor=1.0;
969 1671 : if(!distanceIsMSD) prefactor*=0.5/dist;
970 1671 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
971 1671 : if(!hasDistance)plumed_merror("getDPositionsDerivatives needs to calculate the distance via getDistance first !");
972 1671 : if(!isInitialized)plumed_merror("getDPositionsDerivatives needs to initialize the coreData first!");
973 1671 : Vector csum;
974 1671 : Vector tmp1,tmp2;
975 32688 : for(unsigned iat=0; iat<n; iat++) {
976 31017 : if(alEqDis) {
977 : // there is no need for derivatives of rotation and shift here as it is by construction zero
978 : // (similar to Hellman-Feynman forces)
979 30162 : derivatives[iat]= 2*prefactor*align[iat]*d[iat];
980 : } else {
981 : // these are the derivatives assuming the roto-translation as frozen
982 855 : tmp1=2*displace[iat]*d[iat];
983 855 : derivatives[iat]=tmp1;
984 : // derivative of cpositions
985 855 : ddist_dcpositions+=-tmp1;
986 : // these needed for com corrections
987 855 : tmp2=matmul(ddist_drr01,reference[iat]-creference)*align[iat];
988 855 : derivatives[iat]+=tmp2;
989 855 : csum+=tmp2;
990 : }
991 : }
992 :
993 1671 : if(!alEqDis) for(unsigned iat=0; iat<n; iat++) {derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcpositions-csum)*align[iat]); }
994 :
995 1671 : return derivatives;
996 : }
997 :
998 1 : std::vector<Vector> RMSDCoreData::getDDistanceDReference() {
999 1 : std::vector<Vector> derivatives;
1000 1 : const unsigned n=static_cast<unsigned int>(reference.size());
1001 1 : Vector ddist_dcreference;
1002 1 : derivatives.resize(n);
1003 1 : double prefactor=1.0;
1004 1 : if(!distanceIsMSD) prefactor*=0.5/dist;
1005 1 : Vector csum,tmp1,tmp2;
1006 :
1007 1 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1008 1 : if(!hasDistance)plumed_merror("getDDistanceDReference needs to calculate the distance via getDistance first !");
1009 1 : if(!isInitialized)plumed_merror("getDDistanceDReference to initialize the coreData first!");
1010 : // get the transpose rotation
1011 1 : Tensor t_rotation=rotation.transpose();
1012 1 : Tensor t_ddist_drr01=ddist_drr01.transpose();
1013 :
1014 : // third expensive loop: derivatives
1015 856 : for(unsigned iat=0; iat<n; iat++) {
1016 855 : if(alEqDis) {
1017 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1018 : // (similar to Hellman-Feynman forces)
1019 : //TODO: check this derivative down here
1020 0 : derivatives[iat]= -2*prefactor*align[iat]*matmul(t_rotation,d[iat]);
1021 : } else {
1022 : // these are the derivatives assuming the roto-translation as frozen
1023 855 : tmp1=2*displace[iat]*matmul(t_rotation,d[iat]);
1024 855 : derivatives[iat]= -tmp1;
1025 : // derivative of cpositions
1026 855 : ddist_dcreference+=tmp1;
1027 : // these below are needed for com correction
1028 855 : tmp2=matmul(t_ddist_drr01,positions[iat]-cpositions)*align[iat];
1029 855 : derivatives[iat]+=tmp2;
1030 855 : csum+=tmp2;
1031 : }
1032 : }
1033 :
1034 1 : if(!alEqDis) for(unsigned iat=0; iat<n; iat++) {derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcreference-csum)*align[iat]);}
1035 :
1036 1 : return derivatives;
1037 : }
1038 :
1039 : /// this version does not calculate the derivative of rotation matrix as needed for SOMA
1040 0 : std::vector<Vector> RMSDCoreData::getDDistanceDReferenceSOMA() {
1041 0 : std::vector<Vector> derivatives;
1042 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1043 0 : Vector ddist_dcreference;
1044 0 : derivatives.resize(n);
1045 0 : double prefactor=1.0;
1046 0 : if(!distanceIsMSD) prefactor*=0.5/dist;
1047 0 : Vector csum,tmp1,tmp2;
1048 :
1049 0 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1050 0 : if(!hasDistance)plumed_merror("getDDistanceDReference needs to calculate the distance via getDistance first !");
1051 0 : if(!isInitialized)plumed_merror("getDDistanceDReference to initialize the coreData first!");
1052 : // get the transpose rotation
1053 0 : Tensor t_rotation=rotation.transpose();
1054 :
1055 : // third expensive loop: derivatives
1056 0 : for(unsigned iat=0; iat<n; iat++) {
1057 0 : if(alEqDis) {
1058 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1059 : // (similar to Hellman-Feynman forces)
1060 : //TODO: check this derivative down here
1061 0 : derivatives[iat]= -2*prefactor*align[iat]*matmul(t_rotation,d[iat]);
1062 : } else {
1063 : // these are the derivatives assuming the roto-translation as frozen
1064 0 : tmp1=2*displace[iat]*matmul(t_rotation,d[iat]);
1065 0 : derivatives[iat]= -tmp1;
1066 : // derivative of cpositions
1067 0 : ddist_dcreference+=tmp1;
1068 : }
1069 : }
1070 :
1071 0 : if(!alEqDis) for(unsigned iat=0; iat<n; iat++)derivatives[iat]=prefactor*(derivatives[iat]+ddist_dcreference*align[iat]);
1072 :
1073 0 : return derivatives;
1074 : }
1075 :
1076 :
1077 :
1078 : /*
1079 : This below is the derivative of the rotation matrix that aligns the reference onto the positions
1080 : respect to positions
1081 : note that the this transformation overlap the reference onto position
1082 : if inverseTransform=true then aligns the positions onto reference
1083 : */
1084 1718 : Matrix<std::vector<Vector> > RMSDCoreData::getDRotationDPositions( bool inverseTransform ) {
1085 1718 : const unsigned n=static_cast<unsigned int>(reference.size());
1086 1718 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1087 1718 : if(!isInitialized)plumed_merror("getDRotationDPosition to initialize the coreData first!");
1088 1718 : Matrix<std::vector<Vector> > DRotDPos=Matrix<std::vector<Vector> >(3,3);
1089 : // remember drotation_drr01 is Tensor drotation_drr01[3][3]
1090 : // (3x3 rot) (3x3 components of rr01)
1091 3436 : std::vector<Vector> v(n);
1092 1718 : Vector csum;
1093 : // these below could probably be calculated in the main routine
1094 1718 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1095 1718 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1096 1718 : for(unsigned iat=0; iat<n; iat++) csum+=(reference[iat]-cr)*align[iat];
1097 1718 : for(unsigned iat=0; iat<n; iat++) v[iat]=(reference[iat]-cr-csum)*align[iat];
1098 6872 : for(unsigned a=0; a<3; a++) {
1099 20616 : for(unsigned b=0; b<3; b++) {
1100 15462 : if(inverseTransform) {
1101 15462 : DRotDPos[b][a].resize(n);
1102 289080 : for(unsigned iat=0; iat<n; iat++) {
1103 273618 : DRotDPos[b][a][iat]=matmul(drotation_drr01[a][b],v[iat]);
1104 : }
1105 : } else {
1106 0 : DRotDPos[a][b].resize(n);
1107 0 : for(unsigned iat=0; iat<n; iat++) {
1108 0 : DRotDPos[a][b][iat]=matmul(drotation_drr01[a][b],v[iat]);
1109 : }
1110 : }
1111 : }
1112 : }
1113 3436 : return DRotDPos;
1114 : }
1115 :
1116 : /*
1117 : This below is the derivative of the rotation matrix that aligns the reference onto the positions
1118 : respect to reference
1119 : note that the this transformation overlap the reference onto position
1120 : if inverseTransform=true then aligns the positions onto reference
1121 : */
1122 0 : Matrix<std::vector<Vector> > RMSDCoreData::getDRotationDReference( bool inverseTransform ) {
1123 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1124 0 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1125 0 : if(!isInitialized)plumed_merror("getDRotationDPositions to initialize the coreData first!");
1126 0 : Matrix<std::vector<Vector> > DRotDRef=Matrix<std::vector<Vector> >(3,3);
1127 : // remember drotation_drr01 is Tensor drotation_drr01[3][3]
1128 : // (3x3 rot) (3x3 components of rr01)
1129 0 : std::vector<Vector> v(n);
1130 0 : Vector csum;
1131 : // these below could probably be calculated in the main routine
1132 0 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1133 0 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1134 0 : for(unsigned iat=0; iat<n; iat++) csum+=(positions[iat]-cp)*align[iat];
1135 0 : for(unsigned iat=0; iat<n; iat++) v[iat]=(positions[iat]-cp-csum)*align[iat];
1136 :
1137 0 : for(unsigned a=0; a<3; a++) {
1138 0 : for(unsigned b=0; b<3; b++) {
1139 0 : Tensor t_drotation_drr01=drotation_drr01[a][b].transpose();
1140 0 : if(inverseTransform) {
1141 0 : DRotDRef[b][a].resize(n);
1142 0 : for(unsigned iat=0; iat<n; iat++) {
1143 0 : DRotDRef[b][a][iat]=matmul(t_drotation_drr01,v[iat]);
1144 : }
1145 : } else {
1146 0 : DRotDRef[a][b].resize(n);
1147 0 : for(unsigned iat=0; iat<n; iat++) {
1148 0 : DRotDRef[a][b][iat]=matmul(t_drotation_drr01,v[iat]);
1149 : }
1150 : }
1151 : }
1152 : }
1153 0 : return DRotDRef;
1154 : }
1155 :
1156 :
1157 0 : std::vector<Vector> RMSDCoreData::getAlignedReferenceToPositions() {
1158 0 : std::vector<Vector> alignedref;
1159 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1160 0 : alignedref.resize(n);
1161 0 : if(!isInitialized)plumed_merror("getAlignedReferenceToPostions needs to initialize the coreData first!");
1162 : // avoid to calculate matrix element but use the sum of what you have
1163 0 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1164 0 : for(unsigned iat=0; iat<n; iat++)alignedref[iat]=-d[iat]+positions[iat]-cp;
1165 0 : return alignedref;
1166 : }
1167 1670 : std::vector<Vector> RMSDCoreData::getAlignedPositionsToReference() {
1168 1670 : std::vector<Vector> alignedpos;
1169 1670 : if(!isInitialized)plumed_merror("getAlignedPostionsToReference needs to initialize the coreData first!");
1170 1670 : const unsigned n=static_cast<unsigned int>(positions.size());
1171 1670 : alignedpos.resize(n);
1172 1670 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1173 : // avoid to calculate matrix element but use the sum of what you have
1174 1670 : for(unsigned iat=0; iat<n; iat++)alignedpos[iat]=matmul(rotation.transpose(),positions[iat]-cp);
1175 1670 : return alignedpos;
1176 : }
1177 :
1178 :
1179 1718 : std::vector<Vector> RMSDCoreData::getCenteredPositions() {
1180 1718 : std::vector<Vector> centeredpos;
1181 1718 : const unsigned n=static_cast<unsigned int>(reference.size());
1182 1718 : centeredpos.resize(n);
1183 1718 : if(!isInitialized)plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1184 : // avoid to calculate matrix element but use the sum of what you have
1185 1718 : for(unsigned iat=0; iat<n; iat++)centeredpos[iat]=positions[iat]-cpositions;
1186 1718 : return centeredpos;
1187 : }
1188 :
1189 1670 : std::vector<Vector> RMSDCoreData::getCenteredReference() {
1190 1670 : std::vector<Vector> centeredref;
1191 1670 : const unsigned n=static_cast<unsigned int>(reference.size());
1192 1670 : centeredref.resize(n);
1193 1670 : if(!isInitialized)plumed_merror("getCenteredReference needs to initialize the coreData first!");
1194 : // avoid to calculate matrix element but use the sum of what you have
1195 1670 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1196 1670 : for(unsigned iat=0; iat<n; iat++)centeredref[iat]=reference[iat]-cr;
1197 1670 : return centeredref;
1198 : }
1199 :
1200 :
1201 48 : Vector RMSDCoreData::getPositionsCenter() {
1202 48 : if(!isInitialized)plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1203 48 : return cpositions;
1204 : }
1205 :
1206 0 : Vector RMSDCoreData::getReferenceCenter() {
1207 0 : if(!isInitialized)plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1208 0 : return creference;
1209 : }
1210 :
1211 0 : Tensor RMSDCoreData::getRotationMatrixReferenceToPositions() {
1212 0 : if(!isInitialized)plumed_merror("getRotationMatrixReferenceToPositions needs to initialize the coreData first!");
1213 0 : return rotation;
1214 : }
1215 :
1216 1718 : Tensor RMSDCoreData::getRotationMatrixPositionsToReference() {
1217 1718 : if(!isInitialized)plumed_merror("getRotationMatrixReferenceToPositions needs to initialize the coreData first!");
1218 1718 : return rotation.transpose();
1219 : }
1220 :
1221 :
1222 : template double RMSD::optimalAlignment<true,true>(const std::vector<double> & align,
1223 : const std::vector<double> & displace,
1224 : const std::vector<Vector> & positions,
1225 : const std::vector<Vector> & reference,
1226 : std::vector<Vector> & derivatives, bool squared)const;
1227 : template double RMSD::optimalAlignment<true,false>(const std::vector<double> & align,
1228 : const std::vector<double> & displace,
1229 : const std::vector<Vector> & positions,
1230 : const std::vector<Vector> & reference,
1231 : std::vector<Vector> & derivatives, bool squared)const;
1232 : template double RMSD::optimalAlignment<false,true>(const std::vector<double> & align,
1233 : const std::vector<double> & displace,
1234 : const std::vector<Vector> & positions,
1235 : const std::vector<Vector> & reference,
1236 : std::vector<Vector> & derivatives, bool squared)const;
1237 : template double RMSD::optimalAlignment<false,false>(const std::vector<double> & align,
1238 : const std::vector<double> & displace,
1239 : const std::vector<Vector> & positions,
1240 : const std::vector<Vector> & reference,
1241 : std::vector<Vector> & derivatives, bool squared)const;
1242 :
1243 :
1244 :
1245 2523 : }
|
