modelbuilder3d.java

化学图形处理软件

/*  $RCSfile$
 *  $Author: egonw $
 *  $Date: 2007-02-02 00:15:18 +0100 (Fri, 02 Feb 2007) $
 *  $Revision: 7850 $
 *
 *  Copyright (C) 2005-2007  Christian Hoppe <[e-mail]>
 *
 *  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.modeling.builder3d;

import java.io.IOException;
import java.util.HashMap;
import java.util.Hashtable;
import java.util.Iterator;
import java.util.Map;
import java.util.List;

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

import org.openscience.cdk.CDKConstants;
import org.openscience.cdk.exception.CDKException;
import org.openscience.cdk.geometry.GeometryToolsInternalCoordinates;
import org.openscience.cdk.graph.ConnectivityChecker;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IMolecule;
import org.openscience.cdk.interfaces.IRingSet;
import org.openscience.cdk.layout.AtomPlacer;
import org.openscience.cdk.ringsearch.RingPartitioner;
import org.openscience.cdk.tools.LoggingTool;
import org.openscience.cdk.tools.manipulator.RingSetManipulator;

/**
 *  The main class to generate the 3D coordinates of a molecule ModelBuilder3D.
 *  Its use looks like:
 *  <pre>
 *  ModelBuilder3D mb3d = new ModelBuilder3D();
 *  mb3d.setTemplateHandler();
 *  mb3d.setForceField("mm2");
 *  mb3d.setMolecule(molecule, false);
 *  mb3d.generate3DCoordinates();
 *  Molecule molecule = mb3d.getMolecule();
 *  </pre>
 *
 *  <p>Standing problems:
 *  <ul>
 *    <li>condensed ring systems which are unknown for the template class 
 *    <li>vdWaals clashes
 *    <li>stereochemistry
 *    <li>chains running through ring systems
 *  </ul>
 *
 * @author      cho
 * @author      steinbeck
 * @cdk.created 2004-09-07
 * @cdk.module  builder3d
 * @cdk.keyword 3D coordinates
 * @cdk.keyword coordinate generation, 3D
 * @cdk.bug     1241421
 * @cdk.bug     1315823
 * @cdk.bug     1458647
 */
public class ModelBuilder3D {

	private static Map memyselfandi = new HashMap();
	
	private TemplateHandler3D templateHandler = null;
		
	private Hashtable parameterSet = null;

	private final ForceFieldConfigurator ffc = new ForceFieldConfigurator();

	String forceFieldName = "mm2";
	
	private LoggingTool logger = new LoggingTool(ModelBuilder3D.class);
	
	/**
	 *  Constructor for the ModelBuilder3D object
	 *
	 *@param  molecule         Molecule
	 *@param  templateHandler  templateHandler Object
	 *@param  ffname           name of force field
	 */
	private ModelBuilder3D(TemplateHandler3D templateHandler, String ffname) throws CDKException {
		setTemplateHandler(templateHandler);
		setForceField(ffname);
	}

	public static ModelBuilder3D getInstance(TemplateHandler3D templateHandler, String ffname) throws CDKException {
		if (ffname == null || ffname.length() == 0) throw new CDKException("The given ffname is null or empty!");
		if (templateHandler == null) throw new CDKException("The given template handler is null!");
		
		String builderCode = templateHandler.getClass().getName()+ "#" + ffname;
		if (!memyselfandi.containsKey(builderCode)) {
			ModelBuilder3D builder = new ModelBuilder3D(
				templateHandler, ffname
			);
			memyselfandi.put(builderCode, builder);
			return builder;
		}
		return (ModelBuilder3D)memyselfandi.get(builderCode);
	}

	public static ModelBuilder3D getInstance() throws CDKException {
		return getInstance(TemplateHandler3D.getInstance(), "mm2");
	}

	/**
	 *  gives a list of possible force field types
	 *
	 *@return                the list
	 */
  public String[] getFfTypes(){
    return ffc.getFfTypes();
  }


	/**
	 *  Sets the forceField attribute of the ModelBuilder3D object
	 *
	 *@param  ffname  forceField name
	 * @throws CDKException 
	 */
	private void setForceField(String ffname) throws CDKException {
		if (ffname == null) {
			ffname = "mm2";
		}
		try {
			forceFieldName = ffname;
			ffc.setForceFieldConfigurator(ffname);
			parameterSet = ffc.getParameterSet();
		} catch (Exception ex1) {
			logger.error("Problem with ForceField configuration due to>" + ex1.getMessage());
			logger.debug(ex1);
			throw new CDKException("Problem with ForceField configuration due to>" + ex1.getMessage(), ex1);
		}
	}


