specific_analyses.frame

65 """Class containing the specific methods of the Frame Order theories.""" 66

67 - def __init__(self):

68 """Initialise the class by placing API_common methods into the API.""" 69 70 # Execute the base class __init__ method. 71 super(Frame_order, self).__init__() 72 73 # Place methods into the API. 74 self.overfit_deselect = self._overfit_deselect_dummy 75 self.return_conversion_factor = self._return_no_conversion_factor 76 self.set_param_values = self._set_param_values_global 77 78 # Set up the global parameters. 79 self.PARAMS.add('ave_pos_x', scope='global', units='rad', desc='The average position x translation', py_type=float, set='params', err=True, sim=True) 80 self.PARAMS.add('ave_pos_y', scope='global', units='rad', desc='The average position y translation', py_type=float, set='params', err=True, sim=True) 81 self.PARAMS.add('ave_pos_z', scope='global', units='rad', desc='The average position z translation', py_type=float, set='params', err=True, sim=True) 82 self.PARAMS.add('ave_pos_alpha', scope='global', units='rad', desc='The average position alpha Euler angle', py_type=float, set='params', err=True, sim=True) 83 self.PARAMS.add('ave_pos_beta', scope='global', units='rad', desc='The average position beta Euler angle', py_type=float, set='params', err=True, sim=True) 84 self.PARAMS.add('ave_pos_gamma', scope='global', units='rad', desc='The average position gamma Euler angle', py_type=float, set='params', err=True, sim=True) 85 self.PARAMS.add('eigen_alpha', scope='global', units='rad', desc='The Eigenframe alpha Euler angle', py_type=float, set='params', err=True, sim=True) 86 self.PARAMS.add('eigen_beta', scope='global', units='rad', desc='The Eigenframe beta Euler angle', py_type=float, set='params', err=True, sim=True) 87 self.PARAMS.add('eigen_gamma', scope='global', units='rad', desc='The Eigenframe gamma Euler angle', py_type=float, set='params', err=True, sim=True) 88 self.PARAMS.add('axis_theta', scope='global', units='rad', desc='The cone axis polar angle (for the isotropic cone model)', py_type=float, set='params', err=True, sim=True) 89 self.PARAMS.add('axis_phi', scope='global', units='rad', desc='The cone axis azimuthal angle (for the isotropic cone model)', py_type=float, set='params', err=True, sim=True) 90 self.PARAMS.add('cone_theta_x', scope='global', units='rad', desc='The pseudo-ellipse cone opening half-angle for the x-axis', py_type=float, set='params', err=True, sim=True) 91 self.PARAMS.add('cone_theta_y', scope='global', units='rad', desc='The pseudo-ellipse cone opening half-angle for the y-axis', py_type=float, set='params', err=True, sim=True) 92 self.PARAMS.add('cone_theta', scope='global', units='rad', desc='The isotropic cone opening half-angle', py_type=float, set='params', err=True, sim=True) 93 self.PARAMS.add('cone_s1', scope='global', units='', desc='The isotropic cone order parameter', py_type=float, set='params', err=True, sim=True) 94 self.PARAMS.add('cone_sigma_max', scope='global', units='rad', desc='The torsion angle', py_type=float, set='params', err=True, sim=True) 95 self.PARAMS.add('params', scope='global', desc='The model parameters', py_type=list) 96 97 # Add minimisation structures. 98 self.PARAMS.add_min_data(min_stats_global=True)

99 100

101 - def _assemble_limit_arrays(self):

102 """Assemble and return the limit vectors. 103 104 @return: The lower and upper limit vectors. 105 @rtype: numpy rank-1 array, numpy rank-1 array 106 """ 107 108 # Init. 109 lower = zeros(len(cdp.params), float64) 110 upper = 2.0*pi * ones(len(cdp.params), float64) 111 112 # Return the arrays. 113 return lower, upper

114 115

116 - def _assemble_param_vector(self, sim_index=None):

117 """Assemble and return the parameter vector. 118 119 @return: The parameter vector. 120 @rtype: numpy rank-1 array 121 @keyword sim_index: The Monte Carlo simulation index. 122 @type sim_index: int 123 """ 124 125 # Initialise. 126 param_vect = [] 127 128 # Pivot point. 129 if not self._pivot_fixed(): 130 for i in range(3): 131 param_vect.append(cdp.pivot[i]) 132 133 # Normal values. 134 if sim_index == None: 135 # Average position translation. 136 if not self._translation_fixed(): 137 param_vect.append(cdp.ave_pos_x) 138 param_vect.append(cdp.ave_pos_y) 139 param_vect.append(cdp.ave_pos_z) 140 141 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 142 if cdp.model in ['free rotor', 'iso cone, free rotor']: 143 param_vect.append(cdp.ave_pos_beta) 144 param_vect.append(cdp.ave_pos_gamma) 145 else: 146 param_vect.append(cdp.ave_pos_alpha) 147 param_vect.append(cdp.ave_pos_beta) 148 param_vect.append(cdp.ave_pos_gamma) 149 150 # Frame order eigenframe - the full frame. 151 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 152 param_vect.append(cdp.eigen_alpha) 153 param_vect.append(cdp.eigen_beta) 154 param_vect.append(cdp.eigen_gamma) 155 156 # Frame order eigenframe - the isotropic cone axis. 157 elif cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 158 param_vect.append(cdp.axis_theta) 159 param_vect.append(cdp.axis_phi) 160 161 # Cone parameters - pseudo-elliptic cone parameters. 162 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 163 param_vect.append(cdp.cone_theta_x) 164 param_vect.append(cdp.cone_theta_y) 165 166 # Cone parameters - single isotropic angle or order parameter. 167 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 168 param_vect.append(cdp.cone_theta) 169 elif cdp.model in ['iso cone, free rotor']: 170 param_vect.append(cdp.cone_s1) 171 172 # Cone parameters - torsion angle. 173 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 174 param_vect.append(cdp.cone_sigma_max) 175 176 # Simulation values. 177 else: 178 # Average position translation. 179 if not self._translation_fixed(): 180 param_vect.append(cdp.ave_pos_x_sim[sim_index]) 181 param_vect.append(cdp.ave_pos_y_sim[sim_index]) 182 param_vect.append(cdp.ave_pos_z_sim[sim_index]) 183 184 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 185 if cdp.model in ['free rotor', 'iso cone, free rotor']: 186 param_vect.append(cdp.ave_pos_beta_sim[sim_index]) 187 param_vect.append(cdp.ave_pos_gamma_sim[sim_index]) 188 else: 189 param_vect.append(cdp.ave_pos_alpha_sim[sim_index]) 190 param_vect.append(cdp.ave_pos_beta_sim[sim_index]) 191 param_vect.append(cdp.ave_pos_gamma_sim[sim_index]) 192 193 # Frame order eigenframe - the full frame. 194 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 195 param_vect.append(cdp.eigen_alpha_sim[sim_index]) 196 param_vect.append(cdp.eigen_beta_sim[sim_index]) 197 param_vect.append(cdp.eigen_gamma_sim[sim_index]) 198 199 # Frame order eigenframe - the isotropic cone axis. 200 elif cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 201 param_vect.append(cdp.axis_theta_sim[sim_index]) 202 param_vect.append(cdp.axis_phi_sim[sim_index]) 203 204 # Cone parameters - pseudo-elliptic cone parameters. 205 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 206 param_vect.append(cdp.cone_theta_x_sim[sim_index]) 207 param_vect.append(cdp.cone_theta_y_sim[sim_index]) 208 209 # Cone parameters - single isotropic angle or order parameter. 210 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 211 param_vect.append(cdp.cone_theta_sim[sim_index]) 212 elif cdp.model in ['iso cone, free rotor']: 213 param_vect.append(cdp.cone_s1_sim[sim_index]) 214 215 # Cone parameters - torsion angle. 216 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 217 param_vect.append(cdp.cone_sigma_max_sim[sim_index]) 218 219 # Return as a numpy array. 220 return array(param_vect, float64)

221 222

223 - def _assemble_scaling_matrix(self, data_types=None, scaling=True):

224 """Create and return the scaling matrix. 225 226 @keyword data_types: The base data types used in the optimisation. This list can contain the elements 'rdc', 'pcs' or 'tensor'. 227 @type data_types: list of str 228 @keyword scaling: If False, then the identity matrix will be returned. 229 @type scaling: bool 230 @return: The square and diagonal scaling matrix. 231 @rtype: numpy rank-2 array 232 """ 233 234 # Initialise. 235 scaling_matrix = identity(self._param_num(), float64) 236 237 # Return the identity matrix. 238 if not scaling: 239 return scaling_matrix 240 241 # The pivot point. 242 if not self._pivot_fixed(): 243 for i in range(3): 244 scaling_matrix[i, i] = 1e2 245 246 # Return the matrix. 247 return scaling_matrix

248 249

250 - def _average_position(self, pivot='com', translation=True):

251 """Set up the mechanics of the average domain position. 252 253 @keyword pivot: What to use as the motional pivot. This can be 'com' for the centre of mass of the moving domain, or 'motional' to link the pivot of the motion to the rotation of the average domain position. 254 @type pivot: str 255 @keyword translation: If True, translation to the average domain position will be allowed. If False, then translation will not occur. 256 @type translation: bool 257 """ 258 259 # Test if the current data pipe exists. 260 pipes.test() 261 262 # Check the pivot value. 263 if pivot not in ['com', 'motional']: 264 raise RelaxError("The pivot for the rotation to the average domain position must be either 'com' or 'motional'.") 265 266 # Store the data. 267 cdp.ave_pos_pivot = pivot 268 cdp.ave_pos_translation = translation

269 270

271 - def _base_data_types(self):

272 """Determine all the base data types. 273 274 The base data types can include:: 275 - 'rdc', residual dipolar couplings. 276 - 'pcs', pseudo-contact shifts. 277 - 'noesy', NOE restraints. 278 - 'tensor', alignment tensors. 279 280 @return: A list of all the base data types. 281 @rtype: list of str 282 """ 283 284 # Array of data types. 285 list = [] 286 287 # RDC search. 288 for interatom in interatomic_loop(selection1=self._domain_moving()): 289 if hasattr(interatom, 'rdc'): 290 list.append('rdc') 291 break 292 293 # PCS search. 294 for spin in spin_loop(selection=self._domain_moving()): 295 if hasattr(spin, 'pcs'): 296 list.append('pcs') 297 break 298 299 # Alignment tensor search. 300 if not ('rdc' in list or 'pcs' in list) and hasattr(cdp, 'align_tensors'): 301 list.append('tensor') 302 303 # No data is present. 304 if not list: 305 raise RelaxError("Neither RDCs, PCSs nor alignment tensor data is present.") 306 307 # Return the list. 308 return list

309 310

311 - def _check_rdcs(self, interatom, spin1, spin2):

312 """Check if the RDCs for the given interatomic data container should be used. 313 314 @param interatom: The interatomic data container. 315 @type interatom: InteratomContainer instance 316 @param spin1: The first spin container. 317 @type spin1: SpinContainer instance 318 @param spin2: The second spin container. 319 @type spin2: SpinContainer instance 320 @return: True if the RDCs should be used, False otherwise. 321 """ 322 323 # Skip deselected spins. 324 if not spin1.select or not spin2.select: 325 return False 326 327 # Only use interatomic data containers with RDC data. 328 if not hasattr(interatom, 'rdc'): 329 return False 330 331 # RDC data exists but the interatomic vectors are missing? 332 if not hasattr(interatom, 'vector'): 333 # Throw a warning. 334 warn(RelaxWarning("RDC data exists but the interatomic vectors are missing, skipping the spin pair '%s' and '%s'." % (interatom.spin_id1, interatom.spin_id2))) 335 336 # Jump to the next spin. 337 return False 338 339 # Skip non-Me pseudo-atoms for the first spin. 340 if hasattr(spin1, 'members') and len(spin1.members) != 3: 341 warn(RelaxWarning("Only methyl group pseudo atoms are supported due to their fast rotation, skipping the spin pair '%s' and '%s'." % (interatom.spin_id1, interatom.spin_id2))) 342 return False 343 344 # Skip non-Me pseudo-atoms for the second spin. 345 if hasattr(spin2, 'members') and len(spin2.members) != 3: 346 warn(RelaxWarning("Only methyl group pseudo atoms are supported due to their fast rotation, skipping the spin pair '%s' and '%s'." % (interatom.spin_id1, interatom.spin_id2))) 347 return False 348 349 # Checks. 350 if not hasattr(spin1, 'isotope'): 351 raise RelaxSpinTypeError(interatom.spin_id1) 352 if not hasattr(spin2, 'isotope'): 353 raise RelaxSpinTypeError(interatom.spin_id2) 354 if not hasattr(interatom, 'r'): 355 raise RelaxNoValueError("averaged interatomic distance") 356 357 # Everything is ok. 358 return True

