1# $Id$ 2# 3# Copyright (C) 2003-2006 greg Landrum and Rational Discovery LLC 4# 5# @@ All Rights Reserved @@ 6# This file is part of the RDKit. 7# The contents are covered by the terms of the BSD license 8# which is included in the file license.txt, found at the root 9# of the RDKit source tree. 10# 11""" Calculation of topological/topochemical descriptors. 12 13 14 15""" 16 17from rdkit import Chem 18from rdkit.Chem import Graphs 19from rdkit.Chem import rdchem 20from rdkit.Chem import rdMolDescriptors 21import numpy 22import math 23from rdkit.ML.InfoTheory import entropy 24 25ptable = Chem.GetPeriodicTable() 26 27_log2val = math.log(2) 28 29 30def _VertexDegrees(mat, onlyOnes=0): 31 """ *Internal Use Only* 32 33 this is just a row sum of the matrix... simple, neh? 34 35 """ 36 if not onlyOnes: 37 res = sum(mat) 38 else: 39 res = sum(numpy.equal(mat, 1)) 40 return res 41 42 43def _NumAdjacencies(mol, dMat): 44 """ *Internal Use Only* 45 46 """ 47 res = mol.GetNumBonds() 48 return res 49 50 51def _GetCountDict(arr): 52 """ *Internal Use Only* 53 54 """ 55 res = {} 56 for v in arr: 57 res[v] = res.get(v, 0) + 1 58 return res 59 60# WARNING: this data should probably go somewhere else... 61hallKierAlphas = {'Br': [None, None, 0.48], 62 'C': [-0.22, -0.13, 0.0], 63 'Cl': [None, None, 0.29], 64 'F': [None, None, -0.07], 65 'H': [0.0, 0.0, 0.0], 66 'I': [None, None, 0.73], 67 'N': [-0.29, -0.2, -0.04], 68 'O': [None, -0.2, -0.04], 69 'P': [None, 0.3, 0.43], 70 'S': [None, 0.22, 0.35]} 71 72def _pyHallKierAlpha(m): 73 """ calculate the Hall-Kier alpha value for a molecule 74 75 From equations (58) of Rev. Comp. Chem. vol 2, 367-422, (1991) 76 77 """ 78 alphaSum = 0.0 79 rC = ptable.GetRb0(6) 80 for atom in m.GetAtoms(): 81 atNum = atom.GetAtomicNum() 82 if not atNum: 83 continue 84 symb = atom.GetSymbol() 85 alphaV = hallKierAlphas.get(symb, None) 86 if alphaV is not None: 87 hyb = atom.GetHybridization() - 2 88 if (hyb < len(alphaV)): 89 alpha = alphaV[hyb] 90 if alpha is None: 91 alpha = alphaV[-1] 92 else: 93 alpha = alphaV[-1] 94 else: 95 rA = ptable.GetRb0(atNum) 96 alpha = rA / rC - 1 97 # print(atom.GetIdx(), atom.GetSymbol(), alpha) 98 alphaSum += alpha 99 return alphaSum 100# HallKierAlpha.version="1.0.2" 101 102 103def Ipc(mol, avg=0, dMat=None, forceDMat=0): 104 """This returns the information content of the coefficients of the characteristic 105 polynomial of the adjacency matrix of a hydrogen-suppressed graph of a molecule. 106 107 'avg = 1' returns the information content divided by the total population. 108 109 From D. Bonchev & N. Trinajstic, J. Chem. Phys. vol 67, 4517-4533 (1977) 110 111 """ 112 if forceDMat or dMat is None: 113 if forceDMat: 114 dMat = Chem.GetDistanceMatrix(mol, 0) 115 mol._adjMat = dMat 116 else: 117 try: 118 dMat = mol._adjMat 119 except AttributeError: 120 dMat = Chem.GetDistanceMatrix(mol, 0) 121 mol._adjMat = dMat 122 123 adjMat = numpy.equal(dMat, 1) 124 cPoly = abs(Graphs.CharacteristicPolynomial(mol, adjMat)) 125 if avg: 126 return entropy.InfoEntropy(cPoly) 127 else: 128 return sum(cPoly) * entropy.InfoEntropy(cPoly) 129 130 131Ipc.version = "1.0.0" 132 133 134def _pyKappa1(mol): 135 """ Hall-Kier Kappa1 value 136 137 From equations (58) and (59) of Rev. Comp. Chem. vol 2, 367-422, (1991) 138 139 """ 140 P1 = mol.GetNumBonds(1) 141 A = mol.GetNumHeavyAtoms() 142 alpha = HallKierAlpha(mol) 143 denom = P1 + alpha 144 if denom: 145 kappa = (A + alpha) * (A + alpha - 1)**2 / denom**2 146 else: 147 kappa = 0.0 148 return kappa 149# Kappa1.version="1.0.0" 150 151 152def _pyKappa2(mol): 153 """ Hall-Kier Kappa2 value 154 155 From equations (58) and (60) of Rev. Comp. Chem. vol 2, 367-422, (1991) 156 157 """ 158 P2 = len(Chem.FindAllPathsOfLengthN(mol, 2)) 159 A = mol.GetNumHeavyAtoms() 160 alpha = HallKierAlpha(mol) 161 denom = (P2 + alpha)**2 162 if denom: 163 kappa = (A + alpha - 1) * (A + alpha - 2)**2 / denom 164 else: 165 kappa = 0 166 return kappa 167# Kappa2.version="1.0.0" 168 169 170def _pyKappa3(mol): 171 """ Hall-Kier Kappa3 value 172 173 From equations (58), (61) and (62) of Rev. Comp. Chem. vol 2, 367-422, (1991) 174 175 """ 176 P3 = len(Chem.FindAllPathsOfLengthN(mol, 3)) 177 A = mol.GetNumHeavyAtoms() 178 alpha = HallKierAlpha(mol) 179 denom = (P3 + alpha)**2 180 if denom: 181 if A % 2 == 1: 182 kappa = (A + alpha - 1) * (A + alpha - 3)**2 / denom 183 else: 184 kappa = (A + alpha - 2) * (A + alpha - 3)**2 / denom 185 else: 186 kappa = 0 187 return kappa 188# Kappa3.version="1.0.0" 189 190HallKierAlpha = lambda x: rdMolDescriptors.CalcHallKierAlpha(x) 191HallKierAlpha.version = rdMolDescriptors._CalcHallKierAlpha_version 192Kappa1 = lambda x: rdMolDescriptors.CalcKappa1(x) 193Kappa1.version = rdMolDescriptors._CalcKappa1_version 194Kappa2 = lambda x: rdMolDescriptors.CalcKappa2(x) 195Kappa2.version = rdMolDescriptors._CalcKappa2_version 196Kappa3 = lambda x: rdMolDescriptors.CalcKappa3(x) 197Kappa3.version = rdMolDescriptors._CalcKappa3_version 198 199 200def Chi0(mol): 201 """ From equations (1),(9) and (10) of Rev. Comp. Chem. vol 2, 367-422, (1991) 202 203 """ 204 deltas = [x.GetDegree() for x in mol.GetAtoms()] 205 while 0 in deltas: 206 deltas.remove(0) 207 deltas = numpy.array(deltas, 'd') 208 res = sum(numpy.sqrt(1. / deltas)) 209 return res 210 211 212Chi0.version = "1.0.0" 213 214 215def Chi1(mol): 216 """ From equations (1),(11) and (12) of Rev. Comp. Chem. vol 2, 367-422, (1991) 217 218 """ 219 c1s = [x.GetBeginAtom().GetDegree() * x.GetEndAtom().GetDegree() for x in mol.GetBonds()] 220 while 0 in c1s: 221 c1s.remove(0) 222 c1s = numpy.array(c1s, 'd') 223 res = sum(numpy.sqrt(1. / c1s)) 224 return res 225 226 227Chi1.version = "1.0.0" 228 229 230def _nVal(atom): 231 return ptable.GetNOuterElecs(atom.GetAtomicNum()) - atom.GetTotalNumHs() 232 233 234def _hkDeltas(mol, skipHs=1): 235 global ptable 236 res = [] 237 if hasattr(mol, '_hkDeltas') and mol._hkDeltas is not None: 238 return mol._hkDeltas 239 for atom in mol.GetAtoms(): 240 n = atom.GetAtomicNum() 241 if n > 1: 242 nV = ptable.GetNOuterElecs(n) 243 nHs = atom.GetTotalNumHs() 244 if n <= 10: 245 # first row 246 res.append(float(nV - nHs)) 247 else: 248 # second row and up 249 res.append(float(nV - nHs) / float(n - nV - 1)) 250 elif n == 1: 251 if not skipHs: 252 res.append(0.0) 253 else: 254 res.append(0.0) 255 mol._hkDeltas = res 256 return res 257 258 259def _pyChi0v(mol): 260 """ From equations (5),(9) and (10) of Rev. Comp. Chem. vol 2, 367-422, (1991) 261 262 """ 263 deltas = _hkDeltas(mol) 264 while 0 in deltas: 265 deltas.remove(0) 266 mol._hkDeltas = None 267 res = sum(numpy.sqrt(1. / numpy.array(deltas))) 268 return res 269 270 271def _pyChi1v(mol): 272 """ From equations (5),(11) and (12) of Rev. Comp. Chem. vol 2, 367-422, (1991) 273 274 """ 275 deltas = numpy.array(_hkDeltas(mol, skipHs=0)) 276 res = 0.0 277 for bond in mol.GetBonds(): 278 v = deltas[bond.GetBeginAtomIdx()] * deltas[bond.GetEndAtomIdx()] 279 if v != 0.0: 280 res += numpy.sqrt(1. / v) 281 return res 282 283 284def _pyChiNv_(mol, order=2): 285 """ From equations (5),(15) and (16) of Rev. Comp. Chem. vol 2, 367-422, (1991) 286 287 **NOTE**: because the current path finding code does, by design, 288 detect rings as paths (e.g. in C1CC1 there is *1* atom path of 289 length 3), values of ChiNv with N >= 3 may give results that differ 290 from those provided by the old code in molecules that have rings of 291 size 3. 292 293 """ 294 deltas = numpy.array( 295 [(1. / numpy.sqrt(hkd) if hkd != 0.0 else 0.0) for hkd in _hkDeltas(mol, skipHs=0)]) 296 accum = 0.0 297 for path in Chem.FindAllPathsOfLengthN(mol, order + 1, useBonds=0): 298 accum += numpy.prod(deltas[numpy.array(path)]) 299 return accum 300 301 302def _pyChi2v(mol): 303 """ From equations (5),(15) and (16) of Rev. Comp. Chem. vol 2, 367-422, (1991) 304 305 """ 306 return _pyChiNv_(mol, 2) 307 308 309def _pyChi3v(mol): 310 """ From equations (5),(15) and (16) of Rev. Comp. Chem. vol 2, 367-422, (1991) 311 312 """ 313 return _pyChiNv_(mol, 3) 314 315 316def _pyChi4v(mol): 317 """ From equations (5),(15) and (16) of Rev. Comp. Chem. vol 2, 367-422, (1991) 318 319 **NOTE**: because the current path finding code does, by design, 320 detect rings as paths (e.g. in C1CC1 there is *1* atom path of 321 length 3), values of Chi4v may give results that differ from those 322 provided by the old code in molecules that have 3 rings. 323 324 """ 325 return _pyChiNv_(mol, 4) 326 327 328def _pyChi0n(mol): 329 """ Similar to Hall Kier Chi0v, but uses nVal instead of valence 330 This makes a big difference after we get out of the first row. 331 332 """ 333 deltas = [_nVal(x) for x in mol.GetAtoms()] 334 while deltas.count(0): 335 deltas.remove(0) 336 deltas = numpy.array(deltas, 'd') 337 res = sum(numpy.sqrt(1. / deltas)) 338 return res 339 340 341def _pyChi1n(mol): 342 """ Similar to Hall Kier Chi1v, but uses nVal instead of valence 343 344 """ 345 delts = numpy.array([_nVal(x) for x in mol.GetAtoms()], 'd') 346 res = 0.0 347 for bond in mol.GetBonds(): 348 v = delts[bond.GetBeginAtomIdx()] * delts[bond.GetEndAtomIdx()] 349 if v != 0.0: 350 res += numpy.sqrt(1. / v) 351 return res 352 353 354def _pyChiNn_(mol, order=2): 355 """ Similar to Hall Kier ChiNv, but uses nVal instead of valence 356 This makes a big difference after we get out of the first row. 357 358 **NOTE**: because the current path finding code does, by design, 359 detect rings as paths (e.g. in C1CC1 there is *1* atom path of 360 length 3), values of ChiNn with N >= 3 may give results that differ 361 from those provided by the old code in molecules that have rings of 362 size 3. 363 364 """ 365 nval = [_nVal(x) for x in mol.GetAtoms()] 366 deltas = numpy.array([(1. / numpy.sqrt(x) if x else 0.0) for x in nval]) 367 accum = 0.0 368 for path in Chem.FindAllPathsOfLengthN(mol, order + 1, useBonds=0): 369 accum += numpy.prod(deltas[numpy.array(path)]) 370 return accum 371 372 373def _pyChi2n(mol): 374 """ Similar to Hall Kier Chi2v, but uses nVal instead of valence 375 This makes a big difference after we get out of the first row. 376 377 """ 378 return _pyChiNn_(mol, 2) 379 380 381def _pyChi3n(mol): 382 """ Similar to Hall Kier Chi3v, but uses nVal instead of valence 383 This makes a big difference after we get out of the first row. 384 385 """ 386 return _pyChiNn_(mol, 3) 387 388 389def _pyChi4n(mol): 390 """ Similar to Hall Kier Chi4v, but uses nVal instead of valence 391 This makes a big difference after we get out of the first row. 392 393 394 **NOTE**: because the current path finding code does, by design, 395 detect rings as paths (e.g. in C1CC1 there is *1* atom path of 396 length 3), values of Chi4n may give results that differ from those 397 provided by the old code in molecules that have 3 rings. 398 399 """ 400 return _pyChiNn_(mol, 4) 401 402 403Chi0v = lambda x: rdMolDescriptors.CalcChi0v(x) 404Chi0v.version = rdMolDescriptors._CalcChi0v_version 405Chi1v = lambda x: rdMolDescriptors.CalcChi1v(x) 406Chi1v.version = rdMolDescriptors._CalcChi1v_version 407Chi2v = lambda x: rdMolDescriptors.CalcChi2v(x) 408Chi2v.version = rdMolDescriptors._CalcChi2v_version 409Chi3v = lambda x: rdMolDescriptors.CalcChi3v(x) 410Chi3v.version = rdMolDescriptors._CalcChi3v_version 411Chi4v = lambda x: rdMolDescriptors.CalcChi4v(x) 412Chi4v.version = rdMolDescriptors._CalcChi4v_version 413ChiNv_ = lambda x, y: rdMolDescriptors.CalcChiNv(x, y) 414ChiNv_.version = rdMolDescriptors._CalcChiNv_version 415 416Chi0n = lambda x: rdMolDescriptors.CalcChi0n(x) 417Chi0n.version = rdMolDescriptors._CalcChi0n_version 418Chi1n = lambda x: rdMolDescriptors.CalcChi1n(x) 419Chi1n.version = rdMolDescriptors._CalcChi1n_version 420Chi2n = lambda x: rdMolDescriptors.CalcChi2n(x) 421Chi2n.version = rdMolDescriptors._CalcChi2n_version 422Chi3n = lambda x: rdMolDescriptors.CalcChi3n(x) 423Chi3n.version = rdMolDescriptors._CalcChi3n_version 424Chi4n = lambda x: rdMolDescriptors.CalcChi4n(x) 425Chi4n.version = rdMolDescriptors._CalcChi4n_version 426ChiNn_ = lambda x, y: rdMolDescriptors.CalcChiNn(x, y) 427ChiNn_.version = rdMolDescriptors._CalcChiNn_version 428 429 430def BalabanJ(mol, dMat=None, forceDMat=0): 431 """ Calculate Balaban's J value for a molecule 432 433 **Arguments** 434 435 - mol: a molecule 436 437 - dMat: (optional) a distance/adjacency matrix for the molecule, if this 438 is not provide, one will be calculated 439 440 - forceDMat: (optional) if this is set, the distance/adjacency matrix 441 will be recalculated regardless of whether or not _dMat_ is provided 442 or the molecule already has one 443 444 **Returns** 445 446 - a float containing the J value 447 448 We follow the notation of Balaban's paper: 449 Chem. Phys. Lett. vol 89, 399-404, (1982) 450 451 """ 452 # if no dMat is passed in, calculate one ourselves 453 if forceDMat or dMat is None: 454 if forceDMat: 455 # FIX: should we be using atom weights here or not? 456 dMat = Chem.GetDistanceMatrix(mol, useBO=1, useAtomWts=0, force=1) 457 mol._balabanMat = dMat 458 adjMat = Chem.GetAdjacencyMatrix(mol, useBO=0, emptyVal=0, force=0, prefix="NoBO") 459 mol._adjMat = adjMat 460 else: 461 try: 462 # first check if the molecule already has one 463 dMat = mol._balabanMat 464 except AttributeError: 465 # nope, gotta calculate one 466 dMat = Chem.GetDistanceMatrix(mol, useBO=1, useAtomWts=0, force=0, prefix="Balaban") 467 # now store it 468 mol._balabanMat = dMat 469 try: 470 adjMat = mol._adjMat 471 except AttributeError: 472 adjMat = Chem.GetAdjacencyMatrix(mol, useBO=0, emptyVal=0, force=0, prefix="NoBO") 473 mol._adjMat = adjMat 474 else: 475 adjMat = Chem.GetAdjacencyMatrix(mol, useBO=0, emptyVal=0, force=0, prefix="NoBO") 476 477 s = _VertexDegrees(dMat) 478 q = _NumAdjacencies(mol, dMat) 479 n = mol.GetNumAtoms() 480 mu = q - n + 1 481 482 sum_ = 0. 483 nS = len(s) 484 for i in range(nS): 485 si = s[i] 486 for j in range(i, nS): 487 if adjMat[i, j] == 1: 488 sum_ += 1. / numpy.sqrt(si * s[j]) 489 490 if mu + 1 != 0: 491 J = float(q) / float(mu + 1) * sum_ 492 else: 493 J = 0 494 495 return J 496 497 498BalabanJ.version = "1.0.0" 499 500# ------------------------------------------------------------------------ 501# 502# Start block of BertzCT stuff. 503# 504 505 506def _AssignSymmetryClasses(mol, vdList, bdMat, forceBDMat, numAtoms, cutoff): 507 """ 508 Used by BertzCT 509 510 vdList: the number of neighbors each atom has 511 bdMat: "balaban" distance matrix 512 513 """ 514 if forceBDMat: 515 bdMat = Chem.GetDistanceMatrix(mol, useBO=1, useAtomWts=0, force=1, prefix="Balaban") 516 mol._balabanMat = bdMat 517 518 keysSeen = [] 519 symList = [0] * numAtoms 520 for i in range(numAtoms): 521 tmpList = bdMat[i].tolist() 522 tmpList.sort() 523 theKey = tuple(['%.4f' % x for x in tmpList[:cutoff]]) 524 try: 525 idx = keysSeen.index(theKey) 526 except ValueError: 527 idx = len(keysSeen) 528 keysSeen.append(theKey) 529 symList[i] = idx + 1 530 return tuple(symList) 531 532 533def _LookUpBondOrder(atom1Id, atom2Id, bondDic): 534 """ 535 Used by BertzCT 536 """ 537 if atom1Id < atom2Id: 538 theKey = (atom1Id, atom2Id) 539 else: 540 theKey = (atom2Id, atom1Id) 541 tmp = bondDic[theKey] 542 if tmp == Chem.BondType.AROMATIC: 543 tmp = 1.5 544 else: 545 tmp = float(tmp) 546 # tmp = int(tmp) 547 return tmp 548 549 550def _CalculateEntropies(connectionDict, atomTypeDict, numAtoms): 551 """ 552 Used by BertzCT 553 """ 554 connectionList = list(connectionDict.values()) 555 totConnections = sum(connectionList) 556 connectionIE = totConnections * ( 557 entropy.InfoEntropy(numpy.array(connectionList)) + math.log(totConnections) / _log2val) 558 atomTypeList = list(atomTypeDict.values()) 559 atomTypeIE = numAtoms * entropy.InfoEntropy(numpy.array(atomTypeList)) 560 return atomTypeIE + connectionIE 561 562 563def _CreateBondDictEtc(mol, numAtoms): 564 """ _Internal Use Only_ 565 Used by BertzCT 566 567 """ 568 bondDict = {} 569 nList = [None] * numAtoms 570 vdList = [0] * numAtoms 571 for aBond in mol.GetBonds(): 572 atom1 = aBond.GetBeginAtomIdx() 573 atom2 = aBond.GetEndAtomIdx() 574 if atom1 > atom2: 575 atom2, atom1 = atom1, atom2 576 if not aBond.GetIsAromatic(): 577 bondDict[(atom1, atom2)] = aBond.GetBondType() 578 else: 579 # mark Kekulized systems as aromatic 580 bondDict[(atom1, atom2)] = Chem.BondType.AROMATIC 581 if nList[atom1] is None: 582 nList[atom1] = [atom2] 583 elif atom2 not in nList[atom1]: 584 nList[atom1].append(atom2) 585 if nList[atom2] is None: 586 nList[atom2] = [atom1] 587 elif atom1 not in nList[atom2]: 588 nList[atom2].append(atom1) 589 590 for i, element in enumerate(nList): 591 try: 592 element.sort() 593 vdList[i] = len(element) 594 except Exception: 595 vdList[i] = 0 596 return bondDict, nList, vdList 597 598 599def BertzCT(mol, cutoff=100, dMat=None, forceDMat=1): 600 """ A topological index meant to quantify "complexity" of molecules. 601 602 Consists of a sum of two terms, one representing the complexity 603 of the bonding, the other representing the complexity of the 604 distribution of heteroatoms. 605 606 From S. H. Bertz, J. Am. Chem. Soc., vol 103, 3599-3601 (1981) 607 608 "cutoff" is an integer value used to limit the computational 609 expense. A cutoff value tells the program to consider vertices 610 topologically identical if their distance vectors (sets of 611 distances to all other vertices) are equal out to the "cutoff"th 612 nearest-neighbor. 613 614 **NOTE** The original implementation had the following comment: 615 > this implementation treats aromatic rings as the 616 > corresponding Kekule structure with alternating bonds, 617 > for purposes of counting "connections". 618 Upon further thought, this is the WRONG thing to do. It 619 results in the possibility of a molecule giving two different 620 CT values depending on the kekulization. For example, in the 621 old implementation, these two SMILES: 622 CC2=CN=C1C3=C(C(C)=C(C=N3)C)C=CC1=C2C 623 CC3=CN=C2C1=NC=C(C)C(C)=C1C=CC2=C3C 624 which correspond to differentk kekule forms, yield different 625 values. 626 The new implementation uses consistent (aromatic) bond orders 627 for aromatic bonds. 628 629 THIS MEANS THAT THIS IMPLEMENTATION IS NOT BACKWARDS COMPATIBLE. 630 631 Any molecule containing aromatic rings will yield different 632 values with this implementation. The new behavior is the correct 633 one, so we're going to live with the breakage. 634 635 **NOTE** this barfs if the molecule contains a second (or 636 nth) fragment that is one atom. 637 638 """ 639 atomTypeDict = {} 640 connectionDict = {} 641 numAtoms = mol.GetNumAtoms() 642 if forceDMat or dMat is None: 643 if forceDMat: 644 # nope, gotta calculate one 645 dMat = Chem.GetDistanceMatrix(mol, useBO=0, useAtomWts=0, force=1) 646 mol._adjMat = dMat 647 else: 648 try: 649 dMat = mol._adjMat 650 except AttributeError: 651 dMat = Chem.GetDistanceMatrix(mol, useBO=0, useAtomWts=0, force=1) 652 mol._adjMat = dMat 653 654 if numAtoms < 2: 655 return 0 656 657 bondDict, neighborList, vdList = _CreateBondDictEtc(mol, numAtoms) 658 symmetryClasses = _AssignSymmetryClasses(mol, vdList, dMat, forceDMat, numAtoms, cutoff) 659 # print('Symmm Classes:',symmetryClasses) 660 for atomIdx in range(numAtoms): 661 hingeAtomNumber = mol.GetAtomWithIdx(atomIdx).GetAtomicNum() 662 atomTypeDict[hingeAtomNumber] = atomTypeDict.get(hingeAtomNumber, 0) + 1 663 664 hingeAtomClass = symmetryClasses[atomIdx] 665 numNeighbors = vdList[atomIdx] 666 for i in range(numNeighbors): 667 neighbor_iIdx = neighborList[atomIdx][i] 668 NiClass = symmetryClasses[neighbor_iIdx] 669 bond_i_order = _LookUpBondOrder(atomIdx, neighbor_iIdx, bondDict) 670 # print('\t',atomIdx,i,hingeAtomClass,NiClass,bond_i_order) 671 if (bond_i_order > 1) and (neighbor_iIdx > atomIdx): 672 numConnections = bond_i_order * (bond_i_order - 1) / 2 673 connectionKey = (min(hingeAtomClass, NiClass), max(hingeAtomClass, NiClass)) 674 connectionDict[connectionKey] = connectionDict.get(connectionKey, 0) + numConnections 675 676 for j in range(i + 1, numNeighbors): 677 neighbor_jIdx = neighborList[atomIdx][j] 678 NjClass = symmetryClasses[neighbor_jIdx] 679 bond_j_order = _LookUpBondOrder(atomIdx, neighbor_jIdx, bondDict) 680 numConnections = bond_i_order * bond_j_order 681 connectionKey = (min(NiClass, NjClass), hingeAtomClass, max(NiClass, NjClass)) 682 connectionDict[connectionKey] = connectionDict.get(connectionKey, 0) + numConnections 683 684 if not connectionDict: 685 connectionDict = {'a': 1} 686 687 return _CalculateEntropies(connectionDict, atomTypeDict, numAtoms) 688 689 690BertzCT.version = "2.0.0" 691# Recent Revisions: 692# 1.0.0 -> 2.0.0: 693# - force distance matrix updates properly (Fixed as part of Issue 125) 694# - handle single-atom fragments (Issue 136) 695 696# 697# End block of BertzCT stuff. 698# 699# ------------------------------------------------------------------------ 700