kabschalignment.java

化学图形处理软件

/*
 * Copyright (C) 2004-2007  The Chemistry Development Kit (CDK) project
 * 
 * Contact: [e-mail]
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 
 */

package org.openscience.cdk.geometry.alignment;

import javax.vecmath.Point3d;

import org.openscience.cdk.Atom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.tools.LoggingTool;
import org.openscience.cdk.config.IsotopeFactory;
import org.openscience.cdk.exception.CDKException;

import Jama.EigenvalueDecomposition;
import Jama.Matrix;

/**
 * Aligns two structures to minimize the RMSD using the Kabsch algorithm.
 *
 * <p>This class is an implementation of the Kabsch algorithm ({@cdk.cite KAB76}, {@cdk.cite KAB78})
 * and evaluates the optimal rotation matrix (U) to minimize the RMSD between the two structures.
 * Since the algorithm assumes that the number of points are the same in the two structures
 * it is the job of the caller to pass the proper number of atoms from the two structures. Constructors
 * which take whole <code>AtomContainer</code>'s are provided but they should have the same number
 * of atoms.
 * The algorithm allows for the use of atom weightings and by default all points are given a weight of 1.0
 *
 * <p>Example usage can be:
 * <pre>
 * AtomContainer ac1, ac2;
 *
 * try {
 *    KabschAlignment sa = new KabschAlignment(ac1.getAtoms(),ac2.getAtoms());
 *    sa.align();
 *    System.out.println(sa.getRMSD());
 * } catch (CDKException e){}
 * </pre>
 * In many cases, molecules will be aligned based on some common substructure.
 * In this case the center of masses calculated during alignment refer to these
 * substructures rather than the whole molecules. To superimpose the molecules
 * for display, the second molecule must be rotated and translated by calling
 * <code>rotateAtomContainer</code>. However, since this will also translate the
 * second molecule, the first molecule should also be translated to the center of mass
 * of the substructure specifed for this molecule. This center of mass can be obtained
 * by a call to <code>getCenterOfMass</code> and then manually translating the coordinates.
 * Thus an example would be
 * <pre>
 * AtomContainer ac1, ac2;  // whole molecules
 * Atom[] a1, a2;           // some subset of atoms from the two molecules
 * KabschAlignment sa;
 * 
 * try {
 *    sa = new KabschAlignment(a1,a2);
 *    sa.align();
 * } catch (CDKException e){}
 *
 * Point3d cm1 = sa.getCenterOfMass();
 * for (int i = 0; i &lt; ac1.getAtomCount(); i++) {
 *    Atom a = ac1.getAtomAt(i);
 *    a.setX3d( a.getPoint3d().x - cm1.x );
 *    a.setY3d( a.getPoint3d().y - cm1.y );
 *    a.setY3d( a.getPoint3d().z - cm1.z );
 * }
 * sa.rotateAtomContainer(ac2);
 *
 * // display the two AtomContainer's
 *</pre>
 * 
 * @author           Rajarshi Guha
 * @cdk.created      2004-12-11
 * @cdk.builddepends Jama-1.0.1.jar
 * @cdk.depends      Jama-1.0.1.jar
 * @cdk.dictref      blue-obelisk:alignmentKabsch
 */ 
public class KabschAlignment {

	private LoggingTool logger = new LoggingTool(KabschAlignment.class);
	
    private double[][] U;
    private double rmsd = -1.0;
    private Point3d[] p1,p2,rp; // rp are the rotated coordinates
    private double[] wts;
    private int npoint;
    private Point3d cm1, cm2;
    private double[] atwt1, atwt2;
    
    private Point3d[] getPoint3dArray(org.openscience.cdk.interfaces.IAtom[] a) {
        Point3d[] p = new Point3d[ a.length ];
        for (int i = 0; i < a.length; i++) {
            p[i] = new Point3d( a[i].getPoint3d() );
        }
        return(p);
    }

    private Point3d[] getPoint3dArray(org.openscience.cdk.interfaces.IAtomContainer ac) {
        Point3d[] p = new Point3d[ ac.getAtomCount() ];
        for (int i = 0; i < ac.getAtomCount(); i++) {
            p[i] = new Point3d( ac.getAtom(i).getPoint3d() );
        }
        return(p);
    }
    
    private double[] getAtomicMasses(org.openscience.cdk.interfaces.IAtom[] a) {
        double[] am = new double[a.length];
        IsotopeFactory factory = null;
        try {
            factory = IsotopeFactory.getInstance(a[0].getBuilder());
        } catch (Exception e) {
        	logger.error("Error while instantiating the isotope factory: ",
        		e.getMessage());
            logger.debug(e);
        }

        for (int i = 0; i < a.length; i++) {
            am[i] = factory.getMajorIsotope( a[i].getSymbol() ).getExactMass();
        }
        return(am);
    }
    
    private double[] getAtomicMasses(org.openscience.cdk.interfaces.IAtomContainer ac) {
        double[] am = new double[ac.getAtomCount()];
        IsotopeFactory factory = null;
        try {
            factory = IsotopeFactory.getInstance(ac.getAtom(0).getBuilder());
        } catch (Exception e) {
        	logger.error("Error while instantiating the isotope factory: ",
            	e.getMessage());
        	logger.debug(e);
        }

        for (int i = 0; i < ac.getAtomCount(); i++) {
            am[i] = factory.getMajorIsotope( ac.getAtom(i).getSymbol() ).getExactMass();
        }
        return(am);
    }

    private Point3d getCenterOfMass(Point3d[] p, double[] atwt) {
        double x = 0.;
        double y = 0.;
        double z = 0.;
        double totalmass = 0.;
        for (int i = 0; i  < p.length; i++) {
            x += atwt[i]*p[i].x;
            y += atwt[i]*p[i].y;
            z += atwt[i]*p[i].z;
            totalmass += atwt[i];
        }
        return( new Point3d(x/totalmass, y/totalmass, z/totalmass) );
    }
        

    /**
     * Sets up variables for the alignment algorithm.
     *
     * The algorithm allows for atom weighting and the default is 1.0 for all
     * atoms.
     *
     * @param al1 An array of {@link Atom} objects
     * @param al2 An array of {@link Atom} objects. This array will have its coordinates rotated
     *            so that the RMDS is minimzed to the coordinates of the first array
     * @throws CDKException if the number of Atom's are not the same in the two arrays
     */            
    public KabschAlignment(Atom[] al1, Atom[] al2) throws CDKException {
        if (al1.length != al2.length) {
            throw new CDKException("The Atom[]'s being aligned must have the same numebr of atoms");
        }
        this.npoint = al1.length;
        this.p1 = getPoint3dArray(al1);
        this.p2 = getPoint3dArray(al2);
        this.wts = new double[this.npoint];

        this.atwt1 = getAtomicMasses(al1);
        this.atwt2 = getAtomicMasses(al2);

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

        this.atwt1 = getAtomicMasses(al1);
        this.atwt2 = getAtomicMasses(al2);
    }

    /**
     * Sets up variables for the alignment algorithm.
     *
     * The algorithm allows for atom weighting and the default is 1.0 for all
     * atoms.
     *
     * @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
     * @throws CDKException if the number of atom's are not the same in the two AtomContainer's
     */            
    public KabschAlignment(IAtomContainer ac1, IAtomContainer ac2) throws CDKException {
        if (ac1.getAtomCount() != ac2.getAtomCount()) {
            throw new CDKException("The AtomContainer's being aligned must have the same 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] = 1.0;

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

    /**
     * Sets up variables for the alignment algorithm.