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 if pcs != None: 226 self.num_align = max(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 # The gradient structures. 370 if not self.centre_fixed: 371 self.dpcs_const_theta = zeros((self.num_align, self.num_spins, self.N, 3), float64) 372 self.dr_theta = -eye(3) 373 374 # Set up the paramagnetic info. 375 self.paramag_info() 376 377 # PCS function, gradient, and Hessian matrices. 378 self.deltaij_theta = zeros((self.num_align, self.num_spins), float64) 379 self.ddeltaij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 380 self.d2deltaij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 381 382 # RDC function, gradient, and Hessian matrices. 383 self.Dij_theta = zeros((self.num_align, self.num_interatom), float64) 384 self.dDij_theta = zeros((self.total_num_params, self.num_align, self.num_interatom), float64) 385 self.d2Dij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_interatom), float64) 386 387 # Set the target function, gradient, and Hessian. 388 self.func = self.func_standard 389 self.dfunc = self.dfunc_standard 390 self.d2func = self.d2func_standard 391 392 # Variable probabilities. 393 self.probs_fixed = True 394 if model == 'population': 395 self.probs_fixed = False 396 397 # Fixed probabilities. 398 if model == 'fixed': 399 # The zero Hessian. 400 self.zero_hessian_rdc = zeros(self.num_interatom, float64) 401 self.zero_hessian_pcs = zeros(self.num_spins, float64) 402 403 # The probability array. 404 if probs: 405 self.probs = probs 406 407 # All structures have initial equal probability. 408 else: 409 self.probs = ones(self.N, float64) / self.N

410 411

412 - def func_2domain(self, params):

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

468 469

470 - def func_standard(self, params):

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

689 690

691 - def dfunc_standard(self, params):

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

981 982

983 - def d2func_standard(self, params):

984 """The Hessian function for optimisation of the standard N-state model. 985 986 Description 987 =========== 988 989 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). 990 991 992 Indices 993 ======= 994 995 For this calculation, six indices are looped over and used in the various data structures. These include: 996 - k, the index over all parameters, 997 - i, the index over alignments, 998 - j, the index over spin systems, 999 - c, the index over the N-states (or over the structures), 1000 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1001 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1002 1003 1004 Equations 1005 ========= 1006 1007 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. 1008 1009 1010 The chi-squared Hessian 1011 ----------------------- 1012 1013 The equation is:: 1014 ___ 1015 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \ 1016 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |. 1017 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 1018 ij 1019 1020 where: 1021 - theta is the parameter vector, 1022 - Dij are the measured RDCs or PCSs, 1023 - Dij(theta) are the back calculated RDCs or PCSs, 1024 - sigma_ij are the RDC or PCS errors, 1025 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1026 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k. 1027 1028 1029 The RDC Hessian 1030 --------------- 1031 1032 pc-pd second partial derivatives 1033 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1034 1035 The probability parameter second partial derivative is:: 1036 1037 d2Dij(theta) 1038 ------------ = 0. 1039 dpc.dpd 1040 1041 1042 pc-Anm second partial derivatives 1043 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1044 1045 The probability parameter-tensor element second partial derivative is:: 1046 1047 d2Dij(theta) T dAi 1048 ------------ = dj . mu_jc . ---- . mu_jc. 1049 dpc.dAmn dAmn 1050 1051 1052 Amn-Aop second partial derivatives 1053 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1054 1055 The alignment tensor element second partial derivative is:: 1056 1057 d2Dij(theta) 1058 ------------ = 0. 1059 dAmn.dAop 1060 1061 1062 The PCS Hessian 1063 --------------- 1064 1065 pc-pd second partial derivatives 1066 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1067 1068 The probability parameter second partial derivative is:: 1069 1070 d2delta_ij(theta) 1071 ----------------- = 0. 1072 dpc.dpd 1073 1074 1075 pc-Anm second partial derivatives 1076 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1077 1078 The probability parameter-tensor element second partial derivative is:: 1079 1080 d2delta_ij(theta) T dAi 1081 ----------------- = djc . mu_jc . ---- . mu_jc. 1082 dpc.dAmn dAmn 1083 1084 1085 Amn-Aop second partial derivatives 1086 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1087 1088 The alignment tensor element second partial derivative is:: 1089 1090 d2delta_ij(theta) 1091 ----------------- = 0 1092 dAmn.dAop 1093 1094 1095 The alignment tensor Hessian 1096 ---------------------------- 1097 1098 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:: 1099 1100 d2Ai | 0 0 0 | 1101 --------- = | 0 0 0 |. 1102 dAmn.dAop | 0 0 0 | 1103 1104 1105 @param params: The vector of parameter values. This is unused as it is assumed that func() was called first. 1106 @type params: numpy rank-1 array 1107 @return: The chi-squared or SSE Hessian. 1108 @rtype: numpy rank-2 array 1109 """ 1110 1111 # Initial chi-squared (or SSE) Hessian. 1112 self.d2chi2 = self.d2chi2 * 0.0 1113 1114 # Loop over each alignment. 1115 for align_index in range(self.num_align): 1116 # Construct the pc-Amn second partial derivative Hessian components. 1117 if not self.probs_fixed: 1118 for c in range(self.N - 1): 1119 # Index in the parameter array. 1120 pc_index = self.num_align_params + c 1121 1122 # Calculate the RDC Hessian component. 1123 for j in range(self.num_interatom): 1124 if self.fixed_tensors[align_index] and self.rdc_flag[align_index] and not self.missing_Dij[align_index, j]: 1125 self.d2Dij_theta[pc_index, align_index*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]) 1126 self.d2Dij_theta[pc_index, align_index*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]) 1127 self.d2Dij_theta[pc_index, align_index*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]) 1128 self.d2Dij_theta[pc_index, align_index*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]) 1129 self.d2Dij_theta[pc_index, align_index*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]) 1130 1131 # Calculate the PCS Hessian component. 1132 for j in range(self.num_spins): 1133 if self.fixed_tensors[align_index] and self.pcs_flag[align_index] and not self.missing_deltaij[align_index, j]: 1134 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]) 1135 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]) 1136 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]) 1137 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]) 1138 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]) 1139 1140 # Construct the paramagnetic centre c partial derivative components for the PCS. 1141 if not self.centre_fixed: 1142 raise RelaxError("The Hessian equations for optimising the paramagnetic centre position are not yet implemented.") 1143 1144 # Construct the chi-squared Hessian element for parameters j and k, alignment i. 1145 for j in range(self.total_num_params): 1146 for k in range(self.total_num_params): 1147 # RDC part of the chi-squared gradient. 1148 if self.rdc_flag[align_index]: 1149 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]) 1150 1151 # PCS part of the chi-squared gradient. 1152 if self.pcs_flag[align_index]: 1153 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]) 1154 1155 # Diagonal scaling. 1156 if self.scaling_flag: 1157 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix) 1158 1159 # Return a copy of the Hessian. 1160 return self.d2chi2 * 1.0

1161 1162

1163 - def paramag_info(self):

1164 """Calculate the paramagnetic centre to spin vectors, distances and constants.""" 1165 1166 # Get the vectors and distances. 1167 if rank(self.paramag_centre) == 1: 1168 vectors_single_centre(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1169 else: 1170 vectors_centre_per_state(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1171 1172 # The PCS constants. 1173 for align_index in range(self.num_align): 1174 for j in range(self.num_spins): 1175 for c in range(self.N): 1176 self.pcs_const[align_index, j, c] = pcs_constant(self.temp[align_index], self.frq[align_index], self.paramag_dist[j, c]) 1177 1178 # The PCS constant gradient components. 1179 if not self.centre_fixed: 1180 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 maths_fns.n_state_model