359 360

361 - def _domain_moving(self):

362 """Return the spin ID string corresponding to the moving domain. 363 364 @return: The spin ID string defining the moving domain. 365 @rtype: str 366 """ 367 368 # Check that the domain is defined. 369 if not hasattr(cdp, 'domain'): 370 raise RelaxError("No domains have been defined. Please use the domain user function.") 371 372 # Only support for 2 domains. 373 if len(list(cdp.domain.keys())) > 2: 374 raise RelaxError("Only two domains are supported in the frame order analysis.") 375 376 # Reference domain not set yet, so return nothing. 377 if not hasattr(cdp, 'ref_domain'): 378 return None 379 380 # Loop over the domains. 381 for id in list(cdp.domain.keys()): 382 # Reference domain. 383 if id == cdp.ref_domain: 384 continue 385 386 # Return the ID. 387 return cdp.domain[id]

388 389

390 - def _grid_row(self, incs, lower, upper, dist_type=None, end_point=True):

391 """Set up a row of the grid search for a given parameter. 392 393 @param incs: The number of increments. 394 @type incs: int 395 @param lower: The lower bounds. 396 @type lower: float 397 @param upper: The upper bounds. 398 @type upper: float 399 @keyword dist_type: The spacing or distribution type between grid nodes. If None, then a linear grid row is returned. If 'acos', then an inverse cos distribution of points is returned (e.g. for uniform sampling in angular space). 400 @type dist_type: None or str 401 @keyword end_point: A flag which if False will cause the end point to be removed. 402 @type end_point: bool 403 @return: The row of the grid. 404 @rtype: list of float 405 """ 406 407 # Init. 408 row = [] 409 410 # Linear grid. 411 if dist_type == None: 412 # Loop over the increments. 413 for i in range(incs): 414 # The row. 415 row.append(lower + i * (upper - lower) / (incs - 1.0)) 416 417 # Inverse cos distribution. 418 elif dist_type == 'acos': 419 # Generate the increment values of v from cos(upper) to cos(lower). 420 v = zeros(incs, float64) 421 val = (cos(lower) - cos(upper)) / (incs - 1.0) 422 for i in range(incs): 423 v[-i-1] = cos(upper) + float(i) * val 424 425 # Generate the distribution. 426 row = arccos(v) 427 428 # Remove the last point. 429 if not end_point: 430 row = row[:-1] 431 432 # Return the row (as a list). 433 return list(row)

434 435

436 - def _minimise_setup_atomic_pos(self, sim_index=None):

437 """Set up the atomic position data structures for optimisation using PCSs and PREs as base data sets. 438 439 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 440 @type sim_index: None or int 441 @return: The atomic positions (the first index is the spins, the second is the structures, and the third is the atomic coordinates) and the paramagnetic centre. 442 @rtype: numpy rank-3 array, numpy rank-1 array. 443 """ 444 445 # Initialise. 446 atomic_pos = [] 447 448 # Store the atomic positions. 449 for spin, spin_id in spin_loop(selection=self._domain_moving(), return_id=True): 450 # Skip deselected spins. 451 if not spin.select: 452 continue 453 454 # Only use spins with alignment/paramagnetic data. 455 if not hasattr(spin, 'pcs') and not hasattr(spin, 'pre'): 456 continue 457 458 # A single atomic position. 459 if spin.pos.shape == (3,): 460 atomic_pos.append(spin.pos) 461 462 # A single model (rank-2 array of a single position). 463 elif spin.pos.shape == (1, 3): 464 atomic_pos.append(spin.pos[0]) 465 466 # Average multiple atomic positions. 467 else: 468 # First throw a warning to tell the user what is happening. 469 if sim_index == None: 470 warn(RelaxWarning("Averaging the %s atomic positions for the PCS for the spin '%s'." % (len(spin.pos), spin_id))) 471 472 # The average position. 473 ave_pos = zeros(3, float64) 474 for i in range(len(spin.pos)): 475 ave_pos += spin.pos[i] 476 ave_pos = ave_pos / len(spin.pos) 477 478 # Store. 479 atomic_pos.append(ave_pos) 480 481 # Convert to numpy objects. 482 atomic_pos = array(atomic_pos, float64) 483 484 # The paramagnetic centre. 485 if not hasattr(cdp, 'paramagnetic_centre'): 486 paramag_centre = zeros(3, float64) 487 else: 488 paramag_centre = array(cdp.paramagnetic_centre) 489 490 # Return the data structures. 491 return atomic_pos, paramag_centre

492 493

494 - def _minimise_setup_pcs(self, sim_index=None):

495 """Set up the data structures for optimisation using PCSs as base data sets. 496 497 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 498 @type sim_index: None or int 499 @return: The assembled data structures for using PCSs as the base data for optimisation. These include: 500 - the PCS values. 501 - the unit vectors connecting the paramagnetic centre (the electron spin) to 502 - the PCS weight. 503 - the nuclear spin. 504 - the pseudocontact shift constants. 505 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array, numpy rank-1 array, numpy rank-1 array) 506 """ 507 508 # Data setup tests. 509 if not hasattr(cdp, 'paramagnetic_centre'): 510 raise RelaxError("The paramagnetic centre has not yet been specified.") 511 if not hasattr(cdp, 'temperature'): 512 raise RelaxError("The experimental temperatures have not been set.") 513 if not hasattr(cdp, 'spectrometer_frq'): 514 raise RelaxError("The spectrometer frequencies of the experiments have not been set.") 515 516 # Initialise. 517 pcs = [] 518 pcs_err = [] 519 pcs_weight = [] 520 temp = [] 521 frq = [] 522 523 # The PCS data. 524 for i in range(len(cdp.align_ids)): 525 # Alias the ID. 526 align_id = cdp.align_ids[i] 527 528 # Skip non-optimised data. 529 if not self._opt_uses_align_data(align_id): 530 continue 531 532 # Append empty arrays to the PCS structures. 533 pcs.append([]) 534 pcs_err.append([]) 535 pcs_weight.append([]) 536 537 # Get the temperature for the PCS constant. 538 if align_id in cdp.temperature: 539 temp.append(cdp.temperature[align_id]) 540 541 # The temperature must be given! 542 else: 543 raise RelaxError("The experimental temperature for the alignment ID '%s' has not been set." % align_id) 544 545 # Get the spectrometer frequency in Tesla units for the PCS constant. 546 if align_id in cdp.spectrometer_frq: 547 frq.append(cdp.spectrometer_frq[align_id] * 2.0 * pi / g1H) 548 549 # The frequency must be given! 550 else: 551 raise RelaxError("The spectrometer frequency for the alignment ID '%s' has not been set." % align_id) 552 553 # Spin loop over the domain. 554 j = 0 555 for spin in spin_loop(selection=self._domain_moving()): 556 # Skip deselected spins. 557 if not spin.select: 558 continue 559 560 # Skip spins without PCS data. 561 if not hasattr(spin, 'pcs'): 562 continue 563 564 # Append the PCSs to the list. 565 if align_id in spin.pcs.keys(): 566 if sim_index != None: 567 pcs[-1].append(spin.pcs_sim[align_id][sim_index]) 568 else: 569 pcs[-1].append(spin.pcs[align_id]) 570 else: 571 pcs[-1].append(None) 572 573 # Append the PCS errors. 574 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 575 pcs_err[-1].append(spin.pcs_err[align_id]) 576 else: 577 pcs_err[-1].append(None) 578 579 # Append the weight. 580 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight.keys(): 581 pcs_weight[-1].append(spin.pcs_weight[align_id]) 582 else: 583 pcs_weight[-1].append(1.0) 584 585 # Spin index. 586 j = j + 1 587 588 # Convert to numpy objects. 589 pcs = array(pcs, float64) 590 pcs_err = array(pcs_err, float64) 591 pcs_weight = array(pcs_weight, float64) 592 593 # Convert the PCS from ppm to no units. 594 pcs = pcs * 1e-6 595 pcs_err = pcs_err * 1e-6 596 597 # Return the data structures. 598 return pcs, pcs_err, pcs_weight, temp, frq

599 600

601 - def _minimise_setup_rdcs(self, sim_index=None):

602 """Set up the data structures for optimisation using RDCs as base data sets. 603 604 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 605 @type sim_index: None or int 606 @return: The assembled data structures for using RDCs as the base data for optimisation. These include: 607 - rdc, the RDC values. 608 - rdc_err, the RDC errors. 609 - rdc_weight, the RDC weights. 610 - vectors, the interatomic vectors. 611 - rdc_const, the dipolar constants. 612 - absolute, the absolute value flags (as 1's and 0's). 613 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array, numpy rank-3 array, numpy rank-2 array, numpy rank-2 array) 614 """ 615 616 # Initialise. 617 rdc = [] 618 rdc_err = [] 619 rdc_weight = [] 620 unit_vect = [] 621 rdc_const = [] 622 absolute = [] 623 624 # The unit vectors and RDC constants. 625 for interatom in interatomic_loop(selection1=self._domain_moving()): 626 # Get the spins. 627 spin1 = return_spin(interatom.spin_id1) 628 spin2 = return_spin(interatom.spin_id2) 629 630 # RDC checks. 631 if not self._check_rdcs(interatom, spin1, spin2): 632 continue 633 634 # A single unit vector. 635 if interatom.vector.shape == (3,): 636 unit_vect.append(interatom.vector) 637 638 # Average multiple unit vectors. 639 else: 640 # First throw a warning to tell the user what is happening. 641 if sim_index == None: 642 warn(RelaxWarning("Averaging the %s unit vectors for the RDC for the spin pair '%s' and '%s'." % (len(interatom.vector), interatom.spin_id1, interatom.spin_id2))) 643 644 # The average position. 645 ave_vector = zeros(3, float64) 646 for i in range(len(interatom.vector)): 647 ave_vector += interatom.vector[i] 648 649 # Store. 650 unit_vect.append(ave_vector) 651 652 # Normalise (to be sure). 653 unit_vect[-1] = unit_vect[-1] / norm(unit_vect[-1]) 654 655 # Gyromagnetic ratios. 656 g1 = return_gyromagnetic_ratio(spin1.isotope) 657 g2 = return_gyromagnetic_ratio(spin2.isotope) 658 659 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 660 rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r)) 661 662 # Fix the unit vector data structure. 663 num = None 664 for rdc_index in range(len(unit_vect)): 665 # Number of vectors. 666 if num == None: 667 if unit_vect[rdc_index] != None: 668 num = len(unit_vect[rdc_index]) 669 continue 670 671 # Check. 672 if unit_vect[rdc_index] != None and len(unit_vect[rdc_index]) != num: 673 raise RelaxError("The number of interatomic vectors for all no match:\n%s" % unit_vect) 674 675 # Missing unit vectors. 676 if num == None: 677 raise RelaxError("No interatomic vectors could be found.") 678 679 # Update None entries. 680 for i in range(len(unit_vect)): 681 if unit_vect[i] == None: 682 unit_vect[i] = [[None, None, None]]*num 683 684 # The RDC data. 685 for i in range(len(cdp.align_ids)): 686 # Alias the ID. 687 align_id = cdp.align_ids[i] 688 689 # Skip non-optimised data. 690 if not self._opt_uses_align_data(align_id): 691 continue 692 693 # Append empty arrays to the RDC structures. 694 rdc.append([]) 695 rdc_err.append([]) 696 rdc_weight.append([]) 697 absolute.append([]) 698 699 # Interatom loop over the domain. 700 for interatom in interatomic_loop(self._domain_moving()): 701 # Get the spins. 702 spin1 = return_spin(interatom.spin_id1) 703 spin2 = return_spin(interatom.spin_id2) 704 705 # Skip deselected spins. 706 if not spin1.select or not spin2.select: 707 continue 708 709 # Only use interatomic data containers with RDC and vector data. 710 if not hasattr(interatom, 'rdc') or not hasattr(interatom, 'vector'): 711 continue 712 713 # Defaults of None. 714 value = None 715 error = None 716 717 # Pseudo-atom set up. 718 if (hasattr(spin1, 'members') or hasattr(spin2, 'members')) and align_id in interatom.rdc.keys(): 719 raise RelaxError("Psuedo-atoms are currently not supported for the frame order analysis.") 720 721 # Normal set up. 722 elif align_id in interatom.rdc.keys(): 723 # The RDC. 724 if sim_index != None: 725 value = interatom.rdc_sim[align_id][sim_index] 726 else: 727 value = interatom.rdc[align_id] 728 729 # The error. 730 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 731 error = interatom.rdc_err[align_id] 732 733 # Append the RDCs to the list. 734 rdc[-1].append(value) 735 736 # Append the RDC errors. 737 rdc_err[-1].append(error) 738 739 # Append the weight. 740 if hasattr(interatom, 'rdc_weight') and align_id in interatom.rdc_weight.keys(): 741 rdc_weight[-1].append(interatom.rdc_weight[align_id]) 742 else: 743 rdc_weight[-1].append(1.0) 744 745 # Append the absolute value flag. 746 if hasattr(interatom, 'absolute_rdc') and align_id in interatom.absolute_rdc.keys(): 747 absolute[-1].append(interatom.absolute_rdc[align_id]) 748 else: 749 absolute[-1].append(False) 750 751 # Convert to numpy objects. 752 rdc = array(rdc, float64) 753 rdc_err = array(rdc_err, float64) 754 rdc_weight = array(rdc_weight, float64) 755 unit_vect = array(unit_vect, float64) 756 rdc_const = array(rdc_const, float64) 757 absolute = array(absolute, float64) 758 759 # Return the data structures. 760 return rdc, rdc_err, rdc_weight, unit_vect, rdc_const, absolute

