| 1 | 
gezelter | 
1335 | 
/* | 
| 2 | 
  | 
  | 
 *  ******************************************************************* | 
| 3 | 
  | 
  | 
 * | 
| 4 | 
  | 
  | 
 *  rmsd.h  | 
| 5 | 
  | 
  | 
 *  (c) 2005 Bosco K Ho | 
| 6 | 
  | 
  | 
 *  | 
| 7 | 
  | 
  | 
 *  Implementation of the Kabsch algorithm to find the RMSD, and  | 
| 8 | 
  | 
  | 
 *  the least-squares rotation matrix for a superposition between  | 
| 9 | 
  | 
  | 
 *  two sets of vectors. | 
| 10 | 
  | 
  | 
 * | 
| 11 | 
  | 
  | 
 *  This implementation is completely self-contained. No other dependencies. | 
| 12 | 
  | 
  | 
 * | 
| 13 | 
  | 
  | 
 *  ************************************************************************** | 
| 14 | 
  | 
  | 
 * | 
| 15 | 
  | 
  | 
 *  This program is free software; you can redistribute it and/or modify | 
| 16 | 
  | 
  | 
 *  it under the terms of the GNU Lesser General Public License as published | 
| 17 | 
  | 
  | 
 *  by the Free Software Foundation; either version 2.1 of the License, or (at | 
| 18 | 
  | 
  | 
 *  your option) any later version. | 
| 19 | 
  | 
  | 
 *   | 
| 20 | 
  | 
  | 
 *  This program is distributed in the hope that it will be useful,  but | 
| 21 | 
  | 
  | 
 *  WITHOUT ANY WARRANTY; without even the implied warranty of | 
| 22 | 
  | 
  | 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU | 
| 23 | 
  | 
  | 
 *  Lesser General Public License for more details.  | 
| 24 | 
  | 
  | 
 *   | 
| 25 | 
  | 
  | 
 *  You should have received a copy of the GNU Lesser General Public License  | 
| 26 | 
  | 
  | 
 *  along with this program; if not, write to the Free Software Foundation,  | 
| 27 | 
  | 
  | 
 *  Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. | 
| 28 | 
  | 
  | 
 * | 
| 29 | 
  | 
  | 
 *  ************************************************************************** | 
| 30 | 
  | 
  | 
 *  | 
| 31 | 
  | 
  | 
 */ | 
| 32 | 
  | 
  | 
 | 
| 33 | 
  | 
  | 
#ifndef MATH_RMSD_HPP | 
| 34 | 
  | 
  | 
#define MATH_RMSD_HPP | 
| 35 | 
  | 
  | 
 | 
| 36 | 
  | 
  | 
#include <vector> | 
| 37 | 
  | 
  | 
 | 
| 38 | 
  | 
  | 
#include "config.h" | 
| 39 | 
  | 
  | 
#include "math/SquareMatrix3.hpp" | 
| 40 | 
  | 
  | 
 | 
| 41 | 
gezelter | 
1390 | 
namespace OpenMD { | 
| 42 | 
gezelter | 
1335 | 
 | 
| 43 | 
  | 
  | 
  class RMSD{ | 
| 44 | 
  | 
  | 
  public: | 
| 45 | 
  | 
  | 
    RMSD(); | 
| 46 | 
  | 
  | 
    RMSD(std::vector<Vector3d> ref) { | 
| 47 | 
  | 
  | 
      set_reference_structure(ref); | 
| 48 | 
  | 
  | 
    } | 
| 49 | 
  | 
  | 
    virtual ~RMSD() { }; | 
| 50 | 
  | 
  | 
 | 
| 51 | 
  | 
  | 
    void set_reference_structure(std::vector<Vector3d> ref) { | 
| 52 | 
  | 
  | 
      ref_ = ref; | 
| 53 | 
  | 
  | 
      ref_com = V3Zero; | 
| 54 | 
  | 
  | 
       | 
| 55 | 
gezelter | 
1782 | 
      for (unsigned int n = 0; n < ref_.size(); n++) { | 
| 56 | 
gezelter | 
1335 | 
        ref_com += ref_[n]; | 
| 57 | 
  | 
  | 
      } | 
| 58 | 
  | 
  | 
      ref_com /= (RealType)ref.size();       | 
| 59 | 
  | 
  | 
    } | 
| 60 | 
  | 
  | 
     | 
| 61 | 
  | 
  | 
    /* | 
| 62 | 
  | 
  | 
     * calculate_rmsd() | 
| 63 | 
  | 
  | 
     * | 
| 64 | 
  | 
  | 
     *   given a vector of Vector3 coordinates, constructs | 
| 65 | 
  | 
  | 
     *    - mov_com: the centre of mass of the mov list | 
| 66 | 
  | 
  | 
     *    - mov_to_ref: vector between the com of mov and ref | 
| 67 | 
  | 
  | 
     *    - U: the rotation matrix for least-squares, usage of | 
| 68 | 
  | 
  | 
     * | 
| 69 | 
  | 
  | 
     *   returns | 
| 70 | 
  | 
  | 
     *    - rmsd: measures similarity between the vectors | 
| 71 | 
  | 
  | 
     */ | 
| 72 | 
  | 
  | 
    RealType calculate_rmsd(std::vector<Vector3d> mov, | 
| 73 | 
  | 
  | 
                            Vector3d mov_com, | 
| 74 | 
cli2 | 
1349 | 
                            Vector3d mov_to_ref); | 
| 75 | 
gezelter | 
1335 | 
     | 
| 76 | 
  | 
  | 
    /*  | 
| 77 | 
  | 
  | 
     * optimal_superposition() | 
| 78 | 
  | 
  | 
     * | 
| 79 | 
  | 
  | 
     *   Returns best-fit rotation matrix  | 
| 80 | 
  | 
  | 
     */ | 
| 81 | 
cli2 | 
1349 | 
    RotMat3x3d optimal_superposition(std::vector<Vector3d> mov, | 
| 82 | 
  | 
  | 
                            Vector3d mov_com, | 
| 83 | 
  | 
  | 
                            Vector3d mov_to_ref); | 
| 84 | 
gezelter | 
1335 | 
 | 
| 85 | 
cli2 | 
1349 | 
     | 
| 86 | 
gezelter | 
1335 | 
  protected: | 
| 87 | 
  | 
  | 
    std::vector<Vector3d> ref_; | 
| 88 | 
  | 
  | 
    Vector3d ref_com; | 
| 89 | 
  | 
  | 
 | 
| 90 | 
  | 
  | 
  }; | 
| 91 | 
  | 
  | 
} | 
| 92 | 
  | 
  | 
     | 
| 93 | 
  | 
  | 
 | 
| 94 | 
  | 
  | 
#endif |