maths_fns.n_state

39 """Class containing the target function of the optimisation of the N-state model.""" 40

41 - def __init__(self, model=None, N=None, init_params=None, probs=None, full_tensors=None, red_data=None, red_errors=None, full_in_ref_frame=None, fixed_tensors=None, pcs=None, pcs_errors=None, pcs_weights=None, rdcs=None, rdc_errors=None, rdc_weights=None, rdc_vect=None, temp=None, frq=None, dip_const=None, absolute_rdc=None, atomic_pos=None, paramag_centre=None, scaling_matrix=None, centre_fixed=True):

42 """Set up the class instance for optimisation. 43 44 The N-state models 45 ================== 46 47 All constant data required for the N-state model are initialised here. Depending on the base data used for optimisation, different data structures can be supplied. However a number of structures must be provided for the N-state model. These are: 48 49 - model, the type of N-state model. This can be '2-domain', 'population', or 'fixed'. 50 - N, the number of states (or structures). 51 - init_params, the initial parameter values. 52 - scaling_matrix, the matrix used for parameter scaling during optimisation. 53 54 55 2-domain N-state model 56 ---------------------- 57 58 If the model type is set to '2-domain', then the following data structures should be supplied: 59 60 - full_tensors, the alignment tensors in matrix form. 61 - red_data, the alignment tensors in 5D form in a rank-1 array. 62 - red_errors, the alignment tensor errors in 5D form in a rank-1 array. This data is not obligatory. 63 - full_in_ref_frame, an array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 64 65 66 The population N-state model 67 ============================ 68 69 In this model, populations are optimised for each state. Additionally the alignment tensors for anisotropic data can also be optimised if they have not been supplied (through the full_tensors arg). 70 71 72 PCS base data 73 ------------- 74 75 If pseudocontact shift data is to be used for optimisation, then the following should be supplied: 76 77 - pcs, the pseudocontact shifts. 78 - pcs_errors, the optional pseudocontact shift error values. 79 - temp, the temperatures of the experiments. 80 - frq, the frequencies of the experiments. 81 82 83 PCS and PRE base data 84 --------------------- 85 86 If either pseudocontact shift or PRE data is to be used for optimisation, then the following should be supplied: 87 88 - atomic_pos, the positions of all atoms. 89 - paramag_centre, the paramagnetic centre position. 90 91 92 RDC base data 93 ------------- 94 95 If residual dipolar coupling data is to be used for optimisation, then the following should be supplied: 96 97 - rdcs, the residual dipolar couplings. 98 - rdc_errors, the optional residual dipolar coupling errors. 99 - rdc_vect, the interatomic vectors. 100 - dip_const, the dipolar constants. 101 - absolute_rdc, the flags specifying whether each RDC is signless. 102 103 104 @keyword model: The N-state model type. This can be one of '2-domain', 'population' or 'fixed'. 105 @type model: str 106 @keyword N: The number of states. 107 @type N: int 108 @keyword init_params: The initial parameter values. Optimisation must start at some point! 109 @type init_params: numpy float64 array 110 @keyword probs: The probabilities for each state c. The length of this list should be equal to N. 111 @type probs: list of float 112 @keyword full_tensors: An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all full tensors. The format is [Sxx1, Syy1, Sxy1, Sxz1, Syz1, Sxx2, Syy2, Sxy2, Sxz2, Syz2, ..., Sxxn, Syyn, Sxyn, Sxzn, Syzn] 113 @type full_tensors: list of rank-2, 3D numpy arrays 114 @keyword red_data: An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all reduced tensors. The format is the same as for full_tensors. 115 @type red_data: numpy float64 array 116 @keyword red_errors: An array of the {Sxx, Syy, Sxy, Sxz, Syz} errors for all reduced tensors. The array format is the same as for full_tensors. 117 @type red_errors: numpy float64 array 118 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 119 @type full_in_ref_frame: numpy rank-1 array 120 @keyword fixed_tensors: An array of flags specifying if the tensor is fixed or will be optimised. 121 @type fixed_tensors: list of bool 122 @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. 123 @type pcs: numpy rank-2 array 124 @keyword pcs_errors: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 125 @type pcs_errors: numpy rank-2 array 126 @keyword pcs_weights: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 127 @type pcs_weights: numpy rank-2 array 128 @keyword rdcs: The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin pairs k. 129 @type rdcs: numpy rank-2 array 130 @keyword rdc_errors: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 131 @type rdc_errors: numpy rank-2 array 132 @keyword rdc_weights: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 133 @type rdc_weights: numpy rank-2 array 134 @keyword rdc_vect: The unit vector lists for the RDC. The first index must correspond to the spin pair and the second index to each structure (its size being equal to the number of states). 135 @type rdc_vect: numpy rank-2 array 136 @keyword temp: The temperature of each experiment, used for the PCS. 137 @type temp: numpy rank-1 array 138 @keyword frq: The frequency of each alignment, used for the PCS. 139 @type frq: numpy rank-1 array 140 @keyword dip_const: The dipolar constants for each spin pair. The indices correspond to the spin pairs k. 141 @type dip_const: numpy rank-1 array 142 @keyword absolute_rdc: The signless or absolute value flags for the RDCs. 143 @type absolute_rdc: numpy rank-2 array 144 @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. 145 @type atomic_pos: numpy rank-3 array 146 @keyword paramag_centre: The paramagnetic centre position (or positions). 147 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 148 @keyword scaling_matrix: The square and diagonal scaling matrix. 149 @type scaling_matrix: numpy rank-2 array 150 @keyword centre_fixed: A flag which if False will cause the paramagnetic centre to be optimised. 151 @type centre_fixed: bool 152 """ 153 154 # Store the data inside the class instance namespace. 155 self.N = N 156 self.params = 1.0 * init_params # Force a copy of the data to be stored. 157 self.fixed_tensors = fixed_tensors 158 self.deltaij = pcs 159 self.Dij = rdcs 160 self.dip_vect = rdc_vect 161 self.dip_const = dip_const 162 self.absolute_rdc = absolute_rdc 163 self.temp = temp 164 self.frq = frq 165 self.atomic_pos = atomic_pos 166 self.paramag_centre = paramag_centre 167 self.centre_fixed = centre_fixed 168 self.total_num_params = len(init_params) 169 170 # Initialise the function value, gradient, and Hessian. 171 self.chi2 = 0.0 172 self.dchi2 = zeros((self.total_num_params), float64) 173 self.d2chi2 = zeros((self.total_num_params, self.total_num_params), float64) 174 175 # Scaling initialisation. 176 self.scaling_matrix = scaling_matrix 177 if self.scaling_matrix != None: 178 self.scaling_flag = True 179 else: 180 self.scaling_flag = False 181 182 # The 2-domain N-state model. 183 if model == '2-domain': 184 # Some checks. 185 if red_data == None or not len(red_data): 186 raise RelaxError("The red_data argument " + repr(red_data) + " must be supplied.") 187 if red_errors == None or not len(red_errors): 188 raise RelaxError("The red_errors argument " + repr(red_errors) + " must be supplied.") 189 if full_in_ref_frame == None or not len(full_in_ref_frame): 190 raise RelaxError("The full_in_ref_frame argument " + repr(full_in_ref_frame) + " must be supplied.") 191 192 # Tensor set up. 193 self.full_tensors = array(full_tensors, float64) 194 self.num_tensors = int(len(self.full_tensors) / 5) 195 self.red_data = red_data 196 self.red_errors = red_errors 197 self.full_in_ref_frame = full_in_ref_frame 198 199 # Alignment tensor in rank-2, 3D form. 200 self.A = zeros((self.num_tensors, 3, 3), float64) 201 for align_index in range(self.num_tensors): 202 to_tensor(self.A[align_index], self.full_tensors[5*align_index:5*align_index+5]) 203 204 # Initialise some empty numpy objects for storage of: 205 # R: the transient rotation matricies. 206 # RT: the transposes of the rotation matricies. 207 # red_bc: the back-calculated reduced alignment tensors. 208 # red_bc_vector: the back-calculated reduced alignment tensors in vector form {Sxx, Syy, Sxy, Sxz, Syz}. 209 self.R = zeros((self.N, 3, 3), float64) 210 self.RT = zeros((self.N, 3, 3), float64) 211 self.red_bc = zeros((self.num_tensors, 3, 3), float64) 212 self.red_bc_vector = zeros(self.num_tensors*5, float64) 213 214 # Set the target function. 215 self.func = self.func_2domain 216 self.dfunc = None 217 self.d2func = None 218 219 # The flexible population or equal probability N-state models. 220 elif model in ['population', 'fixed']: 221 # The total number of alignments (must be the same for the RDC and PCS data). 222 self.num_align = 0 223 if rdcs != None: 224 self.num_align = len(rdcs) 225 elif pcs != None: 226 self.num_align = len(pcs) 227 228 # Set the RDC and PCS flags (indicating the presence of data). 229 self.rdc_flag = [True] * self.num_align 230 self.pcs_flag = [True] * self.num_align 231 for align_index in range(self.num_align): 232 if rdcs == None or len(rdcs[align_index]) == 0: 233 self.rdc_flag[align_index] = False 234 if pcs == None or len(pcs[align_index]) == 0: 235 self.pcs_flag[align_index] = False 236 self.rdc_flag_sum = sum(self.rdc_flag) 237 self.pcs_flag_sum = sum(self.pcs_flag) 238 239 # Some checks. 240 if self.rdc_flag_sum and (rdc_vect == None or not len(rdc_vect)): 241 raise RelaxError("The rdc_vect argument " + repr(rdc_vect) + " must be supplied.") 242 if self.pcs_flag_sum and (atomic_pos == None or not len(atomic_pos)): 243 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.") 244 245 # The total number of spins. 246 self.num_spins = 0 247 if self.pcs_flag_sum: 248 self.num_spins = len(pcs[0]) 249 250 # The total number of interatomic connections. 251 self.num_interatom = 0 252 if self.rdc_flag_sum: 253 self.num_interatom = len(rdcs[0]) 254 255 # Alignment tensor function and gradient matrices. 256 self.A = zeros((self.num_align, 3, 3), float64) 257 self.dA = zeros((5, 3, 3), float64) 258 259 # Fixed alignment tensors. 260 if full_tensors != None: 261 # Convert to numpy. 262 self.full_tensors = array(full_tensors, float64) 263 264 # Set up the alignment data. 265 self.num_align_params = 0 266 index = 0 267 for align_index in range(self.num_align): 268 # Fill the alignment tensor object with the fixed tensors. 269 if fixed_tensors[align_index]: 270 to_tensor(self.A[align_index], self.full_tensors[5*index:5*index+5]) 271 index += 1 272 273 # The number of alignment parameters. 274 if not fixed_tensors[align_index]: 275 self.num_align_params += 5 276 277 # Optimised alignment tensors. 278 else: 279 # The alignment tensor gradients don't change, so pre-calculate them. 280 dAi_dAxx(self.dA[0]) 281 dAi_dAyy(self.dA[1]) 282 dAi_dAxy(self.dA[2]) 283 dAi_dAxz(self.dA[3]) 284 dAi_dAyz(self.dA[4]) 285 286 # PCS errors. 287 if self.pcs_flag_sum: 288 err = False 289 for i in range(len(pcs_errors)): 290 for j in range(len(pcs_errors[i])): 291 if not isNaN(pcs_errors[i, j]): 292 err = True 293 if err: 294 self.pcs_sigma_ij = pcs_errors 295 else: 296 # Missing errors (the values need to be small, close to ppm units, so the chi-squared value is comparable to the RDC). 297 self.pcs_sigma_ij = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64) 298 299 # RDC errors. 300 if self.rdc_flag_sum: 301 err = False 302 for i in range(len(rdc_errors)): 303 for j in range(len(rdc_errors[i])): 304 if not isNaN(rdc_errors[i, j]): 305 err = True 306 if err: 307 self.rdc_sigma_ij = rdc_errors 308 else: 309 # Missing errors. 310 self.rdc_sigma_ij = ones((self.num_align, self.num_interatom), float64) 311 312 # Missing data matrices (RDC). 313 if self.rdc_flag_sum: 314 self.missing_Dij = zeros((self.num_align, self.num_interatom), float64) 315 316 # Missing data matrices (PCS). 317 if self.pcs_flag_sum: 318 self.missing_deltaij = zeros((self.num_align, self.num_spins), float64) 319 320 # Clean up problematic data and put the weights into the errors.. 321 if self.rdc_flag_sum or self.pcs_flag_sum: 322 for align_index in range(self.num_align): 323 if self.rdc_flag_sum: 324 for j in range(self.num_interatom): 325 if isNaN(self.Dij[align_index, j]): 326 # Set the flag. 327 self.missing_Dij[align_index, j] = 1 328 329 # Change the NaN to zero. 330 self.Dij[align_index, j] = 0.0 331 332 # Change the error to one (to avoid zero division). 333 self.rdc_sigma_ij[align_index, j] = 1.0 334 335 # Change the weight to one. 336 rdc_weights[align_index, j] = 1.0 337 338 # The RDC weights. 339 self.rdc_sigma_ij[align_index, j] = self.rdc_sigma_ij[align_index, j] / sqrt(rdc_weights[align_index, j]) 340 341 342 if self.pcs_flag_sum: 343 for j in range(self.num_spins): 344 if isNaN(self.deltaij[align_index, j]): 345 # Set the flag. 346 self.missing_deltaij[align_index, j] = 1 347 348 # Change the NaN to zero. 349 self.deltaij[align_index, j] = 0.0 350 351 # Change the error to one (to avoid zero division). 352 self.pcs_sigma_ij[align_index, j] = 1.0 353 354 # Change the weight to one. 355 pcs_weights[align_index, j] = 1.0 356 357 # The PCS weights. 358 self.pcs_sigma_ij[align_index, j] = self.pcs_sigma_ij[align_index, j] / sqrt(pcs_weights[align_index, j]) 359 360 # The paramagnetic centre vectors and distances. 361 if self.pcs_flag_sum: 362 # Initialise the data structures. 363 self.paramag_unit_vect = zeros(atomic_pos.shape, float64) 364 self.paramag_dist = zeros((self.num_spins, self.N), float64) 365 self.pcs_const = zeros((self.num_align, self.num_spins, self.N), float64) 366 if self.paramag_centre == None: 367 self.paramag_centre = zeros(3, float64) 368 369 # Set up the paramagnetic info. 370 self.paramag_info() 371 372 # PCS function, gradient, and Hessian matrices. 373 self.deltaij_theta = zeros((self.num_align, self.num_spins), float64) 374 self.ddeltaij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 375 self.d2deltaij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 376 377 # RDC function, gradient, and Hessian matrices. 378 self.Dij_theta = zeros((self.num_align, self.num_interatom), float64) 379 self.dDij_theta = zeros((self.total_num_params, self.num_align, self.num_interatom), float64) 380 self.d2Dij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_interatom), float64) 381 382 # Set the target function, gradient, and Hessian (paramagnetic centre optimisation). 383 if not self.centre_fixed: 384 self.func = self.func_population 385 self.dfunc = None 386 self.d2func = None 387 388 # Set the standard target function, gradient, and Hessian. 389 else: 390 self.func = self.func_population 391 self.dfunc = self.dfunc_population 392 self.d2func = self.d2func_population 393 394 # Fixed probabilities. 395 if model == 'fixed': 396 # The probs are unpacked by self.func in the population model, so just override that function. 397 self.func = self.func_tensor_opt 398 self.dfunc = self.dfunc_tensor_opt 399 self.d2func = self.d2func_tensor_opt 400 401 # The zero Hessian. 402 self.zero_hessian_rdc = zeros(self.num_interatom, float64) 403 self.zero_hessian_pcs = zeros(self.num_spins, float64) 404 405 # The probability array. 406 if probs: 407 self.probs = probs 408 409 # All structures have initial equal probability. 410 else: 411 self.probs = ones(self.N, float64) / self.N

