specific_analyses.frame

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.pipes import check_pipe 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 check_pipe() 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 check_pipe() 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 # The internal structural selection object. 103 selection = structure.selection(atom_id=domain_moving()) 104 105 # First rotate the moving domain to the average position. 106 R = zeros((3, 3), float64) 107 if hasattr(cdp, 'ave_pos_alpha'): 108 euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 109 else: 110 euler_to_R_zyz(0.0, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 111 if cdp.ave_pos_pivot == 'com': 112 origin = pipe_centre_of_mass(atom_id=domain_moving(), verbosity=0) 113 else: 114 origin = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) 115 structure.rotate(R=R, origin=origin, selection=selection) 116 117 # Then translate the moving domain. 118 if not translation_fixed(): 119 structure.translate(T=[cdp.ave_pos_x, cdp.ave_pos_y, cdp.ave_pos_z], selection=selection) 120 121 # Write out the PDB file. 122 file = open_write_file(file_name=file, dir=dir, force=force) 123 structure.write_pdb(file=file) 124 file.close()

125 126

127 -def pdb_distribution(file=None, dir=None, force=False):

128 """Create a PDB file of a distribution of positions coving the full dynamics of the moving domain. 129 130 @keyword file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model. 131 @type file: str 132 @keyword dir: The name of the directory to place the PDB file into. 133 @type dir: str 134 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 135 @type force: bool 136 """ 137 138 # Printout. 139 subsection(file=sys.stdout, text="Creating a PDB file of a distribution of positions coving the full dynamics of the moving domain.")

140 141

142 -def pdb_geometric_rep(file=None, dir=None, size=30.0, inc=36, force=False, neg_cone=True):