761 762

763 - def _minimise_setup_tensors(self, sim_index=None):

764 """Set up the data structures for optimisation using alignment tensors as base data sets. 765 766 @keyword sim_index: The simulation index. This should be None if normal optimisation is desired. 767 @type sim_index: None or int 768 @return: The assembled data structures for using alignment tensors as the base data for optimisation. These include: 769 - full_tensors, the full tensors as concatenated arrays. 770 - full_err, the full tensor errors as concatenated arrays. 771 - full_in_ref_frame, the flags specifying if the tensor is the full or reduced tensor in the non-moving reference domain. 772 @rtype: tuple of 3 numpy nx5D, rank-1 arrays 773 """ 774 775 # Checks. 776 if not hasattr(cdp, 'ref_domain'): 777 raise RelaxError("The reference domain has not been set up.") 778 if not hasattr(cdp.align_tensors, 'reduction'): 779 raise RelaxError("The tensor reductions have not been specified.") 780 for i, tensor in self._tensor_loop(): 781 if not hasattr(tensor, 'domain'): 782 raise RelaxError("The domain that the '%s' tensor is attached to has not been set" % tensor.name) 783 784 # Initialise. 785 n = len(cdp.align_tensors.reduction) 786 full_tensors = zeros(n*5, float64) 787 full_err = ones(n*5, float64) * 1e-5 788 full_in_ref_frame = zeros(n, float64) 789 790 # Loop over the full tensors. 791 for i, tensor in self._tensor_loop(red=False): 792 # The full tensor. 793 full_tensors[5*i + 0] = tensor.Axx 794 full_tensors[5*i + 1] = tensor.Ayy 795 full_tensors[5*i + 2] = tensor.Axy 796 full_tensors[5*i + 3] = tensor.Axz 797 full_tensors[5*i + 4] = tensor.Ayz 798 799 # The full tensor corresponds to the frame of reference. 800 if cdp.ref_domain == tensor.domain: 801 full_in_ref_frame[i] = 1 802 803 # The full tensor errors. 804 if hasattr(tensor, 'Axx_err'): 805 full_err[5*i + 0] = tensor.Axx_err 806 full_err[5*i + 1] = tensor.Ayy_err 807 full_err[5*i + 2] = tensor.Axy_err 808 full_err[5*i + 3] = tensor.Axz_err 809 full_err[5*i + 4] = tensor.Ayz_err 810 811 # Return the data structures. 812 return full_tensors, full_err, full_in_ref_frame

813 814

815 - def _num_int_pts(self, num=200000):

816 """Set the number of integration points to use in the quasi-random Sobol' sequence. 817 818 @keyword num: The number of integration points. 819 @type num: int 820 """ 821 822 # Test if the current data pipe exists. 823 pipes.test() 824 825 # Store the value. 826 cdp.num_int_pts = num

827 828

829 - def _opt_uses_align_data(self, align_id=None):

830 """Determine if the PCS or RDC data for the given alignment ID is needed for optimisation. 831 832 @keyword align_id: The optional alignment ID string. 833 @type align_id: str 834 @return: True if alignment data is to be used for optimisation, False otherwise. 835 @rtype: bool 836 """ 837 838 # No alignment IDs. 839 if not hasattr(cdp, 'align_ids'): 840 return False 841 842 # Convert the align IDs to an array, or take all IDs. 843 if align_id: 844 align_ids = [align_id] 845 else: 846 align_ids = cdp.align_ids 847 848 # Check the PCS and RDC. 849 for align_id in align_ids: 850 if self._opt_uses_pcs(align_id) or self._opt_uses_rdc(align_id): 851 return True 852 853 # No alignment data is used for optimisation. 854 return False

855 856

857 - def _opt_uses_pcs(self, align_id):

858 """Determine if the PCS data for the given alignment ID is needed for optimisation. 859 860 @param align_id: The alignment ID string. 861 @type align_id: str 862 @return: True if the PCS data is to be used for optimisation, False otherwise. 863 @rtype: bool 864 """ 865 866 # No alignment IDs. 867 if not hasattr(cdp, 'pcs_ids'): 868 return False 869 870 # No PCS data for the alignment. 871 if align_id not in cdp.pcs_ids: 872 return False 873 874 # The PCS data is to be used for optimisation. 875 return True

876 877

878 - def _opt_uses_rdc(self, align_id):

879 """Determine if the RDC data for the given alignment ID is needed for optimisation. 880 881 @param align_id: The alignment ID string. 882 @type align_id: str 883 @return: True if the RDC data is to be used for optimisation, False otherwise. 884 @rtype: bool 885 """ 886 887 # No alignment IDs. 888 if not hasattr(cdp, 'rdc_ids'): 889 return False 890 891 # No RDC data for the alignment. 892 if align_id not in cdp.rdc_ids: 893 return False 894 895 # The RDC data is to be used for optimisation. 896 return True

897 898

899 - def _param_num(self):

900 """Determine the number of parameters in the model. 901 902 @return: The number of model parameters. 903 @rtype: int 904 """ 905 906 # Init. 907 num = 0 908 909 # Update the model if needed. 910 self._update_model() 911 912 # Determine the data type. 913 data_types = self._base_data_types() 914 915 # Average domain position translation. 916 if not self._translation_fixed(): 917 num += 3 918 919 # The pivot point. 920 if not self._pivot_fixed(): 921 num += 3 922 923 # Average domain position parameters. 924 if cdp.model in ['free rotor', 'iso cone, free rotor']: 925 num += 2 926 else: 927 num += 3 928 929 # Frame order eigenframe - the full frame. 930 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 931 num += 3 932 933 # Frame order eigenframe - the isotropic cone axis. 934 elif cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 935 num += 2 936 937 # Cone parameters - pseudo-elliptic cone parameters. 938 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 939 num += 2 940 941 # Cone parameters - single isotropic angle or order parameter. 942 elif cdp.model in ['iso cone', 'iso cone, torsionless', 'iso cone, free rotor']: 943 num += 1 944 945 # Cone parameters - torsion angle. 946 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 947 num += 1 948 949 # Return the number. 950 return num

951 952

953 - def _pdb_ave_pos(self, file=None, dir=None, force=False):

954 """Create a PDB file of the molecule with the moving domains shifted to the average position. 955 956 @keyword file: The name of the file for the average molecule structure. 957 @type file: str 958 @keyword dir: The name of the directory to place the PDB file into. 959 @type dir: str 960 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 961 @type force: bool 962 """ 963 964 # Printout. 965 subsection(file=sys.stdout, text="Creating a PDB file with the moving domains shifted to the average position.") 966 967 # Make a copy of the structural object (so as to preserve the original structure). 968 structure = deepcopy(cdp.structure) 969 970 # First rotate the moving domain to the average position. 971 R = zeros((3, 3), float64) 972 if hasattr(cdp, 'ave_pos_alpha'): 973 euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 974 else: 975 euler_to_R_zyz(0.0, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 976 if cdp.ave_pos_pivot == 'com': 977 origin = pipe_centre_of_mass(atom_id=self._domain_moving(), verbosity=0) 978 else: 979 origin = cdp.pivot 980 structure.rotate(R=R, origin=origin, atom_id=self._domain_moving()) 981 982 # Then translate the moving domain. 983 if not self._translation_fixed(): 984 structure.translate(T=[cdp.ave_pos_x, cdp.ave_pos_y, cdp.ave_pos_z], atom_id=self._domain_moving()) 985 986 # Write out the PDB file. 987 file = open_write_file(file_name=file, dir=dir, force=force) 988 structure.write_pdb(file=file) 989 file.close()

990 991

992 - def _pdb_distribution(self, file=None, dir=None, force=False):

993 """Create a PDB file of a distribution of positions coving the full dynamics of the moving domain. 994 995 @keyword file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model. 996 @type file: str 997 @keyword dir: The name of the directory to place the PDB file into. 998 @type dir: str 999 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 1000 @type force: bool 1001 """ 1002 1003 # Printout. 1004 subsection(file=sys.stdout, text="Creating a PDB file of a distribution of positions coving the full dynamics of the moving domain.")

1005 1006 1007

1008 - def _pdb_geometric_rep(self, file=None, dir=None, size=30.0, inc=36, force=False, neg_cone=True):

