Home | Trees | Indices | Help |
|
---|
|
1 """Collections of classes for the determination of NMR-related properties. 2 3 Classes: 4 5 * OrderParameter : sets up an order parameter analysis. 6 * OrderParameterContactModel : sets up an order parameter analysis using the contact model approach. 7 """ 8 9 # The python distribution modules 10 import copy 11 import math 12 from numpy import linalg 13 import operator 14 import os 15 import sys 16 from time import asctime 17 from timeit import default_timer 18 19 # The ScientificPython modules 20 from Scientific import N 21 from Scientific.Geometry import Vector 22 from Scientific.IO.NetCDF import NetCDFFile 23 24 # The MMTK distribution modules 25 from MMTK.Collections import Collection 26 from MMTK import Units 27 from MMTK.NucleicAcids import NucleotideChain 28 from MMTK.Proteins import PeptideChain, Protein 29 30 # The nMOLDYN modules 31 from nMOLDYN.Analysis.Analysis import Analysis 32 from nMOLDYN.Core.Error import Error 33 from nMOLDYN.Core.IOFiles import convertNetCDFToASCII 34 from nMOLDYN.Core.Logger import LogMessage 35 from nMOLDYN.Core.Mathematics import correlation 36 from nMOLDYN_S2 import computeS2 37 38 ##################################################################################### 39 # ORDER PARAMETER ANALYSIS 40 #####################################################################################42 """Sets up an order parameter analysis. 43 44 A Subclass of nMOLDYN.Analysis.Analysis. 45 46 Constructor: OrderParameter(|parameters| = None) 47 48 Arguments: 49 50 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters. 51 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class. 52 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame 53 number to consider, 'last' is an integer specifying the last frame number to consider and 54 'step' is an integer specifying the step number between two frames. 55 * group -- a selection string specifying the groups of atoms that will define the vectors on which the 56 analysis will be computed. Each group must contain two and only two atoms. 57 * atomorder -- a string of the form 'atom1,atom2,atom3' where 'atom1', 'atom2' and 'atom3' are 58 respectively the MMTK atom names of the atoms in the way they should be ordered. 59 * op -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension 60 instead of the '.nc' extension. 61 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis. 62 63 Running modes: 64 65 - To run the analysis do: a.runAnalysis() where a is the analysis object. 66 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object. 67 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object. 68 69 70 Comments: 71 72 - This code is based on a first implementation made by Vania Calandrini. 73 """ 74 75 # The minimal set of input parameters names necessary to perform the analysis. 76 inputParametersNames = ('trajectory',\ 77 'timeinfo',\ 78 'bond',\ 79 'atomorder',\ 80 'op',\ 81 'pyroserver',\ 82 ) 83 84 shortName = 'OP' 85 86 # Indicate whether this analysis can be estimated or not. 87 canBeEstimated = True 88271 272 ##################################################################################### 273 # ORDER PARAMETER ANALYSIS USING CONTACT MODEL 274 #####################################################################################90 """The constructor. Insures that the class can not be instanciated directly from here. 91 """ 92 raise Error('This class can not be instanciated directly.')9395 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 96 """ 97 98 # The input parameters are parsed. 99 self.parseInputParameters() 100 101 if self.pyroServer != 'monoprocessor': 102 if self.statusBar is not None: 103 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.') 104 self.statusBar = None 105 106 self.buildTimeInfo() 107 108 self.bond = self.groupSelection(self.universe, self.bond) 109 self.bond = [g for g in self.bond if g.numberOfAtoms() == 2] 110 111 if self.atomOrder is None: 112 self.bond = [sorted(g, key = operator.attrgetter('name')) for g in self.bond] 113 114 else: 115 orderedBonds = [] 116 for g in self.bond: 117 gr = [] 118 for atName in self.atomOrder: 119 found = False 120 for at in g: 121 if at.name == atName: 122 found = True 123 gr.append(at) 124 if len(gr) == 2: 125 orderedBonds.append(gr) 126 gr = [] 127 else: 128 if not found: 129 raise Error('No atom corresponding %s could be found.' % atName) 130 131 self.bond = orderedBonds 132 133 if len(self.bond) == 0: 134 raise Error('Your group selection does not contain any bonds.') 135 136 self.group = self.bond 137 self.nBonds = self.nGroups = len(self.bond) 138 139 self.preLoadTrajectory('atom') 140 141 self.referenceDirection = Vector(0.0,0.0,1.0) 142 143 # The results are stored in dictionnary because the order in which the bonds are treated does not 144 # always follow the structure. 145 self.P2 = {} 146 self.S2 = {} 147 148 if self.pyroServer != 'monoprocessor': 149 # The attribute trajectory is removed because it can not be pickled by Pyro. 150 delattr(self, 'trajectory')151153 """Calculates the contribution for one group. 154 155 @param bondIndex: the index of the group in |self.bond| list. 156 @type bondIndex: integer. 157 158 @param trajname: the name of the trajectory file name. 159 @type trajname: string 160 """ 161 162 # The atoms forming the bond. 163 at1, at2 = self.bond[bondIndex] 164 165 if self.preLoadedTraj is None: 166 if self.pyroServer == 'monoprocessor': 167 t = self.trajectory 168 169 else: 170 # Load the whole trajectory set. 171 t = Trajectory(None, trajname, 'r') 172 173 at1Traj = t.readParticleTrajectory(at1, first = self.first, last = self.last, skip = self.skip) 174 at2Traj = t.readParticleTrajectory(at2, first = self.first, last = self.last, skip = self.skip) 175 176 else: 177 at1Traj = self.preLoadedTraj[at1] 178 at2Traj = self.preLoadedTraj[at2] 179 180 costheta = N.zeros((self.nFrames,), N.Float) 181 sinphi = N.zeros((self.nFrames,), N.Float) 182 cosphi = N.zeros((self.nFrames,), N.Float) 183 sintheta = N.zeros((self.nFrames,), N.Float) 184 185 for frameIndex in range(self.nFrames): 186 187 pVect = Vector(self.universe.distanceVector(at1Traj[frameIndex], at2Traj[frameIndex])) 188 pVect = pVect.normal() 189 190 costheta[frameIndex] = pVect * self.referenceDirection 191 sintheta[frameIndex] = pVect.cross(self.referenceDirection).length() 192 cosphi[frameIndex] = (pVect[0]/sintheta[frameIndex]) 193 sinphi[frameIndex] = (pVect[1]/sintheta[frameIndex]) 194 195 tr2 = 3.0*costheta**2 - 1. 196 cos2phi = 2.0*cosphi**2-1. 197 sin2phi = 2.0*sinphi*cosphi 198 cossintheta = costheta*sintheta 199 sintheta_sq = sintheta**2 200 201 # calcul de <P2(cos(theta))> en termes de somme de fonctions de correlation 202 # d'harmoniques spheriques (theoreme d'addition des harmoniques spheriques). 203 P2 = (0.25*correlation(tr2) + 3.00*correlation(cosphi*cossintheta) + \ 204 3.00*correlation(sinphi*cossintheta) + 0.75*correlation(cos2phi*sintheta_sq) + \ 205 0.75*correlation(sin2phi*sintheta_sq)) 206 207 # calcul du parametre d'ordre S^2 (limite pour t->infini de <P2(cos(theta))>). 208 S2 = (0.75 * (N.sum(cos2phi*sintheta_sq)**2 + N.sum(sin2phi*sintheta_sq)**2) + \ 209 3.00 * (N.sum(cosphi*cossintheta)**2 + N.sum(sinphi*cossintheta)**2) + 210 0.25 * N.sum(tr2)**2) / self.nFrames**2 211 212 return P2, S2213215 """ 216 """ 217 # The atoms forming the bond. 218 at1, at2 = self.bond[bondIndex] 219 if isinstance(at1.topLevelChemicalObject(), (Protein, PeptideChain, NucleotideChain)): 220 bondNumber = at1.parent.parent.sequence_number 221 else: 222 bondNumber = bondIndex + 1 223 224 # calcul de <P2(cos(theta))> en termes de somme de fonctions de correlation 225 # d'harmoniques spheriques (theoreme d'addition des harmoniques spheriques). 226 self.P2[bondNumber] = x[0] 227 228 # calcul du parametre d'ordre S^2 (limite pour t->infini de <P2(cos(theta))>). 229 self.S2[bondNumber] = x[1]230232 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 233 """ 234 235 outputFile = NetCDFFile(self.outputOP, 'w') 236 outputFile.title = self.__class__.__name__ 237 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime() 238 239 outputFile.createDimension('NBONDS', self.nBonds) 240 outputFile.createDimension('NTIMES', self.nFrames) 241 242 TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',)) 243 TIMES[:] = self.times 244 TIMES.units = 'ps' 245 246 BONDS = outputFile.createVariable('bond_id', N.Int, ('NBONDS',)) 247 P2 = outputFile.createVariable('p2', N.Float, ('NBONDS','NTIMES')) 248 P2AVG = outputFile.createVariable('p2-bondavg', N.Float, ('NTIMES',)) 249 S2 = outputFile.createVariable('s2', N.Float, ('NBONDS',)) 250 p2Avg = N.zeros((self.nFrames),N.Float) 251 comp = 0 252 for ind in sorted(self.S2.keys()): 253 BONDS[comp] = ind 254 S2[comp] = self.S2[ind] 255 P2[comp,:] = self.P2[ind] 256 N.add(p2Avg,self.P2[ind],p2Avg) 257 comp += 1 258 259 P2AVG[:] = p2Avg/float(self.nBonds) 260 261 asciiVar = sorted(outputFile.variables.keys()) 262 263 outputFile.close() 264 265 self.toPlot = {'netcdf' : self.outputOP, 'xVar' : 'pair', 'yVar' : 'S2'} 266 267 # Creates an ASCII version of the NetCDF output file. 268 convertNetCDFToASCII(inputFile = self.outputOP,\ 269 outputFile = os.path.splitext(self.outputOP)[0] + '.cdl',\ 270 variables = asciiVar)276 """Sets up an order parameter analysis using the contact model . 277 278 A Subclass of nMOLDYN.Analysis.Analysis. 279 280 Constructor: OrderParameterContactModel(|parameters| = None) 281 282 Arguments: 283 284 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters. 285 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class. 286 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame 287 number to consider, 'last' is an integer specifying the last frame number to consider and 288 'step' is an integer specifying the step number between two frames. 289 * opcm -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension 290 instead of the '.nc' extension. 291 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis. 292 293 Running modes: 294 295 - To run the analysis do: a.runAnalysis() where a is the analysis object. 296 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object. 297 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object. 298 299 Comments: 300 301 - This code is adapted from the s2predict code developped by F. Zhang and R. Bruschweiler and available in: 302 http://nmr.clarku.edu/software/S2/s2predict.html 303 304 - For more details about the method: Zhang, F., Bruschweiler, R. J. AM. Chem. Soc. 2002, 124, 12654-12655. 305 """ 306 307 # The minimal set of input parameters names necessary to perform the analysis. 308 inputParametersNames = ('trajectory',\ 309 'timeinfo',\ 310 'opcm',\ 311 'pyroserver',\ 312 ) 313 314 shortName = 'OPCM' 315 316 # Indicate whether this analysis can be estimated or not. 317 canBeEstimated = True 318484320 """The constructor. Insures that the class can not be instanciated directly from here. 321 """ 322 raise Error('This class can not be instanciated directly.')323325 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 326 """ 327 328 # The input parameters are parsed. 329 self.parseInputParameters() 330 331 if self.pyroServer != 'monoprocessor': 332 if self.statusBar is not None: 333 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.') 334 self.statusBar = None 335 336 self.buildTimeInfo() 337 338 self.preLoadTrajectory('frame') 339 340 if self.universe.basisVectors() is None: 341 raise Error('The universe must be periodic for this kind of calculation.') 342 343 self.S2 = {} 344 self.sequence = {} 345 self.hnLookup = {} 346 for obj in self.universe.objectList(): 347 if isinstance(obj, (PeptideChain, Protein)): 348 self.sequence[obj] = [] 349 self.hnLookup[obj] = [] 350 comp = 1 351 for r in obj.residues()[1:]: 352 if r.symbol.lower() != 'pro': 353 self.sequence[obj].append(r.sequence_number) 354 self.hnLookup[obj].append(comp) 355 comp += 1 356 self.S2[obj] = N.zeros((len(self.sequence[obj]),self.nFrames), N.Float) 357 358 # Case where no protein or peptide chain was found in the universe. 359 if not self.sequence: 360 raise Error('The universe must contains at least one peptide chain to perform the analysis.') 361 362 if self.pyroServer != 'monoprocessor': 363 # The attribute trajectory is removed because it can not be pickled by Pyro. 364 delattr(self, 'trajectory')365367 """Calculates the contribution for one group. 368 369 @param frameIndex: the index of the frame in |self.frameIndexes| array. 370 @type frameIndex: integer. 371 372 @param trajname: the name of the trajectory file name. 373 @type trajname: string 374 """ 375 376 frame = self.frameIndexes[frameIndex] 377 378 if self.preLoadedTraj is None: 379 if self.pyroServer == 'monoprocessor': 380 t = self.trajectory 381 382 else: 383 # Load the whole trajectory set. 384 t = Trajectory(None, trajname, 'r') 385 386 conf = t.configuration[frame] 387 self.universe.setConfiguration(conf) 388 389 else: 390 # The configuration is updated to the current trajectory step. 391 self.universe.setConfiguration(self.preLoadedTraj[frameIndex]) 392 393 directCell = N.ravel(N.array([v for v in self.universe.basisVectors()])) 394 reverseCell = N.ravel(N.transpose(N.array([v for v in self.universe.reciprocalBasisVectors()]))) 395 396 S2 = {} 397 for obj in self.universe.objectList(): 398 if isinstance(obj, (PeptideChain, Protein)): 399 S2[obj] = N.zeros((len(self.sequence[obj])), N.Float) 400 401 for obj in self.sequence.keys(): 402 for v in range(len(self.hnLookup[obj])): 403 aa_no = self.hnLookup[obj][v] 404 aa_no_minus = aa_no - 1 405 406 H = obj.backbone()[aa_no].H 407 O = obj.backbone()[aa_no_minus].O 408 409 heavyAtoms = Collection() 410 for m in range(len(obj.residues())): 411 if (m != aa_no) and (m != aa_no_minus): 412 for atom in obj.residues()[m].atomList(): 413 if atom.symbol != 'H': 414 heavyAtoms.addObject(atom) 415 416 # The indexes of the selected atoms. 417 indexes = N.zeros((heavyAtoms.numberOfAtoms(),), typecode = N.Int) 418 comp = 0 419 for at in heavyAtoms.atomList(): 420 indexes[comp] = at.index 421 comp += 1 422 423 val = computeS2(H.index,\ 424 O.index,\ 425 self.universe.contiguousObjectConfiguration().array,\ 426 directCell,\ 427 reverseCell,\ 428 indexes,\ 429 Units.Ang) 430 431 S2[obj][v] = val 432 433 if self.preLoadedTraj is not None: 434 del self.preLoadedTraj[frameIndex] 435 436 return S2437439 """ 440 """ 441 for obj in self.universe.objectList(): 442 if isinstance(obj, (PeptideChain, Protein)): 443 self.S2[obj][:,frameIndex] = x[obj]444446 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 447 """ 448 449 # The NetCDF output file is opened for writing. 450 outputFile = NetCDFFile(self.outputOPCM, 'w') 451 outputFile.title = self.__class__.__name__ 452 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime() 453 454 outputFile.createDimension('TIME', self.nFrames) 455 TIMES = outputFile.createVariable('time', N.Float, ('TIME',)) 456 TIMES[:] = self.times[:] 457 TIMES.units = 'ps' 458 459 for obj in self.sequence: 460 outputFile.createDimension('SEQ%s' % obj.name, len(self.sequence[obj])) 461 462 SEQUENCE = outputFile.createVariable('%s_sequence' % obj.name, N.Int, ('SEQ%s' % obj.name,)) 463 SEQUENCE[:] = N.array(self.sequence[obj]) 464 SEQUENCE.units = 'unitless' 465 466 S2 = outputFile.createVariable('%s_s2' % obj.name, N.Float, ('SEQ%s' % obj.name,'TIME')) 467 S2[:] = self.S2[obj][:,:] 468 S2.units = 'unitless' 469 470 S2AVG = outputFile.createVariable('%s_s2_timeavg' % obj.name, N.Float, ('SEQ%s' % obj.name,)) 471 S2AVG[:] = self.S2[obj][:,:].sum(1)/float(self.nFrames) 472 S2AVG.units = 'unitless' 473 474 asciiVar = sorted(outputFile.variables.keys()) 475 476 outputFile.close() 477 478 self.toPlot = None 479 480 # Creates an ASCII version of the NetCDF output file. 481 convertNetCDFToASCII(inputFile = self.outputOPCM,\ 482 outputFile = os.path.splitext(self.outputOPCM)[0] + '.cdl',\ 483 variables = asciiVar)
Home | Trees | Indices | Help |
|
---|
Generated by Epydoc 3.0.1 on Thu Oct 8 16:57:07 2009 | http://epydoc.sourceforge.net |