Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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 "Exception.h"
25 : #include <cmath>
26 : #include "Tools.h"
27 :
28 : namespace PLMD {
29 :
30 : // this is ugly, I am sorry
31 : template <typename T>
32 780344 : View<T> make_view(const std::vector<T>& d) {
33 780344 : return View<T> {const_cast<T*>(d.data()),d.size()};
34 : }
35 :
36 47141 : RMSD::RMSD() =default;
37 :
38 : ///
39 : /// 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
40 : /// alignment and beta sets the weight that are used for calculating the displacement.
41 : ///
42 1191 : void RMSD::set(const PDB&pdb,
43 : const std::string & mytype,
44 : const bool remove_center,
45 : const bool normalize_weights ) {
46 :
47 1191 : set(pdb.getOccupancy(),
48 : pdb.getBeta(),
49 : pdb.getPositions(),
50 : mytype,
51 : remove_center,
52 : normalize_weights);
53 :
54 1191 : }
55 28872 : void RMSD::set(const std::vector<double>& alignTo,
56 : const std::vector<double>& displaceFrom,
57 : const std::vector<Vector>& ref,
58 : const std::string & mytype,
59 : const bool remove_center,
60 : const bool normalize_weights ) {
61 :
62 28872 : setReference(ref); // this by default remove the com and assumes uniform weights
63 28872 : setAlign(alignTo, normalize_weights, remove_center); // this recalculates the com with weights. If remove_center=false then it restore the center back
64 28872 : setDisplace(displaceFrom, normalize_weights); // this is does not affect any calculation of the weights
65 28872 : setType(mytype);
66 :
67 28872 : }
68 :
69 74064 : void RMSD::setType(const std::string & mytype) {
70 :
71 74064 : alignmentMethod=SIMPLE; // initialize with the simplest case: no rotation
72 74064 : if (mytype=="SIMPLE") {
73 : alignmentMethod=SIMPLE;
74 72885 : } else if (mytype=="OPTIMAL") {
75 59298 : alignmentMethod=OPTIMAL;
76 13587 : } else if (mytype=="OPTIMAL-FAST") {
77 13587 : alignmentMethod=OPTIMAL_FAST;
78 : } else {
79 0 : plumed_merror("unknown RMSD type" + mytype);
80 : }
81 :
82 74064 : }
83 :
84 27684 : void RMSD::clear() {
85 27684 : reference.clear();
86 : reference_center.zero();
87 27684 : reference_center_is_calculated=false;
88 27684 : reference_center_is_removed=false;
89 27684 : align.clear();
90 27684 : displace.clear();
91 : positions_center.zero();
92 27684 : positions_center_is_calculated=false;
93 27684 : positions_center_is_removed=false;
94 27684 : }
95 :
96 52 : std::string RMSD::getMethod() const {
97 : std::string mystring;
98 52 : switch(alignmentMethod) {
99 2 : case SIMPLE:
100 2 : mystring.assign("SIMPLE");
101 : break;
102 32 : case OPTIMAL:
103 32 : mystring.assign("OPTIMAL");
104 : break;
105 18 : case OPTIMAL_FAST:
106 18 : mystring.assign("OPTIMAL-FAST");
107 : break;
108 : }
109 52 : return mystring;
110 : }
111 : ///
112 : /// this calculates the center of mass for the reference and removes it from the reference itself
113 : /// considering uniform weights for alignment
114 : ///
115 74080 : void RMSD::setReference(const std::vector<Vector>& ref) {
116 74080 : reference=ref;
117 74080 : unsigned n=reference.size();
118 74080 : plumed_massert(align.empty(),"you should first clear() an RMSD object, then set a new reference");
119 74080 : plumed_massert(displace.empty(),"you should first clear() an RMSD object, then set a new reference");
120 74080 : align.resize(n,1.0/n);
121 74080 : displace.resize(n,1.0/n);
122 1063284 : for(unsigned i=0; i<n; i++) {
123 989204 : reference_center+=reference[i]*align[i];
124 : }
125 : #pragma omp simd
126 : for(unsigned i=0; i<n; i++) {
127 989204 : reference[i]-=reference_center;
128 : }
129 74080 : reference_center_is_calculated=true;
130 74080 : reference_center_is_removed=true;
131 74080 : }
132 216081 : const std::vector<Vector>& RMSD::getReference() const {
133 216081 : return reference;
134 : }
135 : ///
136 : /// the alignment weights are here normalized to 1 and the center of the reference is removed accordingly
137 : ///
138 28872 : void RMSD::setAlign(const std::vector<double>& alignTo,
139 : const bool normalize_weights,
140 : const bool remove_center) {
141 28872 : const unsigned n=reference.size();
142 : //NOTE:QUESTION HERE: I am looping on the n from the reference, but I am checking that alignTo and align have the same size, not that reference and align have the same size!!!!
143 28872 : plumed_massert(align.size()==alignTo.size(),"mismatch in dimension of align arrays");
144 28872 : align=alignTo;
145 28872 : if(normalize_weights) {
146 : double w=0.0;
147 57656 : #pragma omp simd reduction(+:w)
148 : for(unsigned i=0; i<n; i++) {
149 400968 : w+=align[i];
150 : }
151 28828 : if(w>epsilon) {
152 28826 : double inv=1.0/w;
153 : #pragma omp simd
154 : for(unsigned i=0; i<n; i++) {
155 400938 : align[i]*=inv;
156 : }
157 : } else {
158 2 : double inv=1.0/n;
159 : #pragma omp simd
160 : for(unsigned i=0; i<n; i++) {
161 30 : align[i]=inv;
162 : }
163 : }
164 : }
165 : // recalculate the center anyway
166 : // just remove the center if that is asked
167 : // if the center was removed before, then add it and store the new one
168 28872 : if(reference_center_is_removed) {
169 28872 : plumed_massert(reference_center_is_calculated," seems that the reference center has been removed but not calculated and stored!");
170 28872 : addCenter(reference,reference_center);
171 : }
172 28872 : reference_center=calculateCenter(reference,align);
173 28872 : reference_center_is_calculated=true;
174 28872 : if(remove_center) {
175 28867 : removeCenter(reference,reference_center);
176 28867 : reference_center_is_removed=true;
177 : } else {
178 5 : reference_center_is_removed=false;
179 : }
180 28872 : }
181 0 : std::vector<double> RMSD::getAlign() {
182 0 : return align;
183 : }
184 : ///
185 : /// here the weigth for normalized weighths are normalized and set
186 : ///
187 28872 : void RMSD::setDisplace(const std::vector<double>& displaceFrom,
188 : const bool normalize_weights) {
189 28872 : const unsigned n=reference.size();
190 : //NOTE:QUESTION HERE: I am looping on the n from the reference, but I am checking that displaceFrom and displace have the same size, not that reference and displace have the same size!!!!
191 28872 : plumed_massert(displace.size()==displaceFrom.size(),"mismatch in dimension of displace arrays");
192 28872 : displace=displaceFrom;
193 28872 : if(normalize_weights) {
194 : double w=0.0;
195 57656 : #pragma omp simd reduction(+:w)
196 : for(unsigned i=0; i<n; i++) {
197 400968 : w+=displace[i];
198 : }
199 28828 : if(w>epsilon) {
200 28826 : double inv=1.0/w;
201 : #pragma omp simd
202 : for(unsigned i=0; i<n; i++) {
203 400938 : displace[i]*=inv;
204 : }
205 : } else {
206 2 : double inv=1.0/n;
207 : #pragma omp simd
208 : for(unsigned i=0; i<n; i++) {
209 30 : displace[i]=inv;
210 : }
211 : }
212 : }
213 28872 : }
214 0 : std::vector<double> RMSD::getDisplace() {
215 0 : return displace;
216 : }
217 : ///
218 : /// This is the main workhorse for rmsd that decides to use specific optimal alignment versions
219 : ///
220 776921 : double RMSD::calculate(const View<Vector> positions,
221 : std::vector<Vector>&derivatives,
222 : const bool squared)const {
223 :
224 : double ret=0.;
225 :
226 776921 : switch(alignmentMethod) {
227 : case SIMPLE : {
228 : // do a simple alignment without rotation
229 277 : std::vector<Vector> displacement( derivatives.size() );
230 277 : ret=simpleAlignment(align,displace,positions,reference,derivatives,displacement,squared);
231 : break;
232 277 : }
233 333674 : case OPTIMAL_FAST : {
234 : // this is calling the fastest option:
235 333674 : if(align==displace) {
236 333672 : ret=optimalAlignment<false,true>(align,displace,positions,reference,derivatives,squared);
237 : } else {
238 2 : ret=optimalAlignment<false,false>(align,displace,positions,reference,derivatives,squared);
239 : }
240 : break;
241 :
242 : }
243 442970 : case OPTIMAL : {
244 : // this is the fast routine but in the "safe" mode, which gives less numerical error:
245 442970 : if(align==displace) {
246 442858 : ret=optimalAlignment<true,true>(align,displace,positions,reference,derivatives,squared);
247 : } else {
248 112 : ret=optimalAlignment<true,false>(align,displace,positions,reference,derivatives,squared);
249 : }
250 : break;
251 : }
252 : }
253 :
254 776921 : return ret;
255 :
256 : }
257 :
258 585289 : double RMSD::calculate(const std::vector<Vector>& positions,
259 : std::vector<Vector>&derivatives,
260 : const bool squared)const {
261 585289 : return calculate(make_view(positions),derivatives,squared);
262 : }
263 :
264 : /// convenience method for calculating the standard derivatives and the derivative of the rmsd respect to the reference position
265 1 : double RMSD::calc_DDistDRef(
266 : const std::vector<Vector>& positions,
267 : std::vector<Vector>&derivatives,
268 : std::vector<Vector>& DDistDRef,
269 : const bool squared ) {
270 : double ret=0.;
271 1 : switch(alignmentMethod) {
272 0 : case SIMPLE:
273 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
274 : break;
275 0 : case OPTIMAL_FAST:
276 0 : if(align==displace) {
277 0 : ret=optimalAlignment_DDistDRef<false,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef, squared);
278 : } else {
279 0 : ret=optimalAlignment_DDistDRef<false,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef,squared);
280 : }
281 : break;
282 1 : case OPTIMAL:
283 1 : if(align==displace) {
284 0 : ret=optimalAlignment_DDistDRef<true,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef,squared);
285 : } else {
286 1 : ret=optimalAlignment_DDistDRef<true,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef,squared);
287 : }
288 : break;
289 : }
290 1 : return ret;
291 :
292 : }
293 :
294 : /// convenience method for calculating the standard derivatives and the derivative of the rmsd respect to the reference position without the matrix contribution
295 : /// as required by SOMA
296 0 : double RMSD::calc_SOMA(const std::vector<Vector>& positions,
297 : std::vector<Vector>&derivatives,
298 : std::vector<Vector>& DDistDRef,
299 : const bool squared) {
300 : double ret=0.;
301 0 : switch(alignmentMethod) {
302 0 : case SIMPLE:
303 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
304 : break;
305 0 : case OPTIMAL_FAST:
306 0 : if(align==displace) {
307 0 : ret=optimalAlignment_SOMA<false,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef, squared);
308 : } else {
309 0 : ret=optimalAlignment_SOMA<false,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef,squared);
310 : }
311 : break;
312 0 : case OPTIMAL:
313 0 : if(align==displace) {
314 0 : ret=optimalAlignment_SOMA<true,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef,squared);
315 : } else {
316 0 : ret=optimalAlignment_SOMA<true,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef,squared);
317 : }
318 : break;
319 : }
320 0 : return ret;
321 :
322 : }
323 :
324 0 : double RMSD::calc_DDistDRef_Rot_DRotDPos(
325 : const std::vector<Vector>& positions,
326 : std::vector<Vector>&derivatives,
327 : std::vector<Vector>& DDistDRef,
328 : Tensor & Rot,
329 : Matrix<std::vector<Vector> >&DRotDPos,
330 : const bool squared) {
331 : double ret=0.;
332 0 : switch(alignmentMethod) {
333 0 : case SIMPLE:
334 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
335 : break;
336 0 : case OPTIMAL_FAST:
337 0 : if(align==displace) {
338 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos<false,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
339 : } else {
340 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos<false,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
341 : }
342 : break;
343 0 : case OPTIMAL:
344 0 : if(align==displace) {
345 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos<true,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
346 : } else {
347 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos<true,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, squared);
348 : }
349 : break;
350 : }
351 0 : return ret;
352 : }
353 :
354 0 : double RMSD::calc_DDistDRef_Rot_DRotDPos_DRotDRef(
355 : const std::vector<Vector>& positions,
356 : std::vector<Vector>&derivatives,
357 : std::vector<Vector>& DDistDRef,
358 : Tensor & Rot,
359 : Matrix<std::vector<Vector> >&DRotDPos,
360 : Matrix<std::vector<Vector> >&DRotDRef,
361 : const bool squared) {
362 : double ret=0.;
363 0 : switch(alignmentMethod) {
364 0 : case SIMPLE:
365 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
366 : break;
367 0 : case OPTIMAL_FAST:
368 0 : if(align==displace) {
369 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<false,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
370 : } else {
371 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<false,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
372 : }
373 : break;
374 0 : case OPTIMAL:
375 0 : if(align==displace) {
376 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<true,true>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
377 : } else {
378 0 : ret=optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef<true,false>(align,displace,make_view(positions),reference,derivatives,DDistDRef, Rot, DRotDPos, DRotDRef, squared);
379 : }
380 : break;
381 : }
382 0 : return ret;
383 : }
384 :
385 1092 : double RMSD::calc_Rot_DRotDRr01( const std::vector<Vector>& positions,
386 : Tensor & Rotation,
387 : std::array<std::array<Tensor,3>,3>& DRotDRr01,
388 : const bool squared) {
389 : double ret=0.;
390 1092 : switch(alignmentMethod) {
391 0 : case SIMPLE:
392 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
393 : break;
394 0 : case OPTIMAL_FAST:
395 0 : if(align==displace) {
396 0 : ret=optimalAlignment_Rot_DRotDRr01<false,true>(align,displace,make_view(positions),reference, Rotation, DRotDRr01, squared);
397 : } else {
398 0 : ret=optimalAlignment_Rot_DRotDRr01<false,false>(align,displace,make_view(positions),reference, Rotation, DRotDRr01, squared);
399 : }
400 : break;
401 1092 : case OPTIMAL:
402 1092 : if(align==displace) {
403 1092 : ret=optimalAlignment_Rot_DRotDRr01<true,true>(align,displace,make_view(positions),reference, Rotation, DRotDRr01, squared);
404 : } else {
405 0 : ret=optimalAlignment_Rot_DRotDRr01<true,false>(align,displace,make_view(positions),reference, Rotation, DRotDRr01, squared);
406 : }
407 : break;
408 : }
409 1092 : return ret;
410 : }
411 :
412 672 : double RMSD::calc_Rot( const std::vector<Vector>& positions,
413 : std::vector<Vector>&derivatives,
414 : Tensor & Rotation,
415 : const bool squared) {
416 : double ret=0.;
417 672 : switch(alignmentMethod) {
418 0 : case SIMPLE:
419 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
420 : break;
421 0 : case OPTIMAL_FAST:
422 0 : if(align==displace) {
423 0 : ret=optimalAlignment_Rot<false,true>(align,displace,make_view(positions),reference,derivatives, Rotation, squared);
424 : } else {
425 0 : ret=optimalAlignment_Rot<false,false>(align,displace,make_view(positions),reference,derivatives, Rotation, squared);
426 : }
427 : break;
428 672 : case OPTIMAL:
429 672 : if(align==displace) {
430 672 : ret=optimalAlignment_Rot<true,true>(align,displace,make_view(positions),reference,derivatives, Rotation, squared);
431 : } else {
432 0 : ret=optimalAlignment_Rot<true,false>(align,displace,make_view(positions),reference,derivatives, Rotation, squared);
433 : }
434 : break;
435 : }
436 672 : return ret;
437 : }
438 :
439 45192 : double RMSD::calculateWithCloseStructure( const std::vector<Vector>& positions,
440 : std::vector<Vector>&derivatives,
441 : const Tensor & rotationPosClose,
442 : const Tensor & rotationRefClose,
443 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01,
444 : const bool squared) {
445 : double ret=0.;
446 45192 : switch(alignmentMethod) {
447 0 : case SIMPLE:
448 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
449 : break;
450 0 : case OPTIMAL_FAST:
451 0 : if(align==displace) {
452 0 : ret=optimalAlignmentWithCloseStructure<false,true>(align,displace,make_view(positions),reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
453 : } else {
454 0 : ret=optimalAlignmentWithCloseStructure<false,false>(align,displace,make_view(positions),reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
455 : }
456 : break;
457 45192 : case OPTIMAL:
458 45192 : if(align==displace) {
459 45192 : ret=optimalAlignmentWithCloseStructure<true,true>(align,displace,make_view(positions),reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
460 : } else {
461 0 : ret=optimalAlignmentWithCloseStructure<true,false>(align,displace,make_view(positions),reference,derivatives, rotationPosClose, rotationRefClose, drotationPosCloseDrr01, squared);
462 : }
463 : break;
464 : }
465 45192 : return ret;
466 : }
467 :
468 147888 : double RMSD::calc_PCAelements(const std::vector<Vector>& positions,
469 : std::vector<Vector>&DDistDPos,
470 : Tensor & Rotation,
471 : Matrix<std::vector<Vector> >& DRotDPos,std::vector<Vector>& alignedpositions,
472 : std::vector<Vector>& centeredpositions,
473 : std::vector<Vector>¢eredreference,
474 : const bool squared) const {
475 : double ret=0.;
476 147888 : switch(alignmentMethod) {
477 0 : case SIMPLE:
478 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
479 : break;
480 123667 : case OPTIMAL_FAST:
481 123667 : if(align==displace) {
482 123667 : ret=optimalAlignment_PCA<false,true>(align,displace,make_view(positions),reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
483 : } else {
484 0 : ret=optimalAlignment_PCA<false,false>(align,displace,make_view(positions),reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
485 : }
486 : break;
487 24221 : case OPTIMAL:
488 24221 : if(align==displace) {
489 24182 : ret=optimalAlignment_PCA<true,true>(align,displace,make_view(positions),reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
490 : } else {
491 39 : ret=optimalAlignment_PCA<true,false>(align,displace,make_view(positions),reference, alignedpositions, centeredpositions,centeredreference,Rotation,DDistDPos,DRotDPos,squared);
492 : }
493 : break;
494 : }
495 147888 : return ret;
496 : }
497 :
498 :
499 210 : double RMSD::calc_FitElements(const std::vector<Vector>& positions,
500 : Tensor & Rotation,
501 : Matrix<std::vector<Vector> >& DRotDPos,
502 : std::vector<Vector>& centeredpositions,
503 : Vector ¢er_positions,
504 : const bool squared) {
505 : double ret=0.;
506 210 : switch(alignmentMethod) {
507 0 : case SIMPLE:
508 0 : plumed_merror("derivative of the refreence frame not implemented for SIMPLE alignmentMethod \n");
509 : break;
510 0 : case OPTIMAL_FAST:
511 0 : if(align==displace) {
512 0 : ret=optimalAlignment_Fit<false,true>(align,displace,make_view(positions),reference, Rotation,DRotDPos,centeredpositions,center_positions,squared);
513 : } else {
514 0 : ret=optimalAlignment_Fit<false,false>(align,displace,make_view(positions),reference, Rotation,DRotDPos,centeredpositions,center_positions,squared);
515 : }
516 : break;
517 210 : case OPTIMAL:
518 210 : if(align==displace) {
519 150 : ret=optimalAlignment_Fit<true,true>(align,displace,make_view(positions),reference,Rotation,DRotDPos,centeredpositions,center_positions,squared);
520 : } else {
521 60 : ret=optimalAlignment_Fit<true,false>(align,displace,make_view(positions),reference,Rotation,DRotDPos,centeredpositions,center_positions,squared);
522 : }
523 : break;
524 : }
525 210 : return ret;
526 : }
527 :
528 47846 : double RMSD::simpleAlignment(const std::vector<double>& alignTo,
529 : const std::vector<double>& displaceFrom,
530 : const View<Vector> positions,
531 : const std::vector<Vector>& ref,
532 : std::vector<Vector>& derivatives,
533 : std::vector<Vector>& displacement,
534 : const bool squared)const {
535 :
536 : double dist(0);
537 47846 : unsigned n=ref.size();
538 :
539 : Vector apositions;
540 : Vector areference;
541 : Vector dpositions;
542 : Vector dreference;
543 :
544 669651 : for(unsigned i=0; i<n; i++) {
545 621805 : double aw=alignTo[i];
546 621805 : double dw=displaceFrom[i];
547 : apositions+=positions[i]*aw;
548 : areference+=ref[i]*aw;
549 : dpositions+=positions[i]*dw;
550 : dreference+=ref[i]*dw;
551 : }
552 :
553 : Vector shift=((apositions-areference)-(dpositions-dreference));
554 669651 : for(unsigned i=0; i<n; i++) {
555 621805 : displacement[i]=(positions[i]-apositions)-(ref[i]-areference);
556 621805 : dist+=displaceFrom[i]*displacement[i].modulo2();
557 621805 : derivatives[i]=2*(displaceFrom[i]*displacement[i]+alignTo[i]*shift);
558 : }
559 :
560 47846 : if(!squared) {
561 : // sqrt
562 1382 : dist=std::sqrt(dist);
563 : ///// sqrt on derivatives
564 20097 : for(unsigned i=0; i<n; i++) {
565 18715 : derivatives[i]*=(0.5/dist);
566 : }
567 : }
568 47846 : return dist;
569 : }
570 :
571 : // this below enable the standard case for rmsd where the rmsd is calculated and the derivative of rmsd respect to positions is retrieved
572 : // additionally this assumes that the com of the reference is already subtracted.
573 : #define OLDRMSD
574 : #ifdef OLDRMSD
575 : // notice that in the current implementation the safe argument only makes sense for
576 : // align==displace
577 : template <bool safe,bool alEqDis>
578 776644 : double RMSD::optimalAlignment(const std::vector<double>& alignTo,
579 : const std::vector<double>& displaceFrom,
580 : const View<Vector> positions,
581 : const std::vector<Vector>& ref,
582 : std::vector<Vector>& derivatives, const bool squared)const {
583 776644 : const unsigned n=ref.size();
584 : // This is the trace of positions*positions + ref*reference
585 : double rr00(0);
586 : double rr11(0);
587 : // This is positions*ref
588 776644 : Tensor rr01;
589 :
590 776644 : derivatives.resize(n);
591 :
592 : Vector cpositions;
593 :
594 : // first expensive loop: compute centers
595 28165818 : for(unsigned iat=0; iat<n; iat++) {
596 27389174 : double w=alignTo[iat];
597 : cpositions+=positions[iat]*w;
598 : }
599 :
600 : // second expensive loop: compute second moments wrt centers
601 28165818 : for(unsigned iat=0; iat<n; iat++) {
602 27389174 : double w=alignTo[iat];
603 27386998 : rr00+=dotProduct(positions[iat]-cpositions,positions[iat]-cpositions)*w;
604 27386998 : rr11+=dotProduct(ref[iat],reference[iat])*w;
605 54778348 : rr01+=Tensor(positions[iat]-cpositions,ref[iat])*w;
606 : }
607 :
608 : Tensor4d m;
609 :
610 776644 : m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
611 776644 : m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
612 776644 : m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
613 776644 : m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
614 776644 : m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
615 776644 : m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
616 776644 : m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
617 776644 : m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
618 776644 : m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
619 776644 : m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
620 776644 : m[1][0] = m[0][1];
621 776644 : m[2][0] = m[0][2];
622 776644 : m[2][1] = m[1][2];
623 776644 : m[3][0] = m[0][3];
624 776644 : m[3][1] = m[1][3];
625 776644 : m[3][2] = m[2][3];
626 :
627 114 : Tensor dm_drr01[4][4];
628 : if(!alEqDis) {
629 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);
630 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);
631 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);
632 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);
633 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);
634 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);
635 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);
636 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);
637 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);
638 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);
639 114 : dm_drr01[1][0] = dm_drr01[0][1];
640 114 : dm_drr01[2][0] = dm_drr01[0][2];
641 114 : dm_drr01[2][1] = dm_drr01[1][2];
642 114 : dm_drr01[3][0] = dm_drr01[0][3];
643 114 : dm_drr01[3][1] = dm_drr01[1][3];
644 114 : dm_drr01[3][2] = dm_drr01[2][3];
645 : }
646 :
647 : double dist=0.0;
648 776530 : Vector4d q;
649 :
650 114 : Tensor dq_drr01[4];
651 : if(!alEqDis) {
652 114 : Vector4d eigenvals;
653 114 : Tensor4d eigenvecs;
654 114 : diagMatSym(m, eigenvals, eigenvecs );
655 114 : dist=eigenvals[0]+rr00+rr11;
656 114 : q=Vector4d(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
657 : double dq_dm[4][4][4];
658 570 : for(unsigned i=0; i<4; i++)
659 2280 : for(unsigned j=0; j<4; j++)
660 9120 : for(unsigned k=0; k<4; k++) {
661 : double tmp=0.0;
662 : // perturbation theory for matrix m
663 29184 : for(unsigned l=1; l<4; l++) {
664 21888 : tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
665 : }
666 7296 : dq_dm[i][j][k]=tmp;
667 : }
668 : // propagation to _drr01
669 570 : for(unsigned i=0; i<4; i++) {
670 456 : Tensor tmp;
671 2280 : for(unsigned j=0; j<4; j++)
672 9120 : for(unsigned k=0; k<4; k++) {
673 14592 : tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
674 : }
675 456 : dq_drr01[i]=tmp;
676 : }
677 : } else {
678 776530 : dist = lowestEigenpairSym(m, q) + rr00 + rr11;
679 : }
680 :
681 :
682 : // This is the rotation matrix that brings ref to positions
683 : // i.e. matmul(rotation,ref[iat])+shift is fitted to positions[iat]
684 :
685 : Tensor rotation;
686 776644 : rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
687 776644 : rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
688 776644 : rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
689 776644 : rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
690 776644 : rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
691 776644 : rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
692 776644 : rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
693 776644 : rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
694 776644 : rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
695 :
696 :
697 114 : std::array<std::array<Tensor,3>,3> drotation_drr01;
698 : if(!alEqDis) {
699 228 : 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];
700 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];
701 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];
702 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]));
703 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]));
704 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]));
705 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]));
706 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]));
707 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]));
708 : }
709 :
710 : double prefactor=2.0;
711 :
712 776530 : if(!squared && alEqDis) {
713 233414 : prefactor*=0.5/std::sqrt(dist);
714 : }
715 :
716 : // if "safe", recompute dist here to a better accuracy
717 : if(safe || !alEqDis) {
718 : dist=0.0;
719 : }
720 :
721 : // If safe is set to "false", MSD is taken from the eigenvalue of the M matrix
722 : // If safe is set to "true", MSD is recomputed from the rotational matrix
723 : // For some reason, this last approach leads to less numerical noise but adds an overhead
724 :
725 114 : Tensor ddist_drotation;
726 : Vector ddist_dcpositions;
727 :
728 : // third expensive loop: derivatives
729 28165818 : for(unsigned iat=0; iat<n; iat++) {
730 27389174 : Vector d(positions[iat]-cpositions - matmul(rotation,ref[iat]));
731 : if(alEqDis) {
732 : // there is no need for derivatives of rotation and shift here as it is by construction zero
733 : // (similar to Hellman-Feynman forces)
734 27386998 : derivatives[iat]= prefactor*alignTo[iat]*d;
735 : if(safe) {
736 7454522 : dist+=alignTo[iat]*modulo2(d);
737 : }
738 : } else {
739 : // the case for alignTo != displaceFrom is different, sob:
740 2176 : dist+=displaceFrom[iat]*modulo2(d);
741 : // these are the derivatives assuming the roto-translation as frozen
742 2176 : derivatives[iat]=2*displaceFrom[iat]*d;
743 : // here I accumulate derivatives wrt rotation matrix ..
744 4352 : ddist_drotation+=-2*displaceFrom[iat]*extProduct(d,ref[iat]);
745 : // .. and cpositions
746 : ddist_dcpositions+=-2*displaceFrom[iat]*d;
747 : }
748 : }
749 :
750 : if(!alEqDis) {
751 114 : Tensor ddist_drr01;
752 456 : for(unsigned i=0; i<3; i++)
753 1368 : for(unsigned j=0; j<3; j++) {
754 2052 : ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
755 : }
756 2290 : for(unsigned iat=0; iat<n; iat++) {
757 : // this is propagating to positions.
758 : // I am implicitly using the derivative of rr01 wrt positions here
759 2176 : derivatives[iat]+=matmul(ddist_drr01,ref[iat])*alignTo[iat];
760 2176 : derivatives[iat]+=ddist_dcpositions*alignTo[iat];
761 : }
762 : }
763 776644 : if(!squared) {
764 233528 : dist=std::sqrt(dist);
765 : if(!alEqDis) {
766 114 : double xx=0.5/dist;
767 2290 : for(unsigned iat=0; iat<n; iat++) {
768 2176 : derivatives[iat]*=xx;
769 : }
770 : }
771 : }
772 :
773 776644 : return dist;
774 : }
775 : #else
776 : /// note that this method is intended to be repeatedly invoked
777 : /// when the reference does already have the center subtracted
778 : /// but the position has not calculated center and not subtracted
779 : template <bool safe,bool alEqDis>
780 : double RMSD::optimalAlignment(const std::vector<double>& align,
781 : const std::vector<double>& displace,
782 : const View<Vector> positions,
783 : const std::vector<Vector>& reference,
784 : std::vector<Vector>& derivatives,
785 : const bool squared) const {
786 : RMSDCoreData cd(align,displace,make_view(positions),reference);
787 :
788 : // transfer the settings for the center to let the CoreCalc deal with it
789 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
790 : if(positions_center_is_calculated) {
791 : cd.setPositionsCenter(positions_center);
792 : } else {
793 : cd.calcPositionsCenter();
794 : };
795 :
796 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
797 : if(!reference_center_is_calculated) {
798 : cd.calcReferenceCenter();
799 : } else {
800 : cd.setReferenceCenter(reference_center);
801 : }
802 :
803 : // Perform the diagonalization and all the needed stuff
804 : cd.doCoreCalc(safe,alEqDis);
805 : // make the core calc distance
806 : double dist=cd.getDistance(squared);
807 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
808 : derivatives=cd.getDDistanceDPositions();
809 : return dist;
810 : }
811 : #endif
812 : template <bool safe,bool alEqDis>
813 1 : double RMSD::optimalAlignment_DDistDRef(const std::vector<double>& alignTo,
814 : const std::vector<double>& displaceFrom,
815 : const View<Vector> positions,
816 : const std::vector<Vector>& ref,
817 : std::vector<Vector>& derivatives,
818 : std::vector<Vector>& ddistdref,
819 : const bool squared) const {
820 : //initialize the data into the structure
821 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
822 1 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
823 : // transfer the settings for the center to let the CoreCalc deal with it
824 : // transfer the settings for the center to let the CoreCalc deal with it
825 1 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
826 1 : if(positions_center_is_calculated) {
827 0 : cd.setPositionsCenter(positions_center);
828 : } else {
829 1 : cd.calcPositionsCenter();
830 : };
831 :
832 1 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
833 1 : if(!reference_center_is_calculated) {
834 0 : cd.calcReferenceCenter();
835 : } else {
836 1 : cd.setReferenceCenter(reference_center);
837 : }
838 :
839 : // Perform the diagonalization and all the needed stuff
840 1 : cd.doCoreCalc(safe,alEqDis);
841 : // make the core calc distance
842 1 : double dist=cd.getDistance(squared);
843 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
844 2 : derivatives=cd.getDDistanceDPositions();
845 2 : ddistdref=cd.getDDistanceDReference();
846 1 : return dist;
847 : }
848 :
849 : template <bool safe,bool alEqDis>
850 0 : double RMSD::optimalAlignment_SOMA(const std::vector<double>& alignTo,
851 : const std::vector<double>& displaceFrom,
852 : const View<Vector> positions,
853 : const std::vector<Vector>& ref,
854 : std::vector<Vector>& derivatives,
855 : std::vector<Vector>& ddistdref,
856 : const bool squared) const {
857 : //initialize the data into the structure
858 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
859 0 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
860 : // transfer the settings for the center to let the CoreCalc deal with it
861 : // transfer the settings for the center to let the CoreCalc deal with it
862 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
863 0 : if(positions_center_is_calculated) {
864 0 : cd.setPositionsCenter(positions_center);
865 : } else {
866 0 : cd.calcPositionsCenter();
867 : };
868 :
869 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
870 0 : if(!reference_center_is_calculated) {
871 0 : cd.calcReferenceCenter();
872 : } else {
873 0 : cd.setReferenceCenter(reference_center);
874 : }
875 :
876 : // Perform the diagonalization and all the needed stuff
877 0 : cd.doCoreCalc(safe,alEqDis);
878 : // make the core calc distance
879 0 : double dist=cd.getDistance(squared);
880 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
881 0 : derivatives=cd.getDDistanceDPositions();
882 0 : ddistdref=cd.getDDistanceDReferenceSOMA();
883 0 : return dist;
884 : }
885 :
886 :
887 : template <bool safe,bool alEqDis>
888 0 : double RMSD::optimalAlignment_DDistDRef_Rot_DRotDPos(
889 : const std::vector<double>& alignTo,
890 : const std::vector<double>& displaceFrom,
891 : const View<Vector> positions,
892 : const std::vector<Vector>& ref,
893 : std::vector<Vector>& derivatives,
894 : std::vector<Vector>& ddistdref,
895 : Tensor & Rotation,
896 : Matrix<std::vector<Vector> >&DRotDPos,
897 : const bool squared) const {
898 : //initialize the data into the structure
899 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
900 0 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
901 : // transfer the settings for the center to let the CoreCalc deal with it
902 : // transfer the settings for the center to let the CoreCalc deal with it
903 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
904 0 : if(positions_center_is_calculated) {
905 0 : cd.setPositionsCenter(positions_center);
906 : } else {
907 0 : cd.calcPositionsCenter();
908 : };
909 :
910 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
911 0 : if(!reference_center_is_calculated) {
912 0 : cd.calcReferenceCenter();
913 : } else {
914 0 : cd.setReferenceCenter(reference_center);
915 : }
916 :
917 : // Perform the diagonalization and all the needed stuff
918 0 : cd.doCoreCalc(safe,alEqDis);
919 : // make the core calc distance
920 0 : double dist=cd.getDistance(squared);
921 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
922 0 : derivatives=cd.getDDistanceDPositions();
923 0 : ddistdref=cd.getDDistanceDReference();
924 : // get the rotation matrix
925 0 : Rotation=cd.getRotationMatrixReferenceToPositions();
926 : // get its derivative
927 0 : DRotDPos=cd.getDRotationDPositions();
928 0 : return dist;
929 : }
930 :
931 : template <bool safe,bool alEqDis>
932 0 : double RMSD::optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef(
933 : const std::vector<double>& alignTo,
934 : const std::vector<double>& displaceFrom,
935 : const View<Vector> positions,
936 : const std::vector<Vector>& ref,
937 : std::vector<Vector>& derivatives,
938 : std::vector<Vector>& ddistdref,
939 : Tensor & Rotation,
940 : Matrix<std::vector<Vector> >&DRotDPos,
941 : Matrix<std::vector<Vector> >&DRotDRef,
942 : const bool squared) const {
943 : //initialize the data into the structure
944 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
945 0 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
946 : // transfer the settings for the center to let the CoreCalc deal with it
947 : // transfer the settings for the center to let the CoreCalc deal with it
948 0 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
949 0 : if(positions_center_is_calculated) {
950 0 : cd.setPositionsCenter(positions_center);
951 : } else {
952 0 : cd.calcPositionsCenter();
953 : };
954 :
955 0 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
956 0 : if(!reference_center_is_calculated) {
957 0 : cd.calcReferenceCenter();
958 : } else {
959 0 : cd.setReferenceCenter(reference_center);
960 : }
961 :
962 : // Perform the diagonalization and all the needed stuff
963 0 : cd.doCoreCalc(safe,alEqDis);
964 : // make the core calc distance
965 0 : double dist=cd.getDistance(squared);
966 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
967 0 : derivatives=cd.getDDistanceDPositions();
968 0 : ddistdref=cd.getDDistanceDReference();
969 : // get the rotation matrix
970 0 : Rotation=cd.getRotationMatrixReferenceToPositions();
971 : // get its derivative
972 0 : DRotDPos=cd.getDRotationDPositions();
973 0 : DRotDRef=cd.getDRotationDReference();
974 0 : return dist;
975 : }
976 :
977 : template <bool safe,bool alEqDis>
978 1092 : double RMSD::optimalAlignment_Rot_DRotDRr01(const std::vector<double>& alignTo,
979 : const std::vector<double>& displaceFrom,
980 : const View<Vector> positions,
981 : const std::vector<Vector>& ref,
982 : Tensor & Rotation,
983 : std::array<std::array<Tensor,3>,3>& DRotDRr01,
984 : const bool squared) const {
985 : //initialize the data into the structure
986 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
987 1092 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
988 : // transfer the settings for the center to let the CoreCalc deal with it
989 1092 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
990 1092 : if(positions_center_is_calculated) {
991 0 : cd.setPositionsCenter(positions_center);
992 : } else {
993 1092 : cd.calcPositionsCenter();
994 : };
995 :
996 1092 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
997 1092 : if(!reference_center_is_calculated) {
998 0 : cd.calcReferenceCenter();
999 : } else {
1000 1092 : cd.setReferenceCenter(reference_center);
1001 : }
1002 :
1003 : // Perform the diagonalization and all the needed stuff
1004 1092 : cd.doCoreCalc(safe,alEqDis);
1005 : // make the core calc distance
1006 1092 : double dist=cd.getDistance(squared);
1007 : // get the rotation matrix
1008 1092 : Rotation=cd.getRotationMatrixReferenceToPositions();
1009 : //get detivative w.r.t. rr01
1010 1092 : DRotDRr01=cd.getDRotationDRr01();
1011 1092 : return dist;
1012 : }
1013 :
1014 : template <bool safe,bool alEqDis>
1015 672 : double RMSD::optimalAlignment_Rot(const std::vector<double>& alignTo,
1016 : const std::vector<double>& displaceFrom,
1017 : const View<Vector> positions,
1018 : const std::vector<Vector>& ref,
1019 : std::vector<Vector>& derivatives,
1020 : Tensor & Rotation,
1021 : const bool squared) const {
1022 : //initialize the data into the structure
1023 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
1024 672 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
1025 : // transfer the settings for the center to let the CoreCalc deal with it
1026 672 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
1027 672 : if(positions_center_is_calculated) {
1028 0 : cd.setPositionsCenter(positions_center);
1029 : } else {
1030 672 : cd.calcPositionsCenter();
1031 : };
1032 :
1033 672 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
1034 672 : if(!reference_center_is_calculated) {
1035 0 : cd.calcReferenceCenter();
1036 : } else {
1037 672 : cd.setReferenceCenter(reference_center);
1038 : }
1039 :
1040 : // Perform the diagonalization and all the needed stuff
1041 672 : cd.doCoreCalc(safe,alEqDis);
1042 : // make the core calc distance
1043 672 : double dist=cd.getDistance(squared);
1044 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
1045 1344 : derivatives=cd.getDDistanceDPositions();
1046 : // get the rotation matrix
1047 672 : Rotation=cd.getRotationMatrixReferenceToPositions();
1048 672 : return dist;
1049 : }
1050 :
1051 : template <bool safe,bool alEqDis>
1052 45192 : double RMSD::optimalAlignmentWithCloseStructure(const std::vector<double>& alignTo,
1053 : const std::vector<double>& displaceFrom,
1054 : const View<Vector> positions,
1055 : const std::vector<Vector>& ref,
1056 : std::vector<Vector>& derivatives,
1057 : const Tensor & rotationPosClose,
1058 : const Tensor & rotationRefClose,
1059 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01,
1060 : const bool squared) const {
1061 : //initialize the data into the structure
1062 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
1063 45192 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
1064 : // transfer the settings for the center to let the CoreCalc deal with it
1065 45192 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
1066 45192 : if(positions_center_is_calculated) {
1067 0 : cd.setPositionsCenter(positions_center);
1068 : } else {
1069 45192 : cd.calcPositionsCenter();
1070 : };
1071 :
1072 45192 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
1073 45192 : if(!reference_center_is_calculated) {
1074 0 : cd.calcReferenceCenter();
1075 : } else {
1076 45192 : cd.setReferenceCenter(reference_center);
1077 : }
1078 :
1079 : // instead of diagonalization, approximate with saved rotation matrix
1080 45192 : cd.doCoreCalcWithCloseStructure(safe,alEqDis, rotationPosClose, rotationRefClose, drotationPosCloseDrr01);
1081 : // make the core calc distance
1082 45192 : double dist=cd.getDistance(squared);
1083 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
1084 90384 : derivatives=cd.getDDistanceDPositions();
1085 45192 : return dist;
1086 : }
1087 :
1088 :
1089 : template <bool safe,bool alEqDis>
1090 147888 : double RMSD::optimalAlignment_PCA(const std::vector<double>& alignTo,
1091 : const std::vector<double>& displaceFrom,
1092 : const View<Vector> positions,
1093 : const std::vector<Vector>& ref,
1094 : std::vector<Vector>& alignedpositions,
1095 : std::vector<Vector>& centeredpositions,
1096 : std::vector<Vector>& centeredreference,
1097 : Tensor & Rotation,
1098 : std::vector<Vector>& DDistDPos,
1099 : Matrix<std::vector<Vector> >& DRotDPos,
1100 : const bool squared) const {
1101 : //initialize the data into the structure
1102 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
1103 147888 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
1104 : // transfer the settings for the center to let the CoreCalc deal with it
1105 147888 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
1106 147888 : if(positions_center_is_calculated) {
1107 0 : cd.setPositionsCenter(positions_center);
1108 : } else {
1109 147888 : cd.calcPositionsCenter();
1110 : };
1111 :
1112 147888 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
1113 147888 : if(!reference_center_is_calculated) {
1114 0 : cd.calcReferenceCenter();
1115 : } else {
1116 147888 : cd.setReferenceCenter(reference_center);
1117 : }
1118 :
1119 : // Perform the diagonalization and all the needed stuff
1120 147888 : cd.doCoreCalc(safe,alEqDis);
1121 : // make the core calc distance
1122 147888 : double dist=cd.getDistance(squared);
1123 : // make the derivatives by using pieces calculated in coreCalc (probably the best is just to copy the vector...)
1124 295776 : DDistDPos=cd.getDDistanceDPositions();
1125 : // get the rotation matrix
1126 147888 : Rotation=cd.getRotationMatrixPositionsToReference();
1127 : // get its derivative
1128 147888 : DRotDPos=cd.getDRotationDPositions(true); // this gives back the inverse
1129 : // get aligned positions
1130 295776 : alignedpositions=cd.getAlignedPositionsToReference();
1131 : // get centered positions
1132 295776 : centeredpositions=cd.getCenteredPositions();
1133 : // get centered ref
1134 295776 : centeredreference=cd.getCenteredReference();
1135 147888 : return dist;
1136 : }
1137 :
1138 :
1139 : template <bool safe,bool alEqDis>
1140 210 : double RMSD::optimalAlignment_Fit(const std::vector<double>& alignTo,
1141 : const std::vector<double>& displaceFrom,
1142 : const View<Vector> positions,
1143 : const std::vector<Vector>& ref,
1144 : Tensor& Rotation,
1145 : Matrix<std::vector<Vector>>& DRotDPos,
1146 : std::vector<Vector>& centeredpositions,
1147 : Vector& center_positions,
1148 : const bool squared) {
1149 : //initialize the data into the structure
1150 : // typically the positions do not have the com neither calculated nor subtracted. This layer takes care of this business
1151 210 : RMSDCoreData cd(alignTo,displaceFrom,positions,ref);
1152 : // transfer the settings for the center to let the CoreCalc deal with it
1153 210 : cd.setPositionsCenterIsRemoved(positions_center_is_removed);
1154 210 : if(positions_center_is_calculated) {
1155 0 : cd.setPositionsCenter(positions_center);
1156 : } else {
1157 210 : cd.calcPositionsCenter();
1158 : };
1159 :
1160 210 : cd.setReferenceCenterIsRemoved(reference_center_is_removed);
1161 210 : if(!reference_center_is_calculated) {
1162 0 : cd.calcReferenceCenter();
1163 : } else {
1164 210 : cd.setReferenceCenter(reference_center);
1165 : }
1166 :
1167 : // Perform the diagonalization and all the needed stuff
1168 210 : cd.doCoreCalc(safe,alEqDis);
1169 : // make the core calc distance
1170 210 : double dist=cd.getDistance(squared);
1171 : // get the rotation matrix
1172 210 : Rotation=cd.getRotationMatrixPositionsToReference();
1173 : // get its derivative
1174 210 : DRotDPos=cd.getDRotationDPositions(true); // this gives back the inverse
1175 : // get centered positions
1176 420 : centeredpositions=cd.getCenteredPositions();
1177 : // get center
1178 210 : center_positions=cd.getPositionsCenter();
1179 210 : return dist;
1180 : }
1181 :
1182 : /// This calculates the elements needed by the quaternion to calculate everything that is needed
1183 : /// additional calls retrieve different components
1184 : /// note that this considers that the centers of both reference and positions are already set
1185 : /// but automatically should properly account for non removed components: if not removed then it
1186 : /// removes prior to calculation of the alignment
1187 149863 : void RMSDCoreData::doCoreCalc(const bool setSafe,
1188 : const bool setAlEqDis) {
1189 :
1190 : // transfer this bools to the cd so that this settings will be reflected in the other calls
1191 149863 : alEqDis=setAlEqDis;
1192 149863 : safe=setSafe;
1193 149863 : const unsigned n=static_cast<unsigned int>(reference.size());
1194 :
1195 149863 : plumed_massert(creference_is_calculated,"the center of the reference frame must be already provided at this stage");
1196 149863 : plumed_massert(cpositions_is_calculated,"the center of the positions frame must be already provided at this stage");
1197 :
1198 : // This is the trace of positions*positions + reference*reference
1199 149863 : rr00=0.;
1200 149863 : rr11=0.;
1201 : // This is positions*reference
1202 149863 : Tensor rr01;
1203 : // center of mass managing: must subtract the center from the position or not?
1204 : Vector cp;
1205 : cp.zero();
1206 149863 : if(!cpositions_is_removed) {
1207 149863 : cp=cpositions;
1208 : }
1209 : Vector cr;
1210 : cr.zero();
1211 149863 : if(!creference_is_removed) {
1212 60 : cr=creference;
1213 : }
1214 : // second expensive loop: compute second moments wrt centers
1215 2583539 : for(unsigned iat=0; iat<n; iat++) {
1216 2433676 : double w=align[iat];
1217 2433676 : rr00+=dotProduct(positions[iat]-cp,positions[iat]-cp)*w;
1218 2433676 : rr11+=dotProduct(reference[iat]-cr,reference[iat]-cr)*w;
1219 4867352 : rr01+=Tensor(positions[iat]-cp,reference[iat]-cr)*w;
1220 : }
1221 :
1222 : // the quaternion matrix: this is internal
1223 : Tensor4d m;
1224 :
1225 149863 : m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
1226 149863 : m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
1227 149863 : m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
1228 149863 : m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
1229 149863 : m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
1230 149863 : m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
1231 149863 : m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
1232 149863 : m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
1233 149863 : m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
1234 149863 : m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
1235 149863 : m[1][0] = m[0][1];
1236 149863 : m[2][0] = m[0][2];
1237 149863 : m[2][1] = m[1][2];
1238 149863 : m[3][0] = m[0][3];
1239 149863 : m[3][1] = m[1][3];
1240 149863 : m[3][2] = m[2][3];
1241 :
1242 :
1243 149863 : Tensor dm_drr01[4][4];
1244 149863 : 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);
1245 149863 : 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);
1246 149863 : 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);
1247 149863 : 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);
1248 149863 : 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);
1249 149863 : 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);
1250 149863 : 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);
1251 149863 : 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);
1252 149863 : 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);
1253 149863 : 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);
1254 149863 : dm_drr01[1][0] = dm_drr01[0][1];
1255 149863 : dm_drr01[2][0] = dm_drr01[0][2];
1256 149863 : dm_drr01[2][1] = dm_drr01[1][2];
1257 149863 : dm_drr01[3][0] = dm_drr01[0][3];
1258 149863 : dm_drr01[3][1] = dm_drr01[1][3];
1259 149863 : dm_drr01[3][2] = dm_drr01[2][3];
1260 :
1261 :
1262 : Vector4d q;
1263 :
1264 149863 : Tensor dq_drr01[4];
1265 149863 : diagMatSym(m, eigenvals, eigenvecs );
1266 149863 : q=Vector4d(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
1267 : double dq_dm[4][4][4];
1268 749315 : for(unsigned i=0; i<4; i++)
1269 2997260 : for(unsigned j=0; j<4; j++)
1270 11989040 : for(unsigned k=0; k<4; k++) {
1271 : double tmp=0.0;
1272 : // perturbation theory for matrix m
1273 38364928 : for(unsigned l=1; l<4; l++) {
1274 28773696 : tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
1275 : }
1276 9591232 : dq_dm[i][j][k]=tmp;
1277 : }
1278 : // propagation to _drr01
1279 749315 : for(unsigned i=0; i<4; i++) {
1280 599452 : Tensor tmp;
1281 2997260 : for(unsigned j=0; j<4; j++)
1282 11989040 : for(unsigned k=0; k<4; k++) {
1283 19182464 : tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
1284 : }
1285 599452 : dq_drr01[i]=tmp;
1286 : }
1287 :
1288 : // This is the rotation matrix that brings reference to positions
1289 : // i.e. matmul(rotation,reference[iat])+shift is fitted to positions[iat]
1290 :
1291 149863 : rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
1292 149863 : rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
1293 149863 : rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
1294 149863 : rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
1295 149863 : rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
1296 149863 : rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
1297 149863 : rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
1298 149863 : rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
1299 149863 : rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
1300 :
1301 :
1302 299726 : 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];
1303 149863 : 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];
1304 149863 : 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];
1305 149863 : 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]));
1306 149863 : 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]));
1307 149863 : 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]));
1308 149863 : 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]));
1309 149863 : 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]));
1310 149863 : 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]));
1311 :
1312 149863 : d.resize(n);
1313 :
1314 : // calculate rotation matrix derivatives and components distances needed for components only when align!=displacement
1315 149863 : if(!alEqDis) {
1316 : ddist_drotation.zero();
1317 : }
1318 : // This pragma leads to incorrect results with INTEL compiler.
1319 : // Failures are seen in rt65-rmsd2, rt-close-structure, rt64-pca, and others.
1320 : // Not really clear why. GB
1321 : // #pragma omp simd
1322 2583539 : for(unsigned iat=0; iat<n; iat++) {
1323 : // components differences: this is useful externally
1324 2433676 : d[iat]=positions[iat]-cp - matmul(rotation,reference[iat]-cr);
1325 : }
1326 : // ddist_drotation if needed
1327 2583539 : for (unsigned iat=0; iat<n; iat++) {
1328 4867352 : ddist_drotation+=-2*displace[iat]*extProduct(d[iat],reference[iat]-cr);
1329 : }
1330 :
1331 : ddist_drr01.zero();
1332 599452 : for(unsigned i=0; i<3; i++)
1333 1798356 : for(unsigned j=0; j<3; j++) {
1334 2697534 : ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
1335 : }
1336 149863 : isInitialized=true;
1337 :
1338 149863 : }
1339 : /// just retrieve the distance already calculated
1340 195055 : double RMSDCoreData::getDistance( const bool squared) {
1341 :
1342 195055 : if(!isInitialized) {
1343 0 : plumed_merror("getDistance cannot calculate the distance without being initialized first by doCoreCalc ");
1344 : }
1345 :
1346 : double localDist=0.0;
1347 195055 : const unsigned n=static_cast<unsigned int>(reference.size());
1348 195055 : if(safe || !alEqDis) {
1349 : localDist=0.0;
1350 : } else {
1351 123667 : localDist=eigenvals[0]+rr00+rr11;
1352 : }
1353 195055 : if(alEqDis) {
1354 194955 : if(safe) {
1355 142576 : #pragma omp simd reduction(+:localDist)
1356 : for(unsigned iat=0; iat<n; iat++) {
1357 1411839 : localDist+=align[iat]*modulo2(d[iat]);
1358 : }
1359 : }
1360 : } else {
1361 200 : #pragma omp simd reduction(+:localDist)
1362 : for(unsigned iat=0; iat<n; iat++) {
1363 1662 : localDist+=displace[iat]*modulo2(d[iat]);
1364 : }
1365 : }
1366 195055 : if(!squared) {
1367 23679 : dist=std::sqrt(localDist);
1368 23679 : distanceIsMSD=false;
1369 : } else {
1370 171376 : dist=localDist;
1371 171376 : distanceIsMSD=true;
1372 : }
1373 195055 : hasDistance=true;
1374 195055 : return dist;
1375 : }
1376 :
1377 45192 : void RMSDCoreData::doCoreCalcWithCloseStructure(const bool setSafe,
1378 : const bool setAlEqDis,
1379 : const Tensor & rotationPosClose,
1380 : const Tensor & rotationRefClose,
1381 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01) {
1382 :
1383 45192 : alEqDis=setAlEqDis;
1384 45192 : safe=setSafe;
1385 45192 : unsigned natoms = reference.size();
1386 45192 : Tensor ddist_drxy;
1387 : ddist_drr01.zero();
1388 45192 : d.resize(natoms);
1389 :
1390 : // center of mass managing: must subtract the center from the position or not?
1391 : Vector cp;
1392 : cp.zero();
1393 45192 : if(!cpositions_is_removed) {
1394 45192 : cp=cpositions;
1395 : }
1396 : Vector cr;
1397 : cr.zero();
1398 45192 : if(!creference_is_removed) {
1399 0 : cr=creference;
1400 : }
1401 : //distance = \sum_{n=0}^{N} w_n(x_n-cpos-R_{XY} R_{AY} a_n)^2
1402 :
1403 45192 : Tensor rot = matmul(rotationPosClose, rotationRefClose);
1404 :
1405 : // This pragma leads to incorrect results with INTEL compiler.
1406 : // Failures are seen in rt65-rmsd2, rt-close-structure, rt64-pca, and others.
1407 : // Not really clear why. GB
1408 : // #pragma omp simd
1409 632688 : for (unsigned iat=0; iat<natoms; iat++) {
1410 587496 : d[iat] = positions[iat] - cp - matmul(rot, reference[iat]-cr);
1411 : }
1412 45192 : if (!alEqDis) {
1413 0 : for (unsigned iat=0; iat<natoms; iat++) {
1414 : //dist = \sum w_i(x_i - cpos - R_xy * R_ay * a_i)
1415 0 : ddist_drxy += -2*displace[iat]*extProduct(matmul(d[iat], rotationRefClose), reference[iat]-cr);
1416 : }
1417 0 : for(unsigned i=0; i<3; i++)
1418 0 : for(unsigned j=0; j<3; j++) {
1419 0 : ddist_drr01+=ddist_drxy[i][j]*drotationPosCloseDrr01[i][j];
1420 : }
1421 : }
1422 45192 : isInitialized=true;
1423 45192 : }
1424 :
1425 193753 : std::vector<Vector> RMSDCoreData::getDDistanceDPositions() {
1426 193753 : std::vector<Vector> derivatives;
1427 193753 : const unsigned n=static_cast<unsigned int>(reference.size());
1428 : Vector ddist_dcpositions;
1429 193753 : derivatives.resize(n);
1430 : double prefactor=1.0;
1431 193753 : if(!distanceIsMSD) {
1432 23619 : prefactor*=0.5/dist;
1433 : }
1434 193753 : if(!hasDistance) {
1435 0 : plumed_merror("getDPositionsDerivatives needs to calculate the distance via getDistance first !");
1436 : }
1437 193753 : if(!isInitialized) {
1438 0 : plumed_merror("getDPositionsDerivatives needs to initialize the coreData first!");
1439 : }
1440 : Vector csum;
1441 2714729 : for(unsigned iat=0; iat<n; iat++) {
1442 2520976 : if(alEqDis) {
1443 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1444 : // (similar to Hellman-Feynman forces)
1445 2519614 : derivatives[iat]= 2*prefactor*align[iat]*d[iat];
1446 : } else {
1447 : // these are the derivatives assuming the roto-translation as frozen
1448 1362 : Vector tmp1=2*displace[iat]*d[iat];
1449 1362 : derivatives[iat]=tmp1;
1450 : // derivative of cpositions
1451 : ddist_dcpositions+=-tmp1;
1452 : // these needed for com corrections
1453 1362 : Vector tmp2=matmul(ddist_drr01,reference[iat]-creference)*align[iat];
1454 : derivatives[iat]+=tmp2;
1455 : csum+=tmp2;
1456 : }
1457 : }
1458 :
1459 193753 : if(!alEqDis) {
1460 : const Vector shift=ddist_dcpositions-csum;
1461 : #pragma omp simd
1462 : for(unsigned iat=0; iat<n; iat++) {
1463 1362 : const double a=align[iat];
1464 1362 : derivatives[iat][0]=prefactor*(derivatives[iat][0]+shift[0]*a);
1465 1362 : derivatives[iat][1]=prefactor*(derivatives[iat][1]+shift[1]*a);
1466 1362 : derivatives[iat][2]=prefactor*(derivatives[iat][2]+shift[2]*a);
1467 : }
1468 : }
1469 :
1470 193753 : return derivatives;
1471 0 : }
1472 :
1473 1 : std::vector<Vector> RMSDCoreData::getDDistanceDReference() {
1474 1 : std::vector<Vector> derivatives;
1475 1 : const unsigned n=static_cast<unsigned int>(reference.size());
1476 : Vector ddist_dcreference;
1477 1 : derivatives.resize(n);
1478 : double prefactor=1.0;
1479 1 : if(!distanceIsMSD) {
1480 1 : prefactor*=0.5/dist;
1481 : }
1482 : Vector csum;
1483 :
1484 1 : if(!hasDistance) {
1485 0 : plumed_merror("getDDistanceDReference needs to calculate the distance via getDistance first !");
1486 : }
1487 1 : if(!isInitialized) {
1488 0 : plumed_merror("getDDistanceDReference to initialize the coreData first!");
1489 : }
1490 : // get the transpose rotation
1491 1 : Tensor t_rotation=rotation.transpose();
1492 1 : Tensor t_ddist_drr01=ddist_drr01.transpose();
1493 :
1494 : // third expensive loop: derivatives
1495 856 : for(unsigned iat=0; iat<n; iat++) {
1496 855 : if(alEqDis) {
1497 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1498 : // (similar to Hellman-Feynman forces)
1499 : //TODO: check this derivative down here
1500 0 : derivatives[iat]= -2*prefactor*align[iat]*matmul(t_rotation,d[iat]);
1501 : } else {
1502 : // these are the derivatives assuming the roto-translation as frozen
1503 855 : Vector tmp1=2*displace[iat]*matmul(t_rotation,d[iat]);
1504 855 : derivatives[iat]= -tmp1;
1505 : // derivative of cpositions
1506 : ddist_dcreference+=tmp1;
1507 : // these below are needed for com correction
1508 855 : Vector tmp2=matmul(t_ddist_drr01,positions[iat]-cpositions)*align[iat];
1509 : derivatives[iat]+=tmp2;
1510 : csum+=tmp2;
1511 : }
1512 : }
1513 :
1514 1 : if(!alEqDis) {
1515 : const Vector shift=ddist_dcreference-csum;
1516 : #pragma omp simd
1517 : for(unsigned iat=0; iat<n; iat++) {
1518 855 : const double a=align[iat];
1519 855 : derivatives[iat][0]=prefactor*(derivatives[iat][0]+shift[0]*a);
1520 855 : derivatives[iat][1]=prefactor*(derivatives[iat][1]+shift[1]*a);
1521 855 : derivatives[iat][2]=prefactor*(derivatives[iat][2]+shift[2]*a);
1522 : }
1523 : }
1524 :
1525 1 : return derivatives;
1526 0 : }
1527 :
1528 : /// this version does not calculate the derivative of rotation matrix as needed for SOMA
1529 0 : std::vector<Vector> RMSDCoreData::getDDistanceDReferenceSOMA() {
1530 0 : std::vector<Vector> derivatives;
1531 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1532 : Vector ddist_dcreference;
1533 0 : derivatives.resize(n);
1534 : double prefactor=1.0;
1535 0 : if(!distanceIsMSD) {
1536 0 : prefactor*=0.5/dist;
1537 : }
1538 : Vector tmp1;
1539 :
1540 0 : if(!hasDistance) {
1541 0 : plumed_merror("getDDistanceDReference needs to calculate the distance via getDistance first !");
1542 : }
1543 0 : if(!isInitialized) {
1544 0 : plumed_merror("getDDistanceDReference to initialize the coreData first!");
1545 : }
1546 : // get the transpose rotation
1547 0 : Tensor t_rotation=rotation.transpose();
1548 :
1549 : // third expensive loop: derivatives
1550 0 : for(unsigned iat=0; iat<n; iat++) {
1551 0 : if(alEqDis) {
1552 : // there is no need for derivatives of rotation and shift here as it is by construction zero
1553 : // (similar to Hellman-Feynman forces)
1554 : //TODO: check this derivative down here
1555 0 : derivatives[iat]= -2*prefactor*align[iat]*matmul(t_rotation,d[iat]);
1556 : } else {
1557 : // these are the derivatives assuming the roto-translation as frozen
1558 0 : tmp1=2*displace[iat]*matmul(t_rotation,d[iat]);
1559 0 : derivatives[iat]= -tmp1;
1560 : // derivative of cpositions
1561 : ddist_dcreference+=tmp1;
1562 : }
1563 : }
1564 :
1565 0 : if(!alEqDis)
1566 0 : for(unsigned iat=0; iat<n; iat++) {
1567 0 : derivatives[iat]=prefactor*(derivatives[iat]+ddist_dcreference*align[iat]);
1568 : }
1569 :
1570 0 : return derivatives;
1571 0 : }
1572 :
1573 :
1574 :
1575 : /*
1576 : This below is the derivative of the rotation matrix that aligns the reference onto the positions
1577 : respect to positions
1578 : note that the this transformation overlap the reference onto position
1579 : if inverseTransform=true then aligns the positions onto reference
1580 : */
1581 148098 : Matrix<std::vector<Vector> > RMSDCoreData::getDRotationDPositions( const bool inverseTransform ) {
1582 148098 : const unsigned n=static_cast<unsigned int>(reference.size());
1583 148098 : if(!isInitialized) {
1584 0 : plumed_merror("getDRotationDPosition to initialize the coreData first!");
1585 : }
1586 148098 : Matrix<std::vector<Vector> > DRotDPos=Matrix<std::vector<Vector> >(3,3);
1587 : // remember drotation_drr01 is Tensor drotation_drr01[3][3]
1588 : // (3x3 rot) (3x3 components of rr01)
1589 148098 : std::vector<Vector> v(n);
1590 : Vector csum;
1591 : // these below could probably be calculated in the main routine
1592 : Vector cp;
1593 : cp.zero();
1594 148098 : if(!cpositions_is_removed) {
1595 : cp=cpositions;
1596 : }
1597 : Vector cr;
1598 : cr.zero();
1599 148098 : if(!creference_is_removed) {
1600 60 : cr=creference;
1601 : }
1602 2557987 : for(unsigned iat=0; iat<n; iat++) {
1603 2409889 : csum+=(reference[iat]-cr)*align[iat];
1604 : }
1605 2557987 : for(unsigned iat=0; iat<n; iat++) {
1606 2409889 : v[iat]=(reference[iat]-cr-csum)*align[iat];
1607 : }
1608 592392 : for(unsigned a=0; a<3; a++) {
1609 1777176 : for(unsigned b=0; b<3; b++) {
1610 1332882 : if(inverseTransform) {
1611 1332882 : DRotDPos[b][a].resize(n);
1612 23021883 : for(unsigned iat=0; iat<n; iat++) {
1613 21689001 : DRotDPos[b][a][iat]=matmul(drotation_drr01[a][b],v[iat]);
1614 : }
1615 : } else {
1616 0 : DRotDPos[a][b].resize(n);
1617 0 : for(unsigned iat=0; iat<n; iat++) {
1618 0 : DRotDPos[a][b][iat]=matmul(drotation_drr01[a][b],v[iat]);
1619 : }
1620 : }
1621 : }
1622 : }
1623 148098 : return DRotDPos;
1624 148098 : }
1625 :
1626 : /*
1627 : This below is the derivative of the rotation matrix that aligns the reference onto the positions
1628 : respect to reference
1629 : note that the this transformation overlap the reference onto position
1630 : if inverseTransform=true then aligns the positions onto reference
1631 : */
1632 0 : Matrix<std::vector<Vector> > RMSDCoreData::getDRotationDReference( const bool inverseTransform ) {
1633 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1634 0 : if(!isInitialized) {
1635 0 : plumed_merror("getDRotationDPositions to initialize the coreData first!");
1636 : }
1637 0 : Matrix<std::vector<Vector> > DRotDRef=Matrix<std::vector<Vector> >(3,3);
1638 : // remember drotation_drr01 is Tensor drotation_drr01[3][3]
1639 : // (3x3 rot) (3x3 components of rr01)
1640 0 : std::vector<Vector> v(n);
1641 : Vector csum;
1642 : // these below could probably be calculated in the main routine
1643 : Vector cp;
1644 : cp.zero();
1645 0 : if(!cpositions_is_removed) {
1646 0 : cp=cpositions;
1647 : }
1648 : Vector cr;
1649 : cr.zero();
1650 0 : if(!creference_is_removed) {
1651 : cr=creference;
1652 : }
1653 0 : for(unsigned iat=0; iat<n; iat++) {
1654 0 : csum+=(positions[iat]-cp)*align[iat];
1655 : }
1656 0 : for(unsigned iat=0; iat<n; iat++) {
1657 0 : v[iat]=(positions[iat]-cp-csum)*align[iat];
1658 : }
1659 :
1660 0 : for(unsigned a=0; a<3; a++) {
1661 0 : for(unsigned b=0; b<3; b++) {
1662 0 : Tensor t_drotation_drr01=drotation_drr01[a][b].transpose();
1663 0 : if(inverseTransform) {
1664 0 : DRotDRef[b][a].resize(n);
1665 0 : for(unsigned iat=0; iat<n; iat++) {
1666 0 : DRotDRef[b][a][iat]=matmul(t_drotation_drr01,v[iat]);
1667 : }
1668 : } else {
1669 0 : DRotDRef[a][b].resize(n);
1670 0 : for(unsigned iat=0; iat<n; iat++) {
1671 0 : DRotDRef[a][b][iat]=matmul(t_drotation_drr01,v[iat]);
1672 : }
1673 : }
1674 : }
1675 : }
1676 0 : return DRotDRef;
1677 0 : }
1678 :
1679 :
1680 0 : std::vector<Vector> RMSDCoreData::getAlignedReferenceToPositions() {
1681 0 : std::vector<Vector> alignedref;
1682 0 : const unsigned n=static_cast<unsigned int>(reference.size());
1683 0 : alignedref.resize(n);
1684 0 : if(!isInitialized) {
1685 0 : plumed_merror("getAlignedReferenceToPostions needs to initialize the coreData first!");
1686 : }
1687 : // avoid to calculate matrix element but use the sum of what you have
1688 : Vector cp;
1689 : cp.zero();
1690 0 : if(!cpositions_is_removed) {
1691 0 : cp=cpositions;
1692 : }
1693 0 : for(unsigned iat=0; iat<n; iat++) {
1694 0 : alignedref[iat]=-d[iat]+positions[iat]-cp;
1695 : }
1696 0 : return alignedref;
1697 0 : }
1698 147888 : std::vector<Vector> RMSDCoreData::getAlignedPositionsToReference() {
1699 147888 : std::vector<Vector> alignedpos;
1700 147888 : if(!isInitialized) {
1701 0 : plumed_merror("getAlignedPostionsToReference needs to initialize the coreData first!");
1702 : }
1703 147888 : const unsigned n=static_cast<unsigned int>(positions.size());
1704 147888 : alignedpos.resize(n);
1705 : Vector cp;
1706 : cp.zero();
1707 147888 : if(!cpositions_is_removed) {
1708 147888 : cp=cpositions;
1709 : }
1710 : // avoid to calculate matrix element but use the sum of what you have
1711 2071777 : for(unsigned iat=0; iat<n; iat++) {
1712 1923889 : alignedpos[iat]=matmul(rotation.transpose(),positions[iat]-cp);
1713 : }
1714 147888 : return alignedpos;
1715 0 : }
1716 :
1717 :
1718 148098 : std::vector<Vector> RMSDCoreData::getCenteredPositions() {
1719 148098 : std::vector<Vector> centeredpos;
1720 148098 : const unsigned n=static_cast<unsigned int>(reference.size());
1721 148098 : centeredpos.resize(n);
1722 148098 : if(!isInitialized) {
1723 0 : plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1724 : }
1725 : // avoid to calculate matrix element but use the sum of what you have
1726 2557987 : for(unsigned iat=0; iat<n; iat++) {
1727 2409889 : centeredpos[iat]=positions[iat]-cpositions;
1728 : }
1729 148098 : return centeredpos;
1730 0 : }
1731 :
1732 147888 : std::vector<Vector> RMSDCoreData::getCenteredReference() {
1733 147888 : std::vector<Vector> centeredref;
1734 147888 : const unsigned n=static_cast<unsigned int>(reference.size());
1735 147888 : centeredref.resize(n);
1736 147888 : if(!isInitialized) {
1737 0 : plumed_merror("getCenteredReference needs to initialize the coreData first!");
1738 : }
1739 : // avoid to calculate matrix element but use the sum of what you have
1740 : Vector cr;
1741 : cr.zero();
1742 147888 : if(!creference_is_removed) {
1743 0 : cr=creference;
1744 : }
1745 2071777 : for(unsigned iat=0; iat<n; iat++) {
1746 1923889 : centeredref[iat]=reference[iat]-cr;
1747 : }
1748 147888 : return centeredref;
1749 0 : }
1750 :
1751 :
1752 210 : Vector RMSDCoreData::getPositionsCenter() {
1753 210 : if(!isInitialized) {
1754 0 : plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1755 : }
1756 210 : return cpositions;
1757 : }
1758 :
1759 0 : Vector RMSDCoreData::getReferenceCenter() {
1760 0 : if(!isInitialized) {
1761 0 : plumed_merror("getCenteredPositions needs to initialize the coreData first!");
1762 : }
1763 0 : return creference;
1764 : }
1765 :
1766 1764 : Tensor RMSDCoreData::getRotationMatrixReferenceToPositions() {
1767 1764 : if(!isInitialized) {
1768 0 : plumed_merror("getRotationMatrixReferenceToPositions needs to initialize the coreData first!");
1769 : }
1770 1764 : return rotation;
1771 : }
1772 :
1773 148098 : Tensor RMSDCoreData::getRotationMatrixPositionsToReference() {
1774 148098 : if(!isInitialized) {
1775 0 : plumed_merror("getRotationMatrixReferenceToPositions needs to initialize the coreData first!");
1776 : }
1777 148098 : return rotation.transpose();
1778 : }
1779 :
1780 1092 : const std::array<std::array<Tensor,3>,3>& RMSDCoreData::getDRotationDRr01() const {
1781 1092 : if(!isInitialized) {
1782 0 : plumed_merror("getDRotationDRr01 needs to initialize the coreData first!");
1783 : }
1784 1092 : return drotation_drr01;
1785 : }
1786 :
1787 : template double RMSD::optimalAlignment<true,true>(
1788 : const std::vector<double>& align,
1789 : const std::vector<double>& displaceFrom,
1790 : const View<Vector> positions,
1791 : const std::vector<Vector>& reference,
1792 : std::vector<Vector>& derivatives, bool squared)const;
1793 :
1794 : template double RMSD::optimalAlignment<true,false>(
1795 : const std::vector<double>& align,
1796 : const std::vector<double>& displaceFrom,
1797 : const View<Vector> positions,
1798 : const std::vector<Vector>& reference,
1799 : std::vector<Vector>& derivatives, bool squared)const;
1800 :
1801 : template double RMSD::optimalAlignment<false,true>(
1802 : const std::vector<double>& align,
1803 : const std::vector<double>& displaceFrom,
1804 : const View<Vector> positions,
1805 : const std::vector<Vector>& reference,
1806 : std::vector<Vector>& derivatives, bool squared)const;
1807 :
1808 : template double RMSD::optimalAlignment<false,false>(
1809 : const std::vector<double>& align,
1810 : const std::vector<double>& displaceFrom,
1811 : const View<Vector> positions,
1812 : const std::vector<Vector>& reference,
1813 : std::vector<Vector>& derivatives, bool squared)const;
1814 : }
|