	/**
	 * Generate 3D coordinates with force field information.
	 */
	public IMolecule generate3DCoordinates(IMolecule molecule, boolean clone) throws Exception {
		logger.debug("******** GENERATE COORDINATES ********");
		//CHECK FOR CONNECTIVITY!
		logger.debug("#atoms>"+molecule.getAtomCount());
		if (!ConnectivityChecker.isConnected(molecule)) {
			throw new CDKException("Molecule is NOT connected, could not layout.");
		}
		
		// setup helper classes
		AtomPlacer atomPlacer = new AtomPlacer();
		AtomPlacer3D ap3d = new AtomPlacer3D();
		AtomTetrahedralLigandPlacer3D atlp3d = new AtomTetrahedralLigandPlacer3D();
		ap3d.initilize(parameterSet);
		atlp3d.setParameterSet(parameterSet);
		
		if (clone) molecule = (IMolecule)molecule.clone();
		atomPlacer.setMolecule(molecule);
		
		if (ap3d.numberOfUnplacedHeavyAtoms(molecule) == 1) {
			logger.debug("Only one Heavy Atom");
			molecule.getAtom(0).setPoint3d(new Point3d(0.0, 0.0, 0.0));
			try {
				atlp3d.add3DCoordinatesForSinglyBondedLigands(molecule);
			} catch (Exception ex3) {
				logger.error("PlaceSubstitutensERROR: Cannot place substitutents due to:" + ex3.getMessage());
				logger.debug(ex3);
				throw new CDKException("PlaceSubstitutensERROR: Cannot place substitutents due to:" + ex3.getMessage(), ex3);
			}
			return molecule;
		}
		//Assing Atoms to Rings,Aliphatic and Atomtype
		IRingSet ringSetMolecule = ffc.assignAtomTyps(molecule);
		List ringSystems = null;
		IRingSet largestRingSet = null;
		double NumberOfRingAtoms = 0;

		if (ringSetMolecule.getAtomContainerCount() > 0) {
			if(templateHandler==null){
				throw new CDKException("You are trying to generate coordinates for a molecule with rings, but you have no template handler set. Please do setTemplateHandler() before generation!");
			}
			ringSystems = RingPartitioner.partitionRings(ringSetMolecule);
			largestRingSet = getLargestRingSet(ringSystems);
			IAtomContainer largestRingSetContainer = getAllInOneContainer(largestRingSet);
			NumberOfRingAtoms = (double)largestRingSetContainer.getAtomCount();
			templateHandler.mapTemplates(largestRingSetContainer, NumberOfRingAtoms);
			if (!checkAllRingAtomsHasCoordinates(largestRingSetContainer)) {
				throw new CDKException("RingAtomLayoutError: Not every ring atom is placed! Molecule cannot be layout.");
			}

			setAtomsToPlace(largestRingSetContainer);
			searchAndPlaceBranches(molecule, largestRingSetContainer, ap3d, atlp3d, atomPlacer);
			largestRingSet = null;
		} else {
			//logger.debug("****** Start of handling aliphatic molecule ******");
			IAtomContainer ac = null;

			ac = atomPlacer.getInitialLongestChain(molecule);
			setAtomsToUnVisited(molecule);
			setAtomsToUnPlaced(molecule);
			ap3d.placeAliphaticHeavyChain(molecule, ac);
			//ZMatrixApproach
			ap3d.zmatrixChainToCartesian(molecule, false);
			searchAndPlaceBranches(molecule, ac, ap3d, atlp3d, atomPlacer);
		}
		layoutMolecule(ringSystems, molecule, ap3d, atlp3d, atomPlacer);
		//logger.debug("******* PLACE SUBSTITUENTS ******");
		try {
			atlp3d.add3DCoordinatesForSinglyBondedLigands(molecule);
		} catch (Exception ex3) {
			logger.error("PlaceSubstitutensERROR: Cannot place substitutents due to:" + ex3.getMessage());
			logger.debug(ex3);
		}
		return molecule;
	}



	/**
	 *  Gets the ringSetOfAtom attribute of the ModelBuilder3D object
	 *
	 *@param  ringSystems  Description of the Parameter
	 *@param  atom         Description of the Parameter
	 *@return              The ringSetOfAtom value
	 */
	private IRingSet getRingSetOfAtom(List ringSystems, IAtom atom) {
		IRingSet ringSetOfAtom = null;
		for (int i = 0; i < ringSystems.size(); i++) {
			if (((IRingSet) ringSystems.get(i)).contains(atom)) {
				return (IRingSet) ringSystems.get(i);
			}
		}
		return ringSetOfAtom;
	}


