Package target_functions :: Module frame_order
[hide private]
[frames] | no frames]

Source Code for Module target_functions.frame_order

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2009-2013 Edward d'Auvergne                                   # 
   4  #                                                                             # 
   5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
   6  #                                                                             # 
   7  # This program is free software: you can redistribute it and/or modify        # 
   8  # it under the terms of the GNU General Public License as published by        # 
   9  # the Free Software Foundation, either version 3 of the License, or           # 
  10  # (at your option) any later version.                                         # 
  11  #                                                                             # 
  12  # This program is distributed in the hope that it will be useful,             # 
  13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
  14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
  15  # GNU General Public License for more details.                                # 
  16  #                                                                             # 
  17  # You should have received a copy of the GNU General Public License           # 
  18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
  19  #                                                                             # 
  20  ############################################################################### 
  21   
  22  # Module docstring. 
  23  """Module containing the target functions of the Frame Order theories.""" 
  24   
  25  # Python module imports. 
  26  from copy import deepcopy 
  27  from math import acos, pi, sqrt 
  28  from numpy import array, dot, float32, float64, ones, transpose, uint8, zeros 
  29  from numpy.linalg import norm 
  30   
  31  # relax module imports. 
  32  from extern.sobol.sobol_lib import i4_sobol 
  33  from lib.alignment.alignment_tensor import to_5D, to_tensor 
  34  from lib.alignment.pcs import pcs_tensor 
  35  from lib.alignment.rdc import rdc_tensor 
  36  from lib.errors import RelaxError 
  37  from lib.float import isNaN 
  38  from lib.frame_order.free_rotor import compile_2nd_matrix_free_rotor 
  39  from lib.frame_order.iso_cone import compile_2nd_matrix_iso_cone, pcs_numeric_int_iso_cone, pcs_numeric_int_iso_cone_qrint 
  40  from lib.frame_order.iso_cone_free_rotor import compile_2nd_matrix_iso_cone_free_rotor 
  41  from lib.frame_order.iso_cone_torsionless import compile_2nd_matrix_iso_cone_torsionless, pcs_numeric_int_iso_cone_torsionless, pcs_numeric_int_iso_cone_torsionless_qrint 
  42  from lib.frame_order.matrix_ops import reduce_alignment_tensor 
  43  from lib.frame_order.pseudo_ellipse import compile_2nd_matrix_pseudo_ellipse, pcs_numeric_int_pseudo_ellipse, pcs_numeric_int_pseudo_ellipse_qrint 
  44  from lib.frame_order.pseudo_ellipse_free_rotor import compile_2nd_matrix_pseudo_ellipse_free_rotor 
  45  from lib.frame_order.pseudo_ellipse_torsionless import compile_2nd_matrix_pseudo_ellipse_torsionless, pcs_numeric_int_pseudo_ellipse_torsionless, pcs_numeric_int_pseudo_ellipse_torsionless_qrint 
  46  from lib.frame_order.rotor import compile_2nd_matrix_rotor, pcs_numeric_int_rotor, pcs_numeric_int_rotor_qrint 
  47  from lib.geometry.coord_transform import spherical_to_cartesian 
  48  from lib.geometry.rotations import euler_to_R_zyz, two_vect_to_R 
  49  from lib.linear_algebra.kronecker_product import kron_prod 
  50  from lib.order import order_parameters 
  51  from lib.physical_constants import pcs_constant 
  52  from target_functions.chi2 import chi2 
  53   
  54   
