Package nMOLDYN :: Package Analysis :: Module Analysis
[hide private]
[frames] | no frames]

Source Code for Module nMOLDYN.Analysis.Analysis

   1  """ 
   2  This modules implements the base class for all the analysis available in nMOLDYN. 
   3  """ 
   4   
   5  # The python distribution modules 
   6  from distutils.version import LooseVersion 
   7  import getpass 
   8  import os 
   9  import platform 
  10  from random import sample, uniform 
  11  import re 
  12  import sys 
  13  from time import asctime 
  14  from timeit import default_timer 
  15   
  16  # The ScientificPython modules 
  17  from Scientific import N 
  18  from Scientific.Geometry import Tensor, Vector 
  19   
  20  # The MMTK distribution modules 
  21  from MMTK import Atom, AtomCluster, Molecule 
  22  from MMTK.Collections import Collection 
  23  from MMTK.NucleicAcids import NucleotideChain 
  24  from MMTK.ParticleProperties import ParticleScalar 
  25  from MMTK.PDB import PDBConfiguration 
  26  from MMTK.Proteins import PeptideChain, Protein 
  27  from MMTK.Trajectory import Trajectory, TrajectorySet 
  28  from MMTK import Units 
  29   
  30  # The nMOLDYN modules 
  31  import nMOLDYN 
  32  from nMOLDYN.Core.Chemistry import belongToAnAmine, belongToAHydroxy, belongToAMethyl, belongToAThiol 
  33  from nMOLDYN.Core.Error import Error 
  34  from nMOLDYN.Core.IOFiles import TemporaryFile 
  35  from nMOLDYN.Core.Logger import LogMessage 
  36  from nMOLDYN.Core.Mathematics import randomVector 
  37  from nMOLDYN.Core.Preferences import PREFERENCES 
  38   
  39  # Dictionnary describing the residus contents of each chemical family (e.g. hydrophobic ...)  
  40  residusChemFamily = {'acidic'      : ('Asp','Glu'), 
  41                       'aliphatic'   : ('Ile','Leu','Val'), 
  42                       'aromatic'    : ('His','Phe','Trp','Tyr'), 
  43                       'basic'       : ('Arg','His','Lys'), 
  44                       'charged'     : ('Arg','Asp','Glu','His','Lys'), 
  45                       'hydrophobic' : ('Ala','Cys','Cyx','Gly','His','Ile','Leu','Lys','Met','Phe','Thr','Trp','Tyr','Val'), 
  46                       'polar'       : ('Arg','Asn','Asp','Cys','Gln','Glu','His','Lys','Ser','Thr','Trp','Tyr'), 
  47                       'small'       : ('Ala','Asn','Asp','Cys','Cyx','Gly','Pro','Ser','Thr','Val')} 
  48   
  49  # Package path 
  50  nmoldyn_package_path = os.path.dirname(os.path.split(__file__)[0]) 
  51   