1009 """Create a PDB file containing a geometric object representing the frame order dynamics. 1010 1011 @keyword file: The name of the file of the PDB representation of the frame order dynamics to create. 1012 @type file: str 1013 @keyword dir: The name of the directory to place the PDB file into. 1014 @type dir: str 1015 @keyword size: The size of the geometric object in Angstroms. 1016 @type size: float 1017 @keyword inc: The number of increments for the filling of the cone objects. 1018 @type inc: int 1019 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 1020 @type force: bool 1021 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation. 1022 @type neg_cone: bool 1023 """ 1024 1025 # Printout. 1026 subsection(file=sys.stdout, text="Creating a PDB file containing a geometric object representing the frame order dynamics.") 1027 1028 # Monte Carlo simulation flag. 1029 sim = False 1030 num_sim = 0 1031 if hasattr(cdp, 'sim_number'): 1032 sim = True 1033 num_sim = cdp.sim_number 1034 1035 # The inversion matrix. 1036 inv_mat = -eye(3) 1037 1038 # Create the structural object. 1039 structure = Internal() 1040 1041 # Create model for the positive and negative images. 1042 model = structure.add_model(model=1) 1043 if neg_cone: 1044 model_neg = structure.add_model(model=2) 1045 1046 # The rotor object. 1047 if cdp.model in ['rotor', 'free rotor', 'iso cone', 'iso cone, free rotor', 'pseudo-ellipse']: 1048 # The rotor angle. 1049 if cdp.model in ['free rotor', 'iso cone, free rotor']: 1050 rotor_angle = pi 1051 else: 1052 rotor_angle = cdp.cone_sigma_max 1053 1054 # Generate the rotor axis. 1055 if cdp.model in ['rotor', 'free rotor', 'iso cone', 'iso cone, free rotor']: 1056 axis = zeros(3, float64) 1057 spherical_to_cartesian([1.0, cdp.axis_theta, cdp.axis_phi], axis) 1058 else: 1059 axes = zeros((3, 3), float64) 1060 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes) 1061 axis = axes[:, 2] 1062 1063 # Get the CoM of the entire molecule to use as the centre of the rotor. 1064 com = pipe_centre_of_mass(verbosity=0) 1065 1066 # Add the rotor object to the structure as a new molecule. 1067 rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis, axis_pt=cdp.pivot, centre=com, span=2e-9, blade_length=5e-10, staggered=False) 1068 1069 # FIXME: Temporary write out and exit. 1070 print("\nGenerating the PDB file.") 1071 pdb_file = open_write_file(file, dir, force=force) 1072 structure.write_pdb(pdb_file) 1073 pdb_file.close() 1074 return 1075 1076 # Create the molecule. 1077 structure.add_molecule(name='rest') 1078 1079 # Alias the molecules. 1080 mol = model.mol[0] 1081 if neg_cone: 1082 mol_neg = model_neg.mol[0] 1083 1084 1085 # The pivot point. 1086 ################## 1087 1088 # Add the pivot point. 1089 structure.add_atom(mol_name=cdp.model, pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=cdp.pivot, element='C') 1090 1091 1092 # The axes. 1093 ########### 1094 1095 # The spherical angles. 1096 if cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 1097 # Print out. 1098 print("\nGenerating the z-axis system.") 1099 1100 # The axis. 1101 axis = zeros(3, float64) 1102 spherical_to_cartesian([1.0, getattr(cdp, 'axis_theta'), getattr(cdp, 'axis_phi')], axis) 1103 print(("Central axis: %s." % axis)) 1104 1105 # Rotations and inversions. 1106 axis_pos = axis 1107 axis_neg = dot(inv_mat, axis) 1108 1109 # Simulation central axis. 1110 axis_sim_pos = None 1111 axis_sim_neg = None 1112 if sim: 1113 # Init. 1114 axis_sim = zeros((cdp.sim_number, 3), float64) 1115 1116 # Fill the structure. 1117 for i in range(cdp.sim_number): 1118 spherical_to_cartesian([1.0, getattr(cdp, 'axis_theta_sim')[i], getattr(cdp, 'axis_phi_sim')[i]], axis_sim[i]) 1119 1120 # Inversion. 1121 axis_sim_pos = axis_sim 1122 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos))) 1123 1124 # Generate the axis vectors. 1125 print("\nGenerating the axis vectors.") 1126 res_num = geometric.generate_vector_residues(mol=mol, vector=axis_pos, atom_name='z-ax', res_name_vect='AXE', sim_vectors=axis_sim_pos, res_num=2, origin=cdp.pivot, scale=size) 1127 1128 # The negative. 1129 if neg_cone: 1130 res_num = geometric.generate_vector_residues(mol=mol_neg, vector=axis_neg, atom_name='z-ax', res_name_vect='AXE', sim_vectors=axis_sim_neg, res_num=2, origin=cdp.pivot, scale=size) 1131 1132 # The full axis system. 1133 else: 1134 # Print out. 1135 print("\nGenerating the full axis system.") 1136 1137 # The axis system. 1138 axes = zeros((3, 3), float64) 1139 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes) 1140 print(("Axis system:\n%s" % axes)) 1141 1142 # Rotations and inversions. 1143 axes_pos = axes 1144 axes_neg = dot(inv_mat, axes) 1145 1146 # Simulations 1147 axes_sim_pos = None 1148 axes_sim_neg = None 1149 if sim: 1150 # Init. 1151 axes_sim_pos = zeros((cdp.sim_number, 3, 3), float64) 1152 axes_sim_neg = zeros((cdp.sim_number, 3, 3), float64) 1153 1154 # Fill the structure. 1155 for i in range(cdp.sim_number): 1156 # The positive system. 1157 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_pos[i]) 1158 1159 # The negative system. 1160 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_neg[i]) 1161 axes_sim_neg[i] = dot(inv_mat, axes_sim_neg[i]) 1162 1163 # Generate the axis vectors. 1164 print("\nGenerating the axis vectors.") 1165 label = ['x', 'y', 'z'] 1166 for j in range(len(label)): 1167 # The simulation data. 1168 axis_sim_pos = None 1169 axis_sim_neg = None 1170 if sim: 1171 axis_sim_pos = axes_sim_pos[:,:, j] 1172 axis_sim_neg = axes_sim_neg[:,:, j] 1173 1174 # The vectors. 1175 res_num = geometric.generate_vector_residues(mol=mol, vector=axes_pos[:, j], atom_name='%s-ax'%label[j], res_name_vect='AXE', sim_vectors=axis_sim_pos, res_num=2, origin=cdp.pivot, scale=size) 1176 if neg_cone: 1177 res_num = geometric.generate_vector_residues(mol=mol_neg, vector=axes_neg[:, j], atom_name='%s-ax'%label[j], res_name_vect='AXE', sim_vectors=axis_sim_neg, res_num=2, origin=cdp.pivot, scale=size) 1178 1179 1180 # The cone object. 1181 ################## 1182 1183 # Skip models missing a cone. 1184 if cdp.model not in ['rotor', 'free rotor']: 1185 # The rotation matrix (rotation from the z-axis to the cone axis). 1186 if cdp.model not in ['iso cone', 'iso cone, torsionless', 'iso cone, free rotor']: 1187 R = axes 1188 else: 1189 R = zeros((3, 3), float64) 1190 two_vect_to_R(array([0, 0, 1], float64), axis, R) 1191 1192 # Average position rotation. 1193 R_pos = R 1194 R_neg = dot(inv_mat, R) 1195 1196 # The pseudo-ellipse cone object. 1197 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 1198 cone = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y) 1199 1200 # The isotropic cone object. 1201 else: 1202 # The angle. 1203 if hasattr(cdp, 'cone_theta'): 1204 cone_theta = cdp.cone_theta 1205 elif hasattr(cdp, 'cone_s1'): 1206 cone_theta = order_parameters.iso_cone_S_to_theta(cdp.cone_s1) 1207 1208 # The object. 1209 cone = Iso_cone(cone_theta) 1210 1211 # Create the positive and negative cones. 1212 geometric.create_cone_pdb(mol=mol, cone=cone, start_res=mol.res_num[-1]+1, apex=cdp.pivot, R=R_pos, inc=inc, distribution='regular', axis_flag=False) 1213 1214 # The negative. 1215 if neg_cone: 1216 geometric.create_cone_pdb(mol=mol_neg, cone=cone, start_res=mol_neg.res_num[-1]+1, apex=cdp.pivot, R=R_neg, inc=inc, distribution='regular', axis_flag=False) 1217 1218 1219 # Create the PDB file. 1220 ###################### 1221 1222 # Print out. 1223 print("\nGenerating the PDB file.") 1224 1225 # Write the file. 1226 pdb_file = open_write_file(file, dir, force=force) 1227 structure.write_pdb(pdb_file) 1228 pdb_file.close()

1229 1230

1231 - def _pdb_model(self, ave_pos_file="ave_pos.pdb", rep_file="frame_order.pdb", dist_file="domain_distribution.pdb", dir=None, size=30.0, inc=36, force=False, neg_cone=True):

1232 """Create 3 different PDB files for representing the frame order dynamics of the system. 1233 1234 @keyword ave_pos_file: The name of the file for the average molecule structure. 1235 @type ave_pos_file: str 1236 @keyword rep_file: The name of the file of the PDB representation of the frame order dynamics to create. 1237 @type rep_file: str 1238 @keyword dist_file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model. 1239 @type dist_file: str 1240 @keyword dir: The name of the directory to place the PDB file into. 1241 @type dir: str 1242 @keyword size: The size of the geometric object in Angstroms. 1243 @type size: float 1244 @keyword inc: The number of increments for the filling of the cone objects. 1245 @type inc: int 1246 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 1247 @type force: bool 1248 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation. 1249 @type neg_cone: bool 1250 """ 1251 1252 # Test if the current data pipe exists. 1253 pipes.test() 1254 1255 # Create the average position structure. 1256 self._pdb_ave_pos(file=ave_pos_file, dir=dir, force=force) 1257 1258 # Nothing more to do for the rigid model. 1259 if cdp.model == 'rigid': 1260 return 1261 1262 # Create the geometric representation. 1263 self._pdb_geometric_rep(file=rep_file, dir=dir, size=size, inc=inc, force=force, neg_cone=neg_cone) 1264 1265 # Create the distribution. 1266 self._pdb_distribution(file=dist_file, dir=dir, force=force)

1267 1268

1269 - def _pivot(self, pivot=None, fix=None):

1270 """Set the pivot point for the 2 body motion. 1271 1272 @keyword pivot: The pivot point of the two bodies (domains, etc.) in the structural coordinate system. 1273 @type pivot: list of num 1274 @keyword fix: A flag specifying if the pivot point should be fixed during optimisation. 1275 @type fix: bool 1276 """ 1277 1278 # Test if the current data pipe exists. 1279 pipes.test() 1280 1281 # Set the pivot point and fixed flag. 1282 cdp.pivot = pivot 1283 cdp.pivot_fixed = fix 1284 1285 # Convert to floats. 1286 for i in range(3): 1287 cdp.pivot[i] = float(cdp.pivot[i])

1288 1289

1290 - def _pivot_fixed(self):

1291 """Determine if the pivot is fixed or not. 1292 1293 @return: The answer to the question. 1294 @rtype: bool 1295 """ 1296 1297 # A pivot point is not supported by the model. 1298 if cdp.model in ['rigid']: 1299 return True 1300 1301 # The PCS is loaded. 1302 if 'pcs' in self._base_data_types(): 1303 # The fixed flag is not set. 1304 if hasattr(cdp, 'pivot_fixed') and not cdp.pivot_fixed: 1305 return False 1306 1307 # The point is fixed. 1308 return True

1309 1310

1311 - def _quad_int(self, flag=False):

1312 """Turn the high precision Scipy quadratic numerical integration on or off. 1313 1314 @keyword flag: The flag which if True will perform high precision numerical integration via the scipy.integrate quad(), dblquad() and tplquad() integration methods rather than the rough quasi-random numerical integration. 1315 @type flag: bool 1316 """ 1317 1318 # Test if the current data pipe exists. 1319 pipes.test() 1320 1321 # Store the flag. 1322 cdp.quad_int = flag

1323 1324

1325 - def _ref_domain(self, ref=None):

1326 """Set the reference domain for the frame order, multi-domain models. 1327 1328 @param ref: The reference domain. 1329 @type ref: str 1330 """ 1331 1332 # Test if the current data pipe exists. 1333 pipes.test() 1334 1335 # Check that the domain is defined. 1336 if not hasattr(cdp, 'domain') or ref not in list(cdp.domain.keys()): 1337 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % ref) 1338 1339 # Test if the reference domain exists. 1340 exists = False 1341 for tensor_cont in cdp.align_tensors: 1342 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref: 1343 exists = True 1344 if not exists: 1345 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 1346 1347 # Set the reference domain. 1348 cdp.ref_domain = ref 1349 1350 # Update the model. 1351 if hasattr(cdp, 'model'): 1352 self._update_model()

1353 1354

1355 - def _select_model(self, model=None):

1356 """Select the Frame Order model. 1357 1358 @param model: The Frame Order model. This can be one of 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor'. 1359 @type model: str 1360 """ 1361 1362 # Test if the current data pipe exists. 1363 pipes.test() 1364 1365 # Test if the model name exists. 1366 if not model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor']: 1367 raise RelaxError("The model name " + repr(model) + " is invalid.") 1368 1369 # Set the model 1370 cdp.model = model 1371 1372 # Initialise the list of model parameters. 1373 cdp.params = [] 1374 1375 # Set the integration method if needed. 1376 if not hasattr(cdp, 'quad_int'): 1377 # Scipy quadratic numerical integration. 1378 if cdp.model in []: 1379 cdp.quad_int = True 1380 1381 # Quasi-random numerical integration. 1382 else: 1383 cdp.quad_int = False 1384 1385 # Update the model. 1386 self._update_model()

1387 1388

1389 - def _store_bc_data(self, target_fn):