143 """Create a PDB file containing a geometric object representing the frame order dynamics. 144 145 @keyword file: The name of the file of the PDB representation of the frame order dynamics to create. 146 @type file: str 147 @keyword dir: The name of the directory to place the PDB file into. 148 @type dir: str 149 @keyword size: The size of the geometric object in Angstroms. 150 @type size: float 151 @keyword inc: The number of increments for the filling of the cone objects. 152 @type inc: int 153 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 154 @type force: bool 155 @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. 156 @type neg_cone: bool 157 """ 158 159 # Printout. 160 subsection(file=sys.stdout, text="Creating a PDB file containing a geometric object representing the frame order dynamics.") 161 162 # Monte Carlo simulation flag. 163 sim = False 164 num_sim = 0 165 if hasattr(cdp, 'sim_number'): 166 sim = True 167 num_sim = cdp.sim_number 168 169 # The inversion matrix. 170 inv_mat = -eye(3) 171 172 # Create the structural object. 173 structure = Internal() 174 175 # Create model for the positive and negative images. 176 if cdp.model in ['rotor', 'free rotor']: 177 neg_cone = False 178 model = structure.add_model(model=1) 179 model_num = 1 180 if neg_cone: 181 model_neg = structure.add_model(model=2) 182 model_num = 2 183 184 # The pivot point. 185 pivot = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) 186 187 # The rotor object. 188 if cdp.model in ['rotor', 'free rotor', 'iso cone', 'iso cone, free rotor', 'pseudo-ellipse']: 189 # The rotor angle. 190 if cdp.model in ['free rotor', 'iso cone, free rotor']: 191 rotor_angle = pi 192 else: 193 rotor_angle = cdp.cone_sigma_max 194 195 # Get the CoM of the entire molecule to use as the centre of the rotor. 196 com = pivot 197 if cdp.model in ['rotor', 'free rotor']: 198 com = pipe_centre_of_mass(verbosity=0) 199 200 # Generate the rotor axis. 201 if cdp.model in ['rotor']: 202 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com) 203 elif cdp.model in ['free rotor', 'iso cone', 'iso cone, free rotor']: 204 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi) 205 else: 206 axis = create_rotor_axis_euler(alpha=cdp.eigen_alpha, beta=cdp.eigen_beta, gamma=cdp.eigen_gamma) 207 208 # The size of the rotor, taking the 30 Angstrom cone representation into account. 209 if cdp.model in ['rotor', 'free rotor']: 210 span = 20e-10 211 else: 212 span = 35e-10 213 214 # Add the rotor object to the structure as a new molecule. 215 rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis, axis_pt=pivot, centre=com, span=span, blade_length=5e-10, staggered=True) 216 217 # Create the molecule. 218 structure.add_molecule(name='rest') 219 220 # Alias the molecules. 221 mol = model.mol[0] 222 if neg_cone: 223 mol_neg = model_neg.mol[0] 224 225 226 # The pivot point. 227 ################## 228 229 # Add the pivot point. 230 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') 231 232 233 # The axes. 234 ########### 235 236 # The spherical angles. 237 if cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 238 # Print out. 239 print("\nGenerating the z-axis system.") 240 241 # The axis. 242 if cdp.model in ['rotor']: 243 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com) 244 else: 245 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi) 246 print(("Central axis: %s." % axis)) 247 248 # Rotations and inversions. 249 axis_pos = axis 250 axis_neg = dot(inv_mat, axis) 251 252 # Simulation central axis. 253 axis_sim_pos = None 254 axis_sim_neg = None 255 if sim: 256 # Init. 257 axis_sim = zeros((cdp.sim_number, 3), float64) 258 259 # Fill the structure. 260 for i in range(cdp.sim_number): 261 if cdp.model in ['rotor']: 262 axis_sim[i] = create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[i], pivot=pivot, point=com) 263 else: 264 axis_sim[i] = create_rotor_axis_spherical(theta=cdp.axis_theta_sim[i], phi=cdp.axis_phi_sim[i]) 265 266 # Inversion. 267 axis_sim_pos = axis_sim 268 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos))) 269 270 # Generate the axis vectors. 271 print("\nGenerating the axis vectors.") 272 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) 273 274 # The negative. 275 if neg_cone: 276 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) 277 278 # The full axis system. 279 else: 280 # Print out. 281 print("\nGenerating the full axis system.") 282 283 # The axis system. 284 axes = zeros((3, 3), float64) 285 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes) 286 print(("Axis system:\n%s" % axes)) 287 288 # Rotations and inversions. 289 axes_pos = axes 290 axes_neg = dot(inv_mat, axes) 291 292 # Simulations 293 axes_sim_pos = None 294 axes_sim_neg = None 295 if sim: 296 # Init. 297 axes_sim_pos = zeros((cdp.sim_number, 3, 3), float64) 298 axes_sim_neg = zeros((cdp.sim_number, 3, 3), float64) 299 300 # Fill the structure. 301 for i in range(cdp.sim_number): 302 # The positive system. 303 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_pos[i]) 304 305 # The negative system. 306 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_neg[i]) 307 axes_sim_neg[i] = dot(inv_mat, axes_sim_neg[i]) 308 309 # Generate the axis vectors. 310 print("\nGenerating the axis vectors.") 311 label = ['x', 'y', 'z'] 312 for j in range(len(label)): 313 # The simulation data. 314 axis_sim_pos = None 315 axis_sim_neg = None 316 if sim: 317 axis_sim_pos = axes_sim_pos[:,:, j] 318 axis_sim_neg = axes_sim_neg[:,:, j] 319 320 # The vectors. 321 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) 322 if neg_cone: 323 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) 324 325 326 # The cone object. 327 ################## 328 329 # Skip models missing a cone. 330 if cdp.model not in ['rotor', 'free rotor']: 331 # The rotation matrix (rotation from the z-axis to the cone axis). 332 if cdp.model not in ['iso cone', 'iso cone, torsionless', 'iso cone, free rotor']: 333 R = axes 334 else: 335 R = zeros((3, 3), float64) 336 two_vect_to_R(array([0, 0, 1], float64), axis, R) 337 338 # Average position rotation. 339 R_pos = R 340 R_neg = dot(inv_mat, R) 341 342 # The pseudo-ellipse cone object. 343 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 344 cone_obj = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y) 345 346 # The isotropic cone object. 347 else: 348 # The angle. 349 if hasattr(cdp, 'cone_theta'): 350 cone_theta = cdp.cone_theta 351 elif hasattr(cdp, 'cone_s1'): 352 cone_theta = order_parameters.iso_cone_S_to_theta(cdp.cone_s1) 353 354 # The object. 355 cone_obj = Iso_cone(cone_theta) 356 357 # Create the positive and negative cones. 358 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) 359 360 # The negative. 361 if neg_cone: 362 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) 363 364 365 # Create the PDB file. 366 ###################### 367 368 # Print out. 369 print("\nGenerating the PDB file.") 370 371 # Write the file. 372 pdb_file = open_write_file(file, dir, force=force) 373 structure.write_pdb(pdb_file) 374 pdb_file.close()

375 376

377 -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):

