Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 94576 : 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 1057 : void RMSD::set(const PDB&pdb, const string & mytype, bool remove_center, bool normalize_weights ) {
39 :
40 1057 : set(pdb.getOccupancy(),pdb.getBeta(),pdb.getPositions(),mytype,remove_center,normalize_weights);
41 :
42 1057 : }
43 2056 : 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 2056 : setReference(reference); // this by default remove the com and assumes uniform weights
46 2056 : setAlign(align, normalize_weights, remove_center); // this recalculates the com with weights. If remove_center=false then it restore the center back
47 2056 : setDisplace(displace, normalize_weights); // this is does not affect any calculation of the weights
48 2056 : setType(mytype);
49 :
50 2056 : }
51 :
52 47248 : void RMSD::setType(const string & mytype) {
53 :
54 47248 : alignmentMethod=SIMPLE; // initialize with the simplest case: no rotation
55 47248 : if (mytype=="SIMPLE") {
56 0 : alignmentMethod=SIMPLE;
57 : }
58 47248 : else if (mytype=="OPTIMAL") {
59 47248 : 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 47248 : }
67 :
68 1010 : void RMSD::clear() {
69 : reference.clear();
70 1010 : reference_center.zero();
71 1010 : reference_center_is_calculated=false;
72 1010 : reference_center_is_removed=false;
73 : align.clear();
74 : displace.clear();
75 1010 : positions_center.zero();
76 1010 : positions_center_is_calculated=false;
77 1010 : positions_center_is_removed=false;
78 1010 : }
79 :
80 4 : string RMSD::getMethod() {
81 : string mystring;
82 4 : switch(alignmentMethod) {
83 0 : case SIMPLE: mystring.assign("SIMPLE"); break;
84 4 : case OPTIMAL: mystring.assign("OPTIMAL"); break;
85 0 : case OPTIMAL_FAST: mystring.assign("OPTIMAL-FAST"); break;
86 : }
87 4 : return mystring;
88 : }
89 : ///
90 : /// this calculates the center of mass for the reference and removes it from the reference itself
91 : /// considering uniform weights for alignment
92 : ///
93 47264 : void RMSD::setReference(const vector<Vector> & reference) {
94 47264 : unsigned n=reference.size();
95 47264 : this->reference=reference;
96 47264 : plumed_massert(align.empty(),"you should first clear() an RMSD object, then set a new reference");
97 47264 : plumed_massert(displace.empty(),"you should first clear() an RMSD object, then set a new reference");
98 47264 : align.resize(n,1.0/n);
99 47264 : displace.resize(n,1.0/n);
100 1913936 : for(unsigned i=0; i<n; i++) reference_center+=this->reference[i]*align[i];
101 : #pragma omp simd
102 1244448 : for(unsigned i=0; i<n; i++) this->reference[i]-=reference_center;
103 47264 : reference_center_is_calculated=true;
104 47264 : reference_center_is_removed=true;
105 47264 : }
106 45192 : std::vector<Vector> RMSD::getReference() {
107 45192 : return reference;
108 : }
109 : ///
110 : /// the alignment weights are here normalized to 1 and the center of the reference is removed accordingly
111 : ///
112 2056 : void RMSD::setAlign(const vector<double> & align, bool normalize_weights, bool remove_center) {
113 2056 : unsigned n=reference.size();
114 2056 : plumed_massert(this->align.size()==align.size(),"mismatch in dimension of align/displace arrays");
115 2056 : this->align=align;
116 2056 : if(normalize_weights) {
117 : double w=0.0;
118 4104 : #pragma omp simd reduction(+:w)
119 69000 : for(unsigned i=0; i<n; i++) w+=this->align[i];
120 2052 : plumed_massert(w>epsilon,"It looks like weights used for alignment are zero. Check your reference PDB file.");
121 2052 : double inv=1.0/w;
122 : #pragma omp simd
123 69000 : for(unsigned i=0; i<n; i++) this->align[i]*=inv;
124 : }
125 : // recalculate the center anyway
126 : // just remove the center if that is asked
127 : // if the center was removed before, then add it and store the new one
128 2056 : if(reference_center_is_removed) {
129 2056 : plumed_massert(reference_center_is_calculated," seems that the reference center has been removed but not calculated and stored!");
130 2056 : addCenter(reference,reference_center);
131 : }
132 2056 : reference_center=calculateCenter(reference,this->align);
133 2056 : reference_center_is_calculated=true;
134 2056 : if(remove_center) {
135 2052 : removeCenter(reference,reference_center);
136 2052 : reference_center_is_removed=true;
137 : } else {
138 4 : reference_center_is_removed=false;
139 : }
140 2056 : }
141 0 : std::vector<double> RMSD::getAlign() {
142 0 : return align;
143 : }
144 : ///
145 : /// here the weigth for normalized weighths are normalized and set
146 : ///
147 2056 : void RMSD::setDisplace(const vector<double> & displace, bool normalize_weights) {
148 2056 : unsigned n=reference.size();
149 2056 : plumed_massert(this->displace.size()==displace.size(),"mismatch in dimension of align/displace arrays");
150 2056 : this->displace=displace;
151 : double w=0.0;
152 4112 : #pragma omp simd reduction(+:w)
153 69040 : for(unsigned i=0; i<n; i++) w+=this->displace[i];
154 2056 : plumed_massert(w>epsilon,"It looks like weights used for displacement are zero. Check your reference PDB file.");
155 2056 : double inv=1.0/w;
156 2056 : if(normalize_weights) {
157 : #pragma omp simd
158 69000 : for(unsigned i=0; i<n; i++) this->displace[i]*=inv;
159 : }
160 2056 : }
161 0 : std::vector<double> RMSD::getDisplace() {
162 0 : return displace;
163 : }
164 : ///
165 : /// This is the main workhorse for rmsd that decides to use specific optimal alignment versions
166 : ///
167 330337 : double RMSD::calculate(const std::vector<Vector> & positions,std::vector<Vector> &derivatives, bool squared)const {
168 :
169 : double ret=0.;
170 :
171 330337 : switch(alignmentMethod) {
172 : case SIMPLE : {
173 : // do a simple alignment without rotation
174 0 : std::vector<Vector> displacement( derivatives.size() );
175 0 : ret=simpleAlignment(align,displace,positions,reference,derivatives,displacement,squared);
176 : break;
177 0 : } case OPTIMAL_FAST : {
178 : // this is calling the fastest option:
179 0 : if(align==displace) ret=optimalAlignment<false,true>(align,displace,positions,reference,derivatives,squared);
180 0 : else ret=optimalAlignment<false,false>(align,displace,positions,reference,derivatives,squared);
181 : break;
182 :
183 330337 : } case OPTIMAL : {
184 : // this is the fast routine but in the "safe" mode, which gives less numerical error:
185 330337 : if(align==displace) ret=optimalAlignment<true,true>(align,displace,positions,reference,derivatives,squared);
186 1 : else ret=optimalAlignment<true,false>(align,displace,positions,reference,derivatives,squared);
187 : break;
188 : }
189 : }
190 :
191 330337 : return ret;
192 :
193 : }
194 :
195 :
196 : /// convenience method for calculating the standard derivatives and the derivative of the rmsd respect to the reference position
197 1 : double RMSD::calc_DDistDRef( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, std::vector<Vector>& DDistDRef, const bool squared ) {
198 : double ret=0.;
199 1 : switch(alignmentMethod) {
200 0 : case SIMPLE:
201 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
202 : break;
203 0 : case OPTIMAL_FAST:
204 0 : if(align==displace) ret=optimalAlignment_DDistDRef<false,true>(align,displace,positions,reference,derivatives,DDistDRef, squared);
205 0 : else ret=optimalAlignment_DDistDRef<false,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
206 : break;
207 1 : case OPTIMAL:
208 1 : if(align==displace) ret=optimalAlignment_DDistDRef<true,true>(align,displace,positions,reference,derivatives,DDistDRef,squared);
209 1 : else ret=optimalAlignment_DDistDRef<true,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
210 : break;
211 : }
212 1 : return ret;
213 :
214 : }
215 :
216 : /// convenience method for calculating the standard derivatives and the derivative of the rmsd respect to the reference position without the matrix contribution
217 : /// as required by SOMA
218 0 : double RMSD::calc_SOMA( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, std::vector<Vector>& DDistDRef, const bool squared ) {
219 : double ret=0.;
220 0 : switch(alignmentMethod) {
221 0 : case SIMPLE:
222 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
223 : break;
224 0 : case OPTIMAL_FAST:
225 0 : if(align==displace) ret=optimalAlignment_SOMA<false,true>(align,displace,positions,reference,derivatives,DDistDRef, squared);
226 0 : else ret=optimalAlignment_SOMA<false,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
227 : break;
228 0 : case OPTIMAL:
229 0 : if(align==displace) ret=optimalAlignment_SOMA<true,true>(align,displace,positions,reference,derivatives,DDistDRef,squared);
230 0 : else ret=optimalAlignment_SOMA<true,false>(align,displace,positions,reference,derivatives,DDistDRef,squared);
231 : break;
232 : }
233 0 : return ret;
234 :
235 : }
236 :
237 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 ) {
238 : double ret=0.;
239 0 : switch(alignmentMethod) {
240 0 : case SIMPLE:
241 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
242 : break;
243 0 : case OPTIMAL_FAST:
244 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos<false,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
245 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos<false,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
246 : break;
247 0 : case OPTIMAL:
248 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos<true,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
249 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos<true,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
250 : break;
251 : }
252 0 : return ret;
253 : }
254 :
255 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 ) {
256 : double ret=0.;
257 0 : switch(alignmentMethod) {
258 0 : case SIMPLE:
259 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
260 : break;
261 0 : case OPTIMAL_FAST:
262 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<false,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
263 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<false,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
264 : break;
265 0 : case OPTIMAL:
266 0 : if(align==displace) ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<true,true>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
267 0 : else ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<true,false>(align,displace,positions,reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
268 : break;
269 : }
270 0 : return ret;
271 : }
272 :
273 1092 : double RMSD::calc_Rot_DRotDRr01( const std::vector<Vector>& positions, Tensor & Rotation, std::array<std::array<Tensor,3>,3> & DRotDRr01, const bool squared) {
274 : double ret=0.;
275 1092 : switch(alignmentMethod) {
276 0 : case SIMPLE:
277 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
278 : break;
279 0 : case OPTIMAL_FAST:
280 0 : if(align==displace) ret=optimalAlignment_Rot_DRotDRr01<false,true>(align,displace,positions,reference, Rotation, DRotDRr01, squared);
281 0 : else ret=optimalAlignment_Rot_DRotDRr01<false,false>(align,displace,positions,reference, Rotation, DRotDRr01, squared);
282 : break;
283 1092 : case OPTIMAL:
284 1092 : if(align==displace) ret=optimalAlignment_Rot_DRotDRr01<true,true>(align,displace,positions,reference, Rotation, DRotDRr01, squared);
285 0 : else ret=optimalAlignment_Rot_DRotDRr01<true,false>(align,displace,positions,reference, Rotation, DRotDRr01, squared);
286 : break;
287 : }
288 1092 : return ret;
289 : }
290 :
291 672 : double RMSD::calc_Rot( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, Tensor & Rotation, const bool squared) {
292 : double ret=0.;
293 672 : switch(alignmentMethod) {
294 0 : case SIMPLE:
295 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
296 : break;
297 0 : case OPTIMAL_FAST:
298 0 : if(align==displace) ret=optimalAlignment_Rot<false,true>(align,displace,positions,reference,derivatives, Rotation, squared);
299 0 : else ret=optimalAlignment_Rot<false,false>(align,displace,positions,reference,derivatives, Rotation, squared);
300 : break;
301 672 : case OPTIMAL:
302 672 : if(align==displace) ret=optimalAlignment_Rot<true,true>(align,displace,positions,reference,derivatives, Rotation, squared);
303 0 : else ret=optimalAlignment_Rot<true,false>(align,displace,positions,reference,derivatives, Rotation, squared);
304 : break;
305 : }
306 672 : return ret;
307 : }
308 :
309 45192 : double RMSD::calculateWithCloseStructure( const std::vector<Vector>& positions, std::vector<Vector> &derivatives, Tensor & rotationPosClose, Tensor & rotationRefClose, std::array<std::array<Tensor,3>,3> & drotationPosCloseDrr01, const bool squared) {
310 : double ret=0.;
311 45192 : switch(alignmentMethod) {
312 0 : case SIMPLE:
313 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
314 : break;
315 0 : case OPTIMAL_FAST:
316 0 : if(align==displace) ret=optimalAlignmentWithCloseStructure<false,true>(align,displace,positions,reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
317 0 : else ret=optimalAlignmentWithCloseStructure<false,false>(align,displace,positions,reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
318 : break;
319 45192 : case OPTIMAL:
320 45192 : if(align==displace) ret=optimalAlignmentWithCloseStructure<true,true>(align,displace,positions,reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
321 0 : else ret=optimalAlignmentWithCloseStructure<true,false>(align,displace,positions,reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
322 : break;
323 : }
324 45192 : return ret;
325 : }
326 :
327 4441 : 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 {
328 : double ret=0.;
329 4441 : switch(alignmentMethod) {
330 0 : case SIMPLE:
331 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
332 : break;
333 0 : case OPTIMAL_FAST:
334 0 : if(align==displace) ret=optimalAlignment_PCA<false,true>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
335 0 : else ret=optimalAlignment_PCA<false,false>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
336 : break;
337 4441 : case OPTIMAL:
338 4441 : if(align==displace) ret=optimalAlignment_PCA<true,true>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
339 0 : else ret=optimalAlignment_PCA<true,false>(align,displace,positions,reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
340 : break;
341 : }
342 4441 : return ret;
343 : }
344 :
345 :
346 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 ) {
347 : double ret=0.;
348 48 : switch(alignmentMethod) {
349 0 : case SIMPLE:
350 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
351 : break;
352 0 : case OPTIMAL_FAST:
353 0 : if(align==displace)ret=optimalAlignment_Fit<false,true>(align,displace,positions,reference, Rotation,DRotDPos,centeredpositions,center_positions,squared);
354 0 : else ret=optimalAlignment_Fit<false,false>(align,displace,positions,reference, Rotation,DRotDPos,centeredpositions,center_positions,squared);
355 : break;
356 48 : case OPTIMAL:
357 48 : if(align==displace)ret=optimalAlignment_Fit<true,true>(align,displace,positions,reference,Rotation,DRotDPos,centeredpositions,center_positions,squared);
358 48 : else ret=optimalAlignment_Fit<true,false>(align,displace,positions,reference,Rotation,DRotDPos,centeredpositions,center_positions,squared);
359 : break;
360 : }
361 48 : return ret;
362 : }
363 :
364 :
365 :
366 :
367 :
368 :
369 149 : double RMSD::simpleAlignment(const std::vector<double> & align,
370 : const std::vector<double> & displace,
371 : const std::vector<Vector> & positions,
372 : const std::vector<Vector> & reference,
373 : std::vector<Vector> & derivatives,
374 : std::vector<Vector> & displacement,
375 : bool squared)const {
376 :
377 : double dist(0);
378 149 : unsigned n=reference.size();
379 :
380 149 : Vector apositions;
381 149 : Vector areference;
382 149 : Vector dpositions;
383 149 : Vector dreference;
384 :
385 4693 : for(unsigned i=0; i<n; i++) {
386 4544 : double aw=align[i];
387 2272 : double dw=displace[i];
388 2272 : apositions+=positions[i]*aw;
389 2272 : areference+=reference[i]*aw;
390 2272 : dpositions+=positions[i]*dw;
391 2272 : dreference+=reference[i]*dw;
392 : }
393 :
394 149 : Vector shift=((apositions-areference)-(dpositions-dreference));
395 4693 : for(unsigned i=0; i<n; i++) {
396 9088 : displacement[i]=(positions[i]-apositions)-(reference[i]-areference);
397 4544 : dist+=displace[i]*displacement[i].modulo2();
398 6816 : derivatives[i]=2*(displace[i]*displacement[i]+align[i]*shift);
399 : }
400 :
401 149 : if(!squared) {
402 : // sqrt
403 106 : dist=sqrt(dist);
404 : ///// sqrt on derivatives
405 4192 : for(unsigned i=0; i<n; i++) {derivatives[i]*=(0.5/dist);}
406 : }
407 149 : return dist;
408 : }
409 :
410 : // this below enable the standard case for rmsd where the rmsd is calculated and the derivative of rmsd respect to positions is retrieved
411 : // additionally this assumes that the com of the reference is already subtracted.
412 : #define OLDRMSD
413 : #ifdef OLDRMSD
414 : // notice that in the current implementation the safe argument only makes sense for
415 : // align==displace
416 : template <bool safe,bool alEqDis>
417 546477 : double RMSD::optimalAlignment(const std::vector<double> & align,
418 : const std::vector<double> & displace,
419 : const std::vector<Vector> & positions,
420 : const std::vector<Vector> & reference,
421 : std::vector<Vector> & derivatives, bool squared)const {
422 : double dist(0);
423 546477 : const unsigned n=reference.size();
424 : // This is the trace of positions*positions + reference*reference
425 : double rr00(0);
426 : double rr11(0);
427 : // This is positions*reference
428 546477 : Tensor rr01;
429 :
430 546477 : derivatives.resize(n);
431 :
432 546477 : Vector cpositions;
433 :
434 : // first expensive loop: compute centers
435 44117285 : for(unsigned iat=0; iat<n; iat++) {
436 43570808 : double w=align[iat];
437 21785404 : cpositions+=positions[iat]*w;
438 : }
439 :
440 : // second expensive loop: compute second moments wrt centers
441 44117285 : for(unsigned iat=0; iat<n; iat++) {
442 43570808 : double w=align[iat];
443 43570808 : rr00+=dotProduct(positions[iat]-cpositions,positions[iat]-cpositions)*w;
444 21783228 : rr11+=dotProduct(reference[iat],reference[iat])*w;
445 21785404 : rr01+=Tensor(positions[iat]-cpositions,reference[iat])*w;
446 : }
447 :
448 : Matrix<double> m=Matrix<double>(4,4);
449 546477 : m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
450 546477 : m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
451 546477 : m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
452 546477 : m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
453 546477 : m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
454 546477 : m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
455 546477 : m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
456 546477 : m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
457 546477 : m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
458 546477 : m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
459 546477 : m[1][0] = m[0][1];
460 546477 : m[2][0] = m[0][2];
461 546477 : m[2][1] = m[1][2];
462 546477 : m[3][0] = m[0][3];
463 546477 : m[3][1] = m[1][3];
464 546477 : m[3][2] = m[2][3];
465 :
466 2732385 : Tensor dm_drr01[4][4];
467 : if(!alEqDis) {
468 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);
469 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);
470 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);
471 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);
472 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);
473 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);
474 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);
475 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);
476 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);
477 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);
478 114 : dm_drr01[1][0] = dm_drr01[0][1];
479 114 : dm_drr01[2][0] = dm_drr01[0][2];
480 114 : dm_drr01[2][1] = dm_drr01[1][2];
481 114 : dm_drr01[3][0] = dm_drr01[0][3];
482 114 : dm_drr01[3][1] = dm_drr01[1][3];
483 114 : dm_drr01[3][2] = dm_drr01[2][3];
484 : }
485 :
486 : vector<double> eigenvals;
487 : Matrix<double> eigenvecs;
488 546477 : int diagerror=diagMat(m, eigenvals, eigenvecs );
489 :
490 546477 : if (diagerror!=0) {
491 : string sdiagerror;
492 0 : Tools::convert(diagerror,sdiagerror);
493 0 : string msg="DIAGONALIZATION FAILED WITH ERROR CODE "+sdiagerror;
494 0 : plumed_merror(msg);
495 : }
496 :
497 546363 : dist=eigenvals[0]+rr00+rr11;
498 :
499 : Matrix<double> ddist_dm(4,4);
500 :
501 1092954 : Vector4d q(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
502 :
503 2732385 : Tensor dq_drr01[4];
504 : if(!alEqDis) {
505 : double dq_dm[4][4][4];
506 7866 : for(unsigned i=0; i<4; i++) for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
507 : double tmp=0.0;
508 : // perturbation theory for matrix m
509 72960 : for(unsigned l=1; l<4; l++) tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
510 7296 : dq_dm[i][j][k]=tmp;
511 : }
512 : // propagation to _drr01
513 1026 : for(unsigned i=0; i<4; i++) {
514 456 : Tensor tmp;
515 9576 : for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
516 7296 : tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
517 : }
518 456 : dq_drr01[i]=tmp;
519 : }
520 : }
521 :
522 : // This is the rotation matrix that brings reference to positions
523 : // i.e. matmul(rotation,reference[iat])+shift is fitted to positions[iat]
524 :
525 546477 : Tensor rotation;
526 546477 : rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
527 546477 : rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
528 546477 : rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
529 546477 : rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
530 546477 : rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
531 546477 : rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
532 546477 : rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
533 546477 : rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
534 546477 : rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
535 :
536 :
537 546477 : std::array<std::array<Tensor,3>,3> drotation_drr01;
538 : if(!alEqDis) {
539 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];
540 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];
541 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];
542 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]));
543 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]));
544 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]));
545 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]));
546 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]));
547 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]));
548 : }
549 :
550 : double prefactor=2.0;
551 :
552 546363 : if(!squared && alEqDis) prefactor*=0.5/sqrt(dist);
553 :
554 : // if "safe", recompute dist here to a better accuracy
555 : if(safe || !alEqDis) dist=0.0;
556 :
557 : // If safe is set to "false", MSD is taken from the eigenvalue of the M matrix
558 : // If safe is set to "true", MSD is recomputed from the rotational matrix
559 : // For some reason, this last approach leads to less numerical noise but adds an overhead
560 :
561 546477 : Tensor ddist_drotation;
562 546477 : Vector ddist_dcpositions;
563 :
564 : // third expensive loop: derivatives
565 44117285 : for(unsigned iat=0; iat<n; iat++) {
566 65356212 : Vector d(positions[iat]-cpositions - matmul(rotation,reference[iat]));
567 : if(alEqDis) {
568 : // there is no need for derivatives of rotation and shift here as it is by construction zero
569 : // (similar to Hellman-Feynman forces)
570 43566456 : derivatives[iat]= prefactor*align[iat]*d;
571 4921522 : if(safe) dist+=align[iat]*modulo2(d);
572 : } else {
573 : // the case for align != displace is different, sob:
574 2176 : dist+=displace[iat]*modulo2(d);
575 : // these are the derivatives assuming the roto-translation as frozen
576 4352 : derivatives[iat]=2*displace[iat]*d;
577 : // here I accumulate derivatives wrt rotation matrix ..
578 4352 : ddist_drotation+=-2*displace[iat]*extProduct(d,reference[iat]);
579 : // .. and cpositions
580 2176 : ddist_dcpositions+=-2*displace[iat]*d;
581 : }
582 : }
583 :
584 : if(!alEqDis) {
585 114 : Tensor ddist_drr01;
586 1482 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
587 4466 : for(unsigned iat=0; iat<n; iat++) {
588 : // this is propagating to positions.
589 : // I am implicitly using the derivative of rr01 wrt positions here
590 8704 : derivatives[iat]+=matmul(ddist_drr01,reference[iat])*align[iat];
591 4352 : derivatives[iat]+=ddist_dcpositions*align[iat];
592 : }
593 : }
594 546477 : if(!squared) {
595 74171 : dist=sqrt(dist);
596 : if(!alEqDis) {
597 114 : double xx=0.5/dist;
598 4466 : for(unsigned iat=0; iat<n; iat++) derivatives[iat]*=xx;
599 : }
600 : }
601 :
602 546477 : return dist;
603 : }
604 : #else
605 : /// note that this method is intended to be repeatedly invoked
606 : /// when the reference does already have the center subtracted
607 : /// but the position has not calculated center and not subtracted
608 : template <bool safe,bool alEqDis>
609 : double RMSD::optimalAlignment(const std::vector<double> & align,
610 : const std::vector<double> & displace,
611 : const std::vector<Vector> & positions,
612 : const std::vector<Vector> & reference,
613 : std::vector<Vector> & derivatives,
614 : bool squared) const {
615 : //std::cerr<<"setting up the core data \n";
616 : RMSDCoreData cd(align,displace,positions,reference);
617 :
618 : // transfer the settings for the center to let the CoreCalc deal with it
619 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
620 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
621 : else {cd.calcPositionsCenter();};
622 :
623 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
624 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
625 : else {cd.setReferenceCenter(reference_center);}
626 :
627 : // Perform the diagonalization and all the needed stuff
628 : cd.doCoreCalc(safe,alEqDis);
629 : // make the core calc distance
630 : double dist=cd.getDistance(squared);
631 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
632 : derivatives=cd.getDDistanceDPositions();
633 : return dist;
634 : }
635 : #endif
636 : template <bool safe,bool alEqDis>
637 1 : double RMSD::optimalAlignment_DDistDRef(const std::vector<double> & align,
638 : const std::vector<double> & displace,
639 : const std::vector<Vector> & positions,
640 : const std::vector<Vector> & reference,
641 : std::vector<Vector> & derivatives,
642 : std::vector<Vector> & ddistdref,
643 : bool squared) const {
644 : //initialize the data into the structure
645 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
646 2 : RMSDCoreData cd(align,displace,positions,reference);
647 : // transfer the settings for the center to let the CoreCalc deal with it
648 : // transfer the settings for the center to let the CoreCalc deal with it
649 1 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
650 1 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
651 1 : else {cd.calcPositionsCenter();};
652 :
653 1 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
654 1 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
655 1 : else {cd.setReferenceCenter(reference_center);}
656 :
657 : // Perform the diagonalization and all the needed stuff
658 1 : cd.doCoreCalc(safe,alEqDis);
659 : // make the core calc distance
660 1 : double dist=cd.getDistance(squared);
661 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
662 2 : derivatives=cd.getDDistanceDPositions();
663 2 : ddistdref=cd.getDDistanceDReference();
664 1 : return dist;
665 : }
666 :
667 : template <bool safe,bool alEqDis>
668 0 : double RMSD::optimalAlignment_SOMA(const std::vector<double> & align,
669 : const std::vector<double> & displace,
670 : const std::vector<Vector> & positions,
671 : const std::vector<Vector> & reference,
672 : std::vector<Vector> & derivatives,
673 : std::vector<Vector> & ddistdref,
674 : bool squared) const {
675 : //initialize the data into the structure
676 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
677 0 : RMSDCoreData cd(align,displace,positions,reference);
678 : // transfer the settings for the center to let the CoreCalc deal with it
679 : // transfer the settings for the center to let the CoreCalc deal with it
680 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
681 0 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
682 0 : else {cd.calcPositionsCenter();};
683 :
684 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
685 0 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
686 0 : else {cd.setReferenceCenter(reference_center);}
687 :
688 : // Perform the diagonalization and all the needed stuff
689 0 : cd.doCoreCalc(safe,alEqDis);
690 : // make the core calc distance
691 0 : double dist=cd.getDistance(squared);
692 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
693 0 : derivatives=cd.getDDistanceDPositions();
694 0 : ddistdref=cd.getDDistanceDReferenceSOMA();
695 0 : return dist;
696 : }
697 :
698 :
699 : template <bool safe,bool alEqDis>
700 0 : double RMSD::optimalAlignment_DDistDRef_Rot_DRotDPos(const std::vector<double> & align,
701 : const std::vector<double> & displace,
702 : const std::vector<Vector> & positions,
703 : const std::vector<Vector> & reference,
704 : std::vector<Vector> & derivatives,
705 : std::vector<Vector> & ddistdref,
706 : Tensor & Rotation,
707 : Matrix<std::vector<Vector> > &DRotDPos,
708 : bool squared) const {
709 : //initialize the data into the structure
710 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
711 0 : RMSDCoreData cd(align,displace,positions,reference);
712 : // transfer the settings for the center to let the CoreCalc deal with it
713 : // transfer the settings for the center to let the CoreCalc deal with it
714 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
715 0 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
716 0 : else {cd.calcPositionsCenter();};
717 :
718 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
719 0 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
720 0 : else {cd.setReferenceCenter(reference_center);}
721 :
722 : // Perform the diagonalization and all the needed stuff
723 0 : cd.doCoreCalc(safe,alEqDis);
724 : // make the core calc distance
725 0 : double dist=cd.getDistance(squared);
726 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
727 0 : derivatives=cd.getDDistanceDPositions();
728 0 : ddistdref=cd.getDDistanceDReference();
729 : // get the rotation matrix
730 0 : Rotation=cd.getRotationMatrixReferenceToPositions();
731 : // get its derivative
732 0 : DRotDPos=cd.getDRotationDPositions();
733 0 : return dist;
734 : }
735 :
736 : template <bool safe,bool alEqDis>
737 0 : double RMSD::optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef(const std::vector<double> & align,
738 : const std::vector<double> & displace,
739 : const std::vector<Vector> & positions,
740 : const std::vector<Vector> & reference,
741 : std::vector<Vector> & derivatives,
742 : std::vector<Vector> & ddistdref,
743 : Tensor & Rotation,
744 : Matrix<std::vector<Vector> > &DRotDPos,
745 : Matrix<std::vector<Vector> > &DRotDRef,
746 : bool squared) const {
747 : //initialize the data into the structure
748 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
749 0 : RMSDCoreData cd(align,displace,positions,reference);
750 : // transfer the settings for the center to let the CoreCalc deal with it
751 : // transfer the settings for the center to let the CoreCalc deal with it
752 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
753 0 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
754 0 : else {cd.calcPositionsCenter();};
755 :
756 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
757 0 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
758 0 : else {cd.setReferenceCenter(reference_center);}
759 :
760 : // Perform the diagonalization and all the needed stuff
761 0 : cd.doCoreCalc(safe,alEqDis);
762 : // make the core calc distance
763 0 : double dist=cd.getDistance(squared);
764 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
765 0 : derivatives=cd.getDDistanceDPositions();
766 0 : ddistdref=cd.getDDistanceDReference();
767 : // get the rotation matrix
768 0 : Rotation=cd.getRotationMatrixReferenceToPositions();
769 : // get its derivative
770 0 : DRotDPos=cd.getDRotationDPositions();
771 0 : DRotDRef=cd.getDRotationDReference();
772 0 : return dist;
773 : }
774 :
775 : template <bool safe,bool alEqDis>
776 1092 : double RMSD::optimalAlignment_Rot_DRotDRr01(const std::vector<double> & align,
777 : const std::vector<double> & displace,
778 : const std::vector<Vector> & positions,
779 : const std::vector<Vector> & reference,
780 : Tensor & Rotation,
781 : std::array<std::array<Tensor,3>,3> & DRotDRr01,
782 : bool squared) const {
783 : //initialize the data into the structure
784 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
785 2184 : RMSDCoreData cd(align,displace,positions,reference);
786 : // transfer the settings for the center to let the CoreCalc deal with it
787 1092 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
788 1092 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
789 1092 : else {cd.calcPositionsCenter();};
790 :
791 1092 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
792 1092 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
793 1092 : else {cd.setReferenceCenter(reference_center);}
794 :
795 : // Perform the diagonalization and all the needed stuff
796 1092 : cd.doCoreCalc(safe,alEqDis);
797 : // make the core calc distance
798 1092 : double dist=cd.getDistance(squared);
799 : // get the rotation matrix
800 1092 : Rotation=cd.getRotationMatrixReferenceToPositions();
801 : //get detivative w.r.t. rr01
802 1092 : DRotDRr01=cd.getDRotationDRr01();
803 1092 : return dist;
804 : }
805 :
806 : template <bool safe,bool alEqDis>
807 672 : double RMSD::optimalAlignment_Rot(const std::vector<double> & align,
808 : const std::vector<double> & displace,
809 : const std::vector<Vector> & positions,
810 : const std::vector<Vector> & reference,
811 : std::vector<Vector> & derivatives,
812 : Tensor & Rotation,
813 : bool squared) const {
814 : //initialize the data into the structure
815 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
816 1344 : RMSDCoreData cd(align,displace,positions,reference);
817 : // transfer the settings for the center to let the CoreCalc deal with it
818 672 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
819 672 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
820 672 : else {cd.calcPositionsCenter();};
821 :
822 672 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
823 672 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
824 672 : else {cd.setReferenceCenter(reference_center);}
825 :
826 : // Perform the diagonalization and all the needed stuff
827 672 : cd.doCoreCalc(safe,alEqDis);
828 : // make the core calc distance
829 672 : double dist=cd.getDistance(squared);
830 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
831 1344 : derivatives=cd.getDDistanceDPositions();
832 : // get the rotation matrix
833 672 : Rotation=cd.getRotationMatrixReferenceToPositions();
834 672 : return dist;
835 : }
836 :
837 : template <bool safe,bool alEqDis>
838 45192 : double RMSD::optimalAlignmentWithCloseStructure(const std::vector<double> & align,
839 : const std::vector<double> & displace,
840 : const std::vector<Vector> & positions,
841 : const std::vector<Vector> & reference,
842 : std::vector<Vector> & derivatives,
843 : Tensor & rotationPosClose,
844 : Tensor & rotationRefClose,
845 : std::array<std::array<Tensor,3>,3> & drotationPosCloseDrr01,
846 : bool squared) const {
847 : //initialize the data into the structure
848 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
849 90384 : RMSDCoreData cd(align,displace,positions,reference);
850 : // transfer the settings for the center to let the CoreCalc deal with it
851 45192 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
852 45192 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
853 45192 : else {cd.calcPositionsCenter();};
854 :
855 45192 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
856 45192 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
857 45192 : else {cd.setReferenceCenter(reference_center);}
858 :
859 : // instead of diagonalization, approximate with saved rotation matrix
860 45192 : cd.doCoreCalcWithCloseStructure(safe,alEqDis, rotationPosClose, rotationRefClose, drotationPosCloseDrr01);
861 : // make the core calc distance
862 45192 : double dist=cd.getDistance(squared);
863 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
864 90384 : derivatives=cd.getDDistanceDPositions();
865 45192 : return dist;
866 : }
867 :
868 :
869 : template <bool safe,bool alEqDis>
870 4441 : double RMSD::optimalAlignment_PCA(const std::vector<double> & align,
871 : const std::vector<double> & displace,
872 : const std::vector<Vector> & positions,
873 : const std::vector<Vector> & reference,
874 : std::vector<Vector> & alignedpositions,
875 : std::vector<Vector> & centeredpositions,
876 : std::vector<Vector> & centeredreference,
877 : Tensor & Rotation,
878 : std::vector<Vector> & DDistDPos,
879 : Matrix<std::vector<Vector> > & DRotDPos,
880 : bool squared) const {
881 : //initialize the data into the structure
882 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
883 8882 : RMSDCoreData cd(align,displace,positions,reference);
884 : // transfer the settings for the center to let the CoreCalc deal with it
885 4441 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
886 4441 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
887 4441 : else {cd.calcPositionsCenter();};
888 :
889 4441 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
890 4441 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
891 4441 : else {cd.setReferenceCenter(reference_center);}
892 :
893 : // Perform the diagonalization and all the needed stuff
894 4441 : cd.doCoreCalc(safe,alEqDis);
895 : // make the core calc distance
896 4441 : double dist=cd.getDistance(squared);
897 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
898 8882 : DDistDPos=cd.getDDistanceDPositions();
899 : // get the rotation matrix
900 4441 : Rotation=cd.getRotationMatrixPositionsToReference();
901 : // get its derivative
902 8882 : DRotDPos=cd.getDRotationDPositions(true); // this gives back the inverse
903 : // get aligned positions
904 8882 : alignedpositions=cd.getAlignedPositionsToReference();
905 : // get centered positions
906 8882 : centeredpositions=cd.getCenteredPositions();
907 : // get centered reference
908 8882 : centeredreference=cd.getCenteredReference();
909 4441 : return dist;
910 : }
911 :
912 :
913 : template <bool safe,bool alEqDis>
914 48 : double RMSD::optimalAlignment_Fit(const std::vector<double> & align,
915 : const std::vector<double> & displace,
916 : const std::vector<Vector> & positions,
917 : const std::vector<Vector> & reference,
918 : Tensor & Rotation,
919 : Matrix<std::vector<Vector> > & DRotDPos,
920 : std::vector<Vector> & centeredpositions,
921 : Vector & center_positions,
922 : bool squared) {
923 : //initialize the data into the structure
924 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
925 96 : RMSDCoreData cd(align,displace,positions,reference);
926 : // transfer the settings for the center to let the CoreCalc deal with it
927 48 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
928 48 : if(positions_center_is_calculated) {cd.setPositionsCenter(positions_center);}
929 48 : else {cd.calcPositionsCenter();};
930 :
931 48 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
932 48 : if(!reference_center_is_calculated) {cd.calcReferenceCenter();}
933 48 : else {cd.setReferenceCenter(reference_center);}
934 :
935 : // Perform the diagonalization and all the needed stuff
936 48 : cd.doCoreCalc(safe,alEqDis);
937 : // make the core calc distance
938 48 : double dist=cd.getDistance(squared);
939 : // get the rotation matrix
940 48 : Rotation=cd.getRotationMatrixPositionsToReference();
941 : // get its derivative
942 96 : DRotDPos=cd.getDRotationDPositions(true); // this gives back the inverse
943 : // get centered positions
944 96 : centeredpositions=cd.getCenteredPositions();
945 : // get center
946 48 : center_positions=cd.getPositionsCenter();
947 48 : return dist;
948 : }
949 :
950 :
951 :
952 :
953 :
954 :
955 : /// This calculates the elements needed by the quaternion to calculate everything that is needed
956 : /// additional calls retrieve different components
957 : /// note that this considers that the centers of both reference and positions are already setted
958 : /// but automatically should properly account for non removed components: if not removed then it
959 : /// removes prior to calculation of the alignment
960 6254 : void RMSDCoreData::doCoreCalc(bool safe,bool alEqDis, bool only_rotation) {
961 :
962 6254 : retrieve_only_rotation=only_rotation;
963 12508 : const unsigned n=static_cast<unsigned int>(reference.size());
964 :
965 6254 : plumed_massert(creference_is_calculated,"the center of the reference frame must be already provided at this stage");
966 6254 : plumed_massert(cpositions_is_calculated,"the center of the positions frame must be already provided at this stage");
967 :
968 : // This is the trace of positions*positions + reference*reference
969 6254 : rr00=0.;
970 6254 : rr11=0.;
971 : // This is positions*reference
972 6254 : Tensor rr01;
973 : // center of mass managing: must subtract the center from the position or not?
974 6254 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
975 6254 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
976 : // second expensive loop: compute second moments wrt centers
977 186678 : for(unsigned iat=0; iat<n; iat++) {
978 180424 : double w=align[iat];
979 270636 : rr00+=dotProduct(positions[iat]-cp,positions[iat]-cp)*w;
980 270636 : rr11+=dotProduct(reference[iat]-cr,reference[iat]-cr)*w;
981 270636 : rr01+=Tensor(positions[iat]-cp,reference[iat]-cr)*w;
982 : }
983 :
984 : // the quaternion matrix: this is internal
985 : Matrix<double> m=Matrix<double>(4,4);
986 :
987 6254 : m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
988 6254 : m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
989 6254 : m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
990 6254 : m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
991 6254 : m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
992 6254 : m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
993 6254 : m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
994 6254 : m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
995 6254 : m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
996 6254 : m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
997 6254 : m[1][0] = m[0][1];
998 6254 : m[2][0] = m[0][2];
999 6254 : m[2][1] = m[1][2];
1000 6254 : m[3][0] = m[0][3];
1001 6254 : m[3][1] = m[1][3];
1002 6254 : m[3][2] = m[2][3];
1003 :
1004 :
1005 31270 : Tensor dm_drr01[4][4];
1006 6254 : if(!alEqDis or !retrieve_only_rotation) {
1007 6254 : 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);
1008 6254 : 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);
1009 6254 : 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);
1010 6254 : 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);
1011 6254 : 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);
1012 6254 : 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);
1013 6254 : 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);
1014 6254 : 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);
1015 6254 : 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);
1016 6254 : 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);
1017 6254 : dm_drr01[1][0] = dm_drr01[0][1];
1018 6254 : dm_drr01[2][0] = dm_drr01[0][2];
1019 6254 : dm_drr01[2][1] = dm_drr01[1][2];
1020 6254 : dm_drr01[3][0] = dm_drr01[0][3];
1021 6254 : dm_drr01[3][1] = dm_drr01[1][3];
1022 6254 : dm_drr01[3][2] = dm_drr01[2][3];
1023 : }
1024 :
1025 :
1026 6254 : int diagerror=diagMat(m, eigenvals, eigenvecs );
1027 :
1028 6254 : if (diagerror!=0) {
1029 : std::string sdiagerror;
1030 0 : Tools::convert(diagerror,sdiagerror);
1031 0 : std::string msg="DIAGONALIZATION FAILED WITH ERROR CODE "+sdiagerror;
1032 0 : plumed_merror(msg);
1033 : }
1034 :
1035 12508 : Vector4d q(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
1036 : // cerr<<"EIGENVAL "<<eigenvals[0]<<" "<<eigenvals[1]<<" "<<eigenvals[2]<<" "<<eigenvals[3]<<"\n";
1037 :
1038 31270 : Tensor dq_drr01[4];
1039 6254 : if(!alEqDis or !only_rotation) {
1040 : double dq_dm[4][4][4];
1041 431526 : for(unsigned i=0; i<4; i++) for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
1042 : double tmp=0.0;
1043 : // perturbation theory for matrix m
1044 4002560 : for(unsigned l=1; l<4; l++) tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
1045 400256 : dq_dm[i][j][k]=tmp;
1046 : }
1047 : // propagation to _drr01
1048 56286 : for(unsigned i=0; i<4; i++) {
1049 25016 : Tensor tmp;
1050 525336 : for(unsigned j=0; j<4; j++) for(unsigned k=0; k<4; k++) {
1051 400256 : tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
1052 : }
1053 25016 : dq_drr01[i]=tmp;
1054 : }
1055 : }
1056 :
1057 : // This is the rotation matrix that brings reference to positions
1058 : // i.e. matmul(rotation,reference[iat])+shift is fitted to positions[iat]
1059 :
1060 6254 : rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
1061 6254 : rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
1062 6254 : rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
1063 6254 : rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
1064 6254 : rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
1065 6254 : rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
1066 6254 : rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
1067 6254 : rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
1068 6254 : rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
1069 :
1070 :
1071 6254 : if(!alEqDis or !only_rotation) {
1072 12508 : 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];
1073 12508 : 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];
1074 12508 : 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];
1075 12508 : 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]));
1076 12508 : 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]));
1077 12508 : 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]));
1078 12508 : 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]));
1079 12508 : 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]));
1080 12508 : 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]));
1081 : }
1082 :
1083 6254 : d.resize(n);
1084 :
1085 : // calculate rotation matrix derivatives and components distances needed for components only when align!=displacement
1086 6254 : if(!alEqDis)ddist_drotation.zero();
1087 : #pragma omp simd
1088 : for(unsigned iat=0; iat<n; iat++) {
1089 : // components differences: this is useful externally
1090 360848 : d[iat]=positions[iat]-cp - matmul(rotation,reference[iat]-cr);
1091 : //cerr<<"D "<<iat<<" "<<d[iat][0]<<" "<<d[iat][1]<<" "<<d[iat][2]<<"\n";
1092 : }
1093 : // ddist_drotation if needed
1094 6254 : if(!alEqDis or !only_rotation)
1095 186678 : for (unsigned iat=0; iat<n; iat++)
1096 360848 : ddist_drotation+=-2*displace[iat]*extProduct(d[iat],reference[iat]-cr);
1097 :
1098 6254 : if(!alEqDis or !only_rotation) {
1099 6254 : ddist_drr01.zero();
1100 81302 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
1101 : }
1102 : // transfer this bools to the cd so that this settings will be reflected in the other calls
1103 6254 : this->alEqDis=alEqDis;
1104 6254 : this->safe=safe;
1105 6254 : isInitialized=true;
1106 :
1107 6254 : }
1108 : /// just retrieve the distance already calculated
1109 51446 : double RMSDCoreData::getDistance( bool squared) {
1110 :
1111 51446 : if(!isInitialized)plumed_merror("getDistance cannot calculate the distance without being initialized first by doCoreCalc ");
1112 :
1113 : double localDist=0.0;
1114 102892 : const unsigned n=static_cast<unsigned int>(reference.size());
1115 51446 : if(safe || !alEqDis) localDist=0.0;
1116 : else
1117 0 : localDist=eigenvals[0]+rr00+rr11;
1118 102892 : #pragma omp simd reduction(+:localDist)
1119 : for(unsigned iat=0; iat<n; iat++) {
1120 677708 : if(alEqDis) {
1121 2029839 : if(safe) localDist+=align[iat]*modulo2(d[iat]);
1122 : } else {
1123 3285 : localDist+=displace[iat]*modulo2(d[iat]);
1124 : }
1125 : }
1126 51446 : if(!squared) {
1127 75 : dist=sqrt(localDist);
1128 75 : distanceIsMSD=false;
1129 : } else {
1130 51371 : dist=localDist;
1131 51371 : distanceIsMSD=true;
1132 : }
1133 51446 : hasDistance=true;
1134 51446 : return dist;
1135 : }
1136 :
1137 45192 : void RMSDCoreData::doCoreCalcWithCloseStructure(bool safe,bool alEqDis, Tensor & rotationPosClose, Tensor & rotationRefClose, std::array<std::array<Tensor,3>,3> & drotationPosCloseDrr01) {
1138 :
1139 90384 : unsigned natoms = reference.size();
1140 45192 : Tensor ddist_drxy;
1141 45192 : ddist_drr01.zero();
1142 45192 : d.resize(natoms);
1143 :
1144 : // center of mass managing: must subtract the center from the position or not?
1145 45192 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1146 45192 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1147 : //distance = \sum_{n=0}^{N} w_n(x_n-cpos-R_{XY} R_{AY} a_n)^2
1148 :
1149 45192 : Tensor rotation = matmul(rotationPosClose, rotationRefClose);
1150 :
1151 : #pragma omp simd
1152 : for (unsigned iat=0; iat<natoms; iat++) {
1153 2349984 : d[iat] = positions[iat] - cp - matmul(rotation, reference[iat]-cr);
1154 : }
1155 45192 : if (!alEqDis) {
1156 0 : for (unsigned iat=0; iat<natoms; iat++) {
1157 : //dist = \sum w_i(x_i - cpos - R_xy * R_ay * a_i)
1158 0 : ddist_drxy += -2*displace[iat]*extProduct(matmul(d[iat], rotationRefClose), reference[iat]-cr);
1159 : }
1160 : }
1161 :
1162 45192 : if (!alEqDis) {
1163 0 : for(unsigned i=0; i<3; i++)
1164 0 : for(unsigned j=0; j<3; j++)
1165 0 : ddist_drr01+=ddist_drxy[i][j]*drotationPosCloseDrr01[i][j];
1166 : }
1167 45192 : this->alEqDis=alEqDis;
1168 45192 : this->safe=safe;
1169 45192 : isInitialized=true;
1170 45192 : }
1171 :
1172 50306 : std::vector<Vector> RMSDCoreData::getDDistanceDPositions() {
1173 : std::vector<Vector> derivatives;
1174 100612 : const unsigned n=static_cast<unsigned int>(reference.size());
1175 50306 : Vector ddist_dcpositions;
1176 50306 : derivatives.resize(n);
1177 : double prefactor=1.0;
1178 50306 : if(!distanceIsMSD) prefactor*=0.5/dist;
1179 50306 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1180 50306 : if(!hasDistance)plumed_merror("getDPositionsDerivatives needs to calculate the distance via getDistance first !");
1181 50306 : if(!isInitialized)plumed_merror("getDPositionsDerivatives needs to initialize the coreData first!");
1182 50306 : Vector csum;
1183 1376850 : for(unsigned iat=0; iat<n; iat++) {
1184 663272 : if(alEqDis) {
1185 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1186 : // (similar to Hellman-Feynman forces)
1187 2649668 : derivatives[iat]= 2*prefactor*align[iat]*d[iat];
1188 : } else {
1189 : // these are the derivatives assuming the roto-translation as frozen
1190 2565 : Vector tmp1=2*displace[iat]*d[iat];
1191 855 : derivatives[iat]=tmp1;
1192 : // derivative of cpositions
1193 855 : ddist_dcpositions+=-tmp1;
1194 : // these needed for com corrections
1195 2565 : Vector tmp2=matmul(ddist_drr01,reference[iat]-creference)*align[iat];
1196 855 : derivatives[iat]+=tmp2;
1197 855 : csum+=tmp2;
1198 : }
1199 : }
1200 :
1201 50306 : if(!alEqDis)
1202 : #pragma omp simd
1203 3420 : for(unsigned iat=0; iat<n; iat++) {derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcpositions-csum)*align[iat]); }
1204 :
1205 50306 : return derivatives;
1206 : }
1207 :
1208 1 : std::vector<Vector> RMSDCoreData::getDDistanceDReference() {
1209 : std::vector<Vector> derivatives;
1210 2 : const unsigned n=static_cast<unsigned int>(reference.size());
1211 1 : Vector ddist_dcreference;
1212 1 : derivatives.resize(n);
1213 : double prefactor=1.0;
1214 1 : if(!distanceIsMSD) prefactor*=0.5/dist;
1215 1 : Vector csum;
1216 :
1217 1 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1218 1 : if(!hasDistance)plumed_merror("getDDistanceDReference needs to calculate the distance via getDistance first !");
1219 1 : if(!isInitialized)plumed_merror("getDDistanceDReference to initialize the coreData first!");
1220 : // get the transpose rotation
1221 1 : Tensor t_rotation=rotation.transpose();
1222 1 : Tensor t_ddist_drr01=ddist_drr01.transpose();
1223 :
1224 : // third expensive loop: derivatives
1225 1711 : for(unsigned iat=0; iat<n; iat++) {
1226 855 : if(alEqDis) {
1227 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1228 : // (similar to Hellman-Feynman forces)
1229 : //TODO: check this derivative down here
1230 0 : derivatives[iat]= -2*prefactor*align[iat]*matmul(t_rotation,d[iat]);
1231 : } else {
1232 : // these are the derivatives assuming the roto-translation as frozen
1233 2565 : Vector tmp1=2*displace[iat]*matmul(t_rotation,d[iat]);
1234 855 : derivatives[iat]= -tmp1;
1235 : // derivative of cpositions
1236 855 : ddist_dcreference+=tmp1;
1237 : // these below are needed for com correction
1238 2565 : Vector tmp2=matmul(t_ddist_drr01,positions[iat]-cpositions)*align[iat];
1239 855 : derivatives[iat]+=tmp2;
1240 855 : csum+=tmp2;
1241 : }
1242 : }
1243 :
1244 1 : if(!alEqDis)
1245 : #pragma omp simd
1246 3420 : for(unsigned iat=0; iat<n; iat++) {derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcreference-csum)*align[iat]);}
1247 :
1248 1 : return derivatives;
1249 : }
1250 :
1251 : /// this version does not calculate the derivative of rotation matrix as needed for SOMA
1252 0 : std::vector<Vector> RMSDCoreData::getDDistanceDReferenceSOMA() {
1253 : std::vector<Vector> derivatives;
1254 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1255 0 : Vector ddist_dcreference;
1256 0 : derivatives.resize(n);
1257 : double prefactor=1.0;
1258 0 : if(!distanceIsMSD) prefactor*=0.5/dist;
1259 0 : Vector csum,tmp1,tmp2;
1260 :
1261 0 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1262 0 : if(!hasDistance)plumed_merror("getDDistanceDReference needs to calculate the distance via getDistance first !");
1263 0 : if(!isInitialized)plumed_merror("getDDistanceDReference to initialize the coreData first!");
1264 : // get the transpose rotation
1265 0 : Tensor t_rotation=rotation.transpose();
1266 :
1267 : // third expensive loop: derivatives
1268 0 : for(unsigned iat=0; iat<n; iat++) {
1269 0 : if(alEqDis) {
1270 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1271 : // (similar to Hellman-Feynman forces)
1272 : //TODO: check this derivative down here
1273 0 : derivatives[iat]= -2*prefactor*align[iat]*matmul(t_rotation,d[iat]);
1274 : } else {
1275 : // these are the derivatives assuming the roto-translation as frozen
1276 0 : tmp1=2*displace[iat]*matmul(t_rotation,d[iat]);
1277 0 : derivatives[iat]= -tmp1;
1278 : // derivative of cpositions
1279 0 : ddist_dcreference+=tmp1;
1280 : }
1281 : }
1282 :
1283 0 : if(!alEqDis) for(unsigned iat=0; iat<n; iat++)derivatives[iat]=prefactor*(derivatives[iat]+ddist_dcreference*align[iat]);
1284 :
1285 0 : return derivatives;
1286 : }
1287 :
1288 :
1289 :
1290 : /*
1291 : This below is the derivative of the rotation matrix that aligns the reference onto the positions
1292 : respect to positions
1293 : note that the this transformation overlap the reference onto position
1294 : if inverseTransform=true then aligns the positions onto reference
1295 : */
1296 4489 : Matrix<std::vector<Vector> > RMSDCoreData::getDRotationDPositions( bool inverseTransform ) {
1297 8978 : const unsigned n=static_cast<unsigned int>(reference.size());
1298 4489 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1299 4489 : if(!isInitialized)plumed_merror("getDRotationDPosition to initialize the coreData first!");
1300 : Matrix<std::vector<Vector> > DRotDPos=Matrix<std::vector<Vector> >(3,3);
1301 : // remember drotation_drr01 is Tensor drotation_drr01[3][3]
1302 : // (3x3 rot) (3x3 components of rr01)
1303 4489 : std::vector<Vector> v(n);
1304 4489 : Vector csum;
1305 : // these below could probably be calculated in the main routine
1306 4489 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1307 4489 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1308 203764 : for(unsigned iat=0; iat<n; iat++) csum+=(reference[iat]-cr)*align[iat];
1309 270189 : for(unsigned iat=0; iat<n; iat++) v[iat]=(reference[iat]-cr-csum)*align[iat];
1310 31423 : for(unsigned a=0; a<3; a++) {
1311 94269 : for(unsigned b=0; b<3; b++) {
1312 40401 : if(inverseTransform) {
1313 40401 : DRotDPos[b][a].resize(n);
1314 1236051 : for(unsigned iat=0; iat<n; iat++) {
1315 2391300 : DRotDPos[b][a][iat]=matmul(drotation_drr01[a][b],v[iat]);
1316 : }
1317 : } else {
1318 0 : DRotDPos[a][b].resize(n);
1319 0 : for(unsigned iat=0; iat<n; iat++) {
1320 0 : DRotDPos[a][b][iat]=matmul(drotation_drr01[a][b],v[iat]);
1321 : }
1322 : }
1323 : }
1324 : }
1325 4489 : return DRotDPos;
1326 : }
1327 :
1328 : /*
1329 : This below is the derivative of the rotation matrix that aligns the reference onto the positions
1330 : respect to reference
1331 : note that the this transformation overlap the reference onto position
1332 : if inverseTransform=true then aligns the positions onto reference
1333 : */
1334 0 : Matrix<std::vector<Vector> > RMSDCoreData::getDRotationDReference( bool inverseTransform ) {
1335 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1336 0 : plumed_massert(!retrieve_only_rotation,"You used only_rotation=true in doCoreCalc therefore you cannot retrieve this information now");
1337 0 : if(!isInitialized)plumed_merror("getDRotationDPositions to initialize the coreData first!");
1338 : Matrix<std::vector<Vector> > DRotDRef=Matrix<std::vector<Vector> >(3,3);
1339 : // remember drotation_drr01 is Tensor drotation_drr01[3][3]
1340 : // (3x3 rot) (3x3 components of rr01)
1341 0 : std::vector<Vector> v(n);
1342 0 : Vector csum;
1343 : // these below could probably be calculated in the main routine
1344 0 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1345 0 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1346 0 : for(unsigned iat=0; iat<n; iat++) csum+=(positions[iat]-cp)*align[iat];
1347 0 : for(unsigned iat=0; iat<n; iat++) v[iat]=(positions[iat]-cp-csum)*align[iat];
1348 :
1349 0 : for(unsigned a=0; a<3; a++) {
1350 0 : for(unsigned b=0; b<3; b++) {
1351 0 : Tensor t_drotation_drr01=drotation_drr01[a][b].transpose();
1352 0 : if(inverseTransform) {
1353 0 : DRotDRef[b][a].resize(n);
1354 0 : for(unsigned iat=0; iat<n; iat++) {
1355 0 : DRotDRef[b][a][iat]=matmul(t_drotation_drr01,v[iat]);
1356 : }
1357 : } else {
1358 0 : DRotDRef[a][b].resize(n);
1359 0 : for(unsigned iat=0; iat<n; iat++) {
1360 0 : DRotDRef[a][b][iat]=matmul(t_drotation_drr01,v[iat]);
1361 : }
1362 : }
1363 : }
1364 : }
1365 0 : return DRotDRef;
1366 : }
1367 :
1368 :
1369 0 : std::vector<Vector> RMSDCoreData::getAlignedReferenceToPositions() {
1370 : std::vector<Vector> alignedref;
1371 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1372 0 : alignedref.resize(n);
1373 0 : if(!isInitialized)plumed_merror("getAlignedReferenceToPostions needs to initialize the coreData first!");
1374 : // avoid to calculate matrix element but use the sum of what you have
1375 0 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1376 0 : for(unsigned iat=0; iat<n; iat++)alignedref[iat]=-d[iat]+positions[iat]-cp;
1377 0 : return alignedref;
1378 : }
1379 4441 : std::vector<Vector> RMSDCoreData::getAlignedPositionsToReference() {
1380 : std::vector<Vector> alignedpos;
1381 4441 : if(!isInitialized)plumed_merror("getAlignedPostionsToReference needs to initialize the coreData first!");
1382 8882 : const unsigned n=static_cast<unsigned int>(positions.size());
1383 4441 : alignedpos.resize(n);
1384 4441 : Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
1385 : // avoid to calculate matrix element but use the sum of what you have
1386 202996 : for(unsigned iat=0; iat<n; iat++)alignedpos[iat]=matmul(rotation.transpose(),positions[iat]-cp);
1387 4441 : return alignedpos;
1388 : }
1389 :
1390 :
1391 4489 : std::vector<Vector> RMSDCoreData::getCenteredPositions() {
1392 : std::vector<Vector> centeredpos;
1393 8978 : const unsigned n=static_cast<unsigned int>(reference.size());
1394 4489 : centeredpos.resize(n);
1395 4489 : if(!isInitialized)plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1396 : // avoid to calculate matrix element but use the sum of what you have
1397 137339 : for(unsigned iat=0; iat<n; iat++)centeredpos[iat]=positions[iat]-cpositions;
1398 4489 : return centeredpos;
1399 : }
1400 :
1401 4441 : std::vector<Vector> RMSDCoreData::getCenteredReference() {
1402 : std::vector<Vector> centeredref;
1403 8882 : const unsigned n=static_cast<unsigned int>(reference.size());
1404 4441 : centeredref.resize(n);
1405 4441 : if(!isInitialized)plumed_merror("getCenteredReference needs to initialize the coreData first!");
1406 : // avoid to calculate matrix element but use the sum of what you have
1407 4441 : Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
1408 136811 : for(unsigned iat=0; iat<n; iat++)centeredref[iat]=reference[iat]-cr;
1409 4441 : return centeredref;
1410 : }
1411 :
1412 :
1413 48 : Vector RMSDCoreData::getPositionsCenter() {
1414 48 : if(!isInitialized)plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1415 48 : return cpositions;
1416 : }
1417 :
1418 0 : Vector RMSDCoreData::getReferenceCenter() {
1419 0 : if(!isInitialized)plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1420 0 : return creference;
1421 : }
1422 :
1423 1764 : Tensor RMSDCoreData::getRotationMatrixReferenceToPositions() {
1424 1764 : if(!isInitialized)plumed_merror("getRotationMatrixReferenceToPositions needs to initialize the coreData first!");
1425 1764 : return rotation;
1426 : }
1427 :
1428 4489 : Tensor RMSDCoreData::getRotationMatrixPositionsToReference() {
1429 4489 : if(!isInitialized)plumed_merror("getRotationMatrixReferenceToPositions needs to initialize the coreData first!");
1430 4489 : return rotation.transpose();
1431 : }
1432 :
1433 1092 : const std::array<std::array<Tensor,3>,3> & RMSDCoreData::getDRotationDRr01() const {
1434 1092 : if(!isInitialized)plumed_merror("getDRotationDRr01 needs to initialize the coreData first!");
1435 1092 : return drotation_drr01;
1436 : }
1437 :
1438 :
1439 :
1440 : template double RMSD::optimalAlignment<true,true>(const std::vector<double> & align,
1441 : const std::vector<double> & displace,
1442 : const std::vector<Vector> & positions,
1443 : const std::vector<Vector> & reference,
1444 : std::vector<Vector> & derivatives, bool squared)const;
1445 : template double RMSD::optimalAlignment<true,false>(const std::vector<double> & align,
1446 : const std::vector<double> & displace,
1447 : const std::vector<Vector> & positions,
1448 : const std::vector<Vector> & reference,
1449 : std::vector<Vector> & derivatives, bool squared)const;
1450 : template double RMSD::optimalAlignment<false,true>(const std::vector<double> & align,
1451 : const std::vector<double> & displace,
1452 : const std::vector<Vector> & positions,
1453 : const std::vector<Vector> & reference,
1454 : std::vector<Vector> & derivatives, bool squared)const;
1455 : template double RMSD::optimalAlignment<false,false>(const std::vector<double> & align,
1456 : const std::vector<double> & displace,
1457 : const std::vector<Vector> & positions,
1458 : const std::vector<Vector> & reference,
1459 : std::vector<Vector> & derivatives, bool squared)const;
1460 :
1461 :
1462 :
1463 4839 : }
|
