target_functions.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, T_flags=None, j_couplings=None, rdc_pseudo_flags=None, pcs_pseudo_flags=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 - T_flags, the array of flags specifying if the data for the given alignment is of the T = J+D type. 101 - j_couplings, the J couplings for the case when RDC data is of the type T = J+D. 102 - dip_const, the dipolar constants. 103 - absolute_rdc, the flags specifying whether each RDC is signless. 104 105 106 @keyword model: The N-state model type. This can be one of '2-domain', 'population' or 'fixed'. 107 @type model: str 108 @keyword N: The number of states. 109 @type N: int 110 @keyword init_params: The initial parameter values. Optimisation must start at some point! 111 @type init_params: numpy float64 array 112 @keyword probs: The probabilities for each state c. The length of this list should be equal to N. 113 @type probs: list of float 114 @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] 115 @type full_tensors: list of rank-2, 3D numpy arrays 116 @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. 117 @type red_data: numpy float64 array 118 @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. 119 @type red_errors: numpy float64 array 120 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 121 @type full_in_ref_frame: numpy rank-1 array 122 @keyword fixed_tensors: An array of flags specifying if the tensor is fixed or will be optimised. 123 @type fixed_tensors: list of bool 124 @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. 125 @type pcs: numpy rank-2 array 126 @keyword pcs_errors: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 127 @type pcs_errors: numpy rank-2 array 128 @keyword pcs_weights: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 129 @type pcs_weights: numpy rank-2 array 130 @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. 131 @type rdcs: numpy rank-2 array 132 @keyword rdc_errors: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 133 @type rdc_errors: numpy rank-2 array 134 @keyword rdc_weights: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 135 @type rdc_weights: numpy rank-2 array 136 @keyword rdc_vect: The unit vector lists for the RDC. The first index must correspond to the spin pair, the second index to each structure (its size being equal to the number of states) and the third to each pseudo-atom present. 137 @type rdc_vect: list of numpy rank-3 array 138 @keyword T_flags: The array of flags specifying if the data for the given alignment is of the T = J+D type. 139 @type T_flags: numpy rank-2 int32 array 140 @keyword j_couplings: The J couplings list, used when the RDC data is of the type T = J+D. The number of elements must match the second dimension of the rdcs argument. 141 @type j_couplings: numpy rank-1 array 142 @keyword rdc_pseudo_flags: The array of flags specifying if one of the atoms of the interatomic pair for the RDC are pseudo-atoms. The indices correspond to the interatomic pairs. 143 @type rdc_pseudo_flags: numpy rank-1 int32 array 144 @keyword pcs_pseudo_flags: The array of flags specifying if one of the atoms of the interatomic pair for the PCS are pseudo-atoms. The indices correspond to the interatomic pairs. 145 @type pcs_pseudo_flags: numpy rank-1 int32 array 146 @keyword temp: The temperature of each experiment, used for the PCS. 147 @type temp: numpy rank-1 array 148 @keyword frq: The frequency of each alignment, used for the PCS. 149 @type frq: numpy rank-1 array 150 @keyword dip_const: The dipolar constants for each spin pair. The first index corresponds to the spin pairs k and the second to each pseudo-atom present. 151 @type dip_const: list of lists of floats 152 @keyword absolute_rdc: The signless or absolute value flags for the RDCs. 153 @type absolute_rdc: numpy rank-2 int32 array 154 @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. 155 @type atomic_pos: numpy rank-3 array 156 @keyword paramag_centre: The paramagnetic centre position (or positions). 157 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 158 @keyword scaling_matrix: The square and diagonal scaling matrix. 159 @type scaling_matrix: numpy rank-2 array 160 @keyword centre_fixed: A flag which if False will cause the paramagnetic centre to be optimised. 161 @type centre_fixed: bool 162 """ 163 164 # Store the data inside the class instance namespace. 165 self.N = N 166 self.params = 1.0 * init_params # Force a copy of the data to be stored. 167 self.fixed_tensors = fixed_tensors 168 self.deltaij = pcs 169 self.rdc = rdcs 170 self.dip_vect = rdc_vect 171 self.dip_const = dip_const 172 self.absolute_rdc = absolute_rdc 173 self.T_flags = T_flags 174 self.j_couplings = j_couplings 175 self.rdc_pseudo_flags = rdc_pseudo_flags 176 self.pcs_pseudo_flags = pcs_pseudo_flags 177 self.temp = temp 178 self.frq = frq 179 self.atomic_pos = atomic_pos 180 self.paramag_centre = paramag_centre 181 self.centre_fixed = centre_fixed 182 self.total_num_params = len(init_params) 183 184 # Initialise the function value, gradient, and Hessian. 185 self.chi2 = 0.0 186 self.dchi2 = zeros((self.total_num_params), float64) 187 self.d2chi2 = zeros((self.total_num_params, self.total_num_params), float64) 188 189 # Scaling initialisation. 190 self.scaling_matrix = scaling_matrix 191 if self.scaling_matrix != None: 192 self.scaling_flag = True 193 else: 194 self.scaling_flag = False 195 196 # The 2-domain N-state model. 197 if model == '2-domain': 198 # Some checks. 199 if red_data == None or not len(red_data): 200 raise RelaxError("The red_data argument " + repr(red_data) + " must be supplied.") 201 if red_errors == None or not len(red_errors): 202 raise RelaxError("The red_errors argument " + repr(red_errors) + " must be supplied.") 203 if full_in_ref_frame == None or not len(full_in_ref_frame): 204 raise RelaxError("The full_in_ref_frame argument " + repr(full_in_ref_frame) + " must be supplied.") 205 206 # Tensor set up. 207 self.full_tensors = array(full_tensors, float64) 208 self.num_tensors = int(len(self.full_tensors) / 5) 209 self.red_data = red_data 210 self.red_errors = red_errors 211 self.full_in_ref_frame = full_in_ref_frame 212 213 # Alignment tensor in rank-2, 3D form. 214 self.A = zeros((self.num_tensors, 3, 3), float64) 215 for align_index in range(self.num_tensors): 216 to_tensor(self.A[align_index], self.full_tensors[5*align_index:5*align_index+5]) 217 218 # Initialise some empty numpy objects for storage of: 219 # R: the transient rotation matricies. 220 # RT: the transposes of the rotation matricies. 221 # red_bc: the back-calculated reduced alignment tensors. 222 # red_bc_vector: the back-calculated reduced alignment tensors in vector form {Sxx, Syy, Sxy, Sxz, Syz}. 223 self.R = zeros((self.N, 3, 3), float64) 224 self.RT = zeros((self.N, 3, 3), float64) 225 self.red_bc = zeros((self.num_tensors, 3, 3), float64) 226 self.red_bc_vector = zeros(self.num_tensors*5, float64) 227 228 # Set the target function. 229 self.func = self.func_2domain 230 self.dfunc = None 231 self.d2func = None 232 233 # The flexible population or equal probability N-state models. 234 elif model in ['population', 'fixed']: 235 # The total number of alignments (must be the same for the RDC and PCS data). 236 self.num_align = 0 237 if rdcs != None: 238 self.num_align = len(rdcs) 239 if pcs != None: 240 self.num_align = max(self.num_align, len(pcs)) 241 242 # Set the RDC and PCS flags (indicating the presence of data). 243 self.rdc_flag = [True] * self.num_align 244 self.pcs_flag = [True] * self.num_align 245 for align_index in range(self.num_align): 246 if rdcs == None or len(rdcs[align_index]) == 0: 247 self.rdc_flag[align_index] = False 248 if pcs == None or len(pcs[align_index]) == 0: 249 self.pcs_flag[align_index] = False 250 self.rdc_flag_sum = sum(self.rdc_flag) 251 self.pcs_flag_sum = sum(self.pcs_flag) 252 253 # Some checks. 254 if self.rdc_flag_sum and (rdc_vect == None or not len(rdc_vect)): 255 raise RelaxError("The rdc_vect argument " + repr(rdc_vect) + " must be supplied.") 256 if self.pcs_flag_sum and (atomic_pos == None or not len(atomic_pos)): 257 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.") 258 259 # The total number of spins. 260 self.num_spins = 0 261 if self.pcs_flag_sum: 262 self.num_spins = len(pcs[0]) 263 264 # The total number of interatomic connections. 265 self.num_interatom = 0 266 if self.rdc_flag_sum: 267 self.num_interatom = len(rdcs[0]) 268 269 # Alignment tensor function and gradient matrices. 270 self.A = zeros((self.num_align, 3, 3), float64) 271 self.dA = zeros((5, 3, 3), float64) 272 273 # Fixed alignment tensors. 274 if full_tensors != None: 275 # Convert to numpy. 276 self.full_tensors = array(full_tensors, float64) 277 278 # Set up the alignment data. 279 self.num_align_params = 0 280 index = 0 281 for align_index in range(self.num_align): 282 # Fill the alignment tensor object with the fixed tensors. 283 if fixed_tensors[align_index]: 284 to_tensor(self.A[align_index], self.full_tensors[5*index:5*index+5]) 285 index += 1 286 287 # The number of alignment parameters. 288 if not fixed_tensors[align_index]: 289 self.num_align_params += 5 290 291 # Optimised alignment tensors. 292 else: 293 # The alignment tensor gradients don't change, so pre-calculate them. 294 dAi_dAxx(self.dA[0]) 295 dAi_dAyy(self.dA[1]) 296 dAi_dAxy(self.dA[2]) 297 dAi_dAxz(self.dA[3]) 298 dAi_dAyz(self.dA[4]) 299 300 # PCS errors. 301 if self.pcs_flag_sum: 302 err = False 303 for i in range(len(pcs_errors)): 304 for j in range(len(pcs_errors[i])): 305 if not isNaN(pcs_errors[i, j]): 306 err = True 307 if err: 308 self.pcs_errors = pcs_errors 309 else: 310 # Missing errors (the values need to be small, close to ppm units, so the chi-squared value is comparable to the RDC). 311 self.pcs_errors = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64) 312 313 # RDC errors. 314 if self.rdc_flag_sum: 315 err = False 316 for i in range(len(rdc_errors)): 317 for j in range(len(rdc_errors[i])): 318 if not isNaN(rdc_errors[i, j]): 319 err = True 320 if err: 321 self.rdc_errors = rdc_errors 322 else: 323 # Missing errors. 324 self.rdc_errors = ones((self.num_align, self.num_interatom), float64) 325 326 # Missing data matrices (RDC). 327 if self.rdc_flag_sum: 328 self.missing_rdc = zeros((self.num_align, self.num_interatom), float64) 329 330 # Missing data matrices (PCS). 331 if self.pcs_flag_sum: 332 self.missing_deltaij = zeros((self.num_align, self.num_spins), float64) 333 334 # Clean up problematic data and put the weights into the errors.. 335 if self.rdc_flag_sum or self.pcs_flag_sum: 336 for align_index in range(self.num_align): 337 if self.rdc_flag_sum: 338 for j in range(self.num_interatom): 339 if isNaN(self.rdc[align_index, j]): 340 # Set the flag. 341 self.missing_rdc[align_index, j] = 1 342 343 # Change the NaN to zero. 344 self.rdc[align_index, j] = 0.0 345 346 # Change the error to one (to avoid zero division). 347 self.rdc_errors[align_index, j] = 1.0 348 349 # Change the weight to one. 350 rdc_weights[align_index, j] = 1.0 351 352 # The RDC weights. 353 self.rdc_errors[align_index, j] = self.rdc_errors[align_index, j] / sqrt(rdc_weights[align_index, j]) 354 355 356 if self.pcs_flag_sum: 357 for j in range(self.num_spins): 358 if isNaN(self.deltaij[align_index, j]): 359 # Set the flag. 360 self.missing_deltaij[align_index, j] = 1 361 362 # Change the NaN to zero. 363 self.deltaij[align_index, j] = 0.0 364 365 # Change the error to one (to avoid zero division). 366 self.pcs_errors[align_index, j] = 1.0 367 368 # Change the weight to one. 369 pcs_weights[align_index, j] = 1.0 370 371 # The PCS weights. 372 self.pcs_errors[align_index, j] = self.pcs_errors[align_index, j] / sqrt(pcs_weights[align_index, j]) 373 374 # The paramagnetic centre vectors and distances. 375 if self.pcs_flag_sum: 376 # Initialise the data structures. 377 self.paramag_unit_vect = zeros(atomic_pos.shape, float64) 378 self.paramag_dist = zeros((self.num_spins, self.N), float64) 379 self.pcs_const = zeros((self.num_align, self.num_spins, self.N), float64) 380 if self.paramag_centre == None: 381 self.paramag_centre = zeros(3, float64) 382 383 # The gradient structures. 384 if not self.centre_fixed: 385 self.dpcs_const_theta = zeros((self.num_align, self.num_spins, self.N, 3), float64) 386 self.dr_theta = -eye(3) 387 388 # Set up the paramagnetic info. 389 self.paramag_info() 390 391 # PCS function, gradient, and Hessian matrices. 392 self.deltaij_theta = zeros((self.num_align, self.num_spins), float64) 393 self.ddeltaij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 394 self.d2deltaij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 395 396 # RDC function, gradient, and Hessian matrices. 397 self.rdc_theta = zeros((self.num_align, self.num_interatom), float64) 398 self.drdc_theta = zeros((self.total_num_params, self.num_align, self.num_interatom), float64) 399 self.d2rdc_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_interatom), float64) 400 401 # Set the target function, gradient, and Hessian. 402 self.func = self.func_standard 403 self.dfunc = self.dfunc_standard 404 self.d2func = self.d2func_standard 405 406 # Variable probabilities. 407 self.probs_fixed = True 408 if model == 'population': 409 self.probs_fixed = False 410 411 # Fixed probabilities. 412 if model == 'fixed': 413 # The zero Hessian. 414 self.zero_hessian_rdc = zeros(self.num_interatom, float64) 415 self.zero_hessian_pcs = zeros(self.num_spins, float64) 416 417 # The probability array. 418 if probs: 419 self.probs = probs 420 421 # All structures have initial equal probability. 422 else: 423 self.probs = ones(self.N, float64) / self.N

424 425

426 - def func_2domain(self, params):

427 """The target function for optimisation of the 2-domain N-state model. 428 429 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). 430 431 @param params: The vector of parameter values. 432 @type params: list of float 433 @return: The chi-squared or SSE value. 434 @rtype: float 435 """ 436 437 # Scaling. 438 if self.scaling_flag: 439 params = dot(params, self.scaling_matrix) 440 441 # Reset the back-calculated the reduced tensor structure. 442 self.red_bc = self.red_bc * 0.0 443 444 # Loop over the N states. 445 for c in range(self.N): 446 # The rotation matrix. 447 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]) 448 449 # Its transpose. 450 self.RT[c] = transpose(self.R[c]) 451 452 # The probability of state c. 453 if c < self.N-1: 454 pc = params[c] 455 456 # The probability of state N (1 minus the pc of all other states). 457 else: 458 pc = 1.0 459 for c2 in range(self.N-1): 460 pc = pc - params[c2] 461 462 # Back-calculate the reduced tensors for sum element c and add these to red_bc. 463 for align_index in range(self.num_tensors): 464 # Normal RT.X.R rotation. 465 if self.full_in_ref_frame[align_index]: 466 self.red_bc[align_index] = self.red_bc[align_index] + pc * dot(self.RT[c], dot(self.A[align_index], self.R[c])) 467 468 # Inverse R.X.RT rotation. 469 else: 470 self.red_bc[align_index] = self.red_bc[align_index] + pc * dot(self.R[c], dot(self.A[align_index], self.RT[c])) 471 472 # 5D vectorise the back-calculated tensors (create red_bc_vector from red_bc). 473 for align_index in range(self.num_tensors): 474 self.red_bc_vector[5*align_index] = self.red_bc[align_index, 0, 0] # Sxx. 475 self.red_bc_vector[5*align_index+1] = self.red_bc[align_index, 1, 1] # Syy. 476 self.red_bc_vector[5*align_index+2] = self.red_bc[align_index, 0, 1] # Sxy. 477 self.red_bc_vector[5*align_index+3] = self.red_bc[align_index, 0, 2] # Sxz. 478 self.red_bc_vector[5*align_index+4] = self.red_bc[align_index, 1, 2] # Syz. 479 480 # Return the chi-squared value. 481 return chi2(self.red_data, self.red_bc_vector, self.red_errors)

482 483

484 - def func_standard(self, params):

485 """The target function for optimisation of the standard N-state model. 486 487 Description 488 =========== 489 490 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 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). 491 492 493 Indices 494 ======= 495 496 For this calculation, five indices are looped over and used in the various data structures. These include: 497 - i, the index over alignments, 498 - j, the index over spin systems, 499 - c, the index over the N-states (or over the structures), 500 - n, the index over the first dimension of the alignment tensor n = {x, y, z}, 501 - m, the index over the second dimension of the alignment tensor m = {x, y, z}. 502 503 504 Equations 505 ========= 506 507 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC and PCS equations. 508 509 510 The chi-squared equation 511 ------------------------ 512 513 The equations are:: 514 515 ___ 516 \ (Dij - Dij(theta)) ** 2 517 chi^2(theta) = > ----------------------- , 518 /__ sigma_ij ** 2 519 ij 520 521 ___ 522 \ (delta_ij - delta_ij(theta)) ** 2 523 chi^2(theta) = > --------------------------------- , 524 /__ sigma_ij ** 2 525 ij 526 527 where: 528 - theta is the parameter vector, 529 - Dij are the measured RDCs for alignment i, spin j, 530 - Dij(theta) are the back calculated RDCs for alignment i, spin j, 531 - delta_ij are the measured PCSs for alignment i, spin j, 532 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j, 533 - sigma_ij are the RDC or PCS errors. 534 535 Both chi-squared values sum. 536 537 538 The RDC equation 539 ---------------- 540 541 The RDC equation is:: 542 543 _N_ 544 \ T 545 Dij(theta) = dj > pc . mu_jc . Ai . mu_jc, 546 /__ 547 c=1 548 549 where: 550 - dj is the dipolar constant for spin j, 551 - N is the total number of states or structures, 552 - pc is the weight or probability associated with state c, 553 - mu_jc is the unit vector corresponding to spin j and state c, 554 - Ai is the alignment tensor. 555 556 In the fixed and equal probability case, the equation is:: 557 558 _N_ 559 dj \ T 560 Dij(theta) = -- > mu_jc . Ai . mu_jc, 561 N /__ 562 c=1 563 564 The dipolar constant is henceforth defined as:: 565 566 dj = 3 / (2pi) d', 567 568 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:: 569 570 mu0 gI.gS.h_bar 571 d' = - --- ----------- , 572 4pi r**3 573 574 where: 575 - mu0 is the permeability of free space, 576 - gI and gS are the gyromagnetic ratios of the I and S spins, 577 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 578 - r is the distance between the two spins. 579 580 581 The PCS equation 582 ---------------- 583 584 The PCS equation is:: 585 586 _N_ 587 \ T 588 delta_ij(theta) = > pc . dijc . mu_jc . Ai . mu_jc, 589 /__ 590 c=1 591 592 where: 593 - djci is the PCS constant for spin j, state c and experiment or alignment i, 594 - N is the total number of states or structures, 595 - pc is the weight or probability associated with state c, 596 - mu_jc is the unit vector corresponding to spin j and state c, 597 - Ai is the alignment tensor. 598 599 In the fixed and equal probability case, the equation is:: 600 601 _N_ 602 1 \ T 603 delta_ij(theta) = - > dijc . mu_jc . Ai . mu_jc, 604 N /__ 605 c=1 606 607 The PCS constant is defined as:: 608 609 mu0 15kT 1 610 dijc = --- ----- ---- , 611 4pi Bo**2 r**3 612 613 where: 614 - mu0 is the permeability of free space, 615 - k is Boltzmann's constant, 616 - T is the absolute temperature (different for each experiment), 617 - Bo is the magnetic field strength (different for each experiment), 618 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state). 619 620 621 Stored data structures 622 ====================== 623 624 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 PCSs and the alignment tensors. 625 626 Dij(theta) 627 ---------- 628 629 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}. 630 631 delta_ij(theta) 632 --------------- 633 634 The back calculated PCS. This is a rank-2 tensor with indices {i, j}. 635 636 Ai 637 -- 638 639 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}. 640 641 642 @param params: The vector of parameter values. 643 @type params: numpy rank-1 array 644 @return: The chi-squared or SSE value. 645 @rtype: float 646 """ 647 648 # Scaling. 649 if self.scaling_flag: 650 params = dot(params, self.scaling_matrix) 651 652 # Initial chi-squared (or SSE) value. 653 chi2_sum = 0.0 654 655 # Unpack both the probabilities (when the paramagnetic centre is also optimised). 656 if not self.probs_fixed and not self.centre_fixed: 657 # The probabilities. 658 self.probs = params[-(self.N-1)-3:-3] 659 660 # Unpack the probabilities (located at the end of the parameter array). 661 elif not self.probs_fixed: 662 self.probs = params[-(self.N-1):] 663 664 # Unpack the paramagnetic centre (also update the paramagnetic info). 665 if not self.centre_fixed: 666 self.paramag_centre = params[-3:] 667 self.paramag_info() 668 669 # Loop over each alignment. 670 index = 0 671 for align_index in range(self.num_align): 672 # Create tensor i from the parameters. 673 if not self.fixed_tensors[align_index]: 674 to_tensor(self.A[align_index], params[5*index:5*index + 5]) 675 index += 1 676 677 # The back calculated RDC. 678 if self.rdc_flag[align_index]: 679 # Loop over the spin pairs k. 680 for j in range(self.num_interatom): 681 # Calculate the average RDC. 682 if not self.missing_rdc[align_index, j]: 683 if self.rdc_pseudo_flags[j]: 684 self.rdc_theta[align_index, j] = ave_rdc_tensor_pseudoatom(self.dip_const[j], self.dip_vect[j], self.N, self.A[align_index], weights=self.probs) 685 else: 686 self.rdc_theta[align_index, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[align_index], weights=self.probs) 687 688 # Add the J coupling to convert into the back-calculated T = J+D value. 689 if self.T_flags[align_index, j]: 690 self.rdc_theta[align_index, j] += self.j_couplings[j] 691 692 # Take the absolute value. 693 if self.absolute_rdc[align_index, j]: 694 self.rdc_theta[align_index, j] = abs(self.rdc_theta[align_index, j]) 695 696 # The back calculated PCS. 697 if self.pcs_flag[align_index]: 698 # Loop over the spin systems j. 699 for j in range(self.num_spins): 700 # Calculate the average PCS. 701 if not self.missing_deltaij[align_index, j]: 702 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) 703 704 # Calculate and sum the single alignment chi-squared value (for the RDC). 705 if self.rdc_flag[align_index]: 706 chi2_sum = chi2_sum + chi2(self.rdc[align_index], self.rdc_theta[align_index], self.rdc_errors[align_index]) 707 708 # Calculate and sum the single alignment chi-squared value (for the PCS). 709 if self.pcs_flag[align_index]: 710 chi2_sum = chi2_sum + chi2(self.deltaij[align_index], self.deltaij_theta[align_index], self.pcs_errors[align_index]) 711 712 # Return the chi-squared value. 713 return chi2_sum

