target_functions.frame

58 """Class containing the target function of the optimisation of Frame Order matrix components.""" 59

60 - def __init__(self, model=None, init_params=None, full_tensors=None, full_in_ref_frame=None, rdcs=None, rdc_errors=None, rdc_weights=None, rdc_vect=None, dip_const=None, pcs=None, pcs_errors=None, pcs_weights=None, atomic_pos=None, temp=None, frq=None, paramag_centre=zeros(3), scaling_matrix=None, sobol_max_points=200, sobol_oversample=100, com=None, ave_pos_pivot=zeros(3), pivot=None, pivot_opt=False, quad_int=False):

61 """Set up the target functions for the Frame Order theories. 62 63 @keyword model: The name of the Frame Order model. 64 @type model: str 65 @keyword init_params: The initial parameter values. 66 @type init_params: numpy float64 array 67 @keyword full_tensors: An array of the {Axx, Ayy, Axy, Axz, Ayz} values for all full alignment tensors. The format is [Axx1, Ayy1, Axy1, Axz1, Ayz1, Axx2, Ayy2, Axy2, Axz2, Ayz2, ..., Axxn, Ayyn, Axyn, Axzn, Ayzn]. 68 @type full_tensors: numpy nx5D, rank-1 float64 array 69 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 70 @type full_in_ref_frame: numpy rank-1 array 71 @keyword rdcs: The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 72 @type rdcs: numpy rank-2 array 73 @keyword rdc_errors: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 74 @type rdc_errors: numpy rank-2 array 75 @keyword rdc_weights: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 76 @type rdc_weights: numpy rank-2 array 77 @keyword rdc_vect: The unit XH vector lists corresponding to the RDC values. The first index must correspond to the spin systems and the second index to the x, y, z elements. 78 @type rdc_vect: numpy rank-2 array 79 @keyword dip_const: The dipolar constants for each RDC. The indices correspond to the spin systems j. 80 @type dip_const: numpy rank-1 array 81 @keyword pcs: The PCS lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 82 @type pcs: numpy rank-2 array 83 @keyword pcs_errors: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 84 @type pcs_errors: numpy rank-2 array 85 @keyword pcs_weights: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 86 @type pcs_weights: numpy rank-2 array 87 @keyword atomic_pos: The atomic positions of all spins for the PCS and PRE data. The first index is the spin systems j and the second is the structure or state c. 88 @type atomic_pos: numpy rank-3 array 89 @keyword temp: The temperature of each PCS data set. 90 @type temp: numpy rank-1 array 91 @keyword frq: The frequency of each PCS data set. 92 @type frq: numpy rank-1 array 93 @keyword paramag_centre: The paramagnetic centre position (or positions). 94 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 95 @keyword scaling_matrix: The square and diagonal scaling matrix. 96 @type scaling_matrix: numpy rank-2 array 97 @keyword sobol_max_points: The maximum number of Sobol' points to use for the numerical PCS integration technique. 98 @type sobol_max_points: int 99 @keyword sobol_oversample: The oversampling factor Ov used for the total number of points N * Ov * 10**M, where N is the maximum number of Sobol' points and M is the number of dimensions or torsion-tilt angles for the system. 100 @type sobol_oversample: int 101 @keyword com: The centre of mass of the system. This is used for defining the rotor model systems. 102 @type com: numpy 3D rank-1 array 103 @keyword ave_pos_pivot: The pivot point to rotate all atoms about to the average domain position. In most cases this will be the centre of mass of the moving domain. This pivot is shifted by the translation vector. 104 @type ave_pos_pivot: numpy 3D rank-1 array 105 @keyword pivot: The pivot point for the ball-and-socket joint motion. This is needed if PCS or PRE values are used. 106 @type pivot: numpy rank-1, 3D array or None 107 @keyword pivot_opt: A flag which if True will allow the pivot point of the motion to be optimised. 108 @type pivot_opt: bool 109 @keyword quad_int: A 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. 110 @type quad_int: bool 111 """ 112 113 # Model test. 114 if not model: 115 raise RelaxError("The type of Frame Order model must be specified.") 116 117 # Store the initial parameter (as a copy). 118 self.params = deepcopy(init_params) 119 120 # Store the agrs. 121 self.model = model 122 self.full_tensors = deepcopy(full_tensors) 123 self.full_in_ref_frame = deepcopy(full_in_ref_frame) 124 self.rdc = deepcopy(rdcs) 125 self.rdc_weights = deepcopy(rdc_weights) 126 self.rdc_vect = deepcopy(rdc_vect) 127 self.dip_const = deepcopy(dip_const) 128 self.pcs = deepcopy(pcs) 129 self.pcs_weights = deepcopy(pcs_weights) 130 self.atomic_pos = deepcopy(atomic_pos) 131 self.temp = deepcopy(temp) 132 self.frq = deepcopy(frq) 133 self.total_num_params = len(init_params) 134 self.sobol_max_points = sobol_max_points 135 self.sobol_oversample = sobol_oversample 136 self.com = deepcopy(com) 137 self.pivot_opt = pivot_opt 138 self.quad_int = quad_int 139 140 # Tensor setup. 141 self._init_tensors() 142 143 # Scaling initialisation. 144 self.scaling_matrix = scaling_matrix 145 if self.scaling_matrix is not None: 146 self.scaling_flag = True 147 else: 148 self.scaling_flag = False 149 150 # The total number of alignments. 151 self.num_align = 0 152 if rdcs is not None: 153 self.num_align = len(rdcs) 154 elif pcs is not None: 155 self.num_align = len(pcs) 156 157 # Set the RDC and PCS flags (indicating the presence of data). 158 rdc_flag = [True] * self.num_align 159 pcs_flag = [True] * self.num_align 160 for align_index in range(self.num_align): 161 if rdcs is None or len(rdcs[align_index]) == 0: 162 rdc_flag[align_index] = False 163 if pcs is None or len(pcs[align_index]) == 0: 164 pcs_flag[align_index] = False 165 self.rdc_flag = sum(rdc_flag) 166 self.pcs_flag = sum(pcs_flag) 167 168 # Default translation vector (if not optimised). 169 self._translation_vector = zeros(3, float64) 170 171 # Some checks. 172 if self.rdc_flag and (rdc_vect is None or not len(rdc_vect)): 173 raise RelaxError("The rdc_vect argument " + repr(rdc_vect) + " must be supplied.") 174 if self.pcs_flag and (atomic_pos is None or not len(atomic_pos)): 175 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.") 176 177 # The total number of spins. 178 self.num_spins = 0 179 if self.pcs_flag: 180 self.num_spins = len(pcs[0]) 181 182 # The total number of interatomic connections. 183 self.num_interatom = 0 184 if self.rdc_flag: 185 self.num_interatom = len(rdcs[0]) 186 187 # Create multi-dimensional versions of certain structures for faster calculations. 188 if self.pcs_flag: 189 self.spin_ones_struct = ones(self.num_spins, float64) 190 self.pivot = outer(self.spin_ones_struct, pivot) 191 self.paramag_centre = outer(self.spin_ones_struct, paramag_centre) 192 self.ave_pos_pivot = outer(self.spin_ones_struct, ave_pos_pivot) 193 else: 194 self.pivot = array([pivot]) 195 196 # Set up the alignment data. 197 for align_index in range(self.num_align): 198 to_tensor(self.A_3D[align_index], self.full_tensors[5*align_index:5*align_index+5]) 199 200 # PCS errors. 201 if self.pcs_flag: 202 err = False 203 for i in range(len(pcs_errors)): 204 for j in range(len(pcs_errors[i])): 205 if not isNaN(pcs_errors[i, j]): 206 err = True 207 if err: 208 self.pcs_error = pcs_errors 209 else: 210 # Missing errors (default to 0.1 ppm errors). 211 self.pcs_error = 0.1 * 1e-6 * ones((self.num_align, self.num_spins), float64) 212 213 # RDC errors. 214 if self.rdc_flag: 215 err = False 216 for i in range(len(rdc_errors)): 217 for j in range(len(rdc_errors[i])): 218 if not isNaN(rdc_errors[i, j]): 219 err = True 220 if err: 221 self.rdc_error = rdc_errors 222 else: 223 # Missing errors (default to 1 Hz errors). 224 self.rdc_error = ones((self.num_align, self.num_interatom), float64) 225 226 # Missing data matrices (RDC). 227 if self.rdc_flag: 228 self.missing_rdc = zeros((self.num_align, self.num_interatom), uint8) 229 230 # Missing data matrices (PCS). 231 if self.pcs_flag: 232 self.missing_pcs = zeros((self.num_align, self.num_spins), uint8) 233 234 # Clean up problematic data and put the weights into the errors.. 235 if self.rdc_flag or self.pcs_flag: 236 for align_index in range(self.num_align): 237 # Loop over the RDCs. 238 if self.rdc_flag: 239 for j in range(self.num_interatom): 240 if isNaN(self.rdc[align_index, j]): 241 # Set the flag. 242 self.missing_rdc[align_index, j] = 1 243 244 # Change the NaN to zero. 245 self.rdc[align_index, j] = 0.0 246 247 # Change the error to one (to avoid zero division). 248 self.rdc_error[align_index, j] = 1.0 249 250 # Change the weight to one. 251 rdc_weights[align_index, j] = 1.0 252 253 # The RDC weights. 254 if self.rdc_flag: 255 self.rdc_error[align_index, j] = self.rdc_error[align_index, j] / sqrt(rdc_weights[align_index, j]) 256 257 # Loop over the PCSs. 258 if self.pcs_flag: 259 for j in range(self.num_spins): 260 if isNaN(self.pcs[align_index, j]): 261 # Set the flag. 262 self.missing_pcs[align_index, j] = 1 263 264 # Change the NaN to zero. 265 self.pcs[align_index, j] = 0.0 266 267 # Change the error to one (to avoid zero division). 268 self.pcs_error[align_index, j] = 1.0 269 270 # Change the weight to one. 271 pcs_weights[align_index, j] = 1.0 272 273 # The PCS weights. 274 if self.pcs_flag: 275 self.pcs_error[align_index, j] = self.pcs_error[align_index, j] / sqrt(pcs_weights[align_index, j]) 276 277 # The paramagnetic centre vectors and distances. 278 if self.pcs_flag: 279 # Initialise the data structures. 280 self.paramag_unit_vect = zeros(atomic_pos.shape, float64) 281 self.paramag_dist = zeros(self.num_spins, float64) 282 self.pcs_const = zeros((self.num_align, self.num_spins), float64) 283 self.r_pivot_atom = zeros((self.num_spins, 3), float32) 284 self.r_pivot_atom_rev = zeros((self.num_spins, 3), float32) 285 self.r_ln_pivot = self.pivot - self.paramag_centre 286 287 # Set up the paramagnetic constant (without the interatomic distance and in Angstrom units). 288 for align_index in range(self.num_align): 289 self.pcs_const[align_index] = pcs_constant(self.temp[align_index], self.frq[align_index], 1.0) * 1e30 290 291 # PCS function, gradient, and Hessian matrices. 292 self.pcs_theta = None 293 if self.pcs_flag: 294 self.pcs_theta = zeros((self.num_align, self.num_spins), float64) 295 self.pcs_theta_err = zeros((self.num_align, self.num_spins), float64) 296 self.dpcs_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 297 self.d2pcs_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 298 299 # RDC function, gradient, and Hessian matrices. 300 self.rdc_theta = None 301 if self.rdc_flag: 302 self.rdc_theta = zeros((self.num_align, self.num_interatom), float64) 303 self.drdc_theta = zeros((self.total_num_params, self.num_align, self.num_interatom), float64) 304 self.d2rdc_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_interatom), float64) 305 306 # The target function extension. 307 ext = '_qr_int' 308 if self.quad_int: 309 ext = '_quad_int' 310 311 # Non-numerical models. 312 if model in [MODEL_RIGID]: 313 if model == MODEL_RIGID: 314 self.func = self.func_rigid 315 316 # The Sobol' sequence data and target function aliases. 317 else: 318 if model == MODEL_PSEUDO_ELLIPSE: 319 self.create_sobol_data(dims=['theta', 'phi', 'sigma']) 320 self.func = getattr(self, 'func_pseudo_ellipse'+ext) 321 elif model == MODEL_PSEUDO_ELLIPSE_TORSIONLESS: 322 self.create_sobol_data(dims=['theta', 'phi']) 323 self.func = getattr(self, 'func_pseudo_ellipse_torsionless'+ext) 324 elif model == MODEL_PSEUDO_ELLIPSE_FREE_ROTOR: 325 self.create_sobol_data(dims=['theta', 'phi', 'sigma']) 326 self.func = getattr(self, 'func_pseudo_ellipse_free_rotor'+ext) 327 elif model == MODEL_ISO_CONE: 328 self.create_sobol_data(dims=['theta', 'phi', 'sigma']) 329 self.func = getattr(self, 'func_iso_cone'+ext) 330 elif model == MODEL_ISO_CONE_TORSIONLESS: 331 self.create_sobol_data(dims=['theta', 'phi']) 332 self.func = getattr(self, 'func_iso_cone_torsionless'+ext) 333 elif model == MODEL_ISO_CONE_FREE_ROTOR: 334 self.create_sobol_data(dims=['theta', 'phi', 'sigma']) 335 self.func = getattr(self, 'func_iso_cone_free_rotor'+ext) 336 elif model == MODEL_ROTOR: 337 self.create_sobol_data(dims=['sigma']) 338 self.func = getattr(self, 'func_rotor'+ext) 339 elif model == MODEL_FREE_ROTOR: 340 self.create_sobol_data(dims=['sigma']) 341 self.func = getattr(self, 'func_free_rotor'+ext) 342 elif model == MODEL_DOUBLE_ROTOR: 343 self.create_sobol_data(dims=['sigma', 'sigma2']) 344 self.func = getattr(self, 'func_double_rotor'+ext)

