geometrytoolsinternalcoordinates.java

化学图形处理软件

			for (int j=0;j<secondAtomContainer.getAtomCount();j++){
				tmp=Math.floor(distanceMatrix[i][j]*10);
				//logger.debug(tmp/10+"\t");
			}			
		}*/
		
		double minimumDistance;
        for(int i=0;i<firstAtomContainer.getAtomCount();i++){
            minimumDistance=searchRadius;
            for (int j=0;j<secondAtomContainer.getAtomCount();j++){
                if(distanceMatrix[i][j]< searchRadius && distanceMatrix[i][j]< minimumDistance){
                    //logger.debug("Distance OK "+i+" "+j+":"+distanceMatrix[i][j]+" AtomCheck:"+checkAtomMapping(firstAtomContainer,secondAtomContainer, i, j));
                    //check atom properties
                    if (checkAtomMapping(firstAtomContainer,secondAtomContainer, i, j)){
                        minimumDistance=distanceMatrix[i][j];
                        mappedAtoms.put(new Integer(firstAtomContainer.getAtomNumber(firstAtomContainer.getAtom(i))),new Integer(secondAtomContainer.getAtomNumber(secondAtomContainer.getAtom(j))));
                        //firstAtomContainer.getAtomAt(i).setProperty("MappedAtom",new Integer(secondAtomContainer.getAtomNumber(secondAtomContainer.getAtomAt(j))));
                        //logger.debug("#:"+countMappedAtoms+" Atom:"+i+" is mapped to Atom"+j);
                        //logger.debug(firstAtomContainer.getConnectedAtoms(firstAtomContainer.getAtomAt(i)).length);
                    }
                }
            }
        }
		return mappedAtoms;
	}

    /**
     * Returns a Map with the AtomNumbers, the first number corresponds to the first (or the largest
     * AtomContainer) atomContainer.
     * <p/>
     * Only for similar and aligned molecules with coordinates!
     *
     * @param firstAtomContainer  the (largest) first aligned AtomContainer which is the reference
     * @param secondAtomContainer the second aligned AtomContainer
     * @return a Map of the mapped atoms
     * @throws CDKException if there is an error in the UniversalIsomorphismTester
     */
    public static Map mapAtomsOfAlignedStructures(IAtomContainer firstAtomContainer, IAtomContainer secondAtomContainer, Map mappedAtoms) throws CDKException {
        //logger.debug("**** GT MAP ATOMS ****");
        //Map atoms onto each other
        if (firstAtomContainer.getAtomCount() < 1 & secondAtomContainer.getAtomCount() < 1) {
            return mappedAtoms;
        }
        RMap map;
        org.openscience.cdk.interfaces.IAtom atom1;
        org.openscience.cdk.interfaces.IAtom atom2;
        List list;
        try {
            list = UniversalIsomorphismTester.getSubgraphAtomsMap(firstAtomContainer, secondAtomContainer);
            //logger.debug("ListSize:"+list.size());
            for (int i = 0; i < list.size(); i++) {
                map = (RMap) list.get(i);
                atom1 = firstAtomContainer.getAtom(map.getId1());
                atom2 = secondAtomContainer.getAtom(map.getId2());
                if (checkAtomMapping(firstAtomContainer, secondAtomContainer, firstAtomContainer.getAtomNumber(atom1), secondAtomContainer.getAtomNumber(atom2)))
                {
                    mappedAtoms.put(new Integer(firstAtomContainer.getAtomNumber(atom1)), new Integer(secondAtomContainer.getAtomNumber(atom2)));
                    //logger.debug("#:"+countMappedAtoms+" Atom:"+firstAtomContainer.getAtomNumber(atom1)+" is mapped to Atom:"+secondAtomContainer.getAtomNumber(atom2));
                } else {
                    logger.error("Error: Atoms are not similar !!");
                }
            }
        } catch (CDKException e) {
            throw new CDKException("Error in UniversalIsomorphismTester due to:" + e.toString());
        }
        return mappedAtoms;
    }


    private static void getLargestAtomContainer(IAtomContainer firstAC, IAtomContainer secondAC) {
        if (firstAC.getAtomCount() < secondAC.getAtomCount()){
            IAtomContainer tmp;
            try {
                tmp = (IAtomContainer) firstAC.clone();
                firstAC=(IAtomContainer)secondAC.clone();
                secondAC=(IAtomContainer)tmp.clone();
            } catch (CloneNotSupportedException e) {
                // TODO Auto-generated catch block
                e.printStackTrace();
            }
        }
    }
	
	private static boolean checkAtomMapping(IAtomContainer firstAC, IAtomContainer secondAC, int posFirstAtom, int posSecondAtom){
		IAtom firstAtom=firstAC.getAtom(posFirstAtom);
		IAtom secondAtom=secondAC.getAtom(posSecondAtom);
		if (firstAtom.getSymbol().equals(secondAtom.getSymbol()) && firstAC.getConnectedAtomsList(firstAtom).size() == secondAC.getConnectedAtomsList(secondAtom).size() &&
				firstAtom.getBondOrderSum() == secondAtom.getBondOrderSum() &&
				firstAtom.getMaxBondOrder() == secondAtom.getMaxBondOrder() 
		        ){
			return true;
		}else {
			return false;
		}
	}
	
	private static IAtomContainer setVisitedFlagsToFalse(IAtomContainer atomContainer) {
        for (int i=0;i<atomContainer.getAtomCount();i++){
            atomContainer.getAtom(i).setFlag(CDKConstants.VISITED, false);
        }
        return atomContainer;
    }
	
	/**
	 *  Return the RMSD of bonds length between the 2 aligned molecules.
	 *
	 *@param  firstAtomContainer                the (largest) first aligned AtomContainer which is the reference
	 *@param  secondAtomContainer               the second aligned AtomContainer
	 *@param  mappedAtoms             			Map: a Map of the mapped atoms
	 *@param  Coords3d            			    boolean: true if moecules has 3D coords, false if molecules has 2D coords
	 *@return                   				double: all the RMSD of bonds length
	 *
	 **/
	public static double getBondLengthRMSD(IAtomContainer firstAtomContainer, IAtomContainer secondAtomContainer,Map mappedAtoms, boolean Coords3d) {
        //logger.debug("**** GT getBondLengthRMSD ****");
        Iterator firstAtoms=mappedAtoms.keySet().iterator();
        IAtom centerAtomFirstMolecule;
        IAtom centerAtomSecondMolecule;
        java.util.List connectedAtoms;
        double sum=0;
        double n=0;
        double distance1=0;
        double distance2=0;
        setVisitedFlagsToFalse(firstAtomContainer);
        setVisitedFlagsToFalse(secondAtomContainer);
        while(firstAtoms.hasNext()){
            centerAtomFirstMolecule=firstAtomContainer.getAtom(((Integer)firstAtoms.next()).intValue());
            centerAtomFirstMolecule.setFlag(CDKConstants.VISITED, true);
            centerAtomSecondMolecule=secondAtomContainer.getAtom(((Integer)mappedAtoms.get(new Integer(firstAtomContainer.getAtomNumber(centerAtomFirstMolecule)))).intValue());
            connectedAtoms=firstAtomContainer.getConnectedAtomsList(centerAtomFirstMolecule);
            for (int i=0;i<connectedAtoms.size();i++){
            	IAtom conAtom = (IAtom)connectedAtoms.get(i);
                //this step is built to know if the program has already calculate a bond length (so as not to have duplicate values)
                if(!conAtom.getFlag(CDKConstants.VISITED)){
                    if (Coords3d){
                        distance1=((Point3d)centerAtomFirstMolecule.getPoint3d()).distance(conAtom.getPoint3d());
                        distance2=((Point3d)centerAtomSecondMolecule.getPoint3d()).distance(secondAtomContainer.getAtom(((Integer)mappedAtoms.get(new Integer(firstAtomContainer.getAtomNumber(conAtom)))).intValue()).getPoint3d());
                        sum=sum+Math.pow((distance1-distance2),2);
                        n++;
                    }else{
                        distance1=((Point2d)centerAtomFirstMolecule.getPoint2d()).distance(conAtom.getPoint2d());
                        distance2=((Point2d)centerAtomSecondMolecule.getPoint2d()).distance(secondAtomContainer.getAtom(((Integer)mappedAtoms.get(new Integer(firstAtomContainer.getAtomNumber(conAtom)))).intValue()).getPoint2d());
                        sum=sum+Math.pow((distance1-distance2),2);
                        n++;
                    }
                }
            }
        }
        setVisitedFlagsToFalse(firstAtomContainer);
        setVisitedFlagsToFalse(secondAtomContainer);
        return Math.sqrt(sum/n);
    }
	/**
	 *  Return the variation of each angle value between the 2 aligned molecules.
	 *
	 *@param  firstAtomContainer                the (largest) first aligned AtomContainer which is the reference
	 *@param  secondAtomContainer               the second aligned AtomContainer
	 *@param  mappedAtoms             			Map: a Map of the mapped atoms
	 *@return                   				double: the value of the RMSD
	 *
	 **/
	public static double getAngleRMSD(IAtomContainer firstAtomContainer, IAtomContainer secondAtomContainer, Map mappedAtoms) {
        //logger.debug("**** GT getAngleRMSD ****");
        Iterator firstAtoms=mappedAtoms.keySet().iterator();
        //logger.debug("mappedAtoms:"+mappedAtoms.toString());
        IAtom firstAtomfirstAC;
        IAtom centerAtomfirstAC;
        IAtom firstAtomsecondAC;
        IAtom secondAtomsecondAC;
        IAtom centerAtomsecondAC;
        double angleFirstMolecule;
        double angleSecondMolecule;
        double sum=0;
        double n=0;
        while(firstAtoms.hasNext()){
            int firstAtomNumber=((Integer)firstAtoms.next()).intValue();
            centerAtomfirstAC=firstAtomContainer.getAtom(firstAtomNumber);
            java.util.List connectedAtoms=firstAtomContainer.getConnectedAtomsList(centerAtomfirstAC);
            if (connectedAtoms.size() >1){
                //logger.debug("If "+centerAtomfirstAC.getSymbol()+" is the center atom :");
                for (int i=0; i < connectedAtoms.size()-1;i++){
                    firstAtomfirstAC=(IAtom)connectedAtoms.get(i);
                    for (int j=i+1; j < connectedAtoms.size();j++){
                        angleFirstMolecule=getAngle(centerAtomfirstAC,firstAtomfirstAC,(IAtom)connectedAtoms.get(j));
                        centerAtomsecondAC=secondAtomContainer.getAtom(((Integer)mappedAtoms.get(new Integer(firstAtomContainer.getAtomNumber(centerAtomfirstAC)))).intValue());
                        firstAtomsecondAC=secondAtomContainer.getAtom(((Integer)mappedAtoms.get(new Integer(firstAtomContainer.getAtomNumber(firstAtomfirstAC)))).intValue());
                        secondAtomsecondAC=secondAtomContainer.getAtom(((Integer)mappedAtoms.get(new Integer(firstAtomContainer.getAtomNumber((IAtom)connectedAtoms.get(j))))).intValue());
                        angleSecondMolecule=getAngle(centerAtomsecondAC,firstAtomsecondAC,secondAtomsecondAC);
                        sum=sum+Math.pow(angleFirstMolecule-angleSecondMolecule,2);
                        n++;
                        //logger.debug("Error for the "+firstAtomfirstAC.getSymbol().toLowerCase()+"-"+centerAtomfirstAC.getSymbol()+"-"+connectedAtoms[j].getSymbol().toLowerCase()+" Angle :"+deltaAngle+" degrees");
                    }
                }
            }//if
        }
        return Math.sqrt(sum/n);
    }
	
	private static double getAngle(IAtom atom1, IAtom atom2, IAtom atom3){
		
		Vector3d centerAtom = new Vector3d();
		centerAtom.x=atom1.getPoint3d().x;
		centerAtom.y=atom1.getPoint3d().y;
		centerAtom.z=atom1.getPoint3d().z;
		Vector3d firstAtom = new Vector3d();
		Vector3d secondAtom = new Vector3d();
			
		firstAtom.x=atom2.getPoint3d().x;
		firstAtom.y=atom2.getPoint3d().y;
		firstAtom.z=atom2.getPoint3d().z;
				
		secondAtom.x=atom3.getPoint3d().x;
		secondAtom.y=atom3.getPoint3d().y;
		secondAtom.z=atom3.getPoint3d().z;
				
		firstAtom.sub(centerAtom);
		secondAtom.sub(centerAtom);
				
		return firstAtom.angle(secondAtom);
	}
	
	/**
	 *  Return the RMSD between the 2 aligned molecules.
	 *
	 *@param  firstAtomContainer                the (largest) first aligned AtomContainer which is the reference
	 *@param  secondAtomContainer               the second aligned AtomContainer
	 *@param  mappedAtoms             			Map: a Map of the mapped atoms
	 *@param  Coords3d            			    boolean: true if moecules has 3D coords, false if molecules has 2D coords
	 *@return                   				double: the value of the RMSD 
	 *@exception  CDKException
	 *
	 **/
	public static double getAllAtomRMSD(IAtomContainer firstAtomContainer, IAtomContainer secondAtomContainer, Map mappedAtoms, boolean Coords3d)throws CDKException {
		//logger.debug("**** GT getAllAtomRMSD ****");
		double sum=0;
		double RMSD;
		Iterator firstAtoms=mappedAtoms.keySet().iterator();
		int firstAtomNumber;
		int secondAtomNumber;
		int n=0;
		while(firstAtoms.hasNext()){
			firstAtomNumber=((Integer)firstAtoms.next()).intValue();
			try{
				secondAtomNumber=((Integer)mappedAtoms.get(new Integer(firstAtomNumber))).intValue();
				IAtom firstAtom=firstAtomContainer.getAtom(firstAtomNumber);
				if (Coords3d){
					sum=sum+Math.pow(firstAtom.getPoint3d().distance(secondAtomContainer.getAtom(secondAtomNumber).getPoint3d()),2);
					n++;
				}else{
					sum=sum+Math.pow(firstAtom.getPoint2d().distance(secondAtomContainer.getAtom(secondAtomNumber).getPoint2d()),2);
					n++;
				}
			}catch (Exception ex){
                throw new CDKException(ex.toString());
            }
		}
		RMSD=Math.sqrt(sum/n);
		return RMSD;
	}
	
	/**
	 *  Return the RMSD of the heavy atoms between the 2 aligned molecules.
	 *
	 *@param  firstAtomContainer                the (largest) first aligned AtomContainer which is the reference
	 *@param  secondAtomContainer               the second aligned AtomContainer
	 *@param  mappedAtoms             			Map: a Map of the mapped atoms
	 *@param hetAtomOnly                        boolean: true if only hetero atoms should be considered
	 *@param  Coords3d            			    boolean: true if moecules has 3D coords, false if molecules has 2D coords
	 *@return                   				double: the value of the RMSD 
	 *@exception  CDK
	 *
	 **/
	public static double getHeavyAtomRMSD(IAtomContainer firstAtomContainer, IAtomContainer secondAtomContainer, Map mappedAtoms, boolean hetAtomOnly ,boolean Coords3d)throws CDKException {
		//logger.debug("**** GT getAllAtomRMSD ****");
		double sum=0;
		double RMSD=0;
		Iterator firstAtoms=mappedAtoms.keySet().iterator();
		int firstAtomNumber=0;
		int secondAtomNumber=0;
		int n=0;
		while(firstAtoms.hasNext()){
			firstAtomNumber=((Integer)firstAtoms.next()).intValue();
			try{
				secondAtomNumber=((Integer)mappedAtoms.get(new Integer(firstAtomNumber))).intValue();
				IAtom firstAtom=firstAtomContainer.getAtom(firstAtomNumber);
				if (hetAtomOnly){
					if (!firstAtom.getSymbol().equals("H") && !firstAtom.getSymbol().equals("C")){
						if (Coords3d){
							sum=sum+Math.pow(((Point3d)firstAtom.getPoint3d()).distance(secondAtomContainer.getAtom(secondAtomNumber).getPoint3d()),2);
							n++;
						}else{
							sum=sum+Math.pow(((Point2d)firstAtom.getPoint2d()).distance(secondAtomContainer.getAtom(secondAtomNumber).getPoint2d()),2);
							n++;
						}
					}
				}else{
					if (!firstAtom.getSymbol().equals("H")){
						if (Coords3d){
							sum=sum+Math.pow(((Point3d)firstAtom.getPoint3d()).distance(secondAtomContainer.getAtom(secondAtomNumber).getPoint3d()),2);
							n++;
						}else{
							sum=sum+Math.pow(((Point2d)firstAtom.getPoint2d()).distance(secondAtomContainer.getAtom(secondAtomNumber).getPoint2d()),2);
							n++;
						}
					}
				}
			}catch (Exception ex){
			}
		}
		RMSD=Math.sqrt(sum/n);
		return RMSD;
	}



	/**
	 *  An average of all 3D bond length values is produced, using point3ds in atoms.
	 *  Atom's with no coordinates are disregarded.
	 *
	 *@param  ac  The AtomContainer for which the average bond length is to be
	 *      calculated
	 *@return     the average bond length
	 */
	public static double getBondLengthAverage3D(IAtomContainer ac) {
		double bondLengthSum = 0;
        Iterator bonds = ac.bonds();
        int bondCounter = 0;
        while (bonds.hasNext()) {
            IBond bond = (IBond) bonds.next();
			org.openscience.cdk.interfaces.IAtom atom1 = bond.getAtom(0);
			org.openscience.cdk.interfaces.IAtom atom2 = bond.getAtom(1);
			if (atom1.getPoint3d() != null &&
					atom2.getPoint3d() != null) {
				bondCounter++;
				bondLengthSum += atom1.getPoint3d().distance(atom2.getPoint3d());
			}
		}
		return bondLengthSum / bondCounter;
	}
}