412 413

414 - def func_2domain(self, params):

415 """The target function for optimisation of the 2-domain N-state model. 416 417 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no tensor errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared). 418 419 @param params: The vector of parameter values. 420 @type params: list of float 421 @return: The chi-squared or SSE value. 422 @rtype: float 423 """ 424 425 # Scaling. 426 if self.scaling_flag: 427 params = dot(params, self.scaling_matrix) 428 429 # Reset the back-calculated the reduced tensor structure. 430 self.red_bc = self.red_bc * 0.0 431 432 # Loop over the N states. 433 for c in range(self.N): 434 # The rotation matrix. 435 euler_to_R_zyz(params[self.N-1+3*c], params[self.N-1+3*c+1], params[self.N-1+3*c+2], self.R[c]) 436 437 # Its transpose. 438 self.RT[c] = transpose(self.R[c]) 439 440 # The probability of state c. 441 if c < self.N-1: 442 pc = params[c] 443 444 # The probability of state N (1 minus the pc of all other states). 445 else: 446 pc = 1.0 447 for c2 in range(self.N-1): 448 pc = pc - params[c2] 449 450 # Back-calculate the reduced tensors for sum element c and add these to red_bc. 451 for align_index in range(self.num_tensors): 452 # Normal RT.X.R rotation. 453 if self.full_in_ref_frame[align_index]: 454 self.red_bc[align_index] = self.red_bc[align_index] + pc * dot(self.RT[c], dot(self.A[align_index], self.R[c])) 455 456 # Inverse R.X.RT rotation. 457 else: 458 self.red_bc[align_index] = self.red_bc[align_index] + pc * dot(self.R[c], dot(self.A[align_index], self.RT[c])) 459 460 # 5D vectorise the back-calculated tensors (create red_bc_vector from red_bc). 461 for align_index in range(self.num_tensors): 462 self.red_bc_vector[5*align_index] = self.red_bc[align_index, 0, 0] # Sxx. 463 self.red_bc_vector[5*align_index+1] = self.red_bc[align_index, 1, 1] # Syy. 464 self.red_bc_vector[5*align_index+2] = self.red_bc[align_index, 0, 1] # Sxy. 465 self.red_bc_vector[5*align_index+3] = self.red_bc[align_index, 0, 2] # Sxz. 466 self.red_bc_vector[5*align_index+4] = self.red_bc[align_index, 1, 2] # Syz. 467 468 # Return the chi-squared value. 469 return chi2(self.red_data, self.red_bc_vector, self.red_errors)

