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