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 : {
66 : enum AlignmentMethod {SIMPLE, OPTIMAL, OPTIMAL_FAST};
67 : AlignmentMethod alignmentMethod;
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;
77 : bool reference_center_is_removed;
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;
81 : bool positions_center_is_removed;
82 : // calculates the center from the position provided
83 1575 : Vector calculateCenter(const std::vector<Vector> &p,const std::vector<double> &w) {
84 1575 : plumed_massert(p.size()==w.size(),"mismatch in dimension of position/align arrays while calculating the center");
85 1575 : unsigned n; n=p.size();
86 1575 : Vector c; c.zero();
87 40933 : for(unsigned i=0; i<n; i++)c+=p[i]*w[i];
88 1575 : return c;
89 : };
90 : // removes the center for the position provided
91 3145 : void removeCenter(std::vector<Vector> &p, const Vector &c) {
92 3145 : unsigned n; n=p.size();
93 81836 : for(unsigned i=0; i<n; i++)p[i]-=c;
94 3145 : };
95 : // add center
96 1575 : void addCenter(std::vector<Vector> &p, const Vector &c) {Vector cc=c*-1.; removeCenter(p,cc);};
97 :
98 : public:
99 : /// Constructor
100 : RMSD();
101 : /// clear the structure
102 : void clear();
103 : /// set reference, align and displace from input pdb structure: evtl remove com from the initial structure and normalize the input weights from the pdb
104 : void set(const PDB&,const std::string & mytype, bool remove_center=true, bool normalize_weights=true);
105 : /// set align displace reference and type from input vectors
106 : 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 );
107 : /// set the type of alignment we are doing
108 : void setType(const std::string & mytype);
109 : /// set reference coordinates, remove the com by using uniform weights
110 : void setReference(const std::vector<Vector> & reference);
111 : std::vector<Vector> getReference();
112 : /// set weights and remove the center from reference with normalized weights. If the com has been removed, it resets to the new value
113 : void setAlign(const std::vector<double> & align, bool normalize_weights=true, bool remove_center=true);
114 : std::vector<double> getAlign();
115 : /// set align
116 : void setDisplace(const std::vector<double> & displace, bool normalize_weights=true);
117 : std::vector<double> getDisplace();
118 : ///
119 : std::string getMethod();
120 : /// workhorses
121 : double simpleAlignment(const std::vector<double> & align,
122 : const std::vector<double> & displace,
123 : const std::vector<Vector> & positions,
124 : const std::vector<Vector> & reference,
125 : std::vector<Vector> & derivatives,
126 : std::vector<Vector> & displacement,
127 : bool squared=false)const;
128 : template <bool safe,bool alEqDis>
129 : double optimalAlignment(const std::vector<double> & align,
130 : const std::vector<double> & displace,
131 : const std::vector<Vector> & positions,
132 : const std::vector<Vector> & reference,
133 : std::vector<Vector> & DDistDPos, bool squared=false)const;
134 :
135 : template <bool safe, bool alEqDis>
136 : double optimalAlignmentWithCloseStructure(const std::vector<double> & align,
137 : const std::vector<double> & displace,
138 : const std::vector<Vector> & positions,
139 : const std::vector<Vector> & reference,
140 : std::vector<Vector> & derivatives,
141 : const Tensor & rotationPosClose,
142 : const Tensor & rotationRefClose,
143 : std::array<std::array<Tensor,3>,3> & drotationPosCloseDrr01,
144 : bool squared=false)const;
145 :
146 : template <bool safe, bool alEqDis>
147 : double optimalAlignment_Rot_DRotDRr01(const std::vector<double> & align,
148 : const std::vector<double> & displace,
149 : const std::vector<Vector> & positions,
150 : const std::vector<Vector> & reference,
151 : Tensor & Rotation,
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(const std::vector<double> & align,
157 : const std::vector<double> & displace,
158 : const std::vector<Vector> & positions,
159 : const std::vector<Vector> & reference,
160 : std::vector<Vector> & derivatives,
161 : Tensor & Rotation,
162 : bool squared=false)const;
163 :
164 : template <bool safe,bool alEqDis>
165 : double optimalAlignment_DDistDRef(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> & DDistDPos,
170 : std::vector<Vector> & DDistDRef,
171 : bool squared=false) const;
172 :
173 : template <bool safe,bool alEqDis>
174 : double optimalAlignment_SOMA(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_DDistDRef_Rot_DRotDPos(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 : Tensor & Rotation,
190 : Matrix<std::vector<Vector> > &DRotDPos,
191 : bool squared=false) const;
192 :
193 : template <bool safe,bool alEqDis>
194 : double optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef(const std::vector<double> & align,
195 : const std::vector<double> & displace,
196 : const std::vector<Vector> & positions,
197 : const std::vector<Vector> & reference,
198 : std::vector<Vector> & DDistDPos,
199 : std::vector<Vector> & DDistDRef,
200 : Tensor & Rotation,
201 : Matrix<std::vector<Vector> > &DRotDPos,
202 : Matrix<std::vector<Vector> > &DRotDRef,
203 : bool squared=false) const;
204 :
205 : template <bool safe,bool alEqDis>
206 : double optimalAlignment_PCA(const std::vector<double> & align,
207 : const std::vector<double> & displace,
208 : const std::vector<Vector> & positions,
209 : const std::vector<Vector> & reference,
210 : std::vector<Vector> & alignedpositions,
211 : std::vector<Vector> & centeredpositions,
212 : std::vector<Vector> & centeredreference,
213 : Tensor & Rotation,
214 : std::vector<Vector> & DDistDPos,
215 : Matrix<std::vector<Vector> > & DRotDPos,
216 : bool squared=false) const ;
217 :
218 : template <bool safe,bool alEqDis>
219 : double optimalAlignment_Fit(const std::vector<double> & align,
220 : const std::vector<double> & displace,
221 : const std::vector<Vector> & positions,
222 : const std::vector<Vector> & reference,
223 : Tensor & Rotation,
224 : Matrix<std::vector<Vector> > & DRotDPos,
225 : std::vector<Vector> & centeredpositions,
226 : Vector & center_positions,
227 : bool squared=false);
228 :
229 :
230 : /// 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
231 : double calculate(const std::vector<Vector> & positions,std::vector<Vector> &derivatives, bool squared=false)const;
232 : /// Other convenience methods:
233 : /// calculate the derivative of distance respect to position(DDistDPos) and reference (DDistDPos)
234 : double calc_DDistDRef( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, std::vector<Vector>& DDistDRef, const bool squared=false );
235 : /// calculate the derivative for SOMA (i.e. derivative respect to reference frame without the matrix derivative)
236 : double calc_SOMA( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, std::vector<Vector>& DDistDRef, const bool squared=false );
237 : ///
238 : 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 );
239 : 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 );
240 : /// convenience method to retrieve all the bits and pieces for PCA
241 : 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 ;
242 : /// convenience method to retrieve all the bits and pieces needed by FitToTemplate
243 : 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 );
244 : ///calculate rotation matrix, derivative of rotation matrix w.r.t. positions, derivative of rotation matrix w.r.t. rr01
245 : double calc_Rot_DRotDRr01( const std::vector<Vector>& positions, Tensor & Rotation, std::array<std::array<Tensor,3>,3> & DRotDRr01, const bool squared=false );
246 : ///calculate rotation matrix, derivative of rotation matrix w.r.t. positions
247 : double calc_Rot( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, Tensor & Rotation, const bool squared=false );
248 : ///calculate with close structure, i.e. approximate the RMSD without expensive computation of rotation matrix by reusing saved rotation matrices from previous iterations
249 : 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 );
250 : /// static convenience method to get the matrix i,a from drotdpos (which is a bit tricky)
251 900 : static Tensor getMatrixFromDRot(const Matrix< std::vector<Vector> > &drotdpos, const unsigned &i, const unsigned &a) {
252 900 : Tensor t;
253 900 : t[0][0]=drotdpos(0,0)[i][a]; t[0][1]=drotdpos(0,1)[i][a]; t[0][2]=drotdpos(0,2)[i][a];
254 900 : t[1][0]=drotdpos(1,0)[i][a]; t[1][1]=drotdpos(1,1)[i][a]; t[1][2]=drotdpos(1,2)[i][a];
255 900 : t[2][0]=drotdpos(2,0)[i][a]; t[2][1]=drotdpos(2,1)[i][a]; t[2][2]=drotdpos(2,2)[i][a];
256 900 : return t;
257 : };
258 : };
259 :
260 : /// this is a class which is needed to share information across the various non-threadsafe routines
261 : /// so that the public function of rmsd are threadsafe while the inner core can safely share information
262 52152 : class RMSDCoreData
263 : {
264 : private:
265 : bool alEqDis;
266 : bool distanceIsMSD; // default is RMSD but can deliver the MSD
267 : bool hasDistance; // distance is already calculated
268 : bool isInitialized;
269 : bool safe;
270 :
271 : // this need to be copied and they are small, should not affect the performance
272 : Vector creference;
273 : bool creference_is_calculated;
274 : bool creference_is_removed;
275 : Vector cpositions;
276 : bool cpositions_is_calculated;
277 : bool cpositions_is_removed;
278 : bool retrieve_only_rotation;
279 :
280 : // use reference assignment to speed up instead of copying
281 : const std::vector<Vector> &positions;
282 : const std::vector<Vector> &reference;
283 : const std::vector<double> &align;
284 : const std::vector<double> &displace;
285 :
286 : // the needed stuff for distance and more (one could use eigenvecs components and eigenvals for some reason)
287 : double dist;
288 : Vector4d eigenvals;
289 : Tensor4d eigenvecs;
290 : double rr00; // sum of positions squared (needed for dist calc)
291 : double rr11; // sum of reference squared (needed for dist calc)
292 : Tensor rotation; // rotation derived from the eigenvector having the smallest eigenvalue
293 : std::array<std::array<Tensor,3>,3> drotation_drr01; // derivative of the rotation only available when align!=displace
294 : Tensor ddist_drr01;
295 : Tensor ddist_drotation;
296 : std::vector<Vector> d; // difference of components
297 : public:
298 : /// the constructor (note: only references are passed, therefore is rather fast)
299 : /// note: this aligns the reference onto the positions
300 : ///
301 : /// this method assumes that the centers are already calculated and subtracted
302 : 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 ):
303 : alEqDis(false),distanceIsMSD(false),hasDistance(false),isInitialized(false),safe(false),
304 : creference(cr),creference_is_calculated(true),creference_is_removed(true),
305 : 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) {};
306 :
307 : // this constructor does not assume that the positions and reference have the center subtracted
308 52152 : RMSDCoreData(const std::vector<double> &a,const std::vector<double> &d,const std::vector<Vector> &p, const std::vector<Vector> &r):
309 52152 : alEqDis(false),distanceIsMSD(false),hasDistance(false),isInitialized(false),safe(false),
310 52152 : creference_is_calculated(false),creference_is_removed(false),
311 52152 : 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)
312 52152 : {cpositions.zero(); creference.zero();};
313 :
314 : // set the center on the fly without subtracting
315 52152 : void calcPositionsCenter() {
316 52152 : plumed_massert(!cpositions_is_calculated,"the center was already calculated");
317 1221600 : cpositions.zero(); for(unsigned i=0; i<positions.size(); i++) {cpositions+=positions[i]*align[i];} cpositions_is_calculated=true;
318 52152 : }
319 0 : void calcReferenceCenter() {
320 0 : plumed_massert(!creference_is_calculated,"the center was already calculated");
321 0 : creference.zero(); for(unsigned i=0; i<reference.size(); i++) {creference+=reference[i]*align[i];} creference_is_calculated=true;
322 0 : };
323 : // assume the center is given externally
324 0 : void setPositionsCenter(Vector v) {plumed_massert(!cpositions_is_calculated,"You are setting the center two times!"); cpositions=v; cpositions_is_calculated=true;};
325 52152 : void setReferenceCenter(Vector v) {plumed_massert(!creference_is_calculated,"You are setting the center two times!"); creference=v; creference_is_calculated=true;};
326 : // the center is already removed
327 52152 : void setPositionsCenterIsRemoved(bool t) {cpositions_is_removed=t;};
328 52152 : void setReferenceCenterIsRemoved(bool t) {creference_is_removed=t;};
329 : bool getPositionsCenterIsRemoved() {return cpositions_is_removed;};
330 : bool getReferenceCenterIsRemoved() {return creference_is_removed;};
331 : // does the core calc : first thing to call after the constructor:
332 : // only_rotation=true does not retrieve the derivatives, just retrieve the optimal rotation (the same calc cannot be exploit further)
333 : void doCoreCalc(bool safe,bool alEqDis, bool only_rotation=false);
334 : // do calculation with close structure data structures
335 : void doCoreCalcWithCloseStructure(bool safe,bool alEqDis, const Tensor & rotationPosClose, const Tensor & rotationRefClose, std::array<std::array<Tensor,3>,3> & drotationPosCloseDrr01);
336 : // retrieve the distance if required after doCoreCalc
337 : double getDistance(bool squared);
338 : // retrieve the derivative of the distance respect to the position
339 : std::vector<Vector> getDDistanceDPositions();
340 : // retrieve the derivative of the distance respect to the reference
341 : std::vector<Vector> getDDistanceDReference();
342 : // specific version for SOMA calculation (i.e. does not need derivative respect to rotation matrix)
343 : std::vector<Vector> getDDistanceDReferenceSOMA();
344 : // get aligned reference onto position
345 : std::vector<Vector> getAlignedReferenceToPositions();
346 : // get aligned position onto reference
347 : std::vector<Vector> getAlignedPositionsToReference();
348 : // get centered positions
349 : std::vector<Vector> getCenteredPositions();
350 : // get centered reference
351 : std::vector<Vector> getCenteredReference();
352 : // get center of positions
353 : Vector getPositionsCenter();
354 : // get center of reference
355 : Vector getReferenceCenter();
356 : // get rotation matrix (reference ->positions)
357 : Tensor getRotationMatrixReferenceToPositions();
358 : // get rotation matrix (positions -> reference)
359 : Tensor getRotationMatrixPositionsToReference();
360 : // get the derivative of the rotation matrix respect to positions
361 : // note that the this transformation overlap the reference onto position
362 : // if inverseTransform=true then aligns the positions onto reference
363 : Matrix<std::vector<Vector> > getDRotationDPositions( bool inverseTransform=false );
364 : // get the derivative of the rotation matrix respect to reference
365 : // note that the this transformation overlap the reference onto position
366 : // if inverseTransform=true then aligns the positions onto reference
367 : Matrix<std::vector<Vector> > getDRotationDReference(bool inverseTransform=false );
368 : const std::array<std::array<Tensor,3>,3> & getDRotationDRr01() const;
369 : };
370 :
371 : }
372 :
373 : #endif
374 :
|
