atomtools.java

化学图形处理软件

/* $RCSfile$
 * $Author: egonw $
 * $Date: 2007-01-04 18:46:10 +0100 (Thu, 04 Jan 2007) $
 * $Revision: 7636 $
 * 
 * Copyright (C) 2003-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.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * 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;

import java.util.List;

import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;

import org.openscience.cdk.Atom;
import org.openscience.cdk.AtomContainer;
import org.openscience.cdk.Bond;
import org.openscience.cdk.CDKConstants;
import org.openscience.cdk.interfaces.IAtom;

/**
 * A set of static utility classes for geometric calculations on Atoms.
 *
 * @author Peter Murray-Rust
 * @cdk.created 2003-06-14
 */
public class AtomTools {
    
    public final static double TETRAHEDRAL_ANGLE = 
        2.0 * Math.acos(1.0 / Math.sqrt(3.0));


    /** 
     * Generate coordinates for all atoms which are singly bonded and have
     * no coordinates. This is useful when hydrogens are present but have
     * no coords. It knows about C, O, N, S only and will give tetrahedral or 
     * trigonal geometry elsewhere. Bond lengths are computed from covalent radii
     * if available. Angles are tetrahedral or trigonal
     *
     * @param atomContainer the set of atoms involved
     *
     * @cdk.keyword coordinate calculation
     * @cdk.keyword 3D model
     */
    public static void add3DCoordinates1(AtomContainer atomContainer) {
            // atoms without coordinates
        AtomContainer noCoords = new org.openscience.cdk.AtomContainer();
        // get vector of possible referenceAtoms?
        AtomContainer refAtoms = new org.openscience.cdk.AtomContainer();
        for (int i = 0; i < atomContainer.getAtomCount(); i++) {
        	IAtom atom = atomContainer.getAtom(i);
            // is this atom without 3D coords, and has only one ligand?
            if (atom.getPoint3d() == null) {
            	java.util.List connectedAtoms = atomContainer.getConnectedAtomsList(atom);
                if (connectedAtoms.size() == 1) {
                	IAtom refAtom = (IAtom)connectedAtoms.get(0);;
                    if (refAtom.getPoint3d() != null) {
                        refAtoms.addAtom(refAtom);
                        // store atoms with no coords and ref atoms in a 
                        // single container
                        noCoords.addAtom(atom);
                        noCoords.addAtom(refAtom);
                        // bond is required to extract ligands
                        noCoords.addBond(new Bond(atom, refAtom, CDKConstants.BONDORDER_SINGLE));
                    }
                }
            }
        }
        // now add coordinates to ligands of reference atoms
        // use default length of 1.0, which can be adjusted later
        double length = 1.0;
        double angle = TETRAHEDRAL_ANGLE;
        for (int i = 0; i < refAtoms.getAtomCount(); i++) {
        	IAtom refAtom = refAtoms.getAtom(i);
        	java.util.List noCoordLigands = noCoords.getConnectedAtomsList(refAtom);
            int nLigands = noCoordLigands.size();
            int nwanted = nLigands;
            String elementType = refAtom.getSymbol();
            // try to deal with lone pairs on small hetero
            if (elementType.equals("N") ||
                elementType.equals("O") ||
                elementType.equals("S")) {
                nwanted = 3;
            }
            Point3d[] newPoints = calculate3DCoordinatesForLigands(
                atomContainer, refAtom, nwanted, length, angle);
            for (int j = 0; j < nLigands; j++) {
            	IAtom ligand = (IAtom)noCoordLigands.get(j);
                Point3d newPoint = rescaleBondLength(refAtom, ligand, newPoints[j]);
                ligand.setPoint3d(newPoint);
            }
        }
    }
    
    /** 
     * Rescales Point2 so that length 1-2 is sum of covalent radii.
     * if covalent radii cannot be found, use bond length of 1.0
     *
     * @param  atom1  stationary atom
     * @param  atom2  moveable atom
     * @param  point2 coordinates for atom 2
     * @return        new coords for atom 2
     */
    public static Point3d rescaleBondLength(
    		org.openscience.cdk.interfaces.IAtom atom1, org.openscience.cdk.interfaces.IAtom atom2, Point3d point2) {
        Point3d point1 = atom1.getPoint3d();
        double d1 = atom1.getCovalentRadius();
        double d2 = atom2.getCovalentRadius();
// in case we have no covalent radii, set to 1.0        
        double distance = (d1 < 0.1 || d2 < 0.1) ? 1.0 : 
            atom1.getCovalentRadius() + atom2.getCovalentRadius();
        Vector3d vect = new Vector3d(point2);
        vect.sub(point1);
        vect.normalize();
        vect.scale(distance);
        Point3d newPoint = new Point3d(point1);
        newPoint.add(vect);
        return newPoint;
    }
               
