Package specific_analyses :: Package frame_order :: Module uf
[hide private]
[frames] | no frames]

Source Code for Module specific_analyses.frame_order.uf

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2009-2014 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """Module for all of the frame order specific user functions.""" 
 24   
 25  # Python module imports. 
 26  from copy import deepcopy 
 27  from math import pi 
 28  from numpy import array, dot, eye, float64, transpose, zeros 
 29  import sys 
 30   
 31  # relax module imports. 
 32  from lib.arg_check import is_float_array 
 33  from lib.errors import RelaxError 
 34  from lib.frame_order.rotor_axis import create_rotor_axis_alpha, create_rotor_axis_euler, create_rotor_axis_spherical 
 35  from lib.geometry.rotations import euler_to_R_zyz, two_vect_to_R 
 36  from lib.io import open_write_file 
 37  from lib.order import order_parameters 
 38  from lib.structure.cones import Iso_cone, Pseudo_elliptic 
 39  from lib.structure.geometric import generate_vector_residues 
 40  from lib.structure.internal.object import Internal 
 41  from lib.structure.represent.cone import cone 
 42  from lib.structure.represent.rotor import rotor_pdb 
 43  from lib.text.sectioning import subsection 
 44  from pipe_control import pipes 
 45  from pipe_control.structure.mass import pipe_centre_of_mass 
 46  from specific_analyses.frame_order.data import domain_moving, translation_fixed 
 47  from specific_analyses.frame_order.parameters import update_model 
 48   
 49   