1390 """Store the back-calculated data. 1391 1392 @param target_fn: The frame-order target function class. 1393 @type target_fn: class instance 1394 """ 1395 1396 # Loop over the reduced tensors. 1397 for i, tensor in self._tensor_loop(red=True): 1398 # Store the values. 1399 tensor.set(param='Axx', value=target_fn.A_5D_bc[5*i + 0]) 1400 tensor.set(param='Ayy', value=target_fn.A_5D_bc[5*i + 1]) 1401 tensor.set(param='Axy', value=target_fn.A_5D_bc[5*i + 2]) 1402 tensor.set(param='Axz', value=target_fn.A_5D_bc[5*i + 3]) 1403 tensor.set(param='Ayz', value=target_fn.A_5D_bc[5*i + 4]) 1404 1405 # The RDC data. 1406 for i in range(len(cdp.align_ids)): 1407 # The alignment ID. 1408 align_id = cdp.align_ids[i] 1409 1410 # Data flags 1411 rdc_flag = False 1412 if hasattr(cdp, 'rdc_ids') and align_id in cdp.rdc_ids: 1413 rdc_flag = True 1414 pcs_flag = False 1415 if hasattr(cdp, 'pcs_ids') and align_id in cdp.pcs_ids: 1416 pcs_flag = True 1417 1418 # Spin loop over the domain. 1419 pcs_index = 0 1420 for spin in spin_loop(self._domain_moving()): 1421 # Skip deselected spins. 1422 if not spin.select: 1423 continue 1424 1425 # Spins with PCS data. 1426 if pcs_flag and hasattr(spin, 'pcs'): 1427 # Initialise the data structure. 1428 if not hasattr(spin, 'pcs_bc'): 1429 spin.pcs_bc = {} 1430 1431 # Store the back-calculated value (in ppm). 1432 spin.pcs_bc[align_id] = target_fn.pcs_theta[i, pcs_index] * 1e6 1433 1434 # Increment the index. 1435 pcs_index += 1 1436 1437 # Interatomic data container loop. 1438 rdc_index = 0 1439 for interatom in interatomic_loop(self._domain_moving()): 1440 # Get the spins. 1441 spin1 = return_spin(interatom.spin_id1) 1442 spin2 = return_spin(interatom.spin_id2) 1443 1444 # RDC checks. 1445 if not self._check_rdcs(interatom, spin1, spin2): 1446 continue 1447 1448 # Initialise the data structure. 1449 if not hasattr(interatom, 'rdc_bc'): 1450 interatom.rdc_bc = {} 1451 1452 # Store the back-calculated value. 1453 interatom.rdc_bc[align_id] = target_fn.rdc_theta[i, rdc_index] 1454 1455 # Increment the index. 1456 rdc_index += 1

1457 1458

1459 - def _target_fn_setup(self, sim_index=None, scaling=True):

1460 """Initialise the target function for optimisation or direct calculation. 1461 1462 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 1463 @type sim_index: None or int 1464 @param scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 1465 @type scaling: bool 1466 """ 1467 1468 # Check for the average domain displacement information. 1469 if not hasattr(cdp, 'ave_pos_pivot') or not hasattr(cdp, 'ave_pos_translation'): 1470 raise RelaxError("The mechanics of the displacement to the average domain position have not been set up.") 1471 1472 # Assemble the parameter vector. 1473 param_vector = self._assemble_param_vector(sim_index=sim_index) 1474 1475 # Determine if alignment tensors or RDCs are to be used. 1476 data_types = self._base_data_types() 1477 1478 # Diagonal scaling. 1479 scaling_matrix = None 1480 if len(param_vector): 1481 scaling_matrix = self._assemble_scaling_matrix(data_types=data_types, scaling=scaling) 1482 param_vector = dot(inv(scaling_matrix), param_vector) 1483 1484 # Get the data structures for optimisation using the tensors as base data sets. 1485 full_tensors, full_tensor_err, full_in_ref_frame = self._minimise_setup_tensors(sim_index) 1486 1487 # Get the data structures for optimisation using PCSs as base data sets. 1488 pcs, pcs_err, pcs_weight, temp, frq = None, None, None, None, None 1489 if 'pcs' in data_types: 1490 pcs, pcs_err, pcs_weight, temp, frq = self._minimise_setup_pcs(sim_index=sim_index) 1491 1492 # Get the data structures for optimisation using RDCs as base data sets. 1493 rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, absolute_rdc = None, None, None, None, None, None 1494 if 'rdc' in data_types: 1495 rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, absolute_rdc = self._minimise_setup_rdcs(sim_index=sim_index) 1496 1497 # Data checks. 1498 if pcs != None and not len(pcs): 1499 raise RelaxNoPCSError 1500 if rdcs != None and not len(rdcs): 1501 raise RelaxNoRDCError 1502 1503 # Get the atomic_positions. 1504 atomic_pos, paramag_centre = None, None 1505 if 'pcs' in data_types or 'pre' in data_types: 1506 atomic_pos, paramag_centre = self._minimise_setup_atomic_pos(sim_index=sim_index) 1507 1508 # Average domain translation. 1509 translation_opt = False 1510 if not self._translation_fixed(): 1511 translation_opt = True 1512 1513 # The fixed pivot point. 1514 pivot = None 1515 if hasattr(cdp, 'pivot'): 1516 pivot = cdp.pivot 1517 1518 # Pivot optimisation. 1519 pivot_opt = True 1520 if self._pivot_fixed(): 1521 pivot_opt = False 1522 1523 # The number of integration points. 1524 if not hasattr(cdp, 'num_int_pts'): 1525 cdp.num_int_pts = 200000 1526 1527 # The centre of mass of the moving domain - to use as the centroid for the average domain position rotation. 1528 if cdp.ave_pos_pivot == 'com': 1529 com = pipe_centre_of_mass(atom_id=self._domain_moving(), verbosity=0) 1530 ave_pos_piv_sync = False 1531 1532 # The centre of mass of the moving domain - to use as the centroid for the average domain position rotation. 1533 if cdp.ave_pos_pivot == 'motional': 1534 com = pivot 1535 ave_pos_piv_sync = True 1536 1537 # Print outs. 1538 if sim_index == None: 1539 if cdp.model != 'rigid': 1540 if cdp.quad_int: 1541 sys.stdout.write("Numerical integration via Scipy quadratic integration.\n") 1542 else: 1543 sys.stdout.write("Numerical integration via the quasi-random Sobol' sequence.\n") 1544 sys.stdout.write("Number of integration points: %s\n" % cdp.num_int_pts) 1545 base_data = [] 1546 if rdcs != None and len(rdcs): 1547 base_data.append("RDCs") 1548 if pcs != None and len(pcs): 1549 base_data.append("PCSs") 1550 sys.stdout.write("Base data: %s\n" % repr(base_data)) 1551 1552 # Set up the optimisation function. 1553 target = frame_order.Frame_order(model=cdp.model, init_params=param_vector, full_tensors=full_tensors, full_in_ref_frame=full_in_ref_frame, rdcs=rdcs, rdc_errors=rdc_err, rdc_weights=rdc_weight, rdc_vect=rdc_vect, dip_const=rdc_const, pcs=pcs, pcs_errors=pcs_err, pcs_weights=pcs_weight, atomic_pos=atomic_pos, temp=temp, frq=frq, paramag_centre=paramag_centre, scaling_matrix=scaling_matrix, ave_pos_pivot=com, ave_pos_piv_sync=ave_pos_piv_sync, translation_opt=translation_opt, pivot=pivot, pivot_opt=pivot_opt, num_int_pts=cdp.num_int_pts, quad_int=cdp.quad_int) 1554 1555 # Return the data. 1556 return target, param_vector, data_types, scaling_matrix

1557 1558

1559 - def _tensor_loop(self, red=False):

1560 """Generator method for looping over the full or reduced tensors. 1561 1562 @keyword red: A flag which if True causes the reduced tensors to be returned, and if False 1563 the full tensors are returned. 1564 @type red: bool 1565 @return: The tensor index and the tensor. 1566 @rtype: (int, AlignTensorData instance) 1567 """ 1568 1569 # Number of tensor pairs. 1570 n = len(cdp.align_tensors.reduction) 1571 1572 # Alias. 1573 data = cdp.align_tensors 1574 list = data.reduction 1575 1576 # Full or reduced index. 1577 if red: 1578 index = 1 1579 else: 1580 index = 0 1581 1582 # Loop over the reduction list. 1583 for i in range(n): 1584 yield i, data[list[i][index]]

1585 1586

1587 - def _translation_fixed(self):

1588 """Is the translation of the average domain position fixed? 1589 1590 @return: The answer to the question. 1591 @rtype: bool 1592 """ 1593 1594 # PCS data must be present. 1595 if 'pcs' in self._base_data_types(): 1596 # The fixed flag is not set. 1597 if hasattr(cdp, 'ave_pos_translation') and cdp.ave_pos_translation: 1598 return False 1599 1600 # Default to being fixed. 1601 return True

1602 1603

1604 - def _update_model(self):

1605 """Update the model parameters as necessary.""" 1606 1607 # Re-initialise the list of model parameters. 1608 cdp.params = [] 1609 1610 # The pivot parameters. 1611 if not self._pivot_fixed(): 1612 cdp.params.append('pivot_x') 1613 cdp.params.append('pivot_y') 1614 cdp.params.append('pivot_z') 1615 1616 # The average domain position translation parameters. 1617 if not self._translation_fixed(): 1618 cdp.params.append('ave_pos_x') 1619 cdp.params.append('ave_pos_y') 1620 cdp.params.append('ave_pos_z') 1621 1622 # The tensor rotation, or average domain position. 1623 if cdp.model not in ['free rotor', 'iso cone, free rotor']: 1624 cdp.params.append('ave_pos_alpha') 1625 cdp.params.append('ave_pos_beta') 1626 cdp.params.append('ave_pos_gamma') 1627 1628 # Frame order eigenframe - the full frame. 1629 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 1630 cdp.params.append('eigen_alpha') 1631 cdp.params.append('eigen_beta') 1632 cdp.params.append('eigen_gamma') 1633 1634 # Frame order eigenframe - the isotropic cone axis. 1635 elif cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 1636 cdp.params.append('axis_theta') 1637 cdp.params.append('axis_phi') 1638 1639 # Cone parameters - pseudo-elliptic cone parameters. 1640 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 1641 cdp.params.append('cone_theta_x') 1642 cdp.params.append('cone_theta_y') 1643 1644 # Cone parameters - single isotropic angle or order parameter. 1645 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 1646 cdp.params.append('cone_theta') 1647 elif cdp.model in ['iso cone, free rotor']: 1648 cdp.params.append('cone_s1') 1649 1650 # Cone parameters - torsion angle. 1651 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 1652 cdp.params.append('cone_sigma_max') 1653 1654 # Initialise the parameters in the current data pipe. 1655 for param in cdp.params: 1656 if not param in ['pivot_x', 'pivot_y', 'pivot_z'] and not hasattr(cdp, param): 1657 setattr(cdp, param, 0.0)

1658 1659

1660 - def _unpack_opt_results(self, results, scaling=False, scaling_matrix=None, sim_index=None):

