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, num_int_pts=500, com=None, ave_pos_pivot=zeros(3), ave_pos_piv_sync=True, translation_opt=False, pivot=None, pivot2=None, pivot_opt=False, quad_int=True):

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 num_int_pts: The number of points to use for the numerical integration technique. 98 @type num_int_pts: int 99 @keyword com: The centre of mass of the system. This is used for defining the rotor model systems. 100 @type com: numpy 3D rank-1 array 101 @keyword ave_pos_pivot: The pivot point to rotate all atoms about to the average domain position. For example this can be the centre of mass of the moving domain. 102 @type ave_pos_pivot: numpy 3D rank-1 array 103 @keyword ave_pos_piv_sync: A flag which if True will cause pivot point to rotate to the average domain position to be synchronised with the motional pivot. This will cause ave_pos_pivot argument to be ignored. 104 @type ave_pos_piv_sync: bool 105 @keyword translation_opt: A flag which if True will allow the pivot point of the motion to be optimised. 106 @type translation_opt: bool 107 @keyword pivot: The pivot point for the ball-and-socket joint motion. This is needed if PCS or PRE values are used. 108 @type pivot: numpy rank-1, 3D array or None 109 @keyword pivot2: The second pivot point for the motion. This is needed if PCS or PRE values are used and if a double-motional model is to be optimised. 110 @type pivot2: numpy rank-1, 3D array or None 111 @keyword pivot_opt: A flag which if True will allow the pivot point of the motion to be optimised. 112 @type pivot_opt: bool 113 @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. 114 @type quad_int: bool 115 """ 116 117 # Model test. 118 if not model: 119 raise RelaxError("The type of Frame Order model must be specified.") 120 121 # Store the initial parameter (as a copy). 122 self.params = deepcopy(init_params) 123 124 # Store the agrs. 125 self.model = model 126 self.full_tensors = full_tensors 127 self.full_in_ref_frame = full_in_ref_frame 128 self.rdc = rdcs 129 self.rdc_weights = rdc_weights 130 self.rdc_vect = rdc_vect 131 self.dip_const = dip_const 132 self.pcs = pcs 133 self.pcs_weights = pcs_weights 134 self.atomic_pos = atomic_pos 135 self.temp = temp 136 self.frq = frq 137 self.paramag_centre = paramag_centre 138 self.total_num_params = len(init_params) 139 self.num_int_pts = num_int_pts 140 self.com = com 141 self.ave_pos_pivot = ave_pos_pivot 142 self.ave_pos_piv_sync = ave_pos_piv_sync 143 self.translation_opt = translation_opt 144 self._param_pivot = pivot 145 self._param_pivot2 = pivot2 146 self.pivot_opt = pivot_opt 147 148 # Tensor setup. 149 self._init_tensors() 150 151 # Scaling initialisation. 152 self.scaling_matrix = scaling_matrix 153 if self.scaling_matrix != None: 154 self.scaling_flag = True 155 else: 156 self.scaling_flag = False 157 158 # The total number of alignments. 159 self.num_align = 0 160 if rdcs != None: 161 self.num_align = len(rdcs) 162 elif pcs != None: 163 self.num_align = len(pcs) 164 165 # Set the RDC and PCS flags (indicating the presence of data). 166 self.rdc_flag = [True] * self.num_align 167 self.pcs_flag = [True] * self.num_align 168 for align_index in range(self.num_align): 169 if rdcs == None or len(rdcs[align_index]) == 0: 170 self.rdc_flag[align_index] = False 171 if pcs == None or len(pcs[align_index]) == 0: 172 self.pcs_flag[align_index] = False 173 self.rdc_flag_sum = sum(self.rdc_flag) 174 self.pcs_flag_sum = sum(self.pcs_flag) 175 176 # Default translation vector (if not optimised). 177 self._translation_vector = zeros(3, float64) 178 179 # Some checks. 180 if self.rdc_flag_sum and (rdc_vect == None or not len(rdc_vect)): 181 raise RelaxError("The rdc_vect argument " + repr(rdc_vect) + " must be supplied.") 182 if self.pcs_flag_sum and (atomic_pos == None or not len(atomic_pos)): 183 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.") 184 185 # The total number of spins. 186 self.num_spins = 0 187 if self.pcs_flag_sum: 188 self.num_spins = len(pcs[0]) 189 190 # The total number of interatomic connections. 191 self.num_interatom = 0 192 if self.rdc_flag_sum: 193 self.num_interatom = len(rdcs[0]) 194 195 # Set up the alignment data. 196 for align_index in range(self.num_align): 197 to_tensor(self.A_3D[align_index], self.full_tensors[5*align_index:5*align_index+5]) 198 199 # PCS errors. 200 if self.pcs_flag_sum: 201 err = False 202 for i in range(len(pcs_errors)): 203 for j in range(len(pcs_errors[i])): 204 if not isNaN(pcs_errors[i, j]): 205 err = True 206 if err: 207 self.pcs_error = pcs_errors 208 else: 209 # Missing errors (default to 0.1 ppm errors). 210 self.pcs_error = 0.1 * 1e-6 * ones((self.num_align, self.num_spins), float64) 211 212 # RDC errors. 213 if self.rdc_flag_sum: 214 err = False 215 for i in range(len(rdc_errors)): 216 for j in range(len(rdc_errors[i])): 217 if not isNaN(rdc_errors[i, j]): 218 err = True 219 if err: 220 self.rdc_error = rdc_errors 221 else: 222 # Missing errors (default to 1 Hz errors). 223 self.rdc_error = ones((self.num_align, self.num_interatom), float64) 224 225 # Missing data matrices (RDC). 226 if self.rdc_flag_sum: 227 self.missing_rdc = zeros((self.num_align, self.num_interatom), uint8) 228 229 # Missing data matrices (PCS). 230 if self.pcs_flag_sum: 231 self.missing_pcs = zeros((self.num_align, self.num_spins), uint8) 232 233 # Clean up problematic data and put the weights into the errors.. 234 if self.rdc_flag_sum or self.pcs_flag_sum: 235 for align_index in range(self.num_align): 236 # Loop over the RDCs. 237 if self.rdc_flag_sum: 238 for j in range(self.num_interatom): 239 if isNaN(self.rdc[align_index, j]): 240 # Set the flag. 241 self.missing_rdc[align_index, j] = 1 242 243 # Change the NaN to zero. 244 self.rdc[align_index, j] = 0.0 245 246 # Change the error to one (to avoid zero division). 247 self.rdc_error[align_index, j] = 1.0 248 249 # Change the weight to one. 250 rdc_weights[align_index, j] = 1.0 251 252 # The RDC weights. 253 if self.rdc_flag_sum: 254 self.rdc_error[align_index, j] = self.rdc_error[align_index, j] / sqrt(rdc_weights[align_index, j]) 255 256 # Loop over the PCSs. 257 if self.pcs_flag_sum: 258 for j in range(self.num_spins): 259 if isNaN(self.pcs[align_index, j]): 260 # Set the flag. 261 self.missing_pcs[align_index, j] = 1 262 263 # Change the NaN to zero. 264 self.pcs[align_index, j] = 0.0 265 266 # Change the error to one (to avoid zero division). 267 self.pcs_error[align_index, j] = 1.0 268 269 # Change the weight to one. 270 pcs_weights[align_index, j] = 1.0 271 272 # The PCS weights. 273 if self.pcs_flag_sum: 274 self.pcs_error[align_index, j] = self.pcs_error[align_index, j] / sqrt(pcs_weights[align_index, j]) 275 276 # The paramagnetic centre vectors and distances. 277 if self.pcs_flag_sum: 278 # Initialise the data structures. 279 self.paramag_unit_vect = zeros(atomic_pos.shape, float64) 280 self.paramag_dist = zeros(self.num_spins, float64) 281 self.pcs_const = zeros(self.num_align, float64) 282 self.r_pivot_atom = zeros((3, self.num_spins), float32) 283 self.r_pivot_atom_rev = zeros((3, self.num_spins), float32) 284 self.r_ln_pivot = zeros((3, self.num_spins), float32) 285 for j in range(self.num_spins): 286 self.r_ln_pivot[:, j] = pivot - self.paramag_centre 287 if self.paramag_centre == None: 288 self.paramag_centre = zeros(3, float32) 289 290 # Set up the paramagnetic constant (without the interatomic distance and in Angstrom units). 291 for align_index in range(self.num_align): 292 self.pcs_const[align_index] = pcs_constant(self.temp[align_index], self.frq[align_index], 1.0) * 1e30 293 294 # PCS function, gradient, and Hessian matrices. 295 if self.pcs_flag_sum: 296 self.pcs_theta = zeros((self.num_align, self.num_spins), float64) 297 self.pcs_theta_err = zeros((self.num_align, self.num_spins), float64) 298 self.dpcs_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 299 self.d2pcs_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 300 301 # RDC function, gradient, and Hessian matrices. 302 if self.rdc_flag_sum: 303 self.rdc_theta = zeros((self.num_align, self.num_interatom), float64) 304 self.drdc_theta = zeros((self.total_num_params, self.num_align, self.num_interatom), float64) 305 self.d2rdc_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_interatom), float64) 306 307 # The quasi-random integration. 308 if not quad_int and self.pcs_flag_sum and model not in ['rigid']: 309 # The Sobol' sequence data and target function aliases (quasi-random integration). 310 if model == 'pseudo-ellipse': 311 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 312 self.func = self.func_pseudo_ellipse_qrint 313 elif model == 'pseudo-ellipse, torsionless': 314 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi']) 315 self.func = self.func_pseudo_ellipse_torsionless_qrint 316 elif model == 'pseudo-ellipse, free rotor': 317 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 318 self.func = self.func_pseudo_ellipse_free_rotor_qrint 319 elif model == 'iso cone': 320 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 321 self.func = self.func_iso_cone_qrint 322 elif model == 'iso cone, torsionless': 323 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi']) 324 self.func = self.func_iso_cone_torsionless_qrint 325 elif model == 'iso cone, free rotor': 326 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 327 self.func = self.func_iso_cone_free_rotor_qrint 328 elif model == 'line': 329 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'sigma']) 330 self.func = self.func_line_qrint 331 elif model == 'line, torsionless': 332 self.create_sobol_data(n=self.num_int_pts, dims=['theta']) 333 self.func = self.func_line_torsionless_qrint 334 elif model == 'line, free rotor': 335 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'sigma']) 336 self.func = self.func_line_free_rotor_qrint 337 elif model == 'rotor': 338 self.create_sobol_data(n=self.num_int_pts, dims=['sigma']) 339 self.func = self.func_rotor_qrint 340 elif model == 'rigid': 341 self.func = self.func_rigid 342 elif model == 'free rotor': 343 self.create_sobol_data(n=self.num_int_pts, dims=['sigma']) 344 self.func = self.func_free_rotor_qrint 345 elif model == 'double rotor': 346 self.create_sobol_data(n=self.num_int_pts, dims=['sigma', 'sigma2']) 347 self.func = self.func_double_rotor 348 349 # The target function aliases (Scipy numerical integration). 350 else: 351 if model == 'pseudo-ellipse': 352 self.func = self.func_pseudo_ellipse 353 elif model == 'pseudo-ellipse, torsionless': 354 self.func = self.func_pseudo_ellipse_torsionless 355 elif model == 'pseudo-ellipse, free rotor': 356 self.func = self.func_pseudo_ellipse_free_rotor 357 elif model == 'iso cone': 358 self.func = self.func_iso_cone 359 elif model == 'iso cone, torsionless': 360 self.func = self.func_iso_cone_torsionless 361 elif model == 'iso cone, free rotor': 362 self.func = self.func_iso_cone_free_rotor 363 elif model == 'line': 364 self.func = self.func_line 365 elif model == 'line, torsionless': 366 self.func = self.func_line_torsionless 367 elif model == 'line, free rotor': 368 self.func = self.func_line_free_rotor 369 elif model == 'rotor': 370 self.func = self.func_rotor 371 elif model == 'rigid': 372 self.func = self.func_rigid 373 elif model == 'free rotor': 374 self.func = self.func_free_rotor 375 elif model == 'double rotor': 376 self.func = self.func_double_rotor

377 378

379 - def _init_tensors(self):

380 """Set up isotropic cone optimisation against the alignment tensor data.""" 381 382 # Some checks. 383 if self.full_tensors == None or not len(self.full_tensors): 384 raise RelaxError("The full_tensors argument " + repr(self.full_tensors) + " must be supplied.") 385 if self.full_in_ref_frame == None or not len(self.full_in_ref_frame): 386 raise RelaxError("The full_in_ref_frame argument " + repr(self.full_in_ref_frame) + " must be supplied.") 387 388 # Tensor set up. 389 self.num_tensors = int(len(self.full_tensors) / 5) 390 self.A_3D = zeros((self.num_tensors, 3, 3), float64) 391 self.A_3D_bc = zeros((self.num_tensors, 3, 3), float64) 392 self.A_5D_bc = zeros(self.num_tensors*5, float64) 393 394 # The rotation to the Frame Order eigenframe. 395 self.R_eigen = zeros((3, 3), float64) 396 self.R_eigen_2 = zeros((3, 3), float64) 397 self.R_ave = zeros((3, 3), float64) 398 self.Ri_prime = zeros((3, 3), float64) 399 self.tensor_3D = zeros((3, 3), float64) 400 401 # The cone axis storage and molecular frame z-axis. 402 self.cone_axis = zeros(3, float64) 403 self.z_axis = array([0, 0, 1], float64) 404 405 # The rotor axes. 406 self.rotor_axis = zeros(3, float64) 407 self.rotor_axis_2 = zeros(3, float64) 408 409 # Initialise the Frame Order matrices. 410 self.frame_order_2nd = zeros((9, 9), float64) 411 412 # A rotation matrix for general use. 413 self.R = zeros((3, 3), float64)

414 415

416 - def func_double_rotor(self, params):

417 """Target function for the optimisation of the double rotor frame order model. 418 419 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. 420 421 422 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 423 @type params: list of float 424 @return: The chi-squared or SSE value. 425 @rtype: float 426 """ 427 428 # Scaling. 429 if self.scaling_flag: 430 params = dot(params, self.scaling_matrix) 431 432 # Unpack the parameters. 433 if self.translation_opt and self.pivot_opt: 434 self._param_pivot = params[:3] 435 self._param_pivot2 = params[3:6] 436 self._translation_vector = params[6:9] 437 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, axis_theta_2, axis_phi_2, sigma_max, sigma_max_2 = params[9:] 438 elif self.translation_opt: 439 self._translation_vector = params[:3] 440 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, axis_theta_2, axis_phi_2, sigma_max, sigma_max_2 = params[3:] 441 elif self.pivot_opt: 442 self._param_pivot = params[:3] 443 self._param_pivot2 = params[3:6] 444 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, axis_theta_2, axis_phi_2, sigma_max, sigma_max_2 = params[6:] 445 else: 446 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, axis_theta_2, axis_phi_2, sigma_max, sigma_max_2 = params 447 448 # Generate both rotation axes from the spherical angles. 449 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.rotor_axis) 450 spherical_to_cartesian([1.0, axis_theta_2, axis_phi_2], self.rotor_axis_2) 451 452 # Pre-calculate the eigenframe rotation matrix. 453 two_vect_to_R(self.z_axis, self.rotor_axis, self.R_eigen) 454 two_vect_to_R(self.z_axis, self.rotor_axis_2, self.R_eigen_2) 455 456 # Combine the rotations. 457 R_eigen_full = dot(self.R_eigen_2, self.R_eigen) 458 459 # The Kronecker product of the eigenframe rotation. 460 Rx2_eigen = kron_prod(R_eigen_full, R_eigen_full) 461 462 # Generate the 2nd degree Frame Order super matrix. 463 frame_order_2nd = compile_2nd_matrix_double_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max, sigma_max_2) 464 465 # The average frame rotation matrix (and reduce and rotate the tensors). 466 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 467 468 # Pre-transpose matrices for faster calculations. 469 RT_eigen = transpose(R_eigen_full) 470 RT_ave = transpose(self.R_ave) 471 472 # Pre-calculate all the necessary vectors. 473 if self.pcs_flag: 474 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 475 476 # Initial chi-squared (or SSE) value. 477 chi2_sum = 0.0 478 479 # Loop over each alignment. 480 for align_index in range(self.num_align): 481 # RDCs. 482 if self.rdc_flag[align_index]: 483 # Loop over the RDCs. 484 for j in range(self.num_interatom): 485 # The back calculated RDC. 486 if not self.missing_rdc[align_index, j]: 487 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 488 489 # Calculate and sum the single alignment chi-squared value (for the RDC). 490 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 491 492 # PCS via numerical integration. 493 if self.pcs_flag_sum: 494 # Numerical integration of the PCSs. 495 pcs_numeric_int_double_rotor(points=self.sobol_angles, 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, A=self.A_3D, R_eigen=R_eigen_full, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 496 497 # Calculate and sum the single alignment chi-squared value (for the PCS). 498 for align_index in range(self.num_align): 499 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 500 501 # Return the chi-squared value. 502 return chi2_sum

503 504

505 - def func_free_rotor(self, params):

506 """Target function for free rotor model optimisation. 507 508 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 509 510 511 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi}. 512 @type params: list of float 513 @return: The chi-squared or SSE value. 514 @rtype: float 515 """ 516 517 # Scaling. 518 if self.scaling_flag: 519 params = dot(params, self.scaling_matrix) 520 521 # Unpack the parameters. 522 if self.translation_opt and self.pivot_opt: 523 self._param_pivot = params[:3] 524 self._translation_vector = params[3:6] 525 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[6:] 526 elif self.translation_opt: 527 self._translation_vector = params[:3] 528 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 529 elif self.pivot_opt: 530 self._param_pivot = params[:3] 531 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 532 else: 533 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params 534 535 # Generate the cone axis from the spherical angles. 536 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 537 538 # Pre-calculate the eigenframe rotation matrix. 539 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 540 541 # The Kronecker product of the eigenframe rotation. 542 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 543 544 # Generate the 2nd degree Frame Order super matrix. 545 frame_order_2nd = compile_2nd_matrix_free_rotor(self.frame_order_2nd, Rx2_eigen) 546 547 # Reduce and rotate the tensors. 548 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 549 550 # Pre-transpose matrices for faster calculations. 551 RT_eigen = transpose(self.R_eigen) 552 RT_ave = transpose(self.R_ave) 553 554 # Pre-calculate all the necessary vectors. 555 if self.pcs_flag: 556 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 557 558 # Initial chi-squared (or SSE) value. 559 chi2_sum = 0.0 560 561 # Loop over each alignment. 562 for align_index in range(self.num_align): 563 # RDCs. 564 if self.rdc_flag[align_index]: 565 # Loop over the RDCs. 566 for j in range(self.num_interatom): 567 # The back calculated RDC. 568 if not self.missing_rdc[align_index, j]: 569 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 570 571 # Calculate and sum the single alignment chi-squared value (for the RDC). 572 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 573 574 # PCS. 575 if self.pcs_flag[align_index]: 576 # Loop over the PCSs. 577 for j in range(self.num_spins): 578 # The back calculated PCS. 579 if not self.missing_pcs[align_index, j]: 580 # Forwards and reverse rotations. 581 if self.full_in_ref_frame[align_index]: 582 r_pivot_atom = self.r_pivot_atom[:, j] 583 else: 584 r_pivot_atom = self.r_pivot_atom_rev[:, j] 585 586 # The numerical integration. 587 self.pcs_theta[align_index, j] = pcs_numeric_int_rotor(sigma_max=pi, c=self.pcs_const[align_index], 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) 588 589 # Calculate and sum the single alignment chi-squared value (for the PCS). 590 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 591 592 # Return the chi-squared value. 593 return chi2_sum

594 595

596 - def func_free_rotor_qrint(self, params):

597 """Target function for free rotor model optimisation. 598 599 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 600 601 602 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi}. 603 @type params: list of float 604 @return: The chi-squared or SSE value. 605 @rtype: float 606 """ 607 608 # Scaling. 609 if self.scaling_flag: 610 params = dot(params, self.scaling_matrix) 611 612 # Unpack the parameters. 613 if self.translation_opt and self.pivot_opt: 614 self._param_pivot = params[:3] 615 self._translation_vector = params[3:6] 616 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[6:] 617 elif self.translation_opt: 618 self._translation_vector = params[:3] 619 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 620 elif self.pivot_opt: 621 self._param_pivot = params[:3] 622 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 623 else: 624 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params 625 626 # Generate the cone axis from the spherical angles. 627 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 628 629 # Pre-calculate the eigenframe rotation matrix. 630 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 631 632 # The Kronecker product of the eigenframe rotation. 633 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 634 635 # Generate the 2nd degree Frame Order super matrix. 636 frame_order_2nd = compile_2nd_matrix_free_rotor(self.frame_order_2nd, Rx2_eigen) 637 638 # Reduce and rotate the tensors. 639 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 640 641 # Pre-transpose matrices for faster calculations. 642 RT_eigen = transpose(self.R_eigen) 643 RT_ave = transpose(self.R_ave) 644 645 # Pre-calculate all the necessary vectors. 646 if self.pcs_flag: 647 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 648 649 # Initial chi-squared (or SSE) value. 650 chi2_sum = 0.0 651 652 # Loop over each alignment. 653 for align_index in range(self.num_align): 654 # RDCs. 655 if self.rdc_flag[align_index]: 656 # Loop over the RDCs. 657 for j in range(self.num_interatom): 658 # The back calculated RDC. 659 if not self.missing_rdc[align_index, j]: 660 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 661 662 # Calculate and sum the single alignment chi-squared value (for the RDC). 663 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 664 665 # PCS via numerical integration. 666 if self.pcs_flag_sum: 667 # Numerical integration of the PCSs. 668 pcs_numeric_int_rotor_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 669 670 # Calculate and sum the single alignment chi-squared value (for the PCS). 671 for align_index in range(self.num_align): 672 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 673 674 675 # Return the chi-squared value. 676 return chi2_sum

677 678

679 - def func_iso_cone(self, params):

680 """Target function for isotropic cone model optimisation. 681 682 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 683 684 685 @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. 686 @type params: list of float 687 @return: The chi-squared or SSE value. 688 @rtype: float 689 """ 690 691 # Scaling. 692 if self.scaling_flag: 693 params = dot(params, self.scaling_matrix) 694 695 # Unpack the parameters. 696 if self.translation_opt and self.pivot_opt: 697 self._param_pivot = params[:3] 698 self._translation_vector = params[3:6] 699 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[6:] 700 elif self.translation_opt: 701 self._translation_vector = params[:3] 702 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 703 elif self.pivot_opt: 704 self._param_pivot = params[:3] 705 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 706 else: 707 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params 708 709 # Generate the cone axis from the spherical angles. 710 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 711 712 # Pre-calculate the eigenframe rotation matrix. 713 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 714 715 # The Kronecker product of the eigenframe rotation. 716 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 717 718 # Generate the 2nd degree Frame Order super matrix. 719 frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd, Rx2_eigen, cone_theta, sigma_max) 720 721 # Reduce and rotate the tensors. 722 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 723 724 # Pre-transpose matrices for faster calculations. 725 RT_eigen = transpose(self.R_eigen) 726 RT_ave = transpose(self.R_ave) 727 728 # Pre-calculate all the necessary vectors. 729 if self.pcs_flag: 730 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 731 732 # Initial chi-squared (or SSE) value. 733 chi2_sum = 0.0 734 735 # Loop over each alignment. 736 for align_index in range(self.num_align): 737 # RDCs. 738 if self.rdc_flag[align_index]: 739 # Loop over the RDCs. 740 for j in range(self.num_interatom): 741 # The back calculated RDC. 742 if not self.missing_rdc[align_index, j]: 743 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 744 745 # Calculate and sum the single alignment chi-squared value (for the RDC). 746 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 747 748 # PCS. 749 if self.pcs_flag[align_index]: 750 # Loop over the PCSs. 751 for j in range(self.num_spins): 752 # The back calculated PCS. 753 if not self.missing_pcs[align_index, j]: 754 # Forwards and reverse rotations. 755 if self.full_in_ref_frame[align_index]: 756 r_pivot_atom = self.r_pivot_atom[:, j] 757 else: 758 r_pivot_atom = self.r_pivot_atom_rev[:, j] 759 760 # The numerical integration. 761 self.pcs_theta[align_index, j] = pcs_numeric_int_iso_cone(theta_max=cone_theta, sigma_max=sigma_max, c=self.pcs_const[align_index], 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) 762 763 # Calculate and sum the single alignment chi-squared value (for the PCS). 764 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 765 766 # Return the chi-squared value. 767 return chi2_sum

768 769

770 - def func_iso_cone_qrint(self, params):

771 """Target function for isotropic cone model optimisation. 772 773 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 774 775 776 @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. 777 @type params: list of float 778 @return: The chi-squared or SSE value. 779 @rtype: float 780 """ 781 782 # Scaling. 783 if self.scaling_flag: 784 params = dot(params, self.scaling_matrix) 785 786 # Unpack the parameters. 787 if self.translation_opt and self.pivot_opt: 788 self._param_pivot = params[:3] 789 self._translation_vector = params[3:6] 790 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[6:] 791 elif self.translation_opt: 792 self._translation_vector = params[:3] 793 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 794 elif self.pivot_opt: 795 self._param_pivot = params[:3] 796 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 797 else: 798 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params 799 800 # Generate the cone axis from the spherical angles. 801 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 802 803 # Pre-calculate the eigenframe rotation matrix. 804 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 805 806 # The Kronecker product of the eigenframe rotation. 807 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 808 809 # Generate the 2nd degree Frame Order super matrix. 810 frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd, Rx2_eigen, cone_theta, sigma_max) 811 812 # Reduce and rotate the tensors. 813 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 814 815 # Pre-transpose matrices for faster calculations. 816 RT_eigen = transpose(self.R_eigen) 817 RT_ave = transpose(self.R_ave) 818 819 # Pre-calculate all the necessary vectors. 820 if self.pcs_flag: 821 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 822 823 # Initial chi-squared (or SSE) value. 824 chi2_sum = 0.0 825 826 # Loop over each alignment. 827 for align_index in range(self.num_align): 828 # RDCs. 829 if self.rdc_flag[align_index]: 830 # Loop over the RDCs. 831 for j in range(self.num_interatom): 832 # The back calculated RDC. 833 if not self.missing_rdc[align_index, j]: 834 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 835 836 # Calculate and sum the single alignment chi-squared value (for the RDC). 837 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 838 839 # PCS via numerical integration. 840 if self.pcs_flag_sum: 841 # Numerical integration of the PCSs. 842 pcs_numeric_int_iso_cone_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 843 844 # Calculate and sum the single alignment chi-squared value (for the PCS). 845 for align_index in range(self.num_align): 846 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 847 848 # Return the chi-squared value. 849 return chi2_sum

850 851

852 - def func_iso_cone_free_rotor(self, params):

853 """Target function for free rotor isotropic cone model optimisation. 854 855 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 856 857 858 @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. 859 @type params: list of float 860 @return: The chi-squared or SSE value. 861 @rtype: float 862 """ 863 864 # Scaling. 865 if self.scaling_flag: 866 params = dot(params, self.scaling_matrix) 867 868 # Unpack the parameters. 869 if self.translation_opt and self.pivot_opt: 870 self._param_pivot = params[:3] 871 self._translation_vector = params[3:6] 872 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[6:] 873 elif self.translation_opt: 874 self._translation_vector = params[:3] 875 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 876 elif self.pivot_opt: 877 self._param_pivot = params[:3] 878 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 879 else: 880 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params 881 882 # Generate the cone axis from the spherical angles. 883 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 884 885 # Pre-calculate the eigenframe rotation matrix. 886 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 887 888 # The Kronecker product of the eigenframe rotation. 889 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 890 891 # Calculate the cone angle. 892 theta_max = order_parameters.iso_cone_S_to_theta(cone_s1) 893 894 # Generate the 2nd degree Frame Order super matrix. 895 frame_order_2nd = compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_s1) 896 897 # Reduce and rotate the tensors. 898 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 899 900 # Pre-transpose matrices for faster calculations. 901 RT_eigen = transpose(self.R_eigen) 902 RT_ave = transpose(self.R_ave) 903 904 # Pre-calculate all the necessary vectors. 905 if self.pcs_flag: 906 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 907 908 # Initial chi-squared (or SSE) value. 909 chi2_sum = 0.0 910 911 # Loop over each alignment. 912 for align_index in range(self.num_align): 913 # RDCs. 914 if self.rdc_flag[align_index]: 915 # Loop over the RDCs. 916 for j in range(self.num_interatom): 917 # The back calculated RDC. 918 if not self.missing_rdc[align_index, j]: 919 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 920 921 # Calculate and sum the single alignment chi-squared value (for the RDC). 922 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 923 924 # PCS. 925 if self.pcs_flag[align_index]: 926 # Loop over the PCSs. 927 for j in range(self.num_spins): 928 # The back calculated PCS. 929 if not self.missing_pcs[align_index, j]: 930 # Forwards and reverse rotations. 931 if self.full_in_ref_frame[align_index]: 932 r_pivot_atom = self.r_pivot_atom[:, j] 933 else: 934 r_pivot_atom = self.r_pivot_atom_rev[:, j] 935 936 # The numerical integration. 937 self.pcs_theta[align_index, j] = pcs_numeric_int_iso_cone(theta_max=theta_max, sigma_max=pi, c=self.pcs_const[align_index], 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) 938 939 # Calculate and sum the single alignment chi-squared value (for the PCS). 940 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 941 942 # Return the chi-squared value. 943 return chi2_sum

944 945

946 - def func_iso_cone_free_rotor_qrint(self, params):

947 """Target function for free rotor isotropic cone model optimisation. 948 949 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 950 951 952 @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. 953 @type params: list of float 954 @return: The chi-squared or SSE value. 955 @rtype: float 956 """ 957 958 # Scaling. 959 if self.scaling_flag: 960 params = dot(params, self.scaling_matrix) 961 962 # Unpack the parameters. 963 if self.translation_opt and self.pivot_opt: 964 self._param_pivot = params[:3] 965 self._translation_vector = params[3:6] 966 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[6:] 967 elif self.translation_opt: 968 self._translation_vector = params[:3] 969 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 970 elif self.pivot_opt: 971 self._param_pivot = params[:3] 972 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 973 else: 974 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params 975 976 # Generate the cone axis from the spherical angles. 977 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 978 979 # Pre-calculate the eigenframe rotation matrix. 980 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 981 982 # The Kronecker product of the eigenframe rotation. 983 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 984 985 # Calculate the cone angle. 986 theta_max = order_parameters.iso_cone_S_to_theta(cone_s1) 987 988 # Generate the 2nd degree Frame Order super matrix. 989 frame_order_2nd = compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_s1) 990 991 # Reduce and rotate the tensors. 992 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 993 994 # Pre-transpose matrices for faster calculations. 995 RT_eigen = transpose(self.R_eigen) 996 RT_ave = transpose(self.R_ave) 997 998 # Pre-calculate all the necessary vectors. 999 if self.pcs_flag: 1000 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1001 1002 # Initial chi-squared (or SSE) value. 1003 chi2_sum = 0.0 1004 1005 # Loop over each alignment. 1006 for align_index in range(self.num_align): 1007 # RDCs. 1008 if self.rdc_flag[align_index]: 1009 # Loop over the RDCs. 1010 for j in range(self.num_interatom): 1011 # The back calculated RDC. 1012 if not self.missing_rdc[align_index, j]: 1013 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1014 1015 # Calculate and sum the single alignment chi-squared value (for the RDC). 1016 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1017 1018 # PCS via numerical integration. 1019 if self.pcs_flag_sum: 1020 # Numerical integration of the PCSs. 1021 pcs_numeric_int_iso_cone_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 1022 1023 # Calculate and sum the single alignment chi-squared value (for the PCS). 1024 for align_index in range(self.num_align): 1025 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1026 1027 # Return the chi-squared value. 1028 return chi2_sum

1029 1030

1031 - def func_iso_cone_torsionless(self, params):

1032 """Target function for torsionless isotropic cone model optimisation. 1033 1034 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1035 1036 1037 @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. 1038 @type params: list of float 1039 @return: The chi-squared or SSE value. 1040 @rtype: float 1041 """ 1042 1043 # Scaling. 1044 if self.scaling_flag: 1045 params = dot(params, self.scaling_matrix) 1046 1047 # Unpack the parameters. 1048 if self.translation_opt and self.pivot_opt: 1049 self._param_pivot = params[:3] 1050 self._translation_vector = params[3:6] 1051 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[6:] 1052 elif self.translation_opt: 1053 self._translation_vector = params[:3] 1054 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1055 elif self.pivot_opt: 1056 self._param_pivot = params[:3] 1057 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1058 else: 1059 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params 1060 1061 # Generate the cone axis from the spherical angles. 1062 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1063 1064 # Pre-calculate the eigenframe rotation matrix. 1065 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1066 1067 # The Kronecker product of the eigenframe rotation. 1068 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1069 1070 # Generate the 2nd degree Frame Order super matrix. 1071 frame_order_2nd = compile_2nd_matrix_iso_cone_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta) 1072 1073 # Reduce and rotate the tensors. 1074 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1075 1076 # Pre-transpose matrices for faster calculations. 1077 RT_eigen = transpose(self.R_eigen) 1078 RT_ave = transpose(self.R_ave) 1079 1080 # Pre-calculate all the necessary vectors. 1081 if self.pcs_flag: 1082 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1083 1084 # Initial chi-squared (or SSE) value. 1085 chi2_sum = 0.0 1086 1087 # Loop over each alignment. 1088 for align_index in range(self.num_align): 1089 # RDCs. 1090 if self.rdc_flag[align_index]: 1091 # Loop over the RDCs. 1092 for j in range(self.num_interatom): 1093 # The back calculated RDC. 1094 if not self.missing_rdc[align_index, j]: 1095 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1096 1097 # Calculate and sum the single alignment chi-squared value (for the RDC). 1098 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1099 1100 # PCS. 1101 if self.pcs_flag[align_index]: 1102 # Loop over the PCSs. 1103 for j in range(self.num_spins): 1104 # The back calculated PCS. 1105 if not self.missing_pcs[align_index, j]: 1106 # Forwards and reverse rotations. 1107 if self.full_in_ref_frame[align_index]: 1108 r_pivot_atom = self.r_pivot_atom[:, j] 1109 else: 1110 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1111 1112 # The numerical integration. 1113 self.pcs_theta[align_index, j] = pcs_numeric_int_iso_cone_torsionless(theta_max=cone_theta, c=self.pcs_const[align_index], 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) 1114 1115 # Calculate and sum the single alignment chi-squared value (for the PCS). 1116 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1117 1118 # Return the chi-squared value. 1119 return chi2_sum

1120 1121

1122 - def func_iso_cone_torsionless_qrint(self, params):

1123 """Target function for torsionless isotropic cone model optimisation. 1124 1125 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1126 1127 1128 @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. 1129 @type params: list of float 1130 @return: The chi-squared or SSE value. 1131 @rtype: float 1132 """ 1133 1134 # Scaling. 1135 if self.scaling_flag: 1136 params = dot(params, self.scaling_matrix) 1137 1138 # Unpack the parameters. 1139 if self.translation_opt and self.pivot_opt: 1140 self._param_pivot = params[:3] 1141 self._translation_vector = params[3:6] 1142 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[6:] 1143 elif self.translation_opt: 1144 self._translation_vector = params[:3] 1145 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1146 elif self.pivot_opt: 1147 self._param_pivot = params[:3] 1148 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1149 else: 1150 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params 1151 1152 # Generate the cone axis from the spherical angles. 1153 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1154 1155 # Pre-calculate the eigenframe rotation matrix. 1156 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1157 1158 # The Kronecker product of the eigenframe rotation. 1159 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1160 1161 # Generate the 2nd degree Frame Order super matrix. 1162 frame_order_2nd = compile_2nd_matrix_iso_cone_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta) 1163 1164 # Reduce and rotate the tensors. 1165 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1166 1167 # Pre-transpose matrices for faster calculations. 1168 RT_eigen = transpose(self.R_eigen) 1169 RT_ave = transpose(self.R_ave) 1170 1171 # Pre-calculate all the necessary vectors. 1172 if self.pcs_flag: 1173 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1174 1175 # Initial chi-squared (or SSE) value. 1176 chi2_sum = 0.0 1177 1178 # Loop over each alignment. 1179 for align_index in range(self.num_align): 1180 # RDCs. 1181 if self.rdc_flag[align_index]: 1182 # Loop over the RDCs. 1183 for j in range(self.num_interatom): 1184 # The back calculated RDC. 1185 if not self.missing_rdc[align_index, j]: 1186 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1187 1188 # Calculate and sum the single alignment chi-squared value (for the RDC). 1189 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1190 1191 # PCS via numerical integration. 1192 if self.pcs_flag_sum: 1193 # Numerical integration of the PCSs. 1194 pcs_numeric_int_iso_cone_torsionless_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 1195 1196 # Calculate and sum the single alignment chi-squared value (for the PCS). 1197 for align_index in range(self.num_align): 1198 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1199 1200 # Return the chi-squared value. 1201 return chi2_sum

1202 1203

1204 - def func_pseudo_ellipse(self, params):

1205 """Target function for pseudo-elliptic cone model optimisation. 1206 1207 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1208 1209 1210 @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. 1211 @type params: list of float 1212 @return: The chi-squared or SSE value. 1213 @rtype: float 1214 """ 1215 1216 # Scaling. 1217 if self.scaling_flag: 1218 params = dot(params, self.scaling_matrix) 1219 1220 # Unpack the parameters. 1221 if self.translation_opt and self.pivot_opt: 1222 self._param_pivot = params[:3] 1223 self._translation_vector = params[3:6] 1224 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:] 1225 elif self.translation_opt: 1226 self._translation_vector = params[:3] 1227 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:] 1228 elif self.pivot_opt: 1229 self._param_pivot = params[:3] 1230 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:] 1231 else: 1232 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 1233 1234 # Average position rotation. 1235 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1236 1237 # The Kronecker product of the eigenframe rotation. 1238 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1239 1240 # Generate the 2nd degree Frame Order super matrix. 1241 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y, cone_sigma_max) 1242 1243 # Reduce and rotate the tensors. 1244 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1245 1246 # Pre-transpose matrices for faster calculations. 1247 RT_eigen = transpose(self.R_eigen) 1248 RT_ave = transpose(self.R_ave) 1249 1250 # Pre-calculate all the necessary vectors. 1251 if self.pcs_flag: 1252 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1253 1254 # Initial chi-squared (or SSE) value. 1255 chi2_sum = 0.0 1256 1257 # Loop over each alignment. 1258 for align_index in range(self.num_align): 1259 # RDCs. 1260 if self.rdc_flag[align_index]: 1261 # Loop over the RDCs. 1262 for j in range(self.num_interatom): 1263 # The back calculated RDC. 1264 if not self.missing_rdc[align_index, j]: 1265 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1266 1267 # Calculate and sum the single alignment chi-squared value (for the RDC). 1268 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1269 1270 # PCS. 1271 if self.pcs_flag[align_index]: 1272 # Loop over the PCSs. 1273 for j in range(self.num_spins): 1274 # The back calculated PCS. 1275 if not self.missing_pcs[align_index, j]: 1276 # Forwards and reverse rotations. 1277 if self.full_in_ref_frame[align_index]: 1278 r_pivot_atom = self.r_pivot_atom[:, j] 1279 else: 1280 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1281 1282 # The numerical integration. 1283 self.pcs_theta[align_index, j] = pcs_numeric_int_pseudo_ellipse(theta_x=cone_theta_x, theta_y=cone_theta_y, sigma_max=cone_sigma_max, c=self.pcs_const[align_index], 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) 1284 1285 # Calculate and sum the single alignment chi-squared value (for the PCS). 1286 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1287 1288 # Return the chi-squared value. 1289 return chi2_sum

1290 1291

1292 - def func_pseudo_ellipse_qrint(self, params):

1293 """Target function for pseudo-elliptic cone model optimisation. 1294 1295 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. 1296 1297 1298 @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. 1299 @type params: list of float 1300 @return: The chi-squared or SSE value. 1301 @rtype: float 1302 """ 1303 1304 # Scaling. 1305 if self.scaling_flag: 1306 params = dot(params, self.scaling_matrix) 1307 1308 # Unpack the parameters. 1309 if self.translation_opt and self.pivot_opt: 1310 self._param_pivot = params[:3] 1311 self._translation_vector = params[3:6] 1312 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:] 1313 elif self.translation_opt: 1314 self._translation_vector = params[:3] 1315 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:] 1316 elif self.pivot_opt: 1317 self._param_pivot = params[:3] 1318 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:] 1319 else: 1320 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 1321 1322 # Average position rotation. 1323 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1324 1325 # The Kronecker product of the eigenframe rotation. 1326 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1327 1328 # Generate the 2nd degree Frame Order super matrix. 1329 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y, cone_sigma_max) 1330 1331 # Reduce and rotate the tensors. 1332 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1333 1334 # Pre-transpose matrices for faster calculations. 1335 RT_eigen = transpose(self.R_eigen) 1336 RT_ave = transpose(self.R_ave) 1337 1338 # Pre-calculate all the necessary vectors. 1339 if self.pcs_flag: 1340 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1341 1342 # Initial chi-squared (or SSE) value. 1343 chi2_sum = 0.0 1344 1345 # Loop over each alignment. 1346 for align_index in range(self.num_align): 1347 # RDCs. 1348 if self.rdc_flag[align_index]: 1349 # Loop over the RDCs. 1350 for j in range(self.num_interatom): 1351 # The back calculated RDC. 1352 if not self.missing_rdc[align_index, j]: 1353 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1354 1355 # Calculate and sum the single alignment chi-squared value (for the RDC). 1356 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1357 1358 # PCS via numerical integration. 1359 if self.pcs_flag_sum: 1360 # Numerical integration of the PCSs. 1361 pcs_numeric_int_pseudo_ellipse_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 1362 1363 # Calculate and sum the single alignment chi-squared value (for the PCS). 1364 for align_index in range(self.num_align): 1365 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1366 1367 # Return the chi-squared value. 1368 return chi2_sum

1369 1370

1371 - def func_pseudo_ellipse_free_rotor(self, params):

1372 """Target function for free_rotor pseudo-elliptic cone model optimisation. 1373 1374 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1375 1376 1377 @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. 1378 @type params: list of float 1379 @return: The chi-squared or SSE value. 1380 @rtype: float 1381 """ 1382 1383 # Scaling. 1384 if self.scaling_flag: 1385 params = dot(params, self.scaling_matrix) 1386 1387 # Unpack the parameters. 1388 if self.translation_opt and self.pivot_opt: 1389 self._param_pivot = params[:3] 1390 self._translation_vector = params[3:6] 1391 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1392 elif self.translation_opt: 1393 self._translation_vector = params[:3] 1394 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1395 elif self.pivot_opt: 1396 self._param_pivot = params[:3] 1397 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1398 else: 1399 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1400 1401 # Average position rotation. 1402 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1403 1404 # The Kronecker product of the eigenframe rotation. 1405 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1406 1407 # Generate the 2nd degree Frame Order super matrix. 1408 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1409 1410 # Reduce and rotate the tensors. 1411 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1412 1413 # Pre-transpose matrices for faster calculations. 1414 RT_eigen = transpose(self.R_eigen) 1415 RT_ave = transpose(self.R_ave) 1416 1417 # Pre-calculate all the necessary vectors. 1418 if self.pcs_flag: 1419 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1420 1421 # Initial chi-squared (or SSE) value. 1422 chi2_sum = 0.0 1423 1424 # Loop over each alignment. 1425 for align_index in range(self.num_align): 1426 # RDCs. 1427 if self.rdc_flag[align_index]: 1428 # Loop over the RDCs. 1429 for j in range(self.num_interatom): 1430 # The back calculated RDC. 1431 if not self.missing_rdc[align_index, j]: 1432 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1433 1434 # Calculate and sum the single alignment chi-squared value (for the RDC). 1435 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1436 1437 # PCS. 1438 if self.pcs_flag[align_index]: 1439 # Loop over the PCSs. 1440 for j in range(self.num_spins): 1441 # The back calculated PCS. 1442 if not self.missing_pcs[align_index, j]: 1443 # Forwards and reverse rotations. 1444 if self.full_in_ref_frame[align_index]: 1445 r_pivot_atom = self.r_pivot_atom[:, j] 1446 else: 1447 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1448 1449 # The numerical integration. 1450 self.pcs_theta[align_index, j] = pcs_numeric_int_pseudo_ellipse(theta_x=cone_theta_x, theta_y=cone_theta_y, sigma_max=pi, c=self.pcs_const[align_index], 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) 1451 1452 # Calculate and sum the single alignment chi-squared value (for the PCS). 1453 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1454 1455 # Return the chi-squared value. 1456 return chi2_sum

1457 1458

1459 - def func_pseudo_ellipse_free_rotor_qrint(self, params):

1460 """Target function for free_rotor pseudo-elliptic cone model optimisation. 1461 1462 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1463 1464 1465 @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. 1466 @type params: list of float 1467 @return: The chi-squared or SSE value. 1468 @rtype: float 1469 """ 1470 1471 # Scaling. 1472 if self.scaling_flag: 1473 params = dot(params, self.scaling_matrix) 1474 1475 # Unpack the parameters. 1476 if self.translation_opt and self.pivot_opt: 1477 self._param_pivot = params[:3] 1478 self._translation_vector = params[3:6] 1479 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1480 elif self.translation_opt: 1481 self._translation_vector = params[:3] 1482 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1483 elif self.pivot_opt: 1484 self._param_pivot = params[:3] 1485 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1486 else: 1487 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 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(ave_pos_alpha, 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(self._param_pivot, self.R_ave, RT_ave) 1508 1509 # Initial chi-squared (or SSE) value. 1510 chi2_sum = 0.0 1511 1512 # Loop over each alignment. 1513 for align_index in range(self.num_align): 1514 # RDCs. 1515 if self.rdc_flag[align_index]: 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_sum: 1527 # Numerical integration of the PCSs. 1528 pcs_numeric_int_pseudo_ellipse_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 1529 1530 # Calculate and sum the single alignment chi-squared value (for the PCS). 1531 for align_index in range(self.num_align): 1532 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1533 1534 # Return the chi-squared value. 1535 return chi2_sum

1536 1537

1538 - def func_pseudo_ellipse_torsionless(self, params):

1539 """Target function for torsionless pseudo-elliptic cone model optimisation. 1540 1541 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1542 1543 1544 @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. 1545 @type params: list of float 1546 @return: The chi-squared or SSE value. 1547 @rtype: float 1548 """ 1549 1550 # Scaling. 1551 if self.scaling_flag: 1552 params = dot(params, self.scaling_matrix) 1553 1554 # Unpack the parameters. 1555 if self.translation_opt and self.pivot_opt: 1556 self._param_pivot = params[:3] 1557 self._translation_vector = params[3:6] 1558 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1559 elif self.translation_opt: 1560 self._translation_vector = params[:3] 1561 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1562 elif self.pivot_opt: 1563 self._param_pivot = params[:3] 1564 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1565 else: 1566 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1567 1568 # Average position rotation. 1569 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1570 1571 # The Kronecker product of the eigenframe rotation. 1572 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1573 1574 # Generate the 2nd degree Frame Order super matrix. 1575 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1576 1577 # Reduce and rotate the tensors. 1578 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1579 1580 # Pre-transpose matrices for faster calculations. 1581 RT_eigen = transpose(self.R_eigen) 1582 RT_ave = transpose(self.R_ave) 1583 1584 # Pre-calculate all the necessary vectors. 1585 if self.pcs_flag: 1586 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1587 1588 # Initial chi-squared (or SSE) value. 1589 chi2_sum = 0.0 1590 1591 # Loop over each alignment. 1592 for align_index in range(self.num_align): 1593 # RDCs. 1594 if self.rdc_flag[align_index]: 1595 # Loop over the RDCs. 1596 for j in range(self.num_interatom): 1597 # The back calculated RDC. 1598 if not self.missing_rdc[align_index, j]: 1599 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1600 1601 # Calculate and sum the single alignment chi-squared value (for the RDC). 1602 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1603 1604 # PCS. 1605 if self.pcs_flag[align_index]: 1606 # Loop over the PCSs. 1607 for j in range(self.num_spins): 1608 # The back calculated PCS. 1609 if not self.missing_pcs[align_index, j]: 1610 # Forwards and reverse rotations. 1611 if self.full_in_ref_frame[align_index]: 1612 r_pivot_atom = self.r_pivot_atom[:, j] 1613 else: 1614 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1615 1616 # The numerical integration. 1617 self.pcs_theta[align_index, j] = pcs_numeric_int_pseudo_ellipse_torsionless(theta_x=cone_theta_x, theta_y=cone_theta_y, c=self.pcs_const[align_index], 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) 1618 1619 # Calculate and sum the single alignment chi-squared value (for the PCS). 1620 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1621 1622 # Return the chi-squared value. 1623 return chi2_sum

1624 1625

1626 - def func_pseudo_ellipse_torsionless_qrint(self, params):

1627 """Target function for torsionless pseudo-elliptic cone model optimisation. 1628 1629 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1630 1631 1632 @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. 1633 @type params: list of float 1634 @return: The chi-squared or SSE value. 1635 @rtype: float 1636 """ 1637 1638 # Scaling. 1639 if self.scaling_flag: 1640 params = dot(params, self.scaling_matrix) 1641 1642 # Unpack the parameters. 1643 if self.translation_opt and self.pivot_opt: 1644 self._param_pivot = params[:3] 1645 self._translation_vector = params[3:6] 1646 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1647 elif self.translation_opt: 1648 self._translation_vector = params[:3] 1649 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1650 elif self.pivot_opt: 1651 self._param_pivot = params[:3] 1652 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1653 else: 1654 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1655 1656 # Average position rotation. 1657 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1658 1659 # The Kronecker product of the eigenframe rotation. 1660 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1661 1662 # Generate the 2nd degree Frame Order super matrix. 1663 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1664 1665 # Reduce and rotate the tensors. 1666 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1667 1668 # Pre-transpose matrices for faster calculations. 1669 RT_eigen = transpose(self.R_eigen) 1670 RT_ave = transpose(self.R_ave) 1671 1672 # Pre-calculate all the necessary vectors. 1673 if self.pcs_flag: 1674 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1675 1676 # Initial chi-squared (or SSE) value. 1677 chi2_sum = 0.0 1678 1679 # Loop over each alignment. 1680 for align_index in range(self.num_align): 1681 # RDCs. 1682 if self.rdc_flag[align_index]: 1683 # Loop over the RDCs. 1684 for j in range(self.num_interatom): 1685 # The back calculated RDC. 1686 if not self.missing_rdc[align_index, j]: 1687 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1688 1689 # Calculate and sum the single alignment chi-squared value (for the RDC). 1690 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1691 1692 # PCS via numerical integration. 1693 if self.pcs_flag_sum: 1694 # Numerical integration of the PCSs. 1695 pcs_numeric_int_pseudo_ellipse_torsionless_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 1696 1697 # Calculate and sum the single alignment chi-squared value (for the PCS). 1698 for align_index in range(self.num_align): 1699 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1700 1701 # Return the chi-squared value. 1702 return chi2_sum

1703 1704

1705 - def func_rigid(self, params):

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

1776 1777

1778 - def func_rotor(self, params):

1779 """Target function for rotor model optimisation. 1780 1781 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1782 1783 1784 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 1785 @type params: list of float 1786 @return: The chi-squared or SSE value. 1787 @rtype: float 1788 """ 1789 1790 # Scaling. 1791 if self.scaling_flag: 1792 params = dot(params, self.scaling_matrix) 1793 1794 # Unpack the parameters. 1795 if self.translation_opt and self.pivot_opt: 1796 self._param_pivot = params[:3] 1797 self._translation_vector = params[3:6] 1798 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[6:] 1799 elif self.translation_opt: 1800 self._translation_vector = params[:3] 1801 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[3:] 1802 elif self.pivot_opt: 1803 self._param_pivot = params[:3] 1804 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[3:] 1805 else: 1806 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params 1807 1808 # Generate the cone axis from the spherical angles. 1809 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1810 1811 # Pre-calculate the eigenframe rotation matrix. 1812 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1813 1814 # The Kronecker product of the eigenframe rotation. 1815 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1816 1817 # Generate the 2nd degree Frame Order super matrix. 1818 frame_order_2nd = compile_2nd_matrix_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max) 1819 1820 # The average frame rotation matrix (and reduce and rotate the tensors). 1821 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1822 1823 # Pre-transpose matrices for faster calculations. 1824 RT_eigen = transpose(self.R_eigen) 1825 RT_ave = transpose(self.R_ave) 1826 1827 # Pre-calculate all the necessary vectors. 1828 if self.pcs_flag: 1829 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1830 1831 # Initial chi-squared (or SSE) value. 1832 chi2_sum = 0.0 1833 1834 # Loop over each alignment. 1835 for align_index in range(self.num_align): 1836 # RDCs. 1837 if self.rdc_flag[align_index]: 1838 # Loop over the RDCs. 1839 for j in range(self.num_interatom): 1840 # The back calculated RDC. 1841 if not self.missing_rdc[align_index, j]: 1842 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1843 1844 # Calculate and sum the single alignment chi-squared value (for the RDC). 1845 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1846 1847 # PCS. 1848 if self.pcs_flag[align_index]: 1849 # Loop over the PCSs. 1850 for j in range(self.num_spins): 1851 # The back calculated PCS. 1852 if not self.missing_pcs[align_index, j]: 1853 # Forwards and reverse rotations. 1854 if self.full_in_ref_frame[align_index]: 1855 r_pivot_atom = self.r_pivot_atom[:, j] 1856 else: 1857 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1858 1859 # The numerical integration. 1860 self.pcs_theta[align_index, j] = pcs_numeric_int_rotor(sigma_max=sigma_max, c=self.pcs_const[align_index], 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) 1861 1862 # Calculate and sum the single alignment chi-squared value (for the PCS). 1863 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1864 1865 # Return the chi-squared value. 1866 return chi2_sum

1867 1868

1869 - def func_rotor_qrint(self, params):

1870 """Target function for rotor model optimisation. 1871 1872 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. 1873 1874 1875 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 1876 @type params: list of float 1877 @return: The chi-squared or SSE value. 1878 @rtype: float 1879 """ 1880 1881 # Scaling. 1882 if self.scaling_flag: 1883 params = dot(params, self.scaling_matrix) 1884 1885 # Unpack the parameters. 1886 if self.translation_opt and self.pivot_opt: 1887 self._param_pivot = params[:3] 1888 self._translation_vector = params[3:6] 1889 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[6:] 1890 elif self.translation_opt: 1891 self._translation_vector = params[:3] 1892 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[3:] 1893 elif self.pivot_opt: 1894 self._param_pivot = params[:3] 1895 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params[3:] 1896 else: 1897 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_alpha, sigma_max = params 1898 1899 # Generate the rotor axis. 1900 self.cone_axis = create_rotor_axis_alpha(alpha=axis_alpha, pivot=array(self._param_pivot, float64), point=self.com) 1901 1902 # Pre-calculate the eigenframe rotation matrix. 1903 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1904 1905 # The Kronecker product of the eigenframe rotation. 1906 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1907 1908 # Generate the 2nd degree Frame Order super matrix. 1909 frame_order_2nd = compile_2nd_matrix_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max) 1910 1911 # The average frame rotation matrix (and reduce and rotate the tensors). 1912 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1913 1914 # Pre-transpose matrices for faster calculations. 1915 RT_eigen = transpose(self.R_eigen) 1916 RT_ave = transpose(self.R_ave) 1917 1918 # Pre-calculate all the necessary vectors. 1919 if self.pcs_flag: 1920 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1921 1922 # Initial chi-squared (or SSE) value. 1923 chi2_sum = 0.0 1924 1925 # Loop over each alignment. 1926 for align_index in range(self.num_align): 1927 # RDCs. 1928 if self.rdc_flag[align_index]: 1929 # Loop over the RDCs. 1930 for j in range(self.num_interatom): 1931 # The back calculated RDC. 1932 if not self.missing_rdc[align_index, j]: 1933 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1934 1935 # Calculate and sum the single alignment chi-squared value (for the RDC). 1936 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1937 1938 # PCS via numerical integration. 1939 if self.pcs_flag_sum: 1940 # Numerical integration of the PCSs. 1941 pcs_numeric_int_rotor_qrint(points=self.sobol_angles, 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=self.Ri_prime, pcs_theta=self.pcs_theta, pcs_theta_err=self.pcs_theta_err, missing_pcs=self.missing_pcs, error_flag=False) 1942 1943 # Calculate and sum the single alignment chi-squared value (for the PCS). 1944 for align_index in range(self.num_align): 1945 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1946 1947 # Return the chi-squared value. 1948 return chi2_sum

1949 1950

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

1952 """Calculate the pivot to atom and lanthanide to pivot vectors for the target functions. 1953 1954 @keyword pivot: The pivot point. 1955 @type pivot: numpy rank-1, 3D array 1956 @keyword R_ave: The rotation matrix for rotating from the reference frame to the average position. 1957 @type R_ave: numpy rank-2, 3D array 1958 @keyword RT_ave: The transpose of R_ave. 1959 @type RT_ave: numpy rank-2, 3D array 1960 """ 1961 1962 # The rotational pivot. 1963 if self.ave_pos_piv_sync: 1964 ave_pos_pivot = self._param_pivot 1965 else: 1966 ave_pos_pivot = self.ave_pos_pivot 1967 1968 # The pivot to atom vectors. 1969 for j in range(self.num_spins): 1970 # The lanthanide to pivot vector. 1971 if self.pivot_opt: 1972 self.r_ln_pivot[:, j] = pivot - self.paramag_centre 1973 1974 # Rotate then translate the atomic positions, then calculate the pivot to atom vector. 1975 self.r_pivot_atom[:, j] = dot(R_ave, self.atomic_pos[j] - ave_pos_pivot) + ave_pos_pivot 1976 self.r_pivot_atom[:, j] += self._translation_vector 1977 self.r_pivot_atom[:, j] -= pivot 1978 self.r_pivot_atom_rev[:, j] = dot(RT_ave, self.atomic_pos[j] - ave_pos_pivot) + ave_pos_pivot 1979 self.r_pivot_atom_rev[:, j] += self._translation_vector 1980 self.r_pivot_atom_rev[:, j] -= pivot

1981 1982

1983 - def create_sobol_data(self, n=10000, dims=None):

1984 """Create the Sobol' quasi-random data for numerical integration. 1985 1986 This uses the external sobol_lib module to create the data. The algorithm is that modified by Antonov and Saleev. 1987 1988 1989 @keyword n: The number of points to generate. 1990 @type n: int 1991 @keyword dims: The list of parameters. 1992 @type dims: list of str 1993 """ 1994 1995 # The number of dimensions. 1996 m = len(dims) 1997 1998 # Initialise. 1999 self.sobol_angles = zeros((n, m), float32) 2000 2001 # Loop over the points. 2002 for i in range(n): 2003 # The raw point. 2004 point, seed = i4_sobol(m, i) 2005 2006 # Loop over the dimensions, converting the points to angles. 2007 for j in range(m): 2008 # The tilt angle - the angle of rotation about the x-y plane rotation axis. 2009 if dims[j] in ['theta']: 2010 self.sobol_angles[i, j] = acos(2.0*point[j] - 1.0) 2011 2012 # The angle defining the x-y plane rotation axis. 2013 if dims[j] in ['phi']: 2014 self.sobol_angles[i, j] = 2.0 * pi * point[j] 2015 2016 # The torsion angle - the angle of rotation about the z' axis. 2017 if dims[j] in ['sigma', 'sigma2']: 2018 self.sobol_angles[i, j] = 2.0 * pi * (point[j] - 0.5)