345 346

347 - def _init_tensors(self):

348 """Set up isotropic cone optimisation against the alignment tensor data.""" 349 350 # Some checks. 351 if self.full_tensors is None or not len(self.full_tensors): 352 raise RelaxError("The full_tensors argument " + repr(self.full_tensors) + " must be supplied.") 353 if self.full_in_ref_frame is None or not len(self.full_in_ref_frame): 354 raise RelaxError("The full_in_ref_frame argument " + repr(self.full_in_ref_frame) + " must be supplied.") 355 356 # Tensor set up. 357 self.num_tensors = int(len(self.full_tensors) / 5) 358 self.A_3D = zeros((self.num_tensors, 3, 3), float64) 359 self.A_3D_bc = zeros((self.num_tensors, 3, 3), float64) 360 self.A_5D_bc = zeros(self.num_tensors*5, float64) 361 362 # The rotation to the Frame Order eigenframe. 363 self.R_eigen = zeros((3, 3), float64) 364 self.R_eigen_2 = zeros((3, 3), float64) 365 self.R_ave = zeros((3, 3), float64) 366 self.Ri_prime = zeros((3, 3), float64) 367 self.Ri2_prime = zeros((3, 3), float64) 368 self.tensor_3D = zeros((3, 3), float64) 369 370 # The cone axis storage and molecular frame z-axis. 371 self.cone_axis = zeros(3, float64) 372 self.z_axis = array([0, 0, 1], float64) 373 374 # The rotor axes. 375 self.rotor_axis = zeros(3, float64) 376 self.rotor_axis_2 = zeros(3, float64) 377 378 # Initialise the Frame Order matrices. 379 self.frame_order_2nd = zeros((9, 9), float64) 380 381 # A rotation matrix for general use. 382 self.R = zeros((3, 3), float64)

383 384

385 - def func_double_rotor_qr_int(self, params):

386 """Quasi-random Sobol' integration target function for the double rotor model. 387 388 This function optimises the model parameters using the RDC and PCS base data. Quasi-random, Sobol' sequence based, numerical integration is used for the PCS. 389 390 391 @param params: The vector of parameter values. These can include {pivot_x, pivot_y, pivot_z, pivot_disp, ave_pos_x, ave_pos_y, ave_pos_z, ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_sigma_max, cone_sigma_max_2}. 392 @type params: list of float 393 @return: The chi-squared or SSE value. 394 @rtype: float 395 """ 396 397 # Scaling. 398 if self.scaling_flag: 399 params = dot(params, self.scaling_matrix) 400 401 # Unpack the parameters. 402 if self.pivot_opt: 403 pivot2 = outer(self.spin_ones_struct, params[:3]) 404 param_disp = params[3] 405 self._translation_vector = params[4:7] 406 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, sigma_max, sigma_max_2 = params[7:] 407 else: 408 pivot2 = self.pivot 409 param_disp = params[0] 410 self._translation_vector = params[1:4] 411 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, sigma_max, sigma_max_2 = params[4:] 412 413 # Reconstruct the full eigenframe of the motion. 414 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 415 416 # The Kronecker product of the eigenframe. 417 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 418 419 # Generate the 2nd degree Frame Order super matrix. 420 frame_order_2nd = compile_2nd_matrix_double_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max, sigma_max_2) 421 422 # The average frame rotation matrix (and reduce and rotate the tensors). 423 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 424 425 # Pre-transpose matrices for faster calculations. 426 RT_eigen = transpose(self.R_eigen) 427 RT_ave = transpose(self.R_ave) 428 429 # Pre-calculate all the necessary vectors. 430 if self.pcs_flag: 431 # The 1st pivot point (sum of the 2nd pivot and the displacement along the eigenframe z-axis). 432 pivot = pivot2 + param_disp * self.R_eigen[:, 2] 433 434 # Calculate the vectors. 435 self.calc_vectors(pivot=pivot, pivot2=pivot2, R_ave=self.R_ave, RT_ave=RT_ave) 436 437 # Initial chi-squared (or SSE) value. 438 chi2_sum = 0.0 439 440 # RDCs. 441 if self.rdc_flag: 442 # Loop over each alignment. 443 for align_index in range(self.num_align): 444 # Loop over the RDCs. 445 for j in range(self.num_interatom): 446 # The back calculated RDC. 447 if not self.missing_rdc[align_index, j]: 448 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 449 450 # Calculate and sum the single alignment chi-squared value (for the RDC). 451 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 452 453 # PCS via numerical integration. 454 if self.pcs_flag: 455 # Numerical integration of the PCSs. 456 pcs_numeric_qr_int_double_rotor(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, sigma_max=sigma_max, sigma_max_2=sigma_max_2, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, r_inter_pivot=self.r_inter_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, Ri2_prime=sobol_data.Ri2_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 457 458 # Calculate and sum the single alignment chi-squared value (for the PCS). 459 for align_index in range(self.num_align): 460 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 461 462 # Return the chi-squared value. 463 return chi2_sum

464 465

466 - def func_double_rotor_quad_int(self, params):

