1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2    Copyright (c) 2011-2021 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
calculateCenter(std::vector<Vector> & p,std::vector<double> & w)83   Vector calculateCenter(std::vector<Vector> &p,std::vector<double> &w) {
84     plumed_massert(p.size()==w.size(),"mismatch in dimension of position/align arrays while calculating the center");
85     unsigned n; n=p.size();
86     Vector c; c.zero();
87     for(unsigned i=0; i<n; i++)c+=p[i]*w[i];
88     return c;
89   };
90 // removes the center for the position provided
removeCenter(std::vector<Vector> & p,Vector & c)91   void removeCenter(std::vector<Vector> &p, Vector &c) {
92     unsigned n; n=p.size();
93     for(unsigned i=0; i<n; i++)p[i]-=c;
94   };
95 // add center
addCenter(std::vector<Vector> & p,Vector & c)96   void addCenter(std::vector<Vector> &p, 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       Tensor & rotationPosClose,
142       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> &centeredreference, 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, Tensor & rotationPosClose, 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)
getMatrixFromDRot(Matrix<std::vector<Vector>> & drotdpos,const unsigned & i,const unsigned & a)251   static  Tensor getMatrixFromDRot(Matrix< std::vector<Vector> > &drotdpos, const unsigned &i, const unsigned &a) {
252     Tensor t;
253     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     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     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     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 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
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)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
RMSDCoreData(const std::vector<double> & a,const std::vector<double> & d,const std::vector<Vector> & p,const std::vector<Vector> & r)308   RMSDCoreData(const std::vector<double> &a,const std::vector<double> &d,const std::vector<Vector> &p, const std::vector<Vector> &r):
309     alEqDis(false),distanceIsMSD(false),hasDistance(false),isInitialized(false),safe(false),
310     creference_is_calculated(false),creference_is_removed(false),
311     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   {cpositions.zero(); creference.zero();};
313 
314   // set the center on the fly without subtracting
calcPositionsCenter()315   void calcPositionsCenter() {
316     plumed_massert(!cpositions_is_calculated,"the center was already calculated");
317     cpositions.zero(); for(unsigned i=0; i<positions.size(); i++) {cpositions+=positions[i]*align[i];} cpositions_is_calculated=true;
318   }
calcReferenceCenter()319   void calcReferenceCenter() {
320     plumed_massert(!creference_is_calculated,"the center was already calculated");
321     creference.zero(); for(unsigned i=0; i<reference.size(); i++) {creference+=reference[i]*align[i];} creference_is_calculated=true;
322   };
323   // assume the center is given externally
setPositionsCenter(Vector v)324   void setPositionsCenter(Vector v) {plumed_massert(!cpositions_is_calculated,"You are setting the center two times!"); cpositions=v; cpositions_is_calculated=true;};
setReferenceCenter(Vector v)325   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
setPositionsCenterIsRemoved(bool t)327   void setPositionsCenterIsRemoved(bool t) {cpositions_is_removed=t;};
setReferenceCenterIsRemoved(bool t)328   void setReferenceCenterIsRemoved(bool t) {creference_is_removed=t;};
getPositionsCenterIsRemoved()329   bool getPositionsCenterIsRemoved() {return cpositions_is_removed;};
getReferenceCenterIsRemoved()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, Tensor & rotationPosClose, 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 
375