470 471

472 - def func_population(self, params):

473 """The target function for optimisation of the flexible population N-state model. 474 475 Description 476 =========== 477 478 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared). 479 480 481 Indices 482 ======= 483 484 For this calculation, five indices are looped over and used in the various data structures. These include: 485 - i, the index over alignments, 486 - j, the index over spin systems, 487 - c, the index over the N-states (or over the structures), 488 - n, the index over the first dimension of the alignment tensor n = {x, y, z}, 489 - m, the index over the second dimension of the alignment tensor m = {x, y, z}. 490 491 492 Equations 493 ========= 494 495 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC equation. 496 497 498 The chi-squared equation 499 ------------------------ 500 501 The equations are:: 502 503 ___ 504 \ (Dij - Dij(theta)) ** 2 505 chi^2(theta) = > ----------------------- , 506 /__ sigma_ij ** 2 507 ij 508 509 ___ 510 \ (delta_ij - delta_ij(theta)) ** 2 511 chi^2(theta) = > --------------------------------- , 512 /__ sigma_ij ** 2 513 ij 514 515 where: 516 - theta is the parameter vector, 517 - Dij are the measured RDCs for alignment i, spin j, 518 - Dij(theta) are the back calculated RDCs for alignment i, spin j, 519 - delta_ij are the measured PCSs for alignment i, spin j, 520 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j, 521 - sigma_ij are the RDC or PCS errors. 522 523 Both chi-squared values sum. 524 525 526 The RDC equation 527 ---------------- 528 529 The RDC equation is:: 530 531 _N_ 532 \ T 533 Dij(theta) = dj > pc . mu_jc . Ai . mu_jc, 534 /__ 535 c=1 536 537 where: 538 - dj is the dipolar constant for spin j, 539 - N is the total number of states or structures, 540 - pc is the weight or probability associated with state c, 541 - mu_jc is the unit vector corresponding to spin j and state c, 542 - Ai is the alignment tensor. 543 544 The dipolar constant is henceforth defined as:: 545 546 dj = 3 / (2pi) d', 547 548 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is:: 549 550 mu0 gI.gS.h_bar 551 d' = - --- ----------- , 552 4pi r**3 553 554 where: 555 - mu0 is the permeability of free space, 556 - gI and gS are the gyromagnetic ratios of the I and S spins, 557 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 558 - r is the distance between the two spins. 559 560 561 The PCS equation 562 ---------------- 563 564 The PCS equation is:: 565 566 _N_ 567 \ T 568 delta_ij(theta) = > pc . dijc . mu_jc . Ai . mu_jc, 569 /__ 570 c=1 571 572 where: 573 - djci is the PCS constant for spin j, state c and experiment or alignment i, 574 - N is the total number of states or structures, 575 - pc is the weight or probability associated with state c, 576 - mu_jc is the unit vector corresponding to spin j and state c, 577 - Ai is the alignment tensor. 578 579 The PCS constant is defined as:: 580 581 mu0 15kT 1 582 dijc = --- ----- ---- , 583 4pi Bo**2 r**3 584 585 where: 586 - mu0 is the permeability of free space, 587 - k is Boltzmann's constant, 588 - T is the absolute temperature (different for each experiment), 589 - Bo is the magnetic field strength (different for each experiment), 590 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state). 591 592 593 Stored data structures 594 ====================== 595 596 There are a number of data structures calculated by this function and stored for subsequent use in the gradient and Hessian functions. This include the back calculated RDCs and the alignment tensors. 597 598 Dij(theta) 599 ---------- 600 601 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}. 602 603 delta_ij(theta) 604 --------------- 605 606 The back calculated PCS. This is a rank-2 tensor with indices {i, j}. 607 608 Ai 609 -- 610 611 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}. 612 613 614 @param params: The vector of parameter values. 615 @type params: numpy rank-1 array 616 @return: The chi-squared or SSE value. 617 @rtype: float 618 """ 619 620 # Scaling. 621 if self.scaling_flag: 622 params = dot(params, self.scaling_matrix) 623 624 # Initial chi-squared (or SSE) value. 625 chi2_sum = 0.0 626 627 # Unpack the probabilities (located at the end of the parameter array). 628 if self.N > 1: 629 self.probs = params[-(self.N-1):] 630 631 # Unpack the paramagnetic centre. 632 if not self.centre_fixed: 633 # The position. 634 self.paramag_centre = params[-3:] 635 636 # Update the paramagnetic info. 637 self.paramag_info() 638 639 # Loop over each alignment. 640 index = 0 641 for align_index in range(self.num_align): 642 # Create tensor i from the parameters. 643 if not self.fixed_tensors[align_index]: 644 to_tensor(self.A[align_index], params[5*index:5*index + 5]) 645 index += 1 646 647 # The back calculated RDC. 648 if self.rdc_flag[align_index]: 649 # Loop over the spin pairs k. 650 for j in range(self.num_interatom): 651 # Calculate the average RDC. 652 if not self.missing_Dij[align_index, j]: 653 self.Dij_theta[align_index, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[align_index], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 654 655 # The back calculated PCS. 656 if self.pcs_flag[align_index]: 657 # Loop over the spin systems j. 658 for j in range(self.num_spins): 659 # Calculate the average PCS. 660 if not self.missing_deltaij[align_index, j]: 661 self.deltaij_theta[align_index, j] = ave_pcs_tensor(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.A[align_index], weights=self.probs) 662 663 # Calculate and sum the single alignment chi-squared value (for the RDC). 664 if self.rdc_flag[align_index]: 665 chi2_sum = chi2_sum + chi2(self.Dij[align_index], self.Dij_theta[align_index], self.rdc_sigma_ij[align_index]) 666 667 # Calculate and sum the single alignment chi-squared value (for the PCS). 668 if self.pcs_flag[align_index]: 669 chi2_sum = chi2_sum + chi2(self.deltaij[align_index], self.deltaij_theta[align_index], self.pcs_sigma_ij[align_index]) 670 671 # Return the chi-squared value. 672 return chi2_sum