467 """SciPy quadratic integration target function for the double rotor model. 468 469 This function optimises the model parameters using the RDC and PCS base data. Quasi-random, Sobol' sequence based, numerical integration is used for the PCS. 470 471 472 @param params: The vector of parameter values. These can include {pivot_x, pivot_y, pivot_z, pivot_disp, ave_pos_x, ave_pos_y, ave_pos_z, ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_sigma_max, cone_sigma_max_2}. 473 @type params: list of float 474 @return: The chi-squared or SSE value. 475 @rtype: float 476 """ 477 478 # Scaling. 479 if self.scaling_flag: 480 params = dot(params, self.scaling_matrix) 481 482 # Unpack the parameters. 483 if self.pivot_opt: 484 pivot2 = outer(self.spin_ones_struct, params[:3]) 485 param_disp = params[3] 486 self._translation_vector = params[4:7] 487 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, sigma_max, sigma_max_2 = params[7:] 488 else: 489 pivot2 = self.pivot 490 param_disp = params[0] 491 self._translation_vector = params[1:4] 492 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, sigma_max, sigma_max_2 = params[4:] 493 494 # Reconstruct the full eigenframe of the motion. 495 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 496 497 # The Kronecker product of the eigenframe. 498 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 499 500 # Generate the 2nd degree Frame Order super matrix. 501 frame_order_2nd = compile_2nd_matrix_double_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max, sigma_max_2) 502 503 # The average frame rotation matrix (and reduce and rotate the tensors). 504 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 505 506 # Pre-transpose matrices for faster calculations. 507 RT_eigen = transpose(self.R_eigen) 508 RT_ave = transpose(self.R_ave) 509 510 # Pre-calculate all the necessary vectors. 511 if self.pcs_flag: 512 # The 1st pivot point (sum of the 2nd pivot and the displacement along the eigenframe z-axis). 513 pivot = pivot2 + param_disp * self.R_eigen[:, 2] 514 515 # Calculate the vectors. 516 self.calc_vectors(pivot=pivot, pivot2=pivot2, R_ave=self.R_ave, RT_ave=RT_ave) 517 518 # Initial chi-squared (or SSE) value. 519 chi2_sum = 0.0 520 521 # RDCs. 522 if self.rdc_flag: 523 # Loop over each alignment. 524 for align_index in range(self.num_align): 525 # Loop over the RDCs. 526 for j in range(self.num_interatom): 527 # The back calculated RDC. 528 if not self.missing_rdc[align_index, j]: 529 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 530 531 # Calculate and sum the single alignment chi-squared value (for the RDC). 532 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 533 534 # PCS via numerical integration. 535 if self.pcs_flag: 536 # Loop over each alignment. 537 for align_index in range(self.num_align): 538 # Loop over the PCSs. 539 for j in range(self.num_spins): 540 # The back calculated PCS. 541 if not self.missing_pcs[align_index, j]: 542 # Forwards and reverse rotations. 543 if self.full_in_ref_frame[align_index]: 544 r_pivot_atom = self.r_pivot_atom[j] 545 else: 546 r_pivot_atom = self.r_pivot_atom_rev[j] 547 548 # The numerical integration. 549 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_double_rotor(sigma_max=sigma_max, sigma_max_2=sigma_max_2, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], r_inter_pivot=self.r_inter_pivot[j], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime, Ri2_prime=self.Ri2_prime) 550 551 # Calculate and sum the single alignment chi-squared value (for the PCS). 552 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 553 554 # Return the chi-squared value. 555 return chi2_sum

556 557

558 - def func_free_rotor_qr_int(self, params):

559 """Quasi-random Sobol' integration target function for the free rotor model. 560 561 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 562 563 564 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi}. 565 @type params: list of float 566 @return: The chi-squared or SSE value. 567 @rtype: float 568 """ 569 570 # Scaling. 571 if self.scaling_flag: 572 params = dot(params, self.scaling_matrix) 573 574 # Unpack the parameters. 575 if self.pivot_opt: 576 pivot = outer(self.spin_ones_struct, params[:3]) 577 self._translation_vector = params[3:6] 578 ave_pos_beta, ave_pos_gamma, axis_alpha = params[6:] 579 else: 580 pivot = self.pivot 581 self._translation_vector = params[:3] 582 ave_pos_beta, ave_pos_gamma, axis_alpha = params[3:] 583 584 # Generate the rotor axis. 585 self.cone_axis = create_rotor_axis_alpha(alpha=axis_alpha, pivot=pivot[0], point=self.com) 586 587 # Pre-calculate the eigenframe rotation matrix. 588 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 589 590 # The Kronecker product of the eigenframe rotation. 591 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 592 593 # Generate the 2nd degree Frame Order super matrix. 594 frame_order_2nd = compile_2nd_matrix_free_rotor(self.frame_order_2nd, Rx2_eigen) 595 596 # Reduce and rotate the tensors. 597 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 598 599 # Pre-transpose matrices for faster calculations. 600 RT_eigen = transpose(self.R_eigen) 601 RT_ave = transpose(self.R_ave) 602 603 # Pre-calculate all the necessary vectors. 604 if self.pcs_flag: 605 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 606 607 # Initial chi-squared (or SSE) value. 608 chi2_sum = 0.0 609 610 # RDCs. 611 if self.rdc_flag: 612 # Loop over each alignment. 613 for align_index in range(self.num_align): 614 # Loop over the RDCs. 615 for j in range(self.num_interatom): 616 # The back calculated RDC. 617 if not self.missing_rdc[align_index, j]: 618 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 619 620 # Calculate and sum the single alignment chi-squared value (for the RDC). 621 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 622 623 # PCS via numerical integration. 624 if self.pcs_flag: 625 # Numerical integration of the PCSs. 626 pcs_numeric_qr_int_rotor(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, sigma_max=pi, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 627 628 # Calculate and sum the single alignment chi-squared value (for the PCS). 629 for align_index in range(self.num_align): 630 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 631 632 633 # Return the chi-squared value. 634 return chi2_sum

635 636

637 - def func_free_rotor_quad_int(self, params):

638 """SciPy quadratic integration target function for the free rotor model. 639 640 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 641 642 643 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi}. 644 @type params: list of float 645 @return: The chi-squared or SSE value. 646 @rtype: float 647 """ 648 649 # Scaling. 650 if self.scaling_flag: 651 params = dot(params, self.scaling_matrix) 652 653 # Unpack the parameters. 654 if self.pivot_opt: 655 pivot = outer(self.spin_ones_struct, params[:3]) 656 self._translation_vector = params[3:6] 657 ave_pos_beta, ave_pos_gamma, axis_alpha = params[6:] 658 else: 659 pivot = self.pivot 660 self._translation_vector = params[:3] 661 ave_pos_beta, ave_pos_gamma, axis_alpha = params[3:] 662 663 # Generate the rotor axis. 664 self.cone_axis = create_rotor_axis_alpha(alpha=axis_alpha, pivot=pivot[0], point=self.com) 665 666 # Pre-calculate the eigenframe rotation matrix. 667 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 668 669 # The Kronecker product of the eigenframe rotation. 670 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 671 672 # Generate the 2nd degree Frame Order super matrix. 673 frame_order_2nd = compile_2nd_matrix_free_rotor(self.frame_order_2nd, Rx2_eigen) 674 675 # Reduce and rotate the tensors. 676 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 677 678 # Pre-transpose matrices for faster calculations. 679 RT_eigen = transpose(self.R_eigen) 680 RT_ave = transpose(self.R_ave) 681 682 # Pre-calculate all the necessary vectors. 683 if self.pcs_flag: 684 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 685 686 # Initial chi-squared (or SSE) value. 687 chi2_sum = 0.0 688 689 # RDCs. 690 if self.rdc_flag: 691 # Loop over each alignment. 692 for align_index in range(self.num_align): 693 # Loop over the RDCs. 694 for j in range(self.num_interatom): 695 # The back calculated RDC. 696 if not self.missing_rdc[align_index, j]: 697 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 698 699 # Calculate and sum the single alignment chi-squared value (for the RDC). 700 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 701 702 # PCS via numerical integration. 703 if self.pcs_flag: 704 # Loop over each alignment. 705 for align_index in range(self.num_align): 706 # Loop over the PCSs. 707 for j in range(self.num_spins): 708 # The back calculated PCS. 709 if not self.missing_pcs[align_index, j]: 710 # Forwards and reverse rotations. 711 if self.full_in_ref_frame[align_index]: 712 r_pivot_atom = self.r_pivot_atom[j] 713 else: 714 r_pivot_atom = self.r_pivot_atom_rev[j] 715 716 # The numerical integration. 717 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_rotor(sigma_max=pi, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 718 719 # Calculate and sum the single alignment chi-squared value (for the PCS). 720 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 721 722 # Return the chi-squared value. 723 return chi2_sum

724 725

726 - def func_iso_cone_qr_int(self, params):

727 """Quasi-random Sobol' integration target function for the isotropic cone model. 728 729 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 730 731 732 @param params: The vector of parameter values {alpha, beta, gamma, theta, phi, cone_theta, sigma_max} where the first 3 are the tensor rotation Euler angles, the next two are the polar and azimuthal angles of the cone axis, cone_theta is the cone opening half angle, and sigma_max is the torsion angle. 733 @type params: list of float 734 @return: The chi-squared or SSE value. 735 @rtype: float 736 """ 737 738 # Scaling. 739 if self.scaling_flag: 740 params = dot(params, self.scaling_matrix) 741 742 # Unpack the parameters. 743 if self.pivot_opt: 744 pivot = outer(self.spin_ones_struct, params[:3]) 745 self._translation_vector = params[3:6] 746 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[6:] 747 else: 748 pivot = self.pivot 749 self._translation_vector = params[:3] 750 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 751 752 # Generate the cone axis from the spherical angles. 753 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 754 755 # Pre-calculate the eigenframe rotation matrix. 756 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 757 758 # The Kronecker product of the eigenframe rotation. 759 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 760 761 # Generate the 2nd degree Frame Order super matrix. 762 frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd, Rx2_eigen, cone_theta, sigma_max) 763 764 # Reduce and rotate the tensors. 765 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 766 767 # Pre-transpose matrices for faster calculations. 768 RT_eigen = transpose(self.R_eigen) 769 RT_ave = transpose(self.R_ave) 770 771 # Pre-calculate all the necessary vectors. 772 if self.pcs_flag: 773 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 774 775 # Initial chi-squared (or SSE) value. 776 chi2_sum = 0.0 777 778 # RDCs. 779 if self.rdc_flag: 780 # Loop over each alignment. 781 for align_index in range(self.num_align): 782 # Loop over the RDCs. 783 for j in range(self.num_interatom): 784 # The back calculated RDC. 785 if not self.missing_rdc[align_index, j]: 786 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 787 788 # Calculate and sum the single alignment chi-squared value (for the RDC). 789 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 790 791 # PCS via numerical integration. 792 if self.pcs_flag: 793 # Numerical integration of the PCSs. 794 pcs_numeric_qr_int_iso_cone(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, theta_max=cone_theta, sigma_max=sigma_max, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 795 796 # Calculate and sum the single alignment chi-squared value (for the PCS). 797 for align_index in range(self.num_align): 798 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 799 800 # Return the chi-squared value. 801 return chi2_sum