1661 """Unpack and store the Frame Order optimisation results. 1662 1663 @param results: The results tuple returned by the minfx generic_minimise() function. 1664 @type results: tuple 1665 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 1666 @type scaling: bool 1667 @keyword scaling_matrix: The scaling matrix. 1668 @type scaling_matrix: numpy rank-2 array 1669 @keyword sim_index: The index of the simulation to optimise. This should be None for normal optimisation. 1670 @type sim_index: None or int 1671 """ 1672 1673 # Disassemble the results. 1674 if len(results) == 4: # Grid search. 1675 param_vector, func, iter_count, warning = results 1676 f_count = iter_count 1677 g_count = 0.0 1678 h_count = 0.0 1679 else: 1680 param_vector, func, iter_count, f_count, g_count, h_count, warning = results 1681 1682 # Catch infinite chi-squared values. 1683 if isInf(func): 1684 raise RelaxInfError('chi-squared') 1685 1686 # Catch chi-squared values of NaN. 1687 if isNaN(func): 1688 raise RelaxNaNError('chi-squared') 1689 1690 # Scaling. 1691 if scaling: 1692 param_vector = dot(scaling_matrix, param_vector) 1693 1694 # Pivot point. 1695 if not self._pivot_fixed(): 1696 # Store the pivot. 1697 cdp.pivot = param_vector[:3] 1698 1699 # Then remove it from the params. 1700 param_vector = param_vector[3:] 1701 1702 # Average domain position translation. 1703 ave_pos_x, ave_pos_y, ave_pos_z = None, None, None 1704 if not self._translation_fixed(): 1705 ave_pos_x = param_vector[0] 1706 ave_pos_y = param_vector[1] 1707 ave_pos_z = param_vector[2] 1708 1709 # Remove the parameters from the array. 1710 param_vector = param_vector[3:] 1711 1712 # Unpack the parameters. 1713 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = None, None, None 1714 eigen_alpha, eigen_beta, eigen_gamma = None, None, None 1715 axis_theta, axis_phi = None, None 1716 cone_theta_x, cone_theta_y = None, None 1717 cone_theta = None 1718 cone_s1 = None 1719 cone_sigma_max = None 1720 if cdp.model == 'pseudo-ellipse': 1721 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = param_vector 1722 elif cdp.model in ['pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 1723 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = param_vector 1724 elif cdp.model == 'iso cone': 1725 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, cone_sigma_max = param_vector 1726 elif cdp.model in ['iso cone, torsionless']: 1727 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = param_vector 1728 elif cdp.model in ['iso cone, free rotor']: 1729 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = param_vector 1730 ave_pos_alpha = None 1731 elif cdp.model == 'line': 1732 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector 1733 elif cdp.model in ['line, torsionless', 'line, free rotor']: 1734 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector 1735 elif cdp.model in ['rotor']: 1736 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_sigma_max = param_vector 1737 elif cdp.model in ['free rotor']: 1738 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = param_vector 1739 ave_pos_alpha = None 1740 elif cdp.model == 'rigid': 1741 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = param_vector 1742 1743 # Monte Carlo simulation data structures. 1744 if sim_index != None: 1745 # Average position parameters. 1746 if ave_pos_x != None: 1747 cdp.ave_pos_x_sim[sim_index] = ave_pos_x 1748 if ave_pos_y != None: 1749 cdp.ave_pos_y_sim[sim_index] = ave_pos_y 1750 if ave_pos_z != None: 1751 cdp.ave_pos_z_sim[sim_index] = ave_pos_z 1752 if ave_pos_alpha != None: 1753 cdp.ave_pos_alpha_sim[sim_index] = wrap_angles(ave_pos_alpha, cdp.ave_pos_alpha-pi, cdp.ave_pos_alpha+pi) 1754 if ave_pos_beta != None: 1755 cdp.ave_pos_beta_sim[sim_index] = wrap_angles(ave_pos_beta, cdp.ave_pos_beta-pi, cdp.ave_pos_beta+pi) 1756 if ave_pos_gamma != None: 1757 cdp.ave_pos_gamma_sim[sim_index] = wrap_angles(ave_pos_gamma, cdp.ave_pos_gamma-pi, cdp.ave_pos_gamma+pi) 1758 1759 # Eigenframe parameters. 1760 if eigen_alpha != None: 1761 cdp.eigen_alpha_sim[sim_index] = wrap_angles(eigen_alpha, cdp.eigen_alpha-pi, cdp.eigen_alpha+pi) 1762 if eigen_beta != None: 1763 cdp.eigen_beta_sim[sim_index] = wrap_angles(eigen_beta, cdp.eigen_beta-pi, cdp.eigen_beta+pi) 1764 if eigen_gamma != None: 1765 cdp.eigen_gamma_sim[sim_index] = wrap_angles(eigen_gamma, cdp.eigen_gamma-pi, cdp.eigen_gamma+pi) 1766 if axis_theta != None: 1767 cdp.axis_theta_sim[sim_index] = wrap_angles(axis_theta, cdp.axis_theta-pi, cdp.axis_theta+pi) 1768 if axis_phi != None: 1769 cdp.axis_phi_sim[sim_index] = wrap_angles(axis_phi, cdp.axis_phi-pi, cdp.axis_phi+pi) 1770 1771 # Cone parameters. 1772 if cone_theta != None: 1773 cdp.cone_theta_sim[sim_index] = cone_theta 1774 if cone_s1 != None: 1775 cdp.cone_s1_sim[sim_index] = cone_s1 1776 cdp.cone_theta_sim[sim_index] = order_parameters.iso_cone_S_to_theta(cone_s1) 1777 if cone_theta_x != None: 1778 cdp.cone_theta_x_sim[sim_index] = cone_theta_x 1779 if cone_theta_y != None: 1780 cdp.cone_theta_y_sim[sim_index] = cone_theta_y 1781 if cone_sigma_max != None: 1782 cdp.cone_sigma_max_sim[sim_index] = abs(cone_sigma_max) 1783 1784 # Optimisation stats. 1785 cdp.chi2_sim[sim_index] = func 1786 cdp.iter_sim[sim_index] = iter_count 1787 cdp.f_count_sim[sim_index] = f_count 1788 cdp.g_count_sim[sim_index] = g_count 1789 cdp.h_count_sim[sim_index] = h_count 1790 cdp.warning_sim[sim_index] = warning 1791 1792 # Normal data structures. 1793 else: 1794 # Average position parameters. 1795 if ave_pos_x != None: 1796 cdp.ave_pos_x = ave_pos_x 1797 if ave_pos_y != None: 1798 cdp.ave_pos_y = ave_pos_y 1799 if ave_pos_z != None: 1800 cdp.ave_pos_z = ave_pos_z 1801 if ave_pos_alpha != None: 1802 cdp.ave_pos_alpha = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi) 1803 if ave_pos_beta != None: 1804 cdp.ave_pos_beta = wrap_angles(ave_pos_beta, 0.0, 2.0*pi) 1805 if ave_pos_gamma != None: 1806 cdp.ave_pos_gamma = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi) 1807 1808 # Eigenframe parameters. 1809 if eigen_alpha != None: 1810 cdp.eigen_alpha = wrap_angles(eigen_alpha, 0.0, 2.0*pi) 1811 if eigen_beta != None: 1812 cdp.eigen_beta = wrap_angles(eigen_beta, 0.0, 2.0*pi) 1813 if eigen_gamma != None: 1814 cdp.eigen_gamma = wrap_angles(eigen_gamma, 0.0, 2.0*pi) 1815 if axis_theta != None: 1816 cdp.axis_theta = wrap_angles(axis_theta, 0.0, 2.0*pi) 1817 if axis_phi != None: 1818 cdp.axis_phi = wrap_angles(axis_phi, 0.0, 2.0*pi) 1819 1820 # Cone parameters. 1821 if cone_theta != None: 1822 cdp.cone_theta = cone_theta 1823 if cone_s1 != None: 1824 cdp.cone_s1 = cone_s1 1825 cdp.cone_theta = order_parameters.iso_cone_S_to_theta(cone_s1) 1826 if cone_theta_x != None: 1827 cdp.cone_theta_x = cone_theta_x 1828 if cone_theta_y != None: 1829 cdp.cone_theta_y = cone_theta_y 1830 if cone_sigma_max != None: 1831 cdp.cone_sigma_max = abs(cone_sigma_max) 1832 1833 # Optimisation stats. 1834 cdp.chi2 = func 1835 cdp.iter = iter_count 1836 cdp.f_count = f_count 1837 cdp.g_count = g_count 1838 cdp.h_count = h_count 1839 cdp.warning = warning

1840 1841

1842 - def base_data_loop(self):

1843 """Generator method for looping over the base data - alignment tensors, RDCs, PCSs. 1844 1845 This loop yields the following: 1846 1847 - The RDC identification data for the interatomic data container and alignment. 1848 - The PCS identification data for the spin data container and alignment. 1849 1850 @return: The base data type ('rdc' or 'pcs'), the spin or interatomic data container information (either one or two spin IDs), and the alignment ID string. 1851 @rtype: list of str 1852 """ 1853 1854 # Loop over the interatomic data containers for the moving domain (for the RDC data). 1855 for interatom in interatomic_loop(selection1=self._domain_moving()): 1856 # Get the spins. 1857 spin1 = return_spin(interatom.spin_id1) 1858 spin2 = return_spin(interatom.spin_id2) 1859 1860 # RDC checks. 1861 if not self._check_rdcs(interatom, spin1, spin2): 1862 continue 1863 1864 # Loop over the alignment IDs. 1865 for align_id in cdp.rdc_ids: 1866 # Yield the info set. 1867 yield ['rdc', interatom.spin_id1, interatom.spin_id2, align_id] 1868 1869 # Loop over the spin containers for the moving domain (for the PCS data). 1870 for spin, spin_id in spin_loop(selection=self._domain_moving(), return_id=True): 1871 # Skip deselected spins. 1872 if not spin.select: 1873 continue 1874 1875 # No PCS, so skip. 1876 if not hasattr(spin, 'pcs'): 1877 continue 1878 1879 # Loop over the alignment IDs. 1880 for align_id in cdp.pcs_ids: 1881 # Yield the info set. 1882 yield ['pcs', spin_id, align_id]

1883 1884

1885 - def calculate(self, spin_id=None, verbosity=1, sim_index=None):

1886 """Calculate the chi-squared value for the current parameter values. 1887 1888 @keyword spin_id: The spin identification string (unused). 1889 @type spin_id: None 1890 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 1891 @type verbosity: int 1892 @keyword sim_index: The optional MC simulation index (unused). 1893 @type sim_index: None or int 1894 """ 1895 1896 # Set up the target function for direct calculation. 1897 model, param_vector, data_types, scaling_matrix = self._target_fn_setup(sim_index=sim_index) 1898 1899 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 1900 chi2 = model.func(param_vector) 1901 1902 # Set the chi2. 1903 cdp.chi2 = chi2 1904 1905 # Store the back-calculated data. 1906 self._store_bc_data(model) 1907 1908 # Printout. 1909 print("Chi2: %s" % chi2)

1910 1911

1912 - def create_mc_data(self, data_id=None):

1913 """Create the Monte Carlo data by back calculating the reduced tensor data. 1914 1915 @keyword data_id: The data set as yielded by the base_data_loop() generator method. 1916 @type data_id: list of str 1917 @return: The Monte Carlo simulation data. 1918 @rtype: list of floats 1919 """ 1920 1921 # Initialise the MC data structure. 1922 mc_data = [] 1923 1924 # The RDC data. 1925 if data_id[0] == 'rdc': 1926 # Unpack the set. 1927 data_type, spin_id1, spin_id2, align_id = data_id 1928 1929 # Get the interatomic data container. 1930 interatom = return_interatom(spin_id1, spin_id2) 1931 1932 # Does back-calculated data exist? 1933 if not hasattr(interatom, 'rdc_bc'): 1934 self.calculate() 1935 1936 # The data. 1937 if not hasattr(interatom, 'rdc_bc') or align_id not in interatom.rdc_bc: 1938 data = None 1939 else: 1940 data = interatom.rdc_bc[align_id] 1941 1942 # Append the data. 1943 mc_data.append(data) 1944 1945 # The PCS data. 1946 elif data_id[0] == 'pcs': 1947 # Unpack the set. 1948 data_type, spin_id, align_id = data_id 1949 1950 # Get the spin container. 1951 spin = return_spin(spin_id) 1952 1953 # Does back-calculated data exist? 1954 if not hasattr(spin, 'pcs_bc'): 1955 self.calculate() 1956 1957 # The data. 1958 if not hasattr(spin, 'pcs_bc') or align_id not in spin.pcs_bc: 1959 data = None 1960 else: 1961 data = spin.pcs_bc[align_id] 1962 1963 # Append the data. 1964 mc_data.append(data) 1965 1966 # Return the data. 1967 return mc_data

1968 1969

1970 - def deselect(self, model_info, sim_index=None):

1971 """Deselect models or simulations. 1972 1973 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1974 @type model_info: int 1975 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will. 1976 @type sim_index: None or int 1977 """ 1978 1979 # Set the deselection flag. 1980 cdp.select = False

1981 1982

1983 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):

1984 """Duplicate the data specific to a single frame order data pipe. 1985 1986 @keyword pipe_from: The data pipe to copy the data from. 1987 @type pipe_from: str 1988 @keyword pipe_to: The data pipe to copy the data to. 1989 @type pipe_to: str 1990 @param model_info: The model index from model_loop(). 1991 @type model_info: int 1992 @keyword global_stats: The global statistics flag. 1993 @type global_stats: bool 1994 @keyword verbose: Unused. 1995 @type verbose: bool 1996 """ 1997 1998 # Check that the data pipe does not exist. 1999 if pipes.has_pipe(pipe_to): 2000 raise RelaxError("The data pipe '%s' already exists." % pipe_to) 2001 2002 # Create the pipe_to data pipe by copying. 2003 pipes.copy(pipe_from=pipe_from, pipe_to=pipe_to)

