1""" 2 Hydrogen optimization for PDB2PQR 3 4 This is an module for hydrogen optimization routines. 5 6 ---------------------------- 7 8 PDB2PQR -- An automated pipeline for the setup, execution, and analysis of 9 Poisson-Boltzmann electrostatics calculations 10 11 Copyright (c) 2002-2011, Jens Erik Nielsen, University College Dublin; 12 Nathan A. Baker, Battelle Memorial Institute, Developed at the Pacific 13 Northwest National Laboratory, operated by Battelle Memorial Institute, 14 Pacific Northwest Division for the U.S. Department Energy.; 15 Paul Czodrowski & Gerhard Klebe, University of Marburg. 16 17 All rights reserved. 18 19 Redistribution and use in source and binary forms, with or without modification, 20 are permitted provided that the following conditions are met: 21 22 * Redistributions of source code must retain the above copyright notice, 23 this list of conditions and the following disclaimer. 24 * Redistributions in binary form must reproduce the above copyright notice, 25 this list of conditions and the following disclaimer in the documentation 26 and/or other materials provided with the distribution. 27 * Neither the names of University College Dublin, Battelle Memorial Institute, 28 Pacific Northwest National Laboratory, US Department of Energy, or University 29 of Marburg nor the names of its contributors may be used to endorse or promote 30 products derived from this software without specific prior written permission. 31 32 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 33 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 34 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 35 IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, 36 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 37 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 38 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 39 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE 40 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED 41 OF THE POSSIBILITY OF SUCH DAMAGE. 42 43 ---------------------------- 44""" 45 46import os 47import string 48import math 49 50from definitions import * 51from utilities import * 52from quatfit import * 53from routines import * 54import topology 55 56__date__ = "22 April 2009" 57__author__ = "Todd Dolinsky, Jens Erik Nielsen, Yong Huang" 58 59HDEBUG = 0 60HYDPATH = "dat/HYDROGENS.xml" 61TOPOLOGYPATH = "dat/TOPOLOGY.xml" 62ANGLE_CUTOFF = 20.0 # A - D - H(D) angle 63DIST_CUTOFF = 3.3 # H(D) to A distance 64 65class HydrogenHandler(sax.ContentHandler): 66 """ 67 Extends the SAX XML Parser to parse the Hydrogens.xml 68 class 69 """ 70 def __init__(self): 71 """ 72 Initalize the class. 73 """ 74 75 self.curelement = "" 76 self.curatom = None 77 self.curobj = None 78 self.curholder = None 79 self.map = {} 80 81 def startElement(self, name, attributes): 82 """ 83 Create optimization holder objects or atoms 84 """ 85 if name == "class": 86 obj = OptimizationHolder() 87 self.curholder = obj 88 self.curobj = obj 89 elif name == "atom": 90 obj = DefinitionAtom() 91 self.curatom = obj 92 self.curobj = obj 93 else: 94 self.curelement = name 95 return 96 97 def endElement(self, name): 98 """ 99 Complete whatever object is currently passed in 100 by the name parameter 101 """ 102 if name == "class": # Complete Residue object 103 obj = self.curholder 104 if not isinstance(obj, OptimizationHolder): 105 raise PDBInternalError("Internal error parsing XML!") 106 107 self.map[obj.name] = obj 108 self.curholder = None 109 self.curobj = None 110 111 112 elif name == "atom": # Complete atom object 113 atom = self.curatom 114 if not isinstance(atom, DefinitionAtom): 115 raise PDBInternalError("Internal error parsing XML!") 116 atomname = atom.name 117 if atomname == "": 118 raise PDBInternalError("Atom name not set in XML!") 119 else: 120 self.curholder.map[atomname] = atom 121 self.curatom = None 122 self.curobj = self.curholder 123 124 else: # Just free the current element namespace 125 self.curelement = "" 126 127 return self.map 128 129 def characters(self, text): 130 """ 131 Set a given attribute of the object to the text 132 """ 133 134 if text.isspace(): return 135 136 # If this is a float, make it so 137 try: 138 value = float(str(text)) 139 except ValueError: 140 value = str(text) 141 142 setattr(self.curobj, self.curelement, value) 143 144 145class PotentialBond: 146 """ 147 A small class containing the hbond structure 148 """ 149 def __init__(self, atom1, atom2, dist): 150 """ 151 Initialize the class 152 153 Parameters 154 atom1: The first atom in the potential bond (Atom) 155 atom2: The second atom in the potential bond (Atom) 156 dist: The distance between the two atoms (float) 157 """ 158 self.atom1 = atom1 159 self.atom2 = atom2 160 self.dist = dist 161 162 def __str__(self): 163 """ 164 String for debugging 165 """ 166 txt = "%s %s" % (self.atom1.name, self.atom1.residue) 167 txt += " to " 168 txt += "%s %s" % (self.atom2.name, self.atom2.residue) 169 txt += " (%.2f A)" % self.dist 170 return txt 171 172class hydrogenAmbiguity: 173 """ 174 A class containing information about the ambiguity 175 """ 176 def __init__(self, residue, hdef, routines): 177 """ 178 Initialize the class - if the residue has a rotateable hydrogen, 179 remove it. If it can be flipped, pre-flip the residue by creating 180 all additional atoms. 181 182 Parameters 183 residue: The residue in question (residue) 184 hdef: The hydrogen definition matching the residue 185 routines: Pointer to the general routines object 186 """ 187 self.residue = residue 188 self.hdef = hdef 189 self.routines = routines 190 191 def __str__(self): 192 """ 193 Print the ambiguity for debugging purposes 194 """ 195 text = "%s %i %s (%s)" % (self.residue.name, self.residue.resSeq, \ 196 self.residue.chainID, self.hdef.opttype) 197 return text 198 199class Optimize: 200 """ 201 The holder class for the hydrogen optimization 202 routines. Individual optimization types inherit off of this 203 class. Any functions used by multiple types appear here. 204 """ 205 def __init__(self): 206 """ 207 Initialize the class 208 """ 209 return 210 211 def __str__(self): 212 """ 213 String output for debugging 214 """ 215 txt = "%s (%s)" % (self.residue, self.optinstance.opttype) 216 return txt 217 218 def debug(self, txt): 219 """ 220 Easy way to turn on/off debugging 221 """ 222 if HDEBUG: print txt 223 224 @staticmethod 225 def getHbondangle(atom1, atom2, atom3): 226 """ 227 Get the angle between three atoms 228 229 Parameters 230 atom1: The first atom (atom) 231 atom2: The second (vertex) atom (atom) 232 atom3: The third atom (atom) 233 Returns 234 angle: The angle between the atoms (float) 235 """ 236 atom2Coords = atom2.getCoords() 237 coords1 = subtract(atom3.getCoords(), atom2Coords) 238 coords2 = subtract(atom1.getCoords(), atom2Coords) 239 norm1 = normalize(coords1) 240 norm2 = normalize(coords2) 241 dotted = dot(norm1, norm2) 242 if dotted > 1.0: # If normalized, this is due to rounding error 243 dotted = 1.0 244 elif dotted < -1.0: # If normalized, this is due to rounding error 245 dotted = -1.0 246 rad = abs(math.acos(dotted)) 247 angle = rad*180.0/math.pi 248 if angle > 180.0: 249 angle = 360.0 - angle 250 return angle 251 252 def isHbond(self, donor, acc): 253 """ 254 Determine whether this donor acceptor pair is a 255 hydrogen bond 256 """ 257 for donorhatom in donor.bonds: 258 if not donorhatom.isHydrogen(): continue 259 260 # Check the H(D)-A distance 261 262 dist = distance(donorhatom.getCoords(), acc.getCoords()) 263 264 if dist > DIST_CUTOFF: continue 265 266 # Ensure no conflicts if H(A)s if present 267 268 flag = 1 269 for acchatom in acc.bonds: 270 if not acchatom.isHydrogen(): continue 271 272 flag = 0 273 274 # Check the H(D)-H(A) distance 275 276 hdist = distance(donorhatom.getCoords(), acchatom.getCoords()) 277 if hdist < 1.5: continue 278 279 # Check the H(D)-H(A)-A angle 280 281 angle = self.getHbondangle(donorhatom, acchatom, acc) 282 if angle < 110.0: flag = 1 283 284 if flag == 0: continue 285 286 # Check the A-D-H(D) angle 287 288 angle = self.getHbondangle(acc, donor, donorhatom) 289 290 if angle <= ANGLE_CUTOFF: 291 self.debug("Found HBOND! %.4f %.4f" % (dist, angle)) 292 return 1 293 294 # If we get here, no bond is formed 295 296 return 0 297 298 @staticmethod 299 def getPairEnergy(donor, acceptor): 300 """ 301 Get the energy between two atoms 302 303 Parameters 304 donor: The first atom in the pair (Atom) 305 acceptor: The second atom in the pair (Atom) 306 Returns 307 energy: The energy of the pair (float) 308 """ 309 310 # Initialize some variables 311 bump_energy = 10.0 312 bump_distance = 1.5 313 max_hbond_energy = -10.0 314 max_ele_energy = -1.0 315 ADH_angle_cutoff = ANGLE_CUTOFF 316 DhAhA_angle_cutoff = 110.0 317 max_dha_dist = DIST_CUTOFF 318 max_ele_dist = 5.0 319 energy = 0.0 320 321 if not (donor.hdonor and acceptor.hacceptor): 322 return energy 323 324 # See if hydrogens are presently bonded to the acceptor and donor 325 326 donorhs = (bond for bond in donor.bonds if bond.isHydrogen()) 327 328 acceptorhs = [bond for bond in acceptor.bonds if bond.isHydrogen()] 329 330 for donorhatom in donorhs: 331 332 dist = distance(donorhatom.getCoords(), acceptor.getCoords()) 333 if dist > max_dha_dist and dist < max_ele_dist: 334 energy += max_ele_energy/(dist*dist) 335 continue 336 337 # Case 1: Both donor and acceptor hydrogens are present 338 for acceptorhatom in acceptorhs: 339 340 # Penalize if H(D) is too close to H(A) 341 hdist = distance(donorhatom.getCoords(), acceptorhatom.getCoords()) 342 if hdist < bump_distance: 343 energy += bump_energy 344 continue 345 346 # Assign energies based on angles 347 angle1 = Optimize.getHbondangle(acceptor, donor, donorhatom) 348 if angle1 <= ADH_angle_cutoff: 349 angle2 = Optimize.getHbondangle(donorhatom, acceptorhatom, acceptor) 350 if angle2 < DhAhA_angle_cutoff: 351 angle2 = 1.0 352 else: 353 angle2 = (DhAhA_angle_cutoff - angle2)/DhAhA_angle_cutoff 354 355 angleterm = (ADH_angle_cutoff - angle1)/ADH_angle_cutoff 356 energy += max_hbond_energy/pow(dist,3)*angleterm*angle2 357 358 # Case 2: Only donor hydrogens are present 359 if len(acceptorhs) == 0: 360 # Assign energies based on A-D-H(D) angle alone 361 angle1 = Optimize.getHbondangle(acceptor, donor, donorhatom) 362 if angle1 <= ADH_angle_cutoff: 363 angleterm = (ADH_angle_cutoff - angle1)/ADH_angle_cutoff 364 energy += max_hbond_energy/pow(dist,2)*angleterm 365 366 return energy 367 368 def makeAtomWithNoBonds(self, atom, closeatom, addname): 369 """ 370 Called for water oxygen atoms with no current bonds. 371 Uses the closeatom to place the new atom directly 372 colinear with the atom and the closeatom. 373 374 Parameters 375 atom: The oxygen atom of the water 376 closeatom: The nearby atom (donor/acceptor) 377 addname: The name of the atom to add 378 """ 379 newcoords = [] 380 residue = atom.residue 381 382 # Place along line, 1 A away 383 384 vec = subtract(closeatom.getCoords(), atom.getCoords()) 385 dist = distance(atom.getCoords(), closeatom.getCoords()) 386 387 for i in range(3): 388 newcoords.append(vec[i]/dist + atom.getCoords()[i]) 389 390 residue.createAtom(addname, newcoords) 391 newatom = residue.getAtom(addname) 392 self.routines.cells.addCell(newatom) 393 394 # Set the bonds (since not in reference structure) 395 396 if newatom not in atom.bonds: atom.bonds.append(newatom) 397 if atom not in newatom.bonds: newatom.bonds.append(atom) 398 399 def makeWaterWithOneBond(self, atom, addname): 400 """ 401 Add an atom to a water residue that already has 402 one bond. Uses the water reference structure to 403 align the new atom. 404 """ 405 residue = atom.residue 406 nextatom = atom.bonds[0] 407 coords = [atom.getCoords(), nextatom.getCoords()] 408 refcoords = [residue.reference.map[atom.name].getCoords(), \ 409 residue.reference.map["H1"].getCoords()] 410 refatomcoords = residue.reference.map["H2"].getCoords() 411 412 # Make the atom 413 414 newcoords = findCoordinates(2, coords, refcoords, refatomcoords) 415 residue.createAtom(addname, newcoords) 416 417 # Set the bonds (since not in reference structure) 418 419 newatom = residue.getAtom(addname) 420 if newatom not in atom.bonds: atom.bonds.append(newatom) 421 if atom not in newatom.bonds: newatom.bonds.append(atom) 422 423 424 def makeAtomWithOneBondH(self, atom, addname): 425 """ 426 Add a hydrogen to an alcoholic donor with one 427 existing bond. 428 """ 429 430 residue = atom.residue 431 nextatom = atom.bonds[0] 432 coords = [atom.getCoords(), nextatom.getCoords()] 433 refcoords = [residue.reference.map[atom.name].getCoords(), \ 434 residue.reference.map[nextatom.name].getCoords()] 435 refatomcoords = residue.reference.map[addname].getCoords() 436 437 # Make the atom 438 439 newcoords = findCoordinates(2, coords, refcoords, refatomcoords) 440 residue.createAtom(addname, newcoords) 441 442 def makeAtomWithOneBondLP(self, atom, addname): 443 """ 444 Add a lone pair to an alcoholic donor with one 445 existing bond. 446 """ 447 448 # Initialize some variables 449 450 residue = atom.residue 451 for refname in atom.reference.bonds: 452 if refname.startswith("H"): break 453 454 nextatom = atom.bonds[0] 455 coords = [atom.getCoords(), nextatom.getCoords()] 456 refcoords = [residue.reference.map[atom.name].getCoords(), \ 457 residue.reference.map[nextatom.name].getCoords()] 458 refatomcoords = residue.reference.map[refname].getCoords() 459 460 # Make the atom 461 462 newcoords = findCoordinates(2, coords, refcoords, refatomcoords) 463 residue.createAtom(addname, newcoords) 464 465 # Set the bonds (since not in reference structure) 466 467 newatom = residue.getAtom(addname) 468 if newatom not in atom.bonds: atom.bonds.append(newatom) 469 if atom not in newatom.bonds: newatom.bonds.append(atom) 470 471 def trySingleAlcoholicH(self, donor, acc, newatom): 472 """ 473 After a new bond has been added using 474 makeAtomWithOneBond*, try to find the best orientation 475 by rotating to form a hydrogen bond. If a bond 476 cannot be formed, remove the newatom (thereby 477 returning to a single bond). 478 """ 479 # Initialize some variables 480 481 besten = 999.99 482 bestcoords = [] 483 residue = donor.residue 484 pivot = donor.bonds[0] 485 486 for i in range(72): 487 residue.rotateTetrahedral(pivot, donor, 5.0) 488 if self.isHbond(donor, acc): 489 energy = self.getPairEnergy(donor, acc) 490 if energy < besten: 491 bestcoords = newatom.getCoords() 492 besten = energy 493 494 # If a hydrogen bond was made, set at best coordinates 495 496 if bestcoords != []: 497 newatom.x = bestcoords[0] 498 newatom.y = bestcoords[1] 499 newatom.z = bestcoords[2] 500 self.routines.cells.addCell(newatom) 501 return 1 502 else: 503 residue.removeAtom(newatom.name) 504 return 0 505 506 def trySingleAlcoholicLP(self, acc, donor, newatom): 507 """ 508 After a new bond has been added using 509 makeAtomWithOneBond*, ensure that a 510 hydrogen bond has been made. If so, try to 511 minimze the H(D)-A-LP angle. If that cannot 512 be minimized, ignore the bond and remove the 513 atom. 514 """ 515 516 # Initialize some variables 517 518 residue = acc.residue 519 pivot = acc.bonds[0] 520 bestangle = 180.00 521 bestcoords = [] 522 523 # If a hydrogen bond was made, set at best distance 524 525 if not self.isHbond(donor, acc): 526 residue.removeAtom(newatom.name) 527 return 0 528 529 # Grab the H(D) that caused the bond 530 531 for donorhatom in donor.bonds: 532 if donorhatom.isHydrogen() and \ 533 self.getHbondangle(acc, donor, donorhatom) < ANGLE_CUTOFF: break 534 535 for i in range(72): 536 residue.rotateTetrahedral(pivot, acc, 5.0) 537 angle = abs(self.getHbondangle(donorhatom, acc, newatom)) 538 539 if angle < bestangle: 540 bestangle = angle 541 bestcoords = newatom.getCoords() 542 543 # Remove if geometry does not work 544 545 if bestangle > (ANGLE_CUTOFF * 2.0): 546 self.debug("Removing due to geometry %.2f > %.2f" % (bestangle, ANGLE_CUTOFF*2.0)) 547 residue.removeAtom(newatom.name) 548 return 0 549 550 # Otherwise set to best coordinates 551 552 newatom.x = bestcoords[0] 553 newatom.y = bestcoords[1] 554 newatom.z = bestcoords[2] 555 self.routines.cells.addCell(newatom) 556 557 return 1 558 559 560 def getPositionsWithTwoBonds(self, atom): 561 """ 562 Given a tetrahedral geometry with two 563 existing bonds, return the two potential 564 sets of coordinates that are possible for 565 a new bond. 566 """ 567 568 # Initialize some variables 569 570 residue = atom.residue 571 fixed = atom.bonds[0] 572 rotate = atom.bonds[1] 573 574 # Rotate by 120 degrees twice 575 576 residue.rotateTetrahedral(fixed, atom, 120) 577 loc1 = rotate.getCoords() 578 residue.rotateTetrahedral(fixed, atom, 120) 579 loc2 = rotate.getCoords() 580 581 # Set rotate back to original by one more rotation 582 583 residue.rotateTetrahedral(fixed, atom, 120) 584 585 return loc1, loc2 586 587 def tryPositionsWithTwoBondsH(self, donor, acc, newname, loc1, loc2): 588 """ 589 Try adding a new hydrogen two the two potential 590 locations. If both form hydrogen bonds, place at 591 whatever returns the best bond as determined by 592 getPairEnergy. 593 """ 594 595 # Initialize some variables 596 597 besten = 999.99 598 bestcoords = [] 599 residue = donor.residue 600 601 # Try the first position 602 603 residue.createAtom(newname, loc1) 604 if self.isHbond(donor, acc): 605 besten = self.getPairEnergy(donor, acc) 606 bestcoords = loc1 607 608 # Try the second 609 610 newatom = residue.getAtom(newname) 611 newatom.x = loc2[0] 612 newatom.y = loc2[1] 613 newatom.z = loc2[2] 614 if self.isHbond(donor, acc): 615 energy = self.getPairEnergy(donor, acc) 616 if energy < besten: 617 bestcoords = loc2 618 619 # Set at best coords 620 621 if bestcoords != []: 622 newatom.x = bestcoords[0] 623 newatom.y = bestcoords[1] 624 newatom.z = bestcoords[2] 625 self.routines.cells.addCell(newatom) 626 return 1 627 else: 628 residue.removeAtom(newname) 629 return 0 630 631 def tryPositionsWithTwoBondsLP(self, acc, donor, newname, loc1, loc2): 632 """ 633 Try placing an LP on a tetrahedral geometry with 634 two existing bonds. If this isn't a hydrogen bond 635 it can return - otherwise ensure that the H(D)-A-LP 636 angle is minimized. 637 """ 638 639 # Initialize some variables 640 641 bestangle = 180.00 642 bestcoords = [] 643 residue = acc.residue 644 645 # If the donor/acceptor pair is not an hbond return 646 647 if not self.isHbond(donor, acc): return 0 648 649 # Grab the H(D) that caused the bond 650 651 for donorhatom in donor.bonds: 652 if donorhatom.isHydrogen() and \ 653 self.getHbondangle(acc, donor, donorhatom) < ANGLE_CUTOFF: break 654 655 # Try the first position 656 657 residue.createAtom(newname, loc1) 658 newatom = residue.getAtom(newname) 659 angle = abs(self.getHbondangle(donorhatom, acc, newatom)) 660 if angle < bestangle: 661 bestangle = angle 662 bestcoords = loc1 663 664 # Try the second 665 666 newatom.x = loc2[0] 667 newatom.y = loc2[1] 668 newatom.z = loc2[2] 669 angle = self.getHbondangle(donorhatom, acc, newatom) 670 if angle < bestangle: 671 bestcoords = loc2 672 673 # Remove if geometry does not work 674 675 if bestangle > (ANGLE_CUTOFF * 2.0): 676 residue.removeAtom(newname) 677 return 0 678 679 # Otherwise set at best coords 680 681 newatom.x = bestcoords[0] 682 newatom.y = bestcoords[1] 683 newatom.z = bestcoords[2] 684 self.routines.cells.addCell(newatom) 685 686 # Set the bonds (since not in reference structure) 687 688 if newatom not in acc.bonds: acc.bonds.append(newatom) 689 if acc not in newatom.bonds: newatom.bonds.append(acc) 690 691 return 1 692 693 def getPositionWithThreeBonds(self, atom): 694 """ 695 If there's three bonds in a tetrahedral geometry, 696 there's only one available position. Find that 697 position. 698 """ 699 700 # Initialize some variables 701 702 residue = atom.residue 703 pivot = atom.bonds[0] 704 rot1 = atom.bonds[1] 705 rot2 = atom.bonds[2] 706 707 # Find the two new positions 708 709 residue.rotateTetrahedral(pivot, atom, 120) 710 newcoords1 = rot1.getCoords() 711 residue.rotateTetrahedral(pivot, atom, 120) 712 newcoords2 = rot1.getCoords() 713 residue.rotateTetrahedral(pivot, atom, 120) 714 715 # Determine which is unoccupied 716 717 if distance(rot2.getCoords(), newcoords1) > 0.1: return newcoords1 718 else: return newcoords2 719 720 def tryPositionsWithThreeBondsH(self, donor, acc, newname, loc): 721 """ 722 Try making a hydrogen bond with the lone available 723 position. 724 """ 725 residue = donor.residue 726 residue.createAtom(newname, loc) 727 if self.isHbond(donor, acc): 728 newatom = residue.getAtom(newname) 729 self.routines.cells.addCell(newatom) 730 return 1 731 else: 732 residue.removeAtom(newname) 733 return 0 734 735 736 def tryPositionsWithThreeBondsLP(self, acc, donor, newname, loc): 737 """ 738 Try making a hydrogen bond using the lone 739 available hydrogen position. 740 """ 741 residue = acc.residue 742 if not self.isHbond(donor, acc): return 0 743 744 # Grab the H(D) that caused the bond 745 746 for donorhatom in donor.bonds: 747 if donorhatom.isHydrogen() and \ 748 self.getHbondangle(acc, donor, donorhatom) < ANGLE_CUTOFF: break 749 750 residue.createAtom(newname, loc) 751 newatom = residue.getAtom(newname) 752 753 # Remove if geometry does not work 754 755 angle = abs(self.getHbondangle(donorhatom, acc, newatom)) 756 if angle > (ANGLE_CUTOFF * 2.0): 757 residue.removeAtom(newname) 758 return 0 759 760 # Otherwise keep it 761 762 newatom = residue.getAtom(newname) 763 self.routines.cells.addCell(newatom) 764 765 # Set the bonds (since not in reference structure) 766 767 if newatom not in acc.bonds: acc.bonds.append(newatom) 768 if acc not in newatom.bonds: newatom.bonds.append(acc) 769 770 return 1 771 772class Flip(Optimize): 773 """ 774 The holder for optimization of flippable residues. 775 """ 776 def __init__(self, residue, optinstance, routines): 777 """ 778 Initialize a potential flip. Rather than flipping 779 the given residue back and forth, take each atom 780 that would be flipped and pre-flip it, making a 781 new *FLIP atom in its place. 782 783 Parameters 784 residue: The residue to flip (residue) 785 optinstance: The optimization instance containing 786 information about what to optimize 787 788 """ 789 # Initialize some variables 790 791 self.optinstance = optinstance 792 self.residue = residue 793 self.routines = routines 794 self.atomlist = [] 795 self.hbonds = [] 796 797 map = {} 798 799 # Get all moveable names for this angle/residue pair 800 801 dihedral = optinstance.optangle 802 pivot = dihedral.split()[2] 803 moveablenames = self.routines.getMoveableNames(residue, pivot) 804 # HO in CTERM shouldn't be in the list of flip atoms 805 if residue.isCterm: 806 newmoveablenames = [] 807 for i in range(len(moveablenames)): 808 if moveablenames[i] == "HO": pass 809 else: 810 newmoveablenames.append(moveablenames[i]) 811 moveablenames = newmoveablenames 812 813 # Cache current coordinates 814 815 for name in moveablenames: 816 atom = residue.getAtom(name) 817 map[name] = atom.getCoords() 818 819 # Flip the residue about the angle 820 821 anglenum = residue.reference.dihedrals.index(dihedral) 822 if anglenum == -1: 823 raise PDBInternalError("Unable to find Flip dihedral angle!") 824 825 newangle = 180.0 + residue.dihedrals[anglenum] 826 self.routines.setDihedralAngle(residue, anglenum, newangle) 827 828 # Create new atoms at cached positions 829 830 for name in map: 831 newname = "%sFLIP" % name 832 residue.createAtom(newname, map[name]) 833 newatom = residue.getAtom(newname) 834 self.routines.cells.addCell(newatom) 835 836 # Set the bonds 837 838 newatom.reference = residue.reference.map[name] 839 for bond in newatom.reference.bonds: 840 newbond = "%sFLIP" % bond 841 if residue.hasAtom(newbond): 842 bondatom = residue.map[newbond] 843 if bondatom not in newatom.bonds: newatom.bonds.append(bondatom) 844 if newatom not in bondatom.bonds: bondatom.bonds.append(newatom) 845 846 # And connect back to the existing structure 847 848 newbond = bond 849 if residue.hasAtom(newbond): 850 bondatom = residue.map[newbond] 851 if bondatom not in newatom.bonds: newatom.bonds.append(bondatom) 852 if newatom not in bondatom.bonds: bondatom.bonds.append(newatom) 853 854 residue.setDonorsAndAcceptors() 855 856 # Add to the optimization list 857 858 for name in moveablenames: 859 860 # Get the atom 861 862 atom = residue.getAtom(name) 863 if not atom.isHydrogen() and \ 864 (atom.hdonor or atom.hacceptor): self.atomlist.append(atom) 865 866 # And the FLIP 867 868 atom = residue.getAtom("%sFLIP" % name) 869 if not atom.isHydrogen() and \ 870 (atom.hdonor or atom.hacceptor): self.atomlist.append(atom) 871 872 # Special case: Neutral unassigned HIS can be acceptors 873 874 if isinstance(residue, HIS): 875 if "HIS" == residue.name and \ 876 len(residue.patches) == 1: 877 for atom in self.atomlist: 878 if atom.name.startswith("N"): 879 atom.hacceptor = 1 880 881 def tryBoth(self, donor, acc, accobj): 882 """ 883 Called when both the donor and acceptor are optimizeable; 884 If one is fixed, we only need to try one side. Otherwise 885 first try to satisfy the donor - if that's succesful, 886 try to satisfy the acceptor. An undo may be necessary 887 if the donor is satisfied and the acceptor isn't. 888 """ 889 # If one residue if fixed, use other functions 890 891 if donor.residue.fixed: 892 if accobj.tryAcceptor(acc, donor): return 1 893 else: return 0 894 if acc.residue.fixed: 895 if self.tryDonor(donor, acc): return 1 896 else: return 0 897 898 899 self.debug("Working on %s %s (donor) to %s %s (acceptor)" % \ 900 (donor.residue, donor.name, acc.residue, acc.name)) 901 if self.isHbond(donor, acc): 902 if accobj.tryAcceptor(acc, donor): 903 self.fixFlip(donor) 904 donor.hacceptor = 0 905 self.debug("NET BOND SUCCESSFUL!") 906 return 1 907 else: 908 return 0 909 910 def tryDonor(self, donor, acc): 911 """ 912 The main driver for adding a hydrogen to an 913 optimizeable residue. 914 """ 915 residue = self.residue 916 917 # Do some error checking 918 919 if not acc.hacceptor: return 0 920 921 self.debug("Working on %s %s (donor) to %s %s (acceptor)" % \ 922 (donor.residue, donor.name, acc.residue, acc.name)) 923 924 if self.isHbond(donor, acc): 925 residue.fixed = donor.name 926 self.fixFlip(donor) 927 donor.hacceptor = 0 928 return 1 929 else: 930 return 0 931 932 def tryAcceptor(self, acc, donor): 933 """ 934 The main driver for adding an LP to an optimizeable 935 residue. 936 """ 937 residue = acc.residue 938 939 # Do some error checking 940 941 if not donor.hdonor: return 0 942 943 self.debug("Working on %s %s (acceptor) to %s %s (donor)" % \ 944 (acc.residue, acc.name, donor.residue, donor.name)) 945 if self.isHbond(donor, acc): 946 residue.fixed = acc.name 947 self.fixFlip(acc) 948 acc.hdonor = 0 949 return 1 950 else: return 0 951 952 def fixFlip(self, bondatom): 953 """ 954 Called if a hydrogen bond has been found using 955 the bondatom. If bondatom is *FLIP, remove all * 956 atoms, otherwise remove all *FLIP atoms. 957 """ 958 959 # Initialize some variables 960 961 atomlist = [] 962 residue = bondatom.residue 963 for atom in residue.getAtoms(): atomlist.append(atom) 964 965 # Set a flag to see whether to delete the FLIPs or not 966 967 flag = 0 968 if bondatom.name.endswith("FLIP"): flag = 1 969 970 # Delete the appropriate atoms 971 972 for atom in atomlist: 973 atomname = atom.name 974 if atomname.endswith("FLIP") and flag: # Delete the other list 975 if residue.hasAtom(atomname[:-4]): 976 self.routines.cells.removeCell(residue.getAtom(atomname[:-4])) 977 residue.removeAtom(atomname[:-4]) 978 elif atomname.endswith("FLIP"): # Delete the flip 979 self.routines.cells.removeCell(atom) 980 residue.removeAtom(atomname) 981 else: continue 982 983 residue.fixed = 1 984 985 def finalize(self): 986 """ 987 Finalizes a flippable back to its original state - 988 since the original atoms are now *FLIP, it deletes 989 the * atoms and renames the *FLIP atoms back to *. 990 """ 991 residue = self.residue 992 993 if residue.fixed: return 994 atomlist = [] 995 for atom in residue.getAtoms(): atomlist.append(atom) 996 for atom in atomlist: 997 if atom.name.endswith("FLIP"): 998 self.routines.cells.removeCell(atom) 999 residue.removeAtom(atom.name[:-4]) 1000 residue.renameAtom(atom.name, atom.name[:-4]) 1001 residue.fixed = 1 1002 1003 def complete(self): 1004 """ 1005 Complete the flippable residue optimization. Call the finalize 1006 function, and then rename all FLIP atoms back to their standard 1007 names. 1008 """ 1009 residue = self.residue 1010 1011 self.finalize() 1012 1013 # Rename all *FLIP atoms 1014 for atom in residue.getAtoms(): 1015 atomname = atom.name 1016 if atomname.endswith("FLIP"): 1017 residue.renameAtom(atomname, atomname[:-4]) 1018 1019class Alcoholic(Optimize): 1020 """ 1021 The class for alcoholic residues 1022 """ 1023 1024 def __init__(self, residue, optinstance, routines): 1025 """ 1026 Initialize the alcoholic class by removing 1027 the alcoholic hydrogen if it exists. 1028 """ 1029 1030 # Initialize some variables 1031 1032 self.optinstance = optinstance 1033 self.residue = residue 1034 self.routines = routines 1035 self.atomlist = [] 1036 self.hbonds = [] 1037 1038 name = optinstance.map.keys()[0] 1039 self.hname = name 1040 1041 bondname = residue.reference.getAtom(name).bonds[0] 1042 self.atomlist.append(residue.getAtom(bondname)) 1043 if residue.hasAtom(name): 1044 atom = residue.getAtom(name) 1045 self.routines.cells.removeCell(atom) 1046 residue.removeAtom(name) 1047 1048 def tryBoth(self, donor, acc, accobj): 1049 """ 1050 Called when both the donor and acceptor are optimizeable; 1051 If one is fixed, we only need to try one side. Otherwise 1052 first try to satisfy the donor - if that's succesful, 1053 try to satisfy the acceptor. An undo may be necessary 1054 if the donor is satisfied and the acceptor isn't. 1055 """ 1056 # If one residue if fixed, use other functions 1057 1058 residue = donor.residue 1059 if donor.residue.fixed: 1060 if accobj.tryAcceptor(acc, donor): return 1 1061 else: return 0 1062 if acc.residue.fixed: 1063 if self.tryDonor(donor, acc): return 1 1064 else: return 0 1065 1066 if self.tryDonor(donor, acc): 1067 if accobj.tryAcceptor(acc, donor): 1068 self.debug("NET BOND SUCCESSFUL!") 1069 return 1 1070 else: # We need to remove the added H 1071 residue.removeAtom(self.hname) 1072 self.debug("REMOVED NET HBOND") 1073 return 0 1074 else: 1075 return 0 1076 1077 def tryDonor(self, donor, acc): 1078 """ 1079 The main driver for adding a hydrogen to an 1080 optimizeable residue. 1081 """ 1082 residue = self.residue 1083 1084 # Do some error checking 1085 1086 if not acc.hacceptor: return 0 1087 1088 # Get the name of the atom to add 1089 1090 newname = self.hname 1091 if residue.hasAtom(newname): return 0 1092 1093 self.debug("Working on %s %s (donor) to %s %s (acceptor)" % \ 1094 (donor.residue, donor.name, acc.residue, acc.name)) 1095 1096 # Act depending on the number of bonds 1097 1098 if len(donor.bonds) == 1: # No H or LP attached 1099 self.makeAtomWithOneBondH(donor, newname) 1100 newatom = donor.residue.getAtom(newname) 1101 return self.trySingleAlcoholicH(donor, acc, newatom) 1102 elif len(donor.bonds) == 2: 1103 loc1, loc2 = self.getPositionsWithTwoBonds(donor) 1104 return self.tryPositionsWithTwoBondsH(donor, acc, newname, loc1, loc2) 1105 elif len(donor.bonds) == 3: 1106 loc = self.getPositionWithThreeBonds(donor) 1107 return self.tryPositionsWithThreeBondsH(donor, acc, newname, loc) 1108 1109 1110 def tryAcceptor(self, acc, donor): 1111 """ 1112 The main driver for adding an LP to an optimizeable 1113 residue. 1114 """ 1115 residue = acc.residue 1116 1117 # Do some error checking 1118 1119 if not donor.hdonor: return 0 1120 1121 # Get the name of the LP to add 1122 1123 if residue.hasAtom("LP2"): return 0 1124 elif residue.hasAtom("LP1"): newname = "LP2" 1125 else: newname = "LP1" 1126 1127 self.debug("Working on %s %s (acceptor) to %s %s (donor)" % \ 1128 (acc.residue, acc.name, donor.residue, donor.name)) 1129 1130 # Act depending on the number of bonds 1131 1132 if len(acc.bonds) == 1: # No H or LP attached 1133 self.makeAtomWithOneBondLP(acc, newname) 1134 newatom = acc.residue.getAtom(newname) 1135 return self.trySingleAlcoholicLP(acc, donor, newatom) 1136 elif len(acc.bonds) == 2: 1137 loc1, loc2 = self.getPositionsWithTwoBonds(acc) 1138 return self.tryPositionsWithTwoBondsLP(acc, donor, newname, loc1, loc2) 1139 elif len(acc.bonds) == 3: 1140 loc = self.getPositionWithThreeBonds(acc) 1141 return self.tryPositionsWithThreeBondsLP(acc, donor, newname, loc) 1142 1143 def finalize(self): 1144 """ 1145 Finalize an alcoholic residue. Try to minimize 1146 conflict with nearby atoms by building away 1147 from them. Called when LPs are still present 1148 so as to account for their bonds. 1149 """ 1150 1151 # Initialize some variables 1152 1153 residue = self.residue 1154 atom = self.atomlist[0] 1155 1156 # Conditions for return 1157 1158 addname = self.hname 1159 if residue.fixed: return 1160 if residue.hasAtom(addname): return 1161 1162 if len(atom.bonds) == 1: 1163 1164 # Initialize variables 1165 1166 pivot = atom.bonds[0] 1167 #bestdist = 0.0 1168 bestcoords = [] 1169 bestenergy = 999.99 1170 1171 # Add atom and debump 1172 1173 self.makeAtomWithOneBondH(atom, addname) 1174 newatom = residue.getAtom(addname) 1175 self.routines.cells.addCell(newatom) 1176 1177 for _ in range(18): 1178 residue.rotateTetrahedral(pivot, atom, 20.0) 1179 1180 closeatoms = self.routines.cells.getNearCells(atom) 1181 energy = 0.0 1182 for catom in closeatoms: 1183 energy += self.getPairEnergy(atom, catom) + self.getPairEnergy(catom, atom) 1184 if energy < bestenergy: 1185 bestenergy = energy 1186 bestcoords = newatom.getCoords() 1187 1188 #nearatom = self.routines.getClosestAtom(newatom) 1189 1190 # If there is no closest conflict, return 1191 1192 #if nearatom == None: return 1193 1194 #dist = distance(nearatom.getCoords(), newatom.getCoords()) 1195 #if dist > bestdist: 1196 # bestdist = dist 1197 # bestcoords = newatom.getCoords() 1198 1199 if bestcoords != []: 1200 newatom.x = bestcoords[0] 1201 newatom.y = bestcoords[1] 1202 newatom.z = bestcoords[2] 1203 1204 elif len(atom.bonds) == 2: 1205 loc1, loc2 = self.getPositionsWithTwoBonds(atom) 1206 residue.createAtom(addname, loc1) 1207 newatom = residue.getAtom(addname) 1208 self.routines.cells.addCell(newatom) 1209 1210 # Debump residue if necessary by trying the other location 1211 1212 #nearatom = self.routines.getClosestAtom(newatom) 1213 #if nearatom == None: return 1214 #dist1 = distance(newatom.getCoords(), nearatom.getCoords()) 1215 1216 closeatoms = self.routines.cells.getNearCells(atom) 1217 energy1 = 0.0 1218 for catom in closeatoms: 1219 energy1 += self.getPairEnergy(atom, catom) + self.getPairEnergy(catom, atom) 1220 1221 # Place at other location 1222 1223 self.routines.cells.removeCell(newatom) 1224 newatom.x = loc2[0] 1225 newatom.y = loc2[1] 1226 newatom.z = loc2[2] 1227 self.routines.cells.addCell(newatom) 1228 1229 #nearatom = self.routines.getClosestAtom(newatom) 1230 #if nearatom == None: return 1231 1232 energy2 = 0.0 1233 for catom in closeatoms: 1234 energy2 += self.getPairEnergy(atom, catom) + self.getPairEnergy(catom, atom) 1235 1236 # If this is worse, switch back 1237 if energy2 > energy1: 1238 #if distance(newatom.getCoords(), nearatom.getCoords()) < dist1: 1239 self.routines.cells.removeCell(newatom) 1240 newatom.x = loc1[0] 1241 newatom.y = loc1[1] 1242 newatom.z = loc1[2] 1243 self.routines.cells.addCell(newatom) 1244 1245 elif len(atom.bonds) == 3: 1246 1247 loc = self.getPositionWithThreeBonds(atom) 1248 residue.createAtom(addname, loc) 1249 self.routines.cells.addCell(residue.getAtom(addname)) 1250 1251 def complete(self): 1252 """ 1253 Complete an alcoholic optimization. Call finalize(), and then 1254 remove all extra LP atoms. 1255 """ 1256 # Initialize some variables 1257 1258 residue = self.residue 1259 1260 self.finalize() 1261 residue.fixed = 1 1262 1263 # Remove all LP atoms 1264 1265 atomlist = [] 1266 for atom in residue.getAtoms(): atomlist.append(atom) 1267 for atom in atomlist: 1268 if atom.name.startswith("LP"): residue.removeAtom(atom.name) 1269 1270class Water(Optimize): 1271 """ 1272 The class for water residues 1273 """ 1274 def __init__(self, residue, optinstance, routines): 1275 """ 1276 Initialize the water optimization class 1277 """ 1278 self.optinstance = optinstance 1279 self.residue = residue 1280 self.routines = routines 1281 self.hbonds = [] 1282 1283 oxatom = residue.getAtom("O") 1284 if oxatom == None: 1285 raise PDBInternalError("Unable to find oxygen atom in %s!" % residue) 1286 1287 oxatom.hdonor = 1 1288 oxatom.hacceptor = 1 1289 1290 self.atomlist = [oxatom] 1291 1292 def tryBoth(self, donor, acc, accobj): 1293 """ 1294 Called when both the donor and acceptor are optimizeable; 1295 If one is fixed, we only need to try one side. Otherwise 1296 first try to satisfy the donor - if that's succesful, 1297 try to satisfy the acceptor. An undo may be necessary 1298 if the donor is satisfied and the acceptor isn't. 1299 """ 1300 # If one residue if fixed, use other functions 1301 1302 residue = donor.residue 1303 1304 if donor.residue.fixed: 1305 if accobj.tryAcceptor(acc, donor): return 1 1306 else: return 0 1307 if acc.residue.fixed: 1308 if self.tryDonor(donor, acc): return 1 1309 else: return 0 1310 1311 if self.tryDonor(donor, acc): 1312 if accobj.tryAcceptor(acc, donor): 1313 self.debug("NET BOND SUCCESSFUL!") 1314 return 1 1315 else: 1316 # We need to undo what we did to the donor 1317 1318 self.debug("REMOVED NET HBOND") 1319 if residue.hasAtom("H2"): residue.removeAtom("H2") 1320 elif residue.hasAtom("H1"): residue.removeAtom("H1") 1321 return 0 1322 else: 1323 return 0 1324 1325 def tryAcceptor(self, acc, donor): 1326 """ 1327 The main driver for adding an LP to an optimizeable 1328 residue. 1329 """ 1330 residue = acc.residue 1331 1332 # Do some error checking 1333 1334 if not donor.hdonor: return 0 1335 1336 1337 # Get the name of the LP to add 1338 1339 if residue.hasAtom("LP2"): return 0 1340 elif residue.hasAtom("LP1"): newname = "LP2" 1341 else: newname = "LP1" 1342 1343 self.debug("Working on %s %s (acceptor) to %s %s (donor)" % \ 1344 (acc.residue, acc.name, donor.residue, donor.name)) 1345 1346 # Act depending on the number of bonds 1347 1348 if len(acc.bonds) == 0: 1349 1350 if self.isHbond(donor, acc): 1351 1352 # Find the best donor hydrogen and use that 1353 bestdist = distance(acc.getCoords(), donor.getCoords()) 1354 for donorh in donor.bonds: 1355 dist = distance(acc.getCoords(), donorh.getCoords()) 1356 if dist < bestdist: 1357 bestdist = dist 1358 1359 # Point the LP to the best H 1360 1361 self.makeAtomWithNoBonds(acc, donorh, newname) 1362 self.debug("Added %s to %s" % (newname, acc.residue)) 1363 return 1 1364 1365 else: return 0 1366 1367 elif len(acc.bonds) == 1: # No H or LP attached 1368 self.debug("Trying to add %s to %s with one bond" % (newname, acc.residue)) 1369 self.makeWaterWithOneBond(acc, newname) 1370 newatom = acc.residue.getAtom(newname) 1371 return self.trySingleAlcoholicLP(acc, donor, newatom) 1372 elif len(acc.bonds) == 2: 1373 self.debug("Trying to add %s to %s with two bonds" % (newname, acc.residue)) 1374 loc1, loc2 = self.getPositionsWithTwoBonds(acc) 1375 return self.tryPositionsWithTwoBondsLP(acc, donor, newname, loc1, loc2) 1376 elif len(acc.bonds) == 3: 1377 self.debug("Trying to add %s to %s with three bonds" % (newname, acc.residue)) 1378 loc = self.getPositionWithThreeBonds(acc) 1379 return self.tryPositionsWithThreeBondsLP(acc, donor, newname, loc) 1380 1381 def tryDonor(self, donor, acc): 1382 """ 1383 The main driver for adding a hydrogen to an 1384 optimizeable residue. 1385 """ 1386 residue = self.residue 1387 1388 # Do some error checking 1389 1390 if not acc.hacceptor: return 0 1391 1392 # Get the name of the atom to add 1393 1394 if residue.hasAtom("H2"): return 0 1395 elif residue.hasAtom("H1"): newname = "H2" 1396 else: newname = "H1" 1397 1398 self.debug("Working on %s %s (donor) to %s %s (acceptor)" % \ 1399 (donor.residue, donor.name, acc.residue, acc.name)) 1400 1401 # Act depending on the number of bonds 1402 if len(donor.bonds) == 0: 1403 self.makeAtomWithNoBonds(donor, acc, newname) 1404 if self.isHbond(donor, acc): return 1 1405 else: 1406 self.routines.cells.removeCell(residue.getAtom(newname)) 1407 residue.removeAtom(newname) 1408 return 0 1409 if len(donor.bonds) == 1: 1410 self.makeWaterWithOneBond(donor, newname) 1411 newatom = donor.residue.getAtom(newname) 1412 return self.trySingleAlcoholicH(donor, acc, newatom) 1413 elif len(donor.bonds) == 2: 1414 loc1, loc2 = self.getPositionsWithTwoBonds(donor) 1415 return self.tryPositionsWithTwoBondsH(donor, acc, newname, loc1, loc2) 1416 elif len(donor.bonds) == 3: 1417 loc = self.getPositionWithThreeBonds(donor) 1418 return self.tryPositionsWithThreeBondsH(donor, acc, newname, loc) 1419 1420 def finalize(self): 1421 """ 1422 Finalize a water residue. Try to minimize 1423 conflict with nearby atoms by building away 1424 from them. Called when LPs are still present 1425 so as to account for their bonds. 1426 """ 1427 1428 residue = self.residue 1429 1430 # Conditions for return 1431 1432 if residue.fixed: 1433 self.debug("Residue %s already fixed" % residue) 1434 return 1435 if residue.hasAtom("H2"): 1436 self.debug("Residue %s already has H2" % residue) 1437 return 1438 1439 atom = residue.getAtom("O") 1440 if not residue.hasAtom("H1"): addname = "H1" 1441 else: addname = "H2" 1442 1443 self.debug("Finalizing %s by adding %s (%i current O bonds)" % (residue, addname, len(atom.bonds))) 1444 1445 if len(atom.bonds) == 0: 1446 1447 newcoords = [] 1448 1449 # Build hydrogen away from closest atom 1450 1451 closeatom = self.routines.getClosestAtom(atom) 1452 if closeatom != None: 1453 vec = subtract(atom.getCoords(), closeatom.getCoords()) 1454 dist = distance(atom.getCoords(), closeatom.getCoords()) 1455 1456 for i in range(3): 1457 newcoords.append(vec[i]/dist + atom.getCoords()[i]) 1458 1459 else: 1460 newcoords = add(atom.getCoords(), [1.0, 0.0, 0.0]) 1461 1462 residue.createAtom(addname, newcoords) 1463 self.routines.cells.addCell(residue.getAtom(addname)) 1464 1465 self.finalize() 1466 1467 elif len(atom.bonds) == 1: 1468 1469 # Initialize variables 1470 1471 pivot = atom.bonds[0] 1472 bestdist = 0.0 1473 bestcoords = [] 1474 1475 # Add atom and debump 1476 1477 self.makeWaterWithOneBond(atom, addname) 1478 newatom = residue.getAtom(addname) 1479 self.routines.cells.addCell(newatom) 1480 1481 for i in range(18): 1482 residue.rotateTetrahedral(pivot, atom, 20.0) 1483 nearatom = self.routines.getClosestAtom(newatom) 1484 1485 # If there is no closest conflict, continue 1486 1487 if nearatom == None: continue 1488 1489 dist = distance(nearatom.getCoords(), newatom.getCoords()) 1490 1491 if dist > bestdist: 1492 bestdist = dist 1493 bestcoords = newatom.getCoords() 1494 1495 if bestcoords != []: 1496 newatom.x = bestcoords[0] 1497 newatom.y = bestcoords[1] 1498 newatom.z = bestcoords[2] 1499 1500 if addname == "H1": 1501 self.finalize() 1502 1503 residue.fixed = 1 1504 1505 elif len(atom.bonds) == 2: 1506 1507 loc1, loc2 = self.getPositionsWithTwoBonds(atom) 1508 residue.createAtom(addname, loc1) 1509 newatom = residue.getAtom(addname) 1510 self.routines.cells.addCell(newatom) 1511 1512 # Debump residue if necessary by trying the other location 1513 1514 nearatom = self.routines.getClosestAtom(newatom) 1515 if nearatom != None: 1516 dist1 = distance(newatom.getCoords(), nearatom.getCoords()) 1517 1518 # Place at other location 1519 1520 self.routines.cells.removeCell(atom) 1521 newatom.x = loc2[0] 1522 newatom.y = loc2[1] 1523 newatom.z = loc2[2] 1524 self.routines.cells.addCell(atom) 1525 1526 nearatom = self.routines.getClosestAtom(newatom) 1527 if nearatom != None: 1528 1529 # If this is worse, switch back 1530 1531 if distance(newatom.getCoords(), nearatom.getCoords()) < dist1: 1532 self.routines.cells.removeCell(atom) 1533 newatom.x = loc1[0] 1534 newatom.y = loc1[1] 1535 newatom.z = loc1[2] 1536 self.routines.cells.addCell(atom) 1537 1538 if addname == "H1": 1539 self.finalize() 1540 1541 elif len(atom.bonds) == 3: 1542 1543 loc = self.getPositionWithThreeBonds(atom) 1544 residue.createAtom(addname, loc) 1545 self.routines.cells.addCell(residue.getAtom(addname)) 1546 1547 1548 def complete(self): 1549 """ 1550 Complete the water optimization class 1551 """ 1552 self.finalize() 1553 1554 residue = self.residue 1555 1556 # Remove all LP atoms 1557 1558 atomlist = [] 1559 for atom in residue.getAtoms(): atomlist.append(atom) 1560 for atom in atomlist: 1561 if atom.name.startswith("LP"): residue.removeAtom(atom.name) 1562 1563class Carboxylic(Optimize): 1564 """ 1565 The class for carboxylic residues 1566 """ 1567 def __init__(self, residue, optinstance, routines): 1568 """ 1569 Initialize a case where the lone hydrogen atom 1570 can have four different orientations. Works similar 1571 to initializeFlip by preadding the necessary atoms. 1572 1573 This also takes into account that the carboxyl group 1574 has different bond lengths for the two C-O bonds - 1575 this is probably due to one bond being assigned 1576 as a C=O. As a result hydrogens are only added to 1577 the C-O (longer) bond. 1578 1579 Parameters 1580 residue: The residue to flip (residue) 1581 dihedral: The angle to flip about 1582 hname: The name of one of the hydrogens to add 1583 Returns 1584 optlist: A list of optimizeable donors and 1585 acceptors in the residue (list) 1586 """ 1587 # Initialize some variables 1588 1589 self.optinstance = optinstance 1590 self.residue = residue 1591 self.routines = routines 1592 self.atomlist = [] 1593 self.hbonds = [] 1594 self.hlist = [] 1595 1596 hname2 = "" 1597 hname1 = "" 1598 for name in optinstance.map.keys(): 1599 if name.endswith("2"): hname2 = name 1600 else: hname1 = name 1601 1602 bondatom1 = residue.getAtom(optinstance.map[hname1].bond) 1603 bondatom2 = residue.getAtom(optinstance.map[hname2].bond) 1604 longflag = 0 1605 1606 # If one bond in the group is significantly (0.05 A) 1607 # longer than the other, use that group only 1608 1609 for pivotatom in bondatom1.bonds: 1610 if not pivotatom.isHydrogen(): break 1611 1612 d1 = distance(pivotatom.getCoords(), bondatom1.getCoords()) 1613 d2 = distance(pivotatom.getCoords(), bondatom2.getCoords()) 1614 1615 order = [hname1, hname2] 1616 1617 if d2 > d1 and abs(d1 - d2) > 0.05: 1618 longflag = 1 1619 order = [hname2, hname1] 1620 1621 elif d1 > d2 and abs(d1 - d2) > 0.05: 1622 longflag = 1 1623 order = [hname1, hname2] 1624 1625 1626 for hname in order: 1627 bondatom = residue.getAtom(optinstance.map[hname].bond) 1628 1629 # First mirror the hydrogen about the same donor 1630 1631 for di in residue.reference.dihedrals: 1632 if di.endswith(hname): 1633 break 1634 1635 anglenum = residue.reference.dihedrals.index(di) 1636 if anglenum == -1: 1637 raise PDBInternalError("Unable to find Carboxylic dihedral angle!") 1638 1639 if residue.dihedrals[anglenum] == None: 1640 self.atomlist.append(bondatom) 1641 continue 1642 1643 newangle = 180.0 + residue.dihedrals[anglenum] 1644 self.routines.setDihedralAngle(residue, anglenum, newangle) 1645 1646 hatom = residue.getAtom(hname) 1647 newcoords = hatom.getCoords() 1648 1649 # Flip back to return original atom 1650 1651 newangle = 180.0 + residue.dihedrals[anglenum] 1652 self.routines.setDihedralAngle(residue, anglenum, newangle) 1653 1654 # Rename the original atom and rebuild the new atom 1655 1656 residue.renameAtom(hname, "%s1" % hname) 1657 newname = "%s2" % hname 1658 residue.createAtom(newname, newcoords) 1659 newatom = residue.getAtom(newname) 1660 self.routines.cells.addCell(newatom) 1661 newatom.refdistance = hatom.refdistance 1662 1663 # Set the bonds for the new atom 1664 1665 if bondatom not in newatom.bonds: newatom.bonds.append(bondatom) 1666 if newatom not in bondatom.bonds: bondatom.bonds.append(newatom) 1667 1668 # Break if this is the only atom to add 1669 1670 self.atomlist.append(bondatom) 1671 self.hlist.append(residue.getAtom("%s1" % hname)) 1672 self.hlist.append(residue.getAtom("%s2" % hname)) 1673 1674 if longflag: break 1675 1676 residue.setDonorsAndAcceptors() 1677 1678 1679 def tryBoth(self, donor, acc, accobj): 1680 """ 1681 Called when both the donor and acceptor are optimizeable; 1682 If one is fixed, we only need to try one side. Otherwise 1683 first try to satisfy the donor - if that's succesful, 1684 try to satisfy the acceptor. An undo may be necessary 1685 if the donor is satisfied and the acceptor isn't. 1686 """ 1687 # If one residue if fixed, use other functions 1688 1689 if donor.residue.fixed: 1690 if accobj.tryAcceptor(acc, donor): return 1 1691 else: return 0 1692 if acc.residue.fixed: 1693 if self.tryDonor(donor, acc): return 1 1694 else: return 0 1695 1696 self.debug("Working on %s %s (donor) to %s %s (acceptor)" % \ 1697 (donor.residue, donor.name, acc.residue, acc.name)) 1698 1699 if self.isHbond(donor, acc): 1700 if accobj.tryAcceptor(acc, donor): 1701 self.fix(donor, acc) 1702 self.debug("NET BOND SUCCESSFUL!") 1703 return 1 1704 else: 1705 return 0 1706 1707 def isCarboxylicHbond(self, donor, acc): 1708 """ 1709 Determine whether this donor acceptor pair is a 1710 hydrogen bond 1711 """ 1712 for donorhatom in donor.bonds: 1713 if not donorhatom.isHydrogen(): continue 1714 1715 # Check the H(D)-A distance 1716 1717 dist = distance(donorhatom.getCoords(), acc.getCoords()) 1718 if dist > DIST_CUTOFF: continue 1719 1720 # Check the A-D-H(D) angle 1721 1722 angle = self.getHbondangle(acc, donor, donorhatom) 1723 if angle <= ANGLE_CUTOFF: 1724 self.debug("Found HBOND! %.4f %.4f" % (dist, angle)) 1725 return 1 1726 1727 # If we get here, no bond is formed 1728 1729 return 0 1730 1731 def tryAcceptor(self, acc, donor): 1732 """ 1733 The main driver for adding an LP to an optimizeable 1734 residue. 1735 """ 1736 residue = acc.residue 1737 1738 # Do some error checking 1739 1740 if not donor.hdonor: 1741 return 0 1742 1743 self.debug("Working on %s %s (acceptor) to %s %s (donor)" % \ 1744 (acc.residue, acc.name, donor.residue, donor.name)) 1745 1746 # We want to ignore the Hs on the acceptor 1747 1748 if self.isCarboxylicHbond(donor, acc): 1749 1750 # Eliminate the closer hydrogen 1751 1752 hyds = [] 1753 dist = None 1754 donorhatom = None 1755 for hatom in self.hlist: 1756 if hatom.isHydrogen(): hyds.append(hatom) 1757 1758 if len(hyds) < 2: 1759 return 1 1760 1761 dist = distance(hyds[0].getCoords(), donor.getCoords()) 1762 dist2 = distance(hyds[1].getCoords(), donor.getCoords()) 1763 if dist < dist2: # Eliminate hyds[0] 1764 self.hlist.remove(hyds[0]) 1765 self.routines.cells.removeCell(hyds[0]) 1766 residue.removeAtom(hyds[0].name) 1767 donorhatom = residue.getAtom(hyds[1].name) 1768 elif hyds[1] in self.hlist: 1769 self.hlist.remove(hyds[1]) 1770 self.routines.cells.removeCell(hyds[1]) 1771 residue.removeAtom(hyds[1].name) 1772 if residue.hasAtom(hyds[0].name): 1773 donorhatom = residue.getAtom(hyds[0].name) 1774 elif len(self.hlist) != 0 and residue.hasAtom(self.hlist[0].name): 1775 donorhatom = residue.getAtom(self.hlist[0].name) 1776 1777 # If only one H is left, we're done 1778 1779 if len(self.hlist) == 1: 1780 if donorhatom != None: 1781 self.rename(donorhatom) 1782 residue.fixed = 1 1783 return 1 1784 1785 else: 1786 return 0 1787 1788 1789 def tryDonor(self, donor, acc): 1790 """ 1791 The main driver for adding a hydrogen to an 1792 optimizeable residue. 1793 """ 1794 #residue = self.residue 1795 1796 # Do some error checking 1797 1798 if not acc.hacceptor: 1799 return 0 1800 1801 if self.isHbond(donor, acc): 1802 self.fix(donor, acc) 1803 return 1 1804 else: 1805 return 0 1806 1807 def fix(self, donor, acc): 1808 """ 1809 Fix the carboxylic residue. 1810 """ 1811 1812 self.debug("Fixing residue %s due to %s" % (donor.residue, donor.name)) 1813 1814 residue = donor.residue 1815 1816 # Grab the H(D) that caused the bond 1817 1818 for donorhatom in donor.bonds: 1819 if donorhatom.isHydrogen() and \ 1820 self.getHbondangle(acc, donor, donorhatom) <= ANGLE_CUTOFF: break 1821 1822 # Remove all the other available bonded hydrogens 1823 1824 hydrogens = self.hlist[:] 1825 for atom in hydrogens: 1826 if atom != donorhatom: 1827 self.routines.cells.removeCell(atom) 1828 self.hlist.remove(atom) 1829 residue.removeAtom(atom.name) 1830 1831 # Rename the atoms 1832 1833 self.rename(donorhatom) 1834 1835 residue.fixed = 1 1836 1837 def finalize(self): 1838 """ 1839 Finalize a protontated residue. Try to minimize 1840 conflict with nearby atoms. 1841 """ 1842 1843 # Initialize some variables 1844 1845 hydrogens = [] 1846 #bestdist = 0.0 1847 bestatom = None 1848 residue = self.residue 1849 1850 if residue.fixed: return 1851 1852 # For each atom, get the closest atom 1853 1854 #for hydatom in self.hlist: 1855 # closeatom = self.routines.getClosestAtom(hydatom) 1856 # dist = distance(hydatom.getCoords(), closeatom.getCoords()) 1857 # if dist > bestdist: 1858 # bestdist = dist 1859 # bestatom = hydatom 1860 1861 bestenergy = 999.99 1862 for hydatom in self.hlist: 1863 energy = 0 1864 bondedatom = hydatom.bonds[0] 1865 closeatoms = self.routines.cells.getNearCells(bondedatom) 1866 for catom in closeatoms: 1867 energy += self.getPairEnergy(bondedatom, catom) + self.getPairEnergy(catom, bondedatom) 1868 if energy < bestenergy: 1869 bestenergy = energy 1870 bestatom = hydatom 1871 1872 1873 # Keep the bestatom 1874 1875 for hydatom in self.hlist: 1876 hydrogens.append(hydatom) 1877 1878 for hydatom in hydrogens: 1879 if bestatom != hydatom: 1880 self.hlist.remove(hydatom) 1881 residue.removeAtom(hydatom.name) 1882 1883 # Rename the atoms 1884 1885 if bestatom != None and len(bestatom.name) == 4: 1886 self.rename(bestatom) 1887 else: 1888 pass 1889 residue.fixed = 1 1890 1891 def rename(self, hydatom): 1892 """ 1893 Rename the optimized atoms appropriately. This is done 1894 since the forcefields tend to require that the hydrogen is 1895 linked to a specific oxygen, and this atom may have different 1896 parameter values. 1897 1898 Parameters 1899 hydatom: The hydrogen atom that was added. (atom) 1900 """ 1901 residue = self.residue 1902 optinstance = self.optinstance 1903 1904 # No need to rename if hydatom is not in residue.map 1905 if hydatom.name not in residue.map.keys(): 1906 return 1907 # Take off the extension 1908 if len(hydatom.name) == 4: 1909 hname = hydatom.name[:-1] 1910 residue.renameAtom(hydatom.name, hname) 1911 1912 # PATCHES.xml expects *2 - if it's *1 that left, flip names 1913 1914 if len(self.atomlist) == 2: 1915 if hydatom.name.endswith("1"): 1916 residue.renameAtom(hydatom.name,"%s2" % hydatom.name[:-1]) 1917 bondname0 = self.atomlist[0].name 1918 bondname1 = self.atomlist[1].name 1919 tempname = "FLIP" 1920 residue.renameAtom(self.atomlist[0].name, tempname) 1921 residue.renameAtom(self.atomlist[1].name, bondname0) 1922 residue.renameAtom(tempname, bondname1) 1923 elif hydatom.name.endswith("OXT"): 1924 residue.renameAtom(hydatom.name,"HO") 1925 bondname0 = self.atomlist[0].name 1926 bondname1 = self.atomlist[1].name 1927 tempname = "FLIP" 1928 residue.renameAtom(self.atomlist[0].name, tempname) 1929 residue.renameAtom(self.atomlist[1].name, bondname0) 1930 residue.renameAtom(tempname, bondname1) 1931 1932 elif len(self.atomlist) == 1: 1933 1934 # Appending the other bondatom to self.atomlist 1935 hnames = [hname[:-1] + "1", hname[:-1] + "2"] 1936 for hn in hnames: 1937 bondatom = residue.getAtom(optinstance.map[hn].bond) 1938 if bondatom.name != self.atomlist[0].name: 1939 self.atomlist.append(bondatom) 1940 else: 1941 pass 1942 1943 if hydatom.name.endswith("1"): 1944 if (hydatom.name[:-1] + "2") in residue.map.keys(): 1945 residue.removeAtom("%s2" % hydatom.name[:-1]) 1946 residue.renameAtom(hydatom.name, "%s2" % hydatom.name[:-1]) 1947 bondname0 = self.atomlist[0].name 1948 bondname1 = self.atomlist[1].name 1949 tempname = "FLIP" 1950 residue.renameAtom(self.atomlist[0].name, tempname) 1951 residue.renameAtom(self.atomlist[1].name, bondname0) 1952 residue.renameAtom(tempname, bondname1) 1953 elif hydatom.name.endswith("OXT"): 1954 residue.renameAtom(hydatom.name,"HO") 1955 bondname0 = self.atomlist[0].name 1956 bondname1 = self.atomlist[1].name 1957 tempname = "FLIP" 1958 residue.renameAtom(self.atomlist[0].name, tempname) 1959 residue.renameAtom(self.atomlist[1].name, bondname0) 1960 residue.renameAtom(tempname, bondname1) 1961 1962 def complete(self): 1963 """ 1964 If not already fixed, finalize 1965 """ 1966 if len(self.hlist) == 2 and self.residue.fixed ==1: 1967 self.residue.fixed = 0 1968 if not self.residue.fixed: 1969 self.finalize() 1970 1971class Generic(Optimize): 1972 """ 1973 Generic optimization class 1974 """ 1975 def __init__(self, residue, optinstance, routines): 1976 """ 1977 Initialize the generic optimization class 1978 """ 1979 self.optinstance = optinstance 1980 self.residue = residue 1981 self.routines = routines 1982 self.hbonds = [] 1983 self.atomlist = [] 1984 1985 def finalize(self): 1986 1987 # Initialize some variables 1988 1989# hydrogens = [] 1990# bestdist = 0.0 1991# bestatom = None 1992# residue = self.residue 1993 pass 1994 1995 def complete(self): 1996 """ 1997 If not already fixed, finalize 1998 """ 1999 if not self.residue.fixed: self.finalize() 2000 2001class hydrogenRoutines: 2002 """ 2003 The main routines for hydrogen optimization. This could 2004 potentially be extended from the routines object... 2005 """ 2006 def __init__(self, routines): 2007 """ 2008 Parse the XML hydrogenFile and store the data in a map 2009 """ 2010 self.routines = routines 2011 self.protein = routines.protein 2012 self.optlist = [] 2013 self.atomlist = [] 2014 self.resmap = {} 2015 self.hydrodefs = [] 2016 2017 handler = HydrogenHandler() 2018 sax.make_parser() 2019 2020 defpath = getDatFile(HYDPATH) 2021 if defpath == "": 2022 raise PDBInternalError("Could not find %s!" % HYDPATH) 2023 2024 hydrogenFile = open(defpath) 2025 sax.parseString(hydrogenFile.read(), handler) 2026 hydrogenFile.close() 2027 2028 self.map = handler.map 2029 2030 def debug(self, text): 2031 """ 2032 Print text to stdout for debugging purposes. 2033 2034 Parameters 2035 text: The text to output (string) 2036 """ 2037 if HDEBUG: print text 2038 2039 def switchstate(self, states, amb, stateID): 2040 """ 2041 Switch a residue to a new state by first removing all 2042 hydrogens. 2043 2044 Parameters 2045 states: The list of states (list) 2046 amb : The amibiguity to switch (tuple) 2047 stateID : The state id to switch to (int) 2048 """ 2049 # 2050 # This routine does not remove all hydrogens, merely the titratable 2051 # ones as defined in HYDROGENS.DAT 2052 # 2053 # 2054 # If we do pKa calculations, then we use another routine 2055 # 2056 if states=='pKa': 2057 return self.pKa_switchstate(amb,stateID) 2058 # 2059 # JENS: From here on only for Hbond optimisation - is it's used at all? 2060 # 2061 if stateID > len(states): 2062 raise PDBInternalError("Invalid State ID!") 2063 2064 # First Remove all Hs 2065 residue = getattr(amb,"residue") 2066 hdef = getattr(amb,"hdef") 2067 #type = hdef.type 2068 for conf in hdef.conformations: 2069 hname = conf.hname 2070 boundname = conf.boundatom 2071 if residue.getAtom(hname) != None: 2072 print 'Removing',residue.name,residue.resSeq,hname 2073 residue.removeAtom(hname) 2074 residue.getAtom(boundname).hacceptor = 1 2075 residue.getAtom(boundname).hdonor = 0 2076 # Update the IntraBonds 2077 name = residue.get("name") 2078 # 2079 # Special treatment of ligands 2080 # Not sure if it's working - PC 2081 #if residue.type != 2: 2082 defresidue = self.routines.aadef.getResidue(name) 2083 residue.updateIntraBonds(defresidue) 2084 #else: 2085 ## import ligandclean.defligand ### that's not clever 2086 #x defresidue = ligandclean.defligand.LigandDefinitionResidue(residue) 2087 2088 # Now build appropriate atoms 2089 state = states[stateID] 2090 for conf in state: 2091 print conf 2092 refcoords = [] 2093 defcoords = [] 2094 defatomcoords = [] 2095 if conf == (): 2096 continue # Nothing to add 2097 hname = conf.hname 2098 for atom in conf.atoms: 2099 #print confatoms 2100 atomname = atom.get("name") 2101 2102 resatom = residue.getAtom(atomname) 2103 if atomname == hname: 2104 defatomcoords = atom.getCoords() 2105 elif resatom != None: 2106 refcoords.append(resatom.getCoords()) 2107 defcoords.append(atom.getCoords()) 2108 else: 2109 raise PDBInternalError("Could not find necessary atom!") 2110 2111 newcoords = findCoordinates(3, refcoords, defcoords, defatomcoords) 2112 boundname = conf.boundatom 2113 residue.createAtom(hname, newcoords, "ATOM") 2114 residue.addDebumpAtom(residue.getAtom(hname)) 2115 residue.getAtom(boundname).addIntraBond(hname) 2116 residue.getAtom(boundname).hacceptor = 0 2117 residue.getAtom(boundname).hdonor = 1 2118 residue.getAtom(hname).sybylType='H' # Setting the SybylType for the newly built H 2119 residue.getAtom(hname).formalcharge=0.0 # formal charge for PEOE_PB 2120 residue.getAtom(hname).titratableH=True # flag the added hydrogen 2121 residue.getAtom(hname).addIntraBond(boundname) 2122 return 2123 2124 2125 def pKa_switchstate(self,amb,stateID): 2126 """ 2127 Switch a residue to a new state by first removing all 2128 hydrogens. this routine is used in pKa calculations only! 2129 2130 Parameters 2131 amb : The amibiguity to switch (tuple) 2132 stateID : The state id to switch to (list) 2133 2134 """ 2135 # 2136 # This routine does not remove all hydrogens, merely the titratable 2137 # ones as defined in HYDROGENS.DAT 2138 # 2139 titrationdict = {'ASH1c': '1', 'ASH1t': '2', 'ASH2c': '3', 'ASH2t': '4', 'ASP': '0', 2140 'GLH1c': '1', 'GLH1t': '2', 'GLH2c': '3', 'GLH2t': '4', 'GLU': '0', 2141 'ARG0': '1+2+3+4', 'ARG': '1+2+3+4+5', 2142 'LYS': '1', 'LYS0': '0', 2143 'TYR': '1', 'TYR-': '0', 2144 'HSD': '1', 'HSE': '2', 'HSP': '1+2', 2145 'H3': '1', 'H2': '2', 'H3+H2': '1+2', 2146 'CTR01c': '1', 'CTR01t': '2', 'CTR02c': '3', 'CTR02t': '4', 'CTR-': '0'} 2147 stateID=titrationdict[stateID] 2148 stateID=stateID.split('+') 2149 new_stateID=[] 2150 for i in stateID: 2151 new_stateID.append(int(i)) 2152 # 2153 # First Remove all Hs 2154 # 2155 residue = getattr(amb,"residue") 2156 hdef = getattr(amb,"hdef") 2157 #type = hdef.opttype 2158 for conf in hdef.conformations: 2159 hname = conf.hname 2160 boundname = conf.boundatom 2161 if residue.getAtom(hname) != None: 2162 residue.removeAtom(hname) 2163 residue.getAtom(boundname).hacceptor = 1 2164 residue.getAtom(boundname).hdonor = 0 2165 # 2166 # Update the IntraBonds 2167 # 2168 #name = residue.get("name") 2169 #defresidue = self.routines.protein.referencemap[name] #self.routines.aadef.getResidue(name) 2170 #residue.updateIntraBonds(defresidue) 2171 # 2172 # Now build appropriate atoms 2173 # 2174 for stateID in new_stateID: 2175 if stateID==0: 2176 # 2177 # For state 0 we never build any hydrogens 2178 # 2179 continue 2180 conf=hdef.conformations[stateID-1] 2181 refcoords = [] 2182 defcoords = [] 2183 defatomcoords = [] 2184 if conf == (): 2185 continue # Nothing to add 2186 hname = conf.hname 2187 for atom in conf.atoms: # specail case for N-term PRO 2188 if residue.isNterm and residue.name == "PRO": 2189 if atom.name == "H": 2190 atom.name = "CD" 2191 atom.x = 1.874 2192 atom.y = 0.862 2193 atom.z = 1.306 2194 2195 if not Routines.rebuildTetrahedral(residue, hname): 2196 for atom in conf.atoms: 2197 atomname = atom.get("name") 2198 resatom = residue.getAtom(atomname) 2199 if atomname == hname: 2200 defatomcoords = atom.getCoords() 2201 elif resatom != None: 2202 refcoords.append(resatom.getCoords()) 2203 defcoords.append(atom.getCoords()) 2204 else: 2205 raise PDBInternalError("Could not find necessary atom!") 2206 2207 newcoords = findCoordinates(3, refcoords, defcoords, defatomcoords) 2208 residue.createAtom(hname, newcoords) 2209 2210 boundname = conf.boundatom 2211 residue.getAtom(boundname).hacceptor = 0 2212 residue.getAtom(boundname).hdonor = 1 2213 residue.getAtom(hname).sybylType='H' # Setting the SybylType for the newly built H 2214 residue.getAtom(hname).formalcharge=0.0 # formal charge for PEOE_PB 2215 residue.getAtom(hname).titratableH=True # flag the added hydrogen 2216 # 2217 # Update intrabonds again 2218 # 2219 if residue.isNterm and residue.name == "PRO": 2220 for atom in residue.getAtoms(): 2221 if atom.name == "H": 2222 residue.removeAtom("H") 2223 residue.update_terminus_status() 2224 return 2225 2226 def cleanup(self): 2227 """ 2228 If there are any extra carboxlyic *1 atoms, delete them. 2229 2230 This may occur when no optimization is chosen 2231 """ 2232 for residue in self.routines.protein.getResidues(): 2233 if not isinstance(residue, Amino): 2234 continue 2235 if residue.name == "GLH" or "GLH" in residue.patches: 2236 if residue.hasAtom("HE1") and residue.hasAtom("HE2"): 2237 residue.removeAtom("HE1") 2238 elif residue.name == "ASH" or "ASH" in residue.patches: 2239 if residue.hasAtom("HD1") and residue.hasAtom("HD2"): 2240 residue.removeAtom("HD1") 2241 2242 def isOptimizeable(self, residue): 2243 """ 2244 Check to see if the given residue is optimizeable 2245 There are three ways to identify a residue: 2246 2247 1. By name (i.e. HIS) 2248 2. By reference name - a PDB file HSP has 2249 a HIS reference name 2250 3. By patch - applied by PropKa, terminal selection 2251 2252 Parameters 2253 residue: The residue in question (Residue) 2254 Returns 2255 optinstance: None if not optimizeable, otherwise 2256 the OptimizationHolder instance that 2257 corresponds to the residue. 2258 """ 2259 optinstance = None 2260 if not (isinstance(residue, Amino) or isinstance(residue, WAT)): 2261 return optinstance 2262 2263 if residue.name in self.map: 2264 optinstance = self.map[residue.name] 2265 elif residue.reference.name in self.map: 2266 optinstance = self.map[residue.reference.name] 2267 else: 2268 for patch in residue.patches: 2269 if patch in self.map: 2270 optinstance = self.map[patch] 2271 break 2272 2273 # If alcoholic, make sure the hydrogen is present 2274 2275 if optinstance != None: 2276 if optinstance.opttype == "Alcoholic": 2277 atomname = optinstance.map.keys()[0] 2278 if not residue.reference.hasAtom(atomname): 2279 optinstance = None 2280 2281 return optinstance 2282 2283 def setOptimizeableHydrogens(self): 2284 """ 2285 Set any hydrogen listed in HYDROGENS.xml that 2286 is optimizeable. Used BEFORE hydrogen optimization 2287 to label atoms so that they won't be debumped - i.e. 2288 if SER HG is too close to another atom, don't debump 2289 but wait for optimization. This function should not 2290 be used if full optimization is not taking place. 2291 """ 2292 for residue in self.protein.getResidues(): 2293 optinstance = self.isOptimizeable(residue) 2294 if optinstance == None: continue 2295 for atom in residue.getAtoms(): 2296 if atom.name in optinstance.map: 2297 atom.optimizeable = 1 2298 2299 def initializeFullOptimization(self): 2300 """ 2301 Initialize the full optimization. Detects all 2302 optimizeable donors and acceptors and sets the internal 2303 optlist. 2304 """ 2305 self.routines.write("Initializing full optimization...\n") 2306 2307 # Do some setup 2308 2309 self.routines.cells = Cells(5) 2310 self.routines.cells.assignCells(self.protein) 2311 self.routines.calculateDihedralAngles() 2312 self.routines.setDonorsAndAcceptors() 2313 self.routines.updateInternalBonds() 2314 self.routines.setReferenceDistance() 2315 self.optlist = [] 2316 self.atomlist = [] 2317 2318 # First initialize the various types 2319 2320 for residue in self.protein.getResidues(): 2321 optinstance = self.isOptimizeable(residue) 2322 if isinstance(residue, Amino): 2323 if False in residue.stateboolean.values(): 2324 residue.fixed = 1 2325 else: 2326 residue.fixed = 0 2327 if optinstance == None: 2328 continue 2329 2330 type = optinstance.opttype 2331 command = "%s(residue, optinstance, self.routines)" % type 2332 if residue.fixed == 1: 2333 pass 2334 else: 2335 myobj = eval(command) 2336 self.atomlist += myobj.atomlist 2337 self.optlist.append(myobj) 2338 self.resmap[residue] = myobj 2339 2340 self.routines.write("Done.\n") 2341 2342 def initializeWaterOptimization(self): 2343 """ 2344 Initialize optimization for waters only. Detects all 2345 optimizeable donors and acceptors and sets the internal 2346 optlist. 2347 """ 2348 2349 self.routines.write("Initializing water bonding optimization...\n") 2350 2351 # Do some setup 2352 2353 self.routines.cells = Cells(5) 2354 self.routines.cells.assignCells(self.protein) 2355 self.routines.calculateDihedralAngles() 2356 self.routines.setDonorsAndAcceptors() 2357 self.routines.updateInternalBonds() 2358 self.routines.setReferenceDistance() 2359 self.optlist = [] 2360 2361 # First initialize the various types 2362 2363 for residue in self.protein.getResidues(): 2364 optinstance = self.isOptimizeable(residue) 2365 if optinstance == None: continue 2366 2367 type = optinstance.opttype 2368 if type == "Water": 2369 command = "%s(residue, optinstance, self.routines)" % type 2370 myobj = eval(command) 2371 self.atomlist += myobj.atomlist 2372 self.optlist.append(myobj) 2373 self.resmap[residue] = myobj 2374 2375 self.routines.write("Done.\n") 2376 2377 def optimizeHydrogens(self): 2378 """ 2379 The main driver for the optimization. Should be 2380 called only after the optlist has been initialized. 2381 """ 2382 2383 self.routines.write("Optimization progress:\n") 2384 2385 optlist = self.optlist 2386 resmap = {} 2387 connectivity = {} 2388 2389 # Initialize the detection progress 2390 2391 if len(optlist) == 0: return 2392 2393 2394 self.routines.write(" Detecting potential hydrogen bonds:\n") 2395 self.routines.write("0% | | 100%\n", 1) 2396 self.routines.write(" ", 1) 2397 progress = 0.0 2398 increment = 1.0/len(optlist) 2399 2400 for obj in optlist: 2401 connectivity[obj] = [] 2402 for atom in obj.atomlist: 2403 closeatoms = self.routines.cells.getNearCells(atom) 2404 for closeatom in closeatoms: 2405 2406 # Conditions for continuing 2407 2408 if atom.residue == closeatom.residue: continue 2409 if not (closeatom.hacceptor or closeatom.hdonor): continue 2410 if atom.hdonor and not atom.hacceptor \ 2411 and not closeatom.hacceptor: continue 2412 if atom.hacceptor and not atom.hdonor \ 2413 and not closeatom.hdonor: continue 2414 2415 dist = distance(atom.getCoords(), closeatom.getCoords()) 2416 if dist < 4.3: 2417 residue = atom.residue 2418 hbond = PotentialBond(atom, closeatom, dist) 2419 2420 # Store the potential bond 2421 2422 obj.hbonds.append(hbond) 2423 2424 # Keep track of connectivity 2425 2426 if closeatom in self.atomlist: 2427 closeobj = self.resmap[closeatom.residue] 2428 if closeobj not in connectivity[obj]: 2429 connectivity[obj].append(closeobj) 2430 2431 progress += increment 2432 while progress >= 0.0499: 2433 self.routines.write("*") 2434 progress -= 0.05 2435 2436 if len(optlist) > 0: self.routines.write("\n") 2437 2438 # Some residues might have no nearby hbonds - if so, place at 2439 # default state 2440 2441 for obj in optlist: 2442 if len(obj.hbonds) == 0: 2443 if obj.residue.fixed: continue 2444 self.debug("%s has no nearby partners - fixing." % obj.residue) 2445 obj.finalize() 2446 2447 # Determine the distinct networks 2448 2449 networks = [] 2450 seen = [] 2451 for obj in optlist: 2452 if obj.residue.fixed: continue 2453 if obj in seen: continue 2454 network = analyzeConnectivity(connectivity, obj) 2455 for obj in network: 2456 if obj not in seen: seen.append(obj) 2457 2458 networks.append(network) 2459 2460 # Initialize the output progress 2461 2462 if len(networks) > 0: 2463 self.routines.write(" Optimizing hydrogen bonds:\n") 2464 self.routines.write("0% | | 100%\n", 1) 2465 self.routines.write(" ", 1) 2466 progress = 0.0 2467 increment = 1.0/len(networks) 2468 2469 # Work on the networks 2470 2471 for network in networks: 2472 txt = "" 2473 for obj in network: 2474 txt += "%s, " % obj 2475 self.debug("\nStarting network %s" % txt[:-2]) 2476 2477 ### FIRST: Only optimizeable to backbone atoms 2478 2479 self.debug("* Optimizeable to backbone *") 2480 2481 hbondmap = {} 2482 for obj in network: 2483 for hbond in obj.hbonds: 2484 if hbond.atom2 not in self.atomlist: 2485 hbondmap[hbond] = hbond.dist 2486 2487 hbondlist = sortDictByValue(hbondmap) 2488 hbondlist.reverse() 2489 2490 for hbond in hbondlist: 2491 atom = hbond.atom1 2492 atom2 = hbond.atom2 2493 obj = self.resmap[atom.residue] 2494 2495 if atom.residue.fixed: 2496 continue 2497 if atom.hdonor: 2498 obj.tryDonor(atom, atom2) 2499 if atom.hacceptor: 2500 obj.tryAcceptor(atom, atom2) 2501 2502 ### SECOND: Non-dual water Optimizeable to Optimizeable 2503 2504 self.debug("\n* Optimizeable to optimizeable *") 2505 2506 hbondmap = {} 2507 seenlist = [] 2508 for obj in network: 2509 for hbond in obj.hbonds: 2510 if hbond.atom2 in self.atomlist \ 2511 and not (isinstance(hbond.atom1.residue, WAT) \ 2512 and isinstance(hbond.atom2.residue, WAT)): 2513 2514 # Only get one hbond pair 2515 2516 if not (hbond.atom2, hbond.atom1) in seenlist: 2517 hbondmap[hbond] = hbond.dist 2518 seenlist.append((hbond.atom1, hbond.atom2)) 2519 2520 hbondlist = sortDictByValue(hbondmap) 2521 hbondlist.reverse() 2522 2523 for hbond in hbondlist: 2524 atom = hbond.atom1 2525 atom2 = hbond.atom2 2526 obj1 = self.resmap[atom.residue] 2527 obj2 = self.resmap[atom2.residue] 2528 2529 # Atoms may no longer exist if already optimized 2530 2531 if not atom.residue.hasAtom(atom.name): continue 2532 if not atom2.residue.hasAtom(atom2.name): continue 2533 2534 res = 0 2535 if atom.hdonor and atom2.hacceptor: 2536 res = obj1.tryBoth(atom, atom2, obj2) 2537 2538 if atom.hacceptor and atom2.hdonor and res == 0: 2539 obj2.tryBoth(atom2, atom, obj1) 2540 2541 ### THIRD: All water-water residues 2542 2543 self.debug("\n* Water to Water *") 2544 2545 hbondmap = {} 2546 seenlist = [] 2547 for obj in network: 2548 for hbond in obj.hbonds: 2549 residue = hbond.atom1.residue 2550 if isinstance(residue, WAT) and \ 2551 isinstance(hbond.atom2.residue, WAT): 2552 if not (hbond.atom2, hbond.atom1) in seenlist: 2553 hbondmap[hbond] = hbond.dist 2554 seenlist.append((hbond.atom1, hbond.atom2)) 2555 2556 2557 hbondlist = sortDictByValue(hbondmap) 2558 hbondlist.reverse() 2559 2560 for hbond in hbondlist: 2561 atom = hbond.atom1 2562 atom2 = hbond.atom2 2563 obj1 = self.resmap[atom.residue] 2564 obj2 = self.resmap[atom2.residue] 2565 2566 res = 0 2567 if atom.hdonor and atom2.hacceptor: 2568 res = obj1.tryBoth(atom, atom2, obj2) 2569 2570 if atom.hacceptor and atom2.hdonor and res == 0: 2571 obj2.tryBoth(atom2, atom, obj1) 2572 2573 2574 ### FOURTH: Complete all residues 2575 2576 for obj in network: obj.complete() 2577 2578 # STEP 5: Update progress meter 2579 2580 progress += 100.0 * increment 2581 while progress >= 5.0: 2582 self.routines.write("*") 2583 progress -= 5.0 2584 2585 if len(networks) > 0: self.routines.write("\n") 2586 2587 def parseHydrogen(self, res): 2588 """ 2589 Parse a list of lines in order to make a hydrogen 2590 definition 2591 2592 Parameters 2593 lines: The lines to parse (list) 2594 Returns 2595 mydef: The hydrogen definition object (HydrogenDefinition) 2596 2597 This is the current definition: Name Ttyp A R # Stdconf HT Chi OPTm 2598 """ 2599 # 2600 # ------------------ 2601 # 2602 2603 toppath = getDatFile(TOPOLOGYPATH) 2604 if toppath == "": 2605 raise PDBInternalError("Could not find %s!" % TOPOLOGYPATH) 2606 2607 topfile = open(toppath) 2608 top = topology.Topology(topfile) 2609 topfile.close() 2610 2611 2612 name = self.map[res].name 2613 opttype = self.map[res].opttype 2614 optangle = self.map[res].optangle 2615 map = self.map[res].map 2616 2617 mydef = HydrogenDefinition(name, 2618 opttype, 2619 optangle, 2620 map) 2621 conf = [] 2622 patchmap = [] 2623 refmap = {} 2624 titrationstatemap = {} 2625 tautomermap = {} 2626 conformermap = {} 2627 atommap = {} 2628 2629 # reference map from TOPOLOGY.xml 2630 for res in top.residues: 2631 refmap[res.name] = res.reference 2632 for atom in refmap[res.name].atoms: 2633 atommap[res.name, atom.name] = atom 2634 for titrationstate in res.titrationStates: 2635 titrationstatemap[titrationstate.name] = titrationstate 2636 for tautomer in titrationstate.tautomers: 2637 tautomermap[tautomer.name] = tautomer 2638 for conformer in tautomer.conformers: 2639 conformermap[conformer.name] = conformer 2640 2641 if name == 'CYS': 2642 reference = refmap['CYS'] 2643 atoms = ['HG'] 2644 refatoms = ['SG', 'CB'] 2645 elif name == 'HIS': 2646 reference = refmap['HIS'] 2647 atoms = ['HD1', 'HE2'] 2648 for atom in atoms: 2649 refatoms = ['ND1', 'CG', 'CE1'] 2650 elif name == 'LYS': 2651 reference = self.routines.protein.referencemap[name] 2652 patchmap = self.routines.protein.patchmap['LYN'] 2653 atoms = patchmap.remove 2654 refatoms = ['HZ1', 'HZ2', 'NZ'] 2655 elif name == 'TYR': 2656 reference = self.routines.protein.referencemap[name] 2657 patchmap = self.routines.protein.patchmap['TYM'] 2658 atoms = patchmap.remove 2659 refatoms = ['OH', 'CZ', 'CE2'] 2660 elif name == 'WAT': 2661 reference = self.routines.protein.referencemap[name] 2662 patchmap = self.routines.protein.patchmap['HOH'] 2663 atoms = ['H1','H2'] 2664 refatoms = None 2665 elif name == 'NTR': 2666 ntrmap = {} # map for N-TERM 2667 2668 for tautomer in titrationstatemap["NTER"].tautomers: 2669 for conformer in tautomermap[tautomer.name].conformers: 2670 for conformeradds in conformermap[conformer.name].conformerAdds: 2671 for atom in conformeradds.atoms: 2672 ntrmap[atom.name] = atom 2673 atoms = ['H3', 'H2'] 2674 refatoms = ['CA', 'H', 'N'] 2675 elif name == 'CTR': 2676 hmap = {} # map for h atoms 2677 nonhmap ={} # map for refatoms 2678 conformernames = [] 2679 for tautomer in titrationstatemap["CTER"].tautomers: 2680 for conformer in tautomermap[tautomer.name].conformers: 2681 for conformeradds in conformermap[conformer.name].conformerAdds: 2682 for atom in conformeradds.atoms: 2683 nonhmap[atom.name] = atom 2684 for tautomer in titrationstatemap["CTER0"].tautomers: 2685 for conformer in tautomermap[tautomer.name].conformers: 2686 conformernames.append(conformer.name) 2687 for conformeradds in conformermap[conformer.name].conformerAdds: 2688 for atom in conformeradds.atoms: 2689 hmap[conformer.name, atom.name] = atom 2690 2691 atoms = ['HO'] 2692 refatoms = ['O', 'C', 'OXT'] 2693 elif name in ['SER', 'GLN', 'THR', 'ARG', 'ASN']: 2694 reference = refmap[name] 2695 if name == 'SER': 2696 atoms = ['HG'] 2697 refatoms = ['OG', 'CB'] 2698 elif name == 'GLN': 2699 atoms = ['HE21'] 2700 refatoms = ['NE2'] 2701 elif name == 'THR': 2702 atoms = ['HG1'] 2703 refatoms = ['OG1', 'CB'] 2704 elif name == 'ARG': 2705 atoms = ['HH11', 'HH12', 'HH21', 'HH22', 'HE'] 2706 for atom in atoms: 2707 refatoms = ['NH1', 'NH2', 'CZ'] 2708 elif name == 'ASN': 2709 atoms = ['HD21'] 2710 refatoms = ['ND2'] 2711 elif name == 'ASH': 2712 hmap = {} # map for h atoms 2713 nonhmap = {} # map for refatoms 2714 conformernames = [] 2715 reference = refmap['ASP'] 2716 for tautomer in titrationstatemap["ASH"].tautomers: 2717 for conformer in tautomermap[tautomer.name].conformers: 2718 for conformeradds in conformermap[conformer.name].conformerAdds: 2719 for atom in conformeradds.atoms: 2720 hmap[conformer.name, atom.name] = atom 2721 conformernames.append(conformer.name) 2722 atoms = ['HD1', 'HD2'] 2723 refatoms = ['OD1', 'CG', 'OD2'] 2724 elif name == 'GLH': 2725 hmap = {} # map for h atoms 2726 nonhmap ={} # map for refatoms 2727 conformernames = [] 2728 reference = refmap['GLU'] 2729 for tautomer in titrationstatemap["GLH"].tautomers: 2730 for conformer in tautomermap[tautomer.name].conformers: 2731 for conformeradds in conformermap[conformer.name].conformerAdds: 2732 for atom in conformeradds.atoms: 2733 hmap[conformer.name, atom.name] = atom 2734 conformernames.append(conformer.name) 2735 atoms = ['HE1', 'HE2'] 2736 refatoms = ['OE1', 'CD', 'OE2'] 2737 else: 2738 patchmap = self.routines.protein.patchmap[name] 2739 atoms = patchmap.map.keys() 2740 atoms.sort() 2741 2742 if name in ['NTR']: 2743 for atom in atoms: 2744 hname = atom 2745 x = ntrmap[hname].x 2746 y = ntrmap[hname].y 2747 z = ntrmap[hname].z 2748 bondatom = ntrmap[hname].bonds[0] 2749 bondlength = 1.0 2750 myconf = HydrogenConformation(hname, bondatom, bondlength) 2751 atom = DefinitionAtom(hname, x,y,z) 2752 myconf.addAtom(atom) 2753 2754 for atom in refatoms: 2755 if atom == 'N': 2756 natom = DefinitionAtom(atom, 1.201, 0.847, 0.0) 2757 myconf.addAtom(natom) 2758 elif atom == 'CA': 2759 caatom = DefinitionAtom(atom, 0.0, 0.0, 0.0) 2760 myconf.addAtom(caatom) 2761 elif atom == 'H': 2762 caatom = DefinitionAtom(atom, 1.201, 1.847, 0.000) 2763 myconf.addAtom(caatom) 2764 else: pass 2765 mydef.addConf(myconf) 2766 conf = [] 2767 2768 elif name in ['CTR']: 2769 for conformer in conformernames: 2770 for atom in atoms: 2771 hname = atom 2772 x = hmap[conformer, hname].x 2773 y = hmap[conformer, hname].y 2774 z = hmap[conformer, hname].z 2775 bondatom = hmap[conformer, hname].bonds[0] 2776 bondlength = 1.0 2777 myconf = HydrogenConformation(hname, bondatom, bondlength) 2778 atom = DefinitionAtom(hname, x,y,z) 2779 myconf.addAtom(atom) 2780 2781 for atom in refatoms: 2782 if atom == 'C': 2783 catom = DefinitionAtom(atom, -1.250, 0.881, 0.000) 2784 myconf.addAtom(catom) 2785 else: 2786 atomname = atom 2787 x = nonhmap[atom].x 2788 y = nonhmap[atom].y 2789 z = nonhmap[atom].z 2790 atom = DefinitionAtom(atomname, x,y,z) 2791 myconf.addAtom(atom) 2792 2793 mydef.addConf(myconf) 2794 2795 conf = [] 2796 2797 elif name in ['ASH', 'GLH']: 2798 for conformer in conformernames: 2799 for atom in atoms: 2800 hname = atom 2801 if ('1' in conformer and '1' in atom) or ('2' in conformer and '2' in atom): 2802 x = hmap[conformer, hname].x 2803 y = hmap[conformer, hname].y 2804 z = hmap[conformer, hname].z 2805 bondatom = hmap[conformer, hname].bonds[0] 2806 bondlength = 1.0 2807 myconf = HydrogenConformation(hname, bondatom, bondlength) 2808 atom = DefinitionAtom(hname, x,y,z) 2809 myconf.addAtom(atom) 2810 2811 for atom in refatoms: 2812 atomname = atom 2813 if name == 'ASH': 2814 refresname = 'ASP' 2815 elif name == 'GLH': 2816 refresname = 'GLU' 2817 x = atommap[refresname, atom].x 2818 y = atommap[refresname, atom].y 2819 z = atommap[refresname, atom].z 2820 atom = DefinitionAtom(atomname, x,y,z) 2821 myconf.addAtom(atom) 2822 2823 mydef.addConf(myconf) 2824 2825 conf = [] 2826 elif name in ['WAT']: 2827 pass 2828 else: 2829 for atom in atoms: 2830 hname = atom 2831 x = atommap[name, hname].x 2832 y = atommap[name, hname].y 2833 z = atommap[name, hname].z 2834 bondatom = atommap[name, hname].bonds[0] 2835 bondlength = 1.0 2836 myconf = HydrogenConformation(hname, bondatom, bondlength) 2837 atom = DefinitionAtom(hname, x,y,z) 2838 myconf.addAtom(atom) 2839 2840 if refatoms != None: 2841 if name == 'HIS' and atom.name == 'HE2': 2842 refatoms = ['NE2', 'CE1', 'CD2'] 2843 if name == 'ARG' and atom.name == 'HE': 2844 refatoms = ['NE', 'CZ', 'NH1'] 2845 for atom in refatoms: 2846 atomname = atom 2847 x = atommap[name, atomname].x 2848 y = atommap[name, atomname].y 2849 z = atommap[name, atomname].z 2850 atom = DefinitionAtom(atomname, x,y,z) 2851 myconf.addAtom(atom) 2852 2853 mydef.addConf(myconf) 2854 conf = [] 2855 return mydef 2856 2857 def readHydrogenDefinition(self): 2858 """ 2859 Read the Hydrogen Definition file 2860 2861 Returns 2862 hydrodef: The hydrogen definition () 2863 """ 2864 self.hydrodefs = [] 2865 for m in self.map: 2866 res = m 2867 mydef = self.parseHydrogen(res) 2868 self.hydrodefs.append(mydef) 2869 res = '' 2870 2871class OptimizationHolder: 2872 """ 2873 A holder class for the XML parser. 2874 """ 2875 def __init__(self): 2876 """ 2877 Initialize the class. 2878 """ 2879 self.name = "" 2880 self.map = {} 2881 self.opttype = "" 2882 self.optangle = "" 2883 2884 def __str__(self): 2885 """ 2886 A basic string representation for debugging 2887 """ 2888 text = "%s\n" % self.name 2889 text += "Type: %s\n" % self.opttype 2890 if self.optangle != "": 2891 text += "Optimization Angle: %s\n" % self.optangle 2892 text += "Atoms: \n" 2893 for atomname in self.map: 2894 text += "\t%s\n" % str(self.map[atomname]) 2895 return text 2896 2897class HydrogenDefinition: 2898 """ 2899 HydrogenDefinition class 2900 2901 The HydrogenDefinition class provides information on possible 2902 ambiguities in amino acid hydrogens. It is essentially the hydrogen 2903 definition file in object form. 2904 """ 2905 2906 def __init__(self, name, opttype, optangle, map): 2907 """ 2908 Initialize the object with information from the definition file 2909 2910 Parameters: 2911 name: The name of the grouping (string) 2912 opttype: The optimization type of the grouping (string) 2913 optangle: The optimization angle of the grouping (string) 2914 map: The map of Hydrogens. 2915 2916 See HYDROGENS.XML for more information 2917 """ 2918 self.name = name 2919 self.opttype = opttype 2920 self.optangle = optangle 2921 self.map = map 2922 self.conformations = [] 2923 2924 def __str__(self): 2925 """ 2926 Used for debugging purposes 2927 2928 Returns 2929 output: The information about this definition (string) 2930 """ 2931 output = "Name: %s\n" % self.name 2932 output += "Opttype: %s\n" % self.opttype 2933 output += "Optangle: %s\n" % self.optangle 2934 output += "map: %s\n" % self.map 2935 output += "Conformations:\n" 2936 for conf in self.conformations: 2937 output += "\n%s" % conf 2938 output += "*****************************************\n" 2939 return output 2940 2941 def addConf(self, conf): 2942 """ 2943 Add a HydrogenConformation to the list of conformations 2944 2945 Parameters 2946 conf: The conformation to be added (HydrogenConformation) 2947 """ 2948 self.conformations.append(conf) 2949 2950class HydrogenConformation: 2951 """ 2952 HydrogenConformation class 2953 2954 The HydrogenConformation class contains data about possible 2955 hydrogen conformations as specified in the hydrogen data file. 2956 """ 2957 2958 def __init__(self, hname, boundatom, bondlength): 2959 """ 2960 Initialize the object 2961 2962 Parameters 2963 hname : The hydrogen name (string) 2964 boundatom : The atom the hydrogen is bound to (string) 2965 bondlength : The bond length (float) 2966 """ 2967 self.hname = hname 2968 self.boundatom = boundatom 2969 self.bondlength = bondlength 2970 self.atoms = [] 2971 2972 def __str__(self): 2973 """ 2974 Used for debugging purposes 2975 2976 Returns 2977 output: Information about this conformation (string) 2978 """ 2979 output = "Hydrogen Name: %s\n" % self.hname 2980 output += "Bound Atom: %s\n" % self.boundatom 2981 output += "Bond Length: %.2f\n" % self.bondlength 2982 for atom in self.atoms: 2983 output += "\t%s\n" % atom 2984 return output 2985 2986 def addAtom(self, atom): 2987 """ 2988 Add an atom to the list of atoms 2989 2990 Parameters 2991 atom: The atom to be added (DefinitionAtom) 2992 """ 2993 self.atoms.append(atom) 2994