802 803

804 - def func_iso_cone_quad_int(self, params):

805 """SciPy quadratic integration target function for the isotropic cone model. 806 807 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 808 809 810 @param params: The vector of parameter values {beta, gamma, theta, phi, s1} where the first 2 are the tensor rotation Euler angles, the next two are the polar and azimuthal angles of the cone axis, and s1 is the isotropic cone order parameter. 811 @type params: list of float 812 @return: The chi-squared or SSE value. 813 @rtype: float 814 """ 815 816 # Scaling. 817 if self.scaling_flag: 818 params = dot(params, self.scaling_matrix) 819 820 # Unpack the parameters. 821 if self.pivot_opt: 822 pivot = outer(self.spin_ones_struct, params[:3]) 823 self._translation_vector = params[3:6] 824 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[6:] 825 else: 826 pivot = self.pivot 827 self._translation_vector = params[:3] 828 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 829 830 # Generate the cone axis from the spherical angles. 831 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 832 833 # Pre-calculate the eigenframe rotation matrix. 834 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 835 836 # The Kronecker product of the eigenframe rotation. 837 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 838 839 # Generate the 2nd degree Frame Order super matrix. 840 frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd, Rx2_eigen, cone_theta, sigma_max) 841 842 # Reduce and rotate the tensors. 843 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 844 845 # Pre-transpose matrices for faster calculations. 846 RT_eigen = transpose(self.R_eigen) 847 RT_ave = transpose(self.R_ave) 848 849 # Pre-calculate all the necessary vectors. 850 if self.pcs_flag: 851 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 852 853 # Initial chi-squared (or SSE) value. 854 chi2_sum = 0.0 855 856 # RDCs. 857 if self.rdc_flag: 858 # Loop over each alignment. 859 for align_index in range(self.num_align): 860 # Loop over the RDCs. 861 for j in range(self.num_interatom): 862 # The back calculated RDC. 863 if not self.missing_rdc[align_index, j]: 864 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 865 866 # Calculate and sum the single alignment chi-squared value (for the RDC). 867 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 868 869 # PCS via numerical integration. 870 if self.pcs_flag: 871 # Loop over each alignment. 872 for align_index in range(self.num_align): 873 # Loop over the PCSs. 874 for j in range(self.num_spins): 875 # The back calculated PCS. 876 if not self.missing_pcs[align_index, j]: 877 # Forwards and reverse rotations. 878 if self.full_in_ref_frame[align_index]: 879 r_pivot_atom = self.r_pivot_atom[j] 880 else: 881 r_pivot_atom = self.r_pivot_atom_rev[j] 882 883 # The numerical integration. 884 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_iso_cone(theta_max=cone_theta, sigma_max=sigma_max, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 885 886 # Calculate and sum the single alignment chi-squared value (for the PCS). 887 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 888 889 # Return the chi-squared value. 890 return chi2_sum

891 892

893 - def func_iso_cone_free_rotor_qr_int(self, params):

894 """Quasi-random Sobol' integration target function for the free rotor isotropic cone model. 895 896 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 897 898 899 @param params: The vector of parameter values {beta, gamma, theta, phi, s1} where the first 2 are the tensor rotation Euler angles, the next two are the polar and azimuthal angles of the cone axis, and s1 is the isotropic cone order parameter. 900 @type params: list of float 901 @return: The chi-squared or SSE value. 902 @rtype: float 903 """ 904 905 # Scaling. 906 if self.scaling_flag: 907 params = dot(params, self.scaling_matrix) 908 909 # Unpack the parameters. 910 if self.pivot_opt: 911 pivot = outer(self.spin_ones_struct, params[:3]) 912 self._translation_vector = params[3:6] 913 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, theta_max = params[6:] 914 else: 915 pivot = self.pivot 916 self._translation_vector = params[:3] 917 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, theta_max = params[3:] 918 919 # Generate the cone axis from the spherical angles. 920 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 921 922 # Pre-calculate the eigenframe rotation matrix. 923 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 924 925 # The Kronecker product of the eigenframe rotation. 926 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 927 928 # Generate the 2nd degree Frame Order super matrix. 929 frame_order_2nd = compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, Rx2_eigen, theta_max) 930 931 # Reduce and rotate the tensors. 932 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 933 934 # Pre-transpose matrices for faster calculations. 935 RT_eigen = transpose(self.R_eigen) 936 RT_ave = transpose(self.R_ave) 937 938 # Pre-calculate all the necessary vectors. 939 if self.pcs_flag: 940 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 941 942 # Initial chi-squared (or SSE) value. 943 chi2_sum = 0.0 944 945 # RDCs. 946 if self.rdc_flag: 947 # Loop over each alignment. 948 for align_index in range(self.num_align): 949 # Loop over the RDCs. 950 for j in range(self.num_interatom): 951 # The back calculated RDC. 952 if not self.missing_rdc[align_index, j]: 953 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 954 955 # Calculate and sum the single alignment chi-squared value (for the RDC). 956 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 957 958 # PCS via numerical integration. 959 if self.pcs_flag: 960 # Numerical integration of the PCSs. 961 pcs_numeric_qr_int_iso_cone(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, theta_max=theta_max, sigma_max=pi, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 962 963 # Calculate and sum the single alignment chi-squared value (for the PCS). 964 for align_index in range(self.num_align): 965 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 966 967 # Return the chi-squared value. 968 return chi2_sum

969 970

971 - def func_iso_cone_free_rotor_quad_int(self, params):

972 """SciPy quadratic integration target function for the free rotor isotropic cone model. 973 974 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 975 976 977 @param params: The vector of parameter values {beta, gamma, theta, phi, s1} where the first 2 are the tensor rotation Euler angles, the next two are the polar and azimuthal angles of the cone axis, and s1 is the isotropic cone order parameter. 978 @type params: list of float 979 @return: The chi-squared or SSE value. 980 @rtype: float 981 """ 982 983 # Scaling. 984 if self.scaling_flag: 985 params = dot(params, self.scaling_matrix) 986 987 # Unpack the parameters. 988 if self.pivot_opt: 989 pivot = outer(self.spin_ones_struct, params[:3]) 990 self._translation_vector = params[3:6] 991 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, theta_max = params[6:] 992 else: 993 pivot = self.pivot 994 self._translation_vector = params[:3] 995 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, theta_max = params[3:] 996 997 # Generate the cone axis from the spherical angles. 998 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 999 1000 # Pre-calculate the eigenframe rotation matrix. 1001 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1002 1003 # The Kronecker product of the eigenframe rotation. 1004 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1005 1006 # Generate the 2nd degree Frame Order super matrix. 1007 frame_order_2nd = compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, Rx2_eigen, theta_max) 1008 1009 # Reduce and rotate the tensors. 1010 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1011 1012 # Pre-transpose matrices for faster calculations. 1013 RT_eigen = transpose(self.R_eigen) 1014 RT_ave = transpose(self.R_ave) 1015 1016 # Pre-calculate all the necessary vectors. 1017 if self.pcs_flag: 1018 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1019 1020 # Initial chi-squared (or SSE) value. 1021 chi2_sum = 0.0 1022 1023 # RDCs. 1024 if self.rdc_flag: 1025 # Loop over each alignment. 1026 for align_index in range(self.num_align): 1027 # Loop over the RDCs. 1028 for j in range(self.num_interatom): 1029 # The back calculated RDC. 1030 if not self.missing_rdc[align_index, j]: 1031 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1032 1033 # Calculate and sum the single alignment chi-squared value (for the RDC). 1034 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1035 1036 # PCS via numerical integration. 1037 if self.pcs_flag: 1038 # Loop over each alignment. 1039 for align_index in range(self.num_align): 1040 # Loop over the PCSs. 1041 for j in range(self.num_spins): 1042 # The back calculated PCS. 1043 if not self.missing_pcs[align_index, j]: 1044 # Forwards and reverse rotations. 1045 if self.full_in_ref_frame[align_index]: 1046 r_pivot_atom = self.r_pivot_atom[j] 1047 else: 1048 r_pivot_atom = self.r_pivot_atom_rev[j] 1049 1050 # The numerical integration. 1051 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_iso_cone(theta_max=theta_max, sigma_max=pi, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 1052 1053 # Calculate and sum the single alignment chi-squared value (for the PCS). 1054 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1055 1056 # Return the chi-squared value. 1057 return chi2_sum

1058 1059

1060 - def func_iso_cone_torsionless_qr_int(self, params):

1061 """Quasi-random Sobol' integration target function for the torsionless isotropic cone model. 1062 1063 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1064 1065 1066 @param params: The vector of parameter values {beta, gamma, theta, phi, cone_theta} where the first 2 are the tensor rotation Euler angles, the next two are the polar and azimuthal angles of the cone axis, and cone_theta is cone opening angle. 1067 @type params: list of float 1068 @return: The chi-squared or SSE value. 1069 @rtype: float 1070 """ 1071 1072 # Scaling. 1073 if self.scaling_flag: 1074 params = dot(params, self.scaling_matrix) 1075 1076 # Unpack the parameters. 1077 if self.pivot_opt: 1078 pivot = outer(self.spin_ones_struct, params[:3]) 1079 self._translation_vector = params[3:6] 1080 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[6:] 1081 else: 1082 pivot = self.pivot 1083 self._translation_vector = params[:3] 1084 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1085 1086 # Generate the cone axis from the spherical angles. 1087 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1088 1089 # Pre-calculate the eigenframe rotation matrix. 1090 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1091 1092 # The Kronecker product of the eigenframe rotation. 1093 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1094 1095 # Generate the 2nd degree Frame Order super matrix. 1096 frame_order_2nd = compile_2nd_matrix_iso_cone_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta) 1097 1098 # Reduce and rotate the tensors. 1099 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1100 1101 # Pre-transpose matrices for faster calculations. 1102 RT_eigen = transpose(self.R_eigen) 1103 RT_ave = transpose(self.R_ave) 1104 1105 # Pre-calculate all the necessary vectors. 1106 if self.pcs_flag: 1107 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1108 1109 # Initial chi-squared (or SSE) value. 1110 chi2_sum = 0.0 1111 1112 # RDCs. 1113 if self.rdc_flag: 1114 # Loop over each alignment. 1115 for align_index in range(self.num_align): 1116 # Loop over the RDCs. 1117 for j in range(self.num_interatom): 1118 # The back calculated RDC. 1119 if not self.missing_rdc[align_index, j]: 1120 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1121 1122 # Calculate and sum the single alignment chi-squared value (for the RDC). 1123 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1124 1125 # PCS via numerical integration. 1126 if self.pcs_flag: 1127 # Numerical integration of the PCSs. 1128 pcs_numeric_qr_int_iso_cone_torsionless(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, theta_max=cone_theta, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 1129 1130 # Calculate and sum the single alignment chi-squared value (for the PCS). 1131 for align_index in range(self.num_align): 1132 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1133 1134 # Return the chi-squared value. 1135 return chi2_sum