2004 2005

2006 - def eliminate(self, name, value, model_info, args, sim=None):

2007 """Model elimination method. 2008 2009 @param name: The parameter name. 2010 @type name: str 2011 @param value: The parameter value. 2012 @type value: float 2013 @param model_info: The model index from model_info(). 2014 @type model_info: int 2015 @param args: The elimination constant overrides. 2016 @type args: None or tuple of float 2017 @keyword sim: The Monte Carlo simulation index. 2018 @type sim: int 2019 @return: True if the model is to be eliminated, False otherwise. 2020 @rtype: bool 2021 """ 2022 2023 # Text to print out if a model failure occurs. 2024 text = "The %s parameter of %.5g is %s than %.5g, eliminating " 2025 if sim == None: 2026 text += "the model." 2027 else: 2028 text += "simulation %i." % sim 2029 2030 # Isotropic order parameter out of range. 2031 if name == 'cone_s1' and hasattr(cdp, 'cone_s1'): 2032 if cdp.cone_s1 > 1.0: 2033 print(text % ("cone S1 order", cdp.cone_s1, "greater", 1.0)) 2034 return True 2035 if cdp.cone_s1 < -0.125: 2036 print(text % ("cone S1 order", cdp.cone_s1, "less", -0.125)) 2037 return True 2038 2039 # Isotropic cone angle out of range. 2040 if name == 'cone_theta' and hasattr(cdp, 'cone_theta'): 2041 if cdp.cone_theta >= pi: 2042 print(text % ("cone opening angle theta", cdp.cone_theta, "greater", pi)) 2043 return True 2044 if cdp.cone_theta < 0.0: 2045 print(text % ("cone opening angle theta", cdp.cone_theta, "less", 0)) 2046 return True 2047 2048 # Pseudo-ellipse cone angles out of range (0.001 instead of 0.0 because of truncation in the numerical integration). 2049 if name == 'cone_theta_x' and hasattr(cdp, 'cone_theta_x'): 2050 if cdp.cone_theta_x >= pi: 2051 print(text % ("cone opening angle theta x", cdp.cone_theta_x, "greater", pi)) 2052 return True 2053 if cdp.cone_theta_x < 0.001: 2054 print(text % ("cone opening angle theta x", cdp.cone_theta_x, "less", 0.001)) 2055 return True 2056 if name == 'cone_theta_y' and hasattr(cdp, 'cone_theta_y'): 2057 if cdp.cone_theta_y >= pi: 2058 print(text % ("cone opening angle theta y", cdp.cone_theta_y, "greater", pi)) 2059 return True 2060 if cdp.cone_theta_y < 0.001: 2061 print(text % ("cone opening angle theta y", cdp.cone_theta_y, "less", 0.001)) 2062 return True 2063 2064 # Torsion angle out of range. 2065 if name == 'cone_sigma_max' and hasattr(cdp, 'cone_sigma_max'): 2066 if cdp.cone_sigma_max >= pi: 2067 print(text % ("torsion angle sigma_max", cdp.cone_sigma_max, "greater", pi)) 2068 return True 2069 if cdp.cone_sigma_max < 0.0: 2070 print(text % ("torsion angle sigma_max", cdp.cone_sigma_max, "less", 0.0)) 2071 return True 2072 2073 # No failure. 2074 return False

2075 2076

2077 - def get_param_names(self, model_info=None):

2078 """Return a vector of parameter names. 2079 2080 @keyword model_info: The model index from model_info(). 2081 @type model_info: int 2082 @return: The vector of parameter names. 2083 @rtype: list of str 2084 """ 2085 2086 # First update the model, if needed. 2087 self._update_model() 2088 2089 # Return the parameter list object. 2090 return cdp.params

2091 2092

2093 - def get_param_values(self, model_info=None, sim_index=None):

2094 """Return a vector of parameter values. 2095 2096 @keyword model_info: The model index from model_info(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 2097 @type model_info: int 2098 @keyword sim_index: The Monte Carlo simulation index. 2099 @type sim_index: int 2100 @return: The vector of parameter values. 2101 @rtype: list of str 2102 """ 2103 2104 # Assemble the values and return it. 2105 return self._assemble_param_vector(sim_index=sim_index)

2106 2107

2108 - def grid_search(self, lower=None, upper=None, inc=None, constraints=False, verbosity=0, sim_index=None):

2109 """Perform a grid search. 2110 2111 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 2112 @type lower: list of float 2113 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 2114 @type upper: list of float 2115 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. 2116 @type inc: int or list of int 2117 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 2118 @type constraints: bool 2119 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 2120 @type verbosity: int 2121 @keyword sim_index: The Monte Carlo simulation index. 2122 @type sim_index: None or int 2123 """ 2124 2125 # Test if the Frame Order model has been set up. 2126 if not hasattr(cdp, 'model'): 2127 raise RelaxNoModelError('Frame Order') 2128 2129 # Parameter scaling. 2130 scaling_matrix = self._assemble_scaling_matrix(data_types=self._base_data_types(), scaling=True) 2131 2132 # The number of parameters. 2133 n = self._param_num() 2134 2135 # If inc is an int, convert it into an array of that value. 2136 if isinstance(inc, int): 2137 incs = [inc]*n 2138 else: 2139 incs = inc 2140 2141 # Sanity check. 2142 if len(incs) != n: 2143 raise RelaxError("The size of the increment list %s does not match the number of parameters in %s." % (incs, cdp.params)) 2144 2145 # Initialise the grid increments structures. 2146 lower_list = lower 2147 upper_list = upper 2148 grid = [] 2149 """This structure is a list of lists. The first dimension corresponds to the model 2150 parameter. The second dimension are the grid node positions.""" 2151 2152 # Generate the grid. 2153 for i in range(n): 2154 # Fixed parameter. 2155 if incs[i] == None: 2156 grid.append(None) 2157 continue 2158 2159 # Reset. 2160 dist_type = None 2161 end_point = True 2162 2163 # The pivot point. 2164 if cdp.params[i] == 'pivot_x': 2165 lower = cdp.pivot[0] - 10.0 2166 upper = cdp.pivot[0] + 10.0 2167 elif cdp.params[i] == 'pivot_y': 2168 lower = cdp.pivot[1] - 10.0 2169 upper = cdp.pivot[1] + 10.0 2170 elif cdp.params[i] == 'pivot_z': 2171 lower = cdp.pivot[2] - 10.0 2172 upper = cdp.pivot[2] + 10.0 2173 2174 # Average domain position translation (in a +/- 5 Angstrom box). 2175 if cdp.params[i] in ['ave_pos_x', 'ave_pos_y', 'ave_pos_z']: 2176 lower = -5 2177 upper = 5 2178 2179 # Linear angle grid from 0 to one inc before 2pi. 2180 if cdp.params[i] in ['ave_pos_alpha', 'ave_pos_gamma', 'eigen_alpha', 'eigen_gamma', 'axis_phi']: 2181 lower = 0.0 2182 upper = 2*pi * (1.0 - 1.0/incs[i]) 2183 2184 # Arccos grid from 0 to pi. 2185 if cdp.params[i] in ['ave_pos_beta', 'eigen_beta', 'axis_theta']: 2186 # Change the default increment numbers. 2187 if not isinstance(inc, list): 2188 incs[i] = int(incs[i] / 2) + 1 2189 2190 # The distribution type and end point. 2191 dist_type = 'acos' 2192 end_point = False 2193 2194 # Set the default bounds. 2195 lower = 0.0 2196 upper = pi 2197 2198 # The isotropic cone order parameter. 2199 if cdp.params[i] == 'cone_s1': 2200 lower = -0.125 2201 upper = 1.0 2202 2203 # Linear angle grid from 0 to pi excluding the outer points. 2204 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']: 2205 lower = pi * (1.0/incs[i]) 2206 upper = pi * (1.0 - 1.0/incs[i]) 2207 2208 # Over-ride the bounds. 2209 if lower_list: 2210 lower = lower_list[i] 2211 if upper_list: 2212 upper = upper_list[i] 2213 2214 # Append the grid row. 2215 row = self._grid_row(incs[i], lower, upper, dist_type=dist_type, end_point=end_point) 2216 grid.append(row) 2217 2218 # Remove an inc if the end point has been removed. 2219 if not end_point: 2220 incs[i] -= 1 2221 2222 # Total number of points. 2223 total_pts = 1 2224 for i in range(n): 2225 # Fixed parameter. 2226 if grid[i] == None: 2227 continue 2228 2229 total_pts = total_pts * len(grid[i]) 2230 2231 # Check the number. 2232 max_pts = 50e6 2233 if total_pts > max_pts: 2234 raise RelaxError("The total number of grid points '%s' exceeds the maximum of '%s'." % (total_pts, int(max_pts))) 2235 2236 # Build the points array. 2237 pts = zeros((total_pts, n), float64) 2238 indices = zeros(n, int) 2239 for i in range(total_pts): 2240 # Loop over the dimensions. 2241 for j in range(n): 2242 # Fixed parameter. 2243 if grid[j] == None: 2244 # Get the current parameter value (pivot, assuming the pivot point is always at the start of the parameter array). 2245 if cdp.params[j] in ['pivot_x', 'pivot_y', 'pivot_z']: 2246 pts[i, j] = cdp.pivot[j] / scaling_matrix[j, j] 2247 2248 # Get the current parameter value (normal parameter). 2249 else: 2250 pts[i, j] = getattr(cdp, cdp.params[j]) / scaling_matrix[j, j] 2251 2252 # Add the point coordinate. 2253 else: 2254 pts[i, j] = grid[j][indices[j]] / scaling_matrix[j, j] 2255 2256 # Increment the step positions. 2257 for j in range(n): 2258 if incs[j] != None and indices[j] < incs[j]-1: 2259 indices[j] += 1 2260 break # Exit so that the other step numbers are not incremented. 2261 else: 2262 indices[j] = 0 2263 2264 # Minimisation. 2265 self.minimise(min_algor='grid', min_options=pts, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

2266 2267

2268 - def is_spin_param(self, name):

2269 """State that the parameter is not spin specific. 2270 2271 @param name: The name of the parameter. 2272 @type name: str 2273 @return: False. 2274 @rtype: bool 2275 """ 2276 2277 # Not spin specific! 2278 return False

2279 2280

2281 - def map_bounds(self, param, spin_id=None):

2282 """Create bounds for the OpenDX mapping function. 2283 2284 @param param: The name of the parameter to return the lower and upper bounds of. 2285 @type param: str 2286 @param spin_id: The spin identification string (unused). 2287 @type spin_id: None 2288 @return: The upper and lower bounds of the parameter. 2289 @rtype: list of float 2290 """ 2291 2292 # Average domain position. 2293 if param in ['ave_pos_x', 'ave_pos_y', 'ave_pos_z']: 2294 return [-100.0, 100] 2295 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']: 2296 return [0.0, 2*pi] 2297 2298 # Axis spherical coordinate theta. 2299 if param == 'axis_theta': 2300 return [0.0, pi] 2301 2302 # Axis spherical coordinate phi. 2303 if param == 'axis_phi': 2304 return [0.0, 2*pi] 2305 2306 # Cone angle. 2307 if param == 'cone_theta': 2308 return [0.0, pi]

2309 2310

2311 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