    /**
     * Adds 3D coordinates for singly-bonded ligands of a reference atom (A).
     * Initially designed for hydrogens. The ligands of refAtom are identified
     * and those with 3D coordinates used to generate the new points. (This
     * allows strucures with partially known 3D coordinates to be used, as when
     * groups are added.)
     * "Bent" and "non-planar" groups can be formed by taking a subset of the
     * calculated points. Thus R-NH2 could use 2 of the 3 points calculated 
     * from (1,iii)
     * nomenclature: A is point to which new ones are "attached". 
     *     A may have ligands B, C...
     *     B may have ligands J, K..
     *     points X1, X2... are returned
     * The cases (see individual routines, which use idealised geometry by default):
     * (0) zero ligands of refAtom. The resultant points are randomly oriented:
     *    (i) 1 points  required; +x,0,0
     *    (ii) 2 points: use +x,0,0 and -x,0,0
     *    (iii) 3 points: equilateral triangle in xy plane
     *    (iv) 4 points x,x,x, x,-x,-x, -x,x,-x, -x,-x,x
     * (1a) 1 ligand(B) of refAtom which itself has a ligand (J)
     *    (i) 1 points  required; vector along AB vector
     *    (ii) 2 points: 2 vectors in ABJ plane, staggered and eclipsed wrt J
     *    (iii) 3 points: 1 staggered wrt J, the others +- gauche wrt J
     * (1b) 1 ligand(B) of refAtom which has no other ligands. A random J is 
     * generated and (1a) applied
     * (2) 2 ligands(B, C) of refAtom A
     *    (i) 1 points  required; vector in ABC plane bisecting AB, AC. If ABC is
     *        linear, no points
     *    (ii) 2 points: 2 vectors at angle ang, whose resultant is 2i
     * (3) 3 ligands(B, C, D) of refAtom A
     *    (i) 1 points  required; if A, B, C, D coplanar, no points. 
     *       else vector is resultant of BA, CA, DA
     
     * fails if atom itself has no coordinates or >4 ligands
     *
     * @param atomContainer describing the ligands of refAtom. It could be the
     * whole molecule, or could be a selected subset of ligands
     * @param refAtom (A) to which new ligands coordinates could be added
     * @param length A-X length
     * @param angle B-A-X angle (used in certain cases)
     * @return Point3D[] points calculated. If request could not be fulfilled (e.g.
     * too many atoms, or strange geometry, returns empty array (zero length, 
     * not null)
     *
     * @cdk.keyword coordinate generation
     */
    public static Point3d[] calculate3DCoordinatesForLigands(
        AtomContainer atomContainer, org.openscience.cdk.interfaces.IAtom refAtom, int nwanted, 
        double length, double angle) {
        Point3d newPoints[] = new Point3d[0];
        Point3d aPoint = refAtom.getPoint3d();
        // get ligands
	    List connectedAtoms = atomContainer.getConnectedAtomsList(refAtom);
        if (connectedAtoms == null) {
            return newPoints;
        }
        int nligands = connectedAtoms.size();
        AtomContainer ligandsWithCoords    = new org.openscience.cdk.AtomContainer();
        for (int i = 0; i < nligands; i++) {
            Atom ligand = (Atom) connectedAtoms.get(i);
            if (ligand.getPoint3d() != null) {
                ligandsWithCoords.addAtom(ligand);
            }
        }
        int nwithCoords = ligandsWithCoords.getAtomCount();
// too many ligands at present        
        if (nwithCoords > 3) {
            return newPoints;
        }
        if (nwithCoords == 0) {
            newPoints = calculate3DCoordinates0(refAtom.getPoint3d(), nwanted, length);
        } else if (nwithCoords == 1) {
// ligand on A            
        	org.openscience.cdk.interfaces.IAtom bAtom = ligandsWithCoords.getAtom(0);
            connectedAtoms = ligandsWithCoords.getConnectedAtomsList(bAtom);
// does B have a ligand (other than A)            
            Atom jAtom = null;
            for (int i = 0; i < connectedAtoms.size(); i++) {
                Atom connectedAtom = (Atom) connectedAtoms.get(i);
                if (!connectedAtom.equals(refAtom)) {
                    jAtom = connectedAtom;
                    break;
                }
            }
            newPoints = calculate3DCoordinates1(aPoint, bAtom.getPoint3d(), (jAtom != null) ? jAtom.getPoint3d() : null, nwanted, length, angle);
        } else if (nwithCoords == 2) {
            Point3d bPoint = ligandsWithCoords.getAtom(0).getPoint3d();
            Point3d cPoint = ligandsWithCoords.getAtom(1).getPoint3d();
            newPoints = calculate3DCoordinates2(aPoint, bPoint, cPoint, nwanted, length, angle);
        } else if (nwithCoords == 3) {
            Point3d bPoint = ligandsWithCoords.getAtom(0).getPoint3d();
            Point3d cPoint = ligandsWithCoords.getAtom(1).getPoint3d();