1136 1137

1138 - def func_iso_cone_torsionless_quad_int(self, params):

1139 """SciPy quadratic integration target function for the torsionless isotropic cone model. 1140 1141 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1142 1143 1144 @param params: The vector of parameter values {beta, gamma, theta, phi, cone_theta} where the first 2 are the tensor rotation Euler angles, the next two are the polar and azimuthal angles of the cone axis, and cone_theta is cone opening angle. 1145 @type params: list of float 1146 @return: The chi-squared or SSE value. 1147 @rtype: float 1148 """ 1149 1150 # Scaling. 1151 if self.scaling_flag: 1152 params = dot(params, self.scaling_matrix) 1153 1154 # Unpack the parameters. 1155 if self.pivot_opt: 1156 pivot = outer(self.spin_ones_struct, params[:3]) 1157 self._translation_vector = params[3:6] 1158 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[6:] 1159 else: 1160 pivot = self.pivot 1161 self._translation_vector = params[:3] 1162 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1163 1164 # Generate the cone axis from the spherical angles. 1165 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1166 1167 # Pre-calculate the eigenframe rotation matrix. 1168 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1169 1170 # The Kronecker product of the eigenframe rotation. 1171 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1172 1173 # Generate the 2nd degree Frame Order super matrix. 1174 frame_order_2nd = compile_2nd_matrix_iso_cone_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta) 1175 1176 # Reduce and rotate the tensors. 1177 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1178 1179 # Pre-transpose matrices for faster calculations. 1180 RT_eigen = transpose(self.R_eigen) 1181 RT_ave = transpose(self.R_ave) 1182 1183 # Pre-calculate all the necessary vectors. 1184 if self.pcs_flag: 1185 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1186 1187 # Initial chi-squared (or SSE) value. 1188 chi2_sum = 0.0 1189 1190 # RDCs. 1191 if self.rdc_flag: 1192 # Loop over each alignment. 1193 for align_index in range(self.num_align): 1194 # Loop over the RDCs. 1195 for j in range(self.num_interatom): 1196 # The back calculated RDC. 1197 if not self.missing_rdc[align_index, j]: 1198 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1199 1200 # Calculate and sum the single alignment chi-squared value (for the RDC). 1201 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1202 1203 # PCS via numerical integration. 1204 if self.pcs_flag: 1205 # Loop over each alignment. 1206 for align_index in range(self.num_align): 1207 # Loop over the PCSs. 1208 for j in range(self.num_spins): 1209 # The back calculated PCS. 1210 if not self.missing_pcs[align_index, j]: 1211 # Forwards and reverse rotations. 1212 if self.full_in_ref_frame[align_index]: 1213 r_pivot_atom = self.r_pivot_atom[j] 1214 else: 1215 r_pivot_atom = self.r_pivot_atom_rev[j] 1216 1217 # The numerical integration. 1218 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_iso_cone_torsionless(theta_max=cone_theta, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 1219 1220 # Calculate and sum the single alignment chi-squared value (for the PCS). 1221 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1222 1223 # Return the chi-squared value. 1224 return chi2_sum

1225 1226

1227 - def func_pseudo_ellipse_qr_int(self, params):

1228 """Quasi-random Sobol' integration target function for the pseudo-ellipse model. 1229 1230 This function optimises the model parameters using the RDC and PCS base data. Quasi-random, Sobol' sequence based, numerical integration is used for the PCS. 1231 1232 1233 @param params: The vector of parameter values {alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max} where the first 3 are the average position rotation Euler angles, the next 3 are the Euler angles defining the eigenframe, and the last 3 are the pseudo-elliptic cone geometric parameters. 1234 @type params: list of float 1235 @return: The chi-squared or SSE value. 1236 @rtype: float 1237 """ 1238 1239 # Scaling. 1240 if self.scaling_flag: 1241 params = dot(params, self.scaling_matrix) 1242 1243 # Unpack the parameters. 1244 if self.pivot_opt: 1245 pivot = outer(self.spin_ones_struct, params[:3]) 1246 self._translation_vector = params[3:6] 1247 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = params[6:] 1248 else: 1249 pivot = self.pivot 1250 self._translation_vector = params[:3] 1251 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = params[3:] 1252 1253 # Reconstruct the full eigenframe of the motion. 1254 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1255 1256 # The Kronecker product of the eigenframe. 1257 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1258 1259 # Generate the 2nd degree Frame Order super matrix. 1260 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y, cone_sigma_max) 1261 1262 # Reduce and rotate the tensors. 1263 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1264 1265 # Pre-transpose matrices for faster calculations. 1266 RT_eigen = transpose(self.R_eigen) 1267 RT_ave = transpose(self.R_ave) 1268 1269 # Pre-calculate all the necessary vectors. 1270 if self.pcs_flag: 1271 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1272 1273 # Initial chi-squared (or SSE) value. 1274 chi2_sum = 0.0 1275 1276 # RDCs. 1277 if self.rdc_flag: 1278 # Loop over each alignment. 1279 for align_index in range(self.num_align): 1280 # Loop over the RDCs. 1281 for j in range(self.num_interatom): 1282 # The back calculated RDC. 1283 if not self.missing_rdc[align_index, j]: 1284 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1285 1286 # Calculate and sum the single alignment chi-squared value (for the RDC). 1287 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1288 1289 # PCS via numerical integration. 1290 if self.pcs_flag: 1291 # Numerical integration of the PCSs. 1292 pcs_numeric_qr_int_pseudo_ellipse(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, theta_x=cone_theta_x, theta_y=cone_theta_y, sigma_max=cone_sigma_max, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 1293 1294 # Calculate and sum the single alignment chi-squared value (for the PCS). 1295 for align_index in range(self.num_align): 1296 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1297 1298 # Return the chi-squared value. 1299 return chi2_sum

1300 1301

1302 - def func_pseudo_ellipse_quad_int(self, params):

1303 """SciPy quadratic integration target function for the pseudo-ellipse model. 1304 1305 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1306 1307 1308 @param params: The vector of parameter values {alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max} where the first 3 are the average position rotation Euler angles, the next 3 are the Euler angles defining the eigenframe, and the last 3 are the pseudo-elliptic cone geometric parameters. 1309 @type params: list of float 1310 @return: The chi-squared or SSE value. 1311 @rtype: float 1312 """ 1313 1314 # Scaling. 1315 if self.scaling_flag: 1316 params = dot(params, self.scaling_matrix) 1317 1318 # Unpack the parameters. 1319 if self.pivot_opt: 1320 pivot = outer(self.spin_ones_struct, params[:3]) 1321 self._translation_vector = params[3:6] 1322 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = params[6:] 1323 else: 1324 pivot = self.pivot 1325 self._translation_vector = params[:3] 1326 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = params[3:] 1327 1328 # Average position rotation. 1329 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1330 1331 # The Kronecker product of the eigenframe rotation. 1332 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1333 1334 # Generate the 2nd degree Frame Order super matrix. 1335 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y, cone_sigma_max) 1336 1337 # Reduce and rotate the tensors. 1338 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1339 1340 # Pre-transpose matrices for faster calculations. 1341 RT_eigen = transpose(self.R_eigen) 1342 RT_ave = transpose(self.R_ave) 1343 1344 # Pre-calculate all the necessary vectors. 1345 if self.pcs_flag: 1346 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1347 1348 # Initial chi-squared (or SSE) value. 1349 chi2_sum = 0.0 1350 1351 # RDCs. 1352 if self.rdc_flag: 1353 # Loop over each alignment. 1354 for align_index in range(self.num_align): 1355 # Loop over the RDCs. 1356 for j in range(self.num_interatom): 1357 # The back calculated RDC. 1358 if not self.missing_rdc[align_index, j]: 1359 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1360 1361 # Calculate and sum the single alignment chi-squared value (for the RDC). 1362 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1363 1364 # PCS via numerical integration. 1365 if self.pcs_flag: 1366 # Loop over each alignment. 1367 for align_index in range(self.num_align): 1368 # Loop over the PCSs. 1369 for j in range(self.num_spins): 1370 # The back calculated PCS. 1371 if not self.missing_pcs[align_index, j]: 1372 # Forwards and reverse rotations. 1373 if self.full_in_ref_frame[align_index]: 1374 r_pivot_atom = self.r_pivot_atom[j] 1375 else: 1376 r_pivot_atom = self.r_pivot_atom_rev[j] 1377 1378 # The numerical integration. 1379 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_pseudo_ellipse(theta_x=cone_theta_x, theta_y=cone_theta_y, sigma_max=cone_sigma_max, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 1380 1381 # Calculate and sum the single alignment chi-squared value (for the PCS). 1382 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1383 1384 # Return the chi-squared value. 1385 return chi2_sum

1386 1387

1388 - def func_pseudo_ellipse_free_rotor_qr_int(self, params):