673 674

675 - def func_tensor_opt(self, params):

676 """The target function for optimisation of the alignment tensor from RDC and/or PCS data. 677 678 Description 679 =========== 680 681 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no RDC or PCS errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared). 682 683 684 Indices 685 ======= 686 687 For this calculation, five indices are looped over and used in the various data structures. These include: 688 - i, the index over alignments, 689 - j, the index over spin systems, 690 - c, the index over the N-states (or over the structures), 691 - n, the index over the first dimension of the alignment tensor n = {x, y, z}, 692 - m, the index over the second dimension of the alignment tensor m = {x, y, z}. 693 694 695 Equations 696 ========= 697 698 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC and PCS equations. 699 700 701 The chi-squared equation 702 ------------------------ 703 704 The equations are:: 705 706 ___ 707 \ (Dij - Dij(theta)) ** 2 708 chi^2(theta) = > ----------------------- , 709 /__ sigma_ij ** 2 710 ij 711 712 ___ 713 \ (delta_ij - delta_ij(theta)) ** 2 714 chi^2(theta) = > --------------------------------- , 715 /__ sigma_ij ** 2 716 ij 717 718 where: 719 - theta is the parameter vector, 720 - Dij are the measured RDCs for alignment i, spin j, 721 - Dij(theta) are the back calculated RDCs for alignment i, spin j, 722 - delta_ij are the measured PCSs for alignment i, spin j, 723 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j, 724 - sigma_ij are the RDC or PCS errors. 725 726 Both chi-squared values sum. 727 728 729 The RDC equation 730 ---------------- 731 732 The RDC equation is:: 733 734 _N_ 735 dj \ T 736 Dij(theta) = -- > mu_jc . Ai . mu_jc, 737 N /__ 738 c=1 739 740 where: 741 - dj is the dipolar constant for spin j, 742 - N is the total number of states or structures, 743 - mu_jc is the unit vector corresponding to spin j and state c, 744 - Ai is the alignment tensor. 745 746 The dipolar constant is henceforth defined as:: 747 748 dj = 3 / (2pi) d', 749 750 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is:: 751 752 mu0 gI.gS.h_bar 753 d' = - --- ----------- , 754 4pi r**3 755 756 where: 757 - mu0 is the permeability of free space, 758 - gI and gS are the gyromagnetic ratios of the I and S spins, 759 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 760 - r is the distance between the two spins. 761 762 763 The PCS equation 764 ---------------- 765 766 The PCS equation is:: 767 768 _N_ 769 1 \ T 770 delta_ij(theta) = - > dijc . mu_jc . Ai . mu_jc, 771 N /__ 772 c=1 773 774 where: 775 - djci is the PCS constant for spin j, state c and experiment or alignment i, 776 - N is the total number of states or structures, 777 - mu_jc is the unit vector corresponding to spin j and state c, 778 - Ai is the alignment tensor. 779 780 The PCS constant is defined as:: 781 782 mu0 15kT 1 783 dijc = --- ----- ---- , 784 4pi Bo**2 r**3 785 786 where: 787 - mu0 is the permeability of free space, 788 - k is Boltzmann's constant, 789 - T is the absolute temperature (different for each experiment), 790 - Bo is the magnetic field strength (different for each experiment), 791 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state). 792 793 794 Stored data structures 795 ====================== 796 797 There are a number of data structures calculated by this function and stored for subsequent use in the gradient and Hessian functions. This include the back calculated RDCs and the alignment tensors. 798 799 Dij(theta) 800 ---------- 801 802 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}. 803 804 delta_ij(theta) 805 --------------- 806 807 The back calculated PCS. This is a rank-2 tensor with indices {i, j}. 808 809 Ai 810 -- 811 812 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}. 813 814 815 @param params: The vector of parameter values. 816 @type params: numpy rank-1 array 817 @return: The chi-squared or SSE value. 818 @rtype: float 819 """ 820 821 # Scaling. 822 if self.scaling_flag: 823 params = dot(params, self.scaling_matrix) 824 825 # Initial chi-squared (or SSE) value. 826 chi2_sum = 0.0 827 828 # Unpack the paramagnetic centre. 829 if not self.centre_fixed: 830 # The position. 831 self.paramag_centre = params[-3:] 832 833 # Update the paramagnetic info. 834 self.paramag_info() 835 836 # Loop over each alignment. 837 index = 0 838 for align_index in range(self.num_align): 839 # Create tensor i from the parameters. 840 if not self.fixed_tensors[align_index]: 841 to_tensor(self.A[align_index], params[5*index:5*index + 5]) 842 index += 1 843 844 # The back calculated RDC. 845 if self.rdc_flag[align_index]: 846 # Loop over the spin pairs k. 847 for j in range(self.num_interatom): 848 # Calculate the average RDC. 849 if not self.missing_Dij[align_index, j]: 850 self.Dij_theta[align_index, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[align_index], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 851 852 # The back calculated PCS. 853 if self.pcs_flag[align_index]: 854 # Loop over the spin systems j. 855 for j in range(self.num_spins): 856 # Calculate the average PCS. 857 if not self.missing_deltaij[align_index, j]: 858 self.deltaij_theta[align_index, j] = ave_pcs_tensor(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.A[align_index], weights=self.probs) 859 860 # Calculate and sum the single alignment chi-squared value (for the RDC). 861 if self.rdc_flag[align_index]: 862 chi2_sum = chi2_sum + chi2(self.Dij[align_index], self.Dij_theta[align_index], self.rdc_sigma_ij[align_index]) 863 864 # Calculate and sum the single alignment chi-squared value (for the PCS). 865 if self.pcs_flag[align_index]: 866 chi2_sum = chi2_sum + chi2(self.deltaij[align_index], self.deltaij_theta[align_index], self.pcs_sigma_ij[align_index]) 867 868 # Return the chi-squared value. 869 return chi2_sum

