← Back to branch summary

~ubuntu-branches/ubuntu/maverick/cdk/maverick

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/qsar/descriptors/molecular/MomentOfInertiaDescriptor.java

  • Committer: Bazaar Package Importer
  • Author(s): Paul Cager
  • Date: 2008-04-09 21:17:53 UTC
  • Revision ID: james.westby@ubuntu.com-20080409211753-46lmjw5z8mx5pd8d
Tags: upstream-1.0.2
Import upstream version 1.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/org/openscience/cdk/qsar/descriptors/molecular/MomentOfInertiaDescriptor.java
 
1
/*
 
2
 *  Copyright (C) 2004-2007  The Chemistry Development Kit (CDK) project
 
3
 *
 
4
 *  Contact: cdk-devel@lists.sourceforge.net
 
5
 *
 
6
 *  This program is free software; you can redistribute it and/or
 
7
 *  modify it under the terms of the GNU Lesser General Public License
 
8
 *  as published by the Free Software Foundation; either version 2.1
 
9
 *  of the License, or (at your option) any later version.
 
10
 *
 
11
 *  This program is distributed in the hope that it will be useful,
 
12
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
 *  GNU Lesser General Public License for more details.
 
15
 *
 
16
 *  You should have received a copy of the GNU Lesser General Public License
 
17
 *  along with this program; if not, write to the Free Software
 
18
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 
19
 */
 
20
package org.openscience.cdk.qsar.descriptors.molecular;
 
21
 
 
22
import Jama.EigenvalueDecomposition;
 
23
import Jama.Matrix;
 
24
import org.openscience.cdk.config.IsotopeFactory;
 
25
import org.openscience.cdk.exception.CDKException;
 
26
import org.openscience.cdk.geometry.GeometryToolsInternalCoordinates;
 
27
import org.openscience.cdk.interfaces.IAtomContainer;
 
28
import org.openscience.cdk.qsar.DescriptorSpecification;
 
29
import org.openscience.cdk.qsar.DescriptorValue;
 
30
import org.openscience.cdk.qsar.IMolecularDescriptor;
 
31
import org.openscience.cdk.qsar.result.DoubleArrayResult;
 
32
import org.openscience.cdk.qsar.result.IDescriptorResult;
 
33
import org.openscience.cdk.tools.LoggingTool;
 
34
import org.openscience.cdk.tools.MFAnalyser;
 
35
 
 
36
import javax.vecmath.Point3d;
 
37
 
 
38
/**
 
39
 * A descriptor that calculates the moment of inertia and radius of gyration.
 
40
 * Moment of inertia (MI) values characterize the mass distribution of a molecule.
 
41
 * Related to the MI values, ratios of the MI values along the three principal axes
 
42
 * are also well know modeling variables. This descriptor calculates the MI values
 
43
 * along the X, Y and Z axes as well as the ratio's X/Y, X/Z and Y/Z. Finally it also
 
44
 * calculates the radius of gyration of the molecule.
 
45
 * <p/>
 
46
 * The descriptor generates 7 values in the following order
 
47
 * <ul>
 
48
 * <li>MOMI-X - MI along X axis
 
49
 * <li>MOMI-Y - MI along Y axis
 
50
 * <li>MOMI-Z - MI along Z axis
 
51
 * <li>MOMI-XY - X/Y
 
52
 * <li>MOMI-XZ - X/Z
 
53
 * <li>MOMI-YZ Y/Z
 
54
 * <li>MOMI-R - Radius of gyration
 
55
 * </ul>
 
56
 * One important aspect of the algorithm is that if the eigenvalues of the MI tensor
 
57
 * are below 1e-3, then the ratio's are set to a default of 1000.
 
58
 * <p/>
 
59
 * <p>This descriptor uses these parameters:
 
60
 * <table border="1">
 
61
 * <tr>
 
62
 * <td>Name</td>
 
63
 * <td>Default</td>
 
64
 * <td>Description</td>
 
65
 * </tr>
 
66
 * <tr>
 
67
 * <td></td>
 
68
 * <td></td>
 
69
 * <td>no parameters</td>
 
70
 * </tr>
 
71
 * </table>
 
72
 *
 
73
 * @author           Rajarshi Guha
 
74
 * @cdk.created      2005-02-07
 
75
 * @cdk.builddepends Jama-1.0.1.jar
 
76
 * @cdk.depends      Jama-1.0.1.jar
 
77
 * @cdk.module       qsar
 
78
 * @cdk.set          qsar-descriptors
 
79
 * @cdk.dictref      qsar-descriptors:momentOfInertia
 
80
 * @cdk.keyword      moment of inertia
 
81
 */
 