1389 """Quasi-random Sobol' integration target function for the free-rotor pseudo-ellipse model. 1390 1391 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1392 1393 1394 @param params: The vector of parameter values {alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the first 3 are the average position rotation Euler angles, the next 3 are the Euler angles defining the eigenframe, and the last 2 are the free_rotor pseudo-elliptic cone geometric parameters. 1395 @type params: list of float 1396 @return: The chi-squared or SSE value. 1397 @rtype: float 1398 """ 1399 1400 # Scaling. 1401 if self.scaling_flag: 1402 params = dot(params, self.scaling_matrix) 1403 1404 # Unpack the parameters. 1405 if self.pivot_opt: 1406 pivot = outer(self.spin_ones_struct, params[:3]) 1407 self._translation_vector = params[3:6] 1408 ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1409 else: 1410 pivot = self.pivot 1411 self._translation_vector = params[:3] 1412 ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1413 1414 # Reconstruct the full eigenframe of the motion. 1415 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1416 1417 # The Kronecker product of the eigenframe. 1418 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1419 1420 # Generate the 2nd degree Frame Order super matrix. 1421 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1422 1423 # Reduce and rotate the tensors. 1424 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1425 1426 # Pre-transpose matrices for faster calculations. 1427 RT_eigen = transpose(self.R_eigen) 1428 RT_ave = transpose(self.R_ave) 1429 1430 # Pre-calculate all the necessary vectors. 1431 if self.pcs_flag: 1432 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1433 1434 # Initial chi-squared (or SSE) value. 1435 chi2_sum = 0.0 1436 1437 # RDCs. 1438 if self.rdc_flag: 1439 # Loop over each alignment. 1440 for align_index in range(self.num_align): 1441 # Loop over the RDCs. 1442 for j in range(self.num_interatom): 1443 # The back calculated RDC. 1444 if not self.missing_rdc[align_index, j]: 1445 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1446 1447 # Calculate and sum the single alignment chi-squared value (for the RDC). 1448 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1449 1450 # PCS via numerical integration. 1451 if self.pcs_flag: 1452 # Numerical integration of the PCSs. 1453 pcs_numeric_qr_int_pseudo_ellipse(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, theta_x=cone_theta_x, theta_y=cone_theta_y, sigma_max=pi, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 1454 1455 # Calculate and sum the single alignment chi-squared value (for the PCS). 1456 for align_index in range(self.num_align): 1457 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1458 1459 # Return the chi-squared value. 1460 return chi2_sum

1461 1462

1463 - def func_pseudo_ellipse_free_rotor_quad_int(self, params):

1464 """SciPy quadratic integration target function for the free-rotor pseudo-ellipse model. 1465 1466 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1467 1468 1469 @param params: The vector of parameter values {alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the first 3 are the average position rotation Euler angles, the next 3 are the Euler angles defining the eigenframe, and the last 2 are the free_rotor pseudo-elliptic cone geometric parameters. 1470 @type params: list of float 1471 @return: The chi-squared or SSE value. 1472 @rtype: float 1473 """ 1474 1475 # Scaling. 1476 if self.scaling_flag: 1477 params = dot(params, self.scaling_matrix) 1478 1479 # Unpack the parameters. 1480 if self.pivot_opt: 1481 pivot = outer(self.spin_ones_struct, params[:3]) 1482 self._translation_vector = params[3:6] 1483 ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1484 else: 1485 pivot = self.pivot 1486 self._translation_vector = params[:3] 1487 ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1488 1489 # Average position rotation. 1490 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1491 1492 # The Kronecker product of the eigenframe rotation. 1493 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1494 1495 # Generate the 2nd degree Frame Order super matrix. 1496 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1497 1498 # Reduce and rotate the tensors. 1499 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1500 1501 # Pre-transpose matrices for faster calculations. 1502 RT_eigen = transpose(self.R_eigen) 1503 RT_ave = transpose(self.R_ave) 1504 1505 # Pre-calculate all the necessary vectors. 1506 if self.pcs_flag: 1507 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1508 1509 # Initial chi-squared (or SSE) value. 1510 chi2_sum = 0.0 1511 1512 # RDCs. 1513 if self.rdc_flag: 1514 # Loop over each alignment. 1515 for align_index in range(self.num_align): 1516 # Loop over the RDCs. 1517 for j in range(self.num_interatom): 1518 # The back calculated RDC. 1519 if not self.missing_rdc[align_index, j]: 1520 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1521 1522 # Calculate and sum the single alignment chi-squared value (for the RDC). 1523 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1524 1525 # PCS via numerical integration. 1526 if self.pcs_flag: 1527 # Loop over each alignment. 1528 for align_index in range(self.num_align): 1529 # Loop over the PCSs. 1530 for j in range(self.num_spins): 1531 # The back calculated PCS. 1532 if not self.missing_pcs[align_index, j]: 1533 # Forwards and reverse rotations. 1534 if self.full_in_ref_frame[align_index]: 1535 r_pivot_atom = self.r_pivot_atom[j] 1536 else: 1537 r_pivot_atom = self.r_pivot_atom_rev[j] 1538 1539 # The numerical integration. 1540 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_pseudo_ellipse(theta_x=cone_theta_x, theta_y=cone_theta_y, sigma_max=pi, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 1541 1542 # Calculate and sum the single alignment chi-squared value (for the PCS). 1543 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1544 1545 # Return the chi-squared value. 1546 return chi2_sum

1547 1548

1549 - def func_pseudo_ellipse_torsionless_qr_int(self, params):

1550 """Quasi-random Sobol' integration target function for the torsionless pseudo-ellipse model. 1551 1552 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1553 1554 1555 @param params: The vector of parameter values {alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the first 3 are the average position rotation Euler angles, the next 3 are the Euler angles defining the eigenframe, and the last 2 are the torsionless pseudo-elliptic cone geometric parameters. 1556 @type params: list of float 1557 @return: The chi-squared or SSE value. 1558 @rtype: float 1559 """ 1560 1561 # Scaling. 1562 if self.scaling_flag: 1563 params = dot(params, self.scaling_matrix) 1564 1565 # Unpack the parameters. 1566 if self.pivot_opt: 1567 pivot = outer(self.spin_ones_struct, params[:3]) 1568 self._translation_vector = params[3:6] 1569 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1570 else: 1571 pivot = self.pivot 1572 self._translation_vector = params[:3] 1573 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1574 1575 # Reconstruct the full eigenframe of the motion. 1576 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1577 1578 # The Kronecker product of the eigenframe. 1579 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1580 1581 # Generate the 2nd degree Frame Order super matrix. 1582 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1583 1584 # Reduce and rotate the tensors. 1585 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1586 1587 # Pre-transpose matrices for faster calculations. 1588 RT_eigen = transpose(self.R_eigen) 1589 RT_ave = transpose(self.R_ave) 1590 1591 # Pre-calculate all the necessary vectors. 1592 if self.pcs_flag: 1593 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1594 1595 # Initial chi-squared (or SSE) value. 1596 chi2_sum = 0.0 1597 1598 # RDCs. 1599 if self.rdc_flag: 1600 # Loop over each alignment. 1601 for align_index in range(self.num_align): 1602 # Loop over the RDCs. 1603 for j in range(self.num_interatom): 1604 # The back calculated RDC. 1605 if not self.missing_rdc[align_index, j]: 1606 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1607 1608 # Calculate and sum the single alignment chi-squared value (for the RDC). 1609 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1610 1611 # PCS via numerical integration. 1612 if self.pcs_flag: 1613 # Numerical integration of the PCSs. 1614 pcs_numeric_qr_int_pseudo_ellipse_torsionless(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, theta_x=cone_theta_x, theta_y=cone_theta_y, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 1615 1616 # Calculate and sum the single alignment chi-squared value (for the PCS). 1617 for align_index in range(self.num_align): 1618 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1619 1620 # Return the chi-squared value. 1621 return chi2_sum

1622 1623

1624 - def func_pseudo_ellipse_torsionless_quad_int(self, params):

1625 """SciPy quadratic integration target function for the torsionless pseudo-ellipse model. 1626 1627 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1628 1629 1630 @param params: The vector of parameter values {alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the first 3 are the average position rotation Euler angles, the next 3 are the Euler angles defining the eigenframe, and the last 2 are the torsionless pseudo-elliptic cone geometric parameters. 1631 @type params: list of float 1632 @return: The chi-squared or SSE value. 1633 @rtype: float 1634 """ 1635 1636 # Scaling. 1637 if self.scaling_flag: 1638 params = dot(params, self.scaling_matrix) 1639 1640 # Unpack the parameters. 1641 if self.pivot_opt: 1642 pivot = outer(self.spin_ones_struct, params[:3]) 1643 self._translation_vector = params[3:6] 1644 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1645 else: 1646 pivot = self.pivot 1647 self._translation_vector = params[:3] 1648 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1649 1650 # Average position rotation. 1651 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1652 1653 # The Kronecker product of the eigenframe rotation. 1654 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1655 1656 # Generate the 2nd degree Frame Order super matrix. 1657 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1658 1659 # Reduce and rotate the tensors. 1660 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1661 1662 # Pre-transpose matrices for faster calculations. 1663 RT_eigen = transpose(self.R_eigen) 1664 RT_ave = transpose(self.R_ave) 1665 1666 # Pre-calculate all the necessary vectors. 1667 if self.pcs_flag: 1668 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1669 1670 # Initial chi-squared (or SSE) value. 1671 chi2_sum = 0.0 1672 1673 # RDCs. 1674 if self.rdc_flag: 1675 # Loop over each alignment. 1676 for align_index in range(self.num_align): 1677 # Loop over the RDCs. 1678 for j in range(self.num_interatom): 1679 # The back calculated RDC. 1680 if not self.missing_rdc[align_index, j]: 1681 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1682 1683 # Calculate and sum the single alignment chi-squared value (for the RDC). 1684 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1685 1686 # PCS via numerical integration. 1687 if self.pcs_flag: 1688 # Loop over each alignment. 1689 for align_index in range(self.num_align): 1690 # Loop over the PCSs. 1691 for j in range(self.num_spins): 1692 # The back calculated PCS. 1693 if not self.missing_pcs[align_index, j]: 1694 # Forwards and reverse rotations. 1695 if self.full_in_ref_frame[align_index]: 1696 r_pivot_atom = self.r_pivot_atom[j] 1697 else: 1698 r_pivot_atom = self.r_pivot_atom_rev[j] 1699 1700 # The numerical integration. 1701 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_pseudo_ellipse_torsionless(theta_x=cone_theta_x, theta_y=cone_theta_y, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 1702 1703 # Calculate and sum the single alignment chi-squared value (for the PCS). 1704 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1705 1706 # Return the chi-squared value. 1707 return chi2_sum

