kabschalignment.java

化学图形处理软件

     *
     * @param ac1 An {@link IAtomContainer}
     * @param ac2 An {@link IAtomContainer}. This AtomContainer will have its coordinates rotated
     *            so that the RMDS is minimzed to the coordinates of the first one
     * @param wts A vector atom weights.           
     * @throws CDKException if the number of atom's are not the same in the two AtomContainer's or
     *                         length of the weight vector is not the same as number of atoms.
     */            
    public KabschAlignment(IAtomContainer ac1, IAtomContainer ac2, double[] wts) throws CDKException {
        if (ac1.getAtomCount() != ac2.getAtomCount()) {
            throw new CDKException("The AtomContainer's being aligned must have the same number of atoms");
        }
        if (ac1.getAtomCount() != wts.length) {
            throw new CDKException("Number of weights must equal number of atoms");
        }
        this.npoint = ac1.getAtomCount();
        this.p1 = getPoint3dArray(ac1);
        this.p2 = getPoint3dArray(ac2);
        this.wts = new double[npoint];
        for (int i = 0; i < npoint; i++) this.wts[i] = wts[i];

        this.atwt1 = getAtomicMasses(ac1);
        this.atwt2 = getAtomicMasses(ac2);
    }

    /**
     * Perform an alignment.
     *
     * This method aligns to set of atoms which should have been specified
     * prior to this call
     */
    public void align() {
        
        Matrix tmp;

       // get center of gravity and translate both to 0,0,0 
        this.cm1 = new Point3d();
        this.cm2 = new Point3d();

        this.cm1 = getCenterOfMass(p1, atwt1);
        this.cm2 = getCenterOfMass(p2, atwt2);

        // move the points 
        for (int i = 0; i < this.npoint; i++) {
            p1[i].x = p1[i].x - this.cm1.x;
            p1[i].y = p1[i].y - this.cm1.y;
            p1[i].z = p1[i].z - this.cm1.z;

            p2[i].x = p2[i].x - this.cm2.x;
            p2[i].y = p2[i].y - this.cm2.y;
            p2[i].z = p2[i].z - this.cm2.z;
        }
        
        // get the R matrix 
        double[][] tR = new double[3][3];
        for (int i = 0; i < this.npoint; i++) {
            wts[i] = 1.0;
        }
        for (int i = 0; i < this.npoint; i++) {
            tR[0][0] += p1[i].x * p2[i].x * wts[i];
            tR[0][1] += p1[i].x * p2[i].y * wts[i];
            tR[0][2] += p1[i].x * p2[i].z * wts[i];

            tR[1][0] += p1[i].y * p2[i].x * wts[i];
            tR[1][1] += p1[i].y * p2[i].y * wts[i];
            tR[1][2] += p1[i].y * p2[i].z * wts[i];

            tR[2][0] += p1[i].z * p2[i].x * wts[i];
            tR[2][1] += p1[i].z * p2[i].y * wts[i];
            tR[2][2] += p1[i].z * p2[i].z * wts[i];
        }
        double[][] R = new double[3][3];
        tmp = new Matrix(tR);
        R = tmp.transpose().getArray();


        // now get the RtR (=R'R) matrix 
        double[][] RtR = new double[3][3];
        Matrix jamaR = new Matrix(R);
        tmp = tmp.times(jamaR);
        RtR = tmp.getArray();

        // get eigenvalues of RRt (a's)
        Matrix jamaRtR = new Matrix(RtR);
        EigenvalueDecomposition ed = jamaRtR.eig();
        double[] mu = ed.getRealEigenvalues();
        double[][] a = ed.getV().getArray();


 
        // Jama returns the eigenvalues in increasing order so
        // swap the eigenvalues and vectors
        double tmp2 = mu[2];
        mu[2] = mu[0];
        mu[0] = tmp2;
        for (int i = 0; i < 3; i++) {
            tmp2 = a[i][2];
            a[i][2] = a[i][0];
            a[i][0] = tmp2;
        }

        // make sure that the a3 = a1 x a2
        a[0][2] = (a[1][0]*a[2][1]) - (a[1][1]*a[2][0]);
        a[1][2] = (a[0][1]*a[2][0]) - (a[0][0]*a[2][1]);
        a[2][2] = (a[0][0]*a[1][1]) - (a[0][1]*a[1][0]);

        // lets work out the b vectors
        double[][] b = new double[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                for (int k = 0; k < 3; k++) {
                    b[i][j] += R[i][k] * a[k][j];
                }
                b[i][j] = b[i][j] / Math.sqrt(mu[j]);
            }
        }
        
        // normalize and set b3 = b1 x b2
        double norm1 = 0.;
        double norm2 = 0.;
        for (int i = 0; i < 3; i++) {
            norm1 += b[i][0]*b[i][0];
            norm2 += b[i][1]*b[i][1];
        }
        norm1 = Math.sqrt(norm1);
        norm2 = Math.sqrt(norm2);
        for (int i = 0; i < 3; i++) {
            b[i][0] = b[i][0] / norm1;
            b[i][1] = b[i][1] / norm2;
        }
        b[0][2] = (b[1][0]*b[2][1]) - (b[1][1]*b[2][0]);
        b[1][2] = (b[0][1]*b[2][0]) - (b[0][0]*b[2][1]);
        b[2][2] = (b[0][0]*b[1][1]) - (b[0][1]*b[1][0]);

        // get the rotation matrix
        double[][] tU = new double[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                for (int k = 0; k < 3; k++) {
                    tU[i][j] += b[i][k]*a[j][k];
                }
            }
        }

        // take the transpose
        U = new double[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                U[i][j] = tU[j][i];
            }
        }

        // now eval the RMS error
        // first, rotate the second set of points and ...
        this.rp = new Point3d[ this.npoint ];
        for (int i = 0; i < this.npoint; i++) {
            this.rp[i] = new Point3d(
                    U[0][0]*p2[i].x + U[0][1]*p2[i].y + U[0][2]*p2[i].z,
                    U[1][0]*p2[i].x + U[1][1]*p2[i].y + U[1][2]*p2[i].z,
                    U[2][0]*p2[i].x + U[2][1]*p2[i].y + U[2][2]*p2[i].z
                    );
        }

        // ... then eval rms
        double rms = 0.;
        for (int i = 0; i < this.npoint; i++) {
            rms += (p1[i].x-this.rp[i].x)*(p1[i].x-this.rp[i].x) +
                (p1[i].y-this.rp[i].y)*(p1[i].y-this.rp[i].y) +
                (p1[i].z-this.rp[i].z)*(p1[i].z-this.rp[i].z);
        }
        this.rmsd = Math.sqrt(rms / this.npoint);
    }

    /**
     * Returns the RMSD from the alignment.
     *
     * If align() has not been called the return value is -1.0
     *
     * @return The RMSD for this alignment
     * @see #align
     */
    public double getRMSD() {
        return(this.rmsd);
    }

    /**
     * Returns the rotation matrix (u).
     *
     * @return A double[][] representing the rotation matrix
     * @see #align
     */
    public double[][] getRotationMatrix() {
        return(this.U);
    }

    /**
     * Returns the center of mass for the first molecule or fragment used in the calculation.
     *
     * This method is useful when using this class to align the coordinates
     * of two molecules and them displaying them superimposed. Since the center of
     * mass used during the alignment may not be based on the whole molecule (in 
     * general common substructures are aligned), when preparing molecules for display
     * the first molecule should be translated to the center of mass. Then displaying the
     * first molecule and the rotated version of the second one will result in superimposed
     * structures.
     * 
     * @return A Point3d containing the coordinates of the center of mass
     */
    public Point3d getCenterOfMass() {
        return(this.cm1);
    }

    /**
     * Rotates the {@link IAtomContainer} coordinates by the rotation matrix.
     *
     * In general if you align a subset of atoms in a AtomContainer
     * this function can be applied to the whole AtomContainer to rotate all
     * atoms. This should be called with the second AtomContainer (or Atom[])
     * that was passed to the constructor.
     *
     * Note that the AtomContainer coordinates also get translated such that the
     * center of mass of the original fragment used to calculate the alignment is at the origin.
     *
     * @param ac The {@link IAtomContainer} whose coordinates are to be rotated
     */
    public void rotateAtomContainer(IAtomContainer ac)  {
        Point3d[] p = getPoint3dArray( ac );
        for (int i = 0; i < ac.getAtomCount(); i++) {
            // translate the the origin we have calculated
            p[i].x = p[i].x - this.cm2.x;
            p[i].y = p[i].y - this.cm2.y;
            p[i].z = p[i].z - this.cm2.z;

            // do the actual rotation
            ac.getAtom(i).setPoint3d( 
            	new Point3d(
            		U[0][0]*p[i].x + U[0][1]*p[i].y + U[0][2]*p[i].z,
            		U[1][0]*p[i].x + U[1][1]*p[i].y + U[1][2]*p[i].z,
            		U[2][0]*p[i].x + U[2][1]*p[i].y + U[2][2]*p[i].z
            	)
            );
        }
    }

 }