378 """Create 3 different PDB files for representing the frame order dynamics of the system. 379 380 @keyword ave_pos_file: The name of the file for the average molecule structure. 381 @type ave_pos_file: str or None 382 @keyword rep_file: The name of the file of the PDB representation of the frame order dynamics to create. 383 @type rep_file: str or None 384 @keyword dist_file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model. 385 @type dist_file: str or None 386 @keyword dir: The name of the directory to place the PDB file into. 387 @type dir: str 388 @keyword size: The size of the geometric object in Angstroms. 389 @type size: float 390 @keyword inc: The number of increments for the filling of the cone objects. 391 @type inc: int 392 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 393 @type force: bool 394 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation. 395 @type neg_cone: bool 396 """ 397 398 # Check that at least one PDB file name is given. 399 if not ave_pos_file and not rep_file and not dist_file: 400 raise RelaxError("Minimally one PDB file name must be supplied.") 401 402 # Test if the current data pipe exists. 403 check_pipe() 404 405 # Create the average position structure. 406 if ave_pos_file: 407 pdb_ave_pos(file=ave_pos_file, dir=dir, force=force) 408 409 # Nothing more to do for the rigid model. 410 if cdp.model == 'rigid': 411 return 412 413 # Create the geometric representation. 414 if rep_file: 415 pdb_geometric_rep(file=rep_file, dir=dir, size=size, inc=inc, force=force, neg_cone=neg_cone) 416 417 # Create the distribution. 418 if dist_file: 419 pdb_distribution(file=dist_file, dir=dir, force=force)

420 421

422 -def pivot(pivot=None, order=1, fix=False):

423 """Set the pivot point for the 2 body motion. 424 425 @keyword pivot: The pivot point of the two bodies (domains, etc.) in the structural coordinate system. 426 @type pivot: list of num 427 @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. 428 @type order: int 429 @keyword fix: A flag specifying if the pivot point should be fixed during optimisation. 430 @type fix: bool 431 """ 432 433 # Test if the current data pipe exists. 434 check_pipe() 435 436 # Store the fixed flag. 437 cdp.pivot_fixed = fix 438 439 # No pivot given, so update the model if needed and quit. 440 if pivot == None: 441 if hasattr(cdp, 'model'): 442 update_model() 443 return 444 445 # Convert the pivot to a numpy array. 446 pivot = array(pivot, float64) 447 448 # Check the pivot validity. 449 is_float_array(pivot, name='pivot point', size=3) 450 451 # Store the pivot point and fixed flag. 452 if order == 1: 453 cdp.pivot_x = pivot[0] 454 cdp.pivot_y = pivot[1] 455 cdp.pivot_z = pivot[2] 456 else: 457 # The variable names. 458 name_x = 'pivot_x_%i' % order 459 name_y = 'pivot_y_%i' % order 460 name_z = 'pivot_z_%i' % order 461 462 # Store the variables. 463 setattr(cdp, name_x, pivot[0]) 464 setattr(cdp, name_y, pivot[1]) 465 setattr(cdp, name_z, pivot[2]) 466 467 # Update the model. 468 if hasattr(cdp, 'model'): 469 update_model()

470 471

472 -def quad_int(flag=False):

473 """Turn the high precision Scipy quadratic numerical integration on or off. 474 475 @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. 476 @type flag: bool 477 """ 478 479 # Test if the current data pipe exists. 480 check_pipe() 481 482 # Store the flag. 483 cdp.quad_int = flag

484 485

486 -def ref_domain(ref=None):

487 """Set the reference domain for the frame order, multi-domain models. 488 489 @param ref: The reference domain. 490 @type ref: str 491 """ 492 493 # Test if the current data pipe exists. 494 check_pipe() 495 496 # Check that the domain is defined. 497 if not hasattr(cdp, 'domain') or ref not in cdp.domain: 498 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % ref) 499 500 # Test if the reference domain exists. 501 exists = False 502 for tensor_cont in cdp.align_tensors: 503 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref: 504 exists = True 505 if not exists: 506 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 507 508 # Set the reference domain. 509 cdp.ref_domain = ref 510 511 # Update the model. 512 if hasattr(cdp, 'model'): 513 update_model()

514 515

516 -def select_model(model=None):

517 """Select the Frame Order model. 518 519 @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'. 520 @type model: str 521 """ 522 523 # Test if the current data pipe exists. 524 check_pipe() 525 526 # Test if the model name exists. 527 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']: 528 raise RelaxError("The model name " + repr(model) + " is invalid.") 529 530 # Set the model 531 cdp.model = model 532 533 # Initialise the list of model parameters. 534 cdp.params = [] 535 536 # Set the integration method if needed. 537 if not hasattr(cdp, 'quad_int'): 538 # Scipy quadratic numerical integration. 539 if cdp.model in []: 540 cdp.quad_int = True 541 542 # Quasi-random numerical integration. 543 else: 544 cdp.quad_int = False 545 546 # Update the model. 547 update_model()

548

Source Code for Module specific_analyses.frame_order.uf