1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 from math import sqrt
24 from numpy import array, dot, eye, float64, ones, transpose, zeros
25
26
27 from lib.alignment.alignment_tensor import dAi_dAxx, dAi_dAyy, dAi_dAxy, dAi_dAxz, dAi_dAyz, to_tensor
28 from lib.alignment.paramag_centre import vectors_single_centre, vectors_centre_per_state
29 from lib.alignment.pcs import ave_pcs_tensor, ave_pcs_tensor_ddeltaij_dAmn, ave_pcs_tensor_ddeltaij_dc, pcs_constant_grad, pcs_tensor
30 from lib.alignment.rdc import ave_rdc_tensor, ave_rdc_tensor_dDij_dAmn, ave_rdc_tensor_pseudoatom, ave_rdc_tensor_pseudoatom_dDij_dAmn, rdc_tensor
31 from lib.errors import RelaxError
32 from lib.float import isNaN
33 from lib.geometry.rotations import euler_to_R_zyz
34 from lib.physical_constants import pcs_constant
35 from target_functions.chi2 import chi2, dchi2_element, d2chi2_element
36
37
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
165 self.N = N
166 self.params = 1.0 * init_params
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
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
190 self.scaling_matrix = scaling_matrix
191 if self.scaling_matrix is not None:
192 self.scaling_flag = True
193 else:
194 self.scaling_flag = False
195
196
197 if model == '2-domain':
198
199 if red_data is None or not len(red_data):
200 raise RelaxError("The red_data argument " + repr(red_data) + " must be supplied.")
201 if red_errors is 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 is 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
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
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
219
220
221
222
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
229 self.func = self.func_2domain
230 self.dfunc = None
231 self.d2func = None
232
233
234 elif model in ['population', 'fixed']:
235
236 self.num_align = 0
237 if rdcs is not None:
238 self.num_align = len(rdcs)
239 if pcs is not None:
240 self.num_align = max(self.num_align, len(pcs))
241
242
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 is None or len(rdcs[align_index]) == 0:
247 self.rdc_flag[align_index] = False
248 if pcs is 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
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 is None or not len(atomic_pos)):
257 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.")
258
259
260 self.num_spins = 0
261 if self.pcs_flag_sum:
262 self.num_spins = len(pcs[0])
263
264
265 self.num_interatom = 0
266 if self.rdc_flag_sum:
267 self.num_interatom = len(rdcs[0])
268
269
270 self.A = zeros((self.num_align, 3, 3), float64)
271 self.dA = zeros((5, 3, 3), float64)
272
273
274 if full_tensors is not None:
275
276 self.full_tensors = array(full_tensors, float64)
277
278
279 self.num_align_params = 0
280 index = 0
281 for align_index in range(self.num_align):
282
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
288 if not fixed_tensors[align_index]:
289 self.num_align_params += 5
290
291
292 else:
293
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
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
311 self.pcs_errors = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64)
312
313
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
324 self.rdc_errors = ones((self.num_align, self.num_interatom), float64)
325
326
327 if self.rdc_flag_sum:
328 self.missing_rdc = zeros((self.num_align, self.num_interatom), float64)
329
330
331 if self.pcs_flag_sum:
332 self.missing_deltaij = zeros((self.num_align, self.num_spins), float64)
333
334
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
341 self.missing_rdc[align_index, j] = 1
342
343
344 self.rdc[align_index, j] = 0.0
345
346
347 self.rdc_errors[align_index, j] = 1.0
348
349
350 rdc_weights[align_index, j] = 1.0
351
352
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
360 self.missing_deltaij[align_index, j] = 1
361
362
363 self.deltaij[align_index, j] = 0.0
364
365
366 self.pcs_errors[align_index, j] = 1.0
367
368
369 pcs_weights[align_index, j] = 1.0
370
371
372 self.pcs_errors[align_index, j] = self.pcs_errors[align_index, j] / sqrt(pcs_weights[align_index, j])
373
374
375 if self.pcs_flag_sum:
376
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 is None:
381 self.paramag_centre = zeros(3, float64)
382
383
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
389 self.paramag_info()
390
391
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
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
402 self.func = self.func_standard
403 self.dfunc = self.dfunc_standard
404 self.d2func = self.d2func_standard
405
406
407 self.probs_fixed = True
408 if model == 'population':
409 self.probs_fixed = False
410
411
412 if model == 'fixed':
413
414 self.zero_hessian_rdc = zeros(self.num_interatom, float64)
415 self.zero_hessian_pcs = zeros(self.num_spins, float64)
416
417
418 if probs:
419 self.probs = probs
420
421
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
438 if self.scaling_flag:
439 params = dot(params, self.scaling_matrix)
440
441
442 self.red_bc = self.red_bc * 0.0
443
444
445 for c in range(self.N):
446
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
450 self.RT[c] = transpose(self.R[c])
451
452
453 if c < self.N-1:
454 pc = params[c]
455
456
457 else:
458 pc = 1.0
459 for c2 in range(self.N-1):
460 pc = pc - params[c2]
461
462
463 for align_index in range(self.num_tensors):
464
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
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
473 for align_index in range(self.num_tensors):
474 self.red_bc_vector[5*align_index] = self.red_bc[align_index, 0, 0]
475 self.red_bc_vector[5*align_index+1] = self.red_bc[align_index, 1, 1]
476 self.red_bc_vector[5*align_index+2] = self.red_bc[align_index, 0, 1]
477 self.red_bc_vector[5*align_index+3] = self.red_bc[align_index, 0, 2]
478 self.red_bc_vector[5*align_index+4] = self.red_bc[align_index, 1, 2]
479
480
481 return chi2(self.red_data, self.red_bc_vector, self.red_errors)
482
483
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
649 if self.scaling_flag:
650 params = dot(params, self.scaling_matrix)
651
652
653 chi2_sum = 0.0
654
655
656 if not self.probs_fixed and not self.centre_fixed:
657
658 self.probs = params[-(self.N-1)-3:-3]
659
660
661 elif not self.probs_fixed:
662 self.probs = params[-(self.N-1):]
663
664
665 if not self.centre_fixed:
666 self.paramag_centre = params[-3:]
667 self.paramag_info()
668
669
670 index = 0
671 for align_index in range(self.num_align):
672
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
678 if self.rdc_flag[align_index]:
679
680 for j in range(self.num_interatom):
681
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
689 if self.T_flags[align_index, j]:
690 self.rdc_theta[align_index, j] += self.j_couplings[j]
691
692
693 if self.absolute_rdc[align_index, j]:
694 self.rdc_theta[align_index, j] = abs(self.rdc_theta[align_index, j])
695
696
697 if self.pcs_flag[align_index]:
698
699 for j in range(self.num_spins):
700
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
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
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
713 return chi2_sum
714
715
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
936 self.dchi2 = self.dchi2 * 0.0
937
938
939 for align_index in range(self.num_align):
940
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
958 if self.T_flags[align_index, j]:
959 raise RelaxError("Gradients for T = J+D data have not been implemented yet.")
960
961
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
972 if not self.probs_fixed:
973
974 x = 0
975 if not self.centre_fixed:
976 x = 3
977
978
979 for c in range(self.N - 1 - x):
980
981 param_index = self.num_align_params + c
982
983
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
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
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
1002 for k in range(self.total_num_params):
1003
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
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
1012 if self.scaling_flag:
1013 self.dchi2 = dot(self.dchi2, self.scaling_matrix)
1014
1015
1016 return self.dchi2 * 1.0
1017
1018
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
1148 self.d2chi2 = self.d2chi2 * 0.0
1149
1150
1151 for align_index in range(self.num_align):
1152
1153 if not self.probs_fixed:
1154 for c in range(self.N - 1):
1155
1156 pc_index = self.num_align_params + c
1157
1158
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
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
1177 if not self.centre_fixed:
1178 raise RelaxError("The Hessian equations for optimising the paramagnetic centre position are not yet implemented.")
1179
1180
1181 for j in range(self.total_num_params):
1182 for k in range(self.total_num_params):
1183
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
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
1192 if self.scaling_flag:
1193 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix)
1194
1195
1196 return self.d2chi2 * 1.0
1197
1198
1200 """Calculate the paramagnetic centre to spin vectors, distances and constants."""
1201
1202
1203 if self.paramag_centre.ndim == 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
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
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])
1217