[Cdk-commits] CVS: cdk/src/org/openscience/cdk/modeling/builder3d AtomPlacer3D.java,NONE,1.1 AtomTetrahedralLigandPlacer3D.java,NONE,1.1 ForceFieldConfigurator.java,NONE,1.1 MM2BasedParameterSetReader.java,NONE,1.1 MMFF94BasedParameterSetReader.java,NONE,1.1 MMFF94BasedAtomTypePattern.java,NONE,1.1 ModelBuilder3D.java,NONE,1.1 TemplateHandler3D.java,NONE,1.1

[Christian H. <ch...@us...>]

Update of /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/builder3d
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv4122/src/org/openscience/cdk/modeling/builder3d

Added Files:
	AtomPlacer3D.java AtomTetrahedralLigandPlacer3D.java 
	ForceFieldConfigurator.java MM2BasedParameterSetReader.java 
	MMFF94BasedParameterSetReader.java 
	MMFF94BasedAtomTypePattern.java ModelBuilder3D.java 
	TemplateHandler3D.java 
Log Message:
The 3D Model Builder

--- NEW FILE: AtomPlacer3D.java ---
/*
 *  $RCSfile: AtomPlacer3D.java,v $
 *  $Author: chhoppe $
 *  $Date: 2004/11/30 10:24:48 $
 *   *
 *  Copyright (C) 1997-2004  The Chemistry Development Kit (CDK) project
 *
 *  Contact: cdk...@li...
 *
 *  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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */
package org.openscience.cdk.modeling.builder3d;

import java.util.Hashtable;

import javax.vecmath.Point3d;
import javax.vecmath.Point2d;
import javax.vecmath.Vector3d;
import java.util.Vector;
import org.openscience.cdk.Atom;
import org.openscience.cdk.Bond;
import org.openscience.cdk.AtomContainer;
import org.openscience.cdk.CDKConstants;


/**
 * Place aliphatic chains with Z matrix method.
 *
 * @author         chhoppe
 * @cdk.keyword    AtomPlacer3D
 * @cdk.created    2004-10-8
 */
public class AtomPlacer3D {

	private Hashtable pSet = null;
	private double[] distances;
	private int[] first_atoms = null;
	private double[] angles = null;
	private int[] second_atoms = null;
	private double[] dihedrals = null;
	private int[] third_atoms = null;
	private final static double DIHEDRAL_EXTENDED_CHAIN = (180.0 / 180) * Math.PI;
	private final static double DIHEDRAL_BRANCHED_CHAIN = 0.0;
	private final static double DEFAULT_BOND_LENGTH = 1.5;
	private final static double DEFAULT_ANGLE = 109.471;
	
	AtomPlacer3D(){}
	
	/**
	 *  Initialize the atomPlacer class
	 * 
	 * @param  parameterSet  Force Field parameter as Hashtable
	 */
	public void initilize( Hashtable parameterSet) {
		pSet = parameterSet;
	}

	/**
	 *  Count and find first heavy atom(s) (non Hydrogens) in a chain
	 *
	 * @param  chain  chain to be searched
	 * @return        the atom number of the first heavy atom the number of heavy atoms in the chain
	 */
	public int[] findHeavyAtomsInChain(AtomContainer molecule, AtomContainer chain) {
		int[] heavy = {-1, -1};
		int hc = 0;
		for (int i = 0; i < chain.getAtomCount(); i++) {
			if (!(chain.getAtomAt(i).getSymbol()).equals("H")) {
				if (heavy[0] < 0) {
					heavy[0] = molecule.getAtomNumber(chain.getAtomAt(i));
				}
				hc++;
			}
		}
		heavy[1] = hc;
		return heavy;
	}


	/**
	 *  Mark all atoms in chain as placed (CDKConstant ISPLACED)
	 *
	 * @param  ac  chain
	 * @return     chain all atoms marked as placed
	 */
	public AtomContainer markPlaced(AtomContainer ac) {
		for (int i = 0; i < ac.getAtomCount(); i++) {
			ac.getAtomAt(i).setFlag(CDKConstants.ISPLACED, true);
		}
		return ac;
	}

	
	/**
	 *  Method assigns 3Dcoordinates to the heavy atoms in an aliphatic chain
	 *
	 * @param  chain          the atoms to be assigned, must be connected
	 * @exception  Exception  Description of the Exception
	 */
	public void placeAliphaticHeavyChain(AtomContainer molecule, AtomContainer chain) throws Exception {
		//System.out.println("******** Place aliphatic Chain *********");
		int[] first = new int[2];
		int counter = 1;
		int nextAtomNr = 0;
		String ID1 = "";
		String ID2 = "";
		String ID3 = "";
		first = findHeavyAtomsInChain(molecule,chain);
		distances = new double[first[1]];
		first_atoms = new int[first[1]];
		angles = new double[first[1]];
		second_atoms = new int[first[1]];
		dihedrals = new double[first[1]];
		third_atoms = new int[first[1]];
		first_atoms[0] = first[0];
		molecule.getAtomAt(first_atoms[0]).setFlag(CDKConstants.VISITED, true);

		for (int i = 0; i < chain.getAtomCount(); i++) {
			if (!(chain.getAtomAt(i).getSymbol()).equals("H") &
					!chain.getAtomAt(i).getFlag(CDKConstants.VISITED)) {
				//System.out.print("Counter:" + counter);
				nextAtomNr = molecule.getAtomNumber(chain.getAtomAt(i));
				ID2 = molecule.getAtomAt(first_atoms[counter - 1]).getID();
				ID1 = molecule.getAtomAt(nextAtomNr).getID();
				distances[counter] = getBondLengthValue(ID1, ID2);
				//System.out.print(" Distance:" + distances[counter]);
				first_atoms[counter] = nextAtomNr;
				second_atoms[counter] = first_atoms[counter - 1];
				if (counter > 1) {
					ID3 = molecule.getAtomAt(first_atoms[counter - 2]).getID();
					angles[counter] = getAngleValue(ID1, ID2, ID3);
					third_atoms[counter] = first_atoms[counter - 2];
					//System.out.println(" Angle:" + angles[counter]);
				} else {
					angles[counter] = -1;
					third_atoms[counter] = -1;
				}
				if (counter > 2) {
					//double bond
					try{
						if (getDoubleBondConfiguration2D( molecule.getBond(molecule.getAtomAt(first_atoms[counter-1]),molecule.getAtomAt(first_atoms[counter-2])),
										(molecule.getAtomAt(first_atoms[counter])).getPoint2d(),(molecule.getAtomAt(first_atoms[counter-1])).getPoint2d(),
										(molecule.getAtomAt(first_atoms[counter-2])).getPoint2d(),(molecule.getAtomAt(first_atoms[counter-3])).getPoint2d())
							==5){
							dihedrals[counter] = DIHEDRAL_BRANCHED_CHAIN;
						}else{ dihedrals[counter] = DIHEDRAL_EXTENDED_CHAIN;}
					}catch(Exception ex1){
						dihedrals[counter] = DIHEDRAL_EXTENDED_CHAIN;
					}
				} else {
					dihedrals[counter] = -1;
				}
				counter++;
			}
		}
	}

	/**
	 * Takes the given Z Matrix coordinates and converts them to cartesian coordinates.
	 * The first Atom end up in the origin, the second on on the x axis, and the third
	 * one in the XY plane. The rest is added by applying the Zmatrix distances, angles
	 * and dihedrals. Assign coordinates directly to the atoms.
	 *
	 * @param  flag_branched  marks branched chain
	 * @author                egonw,cho
	 */

