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