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