2019 2020

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

2022 """Reduce and rotate the alignments tensors using the frame order matrix and Euler angles. 2023 2024 @keyword ave_pos_alpha: The alpha Euler angle describing the average domain position, the tensor rotation. 2025 @type ave_pos_alpha: float 2026 @keyword ave_pos_beta: The beta Euler angle describing the average domain position, the tensor rotation. 2027 @type ave_pos_beta: float 2028 @keyword ave_pos_gamma: The gamma Euler angle describing the average domain position, the tensor rotation. 2029 @type ave_pos_gamma: float 2030 @keyword daeg: The 2nd degree frame order matrix. 2031 @type daeg: rank-2, 9D array 2032 """ 2033 2034 # Alignment tensor rotation. 2035 euler_to_R_zyz(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, self.R_ave) 2036 2037 # Back calculate the rotated tensors. 2038 for align_index in range(self.num_tensors): 2039 # Tensor indices. 2040 index1 = align_index*5 2041 index2 = align_index*5+5 2042 2043 # Reduction. 2044 if daeg != None: 2045 # Reduce the tensor. 2046 reduce_alignment_tensor(daeg, self.full_tensors[index1:index2], self.A_5D_bc[index1:index2]) 2047 2048 # Convert the reduced tensor to 3D, rank-2 form. 2049 to_tensor(self.tensor_3D, self.A_5D_bc[index1:index2]) 2050 2051 # No reduction: 2052 else: 2053 # Convert the original tensor to 3D, rank-2 form. 2054 to_tensor(self.tensor_3D, self.full_tensors[index1:index2]) 2055 2056 # Rotate the tensor (normal R.X.RT rotation). 2057 if self.full_in_ref_frame[align_index]: 2058 self.A_3D_bc[align_index] = dot(transpose(self.R_ave), dot(self.tensor_3D, self.R_ave)) 2059 2060 # Rotate the tensor (inverse RT.X.R rotation). 2061 else: 2062 self.A_3D_bc[align_index] = dot(self.R_ave, dot(self.tensor_3D, transpose(self.R_ave))) 2063 2064 # Convert the tensor back to 5D, rank-1 form, as the back-calculated reduced tensor. 2065 to_5D(self.A_5D_bc[index1:index2], self.A_3D_bc[align_index])

Source Code for Module target_functions.frame_order