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 : #ifndef __PLUMED_tools_RMSD_h
23 : #define __PLUMED_tools_RMSD_h
24 :
25 : #include "Vector.h"
26 : #include "Matrix.h"
27 : #include "Tensor.h"
28 : #include "View.h"
29 : #include <vector>
30 : #include <string>
31 : #include <array>
32 :
33 : namespace PLMD {
34 :
35 : class Log;
36 : class PDB;
37 :
38 : /** \ingroup TOOLBOX
39 : A class that implements RMSD calculations
40 : This is a class that implements the various infrastructure to calculate the
41 : RMSD or MSD respect a given frame. It can be done through an optimal alignment scheme
42 : as Kearsley or, more simply, by resetting the center of mass.
43 : This is the class that decides this. A very simple use is
44 : \verbatim
45 : #include "tools/PDB.h"
46 : #include "tools/RMSD.h"
47 : #include "tools/Vector.h"
48 : using namespace PLMD;
49 : RMSD rmsd;
50 : PDB pdb;
51 : // get the pdb (see PDB documentation)
52 : pdb.read("file.pdb",true,1.0);
53 : string type;
54 : type.assign("OPTIMAL");
55 : // set the reference and the type
56 : rmsd.set(pdb,type);
57 : // this calculates the rmsd and the derivatives
58 : vector<Vector> derivs;
59 : double val;
60 : val=rmsd.calculate(getPositions(),derivs,true);
61 : \endverbatim
62 :
63 : **/
64 :
65 : class RMSD {
66 : enum AlignmentMethod {SIMPLE, OPTIMAL, OPTIMAL_FAST};
67 : AlignmentMethod alignmentMethod=SIMPLE;
68 : // Reference coordinates
69 : std::vector<Vector> reference;
70 : // Weights for alignment
71 : std::vector<double> align;
72 : // Weights for deviation
73 : std::vector<double> displace;
74 : // Center for reference and flag for its calculation
75 : Vector reference_center;
76 : bool reference_center_is_calculated=false;
77 : bool reference_center_is_removed=false;
78 : // Center for running position (not used in principle but here to reflect reference/positio symmetry
79 : Vector positions_center;
80 : bool positions_center_is_calculated=false;
81 : bool positions_center_is_removed=false;
82 : // calculates the center from the position provided
83 28872 : Vector calculateCenter(const std::vector<Vector>&p,const std::vector<double> &w) {
84 28872 : plumed_massert(p.size()==w.size(),"mismatch in dimension of position/align arrays while calculating the center");
85 : unsigned n;
86 28872 : n=p.size();
87 28872 : Vector c;
88 : c.zero();
89 430372 : for(unsigned i=0; i<n; i++) {
90 401500 : c+=p[i]*w[i];
91 : }
92 28872 : return c;
93 : };
94 : // removes the center for the position provided
95 57739 : void removeCenter(std::vector<Vector>&p, const Vector &c) {
96 : unsigned n;
97 57739 : n=p.size();
98 860714 : for(unsigned i=0; i<n; i++) {
99 802975 : p[i]-=c;
100 : }
101 57739 : };
102 : // add center
103 28872 : void addCenter(std::vector<Vector>&p, const Vector &c) {
104 : Vector cc=c*-1.;
105 28872 : removeCenter(p,cc);
106 28872 : };
107 :
108 : public:
109 : /// Constructor
110 : RMSD();
111 : /// clear the structure
112 : void clear();
113 : /// set reference, align and displace from input pdb structure: evtl remove com from the initial structure and normalize the input weights from the pdb
114 : void set(const PDB&,
115 : const std::string & mytype,
116 : bool remove_center=true,
117 : bool normalize_weights=true);
118 : /// set align displace reference and type from input vectors
119 : void set(const std::vector<double>& alignTo,
120 : const std::vector<double>& displaceFrom,
121 : const std::vector<Vector>& ref,
122 : const std::string & mytype,
123 : bool remove_center=true,
124 : bool normalize_weights=true );
125 : /// set the type of alignment we are doing
126 : void setType(const std::string & mytype);
127 : /// set reference coordinates, remove the com by using uniform weights
128 : void setReference(const std::vector<Vector>& ref);
129 : const std::vector<Vector>& getReference() const ;
130 : /// set weights and remove the center from reference with normalized weights. If the com has been removed, it resets to the new value
131 : void setAlign(const std::vector<double>& alignTo,
132 : bool normalize_weights=true,
133 : bool remove_center=true);
134 : std::vector<double> getAlign();
135 : /// set align
136 : void setDisplace(const std::vector<double>& displaceFrom, bool normalize_weights=true);
137 : std::vector<double> getDisplace();
138 : ///
139 : std::string getMethod() const ;
140 : /// workhorses
141 : double simpleAlignment(const std::vector<double>& alignTo,
142 : const std::vector<double>& displaceFrom,
143 : const View<Vector> positions,
144 : const std::vector<Vector>& ref,
145 : std::vector<Vector>& derivatives,
146 : std::vector<Vector>& displacement,
147 : bool squared=false)const;
148 : template <bool safe,bool alEqDis>
149 : double optimalAlignment(const std::vector<double>& alignTo,
150 : const std::vector<double>& displaceFrom,
151 : const View<Vector> positions,
152 : const std::vector<Vector>& ref,
153 : std::vector<Vector>& DDistDPos, bool squared=false)const;
154 :
155 : template <bool safe, bool alEqDis>
156 : double optimalAlignmentWithCloseStructure(const std::vector<double>& alignTo,
157 : const std::vector<double>& displaceFrom,
158 : const View<Vector> positions,
159 : const std::vector<Vector>& ref,
160 : std::vector<Vector>& derivatives,
161 : const Tensor & rotationPosClose,
162 : const Tensor & rotationRefClose,
163 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01,
164 : bool squared=false)const;
165 :
166 : template <bool safe, bool alEqDis>
167 : double optimalAlignment_Rot_DRotDRr01(const std::vector<double>& alignTo,
168 : const std::vector<double>& displaceFrom,
169 : const View<Vector> positions,
170 : const std::vector<Vector>& ref,
171 : Tensor & Rotation,
172 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01,
173 : bool squared=false)const;
174 :
175 : template <bool safe, bool alEqDis>
176 : double optimalAlignment_Rot(const std::vector<double>& alignTo,
177 : const std::vector<double>& displaceFrom,
178 : const View<Vector> positions,
179 : const std::vector<Vector>& ref,
180 : std::vector<Vector>& derivatives,
181 : Tensor & Rotation,
182 : bool squared=false)const;
183 :
184 : template <bool safe,bool alEqDis>
185 : double optimalAlignment_DDistDRef(const std::vector<double>& alignTo,
186 : const std::vector<double>& displaceFrom,
187 : const View<Vector> positions,
188 : const std::vector<Vector>& ref,
189 : std::vector<Vector>& DDistDPos,
190 : std::vector<Vector>& DDistDRef,
191 : bool squared=false) const;
192 :
193 : template <bool safe,bool alEqDis>
194 : double optimalAlignment_SOMA(const std::vector<double>& alignTo,
195 : const std::vector<double>& displaceFrom,
196 : const View<Vector> positions,
197 : const std::vector<Vector>& ref,
198 : std::vector<Vector>& DDistDPos,
199 : std::vector<Vector>& DDistDRef,
200 : bool squared=false) const;
201 :
202 : template <bool safe,bool alEqDis>
203 : double optimalAlignment_DDistDRef_Rot_DRotDPos(const std::vector<double>& alignTo,
204 : const std::vector<double>& displaceFrom,
205 : const View<Vector> positions,
206 : const std::vector<Vector>& ref,
207 : std::vector<Vector>& DDistDPos,
208 : std::vector<Vector>& DDistDRef,
209 : Tensor & Rotation,
210 : Matrix<std::vector<Vector> >&DRotDPos,
211 : bool squared=false) const;
212 :
213 : template <bool safe,bool alEqDis>
214 : double optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef(const std::vector<double>& alignTo,
215 : const std::vector<double>& displaceFrom,
216 : const View<Vector> positions,
217 : const std::vector<Vector>& ref,
218 : std::vector<Vector>& DDistDPos,
219 : std::vector<Vector>& DDistDRef,
220 : Tensor & Rotation,
221 : Matrix<std::vector<Vector> >&DRotDPos,
222 : Matrix<std::vector<Vector> >&DRotDRef,
223 : bool squared=false) const;
224 :
225 : template <bool safe,bool alEqDis>
226 : double optimalAlignment_PCA(const std::vector<double>& alignTo,
227 : const std::vector<double>& displaceFrom,
228 : const View<Vector> positions,
229 : const std::vector<Vector>& ref,
230 : std::vector<Vector>& alignedpositions,
231 : std::vector<Vector>& centeredpositions,
232 : std::vector<Vector>& centeredreference,
233 : Tensor & Rotation,
234 : std::vector<Vector>& DDistDPos,
235 : Matrix<std::vector<Vector> >& DRotDPos,
236 : bool squared=false) const ;
237 :
238 : template <bool safe,bool alEqDis>
239 : double optimalAlignment_Fit(const std::vector<double>& alignTo,
240 : const std::vector<double>& displaceFrom,
241 : const View<Vector> positions,
242 : const std::vector<Vector>& ref,
243 : Tensor & Rotation,
244 : Matrix<std::vector<Vector> >& DRotDPos,
245 : std::vector<Vector>& centeredpositions,
246 : Vector & center_positions,
247 : bool squared=false);
248 :
249 :
250 : /// Compute rmsd: note that this is an intermediate layer which is kept in order to evtl expand with more alignment types/user options to be called while keeping the workhorses separated
251 : //#pragma acc routine
252 : double calculate(const View<Vector> positions,
253 : std::vector<Vector>&derivatives,
254 : bool squared=false)const;
255 : double calculate(const std::vector<Vector>& positions,
256 : std::vector<Vector>&derivatives,
257 : bool squared=false)const;
258 : /// Other convenience methods:
259 : /// calculate the derivative of distance respect to position(DDistDPos) and reference (DDistDPos)
260 : double calc_DDistDRef( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, std::vector<Vector>& DDistDRef, const bool squared=false );
261 : /// calculate the derivative for SOMA (i.e. derivative respect to reference frame without the matrix derivative)
262 : double calc_SOMA( const std::vector<Vector>& positions,
263 : std::vector<Vector>&DDistDPos,
264 : std::vector<Vector>& DDistDRef,
265 : const bool squared=false);
266 : ///
267 : double calc_DDistDRef_Rot_DRotDPos( const std::vector<Vector>& positions,
268 : std::vector<Vector>&DDistDPos,
269 : std::vector<Vector>& DDistDRef,
270 : Tensor & Rotation,Matrix<std::vector<Vector> >&DRotDPos,
271 : const bool squared=false);
272 : double calc_DDistDRef_Rot_DRotDPos_DRotDRef( const std::vector<Vector>& positions,
273 : std::vector<Vector>&DDistDPos,
274 : std::vector<Vector>& DDistDRef,
275 : Tensor & Rotation,Matrix<std::vector<Vector> >&DRotDPos,Matrix<std::vector<Vector> > &DRotDRef,
276 : const bool squared=false);
277 : /// convenience method to retrieve all the bits and pieces for PCA
278 : double calc_PCAelements( const std::vector<Vector>& pos,
279 : std::vector<Vector>&DDistDPos,
280 : Tensor & Rotation,
281 : Matrix<std::vector<Vector> >& DRotDPos,std::vector<Vector> & alignedpositions,
282 : std::vector<Vector>& centeredpositions,
283 : std::vector<Vector>¢eredreference,
284 : bool squared=false) const ;
285 : /// convenience method to retrieve all the bits and pieces needed by FitToTemplate
286 : double calc_FitElements( const std::vector<Vector>& pos,
287 : Tensor & Rotation,
288 : Matrix<std::vector<Vector> >& DRotDPos,std::vector<Vector> & centeredpositions,Vector & center_positions,
289 : bool squared=false );
290 : ///calculate rotation matrix, derivative of rotation matrix w.r.t. positions, derivative of rotation matrix w.r.t. rr01
291 : double calc_Rot_DRotDRr01( const std::vector<Vector>& positions,
292 : Tensor & Rotation,
293 : std::array<std::array<Tensor,3>,3>& DRotDRr01,
294 : const bool squared=false );
295 : ///calculate rotation matrix, derivative of rotation matrix w.r.t. positions
296 : double calc_Rot( const std::vector<Vector>& positions,
297 : std::vector<Vector>&DDistDPos,
298 : Tensor & Rotation,
299 : const bool squared=false );
300 : ///calculate with close structure, i.e. approximate the RMSD without expensive computation of rotation matrix by reusing saved rotation matrices from previous iterations
301 : double calculateWithCloseStructure( const std::vector<Vector>& positions,
302 : std::vector<Vector>&DDistDPos,
303 : const Tensor & rotationPosClose,
304 : const Tensor & rotationRefClose,
305 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01,
306 : const bool squared=false );
307 : /// static convenience method to get the matrix i,a from drotdpos (which is a bit tricky)
308 : static Tensor getMatrixFromDRot(const Matrix< std::vector<Vector> >&drotdpos, const unsigned i, const unsigned a) {
309 : Tensor t;
310 : t[0][0]=drotdpos(0,0)[i][a];
311 : t[0][1]=drotdpos(0,1)[i][a];
312 : t[0][2]=drotdpos(0,2)[i][a];
313 : t[1][0]=drotdpos(1,0)[i][a];
314 : t[1][1]=drotdpos(1,1)[i][a];
315 : t[1][2]=drotdpos(1,2)[i][a];
316 : t[2][0]=drotdpos(2,0)[i][a];
317 : t[2][1]=drotdpos(2,1)[i][a];
318 : t[2][2]=drotdpos(2,2)[i][a];
319 : return t;
320 : };
321 : };
322 :
323 : /// this is a class which is needed to share information across the various non-threadsafe routines
324 : /// so that the public function of rmsd are threadsafe while the inner core can safely share information
325 195055 : class RMSDCoreData {
326 : private:
327 : //the parameters that are common to the two constructors
328 : //are inititialized here at declaration,
329 : //to underline differences in setups
330 : bool alEqDis=false;
331 : bool distanceIsMSD=false; // default is RMSD but can deliver the MSD
332 : bool hasDistance=false; // distance is already calculated
333 : bool isInitialized=false;
334 : bool safe=false;
335 :
336 : // this need to be copied and they are small, should not affect the performance
337 : Vector creference;
338 : bool creference_is_calculated;
339 : bool creference_is_removed;
340 : Vector cpositions;
341 : bool cpositions_is_calculated;
342 : bool cpositions_is_removed;
343 :
344 : // use reference assignment to speed up instead of copying
345 : const View<Vector> positions;
346 : const std::vector<Vector>&reference;
347 : const std::vector<double>&align;
348 : const std::vector<double>&displace;
349 :
350 : // the needed stuff for distance and more (one could use eigenvecs components and eigenvals for some reason)
351 : double dist=0.0;
352 : Vector4d eigenvals;
353 : Tensor4d eigenvecs;
354 : double rr00=0.0; // sum of positions squared (needed for dist calc)
355 : double rr11=0.0; // sum of reference squared (needed for dist calc)
356 : Tensor rotation; // rotation derived from the eigenvector having the smallest eigenvalue
357 : std::array<std::array<Tensor,3>,3> drotation_drr01; // derivative of the rotation only available when align!=displace
358 : Tensor ddist_drr01;
359 : Tensor ddist_drotation;
360 : std::vector<Vector> d; // difference of components
361 : public:
362 : /// the constructor (note: only references are passed, therefore is rather fast)
363 : /// note: this aligns the reference onto the positions
364 : ///
365 : /// this method assumes that the centers are already calculated and subtracted
366 : RMSDCoreData(const std::vector<double>&a,
367 : const std::vector<double>&disp,
368 : const View<Vector> p,
369 : const std::vector<Vector>&r,
370 : const Vector &cp,
371 : const Vector &cr )
372 : : creference(cr),
373 : creference_is_calculated(true),
374 : creference_is_removed(true),
375 : cpositions(cp),
376 : cpositions_is_calculated(true),
377 : cpositions_is_removed(true),
378 : positions(p),
379 : reference(r),
380 : align(a),
381 : displace(disp) {}
382 :
383 : // this constructor does not assume that the positions and reference have the center subtracted
384 195055 : RMSDCoreData(const std::vector<double>&a,
385 : const std::vector<double>&disp,
386 : const View<Vector> p,
387 : const std::vector<Vector>&r)
388 195055 : : creference(0.0,0.0,0.0),
389 195055 : creference_is_calculated(false),
390 195055 : creference_is_removed(false),
391 : cpositions(0.0,0.0,0.0),
392 195055 : cpositions_is_calculated(false),
393 195055 : cpositions_is_removed(false),
394 195055 : positions(p),
395 195055 : reference(r),
396 195055 : align(a),
397 195055 : displace(disp) {}
398 :
399 : // set the center on the fly without subtracting
400 195055 : void calcPositionsCenter() {
401 195055 : plumed_massert(!cpositions_is_calculated,"the center was already calculated");
402 : cpositions.zero();
403 3216227 : for(unsigned i=0; i<positions.size(); i++) {
404 3021172 : cpositions+=positions[i]*align[i];
405 : }
406 195055 : cpositions_is_calculated=true;
407 195055 : }
408 :
409 0 : void calcReferenceCenter() {
410 0 : plumed_massert(!creference_is_calculated,"the center was already calculated");
411 : creference.zero();
412 0 : for(unsigned i=0; i<reference.size(); i++) {
413 0 : creference+=reference[i]*align[i];
414 : }
415 0 : creference_is_calculated=true;
416 0 : }
417 : // assume the center is given externally
418 : // cppcheck-suppress passedByValue
419 0 : void setPositionsCenter(Vector v) {
420 0 : plumed_massert(!cpositions_is_calculated,"You are setting the center two times!");
421 0 : cpositions=v;
422 0 : cpositions_is_calculated=true;
423 0 : }
424 : // cppcheck-suppress passedByValue
425 195055 : void setReferenceCenter(Vector v) {
426 195055 : plumed_massert(!creference_is_calculated,"You are setting the center two times!");
427 195055 : creference=v;
428 195055 : creference_is_calculated=true;
429 195055 : }
430 : // the center is already removed
431 : void setPositionsCenterIsRemoved(bool t) {
432 195055 : cpositions_is_removed=t;
433 : }
434 : void setReferenceCenterIsRemoved(bool t) {
435 195055 : creference_is_removed=t;
436 : }
437 : bool getPositionsCenterIsRemoved() {
438 : return cpositions_is_removed;
439 : }
440 : bool getReferenceCenterIsRemoved() {
441 : return creference_is_removed;
442 : }
443 : // does the core calc : first thing to call after the constructor:
444 : void doCoreCalc(bool setSafe,bool setAlEqDis);
445 : // do calculation with close structure data structures
446 : void doCoreCalcWithCloseStructure(bool setSafe,
447 : bool setAlEqDis,
448 : const Tensor & rotationPosClose,
449 : const Tensor & rotationRefClose,
450 : std::array<std::array<Tensor,3>,3>& drotationPosCloseDrr01);
451 : // retrieve the distance if required after doCoreCalc
452 : double getDistance(bool squared);
453 : // retrieve the derivative of the distance respect to the position
454 : std::vector<Vector> getDDistanceDPositions();
455 : // retrieve the derivative of the distance respect to the reference
456 : std::vector<Vector> getDDistanceDReference();
457 : // specific version for SOMA calculation (i.e. does not need derivative respect to rotation matrix)
458 : std::vector<Vector> getDDistanceDReferenceSOMA();
459 : // get aligned reference onto position
460 : std::vector<Vector> getAlignedReferenceToPositions();
461 : // get aligned position onto reference
462 : std::vector<Vector> getAlignedPositionsToReference();
463 : // get centered positions
464 : std::vector<Vector> getCenteredPositions();
465 : // get centered reference
466 : std::vector<Vector> getCenteredReference();
467 : // get center of positions
468 : Vector getPositionsCenter();
469 : // get center of reference
470 : Vector getReferenceCenter();
471 : // get rotation matrix (reference ->positions)
472 : Tensor getRotationMatrixReferenceToPositions();
473 : // get rotation matrix (positions -> reference)
474 : Tensor getRotationMatrixPositionsToReference();
475 : // get the derivative of the rotation matrix respect to positions
476 : // note that the this transformation overlap the reference onto position
477 : // if inverseTransform=true then aligns the positions onto reference
478 : Matrix<std::vector<Vector> > getDRotationDPositions( bool inverseTransform=false );
479 : // get the derivative of the rotation matrix respect to reference
480 : // note that the this transformation overlap the reference onto position
481 : // if inverseTransform=true then aligns the positions onto reference
482 : Matrix<std::vector<Vector> > getDRotationDReference(bool inverseTransform=false );
483 : const std::array<std::array<Tensor,3>,3>& getDRotationDRr01() const;
484 : };
485 :
486 : }
487 :
488 : #endif
489 :
|