	public void zmatrixChainToCartesian(AtomContainer molecule, boolean flag_branched) {
		Point3d result = null;
		for (int index = 0; index < distances.length; index++) {
			if (index == 0) {
				result = new Point3d(0d, 0d, 0d);
			} else if (index == 1) {
				result = new Point3d(distances[1], 0d, 0d);
			} else if (index == 2) {
				result = new Point3d(-Math.cos((angles[2] / 180) * Math.PI) * distances[2] + distances[1],
						Math.sin((angles[2] / 180) * Math.PI) * distances[2],
						0d);
			} else {
				Vector3d cd = new Vector3d();
				cd.sub((molecule.getAtomAt(third_atoms[index])).getPoint3d(), (molecule.getAtomAt(second_atoms[index])).getPoint3d());

				Vector3d bc = new Vector3d();
				bc.sub(molecule.getAtomAt(second_atoms[index]).getPoint3d(), molecule.getAtomAt(first_atoms[index - 3]).getPoint3d());

				Vector3d n1 = new Vector3d();
				n1.cross(cd, bc);
				n1.normalize();

				Vector3d n2 = null;
				if (index == 3 & flag_branched) {
					n2 = AtomTetrahedralLigandPlacer3D.rotate(n1, bc, DIHEDRAL_BRANCHED_CHAIN);
				} else {
					n2 = AtomTetrahedralLigandPlacer3D.rotate(n1, bc, dihedrals[index]);
				}
				n2.normalize();

				Vector3d ba = new Vector3d();
				if (index == 3 & flag_branched) {
					ba = AtomTetrahedralLigandPlacer3D.rotate(cd, n2, (-angles[index] / 180) * Math.PI);
					ba = AtomTetrahedralLigandPlacer3D.rotate(ba, cd, (-angles[index] / 180) * Math.PI);
				} else {
					ba = AtomTetrahedralLigandPlacer3D.rotate(cd, n2, (-angles[index] / 180) * Math.PI);
				}

				ba.normalize();

				Vector3d ban = new Vector3d(ba);
				ban.scale(distances[index]);

				result = new Point3d();
				result.add(molecule.getAtomAt(first_atoms[index - 1]).getPoint3d(), ban);
			}

			if ((molecule.getAtomAt(first_atoms[index]).getPoint3d() == null || !(molecule.getAtomAt(first_atoms[index])).getFlag(CDKConstants.ISPLACED))
					 && !(molecule.getAtomAt(first_atoms[index])).getFlag(CDKConstants.ISINRING)
					 && !(molecule.getAtomAt(first_atoms[index])).getSymbol().equals("H")) {
				molecule.getAtomAt(first_atoms[index]).setPoint3d(result);
				molecule.getAtomAt(first_atoms[index]).setFlag(CDKConstants.ISPLACED, true);
			}
		}
	}


	
	private int getDoubleBondConfiguration2D(Bond bond,Point2d a, Point2d b,Point2d c,Point2d d) throws Exception{
		if (bond.getOrder()<1.5 || bond.getOrder()>2){
			return 0;
		}
		Point2d cb=new Point2d(c.x-b.x,c.y-b.y);
		Point2d xT=new Point2d(cb.x-1,cb.y);
		a.y=a.y-b.y-xT.y;
		d.y=d.y-b.y-xT.y;
		if ((a.y>0 && d.y>0)||(a.y<0 && d.y<0)){
			return 5;
		}else {return 6;}
	}

	/**
	 *  Gets the distanceValue attribute of the parameter set
	 *
	 * @param  id1            atom1 id
	 * @param  id2            atom2 id
	 * @return                The distanceValue value from the force field parameter set
	 * @exception  Exception  Description of the Exception
	 */
	public double getBondLengthValue(String id1, String id2) throws Exception {
		String dkey = "";
		if (pSet.containsKey(("bond" + id1 + ";" + id2))) {
			dkey="bond" + id1 + ";" + id2;
		}else if (pSet.containsKey(("bond" + id2 + ";" + id1))) {
			dkey = "bond" + id2 + ";" + id1;
		} else {
			System.out.println("KEYError:Unknown distance key in pSet: " + id2 + " ;" + id1+" take default bon length:"+DEFAULT_BOND_LENGTH);
			return DEFAULT_BOND_LENGTH;
		}
		return ((Double) (((Vector) pSet.get(dkey)).get(0))).doubleValue();
	}

	/**
	 *  Gets the angleKey attribute of the AtomPlacer3D object
	 *
	 * @param  id1            Description of the Parameter
	 * @param  id2            Description of the Parameter
	 * @param  id3            Description of the Parameter
	 * @return                The angleKey value
	 * @exception  Exception  Description of the Exception
	 */
	public double getAngleValue(String id1, String id2, String id3) throws Exception {
		String akey = "";
		if (pSet.containsKey(("angle" + id1 + ";" + id2 + ";" + id3))) {
			akey = "angle" + id1 + ";" + id2 + ";" + id3;
		} else if (pSet.containsKey(("angle" + id3 + ";" + id2 + ";" + id1))) {
			akey = "angle" + id3 + ";" + id2 + ";" + id1;
		} else if (pSet.containsKey(("angle" + id2 + ";" + id1 + ";" + id3))) {
			akey = "angle" + id2 + ";" + id1 + ";" + id3;
		} else if (pSet.containsKey(("angle" + id1 + ";" + id3 + ";" + id2))) {
			akey = "angle" + id1 + ";" + id3 + ";" + id2;
		} else if (pSet.containsKey(("angle" + id3 + ";" + id1 + ";" + id2))) {
			akey = "angle" + id3 + ";" + id1 + ";" + id2;
		} else if (pSet.containsKey(("angle" + id2 + ";" + id3 + ";" + id1))) {
			akey = "angle" + id2 + ";" + id3 + ";" + id1;
		} else {
			System.out.println("KEYErrorAngle:Unknown angle key in pSet: " +id2 + " ; " + id3 + " ; " + id1+" take default angle:"+DEFAULT_ANGLE);
			return DEFAULT_ANGLE;
		}
		return ((Double) (((Vector) pSet.get(akey)).get(0))).doubleValue();
	}


	/**
	 *  Gets the nextUnplacedHeavyAtomWithAliphaticPlacedNeighbour from an atom container or molecule
	 *
	 * @return    The nextUnplacedHeavyAtomWithAliphaticPlacedNeighbour value
	 * @author    steinbeck,cho
	 */
	public Atom getNextUnplacedHeavyAtomWithAliphaticPlacedNeighbour(AtomContainer molecule) {
		Bond[] bonds = molecule.getBonds();
		for (int i = 0; i < bonds.length; i++) {
			Atom[] atoms = bonds[i].getAtoms();
			if (atoms[0].getFlag(CDKConstants.ISPLACED) & !(atoms[1].getFlag(CDKConstants.ISPLACED))) {
				if (atoms[1].getFlag(CDKConstants.ISALIPHATIC) & !atoms[1].getSymbol().equals("H")) {
					return atoms[1];
				}
			}
			if (atoms[1].getFlag(CDKConstants.ISPLACED) & !(atoms[0].getFlag(CDKConstants.ISPLACED))) {
				if (atoms[0].getFlag(CDKConstants.ISALIPHATIC) & !atoms[0].getSymbol().equals("H")) {
					return atoms[0];
				}
			}
		}
		return null;
	}

	/**
	 *  Gets the nextPlacedHeavyAtomWithAliphaticPlacedNeigbor from an atom container or molecule
	 *
	 * @return    The nextUnplacedHeavyAtomWithUnplacedAliphaticNeigbor
	 * @author    steinbeck,cho
	 */
	public Atom getNextPlacedHeavyAtomWithUnplacedAliphaticNeighbour(AtomContainer molecule) {
		Bond[] bonds = molecule.getBonds();
		for (int i = 0; i < bonds.length; i++) {
			Atom[] atoms = bonds[i].getAtoms();
			if (atoms[0].getFlag(CDKConstants.ISPLACED) & !(atoms[1].getFlag(CDKConstants.ISPLACED))) {
				if (atoms[1].getFlag(CDKConstants.ISALIPHATIC) & !atoms[0].getSymbol().equals("H") & !atoms[1].getSymbol().equals("H")) {
					return atoms[0];
				}
			}
			if (atoms[1].getFlag(CDKConstants.ISPLACED) & !(atoms[0].getFlag(CDKConstants.ISPLACED))) {
				if (atoms[0].getFlag(CDKConstants.ISALIPHATIC) & !atoms[1].getSymbol().equals("H") & !atoms[0].getSymbol().equals("H")) {
					return atoms[1];
				}
			}
		}
		return null;
	}

