Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 :
31 : namespace PLMD {
32 :
33 : class Log;
34 : class PDB;
35 :
36 : /** \ingroup TOOLBOX
37 : A class that implements RMSD calculations
38 : This is a class that implements the various infrastructure to calculate the
39 : RMSD or MSD respect a given frame. It can be done through an optimal alignment scheme
40 : as Kearsley or, more simply, by resetting the center of mass.
41 : This is the class that decides this. A very simple use is
42 : \verbatim
43 : #include "tools/PDB.h"
44 : #include "tools/RMSD.h"
45 : #include "tools/Vector.h"
46 : using namespace PLMD;
47 : RMSD rmsd;
48 : PDB pdb;
49 : // get the pdb (see PDB documentation)
50 : pdb.read("file.pdb",true,1.0);
51 : string type;
52 : type.assign("OPTIMAL");
53 : // set the reference and the type
54 : rmsd.set(pdb,type);
55 : // this calculates the rmsd and the derivatives
56 : vector<Vector> derivs;
57 : double val;
58 : val=rmsd.calculate(getPositions(),derivs,true);
59 : \endverbatim
60 :
61 : **/
62 :
63 6590 : class RMSD
64 : {
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 1877 : Vector calculateCenter(std::vector<Vector> &p,std::vector<double> &w) {
83 1877 : plumed_massert(p.size()==w.size(),"mismatch in dimension of position/align arrays while calculating the center");
84 1877 : unsigned n; n=p.size();
85 1877 : Vector c; c.zero();
86 1877 : for(unsigned i=0; i<n; i++)c+=p[i]*w[i];
87 1877 : return c;
88 : };
89 : // removes the center for the position provided
90 3750 : void removeCenter(std::vector<Vector> &p, Vector &c) {
91 3750 : unsigned n; n=p.size();
92 3750 : for(unsigned i=0; i<n; i++)p[i]-=c;
93 3750 : };
94 : // add center
95 1877 : void addCenter(std::vector<Vector> &p, Vector &c) {Vector cc=c*-1.; removeCenter(p,cc);};
96 :
97 : public:
98 : /// Constructor
99 : RMSD();
100 : /// clear the structure
101 : void clear();
102 : /// set reference, align and displace from input pdb structure: evtl remove com from the initial structure and normalize the input weights from the pdb
103 : void set(const PDB&,const std::string & mytype, bool remove_center=true, bool normalize_weights=true);
104 : /// set align displace reference and type from input vectors
105 : 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 );
106 : /// set the type of alignment we are doing
107 : void setType(const std::string & mytype);
108 : /// set reference coordinates, remove the com by using uniform weights
109 : void setReference(const std::vector<Vector> & reference);
110 : /// set weights and remove the center from reference with normalized weights. If the com has been removed, it resets to the new value
111 : void setAlign(const std::vector<double> & align, bool normalize_weights=true, bool remove_center=true);
112 : /// set align
113 : void setDisplace(const std::vector<double> & displace, bool normalize_weights=true);
114 : ///
115 : std::string getMethod();
116 : /// workhorses
117 : double simpleAlignment(const std::vector<double> & align,
118 : const std::vector<double> & displace,
119 : const std::vector<Vector> & positions,
120 : const std::vector<Vector> & reference,
121 : std::vector<Vector> & derivatives,
122 : std::vector<Vector> & displacement,
123 : bool squared=false)const;
124 : template <bool safe,bool alEqDis>
125 : double optimalAlignment(const std::vector<double> & align,
126 : const std::vector<double> & displace,
127 : const std::vector<Vector> & positions,
128 : const std::vector<Vector> & reference,
129 : std::vector<Vector> & DDistDPos, bool squared=false)const;
130 :
131 : template <bool safe,bool alEqDis>
132 : double optimalAlignment_DDistDRef(const std::vector<double> & align,
133 : const std::vector<double> & displace,
134 : const std::vector<Vector> & positions,
135 : const std::vector<Vector> & reference,
136 : std::vector<Vector> & DDistDPos,
137 : std::vector<Vector> & DDistDRef,
138 : bool squared=false) const;
139 :
140 : template <bool safe,bool alEqDis>
141 : double optimalAlignment_SOMA(const std::vector<double> & align,
142 : const std::vector<double> & displace,
143 : const std::vector<Vector> & positions,
144 : const std::vector<Vector> & reference,
145 : std::vector<Vector> & DDistDPos,
146 : std::vector<Vector> & DDistDRef,
147 : bool squared=false) const;
148 :
149 : template <bool safe,bool alEqDis>
150 : double optimalAlignment_DDistDRef_Rot_DRotDPos(const std::vector<double> & align,
151 : const std::vector<double> & displace,
152 : const std::vector<Vector> & positions,
153 : const std::vector<Vector> & reference,
154 : std::vector<Vector> & DDistDPos,
155 : std::vector<Vector> & DDistDRef,
156 : Tensor & Rotation,
157 : Matrix<std::vector<Vector> > &DRotDPos,
158 : bool squared=false) const;
159 :
160 : template <bool safe,bool alEqDis>
161 : double optimalAlignment_DDistDRef_Rot_DRotDPos_DRotDRef(const std::vector<double> & align,
162 : const std::vector<double> & displace,
163 : const std::vector<Vector> & positions,
164 : const std::vector<Vector> & reference,
165 : std::vector<Vector> & DDistDPos,
166 : std::vector<Vector> & DDistDRef,
167 : Tensor & Rotation,
168 : Matrix<std::vector<Vector> > &DRotDPos,
169 : Matrix<std::vector<Vector> > &DRotDRef,
170 : bool squared=false) const;
171 :
172 : template <bool safe,bool alEqDis>
173 : double optimalAlignment_PCA(const std::vector<double> & align,
174 : const std::vector<double> & displace,
175 : const std::vector<Vector> & positions,
176 : const std::vector<Vector> & reference,
177 : std::vector<Vector> & alignedpositions,
178 : std::vector<Vector> & centeredpositions,
179 : std::vector<Vector> & centeredreference,
180 : Tensor & Rotation,
181 : std::vector<Vector> & DDistDPos,
182 : Matrix<std::vector<Vector> > & DRotDPos,
183 : bool squared=false) const ;
184 :
185 : template <bool safe,bool alEqDis>
186 : double optimalAlignment_Fit(const std::vector<double> & align,
187 : const std::vector<double> & displace,
188 : const std::vector<Vector> & positions,
189 : const std::vector<Vector> & reference,
190 : Tensor & Rotation,
191 : Matrix<std::vector<Vector> > & DRotDPos,
192 : std::vector<Vector> & centeredpositions,
193 : Vector & center_positions,
194 : bool squared=false);
195 :
196 :
197 : /// 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
198 : double calculate(const std::vector<Vector> & positions,std::vector<Vector> &derivatives, bool squared=false)const;
199 : /// Other convenience methods:
200 : /// calculate the derivative of distance respect to position(DDistDPos) and reference (DDistDPos)
201 : double calc_DDistDRef( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, std::vector<Vector>& DDistDRef, const bool squared=false );
202 : /// calculate the derivative for SOMA (i.e. derivative respect to reference frame without the matrix derivative)
203 : double calc_SOMA( const std::vector<Vector>& positions, std::vector<Vector> &DDistDPos, std::vector<Vector>& DDistDRef, const bool squared=false );
204 : ///
205 : 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 );
206 : 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 );
207 : /// convenience method to retrieve all the bits and pieces for PCA
208 : 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 ;
209 : /// convenience method to retrieve all the bits and pieces needed by FitToTemplate
210 : 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 );
211 : /// static convenience method to get the matrix i,a from drotdpos (which is a bit tricky)
212 720 : static Tensor getMatrixFromDRot(Matrix< std::vector<Vector> > &drotdpos, const unsigned &i, const unsigned &a) {
213 720 : Tensor t;
214 720 : t[0][0]=drotdpos[0][0][i][a]; t[0][1]=drotdpos[0][1][i][a]; t[0][2]=drotdpos[0][2][i][a];
215 720 : t[1][0]=drotdpos[1][0][i][a]; t[1][1]=drotdpos[1][1][i][a]; t[1][2]=drotdpos[1][2][i][a];
216 720 : t[2][0]=drotdpos[2][0][i][a]; t[2][1]=drotdpos[2][1][i][a]; t[2][2]=drotdpos[2][2][i][a];
217 720 : return t;
218 : };
219 : };
220 :
221 : /// this is a class which is needed to share information across the various non-threadsafe routines
222 : /// so that the public function of rmsd are threadsafe while the inner core can safely share information
223 1719 : class RMSDCoreData
224 : {
225 : private:
226 : bool alEqDis;
227 : bool distanceIsMSD; // default is RMSD but can deliver the MSD
228 : bool hasDistance; // distance is already calculated
229 : bool isInitialized;
230 : bool safe;
231 :
232 : // this need to be copied and they are small, should not affect the performance
233 : Vector creference;
234 : bool creference_is_calculated;
235 : bool creference_is_removed;
236 : Vector cpositions;
237 : bool cpositions_is_calculated;
238 : bool cpositions_is_removed;
239 : bool retrieve_only_rotation;
240 :
241 : // use reference assignment to speed up instead of copying
242 : const std::vector<Vector> &positions;
243 : const std::vector<Vector> &reference;
244 : const std::vector<double> &align;
245 : const std::vector<double> &displace;
246 :
247 : // the needed stuff for distance and more (one could use eigenvecs components and eigenvals for some reason)
248 : double dist;
249 : std::vector<double> eigenvals;
250 : Matrix<double> eigenvecs;
251 : double rr00; // sum of positions squared (needed for dist calc)
252 : double rr11; // sum of reference squared (needed for dist calc)
253 : Tensor rotation; // rotation derived from the eigenvector having the smallest eigenvalue
254 : Tensor drotation_drr01[3][3]; // derivative of the rotation only available when align!=displace
255 : Tensor ddist_drr01;
256 : Tensor ddist_drotation;
257 : std::vector<Vector> d; // difference of components
258 : public:
259 : /// the constructor (note: only references are passed, therefore is rather fast)
260 : /// note: this aligns the reference onto the positions
261 : ///
262 : /// this method assumes that the centers are already calculated and subtracted
263 : 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 ):
264 : alEqDis(false),distanceIsMSD(false),hasDistance(false),isInitialized(false),safe(false),
265 : creference(cr),creference_is_calculated(true),creference_is_removed(true),
266 : 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) {};
267 :
268 : // this constructor does not assume that the positions and reference have the center subtracted
269 1719 : RMSDCoreData(const std::vector<double> &a,const std::vector<double> &d,const std::vector<Vector> &p, const std::vector<Vector> &r):
270 : alEqDis(false),distanceIsMSD(false),hasDistance(false),isInitialized(false),safe(false),
271 : creference_is_calculated(false),creference_is_removed(false),
272 1719 : 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)
273 1719 : {cpositions.zero(); creference.zero();};
274 :
275 : // set the center on the fly without subtracting
276 1719 : void calcPositionsCenter() {
277 1719 : plumed_massert(!cpositions_is_calculated,"the center was already calculated");
278 1719 : cpositions.zero(); for(unsigned i=0; i<positions.size(); i++) {cpositions+=positions[i]*align[i];} cpositions_is_calculated=true;
279 1719 : }
280 0 : void calcReferenceCenter() {
281 0 : plumed_massert(!creference_is_calculated,"the center was already calculated");
282 0 : creference.zero(); for(unsigned i=0; i<reference.size(); i++) {creference+=reference[i]*align[i];} creference_is_calculated=true;
283 0 : };
284 : // assume the center is given externally
285 0 : void setPositionsCenter(Vector v) {plumed_massert(!cpositions_is_calculated,"You are setting the center two times!"); cpositions=v; cpositions_is_calculated=true;};
286 1719 : void setReferenceCenter(Vector v) {plumed_massert(!creference_is_calculated,"You are setting the center two times!"); creference=v; creference_is_calculated=true;};
287 : // the center is already removed
288 1719 : void setPositionsCenterIsRemoved(bool t) {cpositions_is_removed=t;};
289 1719 : void setReferenceCenterIsRemoved(bool t) {creference_is_removed=t;};
290 : bool getPositionsCenterIsRemoved() {return cpositions_is_removed;};
291 : bool getReferenceCenterIsRemoved() {return creference_is_removed;};
292 : // does the core calc : first thing to call after the constructor:
293 : // only_rotation=true does not retrieve the derivatives, just retrieve the optimal rotation (the same calc cannot be exploit further)
294 : void doCoreCalc(bool safe,bool alEqDis, bool only_rotation=false);
295 : // retrieve the distance if required after doCoreCalc
296 : double getDistance(bool squared);
297 : // retrieve the derivative of the distance respect to the position
298 : std::vector<Vector> getDDistanceDPositions();
299 : // retrieve the derivative of the distance respect to the reference
300 : std::vector<Vector> getDDistanceDReference();
301 : // specific version for SOMA calculation (i.e. does not need derivative respect to rotation matrix)
302 : std::vector<Vector> getDDistanceDReferenceSOMA();
303 : // get aligned reference onto position
304 : std::vector<Vector> getAlignedReferenceToPositions();
305 : // get aligned position onto reference
306 : std::vector<Vector> getAlignedPositionsToReference();
307 : // get centered positions
308 : std::vector<Vector> getCenteredPositions();
309 : // get centered reference
310 : std::vector<Vector> getCenteredReference();
311 : // get center of positions
312 : Vector getPositionsCenter();
313 : // get center of reference
314 : Vector getReferenceCenter();
315 : // get rotation matrix (reference ->positions)
316 : Tensor getRotationMatrixReferenceToPositions();
317 : // get rotation matrix (positions -> reference)
318 : Tensor getRotationMatrixPositionsToReference();
319 : // get the derivative of the rotation matrix respect to positions
320 : // note that the this transformation overlap the reference onto position
321 : // if inverseTransform=true then aligns the positions onto reference
322 : Matrix<std::vector<Vector> > getDRotationDPositions( bool inverseTransform=false );
323 : // get the derivative of the rotation matrix respect to reference
324 : // note that the this transformation overlap the reference onto position
325 : // if inverseTransform=true then aligns the positions onto reference
326 : Matrix<std::vector<Vector> > getDRotationDReference(bool inverseTransform=false );
327 : };
328 :
329 : }
330 :
331 : #endif
332 :
|