1708 1709

1710 - def func_rigid(self, params):

1711 """Target function for rigid model optimisation. 1712 1713 This function optimises the isotropic cone model parameters using the RDC and PCS base data. 1714 1715 1716 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma}. 1717 @type params: list of float 1718 @return: The chi-squared or SSE value. 1719 @rtype: float 1720 """ 1721 1722 # Scaling. 1723 if self.scaling_flag: 1724 params = dot(params, self.scaling_matrix) 1725 1726 # Unpack the parameters. 1727 self._translation_vector = params[:3] 1728 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = params[3:6] 1729 1730 # The average frame rotation matrix (and reduce and rotate the tensors). 1731 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma) 1732 1733 # Pre-transpose matrices for faster calculations. 1734 RT_ave = transpose(self.R_ave) 1735 1736 # Initial chi-squared (or SSE) value. 1737 chi2_sum = 0.0 1738 1739 # RDCs. 1740 if self.rdc_flag: 1741 # Loop over each alignment. 1742 for align_index in range(self.num_align): 1743 # Loop over the RDCs. 1744 for j in range(self.num_interatom): 1745 # The back calculated RDC. 1746 if not self.missing_rdc[align_index, j]: 1747 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1748 1749 # Calculate and sum the single alignment chi-squared value (for the RDC). 1750 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1751 1752 # PCS. 1753 if self.pcs_flag: 1754 # Pre-calculate all the necessary vectors. 1755 self.calc_vectors(pivot=self.pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1756 r_ln_atom = self.r_ln_pivot + self.r_pivot_atom 1757 if min(self.full_in_ref_frame) == 0: 1758 r_ln_atom_rev = self.r_ln_pivot + self.r_pivot_atom_rev 1759 1760 # The vector length (to the inverse 5th power). 1761 length = 1.0 / norm(r_ln_atom, axis=1)**5 1762 if min(self.full_in_ref_frame) == 0: 1763 length_rev = 1.0 / norm(r_ln_atom, axis=1)**5 1764 1765 # Loop over each alignment. 1766 for align_index in range(self.num_align): 1767 # Loop over the PCSs. 1768 for j in range(self.num_spins): 1769 # The back calculated PCS. 1770 if not self.missing_pcs[align_index, j]: 1771 # Forwards and reverse rotations. 1772 if self.full_in_ref_frame[align_index]: 1773 r_ln_atom_i = r_ln_atom[j] 1774 length_i = length[j] 1775 else: 1776 r_ln_atom_i = r_ln_atom_rev[j] 1777 length_i = length_rev[j] 1778 1779 # The PCS calculation. 1780 self.pcs_theta[align_index, j] = pcs_tensor(self.pcs_const[align_index, j] * length_i, r_ln_atom_i, self.A_3D[align_index]) 1781 1782 # Calculate and sum the single alignment chi-squared value (for the PCS). 1783 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1784 1785 # Return the chi-squared value. 1786 return chi2_sum

1787 1788

1789 - def func_rotor_qr_int(self, params):

1790 """Quasi-random Sobol' integration target function for the rotor model. 1791 1792 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Quasi-random, Sobol' sequence based, numerical integration is used for the PCS. 1793 1794 1795 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 1796 @type params: list of float 1797 @return: The chi-squared or SSE value. 1798 @rtype: float 1799 """ 1800 1801 # Scaling. 1802 if self.scaling_flag: 1803 params = dot(params, self.scaling_matrix) 1804 1805 # Unpack the parameters. 1806 if self.pivot_opt: 1807 pivot = outer(self.spin_ones_struct, params[:3]) 1808 self._translation_vector = params[3:6] 1809 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[6:] 1810 else: 1811 pivot = self.pivot 1812 self._translation_vector = params[:3] 1813 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[3:] 1814 1815 # Generate the rotor axis. 1816 self.cone_axis = create_rotor_axis_alpha(alpha=axis_alpha, pivot=pivot[0], point=self.com) 1817 1818 # Pre-calculate the eigenframe rotation matrix. 1819 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1820 1821 # The Kronecker product of the eigenframe rotation. 1822 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1823 1824 # Generate the 2nd degree Frame Order super matrix. 1825 frame_order_2nd = compile_2nd_matrix_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max) 1826 1827 # The average frame rotation matrix (and reduce and rotate the tensors). 1828 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1829 1830 # Pre-transpose matrices for faster calculations. 1831 RT_eigen = transpose(self.R_eigen) 1832 RT_ave = transpose(self.R_ave) 1833 1834 # Pre-calculate all the necessary vectors. 1835 if self.pcs_flag: 1836 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1837 1838 # Initial chi-squared (or SSE) value. 1839 chi2_sum = 0.0 1840 1841 # RDCs. 1842 if self.rdc_flag: 1843 # Loop over each alignment. 1844 for align_index in range(self.num_align): 1845 # Loop over the RDCs. 1846 for j in range(self.num_interatom): 1847 # The back calculated RDC. 1848 if not self.missing_rdc[align_index, j]: 1849 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1850 1851 # Calculate and sum the single alignment chi-squared value (for the RDC). 1852 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1853 1854 # PCS via numerical integration. 1855 if self.pcs_flag: 1856 # Numerical integration of the PCSs. 1857 pcs_numeric_qr_int_rotor(points=sobol_data.sobol_angles, max_points=self.sobol_max_points, sigma_max=sigma_max, c=self.pcs_const, full_in_ref_frame=self.full_in_ref_frame, r_pivot_atom=self.r_pivot_atom, r_pivot_atom_rev=self.r_pivot_atom_rev, r_ln_pivot=self.r_ln_pivot, A=self.A_3D, R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=sobol_data.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs) 1858 1859 # Calculate and sum the single alignment chi-squared value (for the PCS). 1860 for align_index in range(self.num_align): 1861 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1862 1863 # Return the chi-squared value. 1864 return chi2_sum

1865 1866

1867 - def func_rotor_quad_int(self, params):

1868 """SciPy quadratic integration target function for rotor model. 1869 1870 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1871 1872 1873 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 1874 @type params: list of float 1875 @return: The chi-squared or SSE value. 1876 @rtype: float 1877 """ 1878 1879 # Scaling. 1880 if self.scaling_flag: 1881 params = dot(params, self.scaling_matrix) 1882 1883 # Unpack the parameters. 1884 if self.pivot_opt: 1885 pivot = outer(self.spin_ones_struct, params[:3]) 1886 self._translation_vector = params[3:6] 1887 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[6:] 1888 else: 1889 pivot = self.pivot 1890 self._translation_vector = params[:3] 1891 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[3:] 1892 1893 # Generate the rotor axis. 1894 self.cone_axis = create_rotor_axis_alpha(alpha=axis_alpha, pivot=pivot[0], point=self.com) 1895 1896 # Pre-calculate the eigenframe rotation matrix. 1897 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1898 1899 # The Kronecker product of the eigenframe rotation. 1900 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1901 1902 # Generate the 2nd degree Frame Order super matrix. 1903 frame_order_2nd = compile_2nd_matrix_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max) 1904 1905 # The average frame rotation matrix (and reduce and rotate the tensors). 1906 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1907 1908 # Pre-transpose matrices for faster calculations. 1909 RT_eigen = transpose(self.R_eigen) 1910 RT_ave = transpose(self.R_ave) 1911 1912 # Pre-calculate all the necessary vectors. 1913 if self.pcs_flag: 1914 self.calc_vectors(pivot=pivot, R_ave=self.R_ave, RT_ave=RT_ave) 1915 1916 # Initial chi-squared (or SSE) value. 1917 chi2_sum = 0.0 1918 1919 # RDCs. 1920 if self.rdc_flag: 1921 # Loop over each alignment. 1922 for align_index in range(self.num_align): 1923 # Loop over the RDCs. 1924 for j in range(self.num_interatom): 1925 # The back calculated RDC. 1926 if not self.missing_rdc[align_index, j]: 1927 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1928 1929 # Calculate and sum the single alignment chi-squared value (for the RDC). 1930 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1931 1932 # PCS via numerical integration. 1933 if self.pcs_flag: 1934 # Loop over each alignment. 1935 for align_index in range(self.num_align): 1936 # Loop over the PCSs. 1937 for j in range(self.num_spins): 1938 # The back calculated PCS. 1939 if not self.missing_pcs[align_index, j]: 1940 # Forwards and reverse rotations. 1941 if self.full_in_ref_frame[align_index]: 1942 r_pivot_atom = self.r_pivot_atom[j] 1943 else: 1944 r_pivot_atom = self.r_pivot_atom_rev[j] 1945 1946 # The numerical integration. 1947 self.pcs_theta[align_index, j] = pcs_numeric_quad_int_rotor(sigma_max=sigma_max, c=self.pcs_const[align_index, j], r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0], A=self.A_3D[align_index], R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime) 1948 1949 # Calculate and sum the single alignment chi-squared value (for the PCS). 1950 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1951 1952 # Return the chi-squared value. 1953 return chi2_sum

1954 1955

1956 - def calc_vectors(self, pivot=None, pivot2=None, R_ave=None, RT_ave=None):

