1 /********************************************************************** 2 align.h - Align two molecules or vectors of vector3 3 4 Copyright (C) 2010 by Noel M. O'Boyle 5 6 This file is part of the Open Babel project. 7 For more information, see <http://openbabel.org/> 8 9 This program is free software; you can redistribute it and/or modify 10 it under the terms of the GNU General Public License as published by 11 the Free Software Foundation version 2 of the License. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 ***********************************************************************/ 18 19 #ifndef OB_ALIGN_H 20 #define OB_ALIGN_H 21 22 #include <openbabel/mol.h> 23 #include <openbabel/math/vector3.h> 24 #include <openbabel/math/matrix3x3.h> 25 #include <openbabel/isomorphism.h> 26 #include <Eigen/Core> 27 28 namespace OpenBabel 29 { 30 /** 31 * \class OBAlign align.h <openbabel/math/align.h> 32 * \brief Perform a least-squares alignment of two molecules or two vectors of vector3 objects 33 * 34 * This class may be used to perform a least-squares alignment of two OBMol 35 * objects or two vector<vector3> objects. The Kabsch algorithm is used 36 * for the alignment. 37 * 38 * During the alignment, the Target is aligned to the Reference. Note that 39 * mutiple alignments to the same Reference will be much faster than multiple 40 * alignments to the same Target. When carrying out multiple alignments, 41 * a single OBAlign instance should be reused by calling 42 * SetTarget() or SetTargetMol() for each additional Target and then calling 43 * Align(). 44 * 45 * When aligning molecules, the atoms of the two molecules must be in the same 46 * order for the results to be sensible. By default, hydrogens are not 47 * included in the least-squares fitting procedure (set @p includeH to 48 * true if you wish to include them) and so the resulting RMSD is the 49 * heavy-atom RMSD (which is usually what you want). 50 * 51 * By default, symmetry is taken 52 * into account when comparing molecules. For example, if a benzene is flipped 53 * by 180 degrees along one of its 2-fold symmetry axes, it will only have an 54 * RMSD of 0 (with respect to its original orientation) if symmetry is 55 * enabled. To turn off handling of symmetry set @p symmetry to false (this 56 * will speed up the alignment). 57 * 58 * Note that neither the Target nor the Reference 59 * are modified by the algorithm. As a result, to update a TargetMol with the 60 * new coordinates from the alignment, you need to use UpdateCoords(). 61 * 62 * @since version 2.3 63 */ 64 class OBAPI OBAlign { 65 public: 66 ///@name Constructors 67 //@{ 68 /** 69 * If this constructor is used, the Target and Reference must be 70 * set using SetRef/SetRefMol and SetTarget/SetTargetMol before running 71 * the alignment. 72 */ 73 OBAlign(bool includeH=false, bool symmetry=true); 74 /** 75 * Align two molecules. 76 */ 77 OBAlign(const OBMol &refmol, const OBMol &targetmol, bool includeH=false, bool symmetry=true); 78 /** 79 * Align two vectors of vector3 objects. 80 */ 81 OBAlign(const std::vector<vector3> &ref, const std::vector<vector3> &target); 82 //@} 83 84 ///@name Partial Setup 85 //@{ 86 /** 87 * Set the Reference (to which the Target will be aligned) in 88 * terms of a vector of vector3 objects. Note that it is faster 89 * to perform multiple alignments to the same Reference, rather than 90 * multiple alignments to the same Target. 91 */ 92 void SetRef(const std::vector<vector3> &ref); 93 /** 94 * Set the Target (which will be aligned to the Reference) in 95 * terms of a vector of vector3 objects. 96 */ 97 void SetTarget(const std::vector<vector3> &target); 98 /** 99 * Set the Reference Molecule (to which the Target Molecule must 100 * be aligned). Note that is faster to perform multiple alignments 101 * to the same Reference Molecule, rather than multple alignments 102 * to the same Target Molecule. 103 */ 104 void SetRefMol(const OBMol &refmol); 105 /** 106 * Set the Target Molecule (which will be aligned to the 107 * Reference Molecule). 108 */ 109 void SetTargetMol(const OBMol &targetmol); 110 //@} 111 112 ///@name Execute the alignment 113 //@{ 114 /** 115 * Align the Target to the Reference using a least-squares alignment. 116 */ 117 bool Align(); 118 119 enum AlignMethod { 120 Kabsch = 0, // Returns matrix and RMSD 121 QCP = 1 // Returns just RMSD (fast) 122 }; 123 124 void SetMethod(enum AlignMethod method); 125 //@} 126 127 ///@name Access the result of the alignment 128 //@{ 129 /** 130 * Return the root mean squared deviation of the target from 131 * the reference. This function should only 132 * be called after running the alignment (using Align()). 133 */ 134 double GetRMSD(); 135 /** 136 * Return the rotation matrix associated with the least squares 137 * alignment. This function should only 138 * be called after running the alignment (using Align()). 139 * 140 * The following example shows how to use the rotation matrix 141 * to rotate all of the atoms in a molecule. 142 * \code 143 * matrix3x3 rotmatrix = align.GetRotMatrix(); 144 * for (unsigned int i = 1; i <= mol.NumAtoms(); ++i) { 145 * vector3 tmpvec = mol.GetAtom(i)->GetVector(); 146 * tmpvec *= rotmatrix; //apply the rotation 147 * mol.GetAtom(i)->SetVector(tmpvec); 148 * } 149 * \endcode 150 * Note that if you wish to use the rotation matrix to find the 151 * aligned coordinates (that is, the same coordinates returned by 152 * GetAlignment()), you should first translate the set of coordinates 153 * to the origin by subtracting the centroid, then apply the rotation, 154 * and finally add the centroid of the reference coordinates. 155 */ 156 matrix3x3 GetRotMatrix(); 157 /** 158 * Return the actual alignment of the Target to the Reference 159 * in terms of a vector of vector3 objects. If you want an OBMol 160 * with the aligned coordinates, you should use UpdateCoords() instead. 161 * This function should only 162 * be called after running the alignment (using Align()). 163 */ 164 std::vector<vector3> GetAlignment(); 165 /** 166 * Set the coordinates of an OBMol to those from the alignment. 167 * This function should only 168 * be called after running the alignment (using Align()). 169 */ 170 bool UpdateCoords(OBMol* target); 171 //@} 172 173 private: 174 bool _ready; 175 bool _fail; 176 bool _symmetry; 177 bool _includeH; 178 enum AlignMethod _method; 179 double _rmsd; 180 OBBitVec _frag_atoms; 181 Automorphisms _aut; 182 const OBMol* _prefmol; 183 const OBMol* _ptargetmol; 184 Eigen::MatrixXd _rotMatrix; 185 Eigen::Vector3d _ref_centr, _target_centr; 186 const std::vector<vector3> *_pref; 187 const std::vector<vector3> *_ptarget; 188 std::vector<vector3> _refmol_coords; 189 std::vector<vector3> _targetmol_coords; 190 Eigen::MatrixXd _result; 191 Eigen::MatrixXd _mref, _mtarget; 192 void VectorsToMatrix(const std::vector<vector3> *pcoords, Eigen::MatrixXd &coords); 193 Eigen::Vector3d MoveToOrigin(Eigen::MatrixXd &coords); 194 void SimpleAlign(const Eigen::MatrixXd &mtarget); 195 void TheobaldAlign(const Eigen::MatrixXd &mtarget); 196 // Generate a mapping from the permutation map to the index of 197 // correct column in _mtarget. Need to handle the fact that the 198 // permutation group contains non-fragment atoms. 199 // For example, map(213465) will be converted to newidx(102354). 200 // If the atom with Idx=3 is not in the fragment, it will be 201 // converted to newidx(10X243) instead. 202 std::vector<unsigned int> _newidx; 203 }; 204 } 205 206 #endif // OB_ALIGN_H 207 208 /// @file align.h 209 /// @brief Perform a least-squares alignment of two molecules or two vectors of vector3 objects 210