52 -class Analysis:
53 """Base class for all analysis defined in nMOLDYN. 54 55 The class Analysis is an abstract-base-class that defines attributes and methods common to all the analysis 56 available in nMOLDYN. To set up an analysis object, use one of its subclass. 57 """ 58 59
60 - def __init__(self, parameters = None, statusBar = None):
61 """ 62 The constructor. 63 64 @param parameters: a dictionnary that contains parameters of the selected analysis. 65 @type parameters: dict 66 67 @param statusBar: if not None, an instance of nMOLDYN.GUI.Widgets.StatusBar. Will attach a status bar to the 68 selected analysis. 69 @type statusBar: instance of nMOLDYN.GUI.Widgets.StatusBar 70 """ 71 72 # This variable will store the the the time taken for the analysis in seconds. 73 self.chrono = None 74 75 try: 76 self.displayProgressRate = int(PREFERENCES.progress_rate) 77 78 except ValueError: 79 self.displayProgressRate = 10 80 81 if (self.displayProgressRate <= 1) or (self.displayProgressRate >= 100): 82 self.displayProgressRate = 10 83 84 # In case of an anlysis run from the GUI, this variable will be the statusbar displaying hte job progress. 85 if statusBar is None: 86 self.statusBar = None 87 else: 88 self.statusBar = statusBar 89 90 if isinstance(parameters, dict): 91 self.parameters = parameters 92 93 else: 94 self.parameters = None
95
96 - def setInputParameters(self, parameters):
97 """Sets the input parameters dictionnary. 98 """ 99 100 if isinstance(parameters, dict): 101 self.parameters = parameters 102 else: 103 raise Error('Bad input for parameters constructor argument. It must be a dictionnary')
104 105
106 - def parseInputParameters(self):
107 """Parses the input parameters stored in |parameters| dictionnary. 108 109 @return: a dictionnary of the parsed parameters. 110 @rtype: dict 111 """ 112 113 # First check for the occurence of some general parameters. 114 if not self.parameters.has_key('version'): 115 self.parameters['version'] = 'unknown' 116 117 if self.canBeEstimated: 118 if not self.parameters.has_key('estimate'): 119 self.parameters['estimate'] = 'no' 120 else: 121 if self.parameters.has_key('estimate'): 122 if self.parameters['estimate'] == 'yes': 123 raise Error('%s analysis can not be estimated.' % self.__class__.__name__.split('_')[0]) 124 else: 125 self.parameters['estimate'] = 'no' 126 127 # First, checks whether some input parameters are not used by the analysis. If so, just raises a warning. 128 for p in self.parameters.keys(): 129 if p.lower() in ['version','estimate']: 130 continue 131 132 if p.lower() not in [v.lower() for v in self.inputParametersNames]: 133 LogMessage('warning','%s is not a useful parameter for %s analysis.' % (p, self.__class__.__name__),['file','console']) 134 del self.parameters[p] 135 136 # Then, checks whether some analysis necessary parameters are missing from the input parameters. 137 # If so, raises an error. 138 for p in self.inputParametersNames: 139 if p.lower() not in [k.lower() for k in self.parameters.keys()]: 140 raise Error('%s input parameter is missing for %s analysis.' % (p, self.__class__.__name__)) 141 142 # Loop over the parameters names necessary for the analysis. 143 for pName in self.parameters: 144 145 # Their corresponding input value. 146 pValue = self.parameters[pName] 147 pName = pName.lower() 148 149 if pName == 'armodelorder': 150 self.arModelOrder = int(pValue) 151 152 elif pName == 'atomorder': 153 154 if pValue is None: 155 self.atomOrder = None 156 157 elif isinstance(pValue, str): 158 159 if pValue.strip().lower() == 'no': 160 self.atomOrder = None 161 162 else: 163 try: 164 self.atomOrder = [p.strip() for p in pValue.split(',')] 165 if not self.atomOrder: 166 raise 167 168 except: 169 raise Error('Error when parsing % parameter: must be comma-separated strings.' % pName) 170 171 else: 172 raise Error('Error when parsing %s parameter: must be a string.' % pName) 173 174 elif pName == 'bond': 175 176 if pValue is None: 177 self.bond = 'all' 178 else: 179 self.bond = pValue 180 181 elif pName == 'deuteration': 182 183 if pValue is None: 184 self.deuteration = 'no' 185 else: 186 self.deuteration = pValue 187 188 elif pName == 'differentiation': 189 190 try: 191 self.differentiation = int(pValue) 192 193 except: 194 raise Error('Error when parsing %s parameter: must be an integer.' % pName) 195 196 elif pName == 'estimate': 197 try: 198 p = pValue.strip().lower() 199 if p == 'yes': 200 self.estimate = True 201 elif p == 'no': 202 self.estimate = False 203 else: 204 raise 205 206 except: 207 raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName) 208 209 # 'fftWindow' input parameter name. 210 elif pName == 'fftwindow': 211 try: 212 self.fftWindow = float(pValue) 213 if (self.fftWindow <= 0.0) or (self.fftWindow >= 100.0): 214 raise 215 216 except: 217 raise Error('Error when parsing %s parameter: must be a float in ]0.0,100.0[.' % pName) 218 219 # 'filter' input parameter name. 220 elif pName == 'filter': 221 try: 222 self.freqMin, self.freqMax = [float(v) for v in pValue.split(':')] 223 224 if self.freqMin >= self.freqMax: 225 raise 226 227 except: 228 raise Error('Error when parsing %s parameter: must be "float:float".' % pName) 229 230 elif pName == 'group': 231 232 if pValue is None: 233 self.group = 'all' 234 else: 235 self.group = pValue 236 237 elif pName == 'normalize': 238 try: 239 p = pValue.strip().lower() 240 if p == 'yes': 241 self.normalize = True 242 elif p == 'no': 243 self.normalize = False 244 else: 245 raise 246 247 except: 248 raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName) 249 250 elif pName == 'phivalues': 251 252 if not isinstance(pValue, str): 253 raise Error('Error when parsing %s parameter: must be a string.' % pName) 254 255 try: 256 257 self.phiMin, self.phiMax, self.dPhi = [float(v) for v in pValue.split(':')] 258 259 if self.phiMin < -180.0: 260 raise 261 262 if self.phiMax > 180.0: 263 raise 264 265 if self.phiMax <= self.phiMin: 266 raise 267 268 if self.dPhi > (self.phiMax - self.phiMin): 269 raise 270 271 self.phiMin *= Units.deg 272 self.phiMax *= Units.deg 273 self.dPhi *= Units.deg 274 275 phi = [] 276 phi0 = self.phiMin 277 while phi0 <= self.phiMax + 1.0e-6: 278 phi.append(phi0) 279 phi0 += self.dPhi 280 281 if not phi: 282 raise 283 284 self.nPhiBins = len(phi) - 1 285 286 # The rMax is readjusted to fit the last r values of the list. 287 self.phiMax = phi[-1] 288 289 self.phiValues = self.phiMin + self.dPhi/2.0 + self.dPhi*N.arange(self.nPhiBins) 290 291 except: 292 raise Error('Error when parsing %s parameter: wrong format for the string.' % pName) 293 294 # 'projection' input parameter name. 295 # The format must be 'float,float,float' 296 elif pName == 'projection': 297 298 if pValue is None: 299 self.projection = None 300 301 elif isinstance(pValue, str): 302 303 if pValue.strip().lower() == 'no': 304 self.projection = None 305 306 else: 307 308 try: 309 self.projection = Vector([float(v) for v in pValue.split(',')]).normal() 310 311 except: 312 raise Error('Error when parsing % parameter: must be "float,float,float".' % pName) 313 314 else: 315 raise Error('Error when parsing %s parameter: must be a string.' % pName) 316 317 elif pName == 'pyroserver': 318 319 if pValue is None: 320 self.pyroServer = 'monoprocessor' 321 self.pyroNodes = {} 322 323 elif isinstance(pValue, str): 324 325 try: 326 if pValue.strip().lower() in ['monoprocessor', 'no']: 327 self.pyroServer = 'monoprocessor' 328 self.pyroNodes = {} 329 330 elif pValue[:14].strip().lower() == 'multiprocessor': 331 hostName, nAllocatedProcs = pValue.split('::')[1].split(':') 332 self.pyroServer = 'multiprocessor' 333 self.pyroNodes = {hostName.strip() : int(nAllocatedProcs.strip())} 334 335 elif pValue[:7].strip().lower() == 'cluster': 336 raise Error('Pyroserver in cluster configuration not implemented yet.') 337 self.pyroServer = 'cluster' 338 self.pyroNodes = {} 339 self.pyroNodes.update([[vv[0].strip(),int(vv[1].strip())] for vv in [v.split(':') for v in pValue.split('::')[1].split(',')]]) 340 341 except: 342 raise Error('Error when parsing %s parameter string in %s case.' % (pName,pValue.strip().lower())) 343 344 # 'qShellValues' input parameter name. 345 # It can be either directly a list (or tuple) of float either a formatted string from which a list of float will be constructed. 346 # The string format is "qmin1:qmax1:qstep1,qmin2:qmax2:qstep2,qmin3:qmax3:qstep3 ..." 347 elif pName == 'qshellvalues': 348 349 try: 350 351 if not isinstance(pValue, (list, tuple, str)): 352 raise 353 354 if isinstance(pValue,(list,tuple)): 355 self.qShellValues = [float(v) for v in pValue] 356 357 elif isinstance(pValue, str): 358 self.qShellValues = [] 359 for v in pValue.split(';'): 360 qMin, qMax, dq = [float(vv) for vv in v.split(':')] 361 362 if qMin < 0: 363 raise 364 365 if qMax <= qMin: 366 raise 367 368 if dq <= 0.0: 369 raise 370 371 if dq > (qMax - qMin): 372 raise 373 374 q0 = qMin 375 while q0 <= qMax + 1.0e-6: 376 self.qShellValues.append(q0) 377 q0 += dq 378 379 self.qShellValues = sorted(list(set(self.qShellValues))) 380 if not self.qShellValues: 381 raise 382 383 except: 384 raise Error('Error when parsing %s parameter".' % pName) 385 386 elif pName == 'qshellwidth': 387 try: 388 self.qShellWidth = float(pValue) 389 if self.qShellWidth <= 0.0: 390 raise 391 392 except: 393 raise Error('Error when parsing %s parameter: must be a float.' % pName) 394 395 # 'qVectorsDirection' input parameter name. 396 elif pName == 'qvectorsdirection': 397 398 if pValue is None: 399 self.qVectorsDirection = None 400 401 elif isinstance(pValue, str): 402 403 if pValue.strip().lower() == 'no': 404 self.qVectorsDirection = None 405 406 else: 407 try: 408 self.qVectorsDirection = [Vector([float(vv.strip()) for vv in v.strip().split(',')]).normal() for v in pValue.split(';')] 409 except: 410 raise Error('Error when parsing %s parameter: each direction must be of the form "x,y,z".' % pName) 411 412 else: 413 raise Error('Error when parsing %s parameter: must be a string.' % pName) 414 415 416 elif pName == 'qvectorsgenerator': 417 418 try: 419 p = pValue.strip().lower() 420 if p not in ['3d isotropic', '2d isotropic', 'anisotropic']: 421 raise 422 423 self.qVectorsGenerator = p 424 425 except: 426 raise Error('Error when parsing %s parameter: must be one of "3d isotropic", "2d isotropic" or "anisotropic".' % pName) 427 428 elif pName == 'qvectorspershell': 429 430 try: 431 self.qVectorsPerShell = int(pValue) 432 if self.qVectorsPerShell <= 0: 433 raise 434 435 except: 436 raise Error('Error when parsing %s parameter: must be an integer.' % pName) 437 438 elif pName == 'referencedirection': 439 if pValue is not None: 440 try: 441 self.referenceDirection = Vector([float(v) for v in pValue.split(',')]).normal() 442 443 except: 444 raise Error('Error when parsing %s parameter: must be "float,float,float".' % pName) 445 446 else: 447 self.referenceDirection = None 448 449 elif pName == 'referenceframe': 450 self.referenceFrame = int(pValue) - 1 451 452 elif pName == 'removetranslation': 453 try: 454 p = pValue.strip().lower() 455 if p == 'yes': 456 self.removeTranslation = True 457 elif p == 'no': 458 self.removeTranslation = False 459 else: 460 raise 461 462 except: 463 raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName) 464 465 elif pName == 'rvalues': 466 467 if not isinstance(pValue, str): 468 raise Error('Error when parsing %s parameter: must be a string.' % pName) 469 470 try: 471 472 self.rMin, self.rMax, self.dr = [float(v) for v in pValue.split(':')] 473 474 if self.rMin < 0.0: 475 raise 476 477 if self.rMax <= self.rMin: 478 raise 479 480 if self.dr > (self.rMax - self.rMin): 481 raise 482 483 r = [] 484 r0 = self.rMin 485 while r0 <= self.rMax + 1.0e-6: 486 r.append(r0) 487 r0 += self.dr 488 489 if not r: 490 raise 491 492 self.nBins = len(r) - 1 493 494 # The rMax is readjusted to fit the last r values of the list. 495 self.rMax = r[-1] 496 497 self.rValues = self.rMin + self.dr/2.0 + self.dr*N.arange(self.nBins) 498 499 except: 500 raise Error('Error when parsing %s parameter: wrong format for the string.' % pName) 501 502 elif pName == 'stepwiserbt': 503 try: 504 p = pValue.strip().lower() 505 if p == 'yes': 506 self.stepwiseRBT = True 507 elif p == 'no': 508 self.stepwiseRBT = False 509 else: 510 raise 511 512 except: 513 raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName) 514 515 elif pName == 'subset': 516 517 if pValue is None: 518 self.subset = 'all' 519 else: 520 self.subset = pValue 521 522 elif pName == 'thetavalues': 523 524 if not isinstance(pValue, str): 525 raise Error('Error when parsing %s parameter: must be a string.' % pName) 526 527 try: 528 529 self.thetaMin, self.thetaMax, self.dTheta = [float(v) for v in pValue.split(':')] 530 531 if self.thetaMin < 0.0: 532 raise 533 534 if self.thetaMin > 180.0: 535 raise 536 537 if self.thetaMax <= self.thetaMin: 538 raise 539 540 if self.dTheta > (self.thetaMax - self.thetaMin): 541 raise 542 543 self.thetaMin *= Units.deg 544 self.thetaMax *= Units.deg 545 self.dTheta *= Units.deg 546 547 theta = [] 548 theta0 = self.thetaMin 549 while theta0 <= self.thetaMax + 1.0e-6: 550 theta.append(theta0) 551 theta0 += self.dTheta 552 553 if not theta: 554 raise 555 556 self.nThetaBins = len(theta) - 1 557 558 # The rMax is readjusted to fit the last r values of the list. 559 self.thetaMax = theta[-1] 560 561 self.thetaValues = self.thetaMin + self.dTheta/2.0 + self.dTheta*N.arange(self.nThetaBins) 562 563 except: 564 raise Error('Error when parsing %s parameter: wrong format for the string.' % pName) 565 566 elif pName == 'timeinfo': 567 568 try: 569 self.first, self.last, self.skip = [int(v) for v in pValue.split(':')] 570 self.first = self.first - 1 571 572 except: 573 raise Error('Error when parsing %s parameter: must be "int:int:int".' % pName) 574 575 elif pName == 'temperature': 576 try: 577 self.temperature = float(pValue) 578 if self.temperature <= 0.0: 579 raise 580 581 except: 582 raise Error('Error when parsing %s parameter: must be a float > 0.' % pName) 583 584 # 'trajectory' input parameter name. 585 586 elif pName == 'trajectory': 587 # If the value is a string, then load the corresponding MMTK object and stores it. 588 if isinstance(pValue, str): 589 # The extension of the trajectory to load. 590 ext = os.path.splitext(pValue)[1] 591 592 # Trajectory set --> loads a TrajectorySet MMTK object. 593 if ext == '.ncs': 594 try: 595 # The trajectory set is opened for reading and its contents sent to |trajSet| list. 596 trajSetFile = open(pValue, 'r') 597 trajSet = [t.strip() for t in trajSetFile.readlines()] 598 trajSetFile.close() 599 600 # Load the whole trajectory set. 601 self.trajectory = TrajectorySet(None, trajSet) 602 603 except: 604 raise Error('Error when parsing %s parameter: unable to open trajectory set file %s.' % pValue) 605 606 # Trajectory --> load a Trajectory. 607 elif ext == '.nc': 608 try: 609 # Load the whole trajectory set. 610 self.trajectory = Trajectory(None, pValue, 'r') 611 612 except IOError: 613 raise Error('Error when parsing %s parameter: unable to open trajectory file %s.' % pValue) 614 615 # Otherwise throw an error. 616 else: 617 raise Error('Error when parsing %s parameter: not a suitable Trajectory or TrajectorySet file name.' % pValue) 618 619 # If the value is already a Trajectory MMTK object, then stores it. 620 elif isinstance(pValue, Trajectory): 621 self.trajectory = pValue 622 623 else: 624 raise Error('Error when parsing %s parameter: neither a trajectory or trajectory set file name neither an instance of Trajectory class.' % pValue) 625 626 self.trajectoryFilename = self.trajectory.filename 627 628 self.universe = self.trajectory.universe 629 630 # The universe contents. 631 setUniverseContents(self.universe) 632 633 elif pName == 'triplet': 634 635 if pValue is None: 636 self.triplet = 'all' 637 else: 638 self.triplet = pValue 639 640 elif pName == 'version': 641 try: 642 if pName == 'unknown': 643 self.version = LooseVersion(vstring = '0.0') 644 else: 645 self.version = LooseVersion(vstring = self.parameters['version']) 646 647 except: 648 raise Error('Error when parsing %s parameter: must be a string of the form "int.int.int".' % pName) 649 650 elif pName == 'weights': 651 self.weights = pValue 652 653 elif pName == 'wignerindexes': 654 655 try: 656 self.wignerIndexes = [int(v) for v in pValue.split(',')] 657 if len(self.wignerIndexes) != 3: 658 raise 659 660 except: 661 raise Error('Error when parsing %s parameter: must be "int:int:int".' % (pName)) 662 663 elif pName in ['ac', 'ados', 'ara', 'avacf', 'cn', 'comt', 'dcsf', 'dcsfar', 'disf', 'disfar',\ 664 'disfg', 'dos', 'eisf', 'gmft', 'isfg', 'msd', 'op', 'opcm', 'pbft', 'pdf', 'qha',\ 665 'rbt','rcf', 'rmsd', 'rog', 'sd', 'sfa', 'scsf', 'sscsf', 'vacf']: 666 667 setattr(self, 'output'+ pName.upper(), pValue)
668
669 - def buildTimeInfo(self):
670 """Builds some attributes related to the frame selection string. They will be used to define at which 671 times a given analysis should be run. 672 """ 673 674 trajLength = len(self.trajectory) 675 676 if self.first < 0: 677 self.first = 0 678 LogMessage('warning', 'The first step must be an integer >= 1.\n\ 679 It will be set to 1 for the running analysis.',['console']) 680 681 if (self.last <= self.first) or (self.last > trajLength): 682 self.last = len(self.trajectory) 683 LogMessage('warning', 'The last step must be an integer in [1,%s].\n\ 684 It will be set to %s for the running analysis.' % (trajLength, trajLength), ['console']) 685 686 if (self.skip <= 0) or (self.skip >= trajLength): 687 self.skip = 1 688 LogMessage('warning', 'The skip step must be an integer in [1,%s].\n\ 689 It will be set to 1 for the running analysis.' % trajLength,['console']) 690 691 self.frameIndexes = N.arange(self.first, self.last, self.skip) 692 if len(self.frameIndexes) == 0: 693 self.first = 1 694 self.last = trajLength 695 self.skip = 1 696 self.frameIndexes = N.arange(self.first,self.last,self.skip) 697 LogMessage('warning','The frame selection is empty.\ 698 It will be set to %s:%s:%s for the running analysis.' % (1,trajLength,1),['console']) 699 700 if hasattr(self.trajectory,'time'): 701 t = self.trajectory.time[self.first:self.last:self.skip] 702 703 else: 704 t = range(trajLength)[self.first:self.last:self.skip] 705 706 if len(t) == 1: 707 self.dt = 1.0 708 else: 709 if t[0] == t[1]: 710 self.dt = 1.0 711 else: 712 self.dt = t[1] - t[0] 713 714 self.times = self.dt * N.arange(len(t)) 715 716 self.nFrames = len(self.times)
717
718 - def preLoadTrajectory(self, structure, differentiation = 1):
719 """ 720 """ 721 722 if '_estimate' in self.__class__.__name__: 723 self.preLoadedTraj = None 724 return 725 726 try: 727 LogMessage('info', 'Trying to preload the trajectory. Please wait ...', ['console']) 728 self.preLoadedTraj = {} 729 if structure == 'atom': 730 test = N.zeros((self.nAtoms, self.nFrames, 3), N.Float) 731 del test 732 for atom in self.subset: 733 if differentiation == 0: 734 self.preLoadedTraj[atom] = self.trajectory.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip, variable = 'velocities').array 735 else: 736 self.preLoadedTraj[atom] = self.trajectory.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array 737 elif structure == 'frame': 738 test = N.zeros((self.nFrames, self.nAtoms, 3), N.Float) 739 del test 740 for frameIndex in range(self.nFrames): 741 frame = self.frameIndexes[frameIndex] 742 self.preLoadedTraj[frameIndex] = self.trajectory.configuration[frame] 743 LogMessage('info', 'Enough memory to preload the trajectory.', ['console']) 744 745 except: 746 self.preLoadedTraj = None 747 LogMessage('info', 'Not enough memory to preload the trajectory. The analysis will run with runtime %s NetCDF extraction.' % structure, ['console'])
748
749 - def saveAnalysis(self, filename):
750 """Saves the settings of an analysis to an output file. 751 752 @param filename: the name of the output file. If the extension is '.nmi' the output file will be a 753 nMOLDYN input script otherwise the output file will be a nMOLDYN autostart script. 754 @type: filename: string 755 """ 756 757 if filename[-4:].lower() == '.nmi': 758 759 # Opens the python script file for writing. 760 pScript = open(filename, 'w') 761 # Write the input file header. 762 pScript.write('################################################################\n') 763 pScript.write('# This is an automatically generated python-nMOLDYN run script #\n') 764 pScript.write('################################################################\n\n') 765 766 # The line that will create an instance of the analysis child class. 767 pScript.write(('\nanalysis = %s\n\n') % (self.__class__.__name__)) 768 769 pScript.write('# General parameters\n') 770 pScript.write(('version = "%s"\n') % nMOLDYN.__version__) 771 pScript.write(('estimate = "%s"\n\n') % self.parameters['estimate']) 772 773 pScript.write('# Analysis-specific parameters\n') 774 # Loop over the individual widget names. 775 for pName in sorted(self.inputParametersNames): 776 pValue = self.parameters[pName] 777 # Special case of strings. Some quotes are added to the value. 778 if isinstance(pValue, str): 779 pScript.write(('%s = "%s"\n') % (pName, pValue)) 780 else: 781 pScript.write(('%s = %s\n') % (pName, pValue)) 782 783 pScript.close() 784 785 else: 786 787 # Opens the python script file for writing. 788 pScript = open(filename, 'w') 789 pScript.write('#!%s\n' % sys.executable) 790 791 # Write the input file header. 792 pScript.write('################################################################\n') 793 pScript.write('# This is an automatically generated python-nMOLDYN run script #\n') 794 pScript.write('################################################################\n\n') 795 # Write the 'import' lines. 796 pScript.write('from MMTK import *\n') 797 pScript.write(('from %s import %s\n\n') % (self.__module__, self.__class__.__name__)) 798 799 # Writes the line that will initialize the parameters dictionnary. 800 pScript.write('parameters = {}\n\n') 801 802 pScript.write('# General parameters\n') 803 pScript.write('parameters[\'version\'] = "%s"\n' % nMOLDYN.__version__) 804 pScript.write('parameters[\'estimate\'] = "%s"\n\n' % self.parameters['estimate']) 805 806 pScript.write('# Analysis-specific parameters\n') 807 # Loop over the individual widget names. 808 for pName in sorted(self.inputParametersNames): 809 pValue = self.parameters[pName] 810 # Special case of strings. Some quotes are added to the value. 811 if isinstance(pValue, str): 812 pScript.write(('parameters[\'%s\'] = "%s"\n') % (pName, pValue)) 813 else: 814 pScript.write(('parameters[\'%s\'] = %s\n') % (pName, pValue)) 815 816 # The line that will create an instance of the analysis child class. 817 pScript.write(('\nanalysis = %s(parameters)') % (self.__class__.__name__)) 818 # The line that will actually launch the analysis. 819 pScript.write('\nanalysis.runAnalysis()') 820 pScript.close() 821 822 # Make the script executable on posix systems. 823 if os.name == 'posix': 824 os.system('chmod u+x %s' % filename)
825
826 - def runAnalysis(self):
827 """Runs an analysis. 828 829 @return: a dictionnary of the form {'days' : d, 'hours' : h, 'minutes' : m, 'seconds' : s} specifying the 830 time the analysis took in dayx, hours, minutes and seconds. 831 @rtype: dict 832 """ 833 834 # The job progress in percentage. 835 self.jobCounter = 0 836 837 # Initialize the analysis. 838 self.initialize() 839 840 # Displays the informations about the running analysis. 841 self.buildJobInfo() 842 843 # Checks the nMOLDYN version of the running analysis 844 if self.version < LooseVersion(vstring = nMOLDYN.__version__): 845 # The elapsed time for the (estimated) analysis is written in the log file. 846 LogMessage('warning','The nMOLDYN version number of the script is unknown or deprecated.\n\ 847 Some analysis keywords name/values may have changed. Please check in case of troubles.',['file','console']) 848 849 # The analysis is actually processed here. 850 self.internalRun() 851 852 # The ending date of the analysis. 853 LogMessage('info','Job finished on '+asctime()+'.\n',['console','file']) 854 855 # The time taken for the analysis. 856 timeTakenForAnalysis = self.analysisTime(self.chrono) 857 858 # The log message that displays the time taken for the analysis is a little bit different 859 # depending on whether we are in estimate mode or not. 860 if self.estimate: 861 # The elapsed time for the (estimated) analysis is written in the log file. 862 LogMessage('info','Estimated time for analysis: %(days)s days %(hours)s hours %(minutes)s minutes %(seconds)s seconds.\n' % timeTakenForAnalysis,['file','console']) 863 else: 864 # The elapsed time for the (estimated) analysis is written in the log file. 865 LogMessage('info','Elapsed time for analysis: %(days)s days %(hours)s hours %(minutes)s minutes %(seconds)s seconds.\n' % timeTakenForAnalysis,['file','console']) 866 867 return timeTakenForAnalysis
868
869 - def updateJobProgress(self, norm):
870 """Check the progress of the running analysis and displays periodically on the console and the logfile how far is the analysis. 871 Called each time a step of an analysis loop is achieved. 872 873 @param norm: the maximum number of steps of the analysis. 874 @type: norm: integer 875 """ 876 877 # computes the old percentage 878 t = int(100*self.jobCounter/norm) 879 880 oldProgress = (t/self.displayProgressRate)*self.displayProgressRate 881 882 # A new step has been done 883 self.jobCounter += 1 884 885 # computes the new percentage 886 t = int(100*self.jobCounter/norm) 887 newProgress = (t/self.displayProgressRate)*self.displayProgressRate 888 889 if newProgress != oldProgress: 890 LogMessage('info','Progress rate = %3d %%' % newProgress,['file','console']) 891 892 if self.statusBar is not None: 893 self.statusBar.setValue(t)
894
895 - def buildJobInfo(self):
896 """Display on the console and in the log file the main ifnormation about the analysis to run. 897 """ 898 899 self.information = "" 900 901 # The ascii time when the job was launched. 902 jobCreationTime = asctime() 903 904 self.information += '#'*90 +'\n' 905 self.information += 'Job information for %s analysis.\n' % self.__class__.__name__ 906 self.information += '#'*90 +'\n\n' 907 908 # Some information about the analysis that will be written as a header of the NetCDF output file. 909 self.information += 'Job launched on: %s\n\n' % jobCreationTime 910 self.information += 'General informations\n' 911 self.information += '--------------------\n' 912 self.information += 'User: %s\n' % getpass.getuser() 913 self.information += 'OS: %s-%s\n' % (platform.system(), platform.release()) 914 self.information += 'Processor: %s\n' % platform.machine() 915 self.information += 'nMOLDYN version: %s\n' % self.parameters['version'] 916 self.information += 'Estimate run: %s\n\n' % self.parameters['estimate'] 917 918 # Some informations for the logfile an output file. 919 self.information += 'Parameters\n' 920 self.information += '----------\n' 921 self.information += '\n'.join(['%s = %s' % (k,self.parameters[k]) for k in self.parameters if k not in ['version','estimate']]) 922 self.information += '\n\n\n' 923 924 # Some informations about the various selections. 925 if hasattr(self, 'subset'): 926 self.information += 'Subset selection\n' 927 self.information += '----------------\n' 928 self.information += 'Number of atoms selected for analysis: %s\n\n' % len(self.subset) 929 930 if hasattr(self, 'deuteration'): 931 self.information += 'Deuteration selection\n' 932 self.information += '---------------------\n' 933 self.information += 'Number of atoms selected for deuteration: %s\n\n' % len(self.deuteration) 934 935 if hasattr(self, 'group'): 936 self.information += 'Group selection\n' 937 self.information += '---------------------\n' 938 self.information += 'Number of groups selected: %s\n\n' % len(self.group) 939 940 if hasattr(self, 'triplet'): 941 self.information += 'Triplet selection\n' 942 self.information += '---------------------\n' 943 self.information += 'Number of triplets selected: %s\n\n' % len(self.triplet) 944 945 if hasattr(self, 'bond'): 946 self.information += 'Bond selection\n' 947 self.information += '---------------------\n' 948 self.information += 'Number of bonds selected: %s\n\n' % len(self.bond) 949 950 self.information += '\n' 951 self.information += 'Job status\n' 952 self.information += '----------\n' 953 954 [LogMessage('info', l, ['file', 'console']) for l in self.information.splitlines()]
955
956 - def analysisTime(self, time):
957 """Converts a time in second in days, hours, minutes and seconds. 958 959 @param time: the time (in seconds) to convert. 960 @type time: integer. 961 962 @return: a dictionnary of the form {'days' : d, 'hours' : h, 'minutes' : m, 'seconds' : s} 963 where d, h, m and s are integers resulting respectively from the conversion of |time| in days, hours, minutes 964 and seconds. 965 @rtype: dict 966 """ 967 968 if time is None: 969 return None 970 971 else: 972 converted = {'days' : 0, 'hours' : 0, 'minutes' : 0, 'seconds' : 0} 973 974 converted['days'], time = divmod(time, 86400) 975 converted['hours'], time = divmod(time, 3600) 976 converted['minutes'], converted['seconds'] = divmod(time, 60) 977 978 return converted
979
980 - def weightingScheme(self, universe, atoms, deuter, scheme = 'equal'):
981 """Returns the weights of |atoms| MMTK collection of |universe| MMTK universe using the 982 weighting scheme |scheme|. 983 984 @param universe: the MMTK universe. 985 @type universe: instance of MMTK.Universe 986 987 @param atoms: the atoms to take into account when defining the weights. 988 @type atoms: instance of MMTK.Collections.Collection 989 990 @param deuter: the hydrogen atoms that will be parametrized as deuterium atoms. 991 @type deuter: instance of MMTK.Collections.Collection 992 993 @param scheme: a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting 994 scheme to use. 995 @type scheme: string 996 997 @return: the weights of the selected atoms. 998 @rtype: an instance of MMTK.ParticleProperties.ParticledScalar 999 """ 1000 1001 try: 1002 if scheme == 'equal': 1003 # MMTK ParticleScalar object containing 1 for each atom of the universe 1004 weights = ParticleScalar(universe, N.ones(universe.numberOfAtoms(), N.Float)) 1005 1006 elif scheme == 'mass': 1007 # The atomic masses of the universe. 1008 weights = universe.masses() 1009 1010 for atom in deuter: 1011 if abs(atom._mass - atom.mass_deut) < 1.e-10: 1012 LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file']) 1013 else: 1014 weights[atom] = atom.mass_deut 1015 1016 elif scheme == 'atomicNumber': 1017 try: 1018 # MMTK ParticleScalar object containing 1 for each atom of the universe 1019 weights = ParticleScalar(universe) 1020 for at in universe.atomList(): 1021 weights[at] = at.atomic_number 1022 except KeyError: 1023 raise Error('One or more element with unknown atomic number in your selection.') 1024 1025 elif scheme == 'coherent': 1026 # getAtomScalarArray is a wrap for getParticleScalar function of Universe.py module 1027 # MMTK ParticleScalar object containing the bcoh. 1028 weights = universe.getAtomScalarArray('b_coherent') 1029 1030 # The b value of the hydrogen to be deuterated are changed to b value of deuterium. 1031 for atom in deuter: 1032 if abs(atom.b_coherent - atom.b_coherent_deut) < 1.e-10: 1033 LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file']) 1034 else: 1035 weights[atom] = atom.b_coherent_deut 1036 1037 elif scheme == 'incoherent': 1038 # getAtomScalarArray is a wrap for getParticleScalar function of Universe.py module 1039 # MMTK ParticleScalar object containing the bincoh. 1040 weights = universe.getAtomScalarArray('b_incoherent') 1041 1042 # The b value of the hydrogen to be deuterated are changed to b value of deuterium. 1043 for atom in deuter: 1044 if abs(atom.b_incoherent - atom.b_incoherent_deut) < 1.e-10: 1045 LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file']) 1046 else: 1047 weights[atom] = atom.b_incoherent_deut 1048 1049 # b^2 1050 weights = weights*weights 1051 1052 else: 1053 raise Error('%s weighting scheme is unknown for %s analysis.' % (scheme, self.__class__.__name__)) 1054 1055 # If there is an atom selection. A boolean mask is applied on the atoms to select. 1056 if atoms != universe: 1057 weights = weights*atoms.booleanMask() 1058 1059 # The normalization factor is different for the coherent weighting scheme. 1060 if scheme == 'coherent': 1061 weights = weights/N.sqrt((weights*weights).sumOverParticles()) 1062 else: 1063 weights = weights/weights.sumOverParticles() 1064 1065 except: 1066 raise Error('Problem when setting up %s weighting scheme for %s analysis.' % (scheme, self.__class__.__name__)) 1067 else: 1068 return weights
1069
1070 - def subsetSelection(self, universe, selection):
1071 """Returns a MMTK collection of atoms that matches |selection| selection string. Used to apply an analysis 1072 to a subset of atoms. 1073 1074 @param universe: the universe on which the selection will be performed. 1075 @type universe: instance of MMTK.Universe 1076 1077 @param selection: the selection string that will define the atoms to select. 1078 @type selection: string 1079 1080 @return: a MMTK Collection of the atoms that matches |selection| selection string. 1081 @rtype: instance of MMTK.Collections.Collection 1082 """ 1083 1084 if not isinstance(selection, (Collection, list, tuple, str)): 1085 raise Error('Wrong format for the subset selection. Must be a string or a MMTK Collection' % str(selection)) 1086 1087 if isinstance(selection, Collection): 1088 return selection 1089 1090 elif isinstance(selection, (list, tuple)): 1091 try: 1092 subsetSelection = Collection(selection) 1093 return subsetSelection 1094 1095 except: 1096 raise Error('The results of a subset selection must be a MMTK Collection.') 1097 1098 # If selection dictionnary is empty or equal to None then, by default, returns the whole universe as the subset. 1099 elif selection.lower() == 'all': 1100 return Collection(universe.atomList()) 1101 1102 else: 1103 1104 # If the selection string contains "filename name-of-file" then it will be a selection from a file. 1105 if selection.lower().strip()[0:8] == 'filename': 1106 subsetSelection = self.__parseFileSelectionString(universe, selection, 'subset') 1107 1108 elif selection.lower().strip()[0:10] == 'expression': 1109 subsetSelection = self.__parseExpressionSelectionString(universe, selection, 'subset') 1110 1111 # Otherwise the selection will be from an object. 1112 else: 1113 subsetSelection = self.__parseObjectSelectionString(universe, selection) 1114 1115 try: 1116 subsetSelection = Collection(list(subsetSelection)) 1117 1118 except: 1119 raise Error('The results of a subset selection must be a MMTK Collection.') 1120 1121 return subsetSelection
1122
1123 - def deuterationSelection(self, universe, selection):
1124 """Returns a MMTK collection of atoms that matches |selection| selection string. Used to switch the parameters 1125 of a subset (or all) of hydrogen atoms to the parameters of deuterium in order to simulate deuterated system. 1126 1127 @param universe: the universe on which the selection will be performed. 1128 @type universe: instance of MMTK.Universe 1129 1130 @param selection: the selection string that will define the atoms to select. 1131 @type selection: string 1132 1133 @return: a MMTK Collection of the atoms that matches |selection| selection string. 1134 @rtype: instance of MMTK.Collections.Collection 1135 """ 1136 1137 if not isinstance(selection, (Collection, list, tuple, str)): 1138 raise Error('Wrong format for the deuteration selection. Must be a string or a MMTK Collection') 1139 1140 if isinstance(selection, Collection): 1141 return selection 1142 1143 elif isinstance(selection, (list, tuple)): 1144 try: 1145 deuterationSelection = Collection(selection) 1146 return deuterationSelection 1147 1148 except: 1149 raise Error('The results of a deuteration selection must be a MMTK Collection.') 1150 1151 elif selection.lower() == 'no': 1152 return Collection() 1153 1154 elif selection.lower() == 'all': 1155 return Collection([at for at in universe.atomList() if at.type.name == 'hydrogen']) 1156 1157 else: 1158 1159 # If the selection string contains "filename name-of-file" then it will be a selection from a file. 1160 if selection.lower().strip()[0:8] == 'filename': 1161 deuterationSelection = self.__parseFileSelectionString(universe, selection, 'deuteration') 1162 1163 elif selection.lower().strip()[0:10] == 'expression': 1164 deuterationSelection = self.__parseExpressionSelectionString(universe, selection, 'deuteration') 1165 1166 # Otherwise the selection will be from an object. 1167 else: 1168 deuterationSelection = self.__parseObjectSelectionString(universe, selection) 1169 1170 deuterationSelection = [at for at in deuterationSelection if at.type.name == 'hydrogen'] 1171 1172 try: 1173 deuterationSelection = Collection(deuterationSelection) 1174 1175 except: 1176 raise Error('The results of a deuteration selection must be a MMTK Collection.') 1177 1178 return deuterationSelection
1179
1180 - def groupSelection(self, universe, selection):
1181 """Returns a list of MMTK collections where each collection defines a group on which will be applied collectively 1182 an analysis. 1183 1184 @param universe: the universe on which the selection will be performed. 1185 @type universe: instance of MMTK.Universe 1186 1187 @param selection: the selection string that will define the contents of each group. 1188 @type selection: string 1189 1190 @return: a list of MMTK Collection where each collection defines a group.. 1191 @rtype: list 1192 """ 1193 1194 if not isinstance(selection, (list, tuple, str)): 1195 raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.') 1196 1197 if isinstance(selection, (list, tuple)): 1198 for el in selection: 1199 if not isinstance(el, Collection): 1200 raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.') 1201 1202 return list(selection) 1203 1204 elif selection.lower() == 'all': 1205 groupSelection = [] 1206 for obj in universe.objectList(): 1207 if isinstance(obj,(PeptideChain, Protein, NucleotideChain)): 1208 for r in obj.residues(): 1209 groupSelection.append(r.atomList()) 1210 1211 else: 1212 groupSelection.append(obj.atomList()) 1213 1214 return [Collection(list(g)) for g in groupSelection] 1215 1216 else: 1217 1218 groupSelection = [] 1219 # If the selection string contains "filename name-of-file" then it will be a selection from a file. 1220 if selection.lower().strip()[0:8] == 'filename': 1221 groupSelection = self.__parseFileSelectionString(universe, selection, 'group') 1222 1223 elif selection.lower().strip()[0:10] == 'expression': 1224 groupSelection = self.__parseExpressionSelectionString(universe, selection, 'group') 1225 1226 # Otherwise the selection will be from an object. 1227 else: 1228 selectionStringLists = [v.strip() for v in re.split('group\d+', selection) if v != ''] 1229 1230 for selString in selectionStringLists: 1231 1232 groupingLevel, sString = [v.strip() for v in re.split('groupinglevel\s+(\w+)\s+', selString) if v != ''] 1233 1234 grp = self.__parseObjectSelectionString(universe, sString) 1235 1236 self.__buildGroup(grp, groupSelection, groupingLevel) 1237 1238 try: 1239 groupSelection = [Collection(list(g)) for g in groupSelection] 1240 1241 except: 1242 raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.') 1243 1244 return groupSelection
1245
1246 - def __buildGroup(self, g, gSelection, gLevel):
1247 """ 1248 """ 1249 1250 gDict= {} 1251 1252 if gLevel.lower() == 'chain': 1253 for at in g: 1254 chain = at.parent.parent.parent 1255 if gDict.has_key(chain): 1256 gDict[chain].append(at) 1257 else: 1258 gDict[chain] = [at] 1259 1260 elif gLevel.lower() == 'residue': 1261 for at in g: 1262 residue = at.parent.parent 1263 if gDict.has_key(residue): 1264 gDict[residue].append(at) 1265 else: 1266 gDict[residue] = [at] 1267 1268 elif gLevel.lower() in ['molecule', 'cluster', 'nucleicacid', 'protein']: 1269 for at in g: 1270 molecule = at.topLevelChemicalObject() 1271 if gDict.has_key(molecule): 1272 gDict[molecule].append(at) 1273 else: 1274 gDict[molecule] = [at] 1275 1276 elif gLevel.lower() == 'atom': 1277 for at in g: 1278 gDict[at] = at 1279 1280 elif gLevel.lower() == 'amine': 1281 for at in g: 1282 amine = belongToAnAmine(at) 1283 if amine is not None: 1284 if gDict.has_key(amine): 1285 gDict[amine].append(at) 1286 else: 1287 gDict[amine] = [at] 1288 1289 elif gLevel.lower() == 'hydroxy': 1290 for at in g: 1291 hydroxy = belongToAHydroxy(at) 1292 if hydroxy is not None: 1293 if gDict.has_key(hydroxy): 1294 gDict[hydroxy].append(at) 1295 else: 1296 gDict[hydroxy] = [at] 1297 1298 elif gLevel.lower() == 'methyl': 1299 for at in g: 1300 methyl = belongToAMethyl(at) 1301 if methyl is not None: 1302 if gDict.has_key(methyl): 1303 gDict[methyl].append(at) 1304 else: 1305 gDict[methyl] = [at] 1306 1307 elif gLevel.lower() == 'thiol': 1308 for at in g: 1309 thiol = belongToAnHydroxy(at) 1310 if thiol is not None: 1311 if gDict.has_key(thiol): 1312 gDict[thiol].append(at) 1313 else: 1314 gDict[thiol] = [at] 1315 1316 elif gLevel.lower() == 'default': 1317 for at in g: 1318 obj = at.topLevelChemicalObject() 1319 1320 if isinstance(obj, Atom): 1321 gDict[at] = at 1322 elif isinstance(obj, (AtomCluster,Molecule)): 1323 if gDict.has_key(obj): 1324 gDict[obj].append(at) 1325 else: 1326 gDict[obj] = [at] 1327 elif isinstance(obj, (NucleotideChain,PeptideChain,Protein)): 1328 residue = at.parent.parent 1329 if gDict.has_key(residue): 1330 gDict[residue].append(at) 1331 else: 1332 gDict[residue] = [at] 1333 1334 for v in gDict.values(): 1335 gSelection.append(set(v))
1336
1337 - def __parseExpressionSelectionString(self, universe, selectionString, selectionMode):
1338 """Parses a nMOLDYN selection string that starts with 'expression' keyword. 1339 1340 @param universe: the universe on which the selection will be performed. 1341 @type universe: instance of MMTK.Universe 1342 1343 @param selectionString: the selection string to parse. 1344 @type selectionString: string 1345 1346 @param selectionMode: a string being one of 'subset', 'deuteration' or 'group' that will specify which kind of 1347 selection should be performed. 1348 @type selectionMode: string 1349 1350 @return: will depend on |selectionMode|. 1351 @rtype: Python set for |selectionMode| = 'subset' or 'deuteration' and dict for |selectionMode| = 'group' 1352 @return: a set of the atoms that match |selectionString| selection string. 1353 @rtype: set 1354 """ 1355 1356 try: 1357 expression = ' '.join(selectionString.split()[1:]) 1358 exec(expression) 1359 except: 1360 raise Error('Could not parse the expression of the selection string.') 1361 1362 if not isinstance(selection,list): 1363 raise Error('The expression of the selection string must provide a list.') 1364 1365 if selectionMode in ['subset', 'deuteration']: 1366 1367 selection = set(selection) 1368 1369 elif selectionMode == 'group': 1370 pass 1371 1372 if not selection: 1373 raise Error('The selection associated with\n%s\nselection string gave an empty selection.\nSomething might be wrong with it.' % selectionString) 1374 1375 return selection
1376
1377 - def __parseFileSelectionString(self, universe, selectionString, selectionMode):
1378 """Parses a nMOLDYN Selection file (nms file) in order to perform a subset, a deuteration or a group 1379 selection. 1380 1381 @param universe: the universe on which the selection will be performed. 1382 @type universe: instance of MMTK.Universe 1383 1384 @param selectionString: the selection string to parse. 1385 @type selectionString: string 1386 1387 @param selectionMode: a string being one of 'subset', 'deuteration' or 'group' that will specify which kind of 1388 selection should be performed. 1389 @type selectionMode: string 1390 1391 @return: will depend on |selectionMode|. 1392 @rtype: Python set for |selectionMode| = 'subset' or 'deuteration' and dict for |selectionMode| = 'group' 1393 """ 1394 1395 # First try to open the nms file that contains the informations about the atom to select. 1396 try: 1397 filename = selectionString.split()[1].strip() 1398 nmsFile = open(filename, 'r') 1399 exec nmsFile 1400 nmsFile.close() 1401 except: 1402 raise Error('Unable to open the file %s.' % filename) 1403 1404 try: 1405 pdb = pdb.strip() 1406 if not os.path.exists(pdb): 1407 raise Error('The pdb file %s does not exist.' % pdb) 1408 1409 except NameError: 1410 raise Error('The "pdb" field was not set in the %s nms file.' % filename) 1411 1412 # Sets up the output variable depending on the selected selection type. 1413 if selectionMode in ['subset', 'deuteration']: 1414 try: 1415 selectedAtoms = [int(v) for v in eval(selectionMode)] 1416 except: 1417 raise Error('The %s selection field must be of the form %s = [int,int,int...]' % (selectionMode,selectionMode)) 1418 1419 selection = set() 1420 1421 elif selectionMode == 'group': 1422 try: 1423 selectedAtoms = [[int(vv) for vv in v] for v in eval(selectionMode)] 1424 1425 except : 1426 raise Error('The group selection field must be of the form group = [[int,int,...],[int,int,...],[int,int,...]...]') 1427 1428 selection = [set()]*len(selectedAtoms) 1429 1430 else: 1431 raise Error('Unknown selection type: %s.' % selectionMode) 1432 1433 try: 1434 atomPos = [] 1435 # Loop over all the atom of the universe. 1436 for atom in universe.atomList(): 1437 atomPos.append([atom, [round(v,2) for v in self.trajectory.configuration[0][atom]]]) 1438 1439 # The PDB file is opened for reading. 1440 pdb = PDBConfiguration(pdb) 1441 1442 if selectionMode in ['subset','deuteration']: 1443 # Loop over all the residues found in the pdb file. Here residues are defined in a broader term. 1444 # If the system is not a protein, residues referes to a single PDB line. 1445 for res in pdb.residues: 1446 # Loop over all the atoms of |res|. 1447 for at in res.atom_list: 1448 # The atom should be selected. 1449 if at.properties['serial_number'] in selectedAtoms: 1450 pos = [round(v,2) for v in at.position] 1451 # Loop over all the atom of the universe. 1452 for atom in atomPos: 1453 # Find the one that matches the coordinates of atom |at|. 1454 if atom[1] == pos: 1455 # And add it to the selection. 1456 selection.add(atom[0]) 1457 break 1458 1459 elif selectionMode == 'group': 1460 # Loop over all the residues found in the pdb file. Here residues are definer in a broader term. 1461 # If the system is not a protein, residues referes to a single PDB line. 1462 for res in pdb.residues: 1463 # Loop over all the atoms of |res|. 1464 for at in res.atom_list: 1465 # Loop over all the selected atom pairs tuples. 1466 for p in selectedAtoms: 1467 # The pdb atom should be selected. 1468 if at.properties['serial_number'] in p: 1469 pos = [round(v,2) for v in at.position] 1470 # Loop over all the atom of the universe. 1471 for atom in atomPos: 1472 # Find the one that matches the coordinates of atom |at|. 1473 if atom[1] == pos: 1474 selection[selectedAtoms.index(p)].add(atom[0]) 1475 break 1476 1477 except: 1478 raise Error('Could not parse properly for selection the PDB file associated with %s selection file.' % filename) 1479 1480 if not selection: 1481 raise Error('The selection from file %s gave an empty selection. Something might be wrong with it.' % filename) 1482 1483 return selection
1484
1485 - def __parseObjectSelectionString(self, universe, selectionString):
1486 """Parses a selection string. 1487 1488 @param universe: the universe on which the selection will be performed. 1489 @type universe: instance of MMTK.Universe 1490 1491 @param selectionString: the selection strin to parse. 1492 @type selectionString: string 1493 1494 @return: a set of the atoms that match |selectionString| selection string. 1495 @rtype: set 1496 """ 1497 1498 selStr = re.sub('\(',' ( ',selectionString) 1499 selStr = re.sub('\)',' ) ',selStr) 1500 1501 l = [v.strip() for v in re.split('\s+|,',selStr) if v.strip()] 1502 1503 processedList = [] 1504 selectionValue = [] 1505 1506 try: 1507 while True: 1508 1509 if not l: 1510 break 1511 1512 l0 = l[0].lower() 1513 if l0 == 'objectname': 1514 objectName = l[1] 1515 del l[0:2] 1516 objectClass = universe.nmoldyncontents[objectName]['objectclass'].lower() 1517 continue 1518 1519 elif l0 in universe.nmoldyncontents[objectName].keys(): 1520 selectionKeyword = l0 1521 del l[0] 1522 1523 while True: 1524 1525 ll0 = l[0].lower() 1526 if ll0 == 'objectname': 1527 raise 1528 1529 elif ll0 in universe.nmoldyncontents[objectName].keys(): 1530 raise 1531 1532 elif ll0 in ['and','or',')']: 1533 if not selectionValue: 1534 raise 1535 indexes = self.__retrieveAtomIndexes(universe, objectClass, objectName, selectionKeyword, selectionValue) 1536 processedList.append(str(set(indexes))) 1537 selectionValue = [] 1538 break 1539 1540 elif ll0 == '(': 1541 raise 1542 1543 else: 1544 selectionValue.append(ll0) 1545 del l[0] 1546 if not l: 1547 indexes = self.__retrieveAtomIndexes(universe, objectClass, objectName, selectionKeyword, selectionValue) 1548 processedList.append(str(set(indexes))) 1549 selectionValue = [] 1550 break 1551 1552 elif l0 in ['and','or','(']: 1553 processedList.append(l0) 1554 del l[0] 1555 if not l: 1556 raise 1557 1558 elif l0 == ')': 1559 processedList.append(l0) 1560 del l[0] 1561 1562 else: 1563 raise 1564 1565 except: 1566 raise Error('Badly formatted selection string.') 1567 1568 processedString = ''.join(processedList) 1569 processedString = processedString.replace('and', '&') 1570 processedString = processedString.replace('or', '|') 1571 1572 # The MMTK objects corresponding to the selection indexes are returned. 1573 selection = set([universe.indexToAtoms[ind] for ind in eval(processedString)]) 1574 1575 if not selection: 1576 raise Error('The selection associated with\n%s\nselection string gave an empty selection.\nSomething might be wrong with it.' % selectionString) 1577 1578 return selection
1579
1580 - def __retrieveAtomIndexes(self, universe, objectClass, objectName, selectionKeyword, selectionValue):
1581 """Retrieves the MMTK index of the atoms matching |objectClass| MMTK chemical object class, 1582 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 1583 1584 @param universe: the universe on which the selection will be performed. 1585 @type universe: instance of MMTK.Universe 1586 1587 @param objectClass: the MMTK chemical object to match. 1588 @type objectClass: string 1589 1590 @param objectName: the nMOLDYN name to match. 1591 @type objectName: string 1592 1593 @param selectionKeyword: the selection keyword to match. 1594 @type selectionKeyword: string 1595 1596 @param selectionValue: the selection value to match. 1597 @type selectionValue: string 1598 1599 @return: a list of MMTK indexes of the selected atoms. 1600 @rtype: list 1601 """ 1602 1603 if objectClass == 'allclass': 1604 # The atom parser. 1605 indexes = self.__allClassParser(universe, selectionKeyword, selectionValue) 1606 1607 elif objectClass == 'atom': 1608 # The atom parser. 1609 indexes = self.__atomParser(universe, objectName.lower(), selectionKeyword, selectionValue) 1610 1611 elif objectClass == 'atomcluster': 1612 # The atom cluster parser. 1613 indexes = self.__atomClusterParser(universe, objectName.lower(), selectionKeyword, selectionValue) 1614 1615 elif objectClass == 'molecule': 1616 # The molecule parser. 1617 indexes = self.__moleculeParser(universe, objectName.lower(), selectionKeyword, selectionValue) 1618 1619 elif objectClass == 'nucleotidechain': 1620 # The nucleotide chain parser. 1621 indexes = self.__nucleotideChainParser(universe, objectName.lower(), selectionKeyword, selectionValue) 1622 1623 elif objectClass == 'peptidechain': 1624 # The peptide chain parser. 1625 indexes = self.__peptideChainParser(universe, objectName.lower(), selectionKeyword, selectionValue) 1626 1627 elif objectClass == 'protein': 1628 # The protein parser. 1629 indexes = self.__proteinParser(universe, objectName.lower(), selectionKeyword, selectionValue) 1630 1631 else: 1632 raise Error('The chemical class %s is not defined in MMTK.' % objectClass) 1633 1634 return indexes
1635
1636 - def __allClassParser(self, universe, selectionKeyword, selectionValue):
1637 """Retrieves the MMTK indexes of the atoms whose nMOLDYN name is '*' matching |selectionKeyword| 1638 selection keyword and |selectionValue| selection value. 1639 1640 @param universe: the universe on which the selection will be performed. 1641 @type universe: instance of MMTK.Universe 1642 1643 @param selectionKeyword: the selection keyword to match. 1644 @type selectionKeyword: string 1645 1646 @param selectionValue: the selection value to match. 1647 @type selectionValue: string 1648 1649 @return: a list of MMTK indexes of the selected atoms. 1650 @rtype: list 1651 """ 1652 1653 selectedObjects = universe.objectList() 1654 1655 selection = [] 1656 1657 if selectionKeyword == 'atomelement': 1658 if '*' in selectionValue: 1659 for obj in selectedObjects: 1660 selection.extend([at.index for at in obj.atomList()]) 1661 else: 1662 selection = [] 1663 for v in selectionValue: 1664 for obj in selectedObjects: 1665 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v]) 1666 1667 return selection
1668
1669 - def __atomParser(self, universe, objectName, selectionKeyword, selectionValue):
1670 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Atom' matching 1671 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 1672 1673 @param universe: the universe on which the selection will be performed. 1674 @type universe: instance of MMTK.Universe 1675 1676 @param objectName: the nMOLDYN name to match. 1677 @type objectName: string 1678 1679 @param selectionKeyword: the selection keyword to match. 1680 @type selectionKeyword: string 1681 1682 @param selectionValue: the selection value to match. 1683 @type selectionValue: string 1684 1685 @return: a list of MMTK indexes of the selected atoms. 1686 @rtype: list 1687 """ 1688 1689 selection = [] 1690 1691 if selectionKeyword == 'name': 1692 # If the '*' selection keyword is in the 'name'-key values then the selection is made of all the objects whose 1693 # nmoldynname is 'selectionKey'. 1694 if '*' in selectionValue: 1695 selection = [obj.index for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName] 1696 1697 # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword 1698 # processed. 1699 else: 1700 selection = [] 1701 for v in selectionValue: 1702 selection.extend([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v]) 1703 1704 return selection
1705
1706 - def __atomClusterParser(self, universe, objectName, selectionKeyword, selectionValue):
1707 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'AtomCluster' matching 1708 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 1709 1710 @param universe: the universe on which the selection will be performed. 1711 @type universe: instance of MMTK.Universe 1712 1713 @param objectName: the nMOLDYN name to match. 1714 @type objectName: string 1715 1716 @param selectionKeyword: the selection keyword to match. 1717 @type selectionKeyword: string 1718 1719 @param selectionValue: the selection value to match. 1720 @type selectionValue: string 1721 1722 @return: a list of MMTK indexes of the selected atoms. 1723 @rtype: list 1724 """ 1725 1726 selection = [] 1727 1728 # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects 1729 # that will match the selection 'name'-key corresponding values. 1730 if selectionKeyword == 'name': 1731 # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose 1732 # nmoldynname is 'selectionKey'. 1733 if '*' in selectionValue: 1734 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 1735 1736 # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword 1737 # processed. 1738 else: 1739 selectedObjects = set() 1740 for v in selectionValue: 1741 selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v]) 1742 1743 for obj in selectedObjects: 1744 selection.extend([at.index for at in obj.atomList()]) 1745 1746 else: 1747 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 1748 1749 if selectionKeyword == 'atomname': 1750 for v in selectionValue: 1751 if v == '*': 1752 for obj in selectedObjects: 1753 selection.extend([at.index for at in obj.atomList()]) 1754 else: 1755 for obj in selectedObjects: 1756 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v]) 1757 1758 elif selectionKeyword == 'atomelement': 1759 for v in selectionValue: 1760 if v == '*': 1761 for obj in selectedObjects: 1762 selection.extend([at.index for at in obj.atomList()]) 1763 else: 1764 for obj in selectedObjects: 1765 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v]) 1766 1767 return selection
1768
1769 - def __moleculeParser(self, universe, objectName, selectionKeyword, selectionValue):
1770 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Molecule' matching 1771 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 1772 1773 @param universe: the universe on which the selection will be performed. 1774 @type universe: instance of MMTK.Universe 1775 1776 @param objectName: the nMOLDYN name to match. 1777 @type objectName: string 1778 1779 @param selectionKeyword: the selection keyword to match. 1780 @type selectionKeyword: string 1781 1782 @param selectionValue: the selection value to match. 1783 @type selectionValue: string 1784 1785 @return: a list of MMTK indexes of the selected atoms. 1786 @rtype: list 1787 """ 1788 1789 selection = [] 1790 # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects 1791 # that will match the selection 'name'-key corresponding values. 1792 if selectionKeyword == 'name': 1793 # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose 1794 # nmoldynname is 'selectionKey'. 1795 if '*' in selectionValue: 1796 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 1797 # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword 1798 # processed. 1799 1800 else: 1801 selectedObjects = set() 1802 for v in selectionValue: 1803 selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v]) 1804 for obj in selectedObjects: 1805 selection.extend([at.index for at in obj.atomList()]) 1806 1807 else: 1808 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 1809 1810 if selectionKeyword == 'atomname': 1811 for v in selectionValue: 1812 if v == '*': 1813 for obj in selectedObjects: 1814 selection.extend([at.index for at in obj.atomList()]) 1815 else: 1816 for obj in selectedObjects: 1817 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v]) 1818 1819 elif selectionKeyword == 'atomelement': 1820 for v in selectionValue: 1821 if v == '*': 1822 for obj in selectedObjects: 1823 selection.extend([at.index for at in obj.atomList()]) 1824 else: 1825 for obj in selectedObjects: 1826 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v]) 1827 1828 elif selectionKeyword == 'chemfragment': 1829 for v in selectionValue: 1830 1831 if v == 'amine': 1832 for obj in selectedObjects: 1833 nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen'] 1834 for n in nitrogens: 1835 neighbours = n.bondedTo() 1836 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 1837 # The amine. 1838 if len(hydrogens) == 2: 1839 selection.extend([n.index] + [h.index for h in hydrogens]) 1840 # The ammonium. 1841 elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4): 1842 selection.extend([n.index] + [h.index for h in hydrogens]) 1843 1844 elif v == 'hydroxy': 1845 for obj in selectedObjects: 1846 oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen'] 1847 for o in oxygens: 1848 neighbours = o.bondedTo() 1849 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 1850 # The hydroxy. 1851 if len(hydrogens) == 1: 1852 selection.extend([o.index] + [h.index for h in hydrogens]) 1853 # The water. 1854 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3): 1855 selection.extend([o.index] + [h.index for h in hydrogens]) 1856 1857 elif v == 'methyl': 1858 for obj in selectedObjects: 1859 carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon'] 1860 for c in carbons: 1861 neighbours = c.bondedTo() 1862 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 1863 # The methyl 1864 if len(hydrogens) == 3: 1865 selection.extend([c.index] + [h.index for h in hydrogens]) 1866 # The methane 1867 elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5): 1868 selection.extend([c.index] + [h.index for h in hydrogens]) 1869 1870 elif v == 'thiol': 1871 for obj in selectedObjects: 1872 sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur','sulfur']] 1873 for s in sulphurs: 1874 neighbours = s.bondedTo() 1875 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 1876 # The thiol. 1877 if len(hydrogens) == 1: 1878 selection.extend([o.index] + [h.index for h in hydrogens]) 1879 # The SH2. Quite unusual ... 1880 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3): 1881 selection.extend([s.index] + [h.index for h in hydrogens]) 1882 1883 return selection
1884
1885 - def __nucleotideChainParser(self, universe, objectName, selectionKeyword, selectionValue):
1886 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'NucleotideChain' matching 1887 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 1888 1889 @param universe: the universe on which the selection will be performed. 1890 @type universe: instance of MMTK.Universe 1891 1892 @param objectName: the nMOLDYN name to match. 1893 @type objectName: string 1894 1895 @param selectionKeyword: the selection keyword to match. 1896 @type selectionKeyword: string 1897 1898 @param selectionValue: the selection value to match. 1899 @type selectionValue: string 1900 1901 @return: a list of MMTK indexes of the selected atoms. 1902 @rtype: list 1903 """ 1904 1905 selection = [] 1906 1907 # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects 1908 # that will match the selection 'name'-key corresponding values. 1909 if selectionKeyword == 'name': 1910 # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose 1911 # nmoldynname is 'selectionKey'. 1912 if '*' in selectionValue: 1913 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 1914 1915 # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword 1916 # processed. 1917 else: 1918 selectedObjects = set([]) 1919 for v in selectionValue: 1920 selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v]) 1921 1922 for obj in selectedObjects: 1923 selection.extend([at.index for at in obj.atomList()]) 1924 1925 else: 1926 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 1927 1928 if selectionKeyword == 'atomname': 1929 for v in selectionValue: 1930 if v == '*': 1931 for obj in selectedObjects: 1932 selection.extend([at.index for at in obj.atomList()]) 1933 else: 1934 for obj in selectedObjects: 1935 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v]) 1936 1937 elif selectionKeyword == 'atomtype': 1938 for v in selectionValue: 1939 if v == '*': 1940 for obj in selectedObjects: 1941 selection.extend([at.index for at in obj.atomList()]) 1942 1943 else: 1944 for obj in selectedObjects: 1945 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v]) 1946 1947 elif selectionKeyword == 'atomelement': 1948 for v in selectionValue: 1949 if v == '*': 1950 for obj in selectedObjects: 1951 selection.extend([at.index for at in obj.atomList()]) 1952 else: 1953 for obj in selectedObjects: 1954 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v]) 1955 1956 elif selectionKeyword == 'nuclname': 1957 for v in selectionValue: 1958 if v == '*': 1959 for obj in selectedObjects: 1960 selection.extend([at.index for at in obj.atomList()]) 1961 else: 1962 for obj in selectedObjects: 1963 for nucl in obj.residues(): 1964 if nucl.fullName().strip().lower() == v: 1965 selection.extend([at.index for at in nucl.atomList()]) 1966 1967 elif selectionKeyword == 'nucltype': 1968 for v in selectionValue: 1969 if v == '*': 1970 for obj in selectedObjects: 1971 selection.extend([at.index for at in obj.atomList()]) 1972 else: 1973 for obj in selectedObjects: 1974 for nucl in obj.residues(): 1975 if nucl.symbol.strip().lower() == v: 1976 selection.extend([at.index for at in nucl.atomList()]) 1977 1978 elif selectionKeyword == 'misc': 1979 for v in selectionValue: 1980 if v == 'bases': 1981 for obj in selectedObjects: 1982 selection.extend([at.index for at in obj.bases().atomList()]) 1983 1984 elif v == 'backbone': 1985 for obj in selectedObjects: 1986 selection.extend([at.index for at in obj.backbone().atomList()]) 1987 1988 return selection
1989
1990 - def __peptideChainParser(self, universe, objectName, selectionKeyword, selectionValue):
1991 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'PeptideChain' matching 1992 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 1993 1994 @param universe: the universe on which the selection will be performed. 1995 @type universe: instance of MMTK.Universe 1996 1997 @param objectName: the nMOLDYN name to match. 1998 @type objectName: string 1999 2000 @param selectionKeyword: the selection keyword to match. 2001 @type selectionKeyword: string 2002 2003 @param selectionValue: the selection value to match. 2004 @type selectionValue: string 2005 2006 @return: a list of MMTK indexes of the selected atoms. 2007 @rtype: list 2008 """ 2009 2010 selection = [] 2011 # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects 2012 # that will match the selection 'name'-key corresponding values. 2013 if selectionKeyword == 'name': 2014 2015 # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose 2016 # nmoldynname is 'selectionKey'. 2017 if '*' in selectionValue: 2018 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 2019 2020 # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword 2021 # processed. 2022 else: 2023 selectedObjects = set([]) 2024 for v in selectionValue: 2025 selectedObjects.update([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v]) 2026 2027 for obj in selectedObjects: 2028 selection.extend([at.index for at in obj.atomList()]) 2029 else: 2030 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 2031 2032 if selectionKeyword == 'atomname': 2033 for v in selectionValue: 2034 if v == '*': 2035 for obj in selectedObjects: 2036 selection.extend([at.index for at in obj.atomList()]) 2037 else: 2038 for obj in selectedObjects: 2039 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v]) 2040 2041 elif selectionKeyword == 'atomtype': 2042 for v in selectionValue: 2043 if v == '*': 2044 for obj in selectedObjects: 2045 selection.extend([at.index for at in obj.atomList()]) 2046 else: 2047 for obj in selectedObjects: 2048 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v]) 2049 2050 elif selectionKeyword == 'atomelement': 2051 for v in selectionValue: 2052 if v == '*': 2053 for obj in selectedObjects: 2054 selection.extend([at.index for at in obj.atomList()]) 2055 else: 2056 for obj in selectedObjects: 2057 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v]) 2058 2059 elif selectionKeyword == 'chemfragment': 2060 for v in selectionValue: 2061 if v == 'amine': 2062 for obj in selectedObjects: 2063 nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen'] 2064 for n in nitrogens: 2065 neighbours = n.bondedTo() 2066 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2067 # The amine. 2068 if len(hydrogens) == 2: 2069 selection.extend([n.index] + [h.index for h in hydrogens]) 2070 # The ammonium. 2071 elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4): 2072 selection.extend([n.index] + [h.index for h in hydrogens]) 2073 2074 elif v == 'c_alphas': 2075 for obj in selectedObjects: 2076 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == 'c_alpha']) 2077 2078 elif v == 'hydroxy': 2079 for obj in selectedObjects: 2080 oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen'] 2081 for o in oxygens: 2082 neighbours = o.bondedTo() 2083 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2084 # The hydroxy. 2085 if len(hydrogens) == 1: 2086 selection.extend([o.index] + [h.index for h in hydrogens]) 2087 # The water. 2088 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3): 2089 selection.extend([o.index] + [h.index for h in hydrogens]) 2090 2091 elif v == 'methyl': 2092 for obj in selectedObjects: 2093 carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon'] 2094 for c in carbons: 2095 neighbours = c.bondedTo() 2096 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2097 # The methyl 2098 if len(hydrogens) == 3: 2099 selection.extend([c.index] + [h.index for h in hydrogens]) 2100 # The methane 2101 elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5): 2102 selection.extend([c.index] + [h.index for h in hydrogens]) 2103 2104 elif v == 'thiol': 2105 for obj in selectedObjects: 2106 sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur', 'sulfur']] 2107 for s in sulphurs: 2108 neighbours = s.bondedTo() 2109 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2110 # The thiol. 2111 if len(hydrogens) == 1: 2112 selection.extend([o.index] + [h.index for h in hydrogens]) 2113 # The SH2. Quite unusual ... 2114 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3): 2115 selection.extend([s.index] + [h.index for h in hydrogens]) 2116 2117 elif selectionKeyword == 'resname': 2118 for v in selectionValue: 2119 if v == '*': 2120 for obj in selectedObjects: 2121 selection.extend([at.index for at in obj.atomList()]) 2122 else: 2123 for obj in selectedObjects: 2124 for res in obj.residues(): 2125 if res.fullName().strip().lower() == v: 2126 selection.extend([at.index for at in res.atomList()]) 2127 2128 elif selectionKeyword == 'restype': 2129 for v in selectionValue: 2130 if v == '*': 2131 for obj in selectedObjects: 2132 selection.extend([at.index for at in obj.atomList()]) 2133 else: 2134 for obj in selectedObjects: 2135 for res in obj.residues(): 2136 if res.symbol.strip().lower() == v: 2137 selection.extend([at.index for at in res.atomList()]) 2138 2139 elif selectionKeyword == 'resclass': 2140 for v in selectionValue: 2141 for obj in selectedObjects: 2142 for res in obj.residues(): 2143 if res.symbol.strip().lower() in residusChemFamily[v]: 2144 selection.extend([at.index for at in res.atomList()]) 2145 2146 elif selectionKeyword == 'misc': 2147 for v in selectionValue: 2148 if v == 'sidechains': 2149 for obj in selectedObjects: 2150 selection.extend([at.index for at in obj.sidechains().atomList()]) 2151 2152 elif v == 'backbone': 2153 for obj in selectedObjects: 2154 selection.extend([at.index for at in obj.backbone().atomList()]) 2155 2156 return selection
2157
2158 - def __proteinParser(self, universe, objectName, selectionKeyword, selectionValue):
2159 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Protein' matching 2160 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value. 2161 2162 @param universe: the universe on which the selection will be performed. 2163 @type universe: instance of MMTK.Universe 2164 2165 @param objectName: the nMOLDYN name to match. 2166 @type objectName: string 2167 2168 @param selectionKeyword: the selection keyword to match. 2169 @type selectionKeyword: string 2170 2171 @param selectionValue: the selection value to match. 2172 @type selectionValue: string 2173 2174 @return: a list of MMTK indexes of the selected atoms. 2175 @rtype: list 2176 """ 2177 2178 selection = [] 2179 # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects 2180 # that will match the selection 'name'-key corresponding values. 2181 if selectionKeyword == 'name': 2182 # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose 2183 # nmoldynname is 'selectionKey'. 2184 if '*' in selectionValue: 2185 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 2186 # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword 2187 # processed. 2188 else: 2189 selectedObjects = set([]) 2190 for v in selectionValue: 2191 selectedObjects.update([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v]) 2192 2193 for obj in selectedObjects: 2194 selection.extend([at.index for at in obj.atomList()]) 2195 else: 2196 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]) 2197 2198 2199 if selectionKeyword == 'atomname': 2200 for v in selectionValue: 2201 if v == '*': 2202 for obj in selectedObjects: 2203 selection.extend([at.index for at in obj.atomList()]) 2204 else: 2205 for obj in selectedObjects: 2206 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v]) 2207 2208 elif selectionKeyword == 'atomtype': 2209 for v in selectionValue: 2210 if v == '*': 2211 for obj in selectedObjects: 2212 selection.extend([at.index for at in obj.atomList()]) 2213 else: 2214 for obj in selectedObjects: 2215 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v]) 2216 2217 elif selectionKeyword == 'atomelement': 2218 for v in selectionValue: 2219 if v == '*': 2220 for obj in selectedObjects: 2221 selection.extend([at.index for at in obj.atomList()]) 2222 else: 2223 for obj in selectedObjects: 2224 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v]) 2225 2226 elif selectionKeyword == 'chemfragment': 2227 for v in selectionValue: 2228 2229 if v == 'amine': 2230 for obj in selectedObjects: 2231 nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen'] 2232 for n in nitrogens: 2233 neighbours = n.bondedTo() 2234 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2235 # The amine. 2236 if len(hydrogens) == 2: 2237 selection.extend([n.index] + [h.index for h in hydrogens]) 2238 # The ammonium. 2239 elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4): 2240 selection.extend([n.index] + [h.index for h in hydrogens]) 2241 2242 elif v == 'c_alphas': 2243 for obj in selectedObjects: 2244 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == 'c_alpha']) 2245 2246 elif v == 'hydroxy': 2247 for obj in selectedObjects: 2248 oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen'] 2249 for o in oxygens: 2250 neighbours = o.bondedTo() 2251 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2252 # The hydroxy. 2253 if len(hydrogens) == 1: 2254 selection.extend([o.index] + [h.index for h in hydrogens]) 2255 # The water. 2256 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3): 2257 selection.extend([o.index] + [h.index for h in hydrogens]) 2258 2259 elif v == 'methyl': 2260 for obj in selectedObjects: 2261 carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon'] 2262 for c in carbons: 2263 neighbours = c.bondedTo() 2264 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2265 # The methyl 2266 if len(hydrogens) == 3: 2267 selection.extend([c.index] + [h.index for h in hydrogens]) 2268 # The methane 2269 elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5): 2270 selection.extend([c.index] + [h.index for h in hydrogens]) 2271 2272 elif v == 'thiol': 2273 for obj in selectedObjects: 2274 sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur','sulfur']] 2275 for s in sulphurs: 2276 neighbours = s.bondedTo() 2277 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen'] 2278 # The thiol. 2279 if len(hydrogens) == 1: 2280 selection.extend([o.index] + [h.index for h in hydrogens]) 2281 # The SH2. Quite unusual ... 2282 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3): 2283 selection.extend([s.index] + [h.index for h in hydrogens]) 2284 2285 elif selectionKeyword == 'resname': 2286 for v in selectionValue: 2287 if v == '*': 2288 for obj in selectedObjects: 2289 selection.extend([at.index for at in obj.atomList()]) 2290 else: 2291 for obj in selectedObjects: 2292 for res in obj.residues(): 2293 if res.fullName().strip().lower() == v: 2294 selection.extend([at.index for at in res.atomList()]) 2295 2296 elif selectionKeyword == 'restype': 2297 for v in selectionValue: 2298 if v == '*': 2299 for obj in selectedObjects: 2300 selection.extend([at.index for at in obj.atomList()]) 2301 else: 2302 for obj in selectedObjects: 2303 for res in obj.residues(): 2304 if res.symbol.strip().lower() == v: 2305 selection.extend([at.index for at in res.atomList()]) 2306 2307 elif selectionKeyword == 'resclass': 2308 for v in selectionValue: 2309 for obj in selectedObjects: 2310 for res in obj.residues(): 2311 if res.symbol.strip().lower() in residusChemFamily[v]: 2312 selection.extend([at.index for at in res.atomList()]) 2313 2314 elif selectionKeyword == 'chainname': 2315 for v in selectionValue: 2316 if v == '*': 2317 for obj in selectedObjects: 2318 selection.extend([at.index for at in obj.atomList()]) 2319 else: 2320 for obj in selectedObjects: 2321 for chain in obj: 2322 if chain.name.strip().lower() == v: 2323 selection.extend([at.index for at in chain.atomList()]) 2324 2325 elif selectionKeyword == 'misc': 2326 for v in selectionValue: 2327 if v == 'sidechains': 2328 for obj in selectedObjects: 2329 selection.extend([at.index for at in obj.sidechains().atomList()]) 2330 2331 elif v == 'backbone': 2332 for obj in selectedObjects: 2333 selection.extend([at.index for at in obj.backbone().atomList()]) 2334 2335 return selection
2336
2337 -def setUniverseContents(universe):
2338 """Sets the contents of each object found in the universe. 2339 2340 @param universe: the MMTK universe to look in. 2341 @type universe: a instance of MMTK.Universe. 2342 """ 2343 2344 if hasattr(universe, 'nmoldyncontents'): 2345 return 2346 2347 isotope = {} 2348 universe.indexToAtoms = {} 2349 # Loop over the object found in the universe. 2350 for obj in universe.objectList(): 2351 if not isinstance(obj,(Atom, AtomCluster, Molecule, NucleotideChain, PeptideChain, Protein)): 2352 continue 2353 2354 for at in obj.atomList(): 2355 try: 2356 universe.indexToAtoms[at.index] = at 2357 except: 2358 pass 2359 2360 try: 2361 # Check if the object is referenced in the MMTK database. In that case, it has a name. 2362 if obj.type.name: 2363 obj.nmoldynname = obj.type.name 2364 else: 2365 raise 2366 2367 except: 2368 2369 if isinstance(obj, Atom): 2370 obj.nmoldynname = 'A' 2371 elif isinstance(obj, AtomCluster): 2372 obj.nmoldynname = 'AC' 2373 elif isinstance(obj, Molecule): 2374 obj.nmoldynname = 'M' 2375 elif isinstance(obj, NucleotideChain): 2376 obj.nmoldynname = 'NC' 2377 elif isinstance(obj, PeptideChain): 2378 obj.nmoldynname = 'PC' 2379 elif isinstance(obj, Protein): 2380 obj.nmoldynname = 'P' 2381 2382 obj.nmoldynname += str(obj.numberOfAtoms()) 2383 2384 # Case of the isotopes. 2385 obj.nmoldynnamelist = sorted([at.name for at in obj.atomList()]) 2386 if not isotope.has_key(obj.nmoldynname): 2387 isotope[obj.nmoldynname] = [obj.nmoldynnamelist] 2388 2389 else: 2390 if not obj.nmoldynnamelist in isotope[obj.nmoldynname]: 2391 isotope[obj.nmoldynname].append(obj.nmoldynnamelist) 2392 2393 universe.nmoldyncontents = {} 2394 2395 # Loop over the object found in the universe. 2396 for obj in universe.objectList(): 2397 2398 if not isinstance(obj,(Atom, AtomCluster, Molecule, NucleotideChain, PeptideChain, Protein)): 2399 continue 2400 2401 if len(isotope[obj.nmoldynname]) > 1: 2402 LogMessage('info','Isotope found for %s' % obj.nmoldynname,['file','console']) 2403 iso = str(isotope[obj.nmoldynname].index(obj.nmoldynnamelist) + 1) 2404 obj.nmoldynname += '_iso' + iso 2405 obj.writeToFile(os.path.join(PREFERENCES.outputfile_path,obj.nmoldynname + '.pdb'), format = 'pdb') 2406 2407 delattr(obj, 'nmoldynnamelist') 2408 2409 if not obj.name: 2410 obj.name = obj.nmoldynname 2411 2412 # If the objet type name is already a key of the |universeContents| dictionnary then just add the name of 2413 # the object to the 'name' set and increment the 'number' counter. 2414 if universe.nmoldyncontents.has_key(obj.nmoldynname): 2415 universe.nmoldyncontents[obj.nmoldynname]['name'].add(obj.name) 2416 universe.nmoldyncontents[obj.nmoldynname]['number'] += 1 2417 2418 else: 2419 universe.nmoldyncontents[obj.nmoldynname] = {'name' : set(['*', obj.name]) , 'number' : 1} 2420 2421 # Case of an object of |Atom| class. 2422 # No additional selection fields: 2423 if isinstance(obj, Atom): 2424 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Atom' 2425 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [0, ['atom',]] 2426 2427 # Case of an object of |AtomCluster| class. 2428 # Two additional selection fields: 2429 # -atomname: the full name of the atom. 2430 # -atomelement: the element of the atom. 2431 elif isinstance(obj, AtomCluster): 2432 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'AtomCluster' 2433 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [1, ['atom','cluster']] 2434 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = set(['*']) 2435 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*']) 2436 for atom in obj.atomList(): 2437 universe.nmoldyncontents[obj.nmoldynname]['atomname'].add(atom.fullName()) 2438 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name) 2439 2440 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomname']) 2441 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement']) 2442 2443 # Case of an object of |Molecule| class. 2444 # Three additional selection fields: 2445 # -atomname: the full name of the atom. 2446 # -atomelement: the element of the atom. 2447 # -chemfragment: the chemical class of the atom. 2448 elif isinstance(obj, Molecule): 2449 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Molecule' 2450 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom','amine', 'hydroxy', 'methyl', 'thiol', 'molecule']] 2451 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = set(['*']) 2452 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*']) 2453 universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'hydroxy', 'methyl', 'thiol'] 2454 2455 for atom in obj.atomList(): 2456 universe.nmoldyncontents[obj.nmoldynname]['atomname'].add(atom.fullName()) 2457 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name) 2458 2459 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomname']) 2460 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement']) 2461 2462 # Case of an object of |NucleotideChain| class. 2463 # Six additional selection fields: 2464 # -atomname: the full name of the atom. 2465 # -atomtype: the type of the atom. 2466 # -atomelement: the element of the atom. 2467 # -nuclname: the name of the nucleotide. 2468 # -nucltype: the type of the nucleotide. 2469 # -chemfragment: some miscellaneous selection keywords. 2470 elif isinstance(obj, NucleotideChain): 2471 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'NucleotideChain' 2472 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [2, ['atom', 'amine', 'residue', 'nucleicacid']] 2473 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*'] 2474 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*']) 2475 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*']) 2476 universe.nmoldyncontents[obj.nmoldynname]['nuclname'] = ['*'] 2477 universe.nmoldyncontents[obj.nmoldynname]['nucltype'] = set(['*']) 2478 universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'bases'] 2479 2480 # Loop over the residues of the chain. 2481 for nucl in obj.residues(): 2482 universe.nmoldyncontents[obj.nmoldynname]['nuclname'].append(nucl.fullName()) 2483 universe.nmoldyncontents[obj.nmoldynname]['nucltype'].add(nucl.symbol) 2484 2485 # Loop over the atoms of the residue. 2486 for atom in nucl.atomList(): 2487 universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName()) 2488 universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name) 2489 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name) 2490 2491 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype']) 2492 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement']) 2493 2494 universe.nmoldyncontents[obj.nmoldynname]['nucltype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['nucltype']) 2495 2496 # Case of an object of |PeptideChain| class. 2497 # Eight additional selection fields: 2498 # -atomname: the full name of the atom. 2499 # -atomtype: the type of the atom. 2500 # -atomelement: the element of the atom. 2501 # -chemfragment: the chemical class of the atom. 2502 # -resname: the name of the residue. 2503 # -restype: the type of the residue. 2504 # -resclass: the chemical class of the residue. 2505 # -chemfragment: some miscellaneous selection keywords. 2506 elif isinstance(obj, PeptideChain): 2507 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'PeptideChain' 2508 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom', 'amine', 'hydroxy', 'methyl', 'thiol', 'residue','chain']] 2509 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*'] 2510 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*']) 2511 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*']) 2512 universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'c_alphas', 'hydroxy', 'methyl', 'thiol'] 2513 universe.nmoldyncontents[obj.nmoldynname]['resname'] = ['*'] 2514 universe.nmoldyncontents[obj.nmoldynname]['restype'] = set(['*']) 2515 universe.nmoldyncontents[obj.nmoldynname]['resclass'] = sorted(residusChemFamily.keys()) 2516 universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'sidechains'] 2517 2518 # Loop over the residues of the chain. 2519 for residue in obj.residues(): 2520 universe.nmoldyncontents[obj.nmoldynname]['resname'].append(residue.fullName()) 2521 universe.nmoldyncontents[obj.nmoldynname]['restype'].add(residue.symbol) 2522 2523 # Loop over the atoms of the residue. 2524 for atom in residue.atomList(): 2525 universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName()) 2526 universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name) 2527 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name) 2528 2529 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype']) 2530 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement']) 2531 2532 universe.nmoldyncontents[obj.nmoldynname]['restype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['restype']) 2533 2534 # Case of an object of |Protein| class. 2535 # Nine additional selection fields: 2536 # -atomname: the full name of the atom. 2537 # -atomtype: the type of the atom. 2538 # -atomelement: the element of the atom. 2539 # -chemfragment: the chemical class of the atom. 2540 # -resname: the name of the residue. 2541 # -restype: the type of the residue. 2542 # -resclass: the chemical class of the residue. 2543 # -chainname: the name of the chain. 2544 # -chemfragment: some miscellaneous selection keywords. 2545 elif isinstance(obj, Protein): 2546 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Protein' 2547 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom', 'amine', 'hydroxy', 'methyl', 'thiol', 'residue','chain', 'protein']] 2548 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*'] 2549 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*']) 2550 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*']) 2551 universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'c_alphas', 'hydroxy', 'methyl', 'thiol'] 2552 universe.nmoldyncontents[obj.nmoldynname]['resname'] = ['*'] 2553 universe.nmoldyncontents[obj.nmoldynname]['restype'] = set(['*']) 2554 universe.nmoldyncontents[obj.nmoldynname]['resclass'] = sorted(residusChemFamily.keys()) 2555 universe.nmoldyncontents[obj.nmoldynname]['chainname'] = ['*'] 2556 universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'sidechains'] 2557 2558 # Loop over the chains of the protein. 2559 for chain in obj: 2560 # Append the chain name. 2561 universe.nmoldyncontents[obj.nmoldynname]['chainname'].append(chain.fullName()) 2562 2563 # Loop over the residues of the chain. 2564 for residue in chain.residues(): 2565 universe.nmoldyncontents[obj.nmoldynname]['resname'].append(residue.fullName()) 2566 universe.nmoldyncontents[obj.nmoldynname]['restype'].add(residue.symbol) 2567 2568 # Loop over the atoms of the residue. 2569 for atom in residue.atomList(): 2570 universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName()) 2571 universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name) 2572 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name) 2573 2574 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype']) 2575 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement']) 2576 2577 universe.nmoldyncontents[obj.nmoldynname]['restype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['restype']) 2578 2579 # The 'name' set value is converted to a sorted list. 2580 for objTypeName in universe.nmoldyncontents.keys(): 2581 universe.nmoldyncontents[objTypeName]['name'] = sorted(universe.nmoldyncontents[objTypeName]['name']) 2582 2583 # Special case of * key used to selected all the nmoldyn type found in the universe. 2584 universe.nmoldyncontents['*'] = {'atomelement' : set(['*']),\ 2585 'number' : universe.numberOfAtoms(),\ 2586 'objectclass' : 'AllClass',\ 2587 'groupinglevel' : [0,['default',]]} 2588 for atom in universe.atomList(): 2589 universe.nmoldyncontents['*']['atomelement'].add(atom.type.name) 2590 universe.nmoldyncontents['*']['atomelement'] = sorted(universe.nmoldyncontents['*']['atomelement'])
2591
2592 -class QVectors:
2593 """Generates a set of QVectors within a given shell. 2594 """
2595 - def __init__(self, universe, generator, qRadii, dq, qVectorsPerShell, qVectorsDirection = None):
2596 """The constructor. 2597 2598 @param universe: the MMTK universe used to define the reciprocal space. 2599 @type universe: a MMTK.Universe subclass object 2600 2601 @param generator: a string being one of '3d isotropic', '2d isotropic' or 'anistropic' the way the q-vectors 2602 should be generated. 2603 @type generator: string 2604 2605 @param qRadii: a list of floats specifying the radii of the shell in which the q vectors have to be 2606 generated. 2607 @type qRadii: list 2608 2609 @param dq: a float specifying the width of a qhsell defined as [|qRadius| - dq/2,|qRadius|+dq/2]. 2610 @type dq: float 2611 2612 @param qVectorsPerShell: an integer specifying the number of q-vectors to generate for each shell. 2613 @type qVectorsPerShell: integer 2614 2615 @param qVectorsDirection: a list of Scientific.Geometry.Vector objects specifying the directions along which 2616 the q-vectors should be generated. If None, the q-vectors generation will be isotropic. 2617 @type qVectorsDirection: list 2618 """ 2619 2620 self.universe = universe 2621 self.generator = generator 2622 self.dq = dq 2623 self.qVectorsPerShell = qVectorsPerShell 2624 self.qVectorsDirection = qVectorsDirection 2625 2626 if self.generator == '3d isotropic': 2627 self.qVectorsDirection = None 2628 2629 else: 2630 2631 if self.qVectorsDirection is None: 2632 raise Error('Directions must be given for %s q vectors generation.' % self.generator) 2633 2634 if self.generator == '2d isotropic': 2635 if len(self.qVectorsDirection) != 2: 2636 raise Error('Two vectors must be given to generate q vectors in a plane.') 2637 2638 v1, v2 = self.qVectorsDirection 2639 if abs(v1.angle(v2)) < 1.e-10: 2640 raise Error('The two vectors are colinear.') 2641 2642 self.qVectorsDirection = qVectorsDirection 2643 2644 self.qRadii = [] 2645 self.qMin = [] 2646 self.qMax = [] 2647 # Loop over all the sphere radius on which the q vectors have to be generated. 2648 for qRadius in sorted(qRadii): 2649 # The case q = 0 produces untrustable results. So skip it. 2650 if qRadius <= 0.0: 2651 LogMessage('warning','No Q vector could be generated for shell with radius Q = %s.' % qRadius,['file','console']) 2652 continue 2653 else: 2654 self.qRadii.append(qRadius) 2655 self.qMin.append(max(qRadius - 0.5*self.dq,0.0)) 2656 self.qMax.append(qRadius + 0.5*self.dq) 2657 2658 # Get the reciprocal basis of the elementary cell. If universe is non-periodic, basis = None 2659 self.reciprocalBasis = universe.reciprocalBasisVectors() 2660 2661 # Non-periodic universe. 2662 if self.reciprocalBasis is None: 2663 self.directMatrix = None 2664 self.reciprocalCellVolume = None 2665 2666 # Periodic universe. 2667 else: 2668 self.reciprocalBasis = [2.0*N.pi*v for v in self.reciprocalBasis] 2669 self.directMatrix = N.array([N.array(v) for v in self.universe.basisVectors()]) 2670 self.directMatrix = Tensor(self.directMatrix)/(2.*N.pi) 2671 self.reciprocalCellVolume = self.reciprocalBasis[0]*(self.reciprocalBasis[1].cross(self.reciprocalBasis[2])) 2672 2673 self.qVectors = {} 2674 2675 # Loop over all the sphere radius on which the q vectors have to be generated. 2676 for index in range(len(self.qRadii)): 2677 qRadius = self.qRadii[index] 2678 self.qVectors[qRadius] = [] 2679 qMin = self.qMin[index] 2680 qMax = self.qMax[index] 2681 2682 if self.generator in ['2d isotropic', '3d isotropic']: 2683 self.__isotropicGeneration(qRadius, qMin, qMax) 2684 2685 elif self.generator == '2d isotropic': 2686 self.__planarIsotropicGeneration(qRadius, qMin, qMax) 2687 2688 elif self.generator == 'anisotropic': 2689 self.__anisotropicGeneration(qRadius, qMin, qMax) 2690 2691 if len(self.qVectors[qRadius]) == 0: 2692 LogMessage('warning','No Q vector could be generated for shell with radius Q = %s.' % qRadius,['file','console']) 2693 del self.qVectors[qRadius] 2694 continue 2695 2696 if not self.qVectors: 2697 raise Error('No Q Vectors generated. Please check your settings.') 2698 2699 self.statistics = N.zeros((len(self.qVectors), 8), typecode = N.Int) 2700 comp = 0 2701 for qRadius in sorted(self.qVectors.keys()): 2702 qVectors = self.qVectors[qRadius] 2703 for qVect in qVectors: 2704 if qVect[0] > 0: 2705 if qVect[1] > 0: 2706 # X+Y+Z+ 2707 if qVect[2] > 0: 2708 self.statistics[comp, 0] += 1 2709 # X+Y+Z- 2710 else: 2711 self.statistics[comp, 1] += 1 2712 else: 2713 # X+Y-Z+ 2714 if qVect[2] > 0: 2715 self.statistics[comp, 2] += 1 2716 # X+Y-Z- 2717 else: 2718 self.statistics[comp, 3] += 1 2719 else: 2720 if qVect[1] > 0: 2721 # X-Y+Z+ 2722 if qVect[2] > 0: 2723 self.statistics[comp, 4] += 1 2724 # X-Y+Z- 2725 else: 2726 self.statistics[comp, 5] += 1 2727 else: 2728 # X-Y-Z+ 2729 if qVect[2] > 0: 2730 self.statistics[comp, 6] += 1 2731 # X-Y-Z- 2732 else: 2733 self.statistics[comp, 7] += 1 2734 comp += 1
2735
2736 - def __isotropicGeneration(self, qRadius, qMin, qMax):
2737 if self.reciprocalCellVolume is None: 2738 self.__makeRandomQList() 2739 2740 else: 2741 reciprocalShellVolume = 4.0*N.pi*qRadius**2 * (qMax - qMin) 2742 nCells = reciprocalShellVolume/self.reciprocalCellVolume 2743 if nCells > 500: 2744 LogMessage('info', 'Random Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console']) 2745 self.__makeRandomQList(qRadius, qMin, qMax) 2746 else: 2747 LogMessage('info', 'Explicit Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console']) 2748 self.__makeExplicitQList(qRadius, qMin, qMax)
2749
2750 - def __anisotropicGeneration(self, qRadius, qMin, qMax):
2751 LogMessage('info', 'Explicit Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console']) 2752 self.__makeExplicitQList(qRadius, qMin, qMax)
2753
2754 - def __makeRandomQList(self, qRadius, qMin, qMax):
2755 2756 # Generates randomly, 'number' q vectors. 2757 # This counter is used to avoid infinite loop that can occur during the q vector generation. 2758 redundancyCounter = 0 2759 while len(self.qVectors[qRadius]) < self.qVectorsPerShell: 2760 # Generate a normalized random vector. 2761 # The actual q vector. 2762 q = randomVector(self.qVectorsDirection)*uniform(qMin, qMax) 2763 if self.reciprocalBasis is not None: 2764 # Its projection on the reciprocal basis. If the vector is not an exact reciprocal space vector, it is 2765 # rounded to the lower corner of its corresponding reciprocal cell. 2766 q = [round(v) for v in N.array(self.directMatrix*q)] 2767 # The approximated reciprocal space vector. 2768 q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2] 2769 2770 # Checks whether the length of q vector is in the user-defined length interval [lower, upper] 2771 # If not, regenerates a new q vector 2772 if qMin < q.length() <= qMax: 2773 # The generated q vector is already in the q vector list so skip it. 2774 if q in self.qVectors[qRadius]: 2775 redundancyCounter += 1 2776 # After 5000 trials to generate a new q vector, gives up and returns the q vector list in its current state. 2777 if redundancyCounter == 5000: 2778 return 2779 continue 2780 self.qVectors[qRadius].append(q) 2781 # A new q vector could be found, so reinit the redundancy counter. 2782 redundancyCounter = 0 2783 else: 2784 redundancyCounter += 1 2785 # After 5000 trials to generate a new q vector, gives up and returns the q vector list in its current state. 2786 if redundancyCounter == 5000: 2787 return
2788
2789 - def __makeExplicitQList(self, qRadius, qMin, qMax):
2790 2791 # isotropic case 2792 if self.generator == '3d isotropic': 2793 upperLimit = [int(qMax/v.length()) for v in self.reciprocalBasis] 2794 # Computes all the possible combination of basis vectors and keep the ones that falls in the selected shell. 2795 for l in range(-upperLimit[0], upperLimit[0]+1): 2796 for m in range(-upperLimit[1], upperLimit[1]+1): 2797 for n in range(-upperLimit[2], upperLimit[2]+1): 2798 if l == m == n == 0: 2799 continue 2800 2801 q = l*self.reciprocalBasis[0] + m*self.reciprocalBasis[1] + n*self.reciprocalBasis[2] 2802 2803 if qMin < q.length() <= qMax: 2804 if q in self.qVectors[qRadius]: 2805 continue 2806 self.qVectors[qRadius].append(q) 2807 2808 elif self.generator == '2d isotropic': 2809 upperLimit = [int(qMax/v.length()) for v in self.qVectorsDirection] 2810 for l in range(-upperLimit[0], upperLimit[0]+1): 2811 for m in range(-upperLimit[1], upperLimit[1]+1): 2812 if l == m == 0: 2813 continue 2814 2815 # The actual q vector. 2816 q = l*self.qVectorsDirection[0] + m*self.qVectorsDirection[1] 2817 2818 # Its projection on the reciprocal basis vector. If the vector is not an exact reciprocal space vector, it is 2819 # rounded to the lower corner of its corresponding reciprocal cell. 2820 q = [round(v) for v in N.array(self.directMatrix*q)] 2821 # The approximated reciprocal space vector. 2822 q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2] 2823 2824 if qMin < q.length() <= qMax: 2825 if q in self.qVectors[qRadius]: 2826 continue 2827 self.qVectors[qRadius].append(q) 2828 2829 elif self.generator == 'anisotropic': 2830 2831 for qVect in self.qVectorsDirection: 2832 upperLimit = int(qMax/qVect.length()) 2833 for l in range(-upperLimit, upperLimit + 1): 2834 if l == 0: 2835 continue 2836 2837 # Generate a normalized random vector. 2838 # The actual q vector. 2839 q = l*qVect 2840 2841 # Its projection on the reciprocal basis vector. If the vector is not an exact reciprocal space vector, it is 2842 # rounded to the lower corner of its corresponding reciprocal cell. 2843 q = [round(v) for v in N.array(self.directMatrix*q)] 2844 # The approximated reciprocal space vector. 2845 q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2] 2846 2847 if qMin < q.length() <= qMax: 2848 if q in self.qVectors[qRadius]: 2849 continue 2850 self.qVectors[qRadius].append(q) 2851 2852 # If there are more q vectors generated than the user-defined number of q vectors per shell 2853 # then pick up randomly |number| q vectors among the generated q vectors. 2854 if len(self.qVectors[qRadius]) > self.qVectorsPerShell: 2855 LogMessage('warning', '%s qVectors generated for shell with radius Q = %s. Will be reduced to %s.' % (len(self.qVectors[qRadius]), qRadius, self.qVectorsPerShell), ['file','console']) 2856 self.qVectors[qRadius] = sample(self.qVectors[qRadius], self.qVectorsPerShell)
2857