	/**
	 *  Gets the nextPlacedHeavyAtomWithUnplacedRingNeighbour attribute of the AtomPlacer3D object
	 *
	 * @return    The nextPlacedHeavyAtomWithUnplacedRingNeighbour value
	 */
	public Atom getNextPlacedHeavyAtomWithUnplacedRingNeighbour(AtomContainer molecule) {
		Bond[] bonds = molecule.getBonds();
		for (int i = 0; i < bonds.length; i++) {
			Atom[] atoms = bonds[i].getAtoms();
			if (atoms[0].getFlag(CDKConstants.ISPLACED) & !(atoms[1].getFlag(CDKConstants.ISPLACED))) {
				if (atoms[1].getFlag(CDKConstants.ISINRING) & !atoms[0].getSymbol().equals("H") & !atoms[1].getSymbol().equals("H")) {
					return atoms[0];
				}
			}
			if (atoms[1].getFlag(CDKConstants.ISPLACED) & !(atoms[0].getFlag(CDKConstants.ISPLACED))) {
				if (atoms[0].getFlag(CDKConstants.ISINRING) & !atoms[1].getSymbol().equals("H") & !atoms[0].getSymbol().equals("H")) {
					return atoms[1];
				}
			}
		}
		return null;
	}

	/**
	 *  Gets the farthestAtom attribute of the AtomPlacer3D object
	 *
	 * @param  refAtomPoint  Description of the Parameter
	 * @param  ac            Description of the Parameter
	 * @return               The farthestAtom value
	 */
	public Atom getFarthestAtom(Point3d refAtomPoint, AtomContainer ac) {
		double distance = 0;
		Atom atom = null;
		for (int i = 0; i < ac.getAtomCount(); i++) {
			if (ac.getAtomAt(i).getPoint3d() != null) {
				if (Math.abs(refAtomPoint.distance(ac.getAtomAt(i).getPoint3d())) > distance) {
					atom = ac.getAtomAt(i);
					distance = Math.abs(refAtomPoint.distance(ac.getAtomAt(i).getPoint3d()));
				}
			}
		}
		return atom;
	}

	/**
	 *  Gets the unplacedRingHeavyAtom attribute of the AtomPlacer3D object
	 *
	 * @param  atom  Description of the Parameter
	 * @return       The unplacedRingHeavyAtom value
	 */
	public Atom getUnplacedRingHeavyAtom(AtomContainer molecule, Atom atom) {

		Bond[] bonds = molecule.getConnectedBonds(atom);
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atom);
			if (!connectedAtom.getFlag(CDKConstants.ISPLACED) && !connectedAtom.getSymbol().equals("H") && connectedAtom.getFlag(CDKConstants.ISINRING)) {
				return connectedAtom;
			}
		}
		return connectedAtom;
	}

	/**
	 *  Calculates the geometric center of all placed atoms in the atomcontainer
	 *
	 * @return    Point3d the geometric center
	 */
	public Point3d geometricCenterAllPlacedAtoms(AtomContainer molecule) {
		AtomContainer allPlacedAtoms = getAllPlacedAtoms(molecule);
		return allPlacedAtoms.get3DCenter();
	}

	/**
	 *  Returns an unplaced atom connected to a given atom
	 *
	 * @param  atom  The Atom whose unplaced bonding partners are to be returned
	 * @return       an unplaced heavy atom connected to a given atom
	 * @author       steinbeck,cho
	 */
	public Atom getUnplacedHeavyAtom(AtomContainer molecule, Atom atom) {

		Bond[] bonds = molecule.getConnectedBonds(atom);
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atom);
			if (!connectedAtom.getFlag(CDKConstants.ISPLACED) & !connectedAtom.getSymbol().equals("H")) {
				return connectedAtom;
			}
		}
		return connectedAtom;
	}

	/**
	 *  Returns a placed atom connected to a given atom
	 *
	 * @param  atom  The Atom whose placed bonding partners are to be returned
	 * @return       a placed heavy atom connected to a given atom
	 * @author       steinbeck
	 */
	public Atom getPlacedHeavyAtom(AtomContainer molecule, Atom atom) {

		Bond[] bonds = molecule.getConnectedBonds(atom);
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atom);
			if (connectedAtom.getFlag(CDKConstants.ISPLACED) & !connectedAtom.getSymbol().equals("H")) {
				return connectedAtom;
			}
		}
		return connectedAtom;
	}

	/**
	 *  Gets the placedAtom attribute of the AtomPlacer3D object
	 *
	 * @param  atomA  Description of the Parameter
	 * @param  atomB  Description of the Parameter
	 * @return        The placedAtom value
	 */
	public Atom getPlacedAtom(AtomContainer molecule, Atom atomA, Atom atomB) {

		Bond[] bonds = molecule.getConnectedBonds(atomA);
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atomA);
			if (connectedAtom.getFlag(CDKConstants.ISPLACED) && connectedAtom != atomB) {
				return connectedAtom;
			}
		}
		return connectedAtom;
	}

	/**
	 *  Gets the placedHeavyAtom attribute of the AtomPlacer3D object
	 *
	 * @param  atomA  Description of the Parameter
	 * @param  atomB  Description of the Parameter
	 * @return        The placedHeavyAtom value
	 */
	public Atom getPlacedHeavyAtom(AtomContainer molecule, Atom atomA, Atom atomB) {
		Bond[] bonds = molecule.getConnectedBonds(atomA);
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atomA);
			if (connectedAtom.getFlag(CDKConstants.ISPLACED) && !connectedAtom.getSymbol().equals("H")
					 && connectedAtom != atomB) {
				return connectedAtom;
			}
		}
		return connectedAtom;
	}

	/**
	 *  Gets the placedHeavyAtoms attribute of the AtomPlacer3D object
	 *
	 * @param  atom  Description of the Parameter
	 * @return       The placedHeavyAtoms value
	 */
	public AtomContainer getPlacedHeavyAtoms(AtomContainer molecule, Atom atom) {

		Bond[] bonds = molecule.getConnectedBonds(atom);
		AtomContainer connectedAtoms = new AtomContainer();
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atom);
			if (connectedAtom.getFlag(CDKConstants.ISPLACED) & !(connectedAtom.getSymbol().equals("H"))) {
				connectedAtoms.addAtom(connectedAtom);
			}
		}
		return connectedAtoms;
	}

	/**
	 *  Gets all unplacedAtoms around a atom
	 *
	 * @param  molecule AtomContainer
	 * @param  atom  atom: center of search
	 * @return       AtomContainer of unplaced atoms
	 */
	public AtomContainer getUnplacedAtoms(AtomContainer molecule, Atom atom) {

		Bond[] bonds = molecule.getConnectedBonds(atom);
		AtomContainer connectedAtoms = new AtomContainer();
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atom);
			if (!connectedAtom.getFlag(CDKConstants.ISPLACED)) {
				connectedAtoms.addAtom(connectedAtom);
			}
		}
		return connectedAtoms;
	}
	
	/**
	 *  Gets the placedAtoms attribute of the AtomPlacer3D object
	 *
	 * @param  atom  Description of the Parameter
	 * @return       The placedAtoms value
	 */
	public int numberOfUnplacedHeavyAtoms(AtomContainer ac) {
		int nUnplacedHeavyAtoms=0;
		for (int i = 0; i < ac.getAtomCount(); i++) {
			if (!ac.getAtomAt(i).getFlag(CDKConstants.ISPLACED) && !ac.getAtomAt(i).equals("H")) {
				nUnplacedHeavyAtoms+=1;
			}
		}
		return nUnplacedHeavyAtoms;
	}
	
	/**
	 *  Gets the allPlacedAtoms attribute of the AtomPlacer3D object
	 *
	 * @return    The allPlacedAtoms value
	 */
	public AtomContainer getAllPlacedAtoms(AtomContainer molecule) {
		AtomContainer placedAtoms = new AtomContainer();
		for (int i = 0; i < molecule.getAtomCount(); i++) {
			if (molecule.getAtomAt(i).getFlag(CDKConstants.ISPLACED)) {
				placedAtoms.addAtom(molecule.getAtomAt(i));
			}
		}
		return placedAtoms;
	}

	/**
	 *  True is all the atoms in the given AtomContainer have been placed
	 *
	 * @param  ac  The AtomContainer to be searched
	 * @return     True is all the atoms in the given AtomContainer have been placed
	 */
	public boolean allHeavyAtomsPlaced(AtomContainer ac) {
		for (int i = 0; i < ac.getAtomCount(); i++) {
			if (!ac.getAtomAt(i).getFlag(CDKConstants.ISPLACED) & !(ac.getAtomAt(i).getSymbol().equals("H"))) {
				return false;
			}
		}
		return true;
	}
}



