1
/* $Revision: 7895 $ $Author: egonw $ $Date: 2007-02-07 23:20:44 +0100 (Wed, 07 Feb 2007) $
3
* Copyright (C) 2004-2007 Matteo Floris <mfe4@users.sf.net>
4
* Copyright (C) 2006-2007 Federico
6
* Contact: cdk-devel@lists.sourceforge.net
8
* This program is free software; you can redistribute it and/or
9
* modify it under the terms of the GNU Lesser General Public License
10
* as published by the Free Software Foundation; either version 2.1
11
* of the License, or (at your option) any later version.
13
* This program is distributed in the hope that it will be useful,
14
* but WITHOUT ANY WARRANTY; without even the implied warranty of
15
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16
* GNU Lesser General Public License for more details.
18
* You should have received a copy of the GNU Lesser General Public License
19
* along with this program; if not, write to the Free Software
20
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
22
package org.openscience.cdk.qsar.descriptors.atomic;
24
import org.openscience.cdk.AtomContainerSet;
25
import org.openscience.cdk.CDKConstants;
26
import org.openscience.cdk.Molecule;
27
import org.openscience.cdk.Ring;
28
import org.openscience.cdk.aromaticity.HueckelAromaticityDetector;
29
import org.openscience.cdk.charges.GasteigerMarsiliPartialCharges;
30
import org.openscience.cdk.exception.CDKException;
31
import org.openscience.cdk.graph.invariant.ConjugatedPiSystemsDetector;
32
import org.openscience.cdk.interfaces.IAtom;
33
import org.openscience.cdk.interfaces.IAtomContainer;
34
import org.openscience.cdk.interfaces.IBond;
35
import org.openscience.cdk.interfaces.IRingSet;
36
import org.openscience.cdk.qsar.DescriptorSpecification;
37
import org.openscience.cdk.qsar.DescriptorValue;
38
import org.openscience.cdk.qsar.IAtomicDescriptor;
39
import org.openscience.cdk.qsar.result.DoubleArrayResult;
40
import org.openscience.cdk.ringsearch.AllRingsFinder;
41
import org.openscience.cdk.tools.LoggingTool;
43
import javax.vecmath.Point3d;
44
import javax.vecmath.Vector3d;
45
import java.util.ArrayList;
46
import java.util.Iterator;
47
import java.util.List;
50
* This class calculates G3R proton descriptors used in neural networks for H1
54
* This descriptor uses these parameters: <table border="1">
58
* <td>Description</td>
61
* <td>checkAromaticity</td>
63
* <td>True is the aromaticity has to be checked</td>
68
* @cdk.created 2006-12-11
70
* @cdk.set qsar-descriptors
71
* @cdk.dictref qsar-descriptors:rdfProtonCalculatedValues
74
public class RDFProtonDescriptor_G3R implements IAtomicDescriptor {
76
private boolean checkAromaticity = false;
78
private IAtomContainer acold = null;
80
private IRingSet varRingSet = null;
82
private AtomContainerSet varAtomContainerSet = null;
84
private final static LoggingTool logger = new LoggingTool(RDFProtonDescriptor_G3R.class);
86
private static String[] descriptorNames;
89
* Constructor for the RDFProtonDescriptor object
91
public RDFProtonDescriptor_G3R() {
95
* Gets the specification attribute of the RDFProtonDescriptor_G3R object
97
* @return The specification value
99
public DescriptorSpecification getSpecification() {
100
return new DescriptorSpecification(
101
"http://www.blueobelisk.org/ontologies/chemoinformatics-algorithms/#rdfProtonCalculatedValues",
102
this.getClass().getName(),
103
"$Id: RDFProtonDescriptor.java 7032 2006-09-22 15:26:48 +0000 (ven, 22 set 2006) kaihartmann $",
104
"The Chemistry Development Kit");
108
* Sets the parameters attribute of the RDFProtonDescriptor object
111
* Parameters are the proton position and a boolean (true if you
112
* need to detect aromaticity)
113
* @exception CDKException
114
* Possible Exceptions
116
public void setParameters(Object[] params) throws CDKException {
117
if (params.length > 1) {
118
throw new CDKException(
119
"RDFProtonDescriptor only expects one parameters");
121
if (!(params[0] instanceof Boolean)) {
122
throw new CDKException(
123
"The second parameter must be of type Boolean");
125
checkAromaticity = ((Boolean) params[0]).booleanValue();
129
* Gets the parameters attribute of the RDFProtonDescriptor object
131
* @return The parameters value
133
public Object[] getParameters() {
134
// return the parameters as used for the descriptor calculation
135
Object[] params = new Object[1];
136
params[0] = Boolean.valueOf(checkAromaticity);
140
public DescriptorValue calculate(IAtom atom,
141
IAtomContainer varAtomContainerSet) throws CDKException {
142
return (calculate(atom, varAtomContainerSet, null));
145
public DescriptorValue calculate(IAtom atom,
146
IAtomContainer atomContainer, IRingSet precalculatedringset)
147
throws CDKException {
149
IAtomContainer varAtomContainer;
151
varAtomContainer = (IAtomContainer) atomContainer.clone();
152
} catch (CloneNotSupportedException e) {
153
throw new CDKException("Error during clone");
156
int atomPosition = atomContainer.getAtomNumber(atom);
157
IAtom clonedAtom = varAtomContainer.getAtom(atomPosition);
159
final int g3r_desc_length = 13;
161
DoubleArrayResult rdfProtonCalculatedValues = new DoubleArrayResult(
163
if (!atom.getSymbol().equals("H")) {
164
throw new CDKException("You tried calculation on a "
166
+ " atom. This is not allowed! Atom must be a H atom.");
169
// ///////////////////////FIRST SECTION OF MAIN METHOD: DEFINITION OF
171
// ///////////////////////AND AROMATICITY AND PI-SYSTEM AND RINGS
174
Molecule mol = new Molecule(varAtomContainer);
175
if (varAtomContainer != acold) {
176
acold = varAtomContainer;
177
// DETECTION OF pi SYSTEMS
178
varAtomContainerSet = ConjugatedPiSystemsDetector.detect(mol);
179
if (precalculatedringset == null)
180
varRingSet = (new AllRingsFinder())
181
.findAllRings(varAtomContainer);
183
varRingSet = precalculatedringset;
185
GasteigerMarsiliPartialCharges peoe = new GasteigerMarsiliPartialCharges();
186
peoe.assignGasteigerMarsiliSigmaPartialCharges(mol, true);
187
} catch (Exception ex1) {
188
throw new CDKException(
189
"Problems with assignGasteigerMarsiliPartialCharges due to "
190
+ ex1.toString(), ex1);
193
if (checkAromaticity) {
194
HueckelAromaticityDetector.detectAromaticity(varAtomContainer,
199
List ringsWithThisBond;
200
// SET ISINRING FLAGS FOR BONDS
201
// org.openscience.cdk.interfaces.IBond[] bondsInContainer = varAtomContainer.getBonds();
203
Iterator bondsInContainer = varAtomContainer.bonds();
204
while (bondsInContainer.hasNext()) {
205
IBond bond = (IBond) bondsInContainer.next();
206
ringsWithThisBond = varRingSet.getRings(bond);
207
if (ringsWithThisBond.size() > 0) {
208
bond.setFlag(CDKConstants.ISINRING, true);
211
// SET ISINRING FLAGS FOR ATOMS
212
org.openscience.cdk.interfaces.IRingSet ringsWithThisAtom;
214
for (int w = 0; w < varAtomContainer.getAtomCount(); w++) {
215
ringsWithThisAtom = varRingSet
216
.getRings(varAtomContainer.getAtom(w));
217
if (ringsWithThisAtom.getAtomContainerCount() > 0) {
218
varAtomContainer.getAtom(w)
219
.setFlag(CDKConstants.ISINRING, true);
223
IAtomContainer detected = varAtomContainerSet.getAtomContainer(0);
225
// neighboors[0] is the atom joined to the target proton:
226
java.util.List neighboors = mol.getConnectedAtomsList(clonedAtom);
227
IAtom neighbour0 = (IAtom) neighboors.get(0);
229
// 2', 3', 4', 5', 6', and 7' atoms up to the target are detected:
230
List atomsInSecondSphere = mol.getConnectedAtomsList(neighbour0);
231
List atomsInThirdSphere = null;
232
List atomsInFourthSphere = null;
233
List atomsInFifthSphere = null;
234
List atomsInSixthSphere = null;
235
List atomsInSeventhSphere = null;
237
// SOME LISTS ARE CREATED FOR STORING OF INTERESTING ATOMS AND BONDS
239
ArrayList singles = new ArrayList(); // list of any bond not
241
ArrayList doubles = new ArrayList(); // list with only double bonds
242
ArrayList atoms = new ArrayList(); // list with all the atoms in
244
// atoms.add( new Integer( mol.getAtomNumber(neighboors[0]) ) );
245
ArrayList bondsInCycloex = new ArrayList(); // list for bonds in
246
// cycloexane-like rings
248
// 2', 3', 4', 5', 6', and 7' bonds up to the target are detected:
249
IBond secondBond; // (remember that first bond is proton bond)
254
IBond seventhBond; //
256
// definition of some variables used in the main FOR loop for detection
257
// of interesting atoms and bonds:
258
boolean theBondIsInA6MemberedRing; // this is like a flag for bonds
259
// which are in cycloexane-like
260
// rings (rings with more than 4
266
// THIS MAIN FOR LOOP DETECT RIGID BONDS IN 7 SPHERES:
267
for (int a = 0; a < atomsInSecondSphere.size(); a++) {
268
IAtom curAtomSecond = (IAtom) atomsInSecondSphere.get(a);
269
secondBond = mol.getBond(neighbour0, curAtomSecond);
270
if (mol.getAtomNumber(curAtomSecond) != atomPosition
271
&& getIfBondIsNotRotatable(mol, secondBond, detected)) {
273
bondOrder = secondBond.getOrder();
274
bondNumber = mol.getBondNumber(secondBond);
275
theBondIsInA6MemberedRing = false;
276
checkAndStore(bondNumber, bondOrder, singles, doubles,
277
bondsInCycloex, mol.getAtomNumber(curAtomSecond),
278
atoms, sphere, theBondIsInA6MemberedRing);
279
atomsInThirdSphere = mol.getConnectedAtomsList(curAtomSecond);
280
if (atomsInThirdSphere.size() > 0) {
281
for (int b = 0; b < atomsInThirdSphere.size(); b++) {
282
IAtom curAtomThird = (IAtom) atomsInThirdSphere.get(b);
283
thirdBond = mol.getBond(curAtomThird, curAtomSecond);
284
// IF THE ATOMS IS IN THE THIRD SPHERE AND IN A
285
// CYCLOEXANE-LIKE RING, IT IS STORED IN THE PROPER
287
if (mol.getAtomNumber(curAtomThird) != atomPosition
288
&& getIfBondIsNotRotatable(mol, thirdBond,
291
bondOrder = thirdBond.getOrder();
292
bondNumber = mol.getBondNumber(thirdBond);
293
theBondIsInA6MemberedRing = false;
295
// if the bond is in a cyclohexane-like ring (a ring
296
// with 5 or more atoms, not aromatic)
297
// the boolean "theBondIsInA6MemberedRing" is set to
299
if (!thirdBond.getFlag(CDKConstants.ISAROMATIC)) {
300
if (!curAtomThird.equals(neighbour0)) {
301
rsAtom = varRingSet.getRings(thirdBond);
302
for (int f = 0; f < rsAtom.size(); f++) {
303
ring = (Ring) rsAtom.get(f);
304
if (ring.getRingSize() > 4
305
&& ring.contains(thirdBond)) {
306
theBondIsInA6MemberedRing = true;
311
checkAndStore(bondNumber, bondOrder, singles,
312
doubles, bondsInCycloex, mol
313
.getAtomNumber(curAtomThird),
314
atoms, sphere, theBondIsInA6MemberedRing);
315
theBondIsInA6MemberedRing = false;
316
atomsInFourthSphere = mol
317
.getConnectedAtomsList(curAtomThird);
318
if (atomsInFourthSphere.size() > 0) {
319
for (int c = 0; c < atomsInFourthSphere.size(); c++) {
320
IAtom curAtomFourth = (IAtom) atomsInFourthSphere
322
fourthBond = mol.getBond(curAtomThird,
324
if (mol.getAtomNumber(curAtomFourth) != atomPosition
325
&& getIfBondIsNotRotatable(mol,
326
fourthBond, detected)) {
328
bondOrder = fourthBond.getOrder();
330
.getBondNumber(fourthBond);
331
theBondIsInA6MemberedRing = false;
339
.getAtomNumber(curAtomFourth),
341
theBondIsInA6MemberedRing);
342
atomsInFifthSphere = mol
343
.getConnectedAtomsList(curAtomFourth);
344
if (atomsInFifthSphere.size() > 0) {
345
for (int d = 0; d < atomsInFifthSphere
347
IAtom curAtomFifth = (IAtom) atomsInFifthSphere
349
fifthBond = mol.getBond(
353
.getAtomNumber(curAtomFifth) != atomPosition
354
&& getIfBondIsNotRotatable(
358
bondOrder = fifthBond
361
.getBondNumber(fifthBond);
362
theBondIsInA6MemberedRing = false;
370
.getAtomNumber(curAtomFifth),
372
theBondIsInA6MemberedRing);
373
atomsInSixthSphere = mol
374
.getConnectedAtomsList(curAtomFifth);
375
if (atomsInSixthSphere
377
for (int e = 0; e < atomsInSixthSphere
379
IAtom curAtomSixth = (IAtom) atomsInSixthSphere
386
.getAtomNumber(curAtomSixth) != atomPosition
387
&& getIfBondIsNotRotatable(
392
bondOrder = sixthBond
395
.getBondNumber(sixthBond);
396
theBondIsInA6MemberedRing = false;
404
.getAtomNumber(curAtomSixth),
407
theBondIsInA6MemberedRing);
408
atomsInSeventhSphere = mol
409
.getConnectedAtomsList(curAtomSixth);
410
if (atomsInSeventhSphere
412
for (int f = 0; f < atomsInSeventhSphere
414
IAtom curAtomSeventh = (IAtom) atomsInSeventhSphere
421
.getAtomNumber(curAtomSeventh) != atomPosition
422
&& getIfBondIsNotRotatable(
427
bondOrder = seventhBond
430
.getBondNumber(seventhBond);
431
theBondIsInA6MemberedRing = false;
439
.getAtomNumber(curAtomSeventh),
442
theBondIsInA6MemberedRing);
471
// ////////////////////////LAST DESCRIPTOR IS g3(r), FOR PROTONS BONDED
472
// TO LIKE-CYCLOEXANE RINGS:
474
Vector3d a_a = new Vector3d();
475
Vector3d a_b = new Vector3d();
476
Vector3d b_a = new Vector3d();
477
Vector3d b_b = new Vector3d();
480
if (bondsInCycloex.size() > 0) {
481
IAtom cycloexBondAtom0;
482
IAtom cycloexBondAtom1;
483
org.openscience.cdk.interfaces.IBond theInCycloexBond;
487
step = (limitSup - limitInf) / 13;
493
for (double g3r = 0; g3r < limitSup; g3r = g3r + step) {
495
for (int cyc = 0; cyc < bondsInCycloex.size(); cyc++) {
499
Integer thisInCycloexBond = (Integer) bondsInCycloex
501
position = thisInCycloexBond.intValue();
502
theInCycloexBond = mol.getBond(position);
503
cycloexBondAtom0 = theInCycloexBond.getAtom(0);
504
cycloexBondAtom1 = theInCycloexBond.getAtom(1);
506
connAtoms = mol.getConnectedAtomsList(cycloexBondAtom0);
507
for (int g = 0; g < connAtoms.size(); g++) {
508
if (((IAtom) connAtoms.get(g)).equals(neighbour0))
512
if (ya_counter > 0) {
513
a_a.set(cycloexBondAtom1.getPoint3d().x,
514
cycloexBondAtom1.getPoint3d().y,
515
cycloexBondAtom1.getPoint3d().z);
516
a_b.set(cycloexBondAtom0.getPoint3d().x,
517
cycloexBondAtom0.getPoint3d().y,
518
cycloexBondAtom0.getPoint3d().z);
520
a_a.set(cycloexBondAtom0.getPoint3d().x,
521
cycloexBondAtom0.getPoint3d().y,
522
cycloexBondAtom0.getPoint3d().z);
523
a_b.set(cycloexBondAtom1.getPoint3d().x,
524
cycloexBondAtom1.getPoint3d().y,
525
cycloexBondAtom1.getPoint3d().z);
527
b_a.set(neighbour0.getPoint3d().x,
528
neighbour0.getPoint3d().y,
529
neighbour0.getPoint3d().z);
530
b_b.set(atom.getPoint3d().x, atom.getPoint3d().y, atom
533
angle = calculateAngleBetweenTwoLines(a_a, a_b, b_a, b_b);
535
// logger.debug("gcycr ANGLE: " + angle + " "
536
// +mol.getAtomNumber(cycloexBondAtom0) + "
537
// "+mol.getAtomNumber(cycloexBondAtom1));
539
partial = Math.exp(smooth * (Math.pow((g3r - angle), 2)));
542
// g3r_function.add(new Double(sum));
543
rdfProtonCalculatedValues.add(sum);
544
logger.debug("RDF g3r prob.: "+sum+ " at distance "+g3r);
548
for (int i=0; i<g3r_desc_length; i++) rdfProtonCalculatedValues.add(Double.NaN);
550
return new DescriptorValue(
551
getSpecification(), getParameterNames(),
552
getParameters(), rdfProtonCalculatedValues,
557
// Others definitions
559
private boolean getIfBondIsNotRotatable(Molecule mol,
560
org.openscience.cdk.interfaces.IBond bond, IAtomContainer detected) {
561
boolean isBondNotRotatable = false;
563
IAtom atom0 = bond.getAtom(0);
564
IAtom atom1 = bond.getAtom(1);
565
if (detected != null) {
566
if (detected.contains(bond))
569
if (atom0.getFlag(CDKConstants.ISINRING)) {
570
if (atom1.getFlag(CDKConstants.ISINRING)) {
573
if (atom1.getSymbol().equals("H"))
579
if (atom0.getSymbol().equals("N") && atom1.getSymbol().equals("C")) {
580
if (getIfACarbonIsDoubleBondedToAnOxygen(mol, atom1))
583
if (atom0.getSymbol().equals("C") && atom1.getSymbol().equals("N")) {
584
if (getIfACarbonIsDoubleBondedToAnOxygen(mol, atom0))
588
isBondNotRotatable = true;
589
return isBondNotRotatable;
592
private boolean getIfACarbonIsDoubleBondedToAnOxygen(Molecule mol,
594
boolean isDoubleBondedToOxygen = false;
595
java.util.List neighToCarbon = mol.getConnectedAtomsList(carbonAtom);
596
org.openscience.cdk.interfaces.IBond tmpBond;
598
for (int nei = 0; nei < neighToCarbon.size(); nei++) {
599
IAtom neighbour = (IAtom) neighToCarbon.get(nei);
600
if (neighbour.getSymbol().equals("O")) {
601
tmpBond = mol.getBond(neighbour, carbonAtom);
602
if (tmpBond.getOrder() == 2.0)
607
isDoubleBondedToOxygen = true;
608
return isDoubleBondedToOxygen;
611
// this method calculates the angle between two bonds given coordinates of
614
public double calculateAngleBetweenTwoLines(Vector3d a, Vector3d b,
615
Vector3d c, Vector3d d) {
616
Vector3d firstLine = new Vector3d();
618
Vector3d secondLine = new Vector3d();
619
secondLine.sub(c, d);
620
Vector3d firstVec = new Vector3d(firstLine);
621
Vector3d secondVec = new Vector3d(secondLine);
622
return firstVec.angle(secondVec);
625
// this method store atoms and bonds in proper lists:
626
private void checkAndStore(int bondToStore, double bondOrder,
627
ArrayList singleVec, ArrayList doubleVec, ArrayList cycloexVec,
628
int a1, ArrayList atomVec, int sphere, boolean isBondInCycloex) {
629
if (!atomVec.contains(new Integer(a1))) {
631
atomVec.add(new Integer(a1));
633
if (!cycloexVec.contains(new Integer(bondToStore))) {
634
if (isBondInCycloex) {
635
cycloexVec.add(new Integer(bondToStore));
638
if (bondOrder == 2.0) {
639
if (!doubleVec.contains(new Integer(bondToStore)))
640
doubleVec.add(new Integer(bondToStore));
642
if (bondOrder == 1.0) {
643
if (!singleVec.contains(new Integer(bondToStore)))
644
singleVec.add(new Integer(bondToStore));
648
// generic method for calculation of distance btw 2 atoms
649
private double calculateDistanceBetweenTwoAtoms(IAtom atom1, IAtom atom2) {
651
Point3d firstPoint = atom1.getPoint3d();
652
Point3d secondPoint = atom2.getPoint3d();
653
distance = firstPoint.distance(secondPoint);
657
// given a double bond
658
// this method returns a bond bonded to this double bond
659
private int getNearestBondtoAGivenAtom(Molecule mol, IAtom atom,
660
org.openscience.cdk.interfaces.IBond bond) {
664
IAtom atom0 = bond.getAtom(0);
665
IAtom atom1 = bond.getAtom(1);
666
List bondsAtLeft = mol.getConnectedBondsList(atom0);
668
for (int i = 0; i < bondsAtLeft.size(); i++) {
669
IBond curBond = (IBond) bondsAtLeft.get(i);
670
values = calculateDistanceBetweenAtomAndBond(atom, curBond);
671
partial = mol.getBondNumber(curBond);
673
nearestBond = mol.getBondNumber(curBond);
674
distance = values[0];
676
if (values[0] < distance) {
677
nearestBond = partial;
680
* XXX commented this out, because is has no effect
682
* else { nearestBond = nearestBond; }
690
// method which calculated distance btw an atom and the middle point of a
692
// and returns distance and coordinates of middle point
693
private double[] calculateDistanceBetweenAtomAndBond(IAtom proton,
694
org.openscience.cdk.interfaces.IBond theBond) {
695
Point3d middlePoint = theBond.get3DCenter();
696
Point3d protonPoint = proton.getPoint3d();
697
double[] values = new double[4];
698
values[0] = middlePoint.distance(protonPoint);
699
values[1] = middlePoint.x;
700
values[2] = middlePoint.y;
701
values[3] = middlePoint.z;
706
* Gets the parameterNames attribute of the RDFProtonDescriptor object
708
* @return The parameterNames value
710
public String[] getParameterNames() {
711
String[] params = new String[2];
712
params[0] = "atomPosition";
713
params[1] = "checkAromaticity";
718
* Gets the parameterType attribute of the RDFProtonDescriptor object
721
* Description of the Parameter
722
* @return The parameterType value
724
public Object getParameterType(String name) {
725
if (name.equals("atomPosition"))
726
return new Integer(0);
added
removed
src/org/openscience/cdk/qsar/descriptors/atomic/RDFProtonDescriptor_G3R.java
