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> ¢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, 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