--- NEW FILE: AtomTetrahedralLigandPlacer3D.java ---
/*
 *  $RCSfile: AtomTetrahedralLigandPlacer3D.java,v $
 *  $Author: egonw
 *  $Date: 2004/11/30 10:24:48 $
 *  $ (chhoppe) add force field parameters and some methods
 *  $Revision: 1.1 $
 *
 *  Copyright (C) 2003-2004  The Chemistry Development Kit (CDK) project
 *
 *  Contact: cdk...@li...
 *
 *  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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
package org.openscience.cdk.modeling.builder3d;

import java.util.Vector;
import java.util.Hashtable;
import java.io.IOException;

import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
import javax.vecmath.Matrix3d;
import javax.vecmath.AxisAngle4d;

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

/**
 *  A set of static utility classes for geometric calculations on Atoms.
 *
 *@author         Peter Murray-Rust,chhoppe,egonw
 *@cdk.created    2004-09-07
 */
public class AtomTetrahedralLigandPlacer3D {

	private Hashtable pSet = null;
	private final double DEFAULT_BOND_LENGTH_H = 1.0;
	private final double DEFAULT_BOND_LENGTH_HA = 1.3;

	/**
	 *  Description of the Field
	 */
	public final double TETRAHEDRAL_ANGLE =
			2.0 * Math.acos(1.0 / Math.sqrt(3.0));

	private final double SP2_ANGLE = 120 * Math.PI / 180;
	private final double SP_ANGLE = Math.PI;

	final static Vector3d XV = new Vector3d(1.0, 0.0, 0.0);
	final static Vector3d YV = new Vector3d(0.0, 1.0, 0.0);


	/**
	 *  Constructor for the AtomTetrahedralLigandPlacer3D object
	 */
	AtomTetrahedralLigandPlacer3D() { }


	/**
	 *  Constructor for the setParameterSet object
	 *
	 *@param  moleculeParameter  Description of the Parameter
	 */
	public void setParameterSet(Hashtable moleculeParameter) {
		pSet = moleculeParameter;
	}


	/**
	 *  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 or taken
	 *  out of a paramter set if available. Angles are tetrahedral or trigonal
	 *
	 *@param  atomContainer  the set of atoms involved
	 *@exception  Exception  Description of the Exception
	 *@cdk.keyword           coordinate calculation
	 *@cdk.keyword           3D model
	 */
	public void add3DCoordinatesForSinglyBondedLigands(AtomContainer atomContainer) throws Exception {
		AtomContainer noCoords = new AtomContainer();
		AtomContainer withCoords = new AtomContainer();
		Atom refAtom = null;
		Atom atomC = null;
		int nwanted = 0;
		for (int i = 0; i < atomContainer.getAtomCount(); i++) {
			refAtom = atomContainer.getAtomAt(i);
			if (!refAtom.getSymbol().equals("H") && hasUnsetNeighbour(refAtom, atomContainer)) {
				noCoords = getUnsetAtomsInAtomContainer(refAtom, atomContainer);
				withCoords = getPlacedAtomsInAtomContainer(refAtom, atomContainer);
				if (withCoords.getAtomCount() > 0) {
					atomC = getPlacedHeavyAtomInAtomContainer(withCoords.getAtomAt(0), refAtom, atomContainer);
				}
				if (refAtom.getFormalNeighbourCount() == 0 && refAtom.getSymbol().equals("C")) {
					nwanted = noCoords.getAtomCount();
				} else if (refAtom.getFormalNeighbourCount() == 0 && !refAtom.getSymbol().equals("C")) {
					nwanted = 4;
				} else {
					nwanted = refAtom.getFormalNeighbourCount() - withCoords.getAtomCount();
				}
				Point3d[] newPoints = get3DCoordinatesForLigands(refAtom,
						noCoords, withCoords, atomC, nwanted, DEFAULT_BOND_LENGTH_H, -1);

				for (int j = 0; j < noCoords.getAtomCount(); j++) {
					Atom ligand = noCoords.getAtomAt(j);
					Point3d newPoint = rescaleBondLength(refAtom, ligand, newPoints[j]);
					ligand.setPoint3d(newPoint);
					ligand.setFlag(CDKConstants.ISPLACED, true);
				}

				noCoords.removeAllElements();
				withCoords.removeAllElements();
			}
		}
	}