714 715

716 - def dfunc_standard(self, params):

717 """The gradient function for optimisation of the standard N-state model. 718 719 Description 720 =========== 721 722 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 or PCS 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). 723 724 725 Indices 726 ======= 727 728 For this calculation, six indices are looped over and used in the various data structures. These include: 729 - k, the index over all parameters, 730 - i, the index over alignments, 731 - j, the index over spin systems, 732 - c, the index over the N-states (or over the structures), 733 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 734 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 735 736 737 Equations 738 ========= 739 740 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. 741 742 743 The chi-squared gradient 744 ------------------------ 745 746 The equation is:: 747 ___ 748 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \ 749 ------------- = -2 > | ---------------- . ----------- | 750 dthetak /__ \ sigma_ij**2 dthetak / 751 ij 752 753 where: 754 - theta is the parameter vector, 755 - Dij are the measured RDCs or PCSs, 756 - Dij(theta) are the back calculated RDCs or PCSs, 757 - sigma_ij are the RDC or PCS errors, 758 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 759 760 761 The RDC gradient 762 ---------------- 763 764 This gradient is different for the various parameter types. 765 766 pc partial derivative 767 ~~~~~~~~~~~~~~~~~~~~~ 768 769 The population parameter partial derivative is:: 770 771 dDij(theta) T 772 ----------- = dj . mu_jc . Ai . mu_jc, 773 dpc 774 775 where: 776 - dj is the dipolar constant for spin j, 777 - mu_jc is the unit vector corresponding to spin j and state c, 778 - Ai is the alignment tensor. 779 780 Amn partial derivative 781 ~~~~~~~~~~~~~~~~~~~~~~ 782 783 The alignment tensor element partial derivative is:: 784 785 _N_ 786 dDij(theta) \ T dAi 787 ----------- = dj > pc . mu_jc . ---- . mu_jc, 788 dAmn /__ dAmn 789 c=1 790 791 where: 792 - dj is the dipolar constant for spin j, 793 - pc is the weight or probability associated with state c, 794 - mu_jc is the unit vector corresponding to spin j and state c, 795 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 796 797 In the case of fixed and equal populations, the equation is:: 798 799 _N_ 800 dDij(theta) dj \ T dAi 801 ----------- = -- > mu_jc . ---- . mu_jc, 802 dAmn N /__ dAmn 803 c=1 804 805 806 The PCS gradient 807 ---------------- 808 809 This gradient is also different for the various parameter types. 810 811 pc partial derivative 812 ~~~~~~~~~~~~~~~~~~~~~ 813 814 The population parameter partial derivative is:: 815 816 ddeltaij(theta) T 817 --------------- = dijc . mu_jc . Ai . mu_jc, 818 dpc 819 820 where: 821 - djc is the pseudocontact shift constant for spin j and state c, 822 - mu_jc is the unit vector corresponding to spin j and state c, 823 - Ai is the alignment tensor. 824 825 Amn partial derivative 826 ~~~~~~~~~~~~~~~~~~~~~~ 827 828 The alignment tensor element partial derivative is:: 829 830 _N_ 831 ddelta_ij(theta) \ T dAi 832 ---------------- = > pc . djc . mu_jc . ---- . mu_jc, 833 dAmn /__ dAmn 834 c=1 835 836 where: 837 - djc is the pseudocontact shift constant for spin j and state c, 838 - pc is the weight or probability associated with state c, 839 - mu_jc is the unit vector corresponding to spin j and state c, 840 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 841 842 In the case of fixed and equal populations, the equation is:: 843 844 _N_ 845 ddelta_ij(theta) 1 \ T dAi 846 ---------------- = - > djc . mu_jc . ---- . mu_jc, 847 dAmn N /__ dAmn 848 c=1 849 850 xi partial derivative 851 ~~~~~~~~~~~~~~~~~~~~~ 852 853 The paramagnetic position partial derivative is:: 854 855 _N_ 856 ddelta_ij(theta) \ / ddjc dr_jcT dr_jc \ 857 ---------------- = > pc . | ----.r_jcT.Ai.r_jc + djc.------.Ai.r_jc + djc.r_jcT.Ai.----- | , 858 dxi /__ \ dxi dxi dxi / 859 c=1 860 861 where xi are the paramagnetic position coordinates {x0, x1, x2} and the last two terms in the sum are equal due to the symmetry of the alignment tensor, and:: 862 863 ddjc mu0 15kT 5 (si - xi) 864 ---- = --- ----- --------------------------------------------- , 865 dxi 4pi Bo**2 ((sx-x0)**2 + (sy-x1)**2 + (sz-x2)**2)**(7/2) 866 867 and:: 868 869 dr | 1 | dr | 0 | dr | 0 | 870 -- = - | 0 | , -- = - | 1 | , -- = - | 0 | . 871 dx | 0 | dy | 0 | dy | 1 | 872 873 The pseudocontact shift constant is defined here as:: 874 875 mu0 15kT 1 876 djc = --- ----- ------ , 877 4pi Bo**2 rjc**5 878 879 880 The alignment tensor gradient 881 ----------------------------- 882 883 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are:: 884 885 dAi | 1 0 0 | 886 ---- = | 0 0 0 |, 887 dAxx | 0 0 -1 | 888 889 dAi | 0 0 0 | 890 ---- = | 0 1 0 |, 891 dAyy | 0 0 -1 | 892 893 dAi | 0 1 0 | 894 ---- = | 1 0 0 |, 895 dAxy | 0 0 0 | 896 897 dAi | 0 0 1 | 898 ---- = | 0 0 0 |, 899 dAxz | 1 0 0 | 900 901 dAi | 0 0 0 | 902 ---- = | 0 0 1 |. 903 dAyz | 0 1 0 | 904 905 As these are invariant, they can be pre-calculated. 906 907 908 Stored data structures 909 ====================== 910 911 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. 912 913 dDij(theta)/dthetak 914 ------------------- 915 916 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}. 917 918 ddeltaij(theta)/dthetak 919 ----------------------- 920 921 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}. 922 923 dAi/dAmn 924 -------- 925 926 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}. 927 928 929 @param params: The vector of parameter values. This is unused as it is assumed that func() was called first. 930 @type params: numpy rank-1 array 931 @return: The chi-squared or SSE gradient. 932 @rtype: numpy rank-1 array 933 """ 934 935 # Initial chi-squared (or SSE) gradient. 936 self.dchi2 = self.dchi2 * 0.0 937 938 # Loop over each alignment. 939 for align_index in range(self.num_align): 940 # Construct the Amn partial derivative components for the RDC. 941 if not self.fixed_tensors[align_index]: 942 for j in range(self.num_interatom): 943 if self.rdc_flag[align_index] and not self.missing_rdc[align_index, j]: 944 if self.rdc_pseudo_flags[j]: 945 self.drdc_theta[align_index*5, align_index, j] = ave_rdc_tensor_pseudoatom_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs) 946 self.drdc_theta[align_index*5+1, align_index, j] = ave_rdc_tensor_pseudoatom_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs) 947 self.drdc_theta[align_index*5+2, align_index, j] = ave_rdc_tensor_pseudoatom_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs) 948 self.drdc_theta[align_index*5+3, align_index, j] = ave_rdc_tensor_pseudoatom_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs) 949 self.drdc_theta[align_index*5+4, align_index, j] = ave_rdc_tensor_pseudoatom_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs) 950 else: 951 self.drdc_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) 952 self.drdc_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) 953 self.drdc_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) 954 self.drdc_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) 955 self.drdc_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) 956 957 # Add the J coupling to convert into the back-calculated T = J+D value. 958 if self.T_flags[align_index, j]: 959 raise RelaxError("Gradients for T = J+D data have not been implemented yet.") 960 961 # Construct the Amn partial derivative components for the PCS. 962 if not self.fixed_tensors[align_index]: 963 for j in range(self.num_spins): 964 if self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 965 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) 966 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) 967 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) 968 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) 969 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) 970 971 # Construct the pc partial derivative gradient components, looping over each state. 972 if not self.probs_fixed: 973 # Shift the parameter index if the paramagnetic position is optimised. 974 x = 0 975 if not self.centre_fixed: 976 x = 3 977 978 # Loop over each state. 979 for c in range(self.N - 1 - x): 980 # Index in the parameter array. 981 param_index = self.num_align_params + c 982 983 # Calculate the RDC for state c (this is the pc partial derivative). 984 for j in range(self.num_interatom): 985 if self.rdc_flag[align_index] and not self.missing_rdc[align_index, j]: 986 self.drdc_theta[param_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j][c], self.A[align_index]) 987 988 # Calculate the PCS for state c (this is the pc partial derivative). 989 for j in range(self.num_spins): 990 if self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 991 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]) 992 993 # Construct the paramagnetic centre c partial derivative components for the PCS. 994 if not self.centre_fixed: 995 for j in range(self.num_spins): 996 if self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 997 self.ddeltaij_theta[-3, align_index, j] = ave_pcs_tensor_ddeltaij_dc(ddj=self.dpcs_const_theta[align_index, j,:, 0], dj=self.pcs_const[align_index, j], r=self.paramag_dist[j], unit_vect=self.paramag_unit_vect[j], N=self.N, Ai=self.A[align_index], dr_dc=self.dr_theta[0], weights=self.probs) 998 self.ddeltaij_theta[-2, align_index, j] = ave_pcs_tensor_ddeltaij_dc(ddj=self.dpcs_const_theta[align_index, j,:, 1], dj=self.pcs_const[align_index, j], r=self.paramag_dist[j], unit_vect=self.paramag_unit_vect[j], N=self.N, Ai=self.A[align_index], dr_dc=self.dr_theta[1], weights=self.probs) 999 self.ddeltaij_theta[-1, align_index, j] = ave_pcs_tensor_ddeltaij_dc(ddj=self.dpcs_const_theta[align_index, j,:, 2], dj=self.pcs_const[align_index, j], r=self.paramag_dist[j], unit_vect=self.paramag_unit_vect[j], N=self.N, Ai=self.A[align_index], dr_dc=self.dr_theta[2], weights=self.probs) 1000 1001 # Construct the chi-squared gradient element for parameter k, alignment i. 1002 for k in range(self.total_num_params): 1003 # RDC part of the chi-squared gradient. 1004 if self.rdc_flag[align_index]: 1005 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.rdc[align_index], self.rdc_theta[align_index], self.drdc_theta[k, align_index], self.rdc_errors[align_index]) 1006 1007 # PCS part of the chi-squared gradient. 1008 if self.pcs_flag[align_index]: 1009 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_errors[align_index]) 1010 1011 # Diagonal scaling. 1012 if self.scaling_flag: 1013 self.dchi2 = dot(self.dchi2, self.scaling_matrix) 1014 1015 # Return a copy of the gradient. 1016 return self.dchi2 * 1.0