	/**
	 *  Layout the molecule, starts with ring systems and than aliphatic chains
	 *
	 *@param  ringSetMolecule  ringSystems of the molecule
	 * @param atlp3d 
	 * @param atomPlacer 
	 *@exception  Exception    Description of the Exception
	 */
	private void layoutMolecule(List ringSetMolecule, IMolecule molecule, AtomPlacer3D ap3d, AtomTetrahedralLigandPlacer3D atlp3d, AtomPlacer atomPlacer) throws Exception {
		//logger.debug("****** LAYOUT MOLECULE MAIN *******");
		IAtomContainer ac = null;
		int safetyCounter = 0;
		IAtom atom = null;
		//Place rest Chains/Atoms
		do {
			safetyCounter++;
			atom = ap3d.getNextPlacedHeavyAtomWithUnplacedRingNeighbour(molecule);
			if (atom != null) {
				//logger.debug("layout RingSystem...");
				IAtom unplacedAtom = ap3d.getUnplacedRingHeavyAtom(molecule, atom);
				IRingSet ringSetA = getRingSetOfAtom(ringSetMolecule, unplacedAtom);
				IAtomContainer ringSetAContainer = getAllInOneContainer(ringSetA);
				templateHandler.mapTemplates(ringSetAContainer, (double)ringSetAContainer.getAtomCount());

				if (checkAllRingAtomsHasCoordinates(ringSetAContainer)) {
				} else {
					throw new IOException("RingAtomLayoutError: Not every ring atom is placed! Molecule cannot be layout.Sorry");
				}

				Point3d firstAtomOriginalCoord = unplacedAtom.getPoint3d();
				Point3d centerPlacedMolecule = ap3d.geometricCenterAllPlacedAtoms(molecule);

				setBranchAtom(molecule, unplacedAtom, atom, ap3d.getPlacedHeavyAtoms(molecule, atom), ap3d, atlp3d);
				layoutRingSystem(firstAtomOriginalCoord, unplacedAtom, ringSetA, centerPlacedMolecule, atom, ap3d);
				searchAndPlaceBranches(molecule, ringSetAContainer, ap3d, atlp3d, atomPlacer);
				//logger.debug("Ready layout Ring System");
				ringSetA = null;
				unplacedAtom = null;
				firstAtomOriginalCoord = null;
				centerPlacedMolecule = null;
			} else {
				//logger.debug("layout chains...");
				setAtomsToUnVisited(molecule);
				atom = ap3d.getNextPlacedHeavyAtomWithUnplacedAliphaticNeighbour(molecule);
				if (atom != null) {
					ac = new org.openscience.cdk.AtomContainer();
					ac.addAtom(atom);
					searchAndPlaceBranches(molecule, ac, ap3d, atlp3d, atomPlacer);
					ac = null;
				}
			}
		} while (!ap3d.allHeavyAtomsPlaced(molecule) || safetyCounter > molecule.getAtomCount());
	}


	/**
	 *  Layout the ring system, rotate and translate the template
	 *
	 *@param  originalCoord         coordinates of the placedRingAtom from the template 
	 *@param  placedRingAtom        placedRingAtom
	 *@param  ringSet               ring system which placedRingAtom is part of
	 *@param  centerPlacedMolecule  the geometric center of the already placed molecule
	 *@param  atomB                 placed neighbour atom of  placedRingAtom
	 * @param ap3d 
	 */
	private void layoutRingSystem(Point3d originalCoord, IAtom placedRingAtom, IRingSet ringSet, Point3d centerPlacedMolecule, IAtom atomB, AtomPlacer3D ap3d) {
		//logger.debug("****** Layout ring System ******");System.out.println(">around atom:"+molecule.getAtomNumber(placedRingAtom));
		IAtomContainer ac = getAllInOneContainer(ringSet);
		Point3d newCoord = placedRingAtom.getPoint3d();
		Vector3d axis = new Vector3d(atomB.getPoint3d().x - newCoord.x, atomB.getPoint3d().y - newCoord.y, atomB.getPoint3d().z - newCoord.z);
		translateStructure(originalCoord, newCoord, ac);
		//Rotate Ringsystem to farthest possible point
		Vector3d startAtomVector = new Vector3d(newCoord.x - atomB.getPoint3d().x, newCoord.y - atomB.getPoint3d().y, newCoord.z - atomB.getPoint3d().z);
		IAtom farthestAtom = ap3d.getFarthestAtom(placedRingAtom.getPoint3d(), ac);
		Vector3d farthestAtomVector = new Vector3d(farthestAtom.getPoint3d().x - newCoord.x, farthestAtom.getPoint3d().y - newCoord.y, farthestAtom.getPoint3d().z - newCoord.z);
		Vector3d n1 = new Vector3d();
		n1.cross(axis, farthestAtomVector);
		n1.normalize();
		double lengthFarthestAtomVector = farthestAtomVector.length();
		Vector3d farthestVector = new Vector3d(startAtomVector);