	/**
	 *  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
	 *@exception  Exception  Description of the Exception
	 */
	public Point3d rescaleBondLength(Atom atom1, Atom atom2, Point3d point2) throws Exception {
		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();
		if (pSet != null) {
			distance = getDistanceValue(atom1.getID(), atom2.getID());
		}
		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  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)
	 *@param  nwanted        Description of the Parameter
	 *@param  noCoords       Description of the Parameter
	 *@param  withCoords     Description of the Parameter
	 *@param  atomC          Description of the Parameter
	 *@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)
	 *@exception  Exception  Description of the Exception
	 *@cdk.keyword           coordinate generation
	 */

	public Point3d[] get3DCoordinatesForLigands(Atom refAtom,
			AtomContainer noCoords, AtomContainer withCoords, Atom atomC, int nwanted, double length, double angle) throws Exception {
		Point3d newPoints[] = new Point3d[1];

		if (noCoords.getAtomCount() == 0 && withCoords.getAtomCount() == 0) {
			return newPoints;
		}

		// too many ligands at present
		if (withCoords.getAtomCount() > 3) {
			return newPoints;
		}

		if (refAtom.getFormalNeighbourCount() == 1 || refAtom.getMaxBondOrder() > 4) {
		} else if (refAtom.getFormalNeighbourCount() == 2 || refAtom.getMaxBondOrder() == 3) {
			//sp
			newPoints[0] = get3DCoordinatesForSPLigands(refAtom, withCoords, length, angle);
		} else if (refAtom.getFormalNeighbourCount() == 3 || (refAtom.getMaxBondOrder() > 1 && refAtom.getMaxBondOrder() < 3)) {
			//sp2
			try {
				newPoints = get3DCoordinatesForSP2Ligands(refAtom, noCoords, withCoords, atomC, length, angle);
			} catch (Exception ex1) {
				System.out.println("Get3DCoordinatesForLigandsERROR: Cannot place SP2 Ligands due to:" + ex1.toString());
				throw new IOException("Cannot place sp2 substituents");
			}

		} else {
			//sp3
			try {
				newPoints = get3DCoordinatesForSP3Ligands(refAtom, noCoords, withCoords, atomC, nwanted, length, angle);
			} catch (Exception ex1) {
				System.out.println("Get3DCoordinatesForLigandsERROR: Cannot place SP3 Ligands due to:" + ex1.toString());
				throw new IOException("Cannot place sp3 substituents");
			}
		}
		//System.out.println("...Ready "+newPoints.length+" "+newPoints[0].toString());
		return newPoints;
	}


	/**
	 *  Description of the Method
	 *
	 *@param  refAtom     Description of the Parameter
	 *@param  withCoords  Description of the Parameter
	 *@param  length      Description of the Parameter
	 *@param  angle       Description of the Parameter
	 *@return             Description of the Return Value
	 */
	public Point3d get3DCoordinatesForSPLigands(Atom refAtom, AtomContainer withCoords, double length, double angle) {
		//System.out.println(" SP Ligands start "+refAtom.getPoint3d()+" "+(withCoords.getAtomAt(0)).getPoint3d());
		Vector3d ca = new Vector3d(refAtom.getPoint3d());
		ca.sub((withCoords.getAtomAt(0)).getPoint3d());
		ca.normalize();
		ca.scale(length);
		Point3d newPoint = new Point3d(refAtom.getPoint3d());
		newPoint.add(ca);
		return newPoint;
	}


	/**
	 *  Main method for the calculation of the ligand coordinates for sp2 atoms.
	 *  Decides if one or two coordinates should be created
	 *
	 *@param  refAtom            central atom (Atom)
	 *@param  noCoords           Description of the Parameter
	 *@param  withCoords         Description of the Parameter
	 *@param  atomC              Description of the Parameter
	 *@param  length             Description of the Parameter
	 *@param  angle              Description of the Parameter
	 *@return                    coordinates as Points3d []
	 */
	public Point3d[] get3DCoordinatesForSP2Ligands(Atom refAtom, AtomContainer noCoords, AtomContainer withCoords, Atom atomC, double length, double angle) {
		//System.out.println(" SP2 Ligands start");
		Point3d newPoints[] = new Point3d[1];
		if (angle < 0) {
			angle = SP2_ANGLE;
		}
		if (withCoords.getAtomCount() >= 2) {
			//System.out.println("Wanted:1 "+noCoords.getAtomCount());
			newPoints[0] = calculate3DCoordinatesSP2_1(refAtom.getPoint3d(), (withCoords.getAtomAt(0)).getPoint3d(),
					(withCoords.getAtomAt(1)).getPoint3d(), length, -1 * angle);

		} else if (withCoords.getAtomCount() <= 1) {
			//System.out.println("NoCoords 2:"+noCoords.getAtomCount());
			newPoints = calculate3DCoordinatesSP2_2(refAtom.getPoint3d(), (withCoords.getAtomAt(0)).getPoint3d(),
					(atomC != null) ? atomC.getPoint3d() : null, length, angle);
		}
		//System.out.println("Ready SP2");
		return newPoints;
	}


	/**
	 *  Main method for the calculation of the ligand coordinates for sp3 atoms.
	 *  Decides how many coordinates should be created
	 *
	 *@param  refAtom            central atom (Atom)
	 *@param  nwanted            how many ligands should be created
	 *@param  length             bond length
	 *@param  angle              angle in a B-A-(X) system; a=central atom;
	 *      x=ligand with unknown coordinates
	 *@param  noCoords           Description of the Parameter
	 *@param  withCoords         Description of the Parameter
	 *@param  atomC              Description of the Parameter
	 *@return                    Description of the Return Value
	 */
	public Point3d[] get3DCoordinatesForSP3Ligands(Atom refAtom, AtomContainer noCoords, AtomContainer withCoords, Atom atomC, int nwanted, double length, double angle) {
		//System.out.print("SP3 Ligands start ");
		Point3d newPoints[] = new Point3d[0];
		Point3d aPoint = refAtom.getPoint3d();
		int nwithCoords = withCoords.getAtomCount();
		if (angle < 0) {
			angle = TETRAHEDRAL_ANGLE;
		}
		if (nwithCoords == 0) {
			newPoints = calculate3DCoordinates0(refAtom.getPoint3d(), nwanted, length);
		} else if (nwithCoords == 1) {
			newPoints = calculate3DCoordinates1(aPoint, (withCoords.getAtomAt(0)).getPoint3d(), (atomC != null) ? atomC.getPoint3d() : null, nwanted, length, angle);
		} else if (nwithCoords == 2) {
			Point3d bPoint = withCoords.getAtomAt(0).getPoint3d();
			Point3d cPoint = withCoords.getAtomAt(1).getPoint3d();
			newPoints = calculate3DCoordinates2(aPoint, bPoint, cPoint, nwanted, length, angle);
		} else if (nwithCoords == 3) {
			Point3d bPoint = withCoords.getAtomAt(0).getPoint3d();
			Point3d cPoint = withCoords.getAtomAt(1).getPoint3d();
			newPoints = new Point3d[1];
			Point3d dPoint = withCoords.getAtomAt(2).getPoint3d();
			newPoints[0] = calculate3DCoordinates3(aPoint, bPoint, cPoint, dPoint, length);
		}
		//System.out.println("...Ready");
		return newPoints;
	}


	/**
	 *  Calculates substituent points. Calculate substituent points for (0) zero
	 *  ligands of aPoint. 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 where 3x**2 = bond length
	 *
	 *@param  aPoint   to which substituents are added
	 *@param  nwanted  number of points to calculate (1-4)
	 *@param  length   from aPoint
	 *@return          Point3d[] nwanted points (or zero if failed)
	 */
	public Point3d[] calculate3DCoordinates0(Point3d aPoint, int nwanted, double length) {
		Point3d points[] = new Point3d[0];
		if (nwanted == 1) {
			points = new Point3d[1];
			points[0] = new Point3d(aPoint);
			points[0].add(new Vector3d(length, 0.0, 0.0));
		} else if (nwanted == 2) {
			points[0] = new Point3d(aPoint);
			points[0].add(new Vector3d(length, 0.0, 0.0));
			points[1] = new Point3d(aPoint);
			points[1].add(new Vector3d(-length, 0.0, 0.0));
		} else if (nwanted == 3) {
			points[0] = new Point3d(aPoint);
			points[0].add(new Vector3d(length, 0.0, 0.0));
			points[1] = new Point3d(aPoint);
			points[1].add(new Vector3d(-length * 0.5, -length * 0.5 * Math.sqrt(3.0), 0.0f));
			points[2] = new Point3d(aPoint);
			points[2].add(new Vector3d(-length * 0.5, length * 0.5 * Math.sqrt(3.0), 0.0f));
		} else if (nwanted == 4) {
			double dx = length / Math.sqrt(3.0);
			points[0] = new Point3d(aPoint);
			points[0].add(new Vector3d(dx, dx, dx));
			points[1] = new Point3d(aPoint);
			points[1].add(new Vector3d(dx, -dx, -dx));
			points[2] = new Point3d(aPoint);
			points[2].add(new Vector3d(-dx, -dx, dx));
			points[3] = new Point3d(aPoint);
			points[3].add(new Vector3d(-dx, dx, -dx));
		}
		return points;
	}


	/**
	 *  Calculate new point(s) X in a B-A system to form B-A-X. Use C as reference
	 *  for * staggering about the B-A bond (1a) 1 ligand(B) of refAtom (A) which
	 *  itself has a ligand (C) (i) 1 points required; vector along AB vector (ii)
	 *  2 points: 2 vectors in ABC plane, staggered and eclipsed wrt C (iii) 3
	 *  points: 1 staggered wrt C, the others +- gauche wrt C If C is null, a
	 *  random non-colinear C is generated
	 *
	 *@param  aPoint   to which substituents are added
	 *@param  nwanted  number of points to calculate (1-3)
	 *@param  length   A-X length
	 *@param  angle    B-A-X angle
	 *@param  bPoint   Description of the Parameter
	 *@param  cPoint   Description of the Parameter
	 *@return          Point3d[] nwanted points (or zero if failed)
	 */
	public Point3d[] calculate3DCoordinates1(
			Point3d aPoint, Point3d bPoint, Point3d cPoint,
			int nwanted, double length, double angle
			) {
		Point3d points[] = new Point3d[nwanted];
		// BA vector
		Vector3d ba = new Vector3d(aPoint);
		ba.sub(bPoint);
		ba.normalize();
		// if no cPoint, generate a random reference
		if (cPoint == null) {
			Vector3d cVector = getNonColinearVector(ba);
			cPoint = new Point3d(cVector);
		}
		// CB vector
		Vector3d cb = new Vector3d(bPoint);
		cb.sub(cPoint);
		cb.normalize();
		// if A, B, C colinear, replace C by random point
		double cbdotba = cb.dot(ba);
		if (cbdotba > 0.999999) {
			Vector3d cVector = getNonColinearVector(ba);
			cPoint = new Point3d(cVector);
			cb = new Vector3d(bPoint);
			cb.sub(cPoint);
		}
		// cbxba = c x b
		Vector3d cbxba = new Vector3d();
		cbxba.cross(cb, ba);
		cbxba.normalize();
		// create three perp axes ba, cbxba, and ax
		Vector3d ax = new Vector3d();
		ax.cross(cbxba, ba);
		ax.normalize();
		double drot = Math.PI * 2.0 / (double) nwanted;
		for (int i = 0; i < nwanted; i++) {
			double rot = (double) i * drot;
			points[i] = new Point3d(aPoint);
			Vector3d vx = new Vector3d(ba);
			vx.scale(-Math.cos(angle) * length);
			Vector3d vy = new Vector3d(ax);
			vy.scale(Math.cos(rot) * length);
			Vector3d vz = new Vector3d(cbxba);
			vz.scale(Math.sin(rot) * length);
			points[i].add(vx);
			points[i].add(vy);
			points[i].add(vz);
		}
		/*ax = null;
		cbxba = null;
		ba = null;
		cb = null;*/
		return points;
	}


	/**
	 *  Calculate new point(s) X in a B-A-C system. It forms form a B-A(-C)-X
	 *  system. (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
	 *  points X1, X2, X1-A-X2 = angle about 2i vector
	 *
	 *@param  aPoint   to which substituents are added
	 *@param  bPoint   first ligand of A
	 *@param  cPoint   second ligand of A
	 *@param  nwanted  number of points to calculate (1-2)
	 *@param  length   A-X length
	 *@param  angle    B-A-X angle
	 *@return          Point3d[] nwanted points (or zero if failed)
	 */
	public Point3d[] calculate3DCoordinates2(
			Point3d aPoint, Point3d bPoint, Point3d cPoint,
			int nwanted, double length, double angle) {
		System.out.println("3DCoordinates2");
		Point3d newPoints[] = new Point3d[0];
		double ang2 = angle / 2.0;

		Vector3d ba = new Vector3d(aPoint);
		ba.sub(bPoint);
		Vector3d ca = new Vector3d(aPoint);
		ca.sub(cPoint);
		Vector3d baxca = new Vector3d();
		baxca.cross(ba, ca);
		if (baxca.length() < 0.00000001) {
			;
			// linear
		} else if (nwanted == 1) {
			newPoints = new Point3d[1];
			Vector3d ax = new Vector3d(ba);
			ax.add(ca);
			ax.normalize();
			ax.scale(length);
			newPoints[0] = new Point3d(aPoint);
			newPoints[0].add(ax);
		} else if (nwanted >= 2) {
			newPoints = new Point3d[2];
			Vector3d ax = new Vector3d(ba);
			ax.add(ca);
			ax.normalize();
			baxca.normalize();
			baxca.scale(Math.sin(ang2) * length);
			ax.scale(Math.cos(ang2) * length);
			newPoints[0] = new Point3d(aPoint);
			newPoints[0].add(ax);
			newPoints[0].add(baxca);
			newPoints[1] = new Point3d(aPoint);
			newPoints[1].add(ax);
			newPoints[1].sub(baxca);
		}
		baxca = null;
		ba = null;
		ca = null;
		return newPoints;
	}


	/**
	 *  Calculate new point X in a B-A(-D)-C system. It forms a B-A(-D)(-C)-X
	 *  system. (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
	 *
	 *@param  aPoint  to which substituents are added
	 *@param  bPoint  first ligand of A
	 *@param  cPoint  second ligand of A
	 *@param  dPoint  third ligand of A
	 *@param  length  A-X length
	 *@return         Point3d nwanted points (or null if failed (coplanar))
	 */
	public Point3d calculate3DCoordinates3(
			Point3d aPoint, Point3d bPoint, Point3d cPoint, Point3d dPoint,
			double length) {
		System.out.println("3DCoordinates3");
		Vector3d bc = new Vector3d(bPoint);
		bc.sub(cPoint);
		Vector3d dc = new Vector3d(dPoint);
		dc.sub(cPoint);
		Vector3d ca = new Vector3d(cPoint);
		ca.sub(aPoint);

		Vector3d n1 = new Vector3d();
		Vector3d n2 = new Vector3d();

		n1.cross(bc, dc);
		n1.normalize();
		n1.scale(length);

		Vector3d ax = new Vector3d(aPoint);
		ax.add(n1);
		ax.sub(aPoint);

		Vector3d ax2 = new Vector3d(aPoint);
		ax2.add(n2);
		ax2.sub(aPoint);

		Point3d point = new Point3d(aPoint);

		double dotProduct = ca.dot(ax);
		double angle = Math.acos((dotProduct) / (ax.length() * ca.length()));

		if (angle < 1.5) {
			n2.cross(dc, bc);
			n2.normalize();
			n2.scale(length);
			point.add(n2);
		} else {
			point.add(n1);
		}

		bc = null;
		dc = null;
		ca = null;
		n1 = null;
		n2 = null;
		return point;
	}


	/**
	 *  Calculate new point in B-A-C system. It forms B-A(-X)-C system, where A is
	 *  sp2
	 *
	 *@param  aPoint  central point A (Point3d)
	 *@param  bPoint  B (Point3d)
	 *@param  cPoint  C (Point3d)
	 *@param  length  bond length
	 *@param  angle   angle between B(C)-A-X
	 *@return         new Point (Point3d)
	 */
	public Point3d calculate3DCoordinatesSP2_1(Point3d aPoint, Point3d bPoint, Point3d cPoint,
			double length, double angle) {
		System.out.println("3DCoordinatesSP2_1");
		Vector3d ba = new Vector3d(bPoint);
		ba.sub(aPoint);
		Vector3d ca = new Vector3d(cPoint);
		ca.sub(aPoint);

		Vector3d n1 = new Vector3d();
		n1.cross(ba, ca);
		n1.normalize();

		Vector3d n2 = rotate(ba, n1, angle);
		n2.normalize();

		n2.scale(length);
		Point3d point = new Point3d(aPoint);
		point.add(n2);
		n1 = null;
		n2 = null;
		ba = null;
		ca = null;
		return point;
	}


	/**
	 *  Calculate two new points in B-A system. It forms B-A(-X)(-X) system, where
	 *  A is sp2
	 *
	 *@param  aPoint  central point A (Point3d)
	 *@param  bPoint  B (Point3d)
	 *@param  cPoint  C (Point3d)
	 *@param  length  bond length
	 *@param  angle   angle between B(C)-A-X
	 *@return         new Points (Point3d [])
	 */
	public Point3d[] calculate3DCoordinatesSP2_2(Point3d aPoint, Point3d bPoint, Point3d cPoint, double length, double angle) {
		System.out.println("3DCoordinatesSP_2");
		Vector3d ca = new Vector3d();
		Point3d newPoints[] = new Point3d[2];
		Vector3d ba = new Vector3d(bPoint);
		ba.sub(aPoint);

		if (cPoint != null) {
			ca.x = cPoint.x - aPoint.x;
			ca.y = cPoint.y - aPoint.y;
			ca.z = cPoint.z - aPoint.z;
		} else {
			ca.x = -1 * ba.x;
			ca.y = -1 * ba.y;
			ca.z = -1.5 * ba.z;
		}

		Vector3d crossProduct = new Vector3d();
		crossProduct.cross(ba, ca);

		Vector3d n1 = rotate(ba, crossProduct, 2 * angle);
		n1.normalize();
		n1.scale(length);
		newPoints[0] = new Point3d(aPoint);
		newPoints[0].add(n1);

		Vector3d n2 = rotate(n1, ba, Math.PI);
		n2.normalize();
		n2.scale(length);
		newPoints[1] = new Point3d(aPoint);
		newPoints[1].add(n2);
		n1 = null;
		n2 = null;
		ba = null;
		ca = null;
		return newPoints;
	}


	/**
	 *  Gets the nonColinearVector attribute of the AtomLigandPlacer3D class
	 *
	 *@param  ab  Description of the Parameter
	 *@return     The nonColinearVector value
	 */
	private Vector3d getNonColinearVector(Vector3d ab) {
		Vector3d cr = new Vector3d();
		cr.cross(ab, XV);
		if (cr.length() > 0.00001) {
			return XV;
		} else {
			return YV;
		}
	}


	/**
	 *  Rotates a vector around an axis
	 *
	 *@param  vector  vector to be rotated around axis
	 *@param  axis    axis of rotation
	 *@param  angle   angle to vector rotate around
	 *@return         rotated vector
	 *@author         egonw
	 */
	public static Vector3d rotate(Vector3d vector, Vector3d axis, double angle) {
		Matrix3d rotate = new Matrix3d();
		rotate.set(new AxisAngle4d(axis.x, axis.y, axis.z, angle));
		Vector3d result = new Vector3d();
		rotate.transform(vector, result);
		return result;
	}


	/**
	 *  Gets the distance between two atoms out of the parameter set
	 *
	 *@param  id1            id of the paramter set for atom1 (atom1.getID())
	 *@param  id2            id of the paramter set for atom2
	 *@return                The distanceValue value
	 *@exception  Exception  Description of the Exception
	 */
	private double getDistanceValue(String id1, String id2) throws Exception {
		String dkey = "";
		if (pSet.containsKey(("bond" + id1 + ";" + id2))) {
			dkey = "bond" + id1 + ";" + id2;
		} else if (pSet.containsKey(("bond" + id2 + ";" + id1))) {
			dkey = "bond" + id2 + ";" + id1;
		} else {
			System.out.println("DistanceKEYError:pSet has no key:" + id2 + " ; " + id1 + " take default bond length:" + DEFAULT_BOND_LENGTH_H);
			return DEFAULT_BOND_LENGTH_H;
		}
		return ((Double) (((Vector) pSet.get(dkey)).get(0))).doubleValue();
	}


	/**
	 *  Gets the angleKey attribute of the AtomPlacer3D object
	 *
	 *@param  id1            Description of the Parameter
	 *@param  id2            Description of the Parameter
	 *@param  id3            Description of the Parameter
	 *@return                The angleKey value
	 *@exception  Exception  Description of the Exception
	 */
	public double getAngleValue(String id1, String id2, String id3) throws Exception {
		String akey = "";
		if (pSet.containsKey(("angle" + id1 + ";" + id2 + ";" + id3))) {
			akey = "angle" + id1 + ";" + id2 + ";" + id3;
		} else if (pSet.containsKey(("angle" + id3 + ";" + id2 + ";" + id1))) {
			akey = "angle" + id3 + ";" + id2 + ";" + id1;
		} else if (pSet.containsKey(("angle" + id2 + ";" + id1 + ";" + id3))) {
			akey = "angle" + id2 + ";" + id1 + ";" + id3;
		} else if (pSet.containsKey(("angle" + id1 + ";" + id3 + ";" + id2))) {
			akey = "angle" + id1 + ";" + id3 + ";" + id2;
		} else if (pSet.containsKey(("angle" + id3 + ";" + id1 + ";" + id2))) {
			akey = "angle" + id3 + ";" + id1 + ";" + id2;
		} else if (pSet.containsKey(("angle" + id2 + ";" + id3 + ";" + id1))) {
			akey = "angle" + id2 + ";" + id3 + ";" + id1;
		} else {
			System.out.println("AngleKEYError:Unknown angle " + id1 + " " + id2 + " " + id3 + " take default angle:" + TETRAHEDRAL_ANGLE);
			return TETRAHEDRAL_ANGLE;
		}
		return ((Double) (((Vector) pSet.get(akey)).get(0))).doubleValue();
	}


	/**
	 *  set Atoms in respect to stereoinformation
	 *	take placed neighbours to stereocenter 
	 *		create a x b
	 *	     if right handed system (spatproduct >0)
	 *			if unplaced !=up (relativ to stereocenter)
	 *				n=b x a
	 *	     Determine angle between n and possible ligand place points
	 *	     if angle smaller than 90° take this branch point
	 *
	 *@param  atomA         placed Atom - stereocenter
	 *@param  ax            bond between stereocenter and unplaced atom
	 *@param  atomB         neighbour of atomA (in plane created by atomA, atomB and atomC)
	 *@param  atomC         neighbour of atomA
	 *@param  branchPoints  the two possible placement points for unplaced atom (up and down) 
	 *@return               int value of branch point position
	 */
	public int makeStereocenter(Point3d atomA, Bond ax, Point3d atomB, Point3d atomC, Point3d[] branchPoints) {
		
		Vector3d b = new Vector3d((atomB.x - atomA.x), (atomB.y - atomA.y), (atomB.z - atomA.z));
		Vector3d c = new Vector3d((atomC.x - atomA.x), (atomC.y - atomA.y), (atomC.z - atomA.z));
		Vector3d n1 = new Vector3d();
		Vector3d n2 = null;
		n1.cross(b, c);
		n1.normalize();

		if (getSpatproduct(b, c, n1) >= 0) {
			if (ax.getStereo() != CDKConstants.STEREO_BOND_UP_INV) {
				n1.cross(c, b);
				n1.normalize();
			}
		}
		double dotProduct = 0;
		for (int i = 0; i < branchPoints.length; i++) {
			n2 = new Vector3d(branchPoints[0].x, branchPoints[0].y, branchPoints[0].z);
			dotProduct = n1.dot(n2);
			if (Math.acos(dotProduct / (n1.length() * n2.length())) < 1.6) {
				return i;
			}
		}
		return -1;
	}


	/**
	 *  Gets the spatproduct of three vectors
	 *
	 *@param  a  vector a
	 *@param  b  vector b
	 *@param  c  vector c
	 *@return    double value of the spatproduct
	 */
	public double getSpatproduct(Vector3d a, Vector3d b, Vector3d c) {
		return (c.x * (b.y * a.z - b.z * a.y) + c.y * (b.z * a.x - b.x * a.z) + c.z * (b.x * a.y - b.y * a.x));
	}


	/**
	 *  Calculates the torsionAngle of a-b-c-d 
	 *
	 *@param  a  Point3d
	 *@param  b  Point3d
	 *@param  c  Point3d
	 *@param  d  Point3d
	 *@return    The torsionAngle value
	 */
	public double getTorsionAngle(Point3d a, Point3d b, Point3d c, Point3d d) {
		Vector3d ab = new Vector3d(a.x - b.x, a.y - b.y, a.z - b.z);
		Vector3d cb = new Vector3d(c.x - b.x, c.y - b.y, c.z - b.z);
		Vector3d dc = new Vector3d(d.x - c.x, d.y - c.y, d.z - c.z);
		Vector3d bc = new Vector3d(b.x - c.x, b.y - c.y, b.z - c.z);
		Vector3d n1 = new Vector3d();
		Vector3d n2 = new Vector3d();

		n1.cross(ab, cb);
		if (getSpatproduct(ab, cb, n1) > 0) {
			n1.cross(cb, ab);
		}
		n1.normalize();
		n2.cross(dc, bc);
		if (getSpatproduct(dc, bc, n2) < 0) {
			n2.cross(bc, dc);
		}
		n2.normalize();
		return n1.dot(n2);
	}


	/**
	 *  Gets all placed neighbouring atoms of a atom
	 *
	 *@param  atom  central atom (Atom)
	 *@param  ac    the molecul
	 *@return       all connected and placed atoms to the central atom
	 *      ((AtomContainer)
	 */
	public AtomContainer getPlacedAtomsInAtomContainer(Atom atom, AtomContainer ac) {

		Bond[] bonds = ac.getConnectedBonds(atom);
		AtomContainer connectedAtoms = new AtomContainer();
		Atom connectedAtom = null;
		for (int i = 0; i < bonds.length; i++) {
			connectedAtom = bonds[i].getConnectedAtom(atom);
			if (connectedAtom.getFlag(CDKConstants.ISPLACED)) {
				connectedAtoms.addAtom(connectedAtom);
			}
		}
		return connectedAtoms;
	}


	/**
	 *  Gets the unsetAtomsInAtomContainer attribute of the
	 *  AtomTetrahedralLigandPlacer3D object
	 *
	 *@param  atom  Description of the Parameter
	 *@param  ac    Description of the Parameter
	 *@return       The unsetAtomsInAtomContainer value
	 */
	public AtomContainer getUnsetAtomsInAtomContainer(Atom atom, AtomContainer ac) {
		Atom[] atoms = ac.getConnectedAtoms(atom);
		AtomContainer connectedAtoms = new AtomContainer();
		for (int i = 0; i < atoms.length; i++) {
			if (!atoms[i].getFlag(CDKConstants.ISPLACED)){//&& atoms[i].getPoint3d() == null) {
				connectedAtoms.addAtom(atoms[i]);
			}
		}
		return connectedAtoms;
	}


	/**
	 *  Description of the Method
	 *
	 *@param  atom  Description of the Parameter
	 *@param  ac    Description of the Parameter
	 *@return       Description of the Return Value
	 */
	public boolean hasUnsetNeighbour(Atom atom, AtomContainer ac) {
		Atom[] atoms = ac.getConnectedAtoms(atom);
		for (int i = 0; i < atoms.length; i++) {
			if (!atoms[i].getFlag(CDKConstants.ISPLACED)) {//&& atoms[i].getPoint3d() == null) {
				return true;
			}
		}
		return false;
	}


	/**
	 *  Returns a placed neighbouring atom of a central atom atomA, which is not
	 *  atomB
	 *
	 *@param  atomA  central atom (Atom)
	 *@param  atomB  atom connected to atomA (Atom)
	 *@param  ac     molecule
	 *@return        returns a connected atom (Atom)
	 */
	public Atom getPlacedHeavyAtomInAtomContainer(Atom atomA, Atom atomB, AtomContainer ac) {
		Atom[] atoms = ac.getConnectedAtoms(atomA);
		Atom atom=null;
		for (int i = 0; i < atoms.length; i++) {
			if (atoms[i].getFlag(CDKConstants.ISPLACED) && !atoms[i].getSymbol().equals("H")
					 && atoms[i] != atomB) {
				return atoms[i];
			}
		}
		return atom;
	}
}