1017 1018

1019 - def d2func_standard(self, params):

1020 """The Hessian function for optimisation of the standard N-state model. 1021 1022 Description 1023 =========== 1024 1025 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). 1026 1027 1028 Indices 1029 ======= 1030 1031 For this calculation, six indices are looped over and used in the various data structures. These include: 1032 - k, the index over all parameters, 1033 - i, the index over alignments, 1034 - j, the index over spin systems, 1035 - c, the index over the N-states (or over the structures), 1036 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1037 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1038 1039 1040 Equations 1041 ========= 1042 1043 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. 1044 1045 1046 The chi-squared Hessian 1047 ----------------------- 1048 1049 The equation is:: 1050 ___ 1051 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \ 1052 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |. 1053 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 1054 ij 1055 1056 where: 1057 - theta is the parameter vector, 1058 - Dij are the measured RDCs or PCSs, 1059 - Dij(theta) are the back calculated RDCs or PCSs, 1060 - sigma_ij are the RDC or PCS errors, 1061 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1062 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k. 1063 1064 1065 The RDC Hessian 1066 --------------- 1067 1068 pc-pd second partial derivatives 1069 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1070 1071 The probability parameter second partial derivative is:: 1072 1073 d2Dij(theta) 1074 ------------ = 0. 1075 dpc.dpd 1076 1077 1078 pc-Anm second partial derivatives 1079 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1080 1081 The probability parameter-tensor element second partial derivative is:: 1082 1083 d2Dij(theta) T dAi 1084 ------------ = dj . mu_jc . ---- . mu_jc. 1085 dpc.dAmn dAmn 1086 1087 1088 Amn-Aop second partial derivatives 1089 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1090 1091 The alignment tensor element second partial derivative is:: 1092 1093 d2Dij(theta) 1094 ------------ = 0. 1095 dAmn.dAop 1096 1097 1098 The PCS Hessian 1099 --------------- 1100 1101 pc-pd second partial derivatives 1102 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1103 1104 The probability parameter second partial derivative is:: 1105 1106 d2delta_ij(theta) 1107 ----------------- = 0. 1108 dpc.dpd 1109 1110 1111 pc-Anm second partial derivatives 1112 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1113 1114 The probability parameter-tensor element second partial derivative is:: 1115 1116 d2delta_ij(theta) T dAi 1117 ----------------- = djc . mu_jc . ---- . mu_jc. 1118 dpc.dAmn dAmn 1119 1120 1121 Amn-Aop second partial derivatives 1122 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1123 1124 The alignment tensor element second partial derivative is:: 1125 1126 d2delta_ij(theta) 1127 ----------------- = 0 1128 dAmn.dAop 1129 1130 1131 The alignment tensor Hessian 1132 ---------------------------- 1133 1134 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:: 1135 1136 d2Ai | 0 0 0 | 1137 --------- = | 0 0 0 |. 1138 dAmn.dAop | 0 0 0 | 1139 1140 1141 @param params: The vector of parameter values. This is unused as it is assumed that func() was called first. 1142 @type params: numpy rank-1 array 1143 @return: The chi-squared or SSE Hessian. 1144 @rtype: numpy rank-2 array 1145 """ 1146 1147 # Initial chi-squared (or SSE) Hessian. 1148 self.d2chi2 = self.d2chi2 * 0.0 1149 1150 # Loop over each alignment. 1151 for align_index in range(self.num_align): 1152 # Construct the pc-Amn second partial derivative Hessian components. 1153 if not self.probs_fixed: 1154 for c in range(self.N - 1): 1155 # Index in the parameter array. 1156 pc_index = self.num_align_params + c 1157 1158 # Calculate the RDC Hessian component. 1159 for j in range(self.num_interatom): 1160 if self.fixed_tensors[align_index] and self.rdc_flag[align_index] and not self.missing_rdc[align_index, j]: 1161 self.d2rdc_theta[pc_index, align_index*5+0, align_index, j] = self.d2rdc_theta[align_index*5+0, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j][c], self.dA[0]) 1162 self.d2rdc_theta[pc_index, align_index*5+1, align_index, j] = self.d2rdc_theta[align_index*5+1, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j][c], self.dA[1]) 1163 self.d2rdc_theta[pc_index, align_index*5+2, align_index, j] = self.d2rdc_theta[align_index*5+2, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j][c], self.dA[2]) 1164 self.d2rdc_theta[pc_index, align_index*5+3, align_index, j] = self.d2rdc_theta[align_index*5+3, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j][c], self.dA[3]) 1165 self.d2rdc_theta[pc_index, align_index*5+4, align_index, j] = self.d2rdc_theta[align_index*5+4, pc_index, align_index, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j][c], self.dA[4]) 1166 1167 # Calculate the PCS Hessian component. 1168 for j in range(self.num_spins): 1169 if self.fixed_tensors[align_index] and self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 1170 self.d2deltaij_theta[pc_index, align_index*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]) 1171 self.d2deltaij_theta[pc_index, align_index*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]) 1172 self.d2deltaij_theta[pc_index, align_index*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]) 1173 self.d2deltaij_theta[pc_index, align_index*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]) 1174 self.d2deltaij_theta[pc_index, align_index*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]) 1175 1176 # Construct the paramagnetic centre c partial derivative components for the PCS. 1177 if not self.centre_fixed: 1178 raise RelaxError("The Hessian equations for optimising the paramagnetic centre position are not yet implemented.") 1179 1180 # Construct the chi-squared Hessian element for parameters j and k, alignment i. 1181 for j in range(self.total_num_params): 1182 for k in range(self.total_num_params): 1183 # RDC part of the chi-squared gradient. 1184 if self.rdc_flag[align_index]: 1185 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.rdc[align_index], self.rdc_theta[align_index], self.drdc_theta[j, align_index], self.drdc_theta[k, align_index], self.d2rdc_theta[j, k, align_index], self.rdc_errors[align_index]) 1186 1187 # PCS part of the chi-squared gradient. 1188 if self.pcs_flag[align_index]: 1189 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_errors[align_index]) 1190 1191 # Diagonal scaling. 1192 if self.scaling_flag: 1193 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix) 1194 1195 # Return a copy of the Hessian. 1196 return self.d2chi2 * 1.0

1197 1198

1199 - def paramag_info(self):

1200 """Calculate the paramagnetic centre to spin vectors, distances and constants.""" 1201 1202 # Get the vectors and distances. 1203 if rank(self.paramag_centre) == 1: 1204 vectors_single_centre(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1205 else: 1206 vectors_centre_per_state(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1207 1208 # The PCS constants. 1209 for align_index in range(self.num_align): 1210 for j in range(self.num_spins): 1211 for c in range(self.N): 1212 self.pcs_const[align_index, j, c] = pcs_constant(self.temp[align_index], self.frq[align_index], self.paramag_dist[j, c]) 1213 1214 # The PCS constant gradient components. 1215 if not self.centre_fixed: 1216 pcs_constant_grad(T=self.temp[align_index], Bo=self.frq[align_index], r=self.paramag_dist[j, c], unit_vect=self.paramag_unit_vect[j, c], grad=self.dpcs_const_theta[align_index, j, c])

Source Code for Module target_functions.n_state_model