870 871

872 - def dfunc_population(self, params):

873 """The gradient function for optimisation of the flexible population N-state model. 874 875 Description 876 =========== 877 878 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared gradient corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) gradient is returned instead. The chi-squared gradient is simply the SSE gradient normalised to unit variance (the SSE divided by the error squared). 879 880 881 Indices 882 ======= 883 884 For this calculation, six indices are looped over and used in the various data structures. These include: 885 - k, the index over all parameters, 886 - i, the index over alignments, 887 - j, the index over spin systems, 888 - c, the index over the N-states (or over the structures), 889 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 890 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 891 892 893 Equations 894 ========= 895 896 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 897 898 899 The chi-squared gradient 900 ------------------------ 901 902 The equation is:: 903 ___ 904 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \ 905 ------------- = -2 > | ---------------- . ----------- | 906 dthetak /__ \ sigma_ij**2 dthetak / 907 ij 908 909 where: 910 - theta is the parameter vector, 911 - Dij are the measured RDCs or PCSs, 912 - Dij(theta) are the back calculated RDCs or PCSs, 913 - sigma_ij are the RDC or PCS errors, 914 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 915 916 917 The RDC gradient 918 ---------------- 919 920 This gradient is different for the various parameter types. 921 922 pc partial derivative 923 ~~~~~~~~~~~~~~~~~~~~~ 924 925 The population parameter partial derivative is:: 926 927 dDij(theta) T 928 ----------- = dj . mu_jc . Ai . mu_jc, 929 dpc 930 931 where: 932 - dj is the dipolar constant for spin j, 933 - mu_jc is the unit vector corresponding to spin j and state c, 934 - Ai is the alignment tensor. 935 936 Amn partial derivative 937 ~~~~~~~~~~~~~~~~~~~~~~ 938 939 The alignment tensor element partial derivative is:: 940 941 _N_ 942 dDij(theta) \ T dAi 943 ----------- = dj > pc . mu_jc . ---- . mu_jc, 944 dAmn /__ dAmn 945 c=1 946 947 where: 948 - dj is the dipolar constant for spin j, 949 - pc is the weight or probability associated with state c, 950 - mu_jc is the unit vector corresponding to spin j and state c, 951 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 952 953 954 The PCS gradient 955 ---------------- 956 957 This gradient is also different for the various parameter types. 958 959 pc partial derivative 960 ~~~~~~~~~~~~~~~~~~~~~ 961 962 The population parameter partial derivative is:: 963 964 ddeltaij(theta) T 965 --------------- = dijc . mu_jc . Ai . mu_jc, 966 dpc 967 968 where: 969 - djc is the pseudocontact shift constant for spin j and state c, 970 - mu_jc is the unit vector corresponding to spin j and state c, 971 - Ai is the alignment tensor. 972 973 Amn partial derivative 974 ~~~~~~~~~~~~~~~~~~~~~~ 975 976 The alignment tensor element partial derivative is:: 977 978 _N_ 979 ddelta_ij(theta) \ T dAi 980 ---------------- = > pc . djc . mu_jc . ---- . mu_jc, 981 dAmn /__ dAmn 982 c=1 983 984 where: 985 - djc is the pseudocontact shift constant for spin j and state c, 986 - pc is the weight or probability associated with state c, 987 - mu_jc is the unit vector corresponding to spin j and state c, 988 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 989 990 991 The alignment tensor gradient 992 ----------------------------- 993 994 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are:: 995 996 dAi | 1 0 0 | 997 ---- = | 0 0 0 |, 998 dAxx | 0 0 -1 | 999 1000 dAi | 0 0 0 | 1001 ---- = | 0 1 0 |, 1002 dAyy | 0 0 -1 | 1003 1004 dAi | 0 1 0 | 1005 ---- = | 1 0 0 |, 1006 dAxy | 0 0 0 | 1007 1008 dAi | 0 0 1 | 1009 ---- = | 0 0 0 |, 1010 dAxz | 1 0 0 | 1011 1012 dAi | 0 0 0 | 1013 ---- = | 0 0 1 |. 1014 dAyz | 0 1 0 | 1015 1016 As these are invariant, they can be pre-calculated. 1017 1018 1019 Stored data structures 1020 ====================== 1021 1022 There are a number of data structures calculated by this function and stored for subsequent use in the Hessian function. This include the back calculated RDC and PCS gradients and the alignment tensor gradients. 1023 1024 dDij(theta)/dthetak 1025 ------------------- 1026 1027 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}. 1028 1029 ddeltaij(theta)/dthetak 1030 ----------------------- 1031 1032 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}. 1033 1034 dAi/dAmn 1035 -------- 1036 1037 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}. 1038 1039 1040 @param params: The vector of parameter values. This is unused as it is assumed that 1041 func_population() was called first. 1042 @type params: numpy rank-1 array 1043 @return: The chi-squared or SSE gradient. 1044 @rtype: numpy rank-1 array 1045 """ 1046 1047 # Scaling. 1048 if self.scaling_flag: 1049 params = dot(params, self.scaling_matrix) 1050 1051 # Initial chi-squared (or SSE) gradient. 1052 self.dchi2 = self.dchi2 * 0.0 1053 1054 # Loop over each alignment. 1055 index = 0 1056 for align_index in range(self.num_align): 1057 # Fixed tensor, so skip. 1058 if self.fixed_tensors[align_index]: 1059 continue 1060 1061 # Construct the Amn partial derivative components for the RDC. 1062 if self.fixed_tensors[align_index] and self.rdc_flag[align_index] and not self.missing_Dij[align_index, j]: 1063 for j in range(self.num_interatom): 1064 self.dDij_theta[align_index*5, align_index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1065 self.dDij_theta[align_index*5+1, align_index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1066 self.dDij_theta[align_index*5+2, align_index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1067 self.dDij_theta[align_index*5+3, align_index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1068 self.dDij_theta[align_index*5+4, align_index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1069 1070 # Construct the Amn partial derivative components for the PCS. 1071 if self.fixed_tensors[align_index] and self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 1072 for j in range(self.num_spins): 1073 self.ddeltaij_theta[align_index*5, align_index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.dA[0], weights=self.probs) 1074 self.ddeltaij_theta[align_index*5+1, align_index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.dA[1], weights=self.probs) 1075 self.ddeltaij_theta[align_index*5+2, align_index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.dA[2], weights=self.probs) 1076 self.ddeltaij_theta[align_index*5+3, align_index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.dA[3], weights=self.probs) 1077 self.ddeltaij_theta[align_index*5+4, align_index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[align_index, j], self.paramag_unit_vect[j], self.N, self.dA[4], weights=self.probs) 1078 1079 # Construct the pc partial derivative gradient components, looping over each state. 1080 for c in range(self.N - 1): 1081 # Index in the parameter array. 1082 param_index = self.num_align_params + c 1083 1084 # Calculate the RDC for state c (this is the pc partial derivative). 1085 for j in range(self.num_interatom): 1086 if self.rdc_flag[align_index] and not self.missing_Dij[align_index, j]: 1087 self.dDij_theta[param_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.A[align_index], absolute=self.absolute_rdc[align_index, j]) 1088 1089 # Calculate the PCS for state c (this is the pc partial derivative). 1090 for j in range(self.num_spins): 1091 if self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 1092 self.ddeltaij_theta[param_index, align_index, j] = pcs_tensor(self.pcs_const[align_index, j, c], self.paramag_unit_vect[j, c], self.A[align_index]) 1093 1094 # Construct the chi-squared gradient element for parameter k, alignment i. 1095 for k in range(self.total_num_params): 1096 # RDC part of the chi-squared gradient. 1097 if self.rdc_flag[align_index]: 1098 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.Dij[align_index], self.Dij_theta[align_index], self.dDij_theta[k, align_index], self.rdc_sigma_ij[align_index]) 1099 1100 # PCS part of the chi-squared gradient. 1101 if self.pcs_flag[align_index]: 1102 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.deltaij[align_index], self.deltaij_theta[align_index], self.ddeltaij_theta[k, align_index], self.pcs_sigma_ij[align_index]) 1103 1104 # Increment the index. 1105 index += 1 1106 1107 # Diagonal scaling. 1108 if self.scaling_flag: 1109 self.dchi2 = dot(self.dchi2, self.scaling_matrix) 1110 1111 # The gradient. 1112 return self.dchi2