--- NEW FILE: ForceFieldConfigurator.java ---
/*
 *  $RCSfile: ForceFieldConfigurator.java,v $
 *  $Author: chhoppe
 *  $Date: 2004/11/30 10:24:48 $
 *  Copyright (C) 2003-2004  The Chemistry Development Kit (CDK) project
 *
 *  Contact: cdk...@li...
 *
 *  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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
package org.openscience.cdk.modeling.builder3d;

import java.io.*;
import java.util.Vector;
import java.util.Hashtable;
import java.util.regex.*;

import org.openscience.cdk.AtomType;
import org.openscience.cdk.exception.NoSuchAtomTypeException;
import org.openscience.cdk.Atom;
import org.openscience.cdk.PseudoAtom;
import org.openscience.cdk.CDKConstants;

/**
 *  Reads in a force field configuration file, set the atom types into a vector, and the data into a hashtable
 *  Therfore it uses the class MM2BasedParamterSetReader
 *  private Hashtable parameterSet;
 *  key=nameofdatafield+atomid1+;atomid2;atomxid
 *  
 *  MM2 and MMFF94 force field are implemented 
 *  With force field data it configures the cdk atom (assign atomtype, van der Waals radius, charge...)
 *
 * @author     chhoppe
 * @cdk.created    2004-09-07
 */