1957 """Calculate the pivot to atom and lanthanide to pivot vectors for the target functions. 1958 1959 @keyword pivot: The pivot point. 1960 @type pivot: numpy rank-1, 3D array 1961 @keyword pivot2: The 2nd pivot point. 1962 @type pivot2: numpy rank-1, 3D array 1963 @keyword R_ave: The rotation matrix for rotating from the reference frame to the average position. 1964 @type R_ave: numpy rank-2, 3D array 1965 @keyword RT_ave: The transpose of R_ave. 1966 @type RT_ave: numpy rank-2, 3D array 1967 """ 1968 1969 # The lanthanide to pivot vector. 1970 if self.pivot_opt: 1971 subtract(pivot, self.paramag_centre, self.r_ln_pivot) 1972 if pivot2 is not None: 1973 subtract(pivot2, self.paramag_centre, self.r_ln_pivot) 1974 1975 # Calculate the average position pivot point to atomic positions vectors once. 1976 vect = self.atomic_pos - self.ave_pos_pivot 1977 1978 # Rotate then translate the atomic positions, then calculate the pivot to atom vector. 1979 self.r_pivot_atom[:] = dot(vect, RT_ave) 1980 add(self.r_pivot_atom, self.ave_pos_pivot, self.r_pivot_atom) 1981 add(self.r_pivot_atom, self._translation_vector, self.r_pivot_atom) 1982 subtract(self.r_pivot_atom, pivot, self.r_pivot_atom) 1983 1984 # And the reverse vectors. 1985 if min(self.full_in_ref_frame) == 0: 1986 self.r_pivot_atom_rev[:] = dot(vect, R_ave) 1987 add(self.r_pivot_atom_rev, self.ave_pos_pivot, self.r_pivot_atom_rev) 1988 add(self.r_pivot_atom_rev, self._translation_vector, self.r_pivot_atom_rev) 1989 subtract(self.r_pivot_atom_rev, pivot, self.r_pivot_atom_rev) 1990 1991 # Calculate the inter-pivot vector for the double motion models. 1992 if pivot2 is not None: 1993 self.r_inter_pivot = pivot - pivot2

1994 1995

1996 - def create_sobol_data(self, dims=None):

1997 """Create the Sobol' quasi-random data for numerical integration. 1998 1999 This uses the external sobol_lib module to create the data. The algorithm is that modified by Antonov and Saleev. 2000 2001 2002 @keyword dims: The list of parameters. 2003 @type dims: list of str 2004 """ 2005 2006 # Quadratic integration active, so nothing to do here! 2007 if self.quad_int: 2008 return 2009 2010 # The number of dimensions. 2011 m = len(dims) 2012 2013 # The total number of points. 2014 total_num = int(self.sobol_max_points * self.sobol_oversample * 10**m) 2015 2016 # Reuse pre-created data if available. 2017 if total_num == sobol_data.total_num and self.model == sobol_data.model: 2018 return 2019 2020 # Printout (useful to see how long this takes!). 2021 print("Generating the torsion-tilt angle sampling via the Sobol' sequence for numerical PCS integration.") 2022 2023 # Initialise. 2024 sobol_data.model = self.model 2025 sobol_data.total_num = total_num 2026 sobol_data.sobol_angles = zeros((m, total_num), float32) 2027 sobol_data.Ri_prime = zeros((total_num, 3, 3), float32) 2028 sobol_data.Ri2_prime = zeros((total_num, 3, 3), float32) 2029 2030 # The Sobol' points. 2031 points = i4_sobol_generate(m, total_num, 1000) 2032 2033 # Loop over the points. 2034 for i in range(total_num): 2035 # Loop over the dimensions, converting the points to angles. 2036 theta = None 2037 phi = None 2038 sigma = None 2039 for j in range(m): 2040 # The tilt angle - the angle of rotation about the x-y plane rotation axis. 2041 if dims[j] in ['theta']: 2042 theta = acos(2.0*points[j, i] - 1.0) 2043 sobol_data.sobol_angles[j, i] = theta 2044 2045 # The angle defining the x-y plane rotation axis. 2046 if dims[j] in ['phi']: 2047 phi = 2.0 * pi * points[j, i] 2048 sobol_data.sobol_angles[j, i] = phi 2049 2050 # The 1st torsion angle - the angle of rotation about the z' axis (or y' for the double motion models). 2051 if dims[j] in ['sigma']: 2052 sigma = 2.0 * pi * (points[j, i] - 0.5) 2053 sobol_data.sobol_angles[j, i] = sigma 2054 2055 # The 2nd torsion angle - the angle of rotation about the x' axis. 2056 if dims[j] in ['sigma2']: 2057 sigma2 = 2.0 * pi * (points[j, i] - 0.5) 2058 sobol_data.sobol_angles[j, i] = sigma2 2059 2060 # Pre-calculate the rotation matrices for the double motion models. 2061 if 'sigma2' in dims: 2062 # The 1st rotation about the y-axis. 2063 c_sigma = cos(sigma) 2064 s_sigma = sin(sigma) 2065 sobol_data.Ri_prime[i, 0, 0] = c_sigma 2066 sobol_data.Ri_prime[i, 0, 2] = s_sigma 2067 sobol_data.Ri_prime[i, 1, 1] = 1.0 2068 sobol_data.Ri_prime[i, 2, 0] = -s_sigma 2069 sobol_data.Ri_prime[i, 2, 2] = c_sigma 2070 2071 # The 2nd rotation about the x-axis. 2072 c_sigma2 = cos(sigma2) 2073 s_sigma2 = sin(sigma2) 2074 sobol_data.Ri2_prime[i, 0, 0] = 1.0 2075 sobol_data.Ri2_prime[i, 1, 1] = c_sigma2 2076 sobol_data.Ri2_prime[i, 1, 2] = -s_sigma2 2077 sobol_data.Ri2_prime[i, 2, 1] = s_sigma2 2078 sobol_data.Ri2_prime[i, 2, 2] = c_sigma2 2079 2080 # Pre-calculate the rotation matrix for the full tilt-torsion. 2081 elif theta != None and phi != None and sigma != None: 2082 tilt_torsion_to_R(phi, theta, sigma, sobol_data.Ri_prime[i]) 2083 2084 # Pre-calculate the rotation matrix for the torsionless models. 2085 elif sigma == None: 2086 c_theta = cos(theta) 2087 s_theta = sin(theta) 2088 c_phi = cos(phi) 2089 s_phi = sin(phi) 2090 c_phi_c_theta = c_phi * c_theta 2091 s_phi_c_theta = s_phi * c_theta 2092 sobol_data.Ri_prime[i, 0, 0] = c_phi_c_theta*c_phi + s_phi**2 2093 sobol_data.Ri_prime[i, 0, 1] = c_phi_c_theta*s_phi - c_phi*s_phi 2094 sobol_data.Ri_prime[i, 0, 2] = c_phi*s_theta 2095 sobol_data.Ri_prime[i, 1, 0] = s_phi_c_theta*c_phi - c_phi*s_phi 2096 sobol_data.Ri_prime[i, 1, 1] = s_phi_c_theta*s_phi + c_phi**2 2097 sobol_data.Ri_prime[i, 1, 2] = s_phi*s_theta 2098 sobol_data.Ri_prime[i, 2, 0] = -s_theta*c_phi 2099 sobol_data.Ri_prime[i, 2, 1] = -s_theta*s_phi 2100 sobol_data.Ri_prime[i, 2, 2] = c_theta 2101 2102 # Pre-calculate the rotation matrix for the rotor models. 2103 else: 2104 c_sigma = cos(sigma) 2105 s_sigma = sin(sigma) 2106 sobol_data.Ri_prime[i, 0, 0] = c_sigma 2107 sobol_data.Ri_prime[i, 0, 1] = -s_sigma 2108 sobol_data.Ri_prime[i, 1, 0] = s_sigma 2109 sobol_data.Ri_prime[i, 1, 1] = c_sigma 2110 sobol_data.Ri_prime[i, 2, 2] = 1.0 2111 2112 # Printout (useful to see how long this takes!). 2113 print(" Oversampled to %s points." % total_num)

2114 2115

2116 - def reduce_and_rot(self, ave_pos_alpha=None, ave_pos_beta=None, ave_pos_gamma=None, daeg=None):

2117 """Reduce and rotate the alignments tensors using the frame order matrix and Euler angles. 2118 2119 @keyword ave_pos_alpha: The alpha Euler angle describing the average domain position, the tensor rotation. 2120 @type ave_pos_alpha: float 2121 @keyword ave_pos_beta: The beta Euler angle describing the average domain position, the tensor rotation. 2122 @type ave_pos_beta: float 2123 @keyword ave_pos_gamma: The gamma Euler angle describing the average domain position, the tensor rotation. 2124 @type ave_pos_gamma: float 2125 @keyword daeg: The 2nd degree frame order matrix. 2126 @type daeg: rank-2, 9D array 2127 """ 2128 2129 # Alignment tensor rotation. 2130 euler_to_R_zyz(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, self.R_ave) 2131 2132 # Back calculate the rotated tensors. 2133 for align_index in range(self.num_tensors): 2134 # Tensor indices. 2135 index1 = align_index*5 2136 index2 = align_index*5+5 2137 2138 # Reduction. 2139 if daeg is not None: 2140 # Reduce the tensor. 2141 reduce_alignment_tensor(daeg, self.full_tensors[index1:index2], self.A_5D_bc[index1:index2]) 2142 2143 # Convert the reduced tensor to 3D, rank-2 form. 2144 to_tensor(self.tensor_3D, self.A_5D_bc[index1:index2]) 2145 2146 # No reduction: 2147 else: 2148 # Convert the original tensor to 3D, rank-2 form. 2149 to_tensor(self.tensor_3D, self.full_tensors[index1:index2]) 2150 2151 # Rotate the tensor (normal R.X.RT rotation). 2152 if self.full_in_ref_frame[align_index]: 2153 self.A_3D_bc[align_index] = dot(transpose(self.R_ave), dot(self.tensor_3D, self.R_ave)) 2154 2155 # Rotate the tensor (inverse RT.X.R rotation). 2156 else: 2157 self.A_3D_bc[align_index] = dot(self.R_ave, dot(self.tensor_3D, transpose(self.R_ave))) 2158 2159 # Convert the tensor back to 5D, rank-1 form, as the back-calculated reduced tensor. 2160 to_5D(self.A_5D_bc[index1:index2], self.A_3D_bc[align_index])

Source Code for Module target_functions.frame_order