1113 1114

1115 - def dfunc_tensor_opt(self, params):

1116 """The gradient function for optimisation of the alignment tensor from RDC and/or PCS data. 1117 1118 Description 1119 =========== 1120 1121 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared gradient corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) gradient is returned instead. The chi-squared gradient is simply the SSE gradient normalised to unit variance (the SSE divided by the error squared). 1122 1123 Indices 1124 ======= 1125 1126 For this calculation, six indices are looped over and used in the various data structures. These include: 1127 - k, the index over all parameters, 1128 - i, the index over alignments, 1129 - j, the index over spin systems, 1130 - c, the index over the N-states (or over the structures), 1131 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1132 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1133 1134 1135 Equations 1136 ========= 1137 1138 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 1139 1140 1141 The chi-squared gradient 1142 ------------------------ 1143 1144 The equation is:: 1145 ___ 1146 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \ 1147 ------------- = -2 > | ---------------- . ----------- | 1148 dthetak /__ \ sigma_ij**2 dthetak / 1149 ij 1150 1151 where: 1152 - theta is the parameter vector, 1153 - Dij are the measured RDCs or PCSs, 1154 - Dij(theta) are the back calculated RDCs or PCSs, 1155 - sigma_ij are the RDC or PCS errors, 1156 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1157 1158 1159 The RDC gradient 1160 ---------------- 1161 1162 The only parameters are the tensor components. 1163 1164 Amn partial derivative 1165 ~~~~~~~~~~~~~~~~~~~~~~ 1166 1167 The alignment tensor element partial derivative is:: 1168 1169 _N_ 1170 dDij(theta) dj \ T dAi 1171 ----------- = -- > mu_jc . ---- . mu_jc, 1172 dAmn N /__ dAmn 1173 c=1 1174 1175 where: 1176 - dj is the dipolar constant for spin j, 1177 - N is the total number of states or structures, 1178 - mu_jc is the unit vector corresponding to spin j and state c, 1179 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 1180 1181 1182 The PCS gradient 1183 ---------------- 1184 1185 Amn partial derivative 1186 ~~~~~~~~~~~~~~~~~~~~~~ 1187 1188 The alignment tensor element partial derivative is:: 1189 1190 _N_ 1191 ddelta_ij(theta) 1 \ T dAi 1192 ---------------- = - > djc . mu_jc . ---- . mu_jc, 1193 dAmn N /__ dAmn 1194 c=1 1195 1196 where: 1197 - djc is the pseudocontact shift constant for spin j and state c, 1198 - N is the total number of states or structures, 1199 - mu_jc is the unit vector corresponding to spin j and state c, 1200 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 1201 1202 1203 The alignment tensor gradient 1204 ----------------------------- 1205 1206 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are:: 1207 1208 dAi | 1 0 0 | 1209 ---- = | 0 0 0 |, 1210 dAxx | 0 0 -1 | 1211 1212 dAi | 0 0 0 | 1213 ---- = | 0 1 0 |, 1214 dAyy | 0 0 -1 | 1215 1216 dAi | 0 1 0 | 1217 ---- = | 1 0 0 |, 1218 dAxy | 0 0 0 | 1219 1220 dAi | 0 0 1 | 1221 ---- = | 0 0 0 |, 1222 dAxz | 1 0 0 | 1223 1224 dAi | 0 0 0 | 1225 ---- = | 0 0 1 |. 1226 dAyz | 0 1 0 | 1227 1228 As these are invariant, they can be pre-calculated. 1229 1230 1231 Stored data structures 1232 ====================== 1233 1234 There are a number of data structures calculated by this function and stored for subsequent use in the Hessian function. This include the back calculated RDC and PCS gradients and the alignment tensor gradients. 1235 1236 dDij(theta)/dthetak 1237 ------------------- 1238 1239 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}. 1240 1241 ddeltaij(theta)/dthetak 1242 ----------------------- 1243 1244 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}. 1245 1246 dAi/dAmn 1247 -------- 1248 1249 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}. 1250 1251 1252 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first. 1253 @type params: numpy rank-1 array 1254 @return: The chi-squared or SSE gradient. 1255 @rtype: numpy rank-1 array 1256 """ 1257 1258 # Scaling. 1259 if self.scaling_flag: 1260 params = dot(params, self.scaling_matrix) 1261 1262 # Initial chi-squared (or SSE) gradient. 1263 self.dchi2 = self.dchi2 * 0.0 1264 1265 # Loop over each alignment. 1266 index = 0 1267 for align_index in range(self.num_align): 1268 # Fixed tensor, so skip. 1269 if self.fixed_tensors[align_index]: 1270 continue 1271 1272 # Construct the Amn partial derivative components for the RDC. 1273 for j in range(self.num_interatom): 1274 if self.rdc_flag[align_index] and not self.missing_Dij[align_index, j]: 1275 self.dDij_theta[index*5, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1276 self.dDij_theta[index*5+1, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1277 self.dDij_theta[index*5+2, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1278 self.dDij_theta[index*5+3, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1279 self.dDij_theta[index*5+4, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs, absolute=self.absolute_rdc[align_index, j]) 1280 1281 # Construct the Amn partial derivative components for the PCS. 1282 for j in range(self.num_spins): 1283 if self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 1284 self.ddeltaij_theta[index*5, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[0], weights=self.probs) 1285 self.ddeltaij_theta[index*5+1, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[1], weights=self.probs) 1286 self.ddeltaij_theta[index*5+2, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[2], weights=self.probs) 1287 self.ddeltaij_theta[index*5+3, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[3], weights=self.probs) 1288 self.ddeltaij_theta[index*5+4, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[4], weights=self.probs) 1289 1290 # Construct the chi-squared gradient element for parameter k, alignment i. 1291 for k in range(self.total_num_params): 1292 # RDC part of the chi-squared gradient. 1293 if self.rdc_flag[align_index]: 1294 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.Dij[index], self.Dij_theta[index], self.dDij_theta[k, index], self.rdc_sigma_ij[index]) 1295 1296 # PCS part of the chi-squared gradient. 1297 if self.pcs_flag[align_index]: 1298 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.deltaij[index], self.deltaij_theta[index], self.ddeltaij_theta[k, index], self.pcs_sigma_ij[index]) 1299 1300 # Increment the index. 1301 index += 1 1302 1303 # Diagonal scaling. 1304 if self.scaling_flag: 1305 self.dchi2 = dot(self.dchi2, self.scaling_matrix) 1306 1307 # The gradient. 1308 return self.dchi2