public class ForceFieldConfigurator {

	private File f;
	private String ffName = "mm2";
	private Vector atomTypes;
	private Hashtable parameterSet=null;
	private MM2BasedParameterSetReader mm2 = null;
	private MMFF94BasedParameterSetReader mmff94= null;
	private InputStream ins = null;
	private String[] fftypes = {"mm2","mmff94"};
	
	/**
	 *Constructor for the ForceFieldConfigurator object
	 */
	ForceFieldConfigurator() { }

	/**
	 *  Sets the inputStream attribute of the ForceFieldConfigurator object
	 *
	 * @param  ins  The new inputStream value
	 */
	public void setInputStream(InputStream ins) {
		this.ins = ins;
	}

	/**
	 *  Sets the forceFieldType attribute of the ForceFieldConfigurator object
	 *
	 * @param  ffname  The new forceFieldType name
	 */
	public void checkForceFieldType(String ffname) {
		ffName = ffname.toLowerCase();
		boolean check=false;
		for (int i = 0; i <= fftypes.length; i++) {
			if (fftypes[i].equals(ffName)) {
				check=true;
				break;
			} 
		}
		if (!check) {
			System.out.println("FFError:checkForceFieldType> Unknown forcefield:" + ffName);
		}
	}

	
	private void readFile(File f){
		if (f.exists()) {
			try{
				ins = new FileInputStream(f);
			}catch (Exception ex1){
			}
		} else {
			System.out.println("IOError: ConfigFile doesent exist");
		}
	}
	/**
	 *Constructor for the ForceFieldConfigurator object
	 *
	 * @param  ffname  name of the force field data file
	 */
	public void setForceFieldConfigurator(String ffname) throws Exception {
		if (ffname.toLowerCase()==ffName && parameterSet!=null){
		}else{
			this.checkForceFieldType(ffname);
			if (ffName.equals("mm2")) {
				//System.out.println("ForceFieldConfigurator: open Force Field mm2");
				//f = new File(mm2File);
				//readFile(f);
				ins = this.getClass().getClassLoader().getResourceAsStream("data/mm2.prm");
				mm2 = new MM2BasedParameterSetReader();
				mm2.setInputStream(ins);
				this.setMM2Parameters();
			}else if (ffName.equals("mmff94")) {
				//System.out.println("ForceFieldConfigurator: open Force Field mmff94");
				//f = new File(mmff94File);
				//readFile(f);
				ins = this.getClass().getClassLoader().getResourceAsStream("data/mmff94.prm");
				mmff94= new MMFF94BasedParameterSetReader();
				mmff94.setInputStream(ins);
				this.setMMFF94Parameters();
			}
		}
	}


