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