2312 """Minimisation function. 2313 2314 @param min_algor: The minimisation algorithm to use. 2315 @type min_algor: str 2316 @param min_options: An array of options to be used by the minimisation algorithm. 2317 @type min_options: array of str 2318 @param func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 2319 @type func_tol: None or float 2320 @param grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 2321 @type grad_tol: None or float 2322 @param max_iterations: The maximum number of iterations for the algorithm. 2323 @type max_iterations: int 2324 @param constraints: If True, constraints are used during optimisation. 2325 @type constraints: bool 2326 @param scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 2327 @type scaling: bool 2328 @param verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 2329 @type verbosity: int 2330 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 2331 @type sim_index: None or int 2332 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 2333 @type lower: array of numbers 2334 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 2335 @type upper: array of numbers 2336 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search. 2337 @type inc: array of int 2338 """ 2339 2340 # Set up the target function for direct calculation. 2341 model, param_vector, data_types, scaling_matrix = self._target_fn_setup(sim_index=sim_index, scaling=scaling) 2342 2343 # Constraints not implemented yet. 2344 if constraints: 2345 # Turn the constraints off. 2346 constraints = False 2347 2348 # Remove the method of multipliers arg. 2349 if not search('^[Gg]rid', min_algor): 2350 min_algor = min_options[0] 2351 min_options = min_options[1:] 2352 2353 # Throw a warning. 2354 warn(RelaxWarning("Constraints are as of yet not implemented - turning this option off.")) 2355 2356 # Grid search. 2357 if search('^[Gg]rid', min_algor): 2358 results = grid_point_array(func=model.func, args=(), points=min_options, verbosity=verbosity) 2359 2360 # Minimisation. 2361 else: 2362 results = generic_minimise(func=model.func, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, full_output=True, print_flag=verbosity) 2363 2364 # Unpack the results. 2365 self._unpack_opt_results(results, scaling, scaling_matrix, sim_index) 2366 2367 # Store the back-calculated data. 2368 self._store_bc_data(model)

2369 2370

2371 - def model_desc(self, model_info):

2372 """Return a description of the model. 2373 2374 @param model_info: The model index from model_loop(). 2375 @type model_info: int 2376 @return: The model description. 2377 @rtype: str 2378 """ 2379 2380 return ""

2381 2382

2383 - def model_loop(self):

2384 """Dummy generator method. 2385 2386 In this case only a single model per spin system is assumed. Hence the yielded data is the 2387 spin container object. 2388 2389 2390 @return: Information about the model which for this analysis is the spin container. 2391 @rtype: SpinContainer instance 2392 """ 2393 2394 # Don't return anything, just loop once. 2395 yield None

2396 2397

2398 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):

2399 """Return the k, n, and chi2 model statistics. 2400 2401 k - number of parameters. 2402 n - number of data points. 2403 chi2 - the chi-squared value. 2404 2405 2406 @keyword model_info: Unused. 2407 @type model_info: None 2408 @keyword spin_id: The spin identification string (unused). 2409 @type spin_id: None 2410 @keyword global_stats: Unused. 2411 @type global_stats: None 2412 @return: The optimisation statistics, in tuple format, of the number of 2413 parameters (k), the number of data points (n), and the chi-squared 2414 value (chi2). 2415 @rtype: tuple of (int, int, float) 2416 """ 2417 2418 # Count the number of parameters. 2419 param_vector = self._assemble_param_vector() 2420 k = len(param_vector) 2421 2422 # The number of data points. 2423 n = len(cdp.align_tensors.reduction) 2424 2425 # The chi2 value. 2426 if not hasattr(cdp, 'chi2'): 2427 raise RelaxError("Statistics are not available, most likely because the model has not been optimised.") 2428 chi2 = cdp.chi2 2429 2430 # Return the data. 2431 return k, n, chi2

2432 2433

2434 - def model_type(self):

2435 """Return the type of the model, either being 'local' or 'global'. 2436 2437 @return: The model type, one of 'local' or 'global'. 2438 @rtype: str 2439 """ 2440 2441 return 'global'

2442 2443

2444 - def return_error(self, data_id):

2445 """Return the alignment tensor error structure. 2446 2447 @param data_id: The data set as yielded by the base_data_loop() generator method. 2448 @type data_id: list of str 2449 @return: The array of tensor error values. 2450 @rtype: list of float 2451 """ 2452 2453 # Initialise the MC data structure. 2454 mc_errors = [] 2455 2456 # The RDC data. 2457 if data_id[0] == 'rdc': 2458 # Unpack the set. 2459 data_type, spin_id1, spin_id2, align_id = data_id 2460 2461 # Get the interatomic data container. 2462 interatom = return_interatom(spin_id1, spin_id2) 2463 2464 # Do errors exist? 2465 if not hasattr(interatom, 'rdc_err'): 2466 raise RelaxError("The RDC errors are missing for interatomic data container between spins '%s' and '%s'." % (spin_id1, spin_id2)) 2467 2468 # Handle missing data. 2469 if align_id not in interatom.rdc_err: 2470 mc_errors.append(None) 2471 2472 # Append the data. 2473 else: 2474 mc_errors.append(interatom.rdc_err[align_id]) 2475 2476 # The PCS data. 2477 elif data_id[0] == 'pcs': 2478 # Unpack the set. 2479 data_type, spin_id, align_id = data_id 2480 2481 # Get the spin container. 2482 spin = return_spin(spin_id) 2483 2484 # Do errors exist? 2485 if not hasattr(spin, 'pcs_err'): 2486 raise RelaxError("The PCS errors are missing for spin '%s'." % spin_id) 2487 2488 # Handle missing data. 2489 if align_id not in spin.pcs_err: 2490 mc_errors.append(None) 2491 2492 # Append the data. 2493 else: 2494 mc_errors.append(spin.pcs_err[align_id]) 2495 2496 # Return the errors. 2497 return mc_errors

2498 2499

2500 - def return_units(self, param):

2501 """Return a string representing the parameters units. 2502 2503 @param param: The name of the parameter to return the units string for. 2504 @type param: str 2505 @return: The parameter units string. 2506 @rtype: str 2507 """ 2508 2509 # Average domain position. 2510 if param in ['ave_pos_x', 'ave_pos_y', 'ave_pos_z']: 2511 return 'Angstrom' 2512 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']: 2513 return 'rad' 2514 2515 # Eigenframe angles. 2516 if param in ['eigen_alpha', 'eigen_beta', 'eigen_gamma', 'axis_theta', 'axis_phi']: 2517 return 'rad' 2518 2519 # Cone angles. 2520 if param in ['cone_theta_x', 'cone_theta_y', 'cone_sigma_max', 'cone_theta']: 2521 return 'rad'

2522 2523

2524 - def set_error(self, model_info, index, error):

2525 """Set the parameter errors. 2526 2527 @param model_info: The model information originating from model_loop() (unused). 2528 @type model_info: None 2529 @param index: The index of the parameter to set the errors for. 2530 @type index: int 2531 @param error: The error value. 2532 @type error: float 2533 """ 2534 2535 # Parameter increment counter. 2536 inc = 0 2537 2538 # Loop over the residue specific parameters. 2539 for param in self.data_names(set='params'): 2540 # Not a parameter of the model. 2541 if param not in cdp.params: 2542 continue 2543 2544 # Return the parameter array. 2545 if index == inc: 2546 setattr(cdp, param + "_err", error) 2547 2548 # Increment. 2549 inc = inc + 1 2550 2551 # Add some additional parameters. 2552 if cdp.model == 'iso cone, free rotor' and inc == index: 2553 setattr(cdp, 'cone_theta_err', error)

2554 2555 2556

2557 - def set_selected_sim(self, model_info, select_sim):

2558 """Set the simulation selection flag for the spin. 2559 2560 @param model_info: The model information originating from model_loop() (unused). 2561 @type model_info: None 2562 @param select_sim: The selection flag for the simulations. 2563 @type select_sim: bool 2564 """ 2565 2566 # Set the array. 2567 cdp.select_sim = deepcopy(select_sim)

2568 2569

2570 - def sim_init_values(self):

2571 """Initialise the Monte Carlo parameter values.""" 2572 2573 # Get the parameter object names. 2574 param_names = self.data_names(set='params') 2575 2576 # The model parameters. 2577 model_params = deepcopy(cdp.params) 2578 2579 # Add some additional parameters. 2580 if cdp.model == 'iso cone, free rotor': 2581 param_names.append('cone_theta') 2582 model_params.append('cone_theta') 2583 2584 # Get the minimisation statistic object names. 2585 min_names = self.data_names(set='min') 2586 2587 2588 # Test if Monte Carlo parameter values have already been set. 2589 ############################################################# 2590 2591 # Loop over all the parameter names. 2592 for object_name in param_names: 2593 # Not a parameter of the model. 2594 if object_name not in model_params: 2595 continue 2596 2597 # Name for the simulation object. 2598 sim_object_name = object_name + '_sim' 2599 2600 # Test if the simulation object already exists. 2601 if hasattr(cdp, sim_object_name): 2602 raise RelaxError("Monte Carlo parameter values have already been set.") 2603 2604 2605 # Set the Monte Carlo parameter values. 2606 ####################################### 2607 2608 # Loop over all the data names. 2609 for object_name in param_names: 2610 # Not a parameter of the model. 2611 if object_name not in model_params: 2612 continue 2613 2614 # Name for the simulation object. 2615 sim_object_name = object_name + '_sim' 2616 2617 # Create the simulation object. 2618 setattr(cdp, sim_object_name, []) 2619 2620 # Get the simulation object. 2621 sim_object = getattr(cdp, sim_object_name) 2622 2623 # Loop over the simulations. 2624 for j in range(cdp.sim_number): 2625 # Copy and append the data. 2626 sim_object.append(deepcopy(getattr(cdp, object_name))) 2627 2628 # Loop over all the minimisation object names. 2629 for object_name in min_names: 2630 # Name for the simulation object. 2631 sim_object_name = object_name + '_sim' 2632 2633 # Create the simulation object. 2634 setattr(cdp, sim_object_name, []) 2635 2636 # Get the simulation object. 2637 sim_object = getattr(cdp, sim_object_name) 2638 2639 # Loop over the simulations. 2640 for j in range(cdp.sim_number): 2641 # Copy and append the data. 2642 sim_object.append(deepcopy(getattr(cdp, object_name)))

2643 2644

2645 - def sim_pack_data(self, data_id, sim_data):

2646 """Pack the Monte Carlo simulation data. 2647 2648 @param data_id: The data set as yielded by the base_data_loop() generator method. 2649 @type data_id: list of str 2650 @param sim_data: The Monte Carlo simulation data. 2651 @type sim_data: list of float 2652 """ 2653 2654 # The RDC data. 2655 if data_id[0] == 'rdc': 2656 # Unpack the set. 2657 data_type, spin_id1, spin_id2, align_id = data_id 2658 2659 # Get the interatomic data container. 2660 interatom = return_interatom(spin_id1, spin_id2) 2661 2662 # Initialise. 2663 if not hasattr(interatom, 'rdc_sim'): 2664 interatom.rdc_sim = {} 2665 2666 # Store the data structure. 2667 interatom.rdc_sim[align_id] = [] 2668 for i in range(cdp.sim_number): 2669 interatom.rdc_sim[align_id].append(sim_data[i][0]) 2670 2671 # The PCS data. 2672 elif data_id[0] == 'pcs': 2673 # Unpack the set. 2674 data_type, spin_id, align_id = data_id 2675 2676 # Get the spin container. 2677 spin = return_spin(spin_id) 2678 2679 # Initialise. 2680 if not hasattr(spin, 'pcs_sim'): 2681 spin.pcs_sim = {} 2682 2683 # Store the data structure. 2684 spin.pcs_sim[data_id[2]] = [] 2685 for i in range(cdp.sim_number): 2686 spin.pcs_sim[data_id[2]].append(sim_data[i][0])

2687 2688

2689 - def sim_return_param(self, model_info, index):

2690 """Return the array of simulation parameter values. 2691 2692 @param model_info: The model information originating from model_loop(). 2693 @type model_info: unknown 2694 @param index: The index of the parameter to return the array of values for. 2695 @type index: int 2696 """ 2697 2698 # Parameter increment counter. 2699 inc = 0 2700 2701 # Get the parameter object names. 2702 param_names = self.data_names(set='params') 2703 2704 # Loop over the parameters. 2705 for param in param_names: 2706 # Not a parameter of the model. 2707 if param not in cdp.params: 2708 continue 2709 2710 # Return the parameter array. 2711 if index == inc: 2712 return getattr(cdp, param + "_sim") 2713 2714 # Increment. 2715 inc = inc + 1 2716 2717 # Add some additional parameters. 2718 if cdp.model == 'iso cone, free rotor' and inc == index: 2719 return getattr(cdp, 'cone_theta_sim')

2720 2721

2722 - def sim_return_selected(self, model_info):

2723 """Return the array of selected simulation flags for the spin. 2724 2725 @param model_info: The model information originating from model_loop() (unused). 2726 @type model_info: None 2727 @return: The array of selected simulation flags. 2728 @rtype: list of int 2729 """ 2730 2731 # Return the array. 2732 return cdp.select_sim

Source Code for Package specific_analyses.frame_order