superimpose.h

Go to the documentation of this file.

#ifndef SUPERIMPOSE_H
#define SUPERIMPOSE_H

#include <iostream>
#include <string>
#include <vector>
#include <math.h>
#include "Utils/constants.h"

// - BOOST - //
/*
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/banded.hpp> // for diagonal matrix
#include <boost/numeric/bindings/blas/blas3.hpp>
#include "boost/numeric/bindings/traits/ublas_matrix.hpp"
#include "boost/numeric/bindings/traits/ublas_vector.hpp"

namespace ublas  = boost::numeric::ublas;
namespace blas   = boost::numeric::bindings::blas;
*/
#include <Eigen/Dense>
using namespace Eigen;

namespace MTKpp
{

class molecule;
class vector3d;

// ============================================================
// Class : superimpose()
// ------------------------------------------------------------
// ============================================================
class superimpose
{
public:

    superimpose();

    superimpose(int i);

    virtual ~superimpose();

    double    fit(molecule* pMoleculeA, molecule* pMoleculeB);

    double    fit(double coordsA[][3], double coordsB[][3], int n);

    int       fit(molecule* pMoleculeA, molecule* pMoleculeB,
                  std::vector<int> cor, int type);

    double    rmsd(molecule* pMoleculeA, molecule* pMoleculeB);

    double    rmsd(molecule* pMoleculeA, molecule* pMoleculeB,
              std::vector<std::vector<int> > &correspondenceMatrices);

    double    rmsd(molecule* pMoleculeA, std::vector< vector3d > &molBCoords,
              std::vector<std::vector<int> > &correspondenceMatrices, int &cor);

    double    rmsdNoFit(molecule* pMoleculeA,
              std::vector< vector3d > &molBCoords,
              std::vector<std::vector<int> > &correspondenceMatrices);

    void       getRotationMatrix(double centerA[3], double centerB[3],
                                 double rotMat[][3]);

    int       getRotationMatrix(double coordsA[][3], double coordsB[][3],
                                     int n, double rotMat[][3]);

    int       getRotationMatrix(double coordsA[][3], double coordsB[][3],
                                double centerA[3], double centerB[3],
                                int n, double rotMat[][3]);

    double    calculateRMSD(double coordsA[][3], double coordsB[][3]);

    double    minRMSD(molecule* pMoleculeA, molecule* pMoleculeB, int type,
                      int &nFittedAtoms);

    int initializeCorrespondences(molecule* pMol,
                 int type,
                 std::vector<std::vector<int> > &correspondenceMatrices);

    int initializeCorrespondences(molecule* pMolB,
                 molecule* pMolA,
                 int type,
                 std::vector<std::vector<int> > &correspondenceMatrices);

    void      center(double coords[][3], double center[3]);

    void      updateCoords(double coords[][3], int n, double c1[3],
                           double c2[3], double rotMat[][3]);

protected: // functions

    void      coordinateDiff(double coordsA[][3], double coordsB[][3],
                             double itsDXM[][3], double itsDXP[][3]);

/*
    void      buildQuaternion(ublas::matrix<double, ublas::column_major>& quaternion,
                              double itsDXM[][3], double itsDXP[][3]);
*/
    void      buildQuaternion(Eigen::MatrixXd& quaternion,
                              double itsDXM[][3], double itsDXP[][3]);

/*
    void      buildRotation(double t[][3],
              ublas::matrix<double, ublas::column_major> eigenvectors);
*/
    void      buildRotation(double t[][3],
              Eigen::MatrixXd eigenvectors);

protected: // data

    molecule* pMoleculeA;

    molecule* pMoleculeB;

    int       nAtoms;

    double    centerA[3];

    double    centerB[3];

    double    rotMat[3][3];

    double    dRMSD;

    int* genMatchMatrix;

    int correspondenceType;

    int                           nHeavyAtomsA;

    int*                          heavyAtomIndicesA;

    int                           nHeavyAtomsB;

    int*                          heavyAtomIndicesB;

};

} // MTKpp namespace

#endif // SUPERIMPOSE_H