	/**
	 *  Sets the atomTypes attribute of the ForceFieldConfigurator object
	 *
	 * @param  atomtypes  The new atomTypes 
	 */
	public void setAtomTypes(Vector atomtypes) {
		atomTypes = atomtypes;
	}

	/**
	 *  Sets the parameters attribute of the ForceFieldConfigurator object
	 *
	 * @param  parameterset  The new parameter values
	 */
	public void setParameters(Hashtable parameterset) {
		parameterSet = parameterset;
	}

	/**
	 *  Sets the parameters attribute of the ForceFieldConfigurator object, default is mm2 force field
	 */
	public void setMM2Parameters() throws Exception{
		mm2.readParameterSets();
		parameterSet = mm2.getParamterSet();
		atomTypes = mm2.getAtomTypes();
	}
	
	public void setMMFF94Parameters() throws Exception{
		mmff94.readParameterSets();
		parameterSet = mmff94.getParamterSet();
		atomTypes = mmff94.getAtomTypes();
	}
	
	/**
	 *  Gets the atomTypes attribute of the ForceFieldConfigurator object
	 *
	 * @return    The atomTypes vector
	 */
	public Vector getAtomTypes() {
		return atomTypes;
	}

	/**
	 *  Gets the parameterSet attribute of the ForceFieldConfigurator object
	 *
	 * @return    The parameterSet hashtable
	 */
	public Hashtable getParameterSet() {
		return this.parameterSet;
	}

	/**
	 *  Find the atomType for a id
	 *
	 * @param  ID                           Atomtype id of the forcefield
	 * @return                              The atomType 
	 * @exception  NoSuchAtomTypeException  atomType is not known
	 */
	private AtomType getAtomType(String ID) throws NoSuchAtomTypeException {
		AtomType at = null;
		for (int i = 0; i < atomTypes.size(); i++) {
			at = (AtomType) atomTypes.get(i);
			if (at.getID().equals(ID)) {
				return at;
			}
		}
		throw new NoSuchAtomTypeException("AtomType " + ID + " could not be found");
	}

	/**
	 *  Assigns an atom type to an atom
	 *
	 * @param  atom  The atom to be aasigned
	 * @param  ID    the atom type id
	 * @return       the assigned atom
	 */
	private Atom setAtom(Atom atom, String ID) throws Exception {
		AtomType at = null;
		String key = "";
		Vector data = null;
		Double value = null;
		
		at = getAtomType(ID);
		if (atom.getSymbol()==null){
			atom.setSymbol(at.getSymbol());
		}
		atom.setAtomTypeName(at.getAtomTypeName());
		atom.setFormalNeighbourCount(at.getFormalNeighbourCount());
		key = "vdw" + ID;
		data = (Vector) parameterSet.get(key);
		value = (Double) data.firstElement();
		atom.setVanderwaalsRadius(value.doubleValue());
		key = "charge" + ID;
		if (parameterSet.containsKey(key)) {
			data = (Vector) parameterSet.get(key);
			value = (Double) data.firstElement();
			atom.setCharge(value.doubleValue());
		}
		Object color = at.getProperty("org.openscience.cdk.renderer.color");
		if (color != null) {
			atom.setProperty("org.openscience.cdk.renderer.color", color);
		}
		if (at.getAtomicNumber() != 0) {
			atom.setAtomicNumber(at.getAtomicNumber());
		} 
		if (at.getExactMass() > 0.0) {
			atom.setExactMass(at.getExactMass());
		} 
		return atom;
	}
	
	public Atom configureAtom(Atom atom, String hoseCode, boolean _boolean) throws Exception {
		if (ffName.equals("mm2")){
			return configureMM2BasedAtom(atom, hoseCode,_boolean);
		}else i...
 
[truncated message content]