bcutdescriptor.java

化学图形处理软件

    static private class BurdenMatrix {

        static double[][] evalMatrix(IAtomContainer atomContainer, double[] vsd) {
            IAtomContainer local = AtomContainerManipulator.removeHydrogens(atomContainer);

            int natom = local.getAtomCount();
            double[][] matrix = new double[natom][natom];

            /* set the off diagonal entries */
            for (int i = 0; i < natom - 1; i++) {
                for (int j = i + 1; j < natom; j++) {
                    for (int k = 0; k < local.getBondCount(); k++) {
                        org.openscience.cdk.interfaces.IBond bond = local.getBond(k);
                        if (bond.contains(local.getAtom(i)) && bond.contains(local.getAtom(j))) {
                            if (bond.getOrder() == CDKConstants.BONDORDER_SINGLE) matrix[i][j] = 0.1;
                            else if (bond.getOrder() == CDKConstants.BONDORDER_DOUBLE) matrix[i][j] = 0.2;
                            else if (bond.getOrder() == CDKConstants.BONDORDER_TRIPLE) matrix[i][j] = 0.3;
                            else if (bond.getOrder() == CDKConstants.BONDORDER_AROMATIC) matrix[i][j] = 0.15;

                            if (local.getConnectedBondsCount(i) == 1 || local.getConnectedBondsCount(j) == 1) {
                                matrix[i][j] += 0.01;
                            }
                            matrix[j][i] = matrix[i][j];
                        } else {
                            matrix[i][j] = 0.001;
                            matrix[j][i] = 0.001;
                        }
                    }
                }
            }

            /* set the diagonal entries */
            for (int i = 0; i < natom; i++) {
                if (vsd != null) matrix[i][i] = vsd[i];
                else matrix[i][i] = 0.0;
            }
            return (matrix);
        }
    }

    /**
     * Calculates the three classes of BCUT descriptors.
     *
     * @param container Parameter is the atom container.
     * @return An ArrayList containing the descriptors. The default is to return
     *         all calculated eigenvalues of the Burden matrices in the order described
     *         above. If a parameter list was supplied, then only the specified number
     *         of highest and lowest eigenvalues (for each class of BCUT) will be returned.
     * @throws CDKException if the wrong number of eigenvalues are requested (negative or more than the number
     *                      of heavy atoms)
     */
    public DescriptorValue calculate(IAtomContainer container) throws CDKException {
        int counter;
        Molecule molecule;
        try {
            molecule = (Molecule) container.clone();
        } catch (CloneNotSupportedException e) {
            logger.debug("Error during clone");
            throw new CDKException("Error occured during clone "+e);
        }

        // add H's in case they're not present
        HydrogenAdder hydrogenAdder = new HydrogenAdder();
        try {
            hydrogenAdder.addExplicitHydrogensToSatisfyValency(molecule);
        } catch (Exception e) {
            throw new CDKException("Could not add hydrogens: " + e.getMessage(), e);
        }

        // do aromaticity detecttion for calculating polarizability later on
        if (this.checkAromaticity) {
            HueckelAromaticityDetector.detectAromaticity(molecule);
        }

        // find number of heavy atoms
        int nheavy = 0;
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            if (!molecule.getAtom(i).getSymbol().equals("H")) nheavy++;
        }

        if (this.nhigh > nheavy || this.nlow > nheavy) {
            throw new CDKException("Number of negative or positive eigenvalues cannot be more than number of heavy atoms");
        }

        double[] diagvalue = new double[nheavy];

        // get atomic mass weighted BCUT
        counter = 0;
        try {
            for (int i = 0; i < molecule.getAtomCount(); i++) {
                if (molecule.getAtom(i).getSymbol().equals("H")) continue;
                diagvalue[counter] = IsotopeFactory.getInstance(molecule.getBuilder()).
                        getMajorIsotope(molecule.getAtom(i).getSymbol()).getExactMass();
                counter++;
            }
        } catch (Exception e) {
            throw new CDKException("Could not calculate weight: " + e.getMessage(), e);
        }
        double[][] burdenMatrix = BurdenMatrix.evalMatrix(molecule, diagvalue);
        Matrix matrix = new Matrix(burdenMatrix);
        EigenvalueDecomposition eigenDecomposition = new EigenvalueDecomposition(matrix);
        double[] eval1 = eigenDecomposition.getRealEigenvalues();

        // get charge weighted BCUT
        GasteigerMarsiliPartialCharges peoe;
        try {
            peoe = new GasteigerMarsiliPartialCharges();
            peoe.assignGasteigerMarsiliSigmaPartialCharges(molecule, true);
        } catch (Exception e) {
            throw new CDKException("Could not calculate partial charges: " + e.getMessage(), e);
        }
        counter = 0;
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            if (molecule.getAtom(i).getSymbol().equals("H")) continue;
            diagvalue[counter] = molecule.getAtom(i).getCharge();
            counter++;
        }
        burdenMatrix = BurdenMatrix.evalMatrix(molecule, diagvalue);
        matrix = new Matrix(burdenMatrix);
        eigenDecomposition = new EigenvalueDecomposition(matrix);
        double[] eval2 = eigenDecomposition.getRealEigenvalues();

        // get polarizability weighted BCUT
        Polarizability pol = new Polarizability();
        counter = 0;
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            if (molecule.getAtom(i).getSymbol().equals("H")) continue;
            diagvalue[counter] = pol.calculateGHEffectiveAtomPolarizability(molecule, molecule.getAtom(i), 1000);
            counter++;
        }
        burdenMatrix = BurdenMatrix.evalMatrix(molecule, diagvalue);
        matrix = new Matrix(burdenMatrix);
        eigenDecomposition = new EigenvalueDecomposition(matrix);
        double[] eval3 = eigenDecomposition.getRealEigenvalues();

        DoubleArrayResult retval = new DoubleArrayResult(eval1.length + eval2.length + eval3.length);

        String[] names;
        String[] suffix = {"w", "c", "p"};

        if (nhigh == 0 || nlow == 0) {
            // return all calculated eigenvalues
            for (int i = 0; i < eval1.length; i++) retval.add(eval1[i]);
            for (int i = 0; i < eval2.length; i++) retval.add(eval2[i]);
            for (int i = 0; i < eval3.length; i++) retval.add(eval3[i]);

            names = new String[retval.length()];
            for (int j = 0; j < suffix.length; j++) {
                for (int i = 0; i < eval1.length; i++) {
                    names[i] = "BCUT" + suffix[j] + "-" + (i + 1) + "l";
                }
                for (int i = 0; i < eval1.length; i++) {
                    names[i] = "BCUT" + suffix[j] + "-" + (i + 1) + "h";
                }
            }
        } else {
            // return only the n highest & lowest eigenvalues
            for (int i = 0; i < nlow; i++) retval.add(eval1[i]);
            for (int i = 0; i < nhigh; i++) retval.add(eval1[eval1.length - i - 1]);

            for (int i = 0; i < nlow; i++) retval.add(eval2[i]);
            for (int i = 0; i < nhigh; i++) retval.add(eval2[eval2.length - i - 1]);

            for (int i = 0; i < nlow; i++) retval.add(eval3[i]);
            for (int i = 0; i < nhigh; i++) retval.add(eval3[eval3.length - i - 1]);

            names = new String[3 * nhigh + 3 * nlow];
            counter = 0;
            for (int j = 0; j < suffix.length; j++) {
                for (int i = 0; i < nhigh; i++) {
                    names[counter++] = "BCUT" + suffix[j] + "-" + (i + 1) + "l";
                }
                for (int i = 0; i < nlow; i++) {
                    names[counter++] = "BCUT" + suffix[j] + "-" + (i + 1) + "h";
                }
            }
        }


        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval, names);
    }

    /**
     * Returns the specific type of the DescriptorResult object.
     * <p/>
     * The return value from this method really indicates what type of result will
     * be obtained from the {@link org.openscience.cdk.qsar.DescriptorValue} object. Note that the same result
     * can be achieved by interrogating the {@link org.openscience.cdk.qsar.DescriptorValue} object; this method
     * allows you to do the same thing, without actually calculating the descriptor.
     *
     * @return an object that implements the {@link org.openscience.cdk.qsar.result.IDescriptorResult} interface indicating
     *         the actual type of values returned by the descriptor in the {@link org.openscience.cdk.qsar.DescriptorValue} object
     */
    public IDescriptorResult getDescriptorResultType() {
        return new DoubleArrayResult();
    }
}