1 /* $RCSfile$ 2 * $Author: hansonr $ 3 * $Date: 2007-11-23 12:49:25 -0600 (Fri, 23 Nov 2007) $ 4 * $Revision: 8655 $ 5 * 6 * Copyright (C) 2003-2005 The Jmol Development Team 7 * 8 * Contact: jmol-developers@lists.sf.net 9 * 10 * This library is free software; you can redistribute it and/or 11 * modify it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 23 */ 24 25 package org.jmol.minimize.forcefield; 26 27 import java.io.BufferedReader; 28 import java.util.Hashtable; 29 import java.util.Map; 30 31 import org.jmol.api.SmilesMatcherInterface; 32 import org.jmol.minimize.MinAngle; 33 import org.jmol.minimize.MinAtom; 34 import org.jmol.minimize.MinBond; 35 import org.jmol.minimize.MinObject; 36 import org.jmol.minimize.MinTorsion; 37 import org.jmol.minimize.Minimizer; 38 import org.jmol.modelset.Atom; 39 import org.jmol.modelset.Bond; 40 import org.jmol.util.BSUtil; 41 import org.jmol.util.Elements; 42 import org.jmol.util.Escape; 43 import org.jmol.util.Logger; 44 import org.jmol.viewer.JmolAsyncException; 45 46 import javajs.util.AU; 47 import javajs.util.BS; 48 import javajs.util.Lst; 49 import javajs.util.PT; 50 51 /** 52 * MMFF94 implementation 5/14/2012 53 * 54 * - fully validated for atom types and charges 55 * - reasonably well validated for energies (see below) 56 * 57 * - TODO: add UFF for preliminary/backup calculation 58 * 59 * @author Bob Hanson hansonr@stolaf.edu 60 * 61 * Java implementation by Bob Hanson 5/2012 62 * based loosely on chemKit code by Kyle Lutz and OpenBabel code by Tim Vandermeersch 63 * but primarily from what is described in 64 * 65 * T. A. Halgren; "Merck Molecular Force Field. V. Extension of MMFF94 66 * Using Experimental Data, Additional Computational Data, 67 * and Empirical Rules", J. Comp. Chem. 5 & 6 616-641 (1996). 68 * 69 * Parameter files are clipped from the original Wiley FTP site supplemental material: 70 * 71 * ftp://ftp.wiley.com/public/journals/jcc/suppmat/17/490/MMFF-I_AppendixB.ascii 72 * 73 * Original work, as listed at http://towhee.sourceforge.net/forcefields/mmff94.html: 74 * 75 * T. A. Halgren; "Merck Molecular Force Field. I. Basis, Form, Scope, 76 * Parameterization, and Performance of MMFF94", J. Comp. Chem. 5 & 6 490-519 (1996). 77 * T. A. Halgren; "Merck Molecular Force Field. II. MMFF94 van der Waals 78 * and Electrostatic Parameters for Intermolecular Interactions", 79 * J. Comp. Chem. 5 & 6 520-552 (1996). 80 * T. A. Halgren; "Merck Molecular Force Field. III. Molecular Geometries and 81 * Vibrational Frequencies for MMFF94", J. Comp. Chem. 5 & 6 553-586 (1996). 82 * T. A. Halgren; R. B. Nachbar; "Merck Molecular Force Field. IV. 83 * Conformational Energies and Geometries for MMFF94", J. Comp. Chem. 5 & 6 587-615 (1996). 84 * T. A. Halgren; "Merck Molecular Force Field. V. Extension of MMFF94 85 * Using Experimental Data, Additional Computational Data, 86 * and Empirical Rules", J. Comp. Chem. 5 & 6 616-641 (1996). 87 * T. A. Halgren; "MMFF VII. Characterization of MMFF94, MMFF94s, 88 * and Other Widely Available Force Fields for Conformational Energies 89 * and for Intermolecular-Interaction Energies and Geometries", 90 * J. Comp. Chem. 7 730-748 (1999). 91 * 92 * Validation carried out using MMFF94_opti.log and MMFF94_dative.mol2 (or MMFF94_hypervalent.mol2) 93 * including 761 models using org/jmol/minimize/forcefield/mmff/validate/checkmm.spt (checkAllEnergies) 94 * 95 * All typical compounds validate. The following 7 96 * structures do not validate to within 0.1 kcal/mol total energy; 97 * 98 99 version=12.3.26_dev 100 101 # code: adding empirical rules to MMFF94 calculation 102 # 103 # checkmm.spt;checkAllEnergies 104 # 105 # checking calculated energies for 761 models 106 # 1 COMKAQ E= -7.3250003 Eref= -7.6177 diff= 0.2926998 107 # 2 DUVHUX10 E= 64.759995 Eref= 64.082855 diff= 0.6771393 108 # 3 FORJIF E= 35.978 Eref= 35.833878 diff= 0.14412308 109 # 4 JADLIJ E= 25.104 Eref= 24.7038 diff= 0.4001999 110 # 5 PHOSLA10 E= 111.232994 Eref= 112.07078 diff= 0.8377838 111 # 6 PHOSLB10 E= -93.479004 Eref= -92.64081 diff= 0.8381958 112 # 113 # for 761 atoms, 6 have energy differences outside the range -0.1 to 0.1 114 # with a standard deviation of 0.05309403 115 # 116 # a comment about empirical bond parameter calculation: 117 # 118 # // Well, guess what? As far as I can tell, in Eqn 18 on page 625, 119 # // the reduction term and delta are zero. 120 # 121 # // -- at least in the program run that is at the validation site: 122 # // OPTIMOL: Molecular and Macromolecular Optimization Package 17-Nov-98 16:01:23 123 # // SGI double-precision version ... Updated 5/6/98 124 # // 125 # // This calculation is run only for the following three structures. In each case the 126 # // reported validation values and values from Jmol 12.3.26_dev are shown. Clearly 127 # // the r0 calculated and final energies are very good. subtracting off 0.008 from 128 # // r0 would certainly not give the reported values. Something is odd there. 129 # // 130 # // bond red* r0(here/valid) kb(here/valid) Etotal(here/valid) 131 # // --------------------------------------------------------------------------------------- 132 # // OHWM1 H1-O1 0.03 0.978/0.978 7.510/7.51 -21.727/-21.72690 133 # // ERULE_03 Si1-P1 0.0 2.223/2.224 1.614/1.609 -2.983/ -2.93518 134 # // ERULE_06 N1-F1 0.0 1.381/1.379 5.372/5.438 1.582/ 1.58172 135 # // 136 # // *reduction and delta terms not used in Jmol's calculation 137 # 138 # 139 140 141 COMKAQ 142 -- BATCHMIN ignores 1 of 5-membered ring torsions for a 1-oxo-2-oxa-bicyclo[3.2.0]heptane 143 -- MMFF94_bmin.log: WARNING - Conformational Energies May Not Be Accurate 144 145 DUVHUX10 146 -- BATCHMIN ignores 5-membered ring issue for S-S-containing ring 147 -- MMFF94_bmin.log: WARNING - Conformational Energies May Not Be Accurate 148 149 FORJIF 150 -- BATCHMIN misses four standard 5-membered C-C ring bonds 151 -- MMFF94_bmin.log: WARNING - Conformational Energies May Not Be Accurate 152 153 JADLIJ 154 -- BATCHMIN ignores 5-membered ring for S (note, however, this is not the case in BODKOU) 155 -- MMFF94_bmin.log: WARNING - Conformational Energies May Not Be Accurate 156 157 PHOSLA10 158 -- BATCHMIN ignores all 5-membered ring torsions in ring with P 159 -- (note, however, this is not the case in CUVGAB) 160 -- MMFF94_bmin.log: WARNING - Conformational Energies May Not Be Accurate 161 162 PHOSLB10 163 -- BATCHMIN ignores all 5-membered ring torsions in ring with P 164 -- (note, however, this is not the case in CUVGAB) 165 -- MMFF94_bmin.log: WARNING - Conformational Energies May Not Be Accurate 166 167 OHMW1 168 -- H2O complexed with hydroxide OH(-) 169 -- I don't understand (a) why the OH(-) bond has mltb=1, and even with that 170 I am not getting the correct ro/kb for that bond from empirical rules. 171 Still working on that.... 172 173 */ 174 175 176 public class ForceFieldMMFF extends ForceField { 177 178 179 //private boolean useEmpiricalRules = true; 180 181 private static final int A4_VDW = 122; 182 private static final int A4_BNDK = 123; 183 private static final int A4_CHRG = 124; 184 private static final int A4_SB = 125; 185 private static final int A4_SBDEF = 126; 186 private static final int KEY_SBDEF = 0; 187 private static final int KEY_PBCI = 0; 188 private static final int KEY_VDW = 0; 189 private static final int KEY_BNDK = 0; 190 private static final int KEY_OOP = 6; 191 private static final int TYPE_PBCI = 0x1; 192 private static final int TYPE_VDW = 0x11; 193 private static final int TYPE_BNDK = 0x222; 194 private static final int TYPE_CHRG = 0x22; 195 // the following are bit flags indicating in the 0xF range 196 // which atoms might have default parameter values 197 private static final int TYPE_BOND = 0x3; // 0011 198 private static final int TYPE_ANGLE = 0x5; // 0101 199 private static final int TYPE_SB = 0x15; // 01 0101 200 private static final int TYPE_SBDEF = 0x25; // 10 0101 201 private static final int TYPE_TORSION = 0x9; // 00 1001 202 private static final int TYPE_OOP = 0xD; // 00 1101; 203 204 205 private static Lst<AtomType> atomTypes; 206 private static Map<Object, Object> mmffParams; 207 protected Map<Object, Object> ffParams; 208 209 private int[] rawAtomTypes; 210 private int[] rawBondTypes; 211 private float[] rawMMFF94Charges; // calculated here 212 getAtomTypeDescriptions()213 public String[] getAtomTypeDescriptions() { 214 return getAtomTypeDescs(rawAtomTypes); 215 } 216 getPartialCharges()217 public float[] getPartialCharges() { 218 return rawMMFF94Charges; 219 } 220 221 /* 222 * from SMARTS search when calculating partial charges: 223 * 224 * vRings[0] list of 3-membered rings 225 * vRings[1] list of 4-membered rings 226 * vRings[2] list of all 5-membered rings 227 * vRings[3] list of aromatic 5-membered and 6-membered rings 228 */ 229 private Lst<BS>[] vRings; 230 private static Map<Object, Object> mmff2DParams; 231 ForceFieldMMFF(Minimizer m, boolean isQuick)232 public ForceFieldMMFF(Minimizer m, boolean isQuick) 233 throws JmolAsyncException { 234 this.minimizer = m; 235 if (isQuick) { 236 name = "MMFF2D"; 237 ffParams = mmff2DParams; 238 if (ffParams == null) 239 mmff2DParams = ffParams = getParameters(true); 240 } else { 241 name = "MMFF"; 242 ffParams = mmffParams; 243 if (ffParams == null) 244 mmffParams = ffParams = getParameters(false); 245 } 246 } 247 248 @Override clear()249 public void clear() { 250 // not same atoms? 251 // TODO 252 253 } 254 255 @Override setModel(BS bsElements, int elemnoMax)256 public boolean setModel(BS bsElements, int elemnoMax) { 257 Minimizer m = minimizer; 258 if (!setArrays(m.atoms, m.bsAtoms, m.bonds, m.rawBondCount, false, false)) 259 return false; 260 setModelFields(); 261 if (!fixTypes()) 262 return false; 263 calc = new CalculationsMMFF(this, ffParams, minAtoms, minBonds, 264 minAngles, minTorsions, minPositions, minimizer.constraints); 265 calc.setLoggingEnabled(true); 266 return calc.setupCalculations(); 267 } 268 setArrays(Atom[] atoms, BS bsAtoms, Bond[] bonds, int rawBondCount, boolean doRound, boolean allowUnknowns)269 public boolean setArrays(Atom[] atoms, BS bsAtoms, Bond[] bonds, 270 int rawBondCount, boolean doRound, boolean allowUnknowns) { 271 Minimizer m = minimizer; 272 // these are original atom-index-based, not minAtom-index based. 273 274 vRings = AU.createArrayOfArrayList(4); 275 rawAtomTypes = setAtomTypes(atoms, bsAtoms, m.vwr.getSmilesMatcher(), 276 vRings, allowUnknowns); 277 if (rawAtomTypes == null) 278 return false; 279 rawBondTypes = setBondTypes(bonds, rawBondCount, bsAtoms); 280 rawMMFF94Charges = calculatePartialCharges(bonds, rawBondTypes, atoms, 281 rawAtomTypes, bsAtoms, doRound); 282 return true; 283 } 284 private final static String names = "END.BCI.CHG.ANG.NDK.OND.OOP.TBN.FSB.TOR.VDW."; 285 private final static int[] types = {0, TYPE_PBCI, TYPE_CHRG, TYPE_ANGLE, TYPE_BNDK, TYPE_BOND, TYPE_OOP, TYPE_SB, TYPE_SBDEF, TYPE_TORSION, TYPE_VDW }; 286 getParameters(boolean isQuick)287 protected Map<Object, Object> getParameters(boolean isQuick) throws JmolAsyncException { 288 getAtomTypes(); 289 String resourceName = (isQuick ? "mmff94.par.txt" : "mmff94_2d.par.txt"); 290 291 Hashtable<Object, Object> data = new Hashtable<Object, Object>(); 292 if (Logger.debugging) 293 Logger.debug("reading data from " + resourceName); 294 BufferedReader br = null; 295 String line = null; 296 try { 297 br = getBufferedReader(resourceName); 298 int pt = 0; 299 int dataType = 0; 300 while (true) { 301 while ((pt = (line = br.readLine()).indexOf(".PAR")) < 0) {} 302 if ((dataType = types[names.indexOf(line.substring(pt - 3, pt + 1)) / 4]) < 1) 303 break; 304 readParams(br, dataType, data); 305 } 306 br.close(); 307 } catch (JmolAsyncException e) { 308 throw new JmolAsyncException(e.getFileName()); 309 } catch (Exception e) { 310 System.err.println("Exception " + e.toString() + " in getResource " 311 + resourceName + " line=" + line); 312 } finally { 313 try { 314 br.close(); 315 } catch (Exception e) { 316 //ignore 317 } 318 } 319 return data; 320 } 321 322 private String line; 323 readParams(BufferedReader br, int dataType, Map<Object, Object> data)324 private void readParams(BufferedReader br, int dataType, Map<Object, Object> data) throws Exception { 325 // parameters are keyed by a 32-bit Integer 326 // that is composed of four 7-bit atom types and one 4-bit parameter type 327 // in some cases, the last 7-bit atom type (a4) is used for additional parameter typing 328 Object value = null; 329 int a1 = 0, a2 = 127, a3 = 127, a4 = 127; 330 int type = 0; 331 switch (dataType) { 332 case TYPE_BOND: // bond 333 case TYPE_ANGLE: // angle 334 case TYPE_TORSION: // tor 335 break; 336 case TYPE_CHRG: // chrg/bci, identified by a4 = 4 337 a4 = A4_CHRG; 338 break; 339 case TYPE_SB: // stretch bend, identified by a4 340 a4 = A4_SB; 341 break; 342 case TYPE_BNDK: // A4_BNDK identified by a4 343 a4 = A4_BNDK; 344 type = KEY_BNDK; 345 break; 346 case TYPE_OOP: // oop (tor max type is 5) 347 type = KEY_OOP; 348 break; 349 case TYPE_PBCI: // pbci 350 type = KEY_PBCI; 351 break; 352 case TYPE_SBDEF: // default stretch bend, by row; identified by a4 353 a4 = A4_SBDEF; 354 type = KEY_SBDEF; 355 break; 356 case TYPE_VDW: // vdw identified by a4 = 3, not 127 357 a4 = A4_VDW; 358 type = KEY_VDW; 359 break; 360 } 361 while (!br.readLine().startsWith("*")){} 362 while ((line = br.readLine()).startsWith("*")){} 363 do { 364 switch (dataType) { 365 case TYPE_BNDK: 366 case TYPE_OOP: 367 case TYPE_PBCI: 368 case TYPE_SBDEF: 369 break; 370 case TYPE_VDW: 371 if (line.charAt(5) != ' ') 372 continue; // header stuff 373 break; 374 case TYPE_CHRG: 375 if (line.charAt(0) == '4') 376 continue; // I have no idea what type=4 here would mean. It's supposed to be a bond type 377 //$FALL-THROUGH$ 378 case TYPE_ANGLE: 379 case TYPE_BOND: 380 case TYPE_SB: 381 case TYPE_TORSION: 382 type = line.charAt(0) - '0'; 383 break; 384 } 385 switch (dataType) { 386 case TYPE_OOP: 387 case TYPE_TORSION: 388 a4 = ival(18,20); 389 //$FALL-THROUGH$ 390 case TYPE_ANGLE: 391 case TYPE_SB: 392 case TYPE_SBDEF: 393 a3 = ival(13,15); 394 //$FALL-THROUGH$ 395 case TYPE_BNDK: 396 case TYPE_BOND: 397 case TYPE_CHRG: 398 a2 = ival(8,10); 399 //$FALL-THROUGH$ 400 case TYPE_PBCI: 401 case TYPE_VDW: 402 a1 = ival(3,5); 403 break; 404 } 405 switch (dataType) { 406 case TYPE_BNDK: // empirical bond stretch: kb, r0 (reversed in file) 407 value = new double[] { 408 dval(19,25), 409 dval(13,18) }; 410 break; 411 case TYPE_BOND: // bond stretch: kb, r0 412 value = new double[] { 413 dval(14,20), 414 dval(25,31) }; 415 break; 416 case TYPE_ANGLE: // angles: ka, theta0 417 case TYPE_SB: // stretch-bend: kbaIJK, kbaKJI 418 value = new double[] { 419 dval(19,25), 420 dval(28,35) }; 421 break; 422 case TYPE_CHRG: // bond chrg 423 value = Float.valueOf(fval(10,20)); 424 break; 425 case TYPE_OOP: // oop: koop 426 value = new double[] { dval(24,30) }; 427 break; 428 case TYPE_PBCI: 429 value = Float.valueOf(fval(5,15)); 430 break; 431 case TYPE_SBDEF: // default stretch-bend: F(I_J,K),F(K_J,I) 432 double v1 = dval(19,25); 433 double v2 = dval(28,35); 434 value = new double[] { v1, v2 }; 435 Integer key = MinObject.getKey(type, a1, a2, a3, a4); 436 data.put(key, value); 437 value = new double[] { v2, v1 }; 438 int a = a1; 439 a1 = a3; 440 a3 = a; 441 break; 442 case TYPE_TORSION: // tor: v1, v2, v3 443 value = new double[] { 444 dval(22,28), 445 dval(30,36), 446 dval(38,44) 447 }; 448 break; 449 case TYPE_VDW: // vdw alpha-i, N-i, A-i, G-i, DA 450 value = new double[] { 451 dval(10,15), 452 dval(20,25), 453 dval(30,35), 454 dval(40,45), 455 line.charAt(46) // '-', 'A', 'D' 456 }; 457 break; 458 } 459 Integer key = MinObject.getKey(type, a1, a2, a3, a4); 460 data.put(key, value); 461 if (Logger.debugging) 462 Logger.debug(MinObject.decodeKey(key) + " " + (value instanceof Float ? value : Escape.eAD((double[])value))); 463 } while (!(line = br.readLine()).startsWith("$")); 464 } 465 ival(int i, int j)466 private int ival(int i, int j) { 467 return PT.parseInt(line.substring(i,j).trim()); 468 } 469 fval(int i, int j)470 private float fval(int i, int j) { 471 return Float.valueOf(line.substring(i,j).trim()).floatValue(); 472 } 473 dval(int i, int j)474 private double dval(int i, int j) { 475 return Double.valueOf(line.substring(i,j).trim()).doubleValue(); 476 } 477 getAtomTypes()478 private void getAtomTypes() throws JmolAsyncException { 479 String resourceName = "MMFF94-smarts.txt"; 480 Lst<AtomType> types = new Lst<AtomType>(); 481 try { 482 BufferedReader br = getBufferedReader(resourceName); 483 //turns out from the Jar file 484 // it's a sun.net.www.protocol.jar.JarURLConnection$JarURLInputStream 485 // and within Eclipse it's a BufferedInputStream 486 487 AtomType at; 488 types.addLast(new AtomType(0, 0, 0, 0, 1, "H or NOT FOUND", "")); 489 while ((line = br.readLine()) != null) { 490 if (line.length() == 0 || line.startsWith("#")) 491 continue; 492 //0 1 2 3 4 5 6 493 //0123456789012345678901234567890123456789012345678901234567890123456789 494 //Mg 12 99 0 24 0 DIPOSITIVE MAGNESIUM CATI [MgD0] 495 //#AtSym ElemNo mmType HType formalCharge*12 val Desc Smiles 496 int elemNo = ival(3,5); 497 int mmType = ival(6,8); 498 int hType = ival(9,11); 499 float formalCharge = fval(12,15)/12; 500 int val = ival(16,18); 501 String desc = line.substring(19,44).trim(); 502 String smarts = line.substring(45).trim(); 503 types.addLast(at = new AtomType(elemNo, mmType, hType, formalCharge, val, desc, smarts)); 504 setFlags(at); 505 } 506 br.close(); 507 } catch (JmolAsyncException e) { 508 throw new JmolAsyncException(e.getFileName()); 509 } catch (Exception e) { 510 System.err.println("Exception " + e.toString() + " in getResource " 511 + resourceName + " line=" + line); 512 513 } 514 Logger.info((types.size()-1) + " SMARTS-based atom types read"); 515 atomTypes = types; 516 517 } 518 setFlags(AtomType at)519 private static void setFlags(AtomType at) { 520 // fcadj 521 switch (at.mmType) { 522 523 // Note that these are NOT fractional charges based on 524 // number of connected atoms. These are relatively arbitrary 525 // fractions of the formal charge to be shared with other atoms. 526 // That is, it is not significant that 0.5 is 1/2, and 0.25 is 1/4; 527 // they are just numbers. 528 529 case 32: 530 case 35: 531 case 72: 532 // 32 OXYGEN IN CARBOXYLATE ANION 533 // 32 NITRATE ANION OXYGEN 534 // 32 SINGLE TERMINAL OXYGEN ON TETRACOORD SULFUR 535 // 32 TERMINAL O-S IN SULFONES AND SULFONAMIDES 536 // 32 TERMINAL O IN SULFONATES 537 // 35 OXIDE OXYGEN ON SP2 CARBON, NEGATIVELY CHARGED 538 // 72 TERMINAL SULFUR BONDED TO PHOSPHORUS 539 at.fcadj = 0.5f; 540 break; 541 case 62: 542 case 76: 543 // 62 DEPROTONATED SULFONAMIDE N-; FORMAL CHARGE=-1 544 // 76 NEGATIVELY CHARGED N IN, E.G, TRI- OR TETRAZOLE ANION 545 at.fcadj = 0.25f; 546 break; 547 } 548 549 // arom 550 switch (at.mmType) { 551 case 37: 552 case 38: 553 case 39: 554 case 44: 555 case 58: 556 case 59: 557 case 63: 558 case 64: 559 case 65: 560 case 66: 561 case 69: 562 case 78: 563 case 79: 564 case 81: 565 case 82: 566 at.arom = true; 567 } 568 569 // sbmb 570 switch (at.mmType) { 571 case 2: 572 case 3: 573 case 4: 574 case 9: 575 case 30: 576 case 37: 577 case 39: 578 case 54: 579 case 57: 580 case 58: 581 case 63: 582 case 64: 583 case 67: 584 case 75: 585 case 78: 586 case 80: 587 case 81: 588 at.sbmb = true; 589 } 590 591 // pilp 592 switch(at.mmType) { 593 case 6: 594 case 8: 595 case 10: 596 case 11: 597 case 12: 598 case 13: 599 case 14: 600 case 15: 601 case 26: 602 case 32: 603 case 35: 604 case 39: 605 case 40: 606 case 43: 607 case 44: 608 case 59: 609 case 62: 610 case 70: 611 case 72: 612 case 76: 613 at.pilp = true; 614 } 615 616 // mltb: 617 switch (at.mmType) { 618 case 10: 619 case 32: 620 case 35: 621 case 39: 622 case 41: 623 case 44: 624 case 55: 625 case 56: 626 case 58: 627 case 59: 628 case 69: 629 case 72: 630 case 81: 631 case 82: 632 at.mltb = 1; 633 break; 634 case 2: 635 case 3: 636 case 7: 637 case 9: 638 case 16: 639 case 17: 640 case 30: 641 case 37: 642 case 38: 643 case 45: 644 case 46: 645 case 47: 646 case 51: 647 case 53: 648 case 54: 649 case 57: 650 case 63: 651 case 64: 652 case 65: 653 case 66: 654 case 67: 655 case 74: 656 case 75: 657 case 78: 658 case 79: 659 case 80: 660 at.mltb = 2; 661 break; 662 case 4: 663 case 42: 664 case 60: 665 case 61: 666 at.mltb = 3; 667 break; 668 } 669 } 670 /** 671 * assign partial charges ala MMFF94 672 * 673 * @param bonds 674 * @param bTypes 675 * @param atoms 676 * @param aTypes 677 * @param bsAtoms 678 * @param doRound 679 * @return full array of partial charges 680 */ calculatePartialCharges(Bond[] bonds, int[] bTypes, Atom[] atoms, int[] aTypes, BS bsAtoms, boolean doRound)681 public float[] calculatePartialCharges(Bond[] bonds, int[] bTypes, 682 Atom[] atoms, int[] aTypes, BS bsAtoms, 683 boolean doRound) { 684 685 // start with formal charges specified by MMFF94 (not what is in file!) 686 687 float[] partialCharges = new float[atoms.length]; 688 for (int i = bsAtoms.nextSetBit(0); i >= 0; i = bsAtoms.nextSetBit(i + 1)) 689 partialCharges[i] = atomTypes.get(Math.max(0, aTypes[i])).formalCharge; 690 // run through all bonds, adjusting formal charges as necessary 691 Atom a1 = null; 692 for (int i = bTypes.length; --i >= 0;) { 693 a1 = bonds[i].atom1; 694 Atom a2 = bonds[i].atom2; 695 // It's possible that some of our atoms are not in the atom set, 696 // but we don't want both of them to be out of the set. 697 698 boolean ok1 = bsAtoms.get(a1.i); 699 boolean ok2 = bsAtoms.get(a2.i); 700 if (!ok1 && !ok2) 701 continue; 702 int it = aTypes[a1.i]; 703 AtomType at1 = atomTypes.get(Math.max(0, it)); 704 int type1 = (it < 0 ? -it : at1.mmType); 705 it = aTypes[a2.i]; 706 AtomType at2 = atomTypes.get(Math.max(0, it)); 707 int type2 = (it < 0 ? -it : at2.mmType); 708 709 // we are only interested in bonds that are between different atom types 710 711 // if (type1 == type2) 712 // continue; 713 714 // check for bond charge increment 715 716 // The table is created using the key (100 * type1 + type2), 717 // where type1 < type2. In addition, we encode the partial bci values 718 // with key (100 * type) 719 720 float dq = Float.NaN; // the difference in charge to be added or subtracted from the formal charges 721 try { 722 int bondType = bTypes[i]; 723 float bFactor = (type1 < type2 ? -1 : 1); 724 Integer key = MinObject.getKey(bondType, bFactor == 1 ? type2 : type1, 725 bFactor == 1 ? type1 : type2, 127, A4_CHRG); 726 Float bciValue = (Float) ffParams.get(key); 727 float bci = Float.NaN; 728 String msg = (Logger.debugging 729 ? a1 + "/" + a2 + " mmTypes=" + type1 + "/" + type2 730 + " formalCharges=" + at1.formalCharge + "/" + at2.formalCharge 731 + " bci = " 732 : null); 733 if (bciValue == null) { 734 // no bci was found; we have to use partial bond charge increments 735 // a failure here indicates we don't have information 736 Float a, b; 737 if ((a = (Float) ffParams 738 .get(MinObject.getKey(KEY_PBCI, type1, 127, 127, 127))) != null 739 && (b = (Float) ffParams.get( 740 MinObject.getKey(KEY_PBCI, type2, 127, 127, 127))) != null) { 741 float pa = a.floatValue(); 742 float pb = b.floatValue(); 743 bci = pa - pb; 744 if (Logger.debugging) 745 msg += pa + " - " + pb + " = "; 746 } 747 } else { 748 bci = bFactor * bciValue.floatValue(); 749 } 750 if (Logger.debugging) { 751 msg += bci; 752 Logger.debug(msg); 753 } 754 // Here's the way to do this: 755 // 756 // 1) The formal charge on each atom is adjusted both by 757 // taking on an (arbitrary?) fraction of the formal charge on its partner 758 // and by giving up an (arbitrary?) fraction of its own formal charge 759 // Note that the formal charge is the one specified in the MMFF94 parameters, 760 // NOT the model file. The compounds in MMFF94_dative.mol2, for example, do 761 // not indicate formal charges. The only used fractions are 0, 0.25, and 0.5. 762 // 763 // 2) Then the bond charge increment is added. 764 // 765 // Note that this value I call "dq" is added to one atom and subtracted from its partner 766 767 dq = at2.fcadj * at2.formalCharge - at1.fcadj * at1.formalCharge + bci; 768 } catch (Exception e) { 769 dq = Float.NaN; 770 } 771 if (ok1) 772 partialCharges[a1.i] += dq; 773 if (ok2) 774 partialCharges[a2.i] -= dq; 775 } 776 777 // just rounding to 0.001 here: 778 779 if (doRound) { 780 float abscharge = 0; 781 for (int i = partialCharges.length; --i >= 0;) { 782 partialCharges[i] = (Math.round(partialCharges[i] * 1000)) / 1000f; 783 abscharge += Math.abs(partialCharges[i]); 784 } 785 if (abscharge == 0 && a1 != null) { 786 partialCharges[a1.i]= -0.0f; 787 } 788 } 789 return partialCharges; 790 } 791 792 /** 793 * From forcefieldmmff94.cpp (flag BTij) 794 * 795 * a) single bond between atoms i and j, both i and j are not aromatic and 796 * both types have sbmb set in mmffprop.par, or 797 * 798 * b) between two aromatic atoms, but the bond is not aromatic (e.g. 799 * connecting bond in biphenyl) 800 * 801 * (sbmb is 2, 3, 4, 9, 30, 37, 39, 54, 57, 58, 63, 64, 67, 75, 78, 80, 81) 802 * 803 * 804 * @param at1 805 * @param at2 806 * @return 0 or 1 807 */ isSpecialBondType(AtomType at1, AtomType at2)808 private static boolean isSpecialBondType(AtomType at1, AtomType at2) { 809 return at1.sbmb && at2.sbmb || at1.arom && at2.arom; 810 // but what about at1.sbmb && at2.arom? 811 } 812 813 /** 814 * Get the bond type: 815 * 816 * 1 biphenyl or 817 * 818 * 0 any other 819 * 820 * @param bond 821 * @param at1 822 * @param at2 823 * @param index1 824 * @param index2 825 * @return 826 */ getBondType(Bond bond, AtomType at1, AtomType at2, int index1, int index2)827 private int getBondType(Bond bond, AtomType at1, AtomType at2, 828 int index1, int index2) { 829 return (isSpecialBondType(at1, at2) && 830 bond.getCovalentOrder() == 1 831 && !isAromaticBond(index1, index2) ? 1 : 0); 832 } 833 isAromaticBond(int a1, int a2)834 private boolean isAromaticBond(int a1, int a2) { 835 if (vRings[Raromatic] != null) 836 for (int i = vRings[Raromatic].size(); --i >= 0;) { 837 BS bsRing = vRings[Raromatic].get(i); 838 if (bsRing.get(a1) && bsRing.get(a2)) 839 return true; 840 } 841 return false; 842 } 843 getAtomTypeDescs(int[] types)844 public static String[] getAtomTypeDescs(int[] types) { 845 String[] stypes = new String[types.length]; 846 for (int i = types.length; --i >= 0;) { 847 stypes[i] = String.valueOf(types[i] < 0 ? -types[i] : atomTypes.get(types[i]).mmType); 848 } 849 return stypes; 850 } 851 852 ///////////// MMFF94 object typing ////////////// 853 854 /** 855 * The file MMFF94-smarts.txt is derived from MMFF94-smarts.xlsx. This file 856 * contains records for unique atom type/formal charge sharing/H atom type. 857 * For example, the MMFF94 type 6 is distributed over eight AtomTypes, each 858 * with a different SMARTS match. 859 * 860 * H atom types are given in the file as properties of other atom types, not 861 * as their own individual SMARTS searches. H atom types are determined based 862 * on their attached atom's atom type. 863 * 864 * @param atoms 865 * @param bsAtoms 866 * @param smartsMatcher 867 * @param vRings 868 * @param allowUnknowns 869 * @return array of indexes into AtomTypes or, for H, negative of mmType 870 */ setAtomTypes(Atom[] atoms, BS bsAtoms, SmilesMatcherInterface smartsMatcher, Lst<BS>[] vRings, boolean allowUnknowns)871 private static int[] setAtomTypes(Atom[] atoms, BS bsAtoms, 872 SmilesMatcherInterface smartsMatcher, 873 Lst<BS>[] vRings, boolean allowUnknowns) { 874 Lst<BS> bitSets = new Lst<BS>(); 875 String[] smarts = new String[atomTypes.size()]; 876 int[] types = new int[atoms.length]; 877 BS bsElements = new BS(); 878 BS bsHydrogen = new BS(); 879 BS bsConnected = BSUtil.copy(bsAtoms); 880 881 // It may be important to include all attached atoms 882 883 for (int i = bsAtoms.nextSetBit(0); i >= 0; i = bsAtoms.nextSetBit(i + 1)) { 884 Atom a = atoms[i]; 885 Bond[] bonds = a.bonds; 886 if (bonds != null) 887 for (int j = bonds.length; --j >= 0;) 888 if (bonds[j].isCovalent()) 889 bsConnected.set(bonds[j].getOtherAtom(a).i); 890 } 891 892 // we need to identify H atoms and also make a BitSet of all the elements 893 894 for (int i = bsConnected.nextSetBit(0); i >= 0; i = bsConnected 895 .nextSetBit(i + 1)) { 896 int n = atoms[i].getElementNumber(); 897 switch (n) { 898 case 1: 899 bsHydrogen.set(i); 900 break; 901 default: 902 bsElements.set(n); 903 } 904 } 905 906 // generate a list of SMART codes 907 908 int nUsed = 0; 909 for (int i = 1; i < atomTypes.size(); i++) { 910 AtomType at = atomTypes.get(i); 911 if (!bsElements.get(at.elemNo)) 912 continue; 913 smarts[i] = at.smartsCode; 914 nUsed++; 915 } 916 Logger.info(nUsed + " SMARTS matches used"); 917 918 // The SMARTS list is organized from least general to most general 919 // for each atom. So the FIRST occurrence of an atom in the list 920 // identifies that atom's MMFF94 type. 921 922 try { 923 smartsMatcher.getMMFF94AtomTypes(smarts, atoms, atoms.length, 924 bsConnected, bitSets, vRings); 925 } catch (Exception e) { 926 Logger.error(e.toString()); 927 } 928 BS bsDone = new BS(); 929 for (int j = 0; j < bitSets.size(); j++) { 930 BS bs = bitSets.get(j); 931 if (bs == null) 932 continue; 933 // This is a one-pass system. We first exclude 934 // all atoms that are already identified... 935 bs.andNot(bsDone); 936 // then we set the type of what is remaining... 937 for (int i = bs.nextSetBit(0); i >= 0; i = bs.nextSetBit(i + 1)) 938 types[i] = j; 939 // then we include these atoms in the set of atoms already identified 940 bsDone.or(bs); 941 } 942 943 // now we add in the H atom types as the negative of their MMFF94 type 944 // rather than as an index into AtomTypes. 945 946 for (int i = bsHydrogen.nextSetBit(0); i >= 0; i = bsHydrogen 947 .nextSetBit(i + 1)) { 948 Bond[] bonds = atoms[i].bonds; 949 if (bonds != null) { 950 int j = types[bonds[0].getOtherAtom(atoms[i]).i]; 951 if (j != 0) 952 bsDone.set(i); 953 types[i] = -atomTypes.get(j).hType; 954 } 955 } 956 if (Logger.debugging) 957 for (int i = bsConnected.nextSetBit(0); i >= 0; i = bsConnected 958 .nextSetBit(i + 1)) 959 Logger.debug("atom " 960 + atoms[i] 961 + "\ttype " 962 + (types[i] < 0 ? "" + -types[i] : (atomTypes.get(types[i]).mmType 963 + "\t" + atomTypes.get(types[i]).smartsCode + "\t" + atomTypes 964 .get(types[i]).descr))); 965 966 if (!allowUnknowns && bsDone.cardinality() != bsConnected.cardinality()) 967 return null; 968 return types; 969 } 970 setBondTypes(Bond[] bonds, int bondCount, BS bsAtoms)971 private int[] setBondTypes(Bond[] bonds, int bondCount, BS bsAtoms) { 972 int[] bTypes = new int[bondCount]; 973 for (int i = bondCount; --i >= 0;) { 974 Atom a1 = bonds[i].atom1; 975 Atom a2 = bonds[i].atom2; 976 boolean ok1 = bsAtoms.get(a1.i); 977 boolean ok2 = bsAtoms.get(a2.i); 978 if (!ok1 && !ok2) 979 continue; 980 int it = rawAtomTypes[a1.i]; 981 AtomType at1 = atomTypes.get(Math.max(0, it)); 982 it = rawAtomTypes[a2.i]; 983 AtomType at2 = atomTypes.get(Math.max(0, it)); 984 bTypes[i] = getBondType(bonds[i], at1, at2, a1.i, a2.i); 985 } 986 return bTypes; 987 } 988 fixTypes()989 private boolean fixTypes() { 990 // set atom types in minAtoms 991 for (int i = minAtomCount; --i >= 0;) { 992 MinAtom a = minAtoms[i]; 993 int rawIndex = a.atom.i; 994 int it = rawAtomTypes[rawIndex]; 995 a.ffAtomType = atomTypes.get(Math.max(0, it)); 996 int type = (it < 0 ? -it : atomTypes.get(it).mmType); 997 a.ffType = type; 998 a.vdwKey = MinObject.getKey(KEY_VDW, type, 127, 127, A4_VDW); 999 a.partialCharge = rawMMFF94Charges[rawIndex]; 1000 } 1001 1002 for (int i = minBonds.length; --i >= 0;) { 1003 MinBond bond = minBonds[i]; 1004 bond.type = rawBondTypes[bond.rawIndex]; 1005 bond.key = getKey(bond, bond.type, TYPE_BOND); 1006 if (bond.key == null) 1007 return false; 1008 } 1009 1010 for (int i = minAngles.length; --i >= 0;) { 1011 MinAngle angle = minAngles[i]; 1012 angle.key = getKey(angle, angle.type, TYPE_ANGLE); 1013 angle.sbKey = getKey(angle, angle.sbType, TYPE_SB); 1014 } 1015 1016 for (int i = minTorsions.length; --i >= 0;) { 1017 MinTorsion torsion = minTorsions[i]; 1018 torsion.key = getKey(torsion, torsion.type, TYPE_TORSION); 1019 } 1020 return true; 1021 } 1022 1023 private final static int[] sbMap = {0, 1, 3, 5, 4, 6, 8, 9, 11}; 1024 1025 /** 1026 * Get the angle type: 1027 * 1028 * 0 The angle <i>i-j-k</i> is a "normal" bond angle 1029 * 1030 * 1 Either bond <i>i-j</i> or bond <i>j-k</i> has a bond type of 1 1031 * 1032 * 2 Bonds<i> i-j</i> and <i>j-k</i> each have bond types of 1; the sum is 2. 1033 * 1034 * 3 The angle occurs in a three-membered ring 1035 * 1036 * 4 The angle occurs in a four-membered ring 1037 * 1038 * 5 Is in a three-membered ring and the sum of the bond types is 1 1039 * 1040 * 6 Is in a three-membered ring and the sum of the bond types is 2 1041 * 1042 * 7 Is in a four-membered ring and the sum of the bond types is 1 1043 * 1044 * 8 Is in a four-membered ring and the sum of the bond types is 2 1045 */ setAngleType(MinAngle angle)1046 private int setAngleType(MinAngle angle) { 1047 angle.type = minBonds[angle.data[ABI_IJ]].type + minBonds[angle.data[ABI_JK]].type; 1048 if (checkRings(vRings[R3], angle.data, 3)) { 1049 angle.type += (angle.type == 0 ? 3 : 4); 1050 } else if (checkRings(vRings[R4], angle.data, 3)) { 1051 angle.type += (angle.type == 0 ? 4 : 6); 1052 } 1053 1054 /* 1055 SBT AT BT[IJ] BT[JK] 1056 ------------------------------------------------------------- 1057 0 0 0 0 1058 1 1 1 0 1059 2 1 0 1 [error in table] 1060 3 2 1 1 [error in table] 1061 4 4 0 0 1062 5 3 0 0 1063 6 5 1 0 1064 7 5 0 1 1065 8 6 1 1 1066 9 7 1 0 1067 10 7 0 1 1068 11 8 1 1 1069 */ 1070 1071 angle.sbType = sbMap[angle.type]; 1072 switch (angle.type) { 1073 case 1: 1074 case 5: 1075 case 7: 1076 angle.sbType += minBonds[angle.data[ABI_JK]].type; 1077 break; 1078 } 1079 return angle.type; 1080 } 1081 1082 /** 1083 * Get the torsion type for [a,b,c,d], also called "FF class". One of the 1084 * following, determined in this order: 1085 * 1086 * 4: 4-membered ring 1087 * 1088 * 1: ab-cd 1089 * 1090 * 0: biphenyl and not 5-membered ring 1091 * 1092 * 5: 5-membered ring 1093 * 1094 * 2: a-x=x-b 1095 * 1096 * @param t 1097 * @return 1098 */ setTorsionType(MinTorsion t)1099 private int setTorsionType(MinTorsion t) { 1100 if (checkRings(vRings[R4], t.data, 4)) 1101 return (t.type = 4); // in 4-membered ring 1102 t.type = (minBonds[t.data[TBI_BC]].type == 1 ? 1 1103 : minBonds[t.data[TBI_AB]].type == 0 && minBonds[t.data[TBI_CD]].type == 0 ? 0 : 2); 1104 if (t.type == 0 && checkRings(vRings[R5], t.data, 4)) { 1105 t.type = 5; // in 5-membered ring 1106 } 1107 return t.type; 1108 } 1109 typeOf(int iAtom)1110 private int typeOf(int iAtom) { 1111 return minAtoms[iAtom].ffType; 1112 } 1113 checkRings(Lst<BS> v, int[] minlist, int n)1114 private boolean checkRings(Lst<BS> v, int[] minlist, int n) { 1115 if (v != null) 1116 for (int i = v.size(); --i >= 0;) { 1117 BS bs = v.get(i); 1118 if (bs.get(minAtoms[minlist[0]].atom.i) 1119 && bs.get(minAtoms[minlist[1]].atom.i) 1120 && (n < 3 || bs.get(minAtoms[minlist[2]].atom.i)) 1121 && (n < 4 || bs.get(minAtoms[minlist[3]].atom.i))) 1122 return true; 1123 } 1124 return false; 1125 } 1126 getKey(Object obj, int type, int ktype)1127 private Integer getKey(Object obj, int type, int ktype) { 1128 MinObject o = (obj instanceof MinObject ? (MinObject) obj : null); 1129 int[] data = (o == null ? (int[]) obj : o.data); 1130 int n = 4; 1131 switch (ktype) { 1132 case TYPE_BOND: 1133 fixOrder(data, 0, 1); 1134 n = 2; 1135 break; 1136 case TYPE_ANGLE: 1137 if (fixOrder(data, 0, 2) == -1) 1138 swap(data, ABI_IJ, ABI_JK); 1139 type = setAngleType((MinAngle) o); 1140 n = 3; 1141 break; 1142 case TYPE_SB: 1143 n = 3; 1144 break; 1145 case TYPE_TORSION: 1146 switch (fixOrder(data, 1, 2)) { 1147 case 1: 1148 break; 1149 case -1: 1150 swap(data, 0, 3); 1151 swap(data, TBI_AB, TBI_CD); 1152 break; 1153 case 0: 1154 if (fixOrder(data, 0, 3) == -1) 1155 swap(data, TBI_AB, TBI_CD); 1156 break; 1157 } 1158 type = setTorsionType((MinTorsion) o); 1159 } 1160 Integer key = null; 1161 for (int i = 0; i < 4; i++) 1162 typeData[i] = (i < n ? typeOf(data[i]) : 127); 1163 switch (ktype) { 1164 case TYPE_SB: 1165 typeData[3] = A4_SB; 1166 break; 1167 case TYPE_OOP: 1168 sortOop(typeData); 1169 break; 1170 } 1171 key = MinObject.getKey(type, typeData[0], typeData[1], typeData[2], 1172 typeData[3]); 1173 double[] ddata = (double[]) ffParams.get(key); 1174 // default typing 1175 switch (ktype) { 1176 case TYPE_BOND: 1177 // if (!useEmpiricalRules) 1178 // return key; 1179 return (ddata != null && ddata[0] > 0 ? key : applyEmpiricalRules(o, ddata, TYPE_BOND)); 1180 case TYPE_ANGLE: 1181 if (ddata != null && ddata[0] != 0) 1182 return key; 1183 break; 1184 case TYPE_TORSION: 1185 if (ddata == null) { 1186 if (!ffParams.containsKey(key = getTorsionKey(type, 0, 2)) 1187 && !ffParams.containsKey(key = getTorsionKey(type, 2, 0)) 1188 && !ffParams.containsKey(key = getTorsionKey(type, 2, 2))) 1189 key = getTorsionKey(0, 2, 2); 1190 ddata = (double[]) ffParams.get(key); 1191 } 1192 // if (!useEmpiricalRules) 1193 // return key; 1194 return (ddata != null ? key : applyEmpiricalRules(o, ddata, TYPE_TORSION)); 1195 case TYPE_SB: 1196 // use periodic row info 1197 if (ddata != null) 1198 return key; 1199 int r1 = getRowFor(data[0]); 1200 int r2 = getRowFor(data[1]); 1201 int r3 = getRowFor(data[2]); 1202 return MinObject.getKey(KEY_SBDEF, r1, r2, r3, A4_SBDEF); 1203 case TYPE_OOP: 1204 // use periodic row info 1205 if (ddata != null) 1206 return key; 1207 } 1208 // run through equivalent types, really just 3 1209 // if (!useEmpiricalRules && ddata != null) 1210 // return key; 1211 boolean isSwapped = false; 1212 boolean haveKey = false; 1213 for (int i = 0; i < 3 && !haveKey; i++) { 1214 for (int j = 0, bit = 1; j < n; j++, bit <<= 1) 1215 if ((ktype & bit) == bit) 1216 typeData[j] = getEquivalentType(typeOf(data[j]), i); 1217 switch (ktype) { 1218 case TYPE_BOND: 1219 // not really supposed to do this for MMFF94 1220 isSwapped = (fixTypeOrder(typeData, 0, 1)); 1221 break; 1222 case TYPE_ANGLE: 1223 isSwapped = (fixTypeOrder(typeData, 0, 2)); 1224 break; 1225 case TYPE_OOP: 1226 sortOop(typeData); 1227 break; 1228 } 1229 key = MinObject.getKey(type, typeData[0], typeData[1], typeData[2], 1230 typeData[3]); 1231 haveKey = ffParams.containsKey(key); 1232 } 1233 if (haveKey) { 1234 if (isSwapped) 1235 switch (ktype) { 1236 case TYPE_ANGLE: 1237 swap(data, 0, 2); 1238 swap(data, ABI_IJ, ABI_JK); 1239 setAngleType((MinAngle) o); 1240 break; 1241 } 1242 } else if (type != 0 && ktype == TYPE_ANGLE) { 1243 key = Integer.valueOf(key.intValue() ^ 0xFF); 1244 } 1245 1246 // if (!useEmpiricalRules) 1247 // return key; 1248 // 1249 ddata = (double[]) ffParams.get(key); 1250 // default typing 1251 switch (ktype) { 1252 case TYPE_ANGLE: 1253 return (ddata != null && ddata[0] != 0 ? key : 1254 applyEmpiricalRules(o, ddata, TYPE_ANGLE)); 1255 } 1256 return key; 1257 } 1258 getTorsionKey(int type, int i, int j)1259 private Integer getTorsionKey(int type, int i, int j) { 1260 return MinObject.getKey(type, getEquivalentType(typeData[0], i), 1261 typeData[1], typeData[2], getEquivalentType(typeData[3], j)); 1262 } 1263 applyEmpiricalRules(MinObject o, double[] ddata, int ktype)1264 private Integer applyEmpiricalRules(MinObject o, double[] ddata, int ktype) { 1265 double rr, rr2, beta = 0; 1266 MinAtom a, b, c; 1267 switch (ktype) { 1268 case TYPE_BOND: 1269 a = minAtoms[o.data[0]]; 1270 b = minAtoms[o.data[1]]; 1271 int elemno1 = a.atom.getElementNumber(); 1272 int elemno2 = b.atom.getElementNumber(); 1273 Integer key = MinObject.getKey(KEY_BNDK, Math.min(elemno1, elemno2), Math 1274 .max(elemno1, elemno2), 127, A4_BNDK); 1275 ddata = (double[]) ffParams.get(key); 1276 if (ddata == null) 1277 return null; 1278 double kbref = ddata[0]; 1279 double r0ref = ddata[1]; 1280 1281 double r0 = getRuleBondLength(a, b, ((MinBond) o).order, isAromaticBond( 1282 o.data[0], o.data[1])); 1283 if (r0 == 0) 1284 return null; 1285 1286 rr = r0ref / r0; 1287 rr2 = rr * rr; 1288 double rr4 = rr2 * rr2; 1289 double rr6 = rr4 * rr2; 1290 double kb = kbref * rr6; 1291 o.ddata = new double[] { kb, r0 }; 1292 return Integer.valueOf(-1); 1293 case TYPE_ANGLE: 1294 // from OpenBabel 1295 double theta0; 1296 if (ddata == null || (theta0 = ddata[1]) == 0) { 1297 b = minAtoms[o.data[1]]; 1298 Atom atom = b.atom; 1299 int elemno = atom.getElementNumber(); 1300 switch (o.type) { 1301 case 3: 1302 case 5: 1303 case 6: 1304 // 3-membered rings 1305 theta0 = 60; 1306 beta *= 0.05; 1307 break; 1308 case 4: 1309 case 7: 1310 case 8: 1311 // 4-membered rings 1312 theta0 = 90; 1313 break; 1314 default: 1315 // everything else 1316 theta0 = 120; 1317 int crd = atom.getCovalentBondCount(); 1318 switch (crd) { 1319 case 2: 1320 if (MinAtom.isLinear(b)) 1321 theta0 = 180; 1322 else if (elemno == 8) 1323 theta0 = 105; 1324 else if (elemno > 10) 1325 theta0 = 95.0; 1326 break; 1327 case 3: 1328 if (b.ffAtomType.mltb == 0 && b.ffAtomType.val == 3) 1329 theta0 = (elemno == 7 ? 107 : 92); 1330 break; 1331 case 4: 1332 theta0 = 109.45; 1333 break; 1334 } 1335 } 1336 } 1337 1338 beta = 1.75; 1339 switch (o.type) { 1340 case 3: 1341 case 5: 1342 case 6: 1343 // 3-membered rings 1344 beta *= 0.05; 1345 break; 1346 case 4: 1347 case 7: 1348 case 8: 1349 // 4-membered rings 1350 beta *= 0.85; 1351 break; 1352 } 1353 double za = getZParam(minAtoms[o.data[0]].atom.getElementNumber()); 1354 double cb = getCParam(minAtoms[o.data[1]].atom.getElementNumber()); 1355 double zc = getZParam(minAtoms[o.data[2]].atom.getElementNumber()); 1356 double r0ab = getR0(minBonds[o.data[ABI_IJ]]); 1357 double r0bc = getR0(minBonds[o.data[ABI_JK]]); 1358 rr = r0ab + r0bc; 1359 rr2 = rr * rr; 1360 double D = (r0ab - r0bc) / rr2; 1361 double theta2 = theta0 * Calculations.DEG_TO_RAD; 1362 theta2 *= theta2; 1363 double ka = (beta * za * cb * zc * Math.exp(-2 * D)) / (rr * theta2); 1364 o.ddata = new double[] { ka, theta0 }; 1365 return Integer.valueOf(-1); 1366 case TYPE_TORSION: 1367 int ib = o.data[1]; 1368 int ic = o.data[2]; 1369 b = minAtoms[ib]; 1370 c = minAtoms[ic]; 1371 1372 // rule (a) page 631: no linear systems 1373 1374 if (MinAtom.isLinear(b) || MinAtom.isLinear(c)) 1375 return null; 1376 1377 MinBond bondBC = minBonds[o.data[TBI_BC]]; 1378 1379 int elemnoB = b.atom.getElementNumber(); 1380 int elemnoC = c.atom.getElementNumber(); 1381 double ub = getUParam(elemnoB); 1382 double uc = getUParam(elemnoC); 1383 double vb = getVParam(elemnoB); 1384 double vc = getVParam(elemnoC); 1385 1386 double v1 = 0, 1387 v2 = 0, 1388 v3 = 0; 1389 double pi_bc = -1; // Eqn 21 1390 double n_bc = -1; // Eqn 22 1391 double wb = -1, 1392 wc = 0; // Eqn 23 1393 int valB = b.ffAtomType.val; 1394 int valC = c.ffAtomType.val; 1395 boolean pilpB = b.ffAtomType.pilp; 1396 boolean pilpC = c.ffAtomType.pilp; 1397 int mltbB = b.ffAtomType.mltb; 1398 int mltbC = c.ffAtomType.mltb; 1399 1400 out: while (true) { 1401 1402 // rule (b) page 631: aromatics 1403 1404 if (isAromaticBond(ib, ic)) { 1405 pi_bc = (pilpB || pilpC ? 0.3 : 0.5); 1406 beta = (valB + valC == 7 ? 3 : 6); 1407 break out; 1408 } 1409 1410 // rule (c) page 631: full double bonds 1411 1412 if (bondBC.order == 2) { 1413 beta = 6; 1414 pi_bc = (mltbB == 2 && mltbC == 2 ? 1.0 : 0.4); 1415 break out; 1416 } 1417 1418 // single bonds only from here on out: 1419 1420 int crdB = b.atom.getCovalentBondCount(); 1421 int crdC = c.atom.getCovalentBondCount(); 1422 1423 // rule (d) page 632: [XD4]-[XD4] 1424 1425 if (crdB == 4 && crdC == 4) { 1426 vb = getVParam(elemnoB); 1427 vc = getVParam(elemnoC); 1428 n_bc = 9; 1429 break out; 1430 } 1431 1432 // rules (e) and (f) page 632, simplified 1433 1434 if (crdB != 4 && (valB > crdB || mltbB > 0) || crdC != 4 1435 && (valC > crdC || mltbC > 0)) 1436 return null; 1437 1438 // rule (g) page 632 (resonant interactions) 1439 1440 boolean case2 = (pilpB && mltbC > 0); 1441 boolean case3 = (pilpC && mltbB > 0); 1442 1443 if (bondBC.order == 1 && (mltbB > 0 && mltbC > 0 || case2 || case3)) { 1444 1445 if (pilpB && pilpC) 1446 return null; 1447 1448 beta = 6; 1449 1450 if (case2) { 1451 pi_bc = (mltbC == 1 ? 0.5 : elemnoB <= 10 && elemnoC <= 10 ? 0.3 1452 : 0.15); 1453 break out; 1454 } 1455 1456 if (case3) { 1457 pi_bc = (mltbB == 1 ? 0.5 : elemnoB <= 10 && elemnoC <= 10 ? 0.3 1458 : 0.15); 1459 break out; 1460 } 1461 1462 if ((mltbB == 1 || mltbC == 1) && (elemnoB == 6 || elemnoC == 6)) { 1463 pi_bc = 0.4; 1464 break out; 1465 } 1466 1467 pi_bc = 0.15; 1468 break out; 1469 } 1470 1471 // rule (h) page 632 (Eqn 23) [XD2-XD2] 1472 1473 switch (elemnoB << 8 + elemnoC) { 1474 case 0x808: // O, O 1475 wb = wc = 2; 1476 break out; 1477 case 0x810: // O, S 1478 wb = 2; 1479 wc = 8; 1480 break out; 1481 case 0x1008: // S, O 1482 wb = 8; 1483 wc = 2; 1484 break out; 1485 case 0x1010: // S, S 1486 wb = wc = 8; 1487 break out; 1488 } 1489 1490 // everything else -- generic Eqn 22 1491 1492 n_bc = crdB * crdC; 1493 break out; 1494 } 1495 if (pi_bc > 0) 1496 v2 = beta * pi_bc * Math.sqrt(ub * uc); // Eqn 21 1497 else if (n_bc > 0) 1498 v3 = Math.sqrt(vb * vc) / n_bc; // Eqn 22 1499 else if (wb != 0) 1500 v2 = -Math.sqrt(wb * wc); // Eqn 23 1501 o.ddata = new double[] { v1, v2, v3 }; 1502 return Integer.valueOf(-1); 1503 default: 1504 return null; 1505 } 1506 } 1507 getR0(MinBond b)1508 private double getR0(MinBond b) { 1509 return (b.ddata == null ? ((double[]) ffParams.get(b.key)) : b.ddata)[1]; 1510 } 1511 getRowFor(int i)1512 private int getRowFor(int i) { 1513 int elemno = minAtoms[i].atom.getElementNumber(); 1514 return (elemno < 3 ? 0 : elemno < 11 ? 1 : elemno < 19 ? 2 : elemno < 37 ? 3 : 4); 1515 } 1516 1517 private int[] typeData = new int[4]; 1518 getOutOfPlaneParameter(int[] data)1519 double getOutOfPlaneParameter(int[] data) { 1520 double[] ddata = (double[]) ffParams.get(getKey(data, KEY_OOP, TYPE_OOP)); 1521 return (ddata == null ? 0 : ddata[0]); 1522 } 1523 sortOop(int[] typeData)1524 private static void sortOop(int[] typeData) { 1525 fixTypeOrder(typeData, 0, 2); 1526 fixTypeOrder(typeData, 0, 3); 1527 fixTypeOrder(typeData, 2, 3); 1528 } 1529 1530 /** 1531 * 1532 * @param a 1533 * @param i 1534 * @param j 1535 * @return true if swapped; false if not 1536 */ fixTypeOrder(int[] a, int i, int j)1537 private static boolean fixTypeOrder(int[] a, int i, int j) { 1538 if (a[i] > a[j]) { 1539 swap(a, i, j); 1540 return true; 1541 } 1542 return false; 1543 } 1544 1545 /** 1546 * 1547 * @param a 1548 * @param i 1549 * @param j 1550 * @return 1 if in order, 0 if same, -1 if reversed 1551 */ fixOrder(int[] a, int i, int j)1552 private int fixOrder(int[] a, int i, int j) { 1553 int test = typeOf(a[j]) - typeOf(a[i]); 1554 if (test < 0) 1555 swap(a, i, j); 1556 return (test < 0 ? -1 : test > 0 ? 1 : 0); 1557 } 1558 swap(int[] a, int i, int j)1559 private static void swap(int[] a, int i, int j) { 1560 int t = a[i]; 1561 a[i] = a[j]; 1562 a[j] = t; 1563 } 1564 1565 private final static int[] equivalentTypes = { 1566 1, 1, // 1 1567 2, 1, // 2 1568 3, 1, // 3 1569 4, 1, // 4 1570 5, 5, // 5 1571 6, 6, // 6 1572 7, 6, // 7 1573 8, 8, // 8 1574 9, 8, // 9 1575 10, 8, // 10 1576 11, 11, // 11 1577 12, 12, // 12 1578 13, 13, // 13 1579 14, 14, // 14 1580 15, 15, // 15 1581 16, 15, // 16 1582 17, 15, // 17 1583 18, 15, // 18 1584 19, 19, // 19 1585 1, 1, // 20 1586 21, 5, // 21 1587 22, 1, // 22 1588 23, 5, // 23 1589 24, 5, // 24 1590 25, 25, // 25 1591 26, 25, // 26 1592 28, 5, // 27 1593 28, 5, // 28 1594 29, 5, // 29 1595 2, 1, // 30 1596 31, 31, // 31 1597 7, 6, // 32 1598 21, 5, // 33 1599 8, 8, // 34 1600 6, 6, // 35 1601 36, 5, // 36 1602 2, 1, // 37 1603 9, 8, // 38 1604 10, 8, // 39 1605 10, 8, // 40 1606 3, 1, // 41 1607 42, 8, // 42 1608 10, 8, // 43 1609 16, 15, // 44 1610 10, 8, // 45 1611 9, 8, // 46 1612 42, 8, // 47 1613 9, 8, // 48 1614 6, 6, // 49 1615 21, 5, // 50 1616 7, 6, // 51 1617 21, 5, // 52 1618 42, 8, // 53 1619 9, 8, // 54 1620 10, 8, // 55 1621 10, 8, // 56 1622 2, 1, // 57 1623 10, 8, // 58 1624 6, 6, // 59 1625 4, 1, // 60 1626 42, 8, // 61 1627 10, 8, // 62 1628 2, 1, // 63 1629 2, 1, // 64 1630 9, 8, // 65 1631 9, 8, // 66 1632 9, 8, // 67 1633 8, 8, // 68 1634 9, 8, // 69 1635 70, 70, // 70 1636 5, 5, // 71 1637 16, 15, // 72 1638 18, 15, // 73 1639 17, 15, // 74 1640 26, 25, // 75 1641 9, 8, // 76 1642 12, 12, // 77 1643 2, 1, // 78 1644 9, 8, // 79 1645 2, 1, // 80 1646 10, 8, // 81 1647 9, 8, // 82 1648 }; 1649 1650 /** 1651 * equivalent types for OOP and torsions 1652 * 1653 * @param type mmFF94 atom type 1654 * @param level 0, 1, or 2. 1655 * @return equivalent type or 0 1656 * 1657 */ getEquivalentType(int type, int level)1658 private static int getEquivalentType(int type, int level) { 1659 return (type == 0 ? 0 : type == 70 || type > 82 ? type : level == 2 ? 0 : 1660 equivalentTypes[((type - 1) << 1) + level]); 1661 } 1662 1663 getZParam(int elemno)1664 private static double getZParam(int elemno) { 1665 switch (elemno) { 1666 case 1: 1667 return 1.395; 1668 case 6: 1669 return 2.494; 1670 case 7: 1671 return 2.711; 1672 case 8: 1673 return 3.045; 1674 case 9: 1675 return 2.847; 1676 case 14: 1677 return 2.350; 1678 case 15: 1679 return 2.350; 1680 case 16: 1681 return 2.980; 1682 case 17: 1683 return 2.909; 1684 case 35: 1685 return 3.017; 1686 case 53: 1687 return 3.086; 1688 } 1689 return 0.0; 1690 } 1691 getCParam(int elemno)1692 private static double getCParam(int elemno) { 1693 switch (elemno) { 1694 case 5: 1695 return 0.704; 1696 case 6: 1697 return 1.016; 1698 case 7: 1699 return 1.113; 1700 case 8: 1701 return 1.337; 1702 case 14: 1703 return 0.811; 1704 case 15: 1705 return 1.068; 1706 case 16: 1707 return 1.249; 1708 case 17: 1709 return 1.078; 1710 case 33: 1711 return 0.825; 1712 } 1713 return 0.0; 1714 } 1715 getUParam(int elemno)1716 private static double getUParam(int elemno) { 1717 switch (elemno) { 1718 case 6: 1719 case 7: 1720 case 8: 1721 return 2.0; 1722 case 14: 1723 case 15: 1724 case 16: 1725 return 1.25; 1726 } 1727 return 0.0; 1728 } 1729 getVParam(int elemno)1730 private static double getVParam(int elemno) { 1731 switch (elemno) { 1732 case 6: 1733 return 2.12; 1734 case 7: 1735 return 1.5; 1736 case 8: 1737 return 0.2; 1738 case 14: 1739 return 1.22; 1740 case 15: 1741 return 2.4; 1742 case 16: 1743 return 0.49; 1744 } 1745 return 0.0; 1746 } 1747 getCovalentRadius(int elemno)1748 private static double getCovalentRadius(int elemno) { 1749 switch(elemno) { 1750 case 1: 1751 return 0.33; // corrected value from MMFF part V 1752 case 5: 1753 return 0.81; 1754 case 6: 1755 return 0.77; // corrected value from MMFF part V 1756 case 7: 1757 return 0.73; 1758 case 8: 1759 return 0.72; 1760 case 9: 1761 return 0.74; 1762 case 13: 1763 return 1.22; 1764 case 14: 1765 return 1.15; 1766 case 15: 1767 return 1.09; 1768 case 16: 1769 return 1.03; 1770 case 17: 1771 return 1.01; 1772 case 31: 1773 return 1.19; 1774 case 32: 1775 return 1.20; 1776 case 33: 1777 return 1.20; 1778 case 34: 1779 return 1.16; 1780 case 35: 1781 return 1.15; 1782 case 44: 1783 return 1.46; 1784 case 50: 1785 return 1.40; 1786 case 51: 1787 return 1.41; 1788 case 52: 1789 return 1.35; 1790 case 53: 1791 return 1.33; 1792 case 81: 1793 return 1.51; 1794 case 82: 1795 return 1.53; 1796 case 83: 1797 return 1.55; 1798 default: 1799 return Elements.getBondingRadius(elemno, 0); 1800 } 1801 } 1802 1803 // private final static double[] r0reductions = new double[] { 1804 // 0.08, 0.03, 0.10, 0.17, 0.075, 0.04 1805 // }; 1806 getRuleBondLength(MinAtom a, MinAtom b, int boAB, boolean isAromatic)1807 private static double getRuleBondLength(MinAtom a, MinAtom b, int boAB, 1808 boolean isAromatic) { 1809 1810 switch (boAB) { 1811 case 1: 1812 case 2: 1813 case 3: 1814 break; 1815 case 5: // aromatic 1816 break; 1817 default: 1818 return 0; 1819 } 1820 // Eqn 18, p. 625 1821 1822 // r0ij = r0i + r0j - c|xi - xj|^n - delta 1823 1824 int elemnoA = a.atom.getElementNumber(); 1825 int elemnoB = b.atom.getElementNumber(); 1826 double r0a = getCovalentRadius(elemnoA); 1827 double r0b = getCovalentRadius(elemnoB); 1828 double Xa = Elements.getAllredRochowElectroNeg(elemnoA); 1829 double Xb = Elements.getAllredRochowElectroNeg(elemnoB); 1830 1831 double c = (elemnoA == 1 || elemnoB == 1 ? 0.05 : 0.085); 1832 //double delta = 0.008; 1833 double n = 1.4; 1834 1835 double r = r0a + r0b; 1836 1837 if (isAromatic) 1838 boAB = (a.ffAtomType.pilp || b.ffAtomType.pilp ? 5 : 4); 1839 else 1840 switch (a.ffAtomType.mltb << 4 + b.ffAtomType.mltb) { 1841 case 0x11: 1842 boAB = 4; 1843 break; 1844 case 0x12: 1845 case 0x21: 1846 boAB = 5; 1847 break; 1848 } 1849 1850 switch (boAB) { 1851 case 1: 1852 // only single bonds involve hybridization 1853 // Halgren uses a complicated way of addressing this, 1854 // but it comes down to the fact that 1855 1856 switch (a.ffAtomType.mltb) { 1857 case 0: // sp3 "H = 3" 1858 break; 1859 case 1: 1860 case 2: 1861 //red += r0reductions[1]; // sp2 "H = 2" 1862 break; 1863 case 3: 1864 //red += r0reductions[0]; // sp "H = 1" 1865 break; 1866 } 1867 1868 // for some reason mltb for RO- is 1 1869 1870 switch (b.ffAtomType.mltb) { 1871 case 0: // sp3 "H = 3" 1872 break; 1873 case 1: 1874 case 2: 1875 //red += r0reductions[1]; // sp2 "H = 2" 1876 break; 1877 case 3: 1878 //red += r0reductions[0]; // sp "H = 1" 1879 break; 1880 } 1881 break; 1882 default: 1883 //red += 2 * r0reductions[boAB]; 1884 break; 1885 } 1886 r -= c * Math.pow(Math.abs(Xa - Xb), n); 1887 //r -= red; 1888 //r -= delta; 1889 1890 // Well, guess what? Actually red and delta are not used. 1891 1892 // -- at least not in the program run that is at the validation site: 1893 // OPTIMOL: Molecular and Macromolecular Optimization Package 17-Nov-98 16:01:23 1894 // SGI double-precision version ... Updated 5/6/98 1895 // 1896 // This calculation is run for only the following structures: 1897 // 1898 // bond red delta ro(here/valid) kb(here/valid) Etotal(here/valid) 1899 // --------------------------------------------------------------------------------------- 1900 // OHWM1 H1-O1 0.03 0.008 0.978/0.978 7.510/7.51 -21.727/-21.72690 1901 // ERULE_03 Si1-P1 0.0 0.008 2.223/2.224 1.614/1.609 -2.983/ -2.93518 1902 // ERULE_06 N1-F1 0.0 0.008 1.381/1.379 5.372/5.438 1.582/ 1.58172 1903 // 1904 // and in each case, we match the results well, but only without red and delta. 1905 1906 return r; 1907 } 1908 1909 } 1910