1 /* Copyright (C) 2008 Miguel Rojas <miguel.rojas@uni-koeln.de> 2 * 3 * Contact: cdk-devel@lists.sourceforge.net 4 * 5 * This program is free software; you can redistribute it and/or 6 * modify it under the terms of the GNU Lesser General Public License 7 * as published by the Free Software Foundation; either version 2.1 8 * of the License, or (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU Lesser General Public License for more details. 14 * 15 * You should have received a copy of the GNU Lesser General Public License 16 * along with this program; if not, write to the Free Software 17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 18 */ 19 package org.openscience.cdk.tools; 20 21 import java.util.ArrayList; 22 import java.util.List; 23 24 import org.openscience.cdk.CDKConstants; 25 import org.openscience.cdk.charges.Electronegativity; 26 import org.openscience.cdk.charges.GasteigerMarsiliPartialCharges; 27 import org.openscience.cdk.charges.GasteigerPEPEPartialCharges; 28 import org.openscience.cdk.charges.PiElectronegativity; 29 import org.openscience.cdk.charges.Polarizability; 30 import org.openscience.cdk.charges.StabilizationCharges; 31 import org.openscience.cdk.exception.CDKException; 32 import org.openscience.cdk.graph.Cycles; 33 import org.openscience.cdk.interfaces.IAtom; 34 import org.openscience.cdk.interfaces.IAtomContainer; 35 import org.openscience.cdk.interfaces.IAtomContainerSet; 36 import org.openscience.cdk.interfaces.IBond; 37 import org.openscience.cdk.interfaces.IReactionSet; 38 import org.openscience.cdk.interfaces.IRingSet; 39 import org.openscience.cdk.reaction.IReactionProcess; 40 import org.openscience.cdk.reaction.type.ElectronImpactNBEReaction; 41 import org.openscience.cdk.reaction.type.parameters.IParameterReact; 42 import org.openscience.cdk.reaction.type.parameters.SetReactionCenter; 43 import org.openscience.cdk.tools.manipulator.RingSetManipulator; 44 45 /** 46 * <p>This class contains the necessary information to predict ionization 47 * potential energy. It contains the models and the families classification. 48 * It is used as IPAtomicLearningDescriptor as the IPMolecularLearningDescriptor. 49 * <p> 50 * 51 * @author Miguel Rojas 52 * @cdk.created 2008-5-15 53 * @cdk.module ionpot 54 * 55 * @see org.openscience.cdk.qsar.descriptors.atomic.IPAtomicLearningDescriptor 56 * @see org.openscience.cdk.qsar.descriptors.molecular.IPMolecularLearningDescriptor 57 * @deprecated Results depend on the order of atoms in bonds, resolving this would require 58 * retraining the models which are not available. 59 */ 60 @Deprecated 61 public class IonizationPotentialTool { 62 63 /** 64 * Method which is predict the Ionization Potential from given atom. 65 * 66 * @param container The IAtomContainer where is contained the IAtom 67 * @param atom The IAtom to prediction the IP 68 * @return The value in eV 69 */ predictIP(IAtomContainer container, IAtom atom)70 public static double predictIP(IAtomContainer container, IAtom atom) throws CDKException { 71 double value = 0; 72 73 // at least one lone pair orbital is necessary to ionize 74 if (container.getConnectedLonePairsCount(atom) == 0) return value; 75 76 // control if the IAtom belongs in some family 77 if (familyHalogen(atom)) 78 value = getDTHalogenF(getQSARs(container, atom)); 79 else if (familyOxygen(atom)) 80 value = getDTOxygenF(getQSARs(container, atom)); 81 else if (familyNitrogen(atom)) value = getDTNitrogenF(getQSARs(container, atom)); 82 83 return value; 84 } 85 86 /** 87 * Method which is predict the Ionization Potential from given atom. 88 * 89 * @param container The IAtomContainer where is contained the IAtom 90 * @param bond The IBond to prediction the IP 91 * @return The value in eV 92 */ predictIP(IAtomContainer container, IBond bond)93 public static double predictIP(IAtomContainer container, IBond bond) throws CDKException { 94 double value = 0; 95 96 if (bond.getOrder().equals(IBond.Order.SINGLE)) return value; 97 98 //if some of the atoms belongs to some of the heteroatom family than 99 // it can not ionized 100 for (int i = 0; i < 2; i++) { 101 IAtom atom = bond.getAtom(i); 102 if (familyHalogen(atom)) 103 return value; 104 else if (familyOxygen(atom)) 105 return value; 106 else if (familyNitrogen(atom)) return value; 107 } 108 109 if (!familyBond(container, bond)) return value; 110 111 return getDTBondF(getQSARs(container, bond)); 112 } 113 114 /** 115 * Looking if the IAtom belongs to the halogen family. 116 * The IAtoms are F, Cl, Br, I. 117 * 118 * @param atom The IAtom 119 * @return True, if it belongs 120 */ familyHalogen(IAtom atom)121 private static boolean familyHalogen(IAtom atom) { 122 String symbol = atom.getSymbol(); 123 if (symbol.equals("F") || symbol.equals("Cl") || symbol.equals("Br") || symbol.equals("I")) 124 return true; 125 else 126 return false; 127 } 128 129 /** 130 * Looking if the Atom belongs to the oxygen family. 131 * The IAtoms are O, S, Se, Te. 132 * 133 * @param atom The IAtom 134 * @return True, if it belongs 135 */ familyOxygen(IAtom atom)136 private static boolean familyOxygen(IAtom atom) { 137 String symbol = atom.getSymbol(); 138 if (symbol.equals("O") || symbol.equals("S") || symbol.equals("Se") || symbol.equals("Te")) 139 return true; 140 else 141 return false; 142 } 143 144 /** 145 * Looking if the Atom belongs to the nitrogen family. 146 * The IAtoms are N, P, As, Sb. 147 * 148 * @param atom The IAtom 149 * @return True, if it belongs 150 */ familyNitrogen(IAtom atom)151 private static boolean familyNitrogen(IAtom atom) { 152 String symbol = atom.getSymbol(); 153 if (symbol.equals("N") || symbol.equals("P") || symbol.equals("As") || symbol.equals("Sb")) 154 return true; 155 else 156 return false; 157 } 158 159 /** 160 * Looking if the Bond belongs to the bond family. 161 * Not in resosance with other heteroatoms. 162 * 163 * @param container The IAtomContainer 164 * @param bond The IBond 165 * @return True, if it belongs 166 */ familyBond(IAtomContainer container, IBond bond)167 private static boolean familyBond(IAtomContainer container, IBond bond) { 168 List<String> normalAt = new ArrayList<String>(); 169 normalAt.add("C"); 170 normalAt.add("H"); 171 172 if (getDoubleBondNumber(container) > 30) // taking to long 173 return false; 174 175 StructureResonanceGenerator gRN = new StructureResonanceGenerator(); 176 IAtomContainer ac = gRN.getContainer(container, bond); 177 178 if (ac == null) return true; 179 180 if (getDoubleBondNumber(ac) > 15) // taking to long 181 return false; 182 183 for (IAtom atom : container.atoms()) { 184 if (!normalAt.contains(atom.getSymbol())) if (ac.contains(atom)) return false; 185 } 186 187 return true; 188 } 189 190 /** 191 * Extract the number of bond with superior ordre. 192 * 193 * @param container The IAtomContainer 194 * @return The number 195 */ getDoubleBondNumber(IAtomContainer container)196 private static int getDoubleBondNumber(IAtomContainer container) { 197 int doubleNumber = 0; 198 for (IBond bond : container.bonds()) { 199 if (bond.getOrder().equals(IBond.Order.DOUBLE) || bond.getOrder().equals(IBond.Order.TRIPLE)) 200 doubleNumber++; 201 } 202 return doubleNumber; 203 } 204 205 /** 206 * Get the results of 7 qsar descriptors been applied. They are: 207 * Electronegativity, 208 * GasteigerMarsiliPartialCharges, 209 * GasteigerPEPEPartialCharges, 210 * Polarizability, 211 * StabilizationCharge, 212 * Number of Atom in resonance 213 * if the container in resonance is aromatic. 214 * 215 * @param container The IAtomContainer which contain the IAtom 216 * @param atom The IAtom to calculate 217 * @return An Array containing the results 218 */ getQSARs(IAtomContainer container, IAtom atom)219 public static double[] getQSARs(IAtomContainer container, IAtom atom) throws CDKException { 220 Electronegativity electronegativity = new Electronegativity(); 221 PiElectronegativity pielectronegativity = new PiElectronegativity(); 222 GasteigerMarsiliPartialCharges peoe = new GasteigerMarsiliPartialCharges(); 223 GasteigerPEPEPartialCharges pepe = new GasteigerPEPEPartialCharges(); 224 Polarizability pol = new Polarizability(); 225 StabilizationCharges stabil = new StabilizationCharges(); 226 StructureResonanceGenerator gRI = new StructureResonanceGenerator(); 227 228 IAtomContainer product = initiateIonization(container, atom); 229 230 double[] results = new double[8]; 231 // sigmaElectronegativity 232 results[0] = electronegativity.calculateSigmaElectronegativity(container, atom); 233 // piElectronegativity 234 results[1] = pielectronegativity.calculatePiElectronegativity(container, atom); 235 // partialSigmaCharge 236 try { 237 peoe.assignGasteigerMarsiliSigmaPartialCharges(container, true); 238 } catch (Exception e) { 239 // ignored, underlying classes are logging this and this class 240 // is deprecated 241 } 242 results[2] = atom.getCharge(); 243 // partialPiCharge 244 for (int i = 0; i < container.getAtomCount(); i++) 245 container.getAtom(i).setCharge(0.0); 246 try { 247 pepe.assignGasteigerPiPartialCharges(container, true); 248 } catch (Exception e) { 249 e.printStackTrace(); 250 } 251 results[3] = atom.getCharge(); 252 // effectiveAtomicPolarizability 253 results[4] = pol.calculateGHEffectiveAtomPolarizability(container, atom, 100, true); 254 255 int position = container.indexOf(atom); 256 if (product != null) 257 results[5] = stabil.calculatePositive(product, product.getAtom(position)); 258 else 259 results[5] = 0.0; 260 // numberResonance 261 IAtomContainer acR = gRI.getContainer(container, atom); 262 if (acR != null) { 263 results[6] = acR.getAtomCount(); 264 // numberAromaticAtoms 265 // boolean isAromatic = Aromaticity.cdkLegacy().apply(container); 266 IRingSet ringSet = Cycles.sssr(container).toRingSet(); 267 RingSetManipulator.markAromaticRings(ringSet); 268 int aromRingCount = 0; 269 for (IAtomContainer ring : ringSet.atomContainers()) { 270 if (ring.getFlag(CDKConstants.ISAROMATIC)) aromRingCount++; 271 } 272 results[7] = aromRingCount; 273 } else { 274 results[6] = 0; 275 results[7] = 0; 276 } 277 278 return results; 279 } 280 281 /** 282 * Get the results of 7 qsar descriptors been applied. They are: 283 * Electronegativity, 284 * GasteigerMarsiliPartialCharges, 285 * GasteigerPEPEPartialCharges, 286 * Polarizability, 287 * StabilizationCharge, 288 * Number of Atom in resonance 289 * if the container in resonance is aromatic. 290 * 291 * @param container The IAtomContainer which contain the IAtom 292 * @param bond The IBond to calculate 293 * @return An Array containing the results 294 */ getQSARs(IAtomContainer container, IBond bond)295 public static double[] getQSARs(IAtomContainer container, IBond bond) throws CDKException { 296 Electronegativity electronegativity = new Electronegativity(); 297 PiElectronegativity pielectronegativity = new PiElectronegativity(); 298 GasteigerMarsiliPartialCharges peoe = new GasteigerMarsiliPartialCharges(); 299 GasteigerPEPEPartialCharges pepe = new GasteigerPEPEPartialCharges(); 300 Polarizability pol = new Polarizability(); 301 StabilizationCharges stabil = new StabilizationCharges(); 302 StructureResonanceGenerator gRI = new StructureResonanceGenerator(); 303 304 double[] results = new double[7]; 305 306 for (int ia = 0; ia < 2; ia++) { 307 IAtom atom = bond.getAtom(ia); 308 309 IAtomContainer product = initiateIonization(container, atom); 310 311 // sigmaElectronegativity 312 results[0] += electronegativity.calculateSigmaElectronegativity(container, atom); 313 // piElectronegativity 314 results[1] += pielectronegativity.calculatePiElectronegativity(container, atom); 315 // partialSigmaCharge 316 try { 317 peoe.assignGasteigerMarsiliSigmaPartialCharges(container, true); 318 } catch (Exception e) { 319 e.printStackTrace(); 320 } 321 results[2] += atom.getCharge(); 322 // partialPiCharge 323 for (int i = 0; i < container.getAtomCount(); i++) 324 container.getAtom(i).setCharge(0.0); 325 try { 326 pepe.assignGasteigerPiPartialCharges(container, true); 327 } catch (Exception e) { 328 e.printStackTrace(); 329 } 330 results[3] += atom.getCharge(); 331 // effectiveAtomicPolarizability 332 results[4] += pol.calculateGHEffectiveAtomPolarizability(container, atom, 100, true); 333 334 int position = container.indexOf(atom); 335 if (product != null) 336 results[5] += stabil.calculatePositive(product, product.getAtom(position)); 337 else 338 results[5] += 0.0; 339 // numberResonance 340 IAtomContainer acR = gRI.getContainer(container, atom); 341 if (acR != null) { 342 results[6] += acR.getAtomCount(); 343 // numberAromaticAtoms 344 // boolean isAromatic = Aromaticity.cdkLegacy().apply(container); 345 // if(isAromatic) 346 // results[7] += 0.1; 347 } else { 348 results[6] += 0; 349 // results[7] += 0; 350 } 351 } 352 for (int i = 0; i < results.length; i++) 353 results[i] = results[i] / 2; 354 return results; 355 } 356 357 /** 358 * Get the prediction result for the Halogen family given a series of values. 359 * It is based on 167 instances and 9 attributes(descriptors) using the Linear Regression Model 360 * with result of Root mean squared error 0.5817 with a cross validation of 10 folds. 361 * 362 * @param resultsH Array which contains the results of each descriptor 363 * @return The result 364 */ getDTHalogenF(double[] resultsH)365 private static double getDTHalogenF(double[] resultsH) { 366 double result = 0.0; 367 double SE = resultsH[0]; 368 double PSC = resultsH[2]; 369 double PIC = resultsH[3]; 370 double ETP = resultsH[4]; 371 double SPC = resultsH[5]; 372 double COUNTR = resultsH[6]; 373 double COUNTAr = resultsH[7]; 374 375 // System.out.println("SE : "+SE+", PE : "+PE+", PSC : "+PSC+", PIC : "+PIC+", ETP : "+ETP+", SPC : "+SPC+", COUNTR : "+COUNTR+", COUNTAr : "+COUNTAr); 376 //model leastMedSq 377 result = 0.272 * SE + 13.5814 * PSC + -4.4765 * PIC + -0.4937 * ETP + 0.0095 * SPC + -0.3706 * COUNTR + 0.5172 378 * COUNTAr + 12.4183; 379 return result; 380 } 381 382 /** 383 * Get the prediction result for the Oxygen family given a series of values. 384 * It is based on 368 instances and 9 attributes(descriptors) using the Linear Regression Model 385 * with result of Root mean squared error 0.64 with a cross validation of 10 folds. 386 * 387 * @param resultsH Array which contains the results of each descriptor 388 * @return The result 389 */ getDTOxygenF(double[] resultsH)390 private static double getDTOxygenF(double[] resultsH) { 391 double result = 0.0; 392 double SE = resultsH[0]; 393 double PE = resultsH[1]; 394 double PSC = resultsH[2]; 395 double PIC = resultsH[3]; 396 double ETP = resultsH[4]; 397 double SPC = resultsH[5]; 398 double COUNTR = resultsH[6]; 399 400 // System.out.println("SE : "+SE+", PE : "+PE+", PSC : "+PSC+", PIC : "+PIC+", ETP : "+ETP+", SPC : "+SPC+", COUNTR : "+COUNTR+", COUNTAr : "+COUNTAr); 401 result = -0.0118 * SE - 0.1859 * PE - 0.0752 * PSC - 8.1697 * PIC - 0.2278 * ETP - 0.0041 * SPC + 0.0175 402 * COUNTR + 11.4835; 403 404 return result; 405 } 406 407 /** 408 * Get the prediction result for the Nitrogen family given a series of values. 409 * It is based on 244 instances and 9 attributes(descriptors) using the Linear Regression Model 410 * with result of Root mean squared error 0.54 with a cross validation of 10 folds. 411 * 412 * @param resultsH Array which contains the results of each descriptor 413 * @return The result 414 */ getDTNitrogenF(double[] resultsH)415 private static double getDTNitrogenF(double[] resultsH) { 416 double result = 0.0; 417 double SE = resultsH[0]; 418 double PE = resultsH[1]; 419 double PSC = resultsH[2]; 420 double PIC = resultsH[3]; 421 double ETP = resultsH[4]; 422 double SPC = resultsH[5]; 423 double COUNTR = resultsH[6]; 424 425 // System.out.println("SE : "+SE+", PE : "+PE+", PSC : "+PSC+", PIC : "+PIC+", ETP : "+ETP+", SPC : "+SPC+", COUNTR : "+COUNTR+", COUNTAr : "+COUNTAr); 426 result = 0.4634 * SE + 0.0201 * PE + 1.1897 * PSC - 3.598 * PIC - 0.2726 * ETP + 0.0006 * SPC - 0.0527 * COUNTR 427 + 6.5419; 428 return result; 429 } 430 431 /** 432 * Get the desicion-tree result for the Halogen family given a series of values. 433 * It is based in 6 qsar descriptors. 434 * 435 * @param resultsH Array which contains the results of each descriptor 436 * @return The result 437 */ getDTBondF(double[] resultsH)438 private static double getDTBondF(double[] resultsH) { 439 double result = 0.0; 440 double SE = resultsH[0]; 441 double PE = resultsH[1]; 442 double PSC = resultsH[2]; 443 double PIC = resultsH[3]; 444 double ETP = resultsH[4]; 445 double COUNTR = resultsH[6]; 446 447 // System.out.println("SE : "+SE+", PE : "+PE+", PSC : "+PSC+", PIC : "+PIC+", ETP : "+ETP+", SPC : "+SPC+", COUNTR : "+COUNTR); 448 result = 0.1691 * SE + 1.1536 * PE + -6.3049 * PSC + -15.2638 * PIC + -0.2456 * ETP + -0.0139 * COUNTR + 2.114; 449 450 return result; 451 } 452 453 /** 454 * Initiate the reaction ElectronImpactNBE. 455 * 456 * @param container The IAtomContainer 457 * @param atom The IAtom to ionize 458 * @return The product resultant 459 * @throws CDKException 460 */ initiateIonization(IAtomContainer container, IAtom atom)461 private static IAtomContainer initiateIonization(IAtomContainer container, IAtom atom) throws CDKException { 462 IReactionProcess reactionNBE = new ElectronImpactNBEReaction(); 463 464 IAtomContainerSet setOfReactants = container.getBuilder().newInstance(IAtomContainerSet.class); 465 setOfReactants.addAtomContainer(container); 466 467 atom.setFlag(CDKConstants.REACTIVE_CENTER, true); 468 List<IParameterReact> paramList = new ArrayList<IParameterReact>(); 469 IParameterReact param = new SetReactionCenter(); 470 param.setParameter(Boolean.TRUE); 471 paramList.add(param); 472 reactionNBE.setParameterList(paramList); 473 474 /* initiate */ 475 IReactionSet setOfReactions = reactionNBE.initiate(setOfReactants, null); 476 atom.setFlag(CDKConstants.REACTIVE_CENTER, false); 477 if (setOfReactions != null && setOfReactions.getReactionCount() == 1 478 && setOfReactions.getReaction(0).getProducts().getAtomContainerCount() == 1) 479 return setOfReactions.getReaction(0).getProducts().getAtomContainer(0); 480 else 481 return null; 482 } 483 484 } 485