55 -class Frame_order:
56 """Class containing the target function of the optimisation of Frame Order matrix components.""" 57
58 - 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, ave_pos_pivot=zeros(3), ave_pos_piv_sync=True, translation_opt=False, pivot=zeros(3), pivot_opt=False, quad_int=True):
59 """Set up the target functions for the Frame Order theories. 60 61 @keyword model: The name of the Frame Order model. 62 @type model: str 63 @keyword init_params: The initial parameter values. 64 @type init_params: numpy float64 array 65 @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]. 66 @type full_tensors: numpy nx5D, rank-1 float64 array 67 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 68 @type full_in_ref_frame: numpy rank-1 array 69 @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. 70 @type rdcs: numpy rank-2 array 71 @keyword rdc_errors: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 72 @type rdc_errors: numpy rank-2 array 73 @keyword rdc_weights: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 74 @type rdc_weights: numpy rank-2 array 75 @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. 76 @type rdc_vect: numpy rank-2 array 77 @keyword dip_const: The dipolar constants for each RDC. The indices correspond to the spin systems j. 78 @type dip_const: numpy rank-1 array 79 @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. 80 @type pcs: numpy rank-2 array 81 @keyword pcs_errors: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 82 @type pcs_errors: numpy rank-2 array 83 @keyword pcs_weights: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 84 @type pcs_weights: numpy rank-2 array 85 @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. 86 @type atomic_pos: numpy rank-3 array 87 @keyword temp: The temperature of each PCS data set. 88 @type temp: numpy rank-1 array 89 @keyword frq: The frequency of each PCS data set. 90 @type frq: numpy rank-1 array 91 @keyword paramag_centre: The paramagnetic centre position (or positions). 92 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 93 @keyword scaling_matrix: The square and diagonal scaling matrix. 94 @type scaling_matrix: numpy rank-2 array 95 @keyword num_int_pts: The number of points to use for the numerical integration technique. 96 @type num_int_pts: int 97 @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. 98 @type ave_pos_pivot: numpy 3D rank-1 array 99 @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. 100 @type ave_pos_piv_sync: bool 101 @keyword translation_opt: A flag which if True will allow the pivot point of the motion to be optimised. 102 @type translation_opt: bool 103 @keyword pivot: The pivot point for the ball-and-socket joint motion. This is needed if PCS or PRE values are used. 104 @type pivot: numpy rank-1, 3D array or None 105 @keyword pivot_opt: A flag which if True will allow the pivot point of the motion to be optimised. 106 @type pivot_opt: bool 107 @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. 108 @type quad_int: bool 109 """ 110 111 # Model test. 112 if not model: 113 raise RelaxError("The type of Frame Order model must be specified.") 114 115 # Store the initial parameter (as a copy). 116 self.params = deepcopy(init_params) 117 118 # Store the agrs. 119 self.model = model 120 self.full_tensors = full_tensors 121 self.full_in_ref_frame = full_in_ref_frame 122 self.rdc = rdcs 123 self.rdc_weights = rdc_weights 124 self.rdc_vect = rdc_vect 125 self.dip_const = dip_const 126 self.pcs = pcs 127 self.pcs_weights = pcs_weights 128 self.atomic_pos = atomic_pos 129 self.temp = temp 130 self.frq = frq 131 self.paramag_centre = paramag_centre 132 self.total_num_params = len(init_params) 133 self.num_int_pts = num_int_pts 134 self.ave_pos_pivot = ave_pos_pivot 135 self.ave_pos_piv_sync = ave_pos_piv_sync 136 self.translation_opt = translation_opt 137 self._param_pivot = pivot 138 self.pivot_opt = pivot_opt 139 140 # Tensor setup. 141 self._init_tensors() 142 143 # Scaling initialisation. 144 self.scaling_matrix = scaling_matrix 145 if self.scaling_matrix != 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 != None: 153 self.num_align = len(rdcs) 154 elif pcs != None: 155 self.num_align = len(pcs) 156 157 # Set the RDC and PCS flags (indicating the presence of data). 158 self.rdc_flag = [True] * self.num_align 159 self.pcs_flag = [True] * self.num_align 160 for align_index in range(self.num_align): 161 if rdcs == None or len(rdcs[align_index]) == 0: 162 self.rdc_flag[align_index] = False 163 if pcs == None or len(pcs[align_index]) == 0: 164 self.pcs_flag[align_index] = False 165 self.rdc_flag_sum = sum(self.rdc_flag) 166 self.pcs_flag_sum = sum(self.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_sum and (rdc_vect == None or not len(rdc_vect)): 173 raise RelaxError("The rdc_vect argument " + repr(rdc_vect) + " must be supplied.") 174 if self.pcs_flag_sum and (atomic_pos == 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_sum: 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_sum: 185 self.num_interatom = len(rdcs[0]) 186 187 # Set up the alignment data. 188 for align_index in range(self.num_align): 189 to_tensor(self.A_3D[align_index], self.full_tensors[5*align_index:5*align_index+5]) 190 191 # PCS errors. 192 if self.pcs_flag_sum: 193 err = False 194 for i in range(len(pcs_errors)): 195 for j in range(len(pcs_errors[i])): 196 if not isNaN(pcs_errors[i, j]): 197 err = True 198 if err: 199 self.pcs_error = pcs_errors 200 else: 201 # Missing errors (default to 0.1 ppm errors). 202 self.pcs_error = 0.1 * 1e-6 * ones((self.num_align, self.num_spins), float64) 203 204 # RDC errors. 205 if self.rdc_flag_sum: 206 err = False 207 for i in range(len(rdc_errors)): 208 for j in range(len(rdc_errors[i])): 209 if not isNaN(rdc_errors[i, j]): 210 err = True 211 if err: 212 self.rdc_error = rdc_errors 213 else: 214 # Missing errors (default to 1 Hz errors). 215 self.rdc_error = ones((self.num_align, self.num_interatom), float64) 216 217 # Missing data matrices (RDC). 218 if self.rdc_flag_sum: 219 self.missing_rdc = zeros((self.num_align, self.num_interatom), uint8) 220 221 # Missing data matrices (PCS). 222 if self.pcs_flag_sum: 223 self.missing_pcs = zeros((self.num_align, self.num_spins), uint8) 224 225 # Clean up problematic data and put the weights into the errors.. 226 if self.rdc_flag_sum or self.pcs_flag_sum: 227 for align_index in range(self.num_align): 228 # Loop over the RDCs. 229 if self.rdc_flag_sum: 230 for j in range(self.num_interatom): 231 if isNaN(self.rdc[align_index, j]): 232 # Set the flag. 233 self.missing_rdc[align_index, j] = 1 234 235 # Change the NaN to zero. 236 self.rdc[align_index, j] = 0.0 237 238 # Change the error to one (to avoid zero division). 239 self.rdc_error[align_index, j] = 1.0 240 241 # Change the weight to one. 242 rdc_weights[align_index, j] = 1.0 243 244 # The RDC weights. 245 if self.rdc_flag_sum: 246 self.rdc_error[align_index, j] = self.rdc_error[align_index, j] / sqrt(rdc_weights[align_index, j]) 247 248 # Loop over the PCSs. 249 if self.pcs_flag_sum: 250 for j in range(self.num_spins): 251 if isNaN(self.pcs[align_index, j]): 252 # Set the flag. 253 self.missing_pcs[align_index, j] = 1 254 255 # Change the NaN to zero. 256 self.pcs[align_index, j] = 0.0 257 258 # Change the error to one (to avoid zero division). 259 self.pcs_error[align_index, j] = 1.0 260 261 # Change the weight to one. 262 pcs_weights[align_index, j] = 1.0 263 264 # The PCS weights. 265 if self.pcs_flag_sum: 266 self.pcs_error[align_index, j] = self.pcs_error[align_index, j] / sqrt(pcs_weights[align_index, j]) 267 268 # The paramagnetic centre vectors and distances. 269 if self.pcs_flag_sum: 270 # Initialise the data structures. 271 self.paramag_unit_vect = zeros(atomic_pos.shape, float64) 272 self.paramag_dist = zeros(self.num_spins, float64) 273 self.pcs_const = zeros(self.num_align, float64) 274 self.r_pivot_atom = zeros((3, self.num_spins), float32) 275 self.r_pivot_atom_rev = zeros((3, self.num_spins), float32) 276 self.r_ln_pivot = zeros((3, self.num_spins), float32) 277 for j in range(self.num_spins): 278 self.r_ln_pivot[:, j] = pivot - self.paramag_centre 279 if self.paramag_centre == None: 280 self.paramag_centre = zeros(3, float32) 281 282 # Set up the paramagnetic constant (without the interatomic distance and in Angstrom units). 283 for align_index in range(self.num_align): 284 self.pcs_const[align_index] = pcs_constant(self.temp[align_index], self.frq[align_index], 1.0) * 1e30 285 286 # PCS function, gradient, and Hessian matrices. 287 if self.pcs_flag_sum: 288 self.pcs_theta = zeros((self.num_align, self.num_spins), float64) 289 self.pcs_theta_err = zeros((self.num_align, self.num_spins), float64) 290 self.dpcs_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 291 self.d2pcs_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 292 293 # RDC function, gradient, and Hessian matrices. 294 if self.rdc_flag_sum: 295 self.rdc_theta = zeros((self.num_align, self.num_interatom), float64) 296 self.drdc_theta = zeros((self.total_num_params, self.num_align, self.num_interatom), float64) 297 self.d2rdc_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_interatom), float64) 298 299 # The quasi-random integration. 300 if not quad_int and self.pcs_flag_sum and model not in ['rigid']: 301 # The Sobol' sequence data and target function aliases (quasi-random integration). 302 if model == 'pseudo-ellipse': 303 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 304 self.func = self.func_pseudo_ellipse_qrint 305 elif model == 'pseudo-ellipse, torsionless': 306 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi']) 307 self.func = self.func_pseudo_ellipse_torsionless_qrint 308 elif model == 'pseudo-ellipse, free rotor': 309 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 310 self.func = self.func_pseudo_ellipse_free_rotor_qrint 311 elif model == 'iso cone': 312 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 313 self.func = self.func_iso_cone_qrint 314 elif model == 'iso cone, torsionless': 315 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi']) 316 self.func = self.func_iso_cone_torsionless_qrint 317 elif model == 'iso cone, free rotor': 318 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'phi', 'sigma']) 319 self.func = self.func_iso_cone_free_rotor_qrint 320 elif model == 'line': 321 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'sigma']) 322 self.func = self.func_line_qrint 323 elif model == 'line, torsionless': 324 self.create_sobol_data(n=self.num_int_pts, dims=['theta']) 325 self.func = self.func_line_torsionless_qrint 326 elif model == 'line, free rotor': 327 self.create_sobol_data(n=self.num_int_pts, dims=['theta', 'sigma']) 328 self.func = self.func_line_free_rotor_qrint 329 elif model == 'rotor': 330 self.create_sobol_data(n=self.num_int_pts, dims=['sigma']) 331 self.func = self.func_rotor_qrint 332 elif model == 'rigid': 333 self.func = self.func_rigid 334 elif model == 'free rotor': 335 self.create_sobol_data(n=self.num_int_pts, dims=['sigma']) 336 self.func = self.func_free_rotor_qrint 337 338 # The target function aliases (Scipy numerical integration). 339 else: 340 if model == 'pseudo-ellipse': 341 self.func = self.func_pseudo_ellipse 342 elif model == 'pseudo-ellipse, torsionless': 343 self.func = self.func_pseudo_ellipse_torsionless 344 elif model == 'pseudo-ellipse, free rotor': 345 self.func = self.func_pseudo_ellipse_free_rotor 346 elif model == 'iso cone': 347 self.func = self.func_iso_cone 348 elif model == 'iso cone, torsionless': 349 self.func = self.func_iso_cone_torsionless 350 elif model == 'iso cone, free rotor': 351 self.func = self.func_iso_cone_free_rotor 352 elif model == 'line': 353 self.func = self.func_line 354 elif model == 'line, torsionless': 355 self.func = self.func_line_torsionless 356 elif model == 'line, free rotor': 357 self.func = self.func_line_free_rotor 358 elif model == 'rotor': 359 self.func = self.func_rotor 360 elif model == 'rigid': 361 self.func = self.func_rigid 362 elif model == 'free rotor': 363 self.func = self.func_free_rotor
364 365
366 - def _init_tensors(self):
367 """Set up isotropic cone optimisation against the alignment tensor data.""" 368 369 # Some checks. 370 if self.full_tensors == None or not len(self.full_tensors): 371 raise RelaxError("The full_tensors argument " + repr(self.full_tensors) + " must be supplied.") 372 if self.full_in_ref_frame == None or not len(self.full_in_ref_frame): 373 raise RelaxError("The full_in_ref_frame argument " + repr(self.full_in_ref_frame) + " must be supplied.") 374 375 # Tensor set up. 376 self.num_tensors = int(len(self.full_tensors) / 5) 377 self.A_3D = zeros((self.num_tensors, 3, 3), float64) 378 self.A_3D_bc = zeros((self.num_tensors, 3, 3), float64) 379 self.A_5D_bc = zeros(self.num_tensors*5, float64) 380 381 # The rotation to the Frame Order eigenframe. 382 self.R_eigen = zeros((3, 3), float64) 383 self.R_ave = zeros((3, 3), float64) 384 self.Ri_prime = zeros((3, 3), float64) 385 self.tensor_3D = zeros((3, 3), float64) 386 387 # The cone axis storage and molecular frame z-axis. 388 self.cone_axis = zeros(3, float64) 389 self.z_axis = array([0, 0, 1], float64) 390 391 # Initialise the Frame Order matrices. 392 self.frame_order_2nd = zeros((9, 9), float64)
393 394
395 - def func_free_rotor(self, params):
396 """Target function for free rotor model optimisation. 397 398 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 399 400 401 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi}. 402 @type params: list of float 403 @return: The chi-squared or SSE value. 404 @rtype: float 405 """ 406 407 # Scaling. 408 if self.scaling_flag: 409 params = dot(params, self.scaling_matrix) 410 411 # Unpack the parameters. 412 if self.translation_opt and self.pivot_opt: 413 self._param_pivot = params[:3] 414 self._translation_vector = params[3:6] 415 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[6:] 416 elif self.translation_opt: 417 self._translation_vector = params[:3] 418 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 419 elif self.pivot_opt: 420 self._param_pivot = params[:3] 421 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 422 else: 423 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params 424 425 # Generate the cone axis from the spherical angles. 426 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 427 428 # Pre-calculate the eigenframe rotation matrix. 429 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 430 431 # The Kronecker product of the eigenframe rotation. 432 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 433 434 # Generate the 2nd degree Frame Order super matrix. 435 frame_order_2nd = compile_2nd_matrix_free_rotor(self.frame_order_2nd, Rx2_eigen) 436 437 # Reduce and rotate the tensors. 438 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 439 440 # Pre-transpose matrices for faster calculations. 441 RT_eigen = transpose(self.R_eigen) 442 RT_ave = transpose(self.R_ave) 443 444 # Pre-calculate all the necessary vectors. 445 if self.pcs_flag: 446 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 447 448 # Initial chi-squared (or SSE) value. 449 chi2_sum = 0.0 450 451 # Loop over each alignment. 452 for align_index in range(self.num_align): 453 # RDCs. 454 if self.rdc_flag[align_index]: 455 # Loop over the RDCs. 456 for j in range(self.num_interatom): 457 # The back calculated RDC. 458 if not self.missing_rdc[align_index, j]: 459 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 460 461 # Calculate and sum the single alignment chi-squared value (for the RDC). 462 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 463 464 # PCS. 465 if self.pcs_flag[align_index]: 466 # Loop over the PCSs. 467 for j in range(self.num_spins): 468 # The back calculated PCS. 469 if not self.missing_pcs[align_index, j]: 470 # Forwards and reverse rotations. 471 if self.full_in_ref_frame[align_index]: 472 r_pivot_atom = self.r_pivot_atom[:, j] 473 else: 474 r_pivot_atom = self.r_pivot_atom_rev[:, j] 475 476 # The numerical integration. 477 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) 478 479 # Calculate and sum the single alignment chi-squared value (for the PCS). 480 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 481 482 # Return the chi-squared value. 483 return chi2_sum
484 485
486 - def func_free_rotor_qrint(self, params):
487 """Target function for free rotor model optimisation. 488 489 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 490 491 492 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi}. 493 @type params: list of float 494 @return: The chi-squared or SSE value. 495 @rtype: float 496 """ 497 498 # Scaling. 499 if self.scaling_flag: 500 params = dot(params, self.scaling_matrix) 501 502 # Unpack the parameters. 503 if self.translation_opt and self.pivot_opt: 504 self._param_pivot = params[:3] 505 self._translation_vector = params[3:6] 506 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[6:] 507 elif self.translation_opt: 508 self._translation_vector = params[:3] 509 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 510 elif self.pivot_opt: 511 self._param_pivot = params[:3] 512 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params[3:] 513 else: 514 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = params 515 516 # Generate the cone axis from the spherical angles. 517 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 518 519 # Pre-calculate the eigenframe rotation matrix. 520 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 521 522 # The Kronecker product of the eigenframe rotation. 523 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 524 525 # Generate the 2nd degree Frame Order super matrix. 526 frame_order_2nd = compile_2nd_matrix_free_rotor(self.frame_order_2nd, Rx2_eigen) 527 528 # Reduce and rotate the tensors. 529 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 530 531 # Pre-transpose matrices for faster calculations. 532 RT_eigen = transpose(self.R_eigen) 533 RT_ave = transpose(self.R_ave) 534 535 # Pre-calculate all the necessary vectors. 536 if self.pcs_flag: 537 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 538 539 # Initial chi-squared (or SSE) value. 540 chi2_sum = 0.0 541 542 # Loop over each alignment. 543 for align_index in range(self.num_align): 544 # RDCs. 545 if self.rdc_flag[align_index]: 546 # Loop over the RDCs. 547 for j in range(self.num_interatom): 548 # The back calculated RDC. 549 if not self.missing_rdc[align_index, j]: 550 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 551 552 # Calculate and sum the single alignment chi-squared value (for the RDC). 553 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 554 555 # PCS via numerical integration. 556 if self.pcs_flag_sum: 557 # Numerical integration of the PCSs. 558 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) 559 560 # Calculate and sum the single alignment chi-squared value (for the PCS). 561 for align_index in range(self.num_align): 562 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 563 564 565 # Return the chi-squared value. 566 return chi2_sum
567 568
569 - def func_iso_cone(self, params):
570 """Target function for isotropic cone model optimisation. 571 572 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 573 574 575 @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. 576 @type params: list of float 577 @return: The chi-squared or SSE value. 578 @rtype: float 579 """ 580 581 # Scaling. 582 if self.scaling_flag: 583 params = dot(params, self.scaling_matrix) 584 585 # Unpack the parameters. 586 if self.translation_opt and self.pivot_opt: 587 self._param_pivot = params[:3] 588 self._translation_vector = params[3:6] 589 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[6:] 590 elif self.translation_opt: 591 self._translation_vector = params[:3] 592 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 593 elif self.pivot_opt: 594 self._param_pivot = params[:3] 595 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 596 else: 597 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params 598 599 # Generate the cone axis from the spherical angles. 600 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 601 602 # Pre-calculate the eigenframe rotation matrix. 603 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 604 605 # The Kronecker product of the eigenframe rotation. 606 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 607 608 # Generate the 2nd degree Frame Order super matrix. 609 frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd, Rx2_eigen, cone_theta, sigma_max) 610 611 # Reduce and rotate the tensors. 612 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 613 614 # Pre-transpose matrices for faster calculations. 615 RT_eigen = transpose(self.R_eigen) 616 RT_ave = transpose(self.R_ave) 617 618 # Pre-calculate all the necessary vectors. 619 if self.pcs_flag: 620 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 621 622 # Initial chi-squared (or SSE) value. 623 chi2_sum = 0.0 624 625 # Loop over each alignment. 626 for align_index in range(self.num_align): 627 # RDCs. 628 if self.rdc_flag[align_index]: 629 # Loop over the RDCs. 630 for j in range(self.num_interatom): 631 # The back calculated RDC. 632 if not self.missing_rdc[align_index, j]: 633 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 634 635 # Calculate and sum the single alignment chi-squared value (for the RDC). 636 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 637 638 # PCS. 639 if self.pcs_flag[align_index]: 640 # Loop over the PCSs. 641 for j in range(self.num_spins): 642 # The back calculated PCS. 643 if not self.missing_pcs[align_index, j]: 644 # Forwards and reverse rotations. 645 if self.full_in_ref_frame[align_index]: 646 r_pivot_atom = self.r_pivot_atom[:, j] 647 else: 648 r_pivot_atom = self.r_pivot_atom_rev[:, j] 649 650 # The numerical integration. 651 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) 652 653 # Calculate and sum the single alignment chi-squared value (for the PCS). 654 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 655 656 # Return the chi-squared value. 657 return chi2_sum
658 659
660 - def func_iso_cone_qrint(self, params):
661 """Target function for isotropic cone model optimisation. 662 663 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 664 665 666 @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. 667 @type params: list of float 668 @return: The chi-squared or SSE value. 669 @rtype: float 670 """ 671 672 # Scaling. 673 if self.scaling_flag: 674 params = dot(params, self.scaling_matrix) 675 676 # Unpack the parameters. 677 if self.translation_opt and self.pivot_opt: 678 self._param_pivot = params[:3] 679 self._translation_vector = params[3:6] 680 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[6:] 681 elif self.translation_opt: 682 self._translation_vector = params[:3] 683 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 684 elif self.pivot_opt: 685 self._param_pivot = params[:3] 686 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params[3:] 687 else: 688 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta, sigma_max = params 689 690 # Generate the cone axis from the spherical angles. 691 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 692 693 # Pre-calculate the eigenframe rotation matrix. 694 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 695 696 # The Kronecker product of the eigenframe rotation. 697 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 698 699 # Generate the 2nd degree Frame Order super matrix. 700 frame_order_2nd = compile_2nd_matrix_iso_cone(self.frame_order_2nd, Rx2_eigen, cone_theta, sigma_max) 701 702 # Reduce and rotate the tensors. 703 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 704 705 # Pre-transpose matrices for faster calculations. 706 RT_eigen = transpose(self.R_eigen) 707 RT_ave = transpose(self.R_ave) 708 709 # Pre-calculate all the necessary vectors. 710 if self.pcs_flag: 711 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 712 713 # Initial chi-squared (or SSE) value. 714 chi2_sum = 0.0 715 716 # Loop over each alignment. 717 for align_index in range(self.num_align): 718 # RDCs. 719 if self.rdc_flag[align_index]: 720 # Loop over the RDCs. 721 for j in range(self.num_interatom): 722 # The back calculated RDC. 723 if not self.missing_rdc[align_index, j]: 724 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 725 726 # Calculate and sum the single alignment chi-squared value (for the RDC). 727 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 728 729 # PCS via numerical integration. 730 if self.pcs_flag_sum: 731 # Numerical integration of the PCSs. 732 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) 733 734 # Calculate and sum the single alignment chi-squared value (for the PCS). 735 for align_index in range(self.num_align): 736 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 737 738 # Return the chi-squared value. 739 return chi2_sum
740 741
742 - def func_iso_cone_free_rotor(self, params):
743 """Target function for free rotor isotropic cone model optimisation. 744 745 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 746 747 748 @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. 749 @type params: list of float 750 @return: The chi-squared or SSE value. 751 @rtype: float 752 """ 753 754 # Scaling. 755 if self.scaling_flag: 756 params = dot(params, self.scaling_matrix) 757 758 # Unpack the parameters. 759 if self.translation_opt and self.pivot_opt: 760 self._param_pivot = params[:3] 761 self._translation_vector = params[3:6] 762 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[6:] 763 elif self.translation_opt: 764 self._translation_vector = params[:3] 765 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 766 elif self.pivot_opt: 767 self._param_pivot = params[:3] 768 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 769 else: 770 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params 771 772 # Generate the cone axis from the spherical angles. 773 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 774 775 # Pre-calculate the eigenframe rotation matrix. 776 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 777 778 # The Kronecker product of the eigenframe rotation. 779 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 780 781 # Calculate the cone angle. 782 theta_max = order_parameters.iso_cone_S_to_theta(cone_s1) 783 784 # Generate the 2nd degree Frame Order super matrix. 785 frame_order_2nd = compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_s1) 786 787 # Reduce and rotate the tensors. 788 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 789 790 # Pre-transpose matrices for faster calculations. 791 RT_eigen = transpose(self.R_eigen) 792 RT_ave = transpose(self.R_ave) 793 794 # Pre-calculate all the necessary vectors. 795 if self.pcs_flag: 796 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 797 798 # Initial chi-squared (or SSE) value. 799 chi2_sum = 0.0 800 801 # Loop over each alignment. 802 for align_index in range(self.num_align): 803 # RDCs. 804 if self.rdc_flag[align_index]: 805 # Loop over the RDCs. 806 for j in range(self.num_interatom): 807 # The back calculated RDC. 808 if not self.missing_rdc[align_index, j]: 809 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 810 811 # Calculate and sum the single alignment chi-squared value (for the RDC). 812 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 813 814 # PCS. 815 if self.pcs_flag[align_index]: 816 # Loop over the PCSs. 817 for j in range(self.num_spins): 818 # The back calculated PCS. 819 if not self.missing_pcs[align_index, j]: 820 # Forwards and reverse rotations. 821 if self.full_in_ref_frame[align_index]: 822 r_pivot_atom = self.r_pivot_atom[:, j] 823 else: 824 r_pivot_atom = self.r_pivot_atom_rev[:, j] 825 826 # The numerical integration. 827 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) 828 829 # Calculate and sum the single alignment chi-squared value (for the PCS). 830 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 831 832 # Return the chi-squared value. 833 return chi2_sum
834 835
836 - def func_iso_cone_free_rotor_qrint(self, params):
837 """Target function for free rotor isotropic cone model optimisation. 838 839 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 840 841 842 @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. 843 @type params: list of float 844 @return: The chi-squared or SSE value. 845 @rtype: float 846 """ 847 848 # Scaling. 849 if self.scaling_flag: 850 params = dot(params, self.scaling_matrix) 851 852 # Unpack the parameters. 853 if self.translation_opt and self.pivot_opt: 854 self._param_pivot = params[:3] 855 self._translation_vector = params[3:6] 856 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[6:] 857 elif self.translation_opt: 858 self._translation_vector = params[:3] 859 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 860 elif self.pivot_opt: 861 self._param_pivot = params[:3] 862 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params[3:] 863 else: 864 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = params 865 866 # Generate the cone axis from the spherical angles. 867 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 868 869 # Pre-calculate the eigenframe rotation matrix. 870 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 871 872 # The Kronecker product of the eigenframe rotation. 873 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 874 875 # Calculate the cone angle. 876 theta_max = order_parameters.iso_cone_S_to_theta(cone_s1) 877 878 # Generate the 2nd degree Frame Order super matrix. 879 frame_order_2nd = compile_2nd_matrix_iso_cone_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_s1) 880 881 # Reduce and rotate the tensors. 882 self.reduce_and_rot(0.0, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 883 884 # Pre-transpose matrices for faster calculations. 885 RT_eigen = transpose(self.R_eigen) 886 RT_ave = transpose(self.R_ave) 887 888 # Pre-calculate all the necessary vectors. 889 if self.pcs_flag: 890 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 891 892 # Initial chi-squared (or SSE) value. 893 chi2_sum = 0.0 894 895 # Loop over each alignment. 896 for align_index in range(self.num_align): 897 # RDCs. 898 if self.rdc_flag[align_index]: 899 # Loop over the RDCs. 900 for j in range(self.num_interatom): 901 # The back calculated RDC. 902 if not self.missing_rdc[align_index, j]: 903 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 904 905 # Calculate and sum the single alignment chi-squared value (for the RDC). 906 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 907 908 # PCS via numerical integration. 909 if self.pcs_flag_sum: 910 # Numerical integration of the PCSs. 911 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) 912 913 # Calculate and sum the single alignment chi-squared value (for the PCS). 914 for align_index in range(self.num_align): 915 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 916 917 # Return the chi-squared value. 918 return chi2_sum
919 920
921 - def func_iso_cone_torsionless(self, params):
922 """Target function for torsionless isotropic cone model optimisation. 923 924 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 925 926 927 @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. 928 @type params: list of float 929 @return: The chi-squared or SSE value. 930 @rtype: float 931 """ 932 933 # Scaling. 934 if self.scaling_flag: 935 params = dot(params, self.scaling_matrix) 936 937 # Unpack the parameters. 938 if self.translation_opt and self.pivot_opt: 939 self._param_pivot = params[:3] 940 self._translation_vector = params[3:6] 941 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[6:] 942 elif self.translation_opt: 943 self._translation_vector = params[:3] 944 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 945 elif self.pivot_opt: 946 self._param_pivot = params[:3] 947 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 948 else: 949 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params 950 951 # Generate the cone axis from the spherical angles. 952 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 953 954 # Pre-calculate the eigenframe rotation matrix. 955 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 956 957 # The Kronecker product of the eigenframe rotation. 958 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 959 960 # Generate the 2nd degree Frame Order super matrix. 961 frame_order_2nd = compile_2nd_matrix_iso_cone_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta) 962 963 # Reduce and rotate the tensors. 964 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 965 966 # Pre-transpose matrices for faster calculations. 967 RT_eigen = transpose(self.R_eigen) 968 RT_ave = transpose(self.R_ave) 969 970 # Pre-calculate all the necessary vectors. 971 if self.pcs_flag: 972 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 973 974 # Initial chi-squared (or SSE) value. 975 chi2_sum = 0.0 976 977 # Loop over each alignment. 978 for align_index in range(self.num_align): 979 # RDCs. 980 if self.rdc_flag[align_index]: 981 # Loop over the RDCs. 982 for j in range(self.num_interatom): 983 # The back calculated RDC. 984 if not self.missing_rdc[align_index, j]: 985 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 986 987 # Calculate and sum the single alignment chi-squared value (for the RDC). 988 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 989 990 # PCS. 991 if self.pcs_flag[align_index]: 992 # Loop over the PCSs. 993 for j in range(self.num_spins): 994 # The back calculated PCS. 995 if not self.missing_pcs[align_index, j]: 996 # Forwards and reverse rotations. 997 if self.full_in_ref_frame[align_index]: 998 r_pivot_atom = self.r_pivot_atom[:, j] 999 else: 1000 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1001 1002 # The numerical integration. 1003 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) 1004 1005 # Calculate and sum the single alignment chi-squared value (for the PCS). 1006 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1007 1008 # Return the chi-squared value. 1009 return chi2_sum
1010 1011
1012 - def func_iso_cone_torsionless_qrint(self, params):
1013 """Target function for torsionless isotropic cone model optimisation. 1014 1015 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1016 1017 1018 @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. 1019 @type params: list of float 1020 @return: The chi-squared or SSE value. 1021 @rtype: float 1022 """ 1023 1024 # Scaling. 1025 if self.scaling_flag: 1026 params = dot(params, self.scaling_matrix) 1027 1028 # Unpack the parameters. 1029 if self.translation_opt and self.pivot_opt: 1030 self._param_pivot = params[:3] 1031 self._translation_vector = params[3:6] 1032 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[6:] 1033 elif self.translation_opt: 1034 self._translation_vector = params[:3] 1035 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1036 elif self.pivot_opt: 1037 self._param_pivot = params[:3] 1038 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params[3:] 1039 else: 1040 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = params 1041 1042 # Generate the cone axis from the spherical angles. 1043 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1044 1045 # Pre-calculate the eigenframe rotation matrix. 1046 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1047 1048 # The Kronecker product of the eigenframe rotation. 1049 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1050 1051 # Generate the 2nd degree Frame Order super matrix. 1052 frame_order_2nd = compile_2nd_matrix_iso_cone_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta) 1053 1054 # Reduce and rotate the tensors. 1055 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1056 1057 # Pre-transpose matrices for faster calculations. 1058 RT_eigen = transpose(self.R_eigen) 1059 RT_ave = transpose(self.R_ave) 1060 1061 # Pre-calculate all the necessary vectors. 1062 if self.pcs_flag: 1063 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1064 1065 # Initial chi-squared (or SSE) value. 1066 chi2_sum = 0.0 1067 1068 # Loop over each alignment. 1069 for align_index in range(self.num_align): 1070 # RDCs. 1071 if self.rdc_flag[align_index]: 1072 # Loop over the RDCs. 1073 for j in range(self.num_interatom): 1074 # The back calculated RDC. 1075 if not self.missing_rdc[align_index, j]: 1076 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1077 1078 # Calculate and sum the single alignment chi-squared value (for the RDC). 1079 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1080 1081 # PCS via numerical integration. 1082 if self.pcs_flag_sum: 1083 # Numerical integration of the PCSs. 1084 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) 1085 1086 # Calculate and sum the single alignment chi-squared value (for the PCS). 1087 for align_index in range(self.num_align): 1088 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1089 1090 # Return the chi-squared value. 1091 return chi2_sum
1092 1093
1094 - def func_pseudo_ellipse(self, params):
1095 """Target function for pseudo-elliptic cone model optimisation. 1096 1097 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1098 1099 1100 @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. 1101 @type params: list of float 1102 @return: The chi-squared or SSE value. 1103 @rtype: float 1104 """ 1105 1106 # Scaling. 1107 if self.scaling_flag: 1108 params = dot(params, self.scaling_matrix) 1109 1110 # Unpack the parameters. 1111 if self.translation_opt and self.pivot_opt: 1112 self._param_pivot = params[:3] 1113 self._translation_vector = params[3:6] 1114 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:] 1115 elif self.translation_opt: 1116 self._translation_vector = params[:3] 1117 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:] 1118 elif self.pivot_opt: 1119 self._param_pivot = params[:3] 1120 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:] 1121 else: 1122 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 1123 1124 # Average position rotation. 1125 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1126 1127 # The Kronecker product of the eigenframe rotation. 1128 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1129 1130 # Generate the 2nd degree Frame Order super matrix. 1131 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y, cone_sigma_max) 1132 1133 # Reduce and rotate the tensors. 1134 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1135 1136 # Pre-transpose matrices for faster calculations. 1137 RT_eigen = transpose(self.R_eigen) 1138 RT_ave = transpose(self.R_ave) 1139 1140 # Pre-calculate all the necessary vectors. 1141 if self.pcs_flag: 1142 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1143 1144 # Initial chi-squared (or SSE) value. 1145 chi2_sum = 0.0 1146 1147 # Loop over each alignment. 1148 for align_index in range(self.num_align): 1149 # RDCs. 1150 if self.rdc_flag[align_index]: 1151 # Loop over the RDCs. 1152 for j in range(self.num_interatom): 1153 # The back calculated RDC. 1154 if not self.missing_rdc[align_index, j]: 1155 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1156 1157 # Calculate and sum the single alignment chi-squared value (for the RDC). 1158 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1159 1160 # PCS. 1161 if self.pcs_flag[align_index]: 1162 # Loop over the PCSs. 1163 for j in range(self.num_spins): 1164 # The back calculated PCS. 1165 if not self.missing_pcs[align_index, j]: 1166 # Forwards and reverse rotations. 1167 if self.full_in_ref_frame[align_index]: 1168 r_pivot_atom = self.r_pivot_atom[:, j] 1169 else: 1170 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1171 1172 # The numerical integration. 1173 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) 1174 1175 # Calculate and sum the single alignment chi-squared value (for the PCS). 1176 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1177 1178 # Return the chi-squared value. 1179 return chi2_sum
1180 1181
1182 - def func_pseudo_ellipse_qrint(self, params):
1183 """Target function for pseudo-elliptic cone model optimisation. 1184 1185 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. 1186 1187 1188 @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. 1189 @type params: list of float 1190 @return: The chi-squared or SSE value. 1191 @rtype: float 1192 """ 1193 1194 # Scaling. 1195 if self.scaling_flag: 1196 params = dot(params, self.scaling_matrix) 1197 1198 # Unpack the parameters. 1199 if self.translation_opt and self.pivot_opt: 1200 self._param_pivot = params[:3] 1201 self._translation_vector = params[3:6] 1202 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:] 1203 elif self.translation_opt: 1204 self._translation_vector = params[:3] 1205 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:] 1206 elif self.pivot_opt: 1207 self._param_pivot = params[:3] 1208 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:] 1209 else: 1210 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 1211 1212 # Average position rotation. 1213 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1214 1215 # The Kronecker product of the eigenframe rotation. 1216 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1217 1218 # Generate the 2nd degree Frame Order super matrix. 1219 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y, cone_sigma_max) 1220 1221 # Reduce and rotate the tensors. 1222 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1223 1224 # Pre-transpose matrices for faster calculations. 1225 RT_eigen = transpose(self.R_eigen) 1226 RT_ave = transpose(self.R_ave) 1227 1228 # Pre-calculate all the necessary vectors. 1229 if self.pcs_flag: 1230 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1231 1232 # Initial chi-squared (or SSE) value. 1233 chi2_sum = 0.0 1234 1235 # Loop over each alignment. 1236 for align_index in range(self.num_align): 1237 # RDCs. 1238 if self.rdc_flag[align_index]: 1239 # Loop over the RDCs. 1240 for j in range(self.num_interatom): 1241 # The back calculated RDC. 1242 if not self.missing_rdc[align_index, j]: 1243 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1244 1245 # Calculate and sum the single alignment chi-squared value (for the RDC). 1246 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1247 1248 # PCS via numerical integration. 1249 if self.pcs_flag_sum: 1250 # Numerical integration of the PCSs. 1251 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) 1252 1253 # Calculate and sum the single alignment chi-squared value (for the PCS). 1254 for align_index in range(self.num_align): 1255 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1256 1257 # Return the chi-squared value. 1258 return chi2_sum
1259 1260
1261 - def func_pseudo_ellipse_free_rotor(self, params):
1262 """Target function for free_rotor pseudo-elliptic cone model optimisation. 1263 1264 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1265 1266 1267 @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. 1268 @type params: list of float 1269 @return: The chi-squared or SSE value. 1270 @rtype: float 1271 """ 1272 1273 # Scaling. 1274 if self.scaling_flag: 1275 params = dot(params, self.scaling_matrix) 1276 1277 # Unpack the parameters. 1278 if self.translation_opt and self.pivot_opt: 1279 self._param_pivot = params[:3] 1280 self._translation_vector = params[3:6] 1281 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1282 elif self.translation_opt: 1283 self._translation_vector = params[:3] 1284 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1285 elif self.pivot_opt: 1286 self._param_pivot = params[:3] 1287 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1288 else: 1289 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1290 1291 # Average position rotation. 1292 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1293 1294 # The Kronecker product of the eigenframe rotation. 1295 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1296 1297 # Generate the 2nd degree Frame Order super matrix. 1298 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1299 1300 # Reduce and rotate the tensors. 1301 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1302 1303 # Pre-transpose matrices for faster calculations. 1304 RT_eigen = transpose(self.R_eigen) 1305 RT_ave = transpose(self.R_ave) 1306 1307 # Pre-calculate all the necessary vectors. 1308 if self.pcs_flag: 1309 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1310 1311 # Initial chi-squared (or SSE) value. 1312 chi2_sum = 0.0 1313 1314 # Loop over each alignment. 1315 for align_index in range(self.num_align): 1316 # RDCs. 1317 if self.rdc_flag[align_index]: 1318 # Loop over the RDCs. 1319 for j in range(self.num_interatom): 1320 # The back calculated RDC. 1321 if not self.missing_rdc[align_index, j]: 1322 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1323 1324 # Calculate and sum the single alignment chi-squared value (for the RDC). 1325 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1326 1327 # PCS. 1328 if self.pcs_flag[align_index]: 1329 # Loop over the PCSs. 1330 for j in range(self.num_spins): 1331 # The back calculated PCS. 1332 if not self.missing_pcs[align_index, j]: 1333 # Forwards and reverse rotations. 1334 if self.full_in_ref_frame[align_index]: 1335 r_pivot_atom = self.r_pivot_atom[:, j] 1336 else: 1337 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1338 1339 # The numerical integration. 1340 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) 1341 1342 # Calculate and sum the single alignment chi-squared value (for the PCS). 1343 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1344 1345 # Return the chi-squared value. 1346 return chi2_sum
1347 1348
1349 - def func_pseudo_ellipse_free_rotor_qrint(self, params):
1350 """Target function for free_rotor pseudo-elliptic cone model optimisation. 1351 1352 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1353 1354 1355 @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. 1356 @type params: list of float 1357 @return: The chi-squared or SSE value. 1358 @rtype: float 1359 """ 1360 1361 # Scaling. 1362 if self.scaling_flag: 1363 params = dot(params, self.scaling_matrix) 1364 1365 # Unpack the parameters. 1366 if self.translation_opt and self.pivot_opt: 1367 self._param_pivot = params[:3] 1368 self._translation_vector = params[3:6] 1369 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1370 elif self.translation_opt: 1371 self._translation_vector = params[:3] 1372 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1373 elif self.pivot_opt: 1374 self._param_pivot = params[:3] 1375 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1376 else: 1377 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1378 1379 # Average position rotation. 1380 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1381 1382 # The Kronecker product of the eigenframe rotation. 1383 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1384 1385 # Generate the 2nd degree Frame Order super matrix. 1386 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1387 1388 # Reduce and rotate the tensors. 1389 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1390 1391 # Pre-transpose matrices for faster calculations. 1392 RT_eigen = transpose(self.R_eigen) 1393 RT_ave = transpose(self.R_ave) 1394 1395 # Pre-calculate all the necessary vectors. 1396 if self.pcs_flag: 1397 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1398 1399 # Initial chi-squared (or SSE) value. 1400 chi2_sum = 0.0 1401 1402 # Loop over each alignment. 1403 for align_index in range(self.num_align): 1404 # RDCs. 1405 if self.rdc_flag[align_index]: 1406 # Loop over the RDCs. 1407 for j in range(self.num_interatom): 1408 # The back calculated RDC. 1409 if not self.missing_rdc[align_index, j]: 1410 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1411 1412 # Calculate and sum the single alignment chi-squared value (for the RDC). 1413 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1414 1415 # PCS via numerical integration. 1416 if self.pcs_flag_sum: 1417 # Numerical integration of the PCSs. 1418 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) 1419 1420 # Calculate and sum the single alignment chi-squared value (for the PCS). 1421 for align_index in range(self.num_align): 1422 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1423 1424 # Return the chi-squared value. 1425 return chi2_sum
1426 1427
1428 - def func_pseudo_ellipse_torsionless(self, params):
1429 """Target function for torsionless pseudo-elliptic cone model optimisation. 1430 1431 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1432 1433 1434 @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. 1435 @type params: list of float 1436 @return: The chi-squared or SSE value. 1437 @rtype: float 1438 """ 1439 1440 # Scaling. 1441 if self.scaling_flag: 1442 params = dot(params, self.scaling_matrix) 1443 1444 # Unpack the parameters. 1445 if self.translation_opt and self.pivot_opt: 1446 self._param_pivot = params[:3] 1447 self._translation_vector = params[3:6] 1448 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1449 elif self.translation_opt: 1450 self._translation_vector = params[:3] 1451 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1452 elif self.pivot_opt: 1453 self._param_pivot = params[:3] 1454 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1455 else: 1456 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1457 1458 # Average position rotation. 1459 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1460 1461 # The Kronecker product of the eigenframe rotation. 1462 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1463 1464 # Generate the 2nd degree Frame Order super matrix. 1465 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1466 1467 # Reduce and rotate the tensors. 1468 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1469 1470 # Pre-transpose matrices for faster calculations. 1471 RT_eigen = transpose(self.R_eigen) 1472 RT_ave = transpose(self.R_ave) 1473 1474 # Pre-calculate all the necessary vectors. 1475 if self.pcs_flag: 1476 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1477 1478 # Initial chi-squared (or SSE) value. 1479 chi2_sum = 0.0 1480 1481 # Loop over each alignment. 1482 for align_index in range(self.num_align): 1483 # RDCs. 1484 if self.rdc_flag[align_index]: 1485 # Loop over the RDCs. 1486 for j in range(self.num_interatom): 1487 # The back calculated RDC. 1488 if not self.missing_rdc[align_index, j]: 1489 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1490 1491 # Calculate and sum the single alignment chi-squared value (for the RDC). 1492 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1493 1494 # PCS. 1495 if self.pcs_flag[align_index]: 1496 # Loop over the PCSs. 1497 for j in range(self.num_spins): 1498 # The back calculated PCS. 1499 if not self.missing_pcs[align_index, j]: 1500 # Forwards and reverse rotations. 1501 if self.full_in_ref_frame[align_index]: 1502 r_pivot_atom = self.r_pivot_atom[:, j] 1503 else: 1504 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1505 1506 # The numerical integration. 1507 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) 1508 1509 # Calculate and sum the single alignment chi-squared value (for the PCS). 1510 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1511 1512 # Return the chi-squared value. 1513 return chi2_sum
1514 1515
1516 - def func_pseudo_ellipse_torsionless_qrint(self, params):
1517 """Target function for torsionless pseudo-elliptic cone model optimisation. 1518 1519 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Simple numerical integration is used for the PCS. 1520 1521 1522 @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. 1523 @type params: list of float 1524 @return: The chi-squared or SSE value. 1525 @rtype: float 1526 """ 1527 1528 # Scaling. 1529 if self.scaling_flag: 1530 params = dot(params, self.scaling_matrix) 1531 1532 # Unpack the parameters. 1533 if self.translation_opt and self.pivot_opt: 1534 self._param_pivot = params[:3] 1535 self._translation_vector = params[3:6] 1536 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[6:] 1537 elif self.translation_opt: 1538 self._translation_vector = params[:3] 1539 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1540 elif self.pivot_opt: 1541 self._param_pivot = params[:3] 1542 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params[3:] 1543 else: 1544 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = params 1545 1546 # Average position rotation. 1547 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen) 1548 1549 # The Kronecker product of the eigenframe rotation. 1550 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1551 1552 # Generate the 2nd degree Frame Order super matrix. 1553 frame_order_2nd = compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, Rx2_eigen, cone_theta_x, cone_theta_y) 1554 1555 # Reduce and rotate the tensors. 1556 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1557 1558 # Pre-transpose matrices for faster calculations. 1559 RT_eigen = transpose(self.R_eigen) 1560 RT_ave = transpose(self.R_ave) 1561 1562 # Pre-calculate all the necessary vectors. 1563 if self.pcs_flag: 1564 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1565 1566 # Initial chi-squared (or SSE) value. 1567 chi2_sum = 0.0 1568 1569 # Loop over each alignment. 1570 for align_index in range(self.num_align): 1571 # RDCs. 1572 if self.rdc_flag[align_index]: 1573 # Loop over the RDCs. 1574 for j in range(self.num_interatom): 1575 # The back calculated RDC. 1576 if not self.missing_rdc[align_index, j]: 1577 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1578 1579 # Calculate and sum the single alignment chi-squared value (for the RDC). 1580 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1581 1582 # PCS via numerical integration. 1583 if self.pcs_flag_sum: 1584 # Numerical integration of the PCSs. 1585 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) 1586 1587 # Calculate and sum the single alignment chi-squared value (for the PCS). 1588 for align_index in range(self.num_align): 1589 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1590 1591 # Return the chi-squared value. 1592 return chi2_sum
1593 1594
1595 - def func_rigid(self, params):
1596 """Target function for rigid model optimisation. 1597 1598 This function optimises the isotropic cone model parameters using the RDC and PCS base data. 1599 1600 1601 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma}. 1602 @type params: list of float 1603 @return: The chi-squared or SSE value. 1604 @rtype: float 1605 """ 1606 1607 # Scaling. 1608 if self.scaling_flag: 1609 params = dot(params, self.scaling_matrix) 1610 1611 # Unpack the parameters. 1612 if self.translation_opt: 1613 self._translation_vector = params[:3] 1614 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = params[3:6] 1615 else: 1616 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = params 1617 1618 # The average frame rotation matrix (and reduce and rotate the tensors). 1619 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma) 1620 1621 # Pre-transpose matrices for faster calculations. 1622 RT_ave = transpose(self.R_ave) 1623 1624 # Pre-calculate all the necessary vectors. 1625 if self.pcs_flag: 1626 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1627 1628 # Initial chi-squared (or SSE) value. 1629 chi2_sum = 0.0 1630 1631 # Loop over each alignment. 1632 for align_index in range(self.num_align): 1633 # RDCs. 1634 if self.rdc_flag[align_index]: 1635 # Loop over the RDCs. 1636 for j in range(self.num_interatom): 1637 # The back calculated RDC. 1638 if not self.missing_rdc[align_index, j]: 1639 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1640 1641 # Calculate and sum the single alignment chi-squared value (for the RDC). 1642 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1643 1644 # PCS. 1645 if self.pcs_flag[align_index]: 1646 # Loop over the PCSs. 1647 for j in range(self.num_spins): 1648 # The back calculated PCS. 1649 if not self.missing_pcs[align_index, j]: 1650 # Forwards and reverse rotations. 1651 if self.full_in_ref_frame[align_index]: 1652 r_pivot_atom = self.r_pivot_atom[:, j] 1653 else: 1654 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1655 1656 # The PCS calculation. 1657 vect = self.r_ln_pivot[:, 0] + r_pivot_atom 1658 length = norm(vect) 1659 self.pcs_theta[align_index, j] = pcs_tensor(self.pcs_const[align_index] / length**5, vect, self.A_3D[align_index]) 1660 1661 # Calculate and sum the single alignment chi-squared value (for the PCS). 1662 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1663 1664 # Return the chi-squared value. 1665 return chi2_sum
1666 1667
1668 - def func_rotor(self, params):
1669 """Target function for rotor model optimisation. 1670 1671 This function optimises the isotropic cone model parameters using the RDC and PCS base data. Scipy quadratic integration is used for the PCS. 1672 1673 1674 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 1675 @type params: list of float 1676 @return: The chi-squared or SSE value. 1677 @rtype: float 1678 """ 1679 1680 # Scaling. 1681 if self.scaling_flag: 1682 params = dot(params, self.scaling_matrix) 1683 1684 # Unpack the parameters. 1685 if self.translation_opt and self.pivot_opt: 1686 self._param_pivot = params[:3] 1687 self._translation_vector = params[3:6] 1688 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[6:] 1689 elif self.translation_opt: 1690 self._translation_vector = params[:3] 1691 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[3:] 1692 elif self.pivot_opt: 1693 self._param_pivot = params[:3] 1694 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[3:] 1695 else: 1696 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params 1697 1698 # Generate the cone axis from the spherical angles. 1699 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1700 1701 # Pre-calculate the eigenframe rotation matrix. 1702 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1703 1704 # The Kronecker product of the eigenframe rotation. 1705 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1706 1707 # Generate the 2nd degree Frame Order super matrix. 1708 frame_order_2nd = compile_2nd_matrix_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max) 1709 1710 # The average frame rotation matrix (and reduce and rotate the tensors). 1711 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1712 1713 # Pre-transpose matrices for faster calculations. 1714 RT_eigen = transpose(self.R_eigen) 1715 RT_ave = transpose(self.R_ave) 1716 1717 # Pre-calculate all the necessary vectors. 1718 if self.pcs_flag: 1719 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1720 1721 # Initial chi-squared (or SSE) value. 1722 chi2_sum = 0.0 1723 1724 # Loop over each alignment. 1725 for align_index in range(self.num_align): 1726 # RDCs. 1727 if self.rdc_flag[align_index]: 1728 # Loop over the RDCs. 1729 for j in range(self.num_interatom): 1730 # The back calculated RDC. 1731 if not self.missing_rdc[align_index, j]: 1732 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1733 1734 # Calculate and sum the single alignment chi-squared value (for the RDC). 1735 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1736 1737 # PCS. 1738 if self.pcs_flag[align_index]: 1739 # Loop over the PCSs. 1740 for j in range(self.num_spins): 1741 # The back calculated PCS. 1742 if not self.missing_pcs[align_index, j]: 1743 # Forwards and reverse rotations. 1744 if self.full_in_ref_frame[align_index]: 1745 r_pivot_atom = self.r_pivot_atom[:, j] 1746 else: 1747 r_pivot_atom = self.r_pivot_atom_rev[:, j] 1748 1749 # The numerical integration. 1750 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) 1751 1752 # Calculate and sum the single alignment chi-squared value (for the PCS). 1753 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1754 1755 # Return the chi-squared value. 1756 return chi2_sum
1757 1758
1759 - def func_rotor_qrint(self, params):
1760 """Target function for rotor model optimisation. 1761 1762 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. 1763 1764 1765 @param params: The vector of parameter values. These are the tensor rotation angles {alpha, beta, gamma, theta, phi, sigma_max}. 1766 @type params: list of float 1767 @return: The chi-squared or SSE value. 1768 @rtype: float 1769 """ 1770 1771 # Scaling. 1772 if self.scaling_flag: 1773 params = dot(params, self.scaling_matrix) 1774 1775 # Unpack the parameters. 1776 if self.translation_opt and self.pivot_opt: 1777 self._param_pivot = params[:3] 1778 self._translation_vector = params[3:6] 1779 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[6:] 1780 elif self.translation_opt: 1781 self._translation_vector = params[:3] 1782 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[3:] 1783 elif self.pivot_opt: 1784 self._param_pivot = params[:3] 1785 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params[3:] 1786 else: 1787 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, sigma_max = params 1788 1789 # Generate the cone axis from the spherical angles. 1790 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis) 1791 1792 # Pre-calculate the eigenframe rotation matrix. 1793 two_vect_to_R(self.z_axis, self.cone_axis, self.R_eigen) 1794 1795 # The Kronecker product of the eigenframe rotation. 1796 Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen) 1797 1798 # Generate the 2nd degree Frame Order super matrix. 1799 frame_order_2nd = compile_2nd_matrix_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max) 1800 1801 # The average frame rotation matrix (and reduce and rotate the tensors). 1802 self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, frame_order_2nd) 1803 1804 # Pre-transpose matrices for faster calculations. 1805 RT_eigen = transpose(self.R_eigen) 1806 RT_ave = transpose(self.R_ave) 1807 1808 # Pre-calculate all the necessary vectors. 1809 if self.pcs_flag: 1810 self.calc_vectors(self._param_pivot, self.R_ave, RT_ave) 1811 1812 # Initial chi-squared (or SSE) value. 1813 chi2_sum = 0.0 1814 1815 # Loop over each alignment. 1816 for align_index in range(self.num_align): 1817 # RDCs. 1818 if self.rdc_flag[align_index]: 1819 # Loop over the RDCs. 1820 for j in range(self.num_interatom): 1821 # The back calculated RDC. 1822 if not self.missing_rdc[align_index, j]: 1823 self.rdc_theta[align_index, j] = rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index]) 1824 1825 # Calculate and sum the single alignment chi-squared value (for the RDC). 1826 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_error[align_index]) 1827 1828 # PCS via numerical integration. 1829 if self.pcs_flag_sum: 1830 # Numerical integration of the PCSs. 1831 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) 1832 1833 # Calculate and sum the single alignment chi-squared value (for the PCS). 1834 for align_index in range(self.num_align): 1835 chi2_sum = chi2_sum + chi2(self.pcs[align_index], self.pcs_theta[align_index], self.pcs_error[align_index]) 1836 1837 # Return the chi-squared value. 1838 return chi2_sum
1839 1840
1841 - def calc_vectors(self, pivot=None, R_ave=None, RT_ave=None):
1842 """Calculate the pivot to atom and lanthanide to pivot vectors for the target functions. 1843 1844 @keyword pivot: The pivot point. 1845 @type pivot: numpy rank-1, 3D array 1846 @keyword R_ave: The rotation matrix for rotating from the reference frame to the average position. 1847 @type R_ave: numpy rank-2, 3D array 1848 @keyword RT_ave: The transpose of R_ave. 1849 @type RT_ave: numpy rank-2, 3D array 1850 """ 1851 1852 # The rotational pivot. 1853 if self.ave_pos_piv_sync: 1854 ave_pos_pivot = self._param_pivot 1855 else: 1856 ave_pos_pivot = self.ave_pos_pivot 1857 1858 # The pivot to atom vectors. 1859 for j in range(self.num_spins): 1860 # The lanthanide to pivot vector. 1861 if self.pivot_opt: 1862 self.r_ln_pivot[:, j] = pivot - self.paramag_centre 1863 1864 # Rotate then translate the atomic positions, then calculate the pivot to atom vector. 1865 self.r_pivot_atom[:, j] = dot(R_ave, self.atomic_pos[j] - ave_pos_pivot) + ave_pos_pivot 1866 self.r_pivot_atom[:, j] += self._translation_vector 1867 self.r_pivot_atom[:, j] -= pivot 1868 self.r_pivot_atom_rev[:, j] = dot(RT_ave, self.atomic_pos[j] - ave_pos_pivot) + ave_pos_pivot 1869 self.r_pivot_atom_rev[:, j] += self._translation_vector 1870 self.r_pivot_atom_rev[:, j] -= pivot
1871 1872
1873 - def create_sobol_data(self, n=10000, dims=None):
1874 """Create the Sobol' quasi-random data for numerical integration. 1875 1876 This uses the external sobol_lib module to create the data. The algorithm is that modified by Antonov and Saleev. 1877 1878 1879 @keyword n: The number of points to generate. 1880 @type n: int 1881 @keyword dims: The list of parameters. 1882 @type dims: list of str 1883 """ 1884 1885 # The number of dimensions. 1886 m = len(dims) 1887 1888 # Initialise. 1889 self.sobol_angles = zeros((n, m), float32) 1890 1891 # Loop over the points. 1892 for i in range(n): 1893 # The raw point. 1894 point, seed = i4_sobol(m, i) 1895 1896 # Loop over the dimensions, converting the points to angles. 1897 for j in range(m): 1898 # The tilt angle - the angle of rotation about the x-y plane rotation axis. 1899 if dims[j] in ['theta']: 1900 self.sobol_angles[i, j] = acos(2.0*point[j] - 1.0) 1901 1902 # The angle defining the x-y plane rotation axis. 1903 if dims[j] in ['phi']: 1904 self.sobol_angles[i, j] = 2.0 * pi * point[j] 1905 1906 # The torsion angle - the angle of rotation about the z' axis. 1907 if dims[j] in ['sigma']: 1908 self.sobol_angles[i, j] = 2.0 * pi * (point[j] - 0.5)
1909 1910
1911 - def reduce_and_rot(self, ave_pos_alpha=None, ave_pos_beta=None, ave_pos_gamma=None, daeg=None):
1912 """Reduce and rotate the alignments tensors using the frame order matrix and Euler angles. 1913 1914 @keyword ave_pos_alpha: The alpha Euler angle describing the average domain position, the tensor rotation. 1915 @type ave_pos_alpha: float 1916 @keyword ave_pos_beta: The beta Euler angle describing the average domain position, the tensor rotation. 1917 @type ave_pos_beta: float 1918 @keyword ave_pos_gamma: The gamma Euler angle describing the average domain position, the tensor rotation. 1919 @type ave_pos_gamma: float 1920 @keyword daeg: The 2nd degree frame order matrix. 1921 @type daeg: rank-2, 9D array 1922 """ 1923 1924 # Alignment tensor rotation. 1925 euler_to_R_zyz(ave_pos_alpha, ave_pos_beta, ave_pos_gamma, self.R_ave) 1926 1927 # Back calculate the rotated tensors. 1928 for align_index in range(self.num_tensors): 1929 # Tensor indices. 1930 index1 = align_index*5 1931 index2 = align_index*5+5 1932 1933 # Reduction. 1934 if daeg != None: 1935 # Reduce the tensor. 1936 reduce_alignment_tensor(daeg, self.full_tensors[index1:index2], self.A_5D_bc[index1:index2]) 1937 1938 # Convert the reduced tensor to 3D, rank-2 form. 1939 to_tensor(self.tensor_3D, self.A_5D_bc[index1:index2]) 1940 1941 # No reduction: 1942 else: 1943 # Convert the original tensor to 3D, rank-2 form. 1944 to_tensor(self.tensor_3D, self.full_tensors[index1:index2]) 1945 1946 # Rotate the tensor (normal R.X.RT rotation). 1947 if self.full_in_ref_frame[align_index]: 1948 self.A_3D_bc[align_index] = dot(transpose(self.R_ave), dot(self.tensor_3D, self.R_ave)) 1949 1950 # Rotate the tensor (inverse RT.X.R rotation). 1951 else: 1952 self.A_3D_bc[align_index] = dot(self.R_ave, dot(self.tensor_3D, transpose(self.R_ave))) 1953 1954 # Convert the tensor back to 5D, rank-1 form, as the back-calculated reduced tensor. 1955 to_5D(self.A_5D_bc[index1:index2], self.A_3D_bc[align_index])
1956