50 -def average_position(pivot='com', translation=True):
51 """Set up the mechanics of the average domain position. 52 53 @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. 54 @type pivot: str 55 @keyword translation: If True, translation to the average domain position will be allowed. If False, then translation will not occur. 56 @type translation: bool 57 """ 58 59 # Test if the current data pipe exists. 60 pipes.test() 61 62 # Check the pivot value. 63 if pivot not in ['com', 'motional']: 64 raise RelaxError("The pivot for the rotation to the average domain position must be either 'com' or 'motional'.") 65 66 # Store the data. 67 cdp.ave_pos_pivot = pivot 68 cdp.ave_pos_translation = translation
69 70
71 -def num_int_pts(num=200000):
72 """Set the number of integration points to use in the quasi-random Sobol' sequence. 73 74 @keyword num: The number of integration points. 75 @type num: int 76 """ 77 78 # Test if the current data pipe exists. 79 pipes.test() 80 81 # Store the value. 82 cdp.num_int_pts = num
83 84
85 -def pdb_ave_pos(file=None, dir=None, force=False):
86 """Create a PDB file of the molecule with the moving domains shifted to the average position. 87 88 @keyword file: The name of the file for the average molecule structure. 89 @type file: str 90 @keyword dir: The name of the directory to place the PDB file into. 91 @type dir: str 92 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 93 @type force: bool 94 """ 95 96 # Printout. 97 subsection(file=sys.stdout, text="Creating a PDB file with the moving domains shifted to the average position.") 98 99 # Make a copy of the structural object (so as to preserve the original structure). 100 structure = deepcopy(cdp.structure) 101 102 # First rotate the moving domain to the average position. 103 R = zeros((3, 3), float64) 104 if hasattr(cdp, 'ave_pos_alpha'): 105 euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 106 else: 107 euler_to_R_zyz(0.0, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 108 if cdp.ave_pos_pivot == 'com': 109 origin = pipe_centre_of_mass(atom_id=domain_moving(), verbosity=0) 110 else: 111 origin = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) 112 structure.rotate(R=R, origin=origin, atom_id=domain_moving()) 113 114 # Then translate the moving domain. 115 if not translation_fixed(): 116 structure.translate(T=[cdp.ave_pos_x, cdp.ave_pos_y, cdp.ave_pos_z], atom_id=domain_moving()) 117 118 # Write out the PDB file. 119 file = open_write_file(file_name=file, dir=dir, force=force) 120 structure.write_pdb(file=file) 121 file.close()
122 123
124 -def pdb_distribution(file=None, dir=None, force=False):
125 """Create a PDB file of a distribution of positions coving the full dynamics of the moving domain. 126 127 @keyword file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model. 128 @type file: str 129 @keyword dir: The name of the directory to place the PDB file into. 130 @type dir: str 131 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 132 @type force: bool 133 """ 134 135 # Printout. 136 subsection(file=sys.stdout, text="Creating a PDB file of a distribution of positions coving the full dynamics of the moving domain.")
137 138
139 -def pdb_geometric_rep(file=None, dir=None, size=30.0, inc=36, force=False, neg_cone=True):
140 """Create a PDB file containing a geometric object representing the frame order dynamics. 141 142 @keyword file: The name of the file of the PDB representation of the frame order dynamics to create. 143 @type file: str 144 @keyword dir: The name of the directory to place the PDB file into. 145 @type dir: str 146 @keyword size: The size of the geometric object in Angstroms. 147 @type size: float 148 @keyword inc: The number of increments for the filling of the cone objects. 149 @type inc: int 150 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 151 @type force: bool 152 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation. This is ignored for the rotor models. 153 @type neg_cone: bool 154 """ 155 156 # Printout. 157 subsection(file=sys.stdout, text="Creating a PDB file containing a geometric object representing the frame order dynamics.") 158 159 # Monte Carlo simulation flag. 160 sim = False 161 num_sim = 0 162 if hasattr(cdp, 'sim_number'): 163 sim = True 164 num_sim = cdp.sim_number 165 166 # The inversion matrix. 167 inv_mat = -eye(3) 168 169 # Create the structural object. 170 structure = Internal() 171 172 # Create model for the positive and negative images. 173 if cdp.model in ['rotor', 'free rotor']: 174 neg_cone = False 175 model = structure.add_model(model=1) 176 model_num = 1 177 if neg_cone: 178 model_neg = structure.add_model(model=2) 179 model_num = 2 180 181 # The pivot point. 182 pivot = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) 183 184 # The rotor object. 185 if cdp.model in ['rotor', 'free rotor', 'iso cone', 'iso cone, free rotor', 'pseudo-ellipse']: 186 # The rotor angle. 187 if cdp.model in ['free rotor', 'iso cone, free rotor']: 188 rotor_angle = pi 189 else: 190 rotor_angle = cdp.cone_sigma_max 191 192 # Get the CoM of the entire molecule to use as the centre of the rotor. 193 com = pivot 194 if cdp.model in ['rotor', 'free rotor']: 195 com = pipe_centre_of_mass(verbosity=0) 196 197 # Generate the rotor axis. 198 if cdp.model in ['rotor']: 199 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com) 200 elif cdp.model in ['free rotor', 'iso cone', 'iso cone, free rotor']: 201 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi) 202 else: 203 axis = create_rotor_axis_euler(alpha=cdp.eigen_alpha, beta=cdp.eigen_beta, gamma=cdp.eigen_gamma) 204 205 # The size of the rotor, taking the 30 Angstrom cone representation into account. 206 if cdp.model in ['rotor', 'free rotor']: 207 span = 20e-10 208 else: 209 span = 35e-10 210 211 # Add the rotor object to the structure as a new molecule. 212 rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis, axis_pt=pivot, centre=com, span=span, blade_length=5e-10, staggered=True) 213 214 # Create the molecule. 215 structure.add_molecule(name='rest') 216 217 # Alias the molecules. 218 mol = model.mol[0] 219 if neg_cone: 220 mol_neg = model_neg.mol[0] 221 222 223 # The pivot point. 224 ################## 225 226 # Add the pivot point. 227 structure.add_atom(mol_name=cdp.model, pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=pivot, element='C') 228 229 230 # The axes. 231 ########### 232 233 # The spherical angles. 234 if cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 235 # Print out. 236 print("\nGenerating the z-axis system.") 237 238 # The axis. 239 if cdp.model in ['rotor']: 240 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com) 241 else: 242 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi) 243 print(("Central axis: %s." % axis)) 244 245 # Rotations and inversions. 246 axis_pos = axis 247 axis_neg = dot(inv_mat, axis) 248 249 # Simulation central axis. 250 axis_sim_pos = None 251 axis_sim_neg = None 252 if sim: 253 # Init. 254 axis_sim = zeros((cdp.sim_number, 3), float64) 255 256 # Fill the structure. 257 for i in range(cdp.sim_number): 258 if cdp.model in ['rotor']: 259 axis_sim[i] = create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[i], pivot=pivot, point=com) 260 else: 261 axis_sim[i] = create_rotor_axis_spherical(theta=cdp.axis_theta_sim[i], phi=cdp.axis_phi_sim[i]) 262 263 # Inversion. 264 axis_sim_pos = axis_sim 265 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos))) 266 267 # Generate the axis vectors. 268 print("\nGenerating the axis vectors.") 269 res_num = 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=pivot, scale=size) 270 271 # The negative. 272 if neg_cone: 273 res_num = 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=pivot, scale=size) 274 275 # The full axis system. 276 else: 277 # Print out. 278 print("\nGenerating the full axis system.") 279 280 # The axis system. 281 axes = zeros((3, 3), float64) 282 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes) 283 print(("Axis system:\n%s" % axes)) 284 285 # Rotations and inversions. 286 axes_pos = axes 287 axes_neg = dot(inv_mat, axes) 288 289 # Simulations 290 axes_sim_pos = None 291 axes_sim_neg = None 292 if sim: 293 # Init. 294 axes_sim_pos = zeros((cdp.sim_number, 3, 3), float64) 295 axes_sim_neg = zeros((cdp.sim_number, 3, 3), float64) 296 297 # Fill the structure. 298 for i in range(cdp.sim_number): 299 # The positive system. 300 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_pos[i]) 301 302 # The negative system. 303 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_neg[i]) 304 axes_sim_neg[i] = dot(inv_mat, axes_sim_neg[i]) 305 306 # Generate the axis vectors. 307 print("\nGenerating the axis vectors.") 308 label = ['x', 'y', 'z'] 309 for j in range(len(label)): 310 # The simulation data. 311 axis_sim_pos = None 312 axis_sim_neg = None 313 if sim: 314 axis_sim_pos = axes_sim_pos[:,:, j] 315 axis_sim_neg = axes_sim_neg[:,:, j] 316 317 # The vectors. 318 res_num = 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=pivot, scale=size) 319 if neg_cone: 320 res_num = 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=pivot, scale=size) 321 322 323 # The cone object. 324 ################## 325 326 # Skip models missing a cone. 327 if cdp.model not in ['rotor', 'free rotor']: 328 # The rotation matrix (rotation from the z-axis to the cone axis). 329 if cdp.model not in ['iso cone', 'iso cone, torsionless', 'iso cone, free rotor']: 330 R = axes 331 else: 332 R = zeros((3, 3), float64) 333 two_vect_to_R(array([0, 0, 1], float64), axis, R) 334 335 # Average position rotation. 336 R_pos = R 337 R_neg = dot(inv_mat, R) 338 339 # The pseudo-ellipse cone object. 340 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 341 cone_obj = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y) 342 343 # The isotropic cone object. 344 else: 345 # The angle. 346 if hasattr(cdp, 'cone_theta'): 347 cone_theta = cdp.cone_theta 348 elif hasattr(cdp, 'cone_s1'): 349 cone_theta = order_parameters.iso_cone_S_to_theta(cdp.cone_s1) 350 351 # The object. 352 cone_obj = Iso_cone(cone_theta) 353 354 # Create the positive and negative cones. 355 cone(mol=mol, cone_obj=cone_obj, start_res=mol.res_num[-1]+1, apex=pivot, R=R_pos, inc=inc, distribution='regular', axis_flag=False) 356 357 # The negative. 358 if neg_cone: 359 cone(mol=mol_neg, cone_obj=cone_obj, start_res=mol_neg.res_num[-1]+1, apex=pivot, R=R_neg, inc=inc, distribution='regular', axis_flag=False) 360 361 362 # Create the PDB file. 363 ###################### 364 365 # Print out. 366 print("\nGenerating the PDB file.") 367 368 # Write the file. 369 pdb_file = open_write_file(file, dir, force=force) 370 structure.write_pdb(pdb_file) 371 pdb_file.close()
372 373
374 -def pdb_model(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):
375 """Create 3 different PDB files for representing the frame order dynamics of the system. 376 377 @keyword ave_pos_file: The name of the file for the average molecule structure. 378 @type ave_pos_file: str or None 379 @keyword rep_file: The name of the file of the PDB representation of the frame order dynamics to create. 380 @type rep_file: str or None 381 @keyword dist_file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model. 382 @type dist_file: str or None 383 @keyword dir: The name of the directory to place the PDB file into. 384 @type dir: str 385 @keyword size: The size of the geometric object in Angstroms. 386 @type size: float 387 @keyword inc: The number of increments for the filling of the cone objects. 388 @type inc: int 389 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 390 @type force: bool 391 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation. 392 @type neg_cone: bool 393 """ 394 395 # Check that at least one PDB file name is given. 396 if not ave_pos_file and not rep_file and not dist_file: 397 raise RelaxError("Minimally one PDB file name must be supplied.") 398 399 # Test if the current data pipe exists. 400 pipes.test() 401 402 # Create the average position structure. 403 if ave_pos_file: 404 pdb_ave_pos(file=ave_pos_file, dir=dir, force=force) 405 406 # Nothing more to do for the rigid model. 407 if cdp.model == 'rigid': 408 return 409 410 # Create the geometric representation. 411 if rep_file: 412 pdb_geometric_rep(file=rep_file, dir=dir, size=size, inc=inc, force=force, neg_cone=neg_cone) 413 414 # Create the distribution. 415 if dist_file: 416 pdb_distribution(file=dist_file, dir=dir, force=force)
417 418
419 -def pivot(pivot=None, order=1, fix=False):
420 """Set the pivot point for the 2 body motion. 421 422 @keyword pivot: The pivot point of the two bodies (domains, etc.) in the structural coordinate system. 423 @type pivot: list of num 424 @keyword order: The ordinal number of the pivot point. The value of 1 is for the first pivot point, the value of 2 for the second pivot point, and so on. 425 @type order: int 426 @keyword fix: A flag specifying if the pivot point should be fixed during optimisation. 427 @type fix: bool 428 """ 429 430 # Test if the current data pipe exists. 431 pipes.test() 432 433 # Store the fixed flag. 434 cdp.pivot_fixed = fix 435 436 # No pivot given, so update the model if needed and quit. 437 if pivot == None: 438 if hasattr(cdp, 'model'): 439 update_model() 440 return 441 442 # Convert the pivot to a numpy array. 443 pivot = array(pivot, float64) 444 445 # Check the pivot validity. 446 is_float_array(pivot, name='pivot point', size=3) 447 448 # Store the pivot point and fixed flag. 449 if order == 1: 450 cdp.pivot_x = pivot[0] 451 cdp.pivot_y = pivot[1] 452 cdp.pivot_z = pivot[2] 453 else: 454 # The variable names. 455 name_x = 'pivot_x_%i' % order 456 name_y = 'pivot_y_%i' % order 457 name_z = 'pivot_z_%i' % order 458 459 # Store the variables. 460 setattr(cdp, name_x, pivot[0]) 461 setattr(cdp, name_y, pivot[1]) 462 setattr(cdp, name_z, pivot[2]) 463 464 # Update the model. 465 if hasattr(cdp, 'model'): 466 update_model()
467 468
469 -def quad_int(flag=False):
470 """Turn the high precision Scipy quadratic numerical integration on or off. 471 472 @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. 473 @type flag: bool 474 """ 475 476 # Test if the current data pipe exists. 477 pipes.test() 478 479 # Store the flag. 480 cdp.quad_int = flag
481 482
483 -def ref_domain(ref=None):
484 """Set the reference domain for the frame order, multi-domain models. 485 486 @param ref: The reference domain. 487 @type ref: str 488 """ 489 490 # Test if the current data pipe exists. 491 pipes.test() 492 493 # Check that the domain is defined. 494 if not hasattr(cdp, 'domain') or ref not in list(cdp.domain.keys()): 495 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % ref) 496 497 # Test if the reference domain exists. 498 exists = False 499 for tensor_cont in cdp.align_tensors: 500 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref: 501 exists = True 502 if not exists: 503 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 504 505 # Set the reference domain. 506 cdp.ref_domain = ref 507 508 # Update the model. 509 if hasattr(cdp, 'model'): 510 update_model()
511 512
513 -def select_model(model=None):
514 """Select the Frame Order model. 515 516 @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', 'double rotor'. 517 @type model: str 518 """ 519 520 # Test if the current data pipe exists. 521 pipes.test() 522 523 # Test if the model name exists. 524 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', 'double rotor']: 525 raise RelaxError("The model name " + repr(model) + " is invalid.") 526 527 # Set the model 528 cdp.model = model 529 530 # Initialise the list of model parameters. 531 cdp.params = [] 532 533 # Set the integration method if needed. 534 if not hasattr(cdp, 'quad_int'): 535 # Scipy quadratic numerical integration. 536 if cdp.model in []: 537 cdp.quad_int = True 538 539 # Quasi-random numerical integration. 540 else: 541 cdp.quad_int = False 542 543 # Update the model. 544 update_model()
545