1309 1310

1311 - def d2func_population(self, params):

1312 """The Hessian function for optimisation of the flexible population N-state model. 1313 1314 Description 1315 =========== 1316 1317 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared Hessian corresponding to that coordinate in the parameter space. If no RDC/PCS errors are supplied, then the SSE (the sum of squares error) Hessian is returned instead. The chi-squared Hessian is simply the SSE Hessian normalised to unit variance (the SSE divided by the error squared). 1318 1319 1320 Indices 1321 ======= 1322 1323 For this calculation, six indices are looped over and used in the various data structures. These include: 1324 - k, the index over all parameters, 1325 - i, the index over alignments, 1326 - j, the index over spin systems, 1327 - c, the index over the N-states (or over the structures), 1328 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1329 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1330 1331 1332 Equations 1333 ========= 1334 1335 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 1336 1337 1338 The chi-squared Hessian 1339 ----------------------- 1340 1341 The equation is:: 1342 ___ 1343 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \ 1344 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |. 1345 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 1346 ij 1347 1348 where: 1349 - theta is the parameter vector, 1350 - Dij are the measured RDCs or PCSs, 1351 - Dij(theta) are the back calculated RDCs or PCSs, 1352 - sigma_ij are the RDC or PCS errors, 1353 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1354 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k. 1355 1356 1357 The RDC Hessian 1358 --------------- 1359 1360 pc-pd second partial derivatives 1361 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1362 1363 The probability parameter second partial derivative is:: 1364 1365 d2Dij(theta) 1366 ------------ = 0. 1367 dpc.dpd 1368 1369 1370 pc-Anm second partial derivatives 1371 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1372 1373 The probability parameter-tensor element second partial derivative is:: 1374 1375 d2Dij(theta) T dAi 1376 ------------ = dj . mu_jc . ---- . mu_jc. 1377 dpc.dAmn dAmn 1378 1379 1380 Amn-Aop second partial derivatives 1381 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1382 1383 The alignment tensor element second partial derivative is:: 1384 1385 d2Dij(theta) 1386 ------------ = 0. 1387 dAmn.dAop 1388 1389 1390 The PCS Hessian 1391 --------------- 1392 1393 pc-pd second partial derivatives 1394 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1395 1396 The probability parameter second partial derivative is:: 1397 1398 d2delta_ij(theta) 1399 ----------------- = 0. 1400 dpc.dpd 1401 1402 1403 pc-Anm second partial derivatives 1404 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1405 1406 The probability parameter-tensor element second partial derivative is:: 1407 1408 d2delta_ij(theta) T dAi 1409 ----------------- = djc . mu_jc . ---- . mu_jc. 1410 dpc.dAmn dAmn 1411 1412 1413 Amn-Aop second partial derivatives 1414 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1415 1416 The alignment tensor element second partial derivative is:: 1417 1418 d2delta_ij(theta) 1419 ----------------- = 0 1420 dAmn.dAop 1421 1422 1423 The alignment tensor Hessian 1424 ---------------------------- 1425 1426 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:: 1427 1428 d2Ai | 0 0 0 | 1429 --------- = | 0 0 0 |. 1430 dAmn.dAop | 0 0 0 | 1431 1432 1433 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first. 1434 @type params: numpy rank-1 array 1435 @return: The chi-squared or SSE Hessian. 1436 @rtype: numpy rank-2 array 1437 """ 1438 1439 # Scaling. 1440 if self.scaling_flag: 1441 params = dot(params, self.scaling_matrix) 1442 1443 # Initial chi-squared (or SSE) Hessian. 1444 self.d2chi2 = self.d2chi2 * 0.0 1445 1446 # Loop over each alignment. 1447 index = 0 1448 for align_index in range(self.num_align): 1449 # Fixed tensor, so skip. 1450 if self.fixed_tensors[align_index]: 1451 continue 1452 1453 # Construct the pc-Amn second partial derivative Hessian components. 1454 for c in range(self.N - 1): 1455 # Index in the parameter array. 1456 pc_index = self.num_align_params + c 1457 1458 # Calculate the RDC Hessian component. 1459 for j in range(self.num_interatom): 1460 if self.fixed_tensors[align_index] and self.rdc_flag[align_index] and not self.missing_Dij[align_index, j]: 1461 self.d2Dij_theta[pc_index, i*5+0, align_index, j] = self.d2Dij_theta[align_index*5+0, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[0], absolute=self.absolute_rdc[align_index, j]) 1462 self.d2Dij_theta[pc_index, i*5+1, align_index, j] = self.d2Dij_theta[align_index*5+1, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[1], absolute=self.absolute_rdc[align_index, j]) 1463 self.d2Dij_theta[pc_index, i*5+2, align_index, j] = self.d2Dij_theta[align_index*5+2, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[2], absolute=self.absolute_rdc[align_index, j]) 1464 self.d2Dij_theta[pc_index, i*5+3, align_index, j] = self.d2Dij_theta[align_index*5+3, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[3], absolute=self.absolute_rdc[align_index, j]) 1465 self.d2Dij_theta[pc_index, i*5+4, align_index, j] = self.d2Dij_theta[align_index*5+4, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[4], absolute=self.absolute_rdc[align_index, j]) 1466 1467 # Calculate the PCS Hessian component. 1468 for j in range(self.num_spins): 1469 if self.fixed_tensors[align_index] and self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 1470 self.d2deltaij_theta[pc_index, i*5+0, align_index, j] = self.d2deltaij_theta[align_index*5+0, pc_index, align_index, j] = pcs_tensor(self.pcs_const[align_index, j, c], self.paramag_unit_vect[j, c], self.dA[0]) 1471 self.d2deltaij_theta[pc_index, i*5+1, align_index, j] = self.d2deltaij_theta[align_index*5+1, pc_index, align_index, j] = pcs_tensor(self.pcs_const[align_index, j, c], self.paramag_unit_vect[j, c], self.dA[1]) 1472 self.d2deltaij_theta[pc_index, i*5+2, align_index, j] = self.d2deltaij_theta[align_index*5+2, pc_index, align_index, j] = pcs_tensor(self.pcs_const[align_index, j, c], self.paramag_unit_vect[j, c], self.dA[2]) 1473 self.d2deltaij_theta[pc_index, i*5+3, align_index, j] = self.d2deltaij_theta[align_index*5+3, pc_index, align_index, j] = pcs_tensor(self.pcs_const[align_index, j, c], self.paramag_unit_vect[j, c], self.dA[3]) 1474 self.d2deltaij_theta[pc_index, i*5+4, align_index, j] = self.d2deltaij_theta[align_index*5+4, pc_index, align_index, j] = pcs_tensor(self.pcs_const[align_index, j, c], self.paramag_unit_vect[j, c], self.dA[4]) 1475 1476 # Increment the index. 1477 index += 1 1478 1479 # Loop over each alignment. 1480 for align_index in range(self.num_align): 1481 # Construct the chi-squared Hessian element for parameters j and k, alignment i. 1482 for j in range(self.total_num_params): 1483 for k in range(self.total_num_params): 1484 # RDC part of the chi-squared gradient. 1485 if self.fixed_tensors[align_index] and self.rdc_flag[align_index]: 1486 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.Dij[align_index], self.Dij_theta[align_index], self.dDij_theta[j, align_index], self.dDij_theta[k, align_index], self.d2Dij_theta[j, k, align_index], self.rdc_sigma_ij[align_index]) 1487 1488 # PCS part of the chi-squared gradient. 1489 if self.fixed_tensors[align_index] and self.pcs_flag[align_index]: 1490 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.deltaij[align_index], self.deltaij_theta[align_index], self.ddeltaij_theta[j, align_index], self.ddeltaij_theta[k, align_index], self.d2deltaij_theta[j, k, align_index], self.pcs_sigma_ij[align_index]) 1491 1492 # Diagonal scaling. 1493 if self.scaling_flag: 1494 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix) 1495 1496 # The gradient. 1497 return self.d2chi2

