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