82
public class MomentOfInertiaDescriptor implements IMolecularDescriptor {
 
83
 
 
84
    private LoggingTool logger;
 
85
 
 
86
    public MomentOfInertiaDescriptor() {
 
87
        logger = new LoggingTool(this);
 
88
    }
 
89
 
 
90
    public DescriptorSpecification getSpecification() {
 
91
        return new DescriptorSpecification(
 
92
                "http://www.blueobelisk.org/ontologies/chemoinformatics-algorithms/#momentOfInertia",
 
93
                this.getClass().getName(),
 
94
                "$Id: MomentOfInertiaDescriptor.java 9056 2007-10-14 18:32:26Z egonw $",
 
95
                "The Chemistry Development Kit");
 
96
    }
 
97
 
 
98
    /**
 
99
     * Sets the parameters attribute of the MomentOfInertiaDescriptor object.
 
100
     *
 
101
     * @param params The new parameters value
 
102
     * @throws CDKException Description of the Exception
 
103
     * @see #getParameters
 
104
     */
 
105
    public void setParameters(Object[] params) throws CDKException {
 
106
        // no parameters for this descriptor
 
107
    }
 
108
 
 
109
    /**
 
110
     * Gets the parameters attribute of the MomentOfInertiaDescriptor object.
 
111
     *
 
112
     * @return The parameters value
 
113
     * @see #setParameters
 
114
     */
 
115
    public Object[] getParameters() {
 
116
        // no parameters to return
 
117
        return (null);
 
118
    }
 
119
 
 
120
    /**
 
121
     * Gets the parameterNames attribute of the MomentOfInertiaDescriptor object.
 
122
     *
 
123
     * @return The parameterNames value
 
124
     */
 
125
    public String[] getParameterNames() {
 
126
        // no param names to return
 
127
        return (null);
 
128
    }
 
129
 
 
130
 
 
131
    /**
 
132
     * Gets the parameterType attribute of the MomentOfInertiaDescriptor object.
 
133
     *
 
134
     * @param name Description of the Parameter
 
135
     * @return The parameterType value
 
136
     */
 
137
    public Object getParameterType(String name) {
 
138
        return (null);
 
139
    }
 
140
 
 
141
    /**
 
142
     * Calculates the 3 MI's, 3 ration and the R_gyr value.
 
143
     *
 
144
     * The molecule should have hydrogens
 
145
     *
 
146
     * @param container Parameter is the atom container.
 
147
     * @return An ArrayList containing 7 elements in the order described above
 
148
     * @throws CDKException if the supplied AtomContainer does not contain 3D coordinates
 
149
     */
 
150
 
 
151
    public DescriptorValue calculate(IAtomContainer container) throws CDKException {
 
152
        IsotopeFactory factory = null;
 
153
        try {
 
154
            factory = IsotopeFactory.getInstance(container.getBuilder());
 
155
        } catch (Exception e) {
 
156
            logger.debug(e);
 
157
        }
 
158
        factory.configureAtoms(container);
 
159
 
 
160
        DoubleArrayResult retval = new DoubleArrayResult(7);
 
161
 
 
162
        double ccf = 1.000138;
 
163
        double eps = 1e-5;
 
164
 
 
165
        double[][] imat = new double[3][3];
 
166
        Point3d centerOfMass = GeometryToolsInternalCoordinates.get3DCentreOfMass(container);
 
167
 
 
168
        double xdif;
 
169
        double ydif;
 
170
        double zdif;
 
171
        double xsq;
 
172
        double ysq;
 
173
        double zsq;
 
174
        for (int i = 0; i < container.getAtomCount(); i++) {
 
175
            org.openscience.cdk.interfaces.IAtom currentAtom = container.getAtom(i);
 
176
            if (currentAtom.getPoint3d() == null) {
 
177
                throw new CDKException("Atom " + i + " did not have any 3D coordinates. These are required");
 
178
            }
 
179
 
 
180
            double mass = factory.getMajorIsotope(currentAtom.getSymbol()).getExactMass();
 
181
 
 
182
            xdif = currentAtom.getPoint3d().x - centerOfMass.x;
 
183
            ydif = currentAtom.getPoint3d().y - centerOfMass.y;
 
184
            zdif = currentAtom.getPoint3d().z - centerOfMass.z;
 
185
            xsq = xdif * xdif;
 
186
            ysq = ydif * ydif;
 
187
            zsq = zdif * zdif;
 
188
 
 
189
            imat[0][0] += mass * (ysq + zsq);
 
190
            imat[1][1] += mass * (xsq + zsq);
 
191
            imat[2][2] += mass * (xsq + ysq);
 
192
 
 
193
            imat[1][0] += -1 * mass * ydif * xdif;
 
194
            imat[0][1] = imat[1][0];
 
195
 
 
196
            imat[2][0] += -1 * mass * xdif * zdif;
 
197
            imat[0][2] = imat[2][0];
 
198
 
 
199
            imat[2][1] += -1 * mass * ydif * zdif;
 
200
            imat[1][2] = imat[2][1];
 
201
        }
 
202
 
 
203
        // diagonalize the MI tensor
 
204
        Matrix tmp = new Matrix(imat);
 
205
        EigenvalueDecomposition eigenDecomp = tmp.eig();
 
206
        double[] eval = eigenDecomp.getRealEigenvalues();
 
207
 
 
208
        retval.add(eval[2]);
 
209
        retval.add(eval[1]);
 
210
        retval.add(eval[0]);
 
211
 
 
212
        double etmp = eval[0];
 
213
        eval[0] = eval[2];
 
214
        eval[2] = etmp;
 
215
 
 
216
        if (Math.abs(eval[1]) > 1e-3) retval.add(eval[0] / eval[1]);
 
217
        else retval.add(1000);
 
218
 
 
219
        if (Math.abs(eval[2]) > 1e-3) {
 
220
            retval.add(eval[0] / eval[2]);
 
221
            retval.add(eval[1] / eval[2]);
 
222
        } else {
 
223
            retval.add(1000);
 
224
            retval.add(1000);
 
225
        }
 
226
 
 
227
        // finally get the radius of gyration
 
228
        double pri = 0.0;
 
229
        MFAnalyser mfa = new MFAnalyser(container);
 
230
        if (Math.abs(eval[2]) > eps) pri = Math.pow(eval[0] * eval[1] * eval[2], 1.0 / 3.0);
 
231
        else pri = Math.sqrt(eval[0] * ccf / mfa.getMass());
 
232
        retval.add(Math.sqrt(Math.PI * 2 * pri * ccf / mfa.getMass()));
 
233
 
 
234
        String[] names = {
 
235
                "MOMI-X", "MOMI-Y", "MOMI-Z",
 
236
                "MOMI-XY", "MOMI-XZ", "MOMI-YZ", "MOMI-R"
 
237
        };
 
238
 
 
239
        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval, names);
 
240
    }
 
241
 
 
242
    /**
 
243
     * Returns the specific type of the DescriptorResult object.
 
244
     * <p/>
 
245
     * The return value from this method really indicates what type of result will
 
246
     * be obtained from the {@link org.openscience.cdk.qsar.DescriptorValue} object. Note that the same result
 
247
     * can be achieved by interrogating the {@link org.openscience.cdk.qsar.DescriptorValue} object; this method
 
248
     * allows you to do the same thing, without actually calculating the descriptor.
 
249
     *
 
250
     * @return an object that implements the {@link org.openscience.cdk.qsar.result.IDescriptorResult} interface indicating
 
251
     *         the actual type of values returned by the descriptor in the {@link org.openscience.cdk.qsar.DescriptorValue} object
 
252
     */
 
253
    public IDescriptorResult getDescriptorResultType() {
 
254
        return new DoubleArrayResult();
 
255
    }
 
256
}
 
257
    
 
258