1498 1499

1500 - def d2func_tensor_opt(self, params):

1501 """The Hessian function for optimisation of the alignment tensor from RDC and/or PCS data. 1502 1503 Description 1504 =========== 1505 1506 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared Hessian corresponding to that coordinate in the parameter space. If no RDC/PCS errors are supplied, then the SSE (the sum of squares error) Hessian is returned instead. The chi-squared Hessian is simply the SSE Hessian normalised to unit variance (the SSE divided by the error squared). 1507 1508 Indices 1509 ======= 1510 1511 For this calculation, six indices are looped over and used in the various data structures. These include: 1512 - k, the index over all parameters, 1513 - i, the index over alignments, 1514 - j, the index over spin systems, 1515 - c, the index over the N-states (or over the structures), 1516 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1517 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1518 1519 1520 Equations 1521 ========= 1522 1523 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 1524 1525 1526 The chi-squared Hessian 1527 ----------------------- 1528 1529 The equation is:: 1530 ___ 1531 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \ 1532 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |. 1533 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 1534 ij 1535 1536 where: 1537 - theta is the parameter vector, 1538 - Dij are the measured RDCs or PCSs, 1539 - Dij(theta) are the back calculated RDCs or PCSs, 1540 - sigma_ij are the RDC or PCS errors, 1541 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1542 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k. 1543 1544 1545 The RDC Hessian 1546 --------------- 1547 1548 The only parameters are the tensor components. 1549 1550 1551 Amn-Aop second partial derivatives 1552 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1553 1554 The alignment tensor element second partial derivative is:: 1555 1556 d2Dij(theta) 1557 ------------ = 0. 1558 dAmn.dAop 1559 1560 1561 The PCS Hessian 1562 --------------- 1563 1564 Amn-Aop second partial derivatives 1565 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1566 1567 The alignment tensor element second partial derivative is:: 1568 1569 d2delta_ij(theta) 1570 ----------------- = 0 1571 dAmn.dAop 1572 1573 1574 The alignment tensor Hessian 1575 ---------------------------- 1576 1577 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:: 1578 1579 d2Ai | 0 0 0 | 1580 --------- = | 0 0 0 |. 1581 dAmn.dAop | 0 0 0 | 1582 1583 1584 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first. 1585 @type params: numpy rank-1 array 1586 @return: The chi-squared or SSE Hessian. 1587 @rtype: numpy rank-2 array 1588 """ 1589 1590 # Scaling. 1591 if self.scaling_flag: 1592 params = dot(params, self.scaling_matrix) 1593 1594 # Initial chi-squared (or SSE) Hessian. 1595 self.d2chi2 = self.d2chi2 * 0.0 1596 1597 # Loop over each alignment. 1598 index = 0 1599 for align_index in range(self.num_align): 1600 # Fixed tensor, so skip. 1601 if self.fixed_tensors[align_index]: 1602 continue 1603 1604 # Construct the chi-squared Hessian element for parameters j and k, alignment i. 1605 for j in range(self.total_num_params): 1606 for k in range(self.total_num_params): 1607 # RDC part of the chi-squared gradient. 1608 if self.rdc_flag[align_index]: 1609 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.Dij[align_index], self.Dij_theta[align_index], self.dDij_theta[j, align_index], self.dDij_theta[k, align_index], self.zero_hessian_rdc, self.rdc_sigma_ij[align_index]) 1610 1611 # PCS part of the chi-squared gradient. 1612 if self.pcs_flag[align_index]: 1613 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.deltaij[align_index], self.deltaij_theta[align_index], self.ddeltaij_theta[j, align_index], self.ddeltaij_theta[k, align_index], self.zero_hessian_pcs, self.pcs_sigma_ij[align_index]) 1614 1615 # Increment the index. 1616 index += 1 1617 1618 # Diagonal scaling. 1619 if self.scaling_flag: 1620 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix) 1621 1622 # The gradient. 1623 return self.d2chi2

1624 1625

1626 - def paramag_info(self):

1627 """Calculate the paramagnetic centre to spin vectors, distances and constants.""" 1628 1629 # Get the vectors and distances. 1630 if rank(self.paramag_centre) == 1: 1631 vectors_single_centre(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1632 else: 1633 vectors_centre_per_state(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1634 1635 # The PCS constants. 1636 for align_index in range(self.num_align): 1637 for j in range(self.num_spins): 1638 for c in range(self.N): 1639 self.pcs_const[align_index, j, c] = pcs_constant(self.temp[align_index], self.frq[align_index], self.paramag_dist[j, c])

Source Code for Module maths_fns.n_state_model