1   
  2   
  3   
  4   
  5   
  6   
  7   
  8   
  9   
 10   
 11   
 12   
 13   
 14   
 15   
 16   
 17   
 18   
 19   
 20   
 21   
 22   
 23   
 24  from math import sqrt 
 25  from numpy import outer 
 26   
 27   
 28   
 29   
 30   
 31   
 32   
 33   
 34   
 35   
 37      """Weight for spherical diffusion. 
 38   
 39      The equation is:: 
 40   
 41          c0 = 1. 
 42      """ 
 43   
 44      data.ci[0] = 1.0 
  45   
 46   
 47   
 48   
 49   
 50   
 51   
 52   
 53   
 54   
 55   
 56   
 58      """Weights for spheroidal diffusion. 
 59   
 60      The equations are:: 
 61   
 62          c-1 = 1/4 (3dz**2 - 1)**2, 
 63   
 64          c0  = 3dz**2 (1 - dz**2), 
 65   
 66          c1  = 3/4 (dz**2 - 1)**2, 
 67   
 68      where dz is the direction cosine of the unit bond vector along the z-axis of the diffusion 
 69      tensor which is calculated as the dot product of the unit bond vector with a unit vector along 
 70      Dpar. 
 71      """ 
 72   
 73       
 74      data.three_dz2_one = 3.0 * data.dz**2 - 1.0 
 75      data.one_two_dz2 = 1.0 - 2.0 * data.dz**2 
 76      data.one_dz2 = 1.0 - data.dz**2 
 77      data.dz2_one = -data.one_dz2 
 78   
 79       
 80      data.ci[0] = 0.25 * data.three_dz2_one**2 
 81      data.ci[1] = 3.0 * data.dz**2 * data.one_dz2 
 82      data.ci[2] = 0.75 * data.dz2_one**2 
  83   
 84   
 85   
 86   
 87   
 88   
 90      """Weight gradient for spheroidal diffusion. 
 91   
 92      The equations are:: 
 93   
 94          dc-1                     ddz 
 95          ---- =  3dz (3dz**2 - 1) --- , 
 96          dOi                      dOi 
 97   
 98          dc0                      ddz 
 99          ---  =  6dz (1 - 2dz**2) --- , 
100          dOi                      dOi 
101   
102          dc1                     ddz 
103          ---  =  3dz (dz**2 - 1) --- , 
104          dOi                     dOi 
105   
106      where the orientation parameter set O is {theta, phi}. 
107      """ 
108   
109       
110      data.dci[2:, 0] = 3.0 * data.dz * data.three_dz2_one * data.ddz_dO 
111      data.dci[2:, 1] = 6.0 * data.dz * data.one_two_dz2 * data.ddz_dO 
112      data.dci[2:, 2] = 3.0 * data.dz * data.dz2_one * data.ddz_dO 
 113   
114   
115   
116   
117   
118   
120      """Weight Hessian for spheroidal diffusion. 
121   
122      The equations are:: 
123   
124           d2c-1        /             ddz   ddz                      d2dz   \  
125          -------  =  3 |(9dz**2 - 1) --- . ---  +  dz (3dz**2 - 1) ------- | , 
126          dOi.dOj       \             dOi   dOj                     dOi.dOj / 
127   
128           d2c0         /             ddz   ddz                      d2dz   \  
129          -------  =  6 |(1 - 6dz**2) --- . ---  +  dz (1 - 2dz**2) ------- | , 
130          dOi.dOj       \             dOi   dOj                     dOi.dOj / 
131   
132           d2c1         /             ddz   ddz                     d2dz   \  
133          -------  =  3 |(3dz**2 - 1) --- . ---  +  dz (dz**2 - 1) ------- | , 
134          dOi.dOj       \             dOi   dOj                    dOi.dOj / 
135   
136      where the orientation parameter set O is {theta, phi}. 
137      """ 
138   
139       
140      op = outer(data.ddz_dO, data.ddz_dO) 
141   
142       
143      data.d2ci[2:, 2:, 0] = 3.0 * ((9.0 * data.dz**2 - 1.0) * op  +  data.dz * data.three_dz2_one * data.d2dz_dO2) 
144      data.d2ci[2:, 2:, 1] = 6.0 * ((1.0 - 6.0*data.dz**2) * op  +  data.dz * data.one_two_dz2 * data.d2dz_dO2) 
145      data.d2ci[2:, 2:, 2] = 3.0 * (data.three_dz2_one * op  +  data.dz * data.dz2_one * data.d2dz_dO2) 
 146   
147   
148   
149   
150   
151   
152   
153   
154   
155   
156   
157   
159      """Weight equations for ellipsoidal diffusion. 
160   
161      The equations are:: 
162   
163          c-2 = 1/4 (d - e), 
164   
165          c-1 = 3dy**2.dz**2, 
166   
167          c0  = 3dx**2.dz**2, 
168   
169          c1  = 3dx**2.dy**2, 
170   
171          c2  = 1/4 (d + e), 
172   
173      where:: 
174   
175          d  = 3(dx**4 + dy**4 + dz**4) - 1, 
176   
177          e  =  1/R [(1 + 3Dr)(dx**4 + 2dy**2.dz**2) + (1 - 3Dr)(dy**4 + 2dx**2.dz**2) 
178                     - 2(dz**4 + 2dx**2.dy**2)], 
179   
180      and where the factor R is defined as:: 
181               ___________ 
182          R = V 1 + 3Dr**2. 
183   
184      dx, dy, and dz are the direction cosines of the XH bond vector along the x, y, and z-axes of the 
185      diffusion tensor, calculated as the dot product of the unit bond vector and the unit vectors 
186      along Dx, Dy, and Dz respectively. 
187      """ 
188   
189       
190      data.R = sqrt(1 + 3.0*diff_data.params[2]**2) 
191      data.inv_R = 1.0 / data.R 
192   
193       
194      data.one_3Dr = 1.0 + 3.0 * diff_data.params[2] 
195      data.one_m3Dr = 1.0 - 3.0 * diff_data.params[2] 
196      data.dx_sqrd = data.dx**2 
197      data.dy_sqrd = data.dy**2 
198      data.dz_sqrd = data.dz**2 
199      data.dx_cubed = data.dx**3 
200      data.dy_cubed = data.dy**3 
201      data.dz_cubed = data.dz**3 
202      data.dx_quar = data.dx**4 
203      data.dy_quar = data.dy**4 
204      data.dz_quar = data.dz**4 
205   
206       
207      data.ex = data.dx_quar + 2.0 * data.dy_sqrd * data.dz_sqrd 
208      data.ey = data.dy_quar + 2.0 * data.dx_sqrd * data.dz_sqrd 
209      data.ez = data.dz_quar + 2.0 * data.dx_sqrd * data.dy_sqrd 
210   
211       
212      d = 3.0 * (data.dx_quar + data.dy_quar + data.dz_quar) - 1.0 
213   
214       
215      e = data.inv_R * (data.one_3Dr * data.ex  +  data.one_m3Dr * data.ey  -  2.0 * data.ez) 
216   
217       
218      data.ci[0] = 0.25 * (d - e) 
219   
220       
221      data.ci[1] = 3.0 * data.dy**2 * data.dz**2 
222   
223       
224      data.ci[2] = 3.0 * data.dx**2 * data.dz**2 
225   
226       
227      data.ci[3] = 3.0 * data.dx**2 * data.dy**2 
228   
229       
230      data.ci[4] = 0.25 * (d + e) 
 231   
232   
233   
234   
235   
236   
238      """Weight gradient for ellipsoidal diffusion. 
239   
240      Oi partial derivatives 
241      ====================== 
242   
243      The equations are:: 
244   
245          dc-2       /       ddx           ddy           ddz \     de 
246          ----  =  3 | dx**3 ---  +  dy**3 ---  +  dz**3 --- |  -  --- , 
247          dOi        \       dOi           dOi           dOi /     dOi 
248   
249   
250          dc-1            /    ddz        ddy \  
251          ----  =  6dy.dz | dy ---  +  dz --- | , 
252          dOi             \    dOi        dOi / 
253   
254   
255          dc0            /    ddz        ddx \  
256          ---  =  6dx.dz | dx ---  +  dz --- | , 
257          dOi            \    dOi        dOi / 
258   
259   
260          dc1            /    ddy        ddx \  
261          ---  =  6dx.dy | dx ---  +  dy --- | , 
262          dOi            \    dOi        dOi / 
263   
264   
265          dc2       /       ddx           ddy           ddz \     de 
266          ---  =  3 | dx**3 ---  +  dy**3 ---  +  dz**3 --- |  +  --- , 
267          dOi       \       dOi           dOi           dOi /     dOi 
268   
269   
270      where:: 
271   
272          de      1 /           /      ddx           /    ddz        ddy \ \  
273          ---  =  - | (1 + 3Dr) |dx**3 ---  +  dy.dz | dy ---  +  dz --- | | 
274          dOi     R \           \      dOi           \    dOi        dOi / / 
275   
276                                /       ddy           /    ddz        ddx \ \  
277                    + (1 - 3Dr) | dy**3 ---  +  dx.dz | dx ---  +  dz --- | | 
278                                \       dOi           \    dOi        dOi / / 
279   
280                        /       ddz           /    ddy        ddx \ \ \  
281                    - 2 | dz**3 ---  +  dx.dy | dx ---  +  dy --- | | | , 
282                        \       dOi           \    dOi        dOi / / / 
283   
284   
285      and where the orietation parameter set O is:: 
286   
287          O = {alpha, beta, gamma}. 
288   
289   
290      tm partial derivatives 
291      ====================== 
292   
293      The equations are:: 
294   
295          dc-2 
296          ----  =  0, 
297          dtm 
298   
299          dc-1 
300          ----  =  0, 
301          dtm 
302   
303          dc0 
304          ---   =  0, 
305          dtm 
306   
307          dc1 
308          ---   =  0, 
309          dtm 
310   
311          dc2 
312          ---   =  0. 
313          dtm 
314   
315   
316      Da partial derivatives 
317      ====================== 
318   
319      The equations are:: 
320   
321          dc-2 
322          ----  =  0, 
323          dDa 
324   
325          dc-1 
326          ----  =  0, 
327          dDa 
328   
329          dc0 
330          ---   =  0, 
331          dDa 
332   
333          dc1 
334          ---   =  0, 
335          dDa 
336   
337          dc2 
338          ---   =  0. 
339          dDa 
340   
341   
342   
343      Dr partial derivatives 
344      ====================== 
345   
346      The equations are:: 
347   
348          dc-2       3 de 
349          ----  =  - - ---, 
350          dDr        4 dDr 
351   
352          dc-1 
353          ----  =  0, 
354          dDr 
355   
356          dc0 
357          ---   =  0, 
358          dDr 
359   
360          dc1 
361          ---   =  0, 
362          dDr 
363   
364          dc2      3 de 
365          ---   =  - ---, 
366          dDr      4 dDr 
367   
368      where:: 
369   
370          de       1   /                                                                                                \  
371          ---  =  ---- | (1 - Dr) (dx**4 + 2dy**2.dz**2) - (1 + Dr) (dy**4 + 2dx**2.dz**2) + 2Dr (dz**4 + 2dx**2.dy**2) | . 
372          dDr     R**3 \                                                                                                / 
373   
374   
375      """ 
376   
377       
378      data.inv_R_cubed = data.inv_R ** 3 
379      data.one_Dr = 1.0 + diff_data.params[2] 
380      data.one_mDr = 1.0 - diff_data.params[2] 
381   
382   
383       
384       
385   
386       
387      data.ci_xy = data.dx * data.dy * (data.dx * data.ddy_dO  +  data.dy * data.ddx_dO) 
388      data.ci_xz = data.dx * data.dz * (data.dx * data.ddz_dO  +  data.dz * data.ddx_dO) 
389      data.ci_yz = data.dy * data.dz * (data.dy * data.ddz_dO  +  data.dz * data.ddy_dO) 
390   
391      data.ci_x = data.dx_cubed * data.ddx_dO 
392      data.ci_y = data.dy_cubed * data.ddy_dO 
393      data.ci_z = data.dz_cubed * data.ddz_dO 
394   
395      data.ci_X = data.ci_x + data.ci_yz 
396      data.ci_Y = data.ci_y + data.ci_xz 
397      data.ci_Z = data.ci_z + data.ci_xy 
398   
399       
400      dd_dOi = 3.0 * (data.ci_x + data.ci_y + data.ci_z) 
401   
402       
403      de_dOi = data.inv_R * (data.one_3Dr * data.ci_X  +  data.one_m3Dr * data.ci_Y  -  2.0 * data.ci_Z) 
404   
405       
406      data.dci[3:, 0] = dd_dOi - de_dOi 
407   
408       
409      data.dci[3:, 1] = 6.0 * data.ci_yz 
410   
411       
412      data.dci[3:, 2] = 6.0 * data.ci_xz 
413   
414       
415      data.dci[3:, 3] = 6.0 * data.ci_xy 
416   
417       
418      data.dci[3:, 4] = dd_dOi + de_dOi 
419   
420   
421       
422       
423   
424       
425      de_dDr = data.inv_R_cubed * (data.one_mDr * data.ex  -  data.one_Dr * data.ey  +  2.0 * diff_data.params[2] * data.ez) 
426   
427       
428      data.dci[2, 0] = -0.75 * de_dDr 
429   
430       
431      data.dci[2, 4] = 0.75 * de_dDr 
 432   
433   
434   
435   
436   
437   
439      """Weight Hessian for ellipsoidal diffusion. 
440   
441      Oi-Oj partial derivatives 
442      ========================= 
443   
444      The equations are:: 
445   
446           d2c-2        /       /     d2dx       ddx   ddx \           /     d2dy       ddy   ddy \           /     d2dz       ddz   ddz \ \       d2e 
447          -------  =  3 | dx**2 | dx ------- + 3 --- . --- |  +  dy**2 | dy ------- + 3 --- . --- |  +  dz**2 | dz ------- + 3 --- . --- | |  -  ------- , 
448          dOi.dOj       \       \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj / /     dOi.dOj 
449   
450   
451           d2c-1              /     d2dz       ddz   ddz \              / ddy   ddz     ddz   ddy \             /     d2dy       ddy   ddy \  
452          -------  =  6 dy**2 | dz -------  +  --- . --- |  +  12 dy.dz | --- . ---  +  --- . --- |  +  6 dz**2 | dy -------  +  --- . --- | , 
453          dOi.dOj             \    dOi.dOj     dOi   dOj /              \ dOi   dOj     dOi   dOj /             \    dOi.dOj     dOi   dOj / 
454   
455   
456           d2c0               /     d2dz       ddz   ddz \              / ddx   ddz     ddz   ddx \             /     d2dx       ddx   ddx \  
457          -------  =  6 dx**2 | dz -------  +  --- . --- |  +  12 dx.dz | --- . ---  +  --- . --- |  +  6 dz**2 | dx -------  +  --- . --- | , 
458          dOi.dOj             \    dOi.dOj     dOi   dOj /              \ dOi   dOj     dOi   dOj /             \    dOi.dOj     dOi   dOj / 
459   
460   
461           d2c1               /     d2dy       ddy   ddy \              / ddx   ddy     ddy   ddx \             /     d2dx       ddx   ddx \  
462          -------  =  6 dx**2 | dy -------  +  --- . --- |  +  12 dx.dy | --- . ---  +  --- . --- |  +  6 dy**2 | dx -------  +  --- . --- | , 
463          dOi.dOj             \    dOi.dOj     dOi   dOj /              \ dOi   dOj     dOi   dOj /             \    dOi.dOj     dOi   dOj / 
464   
465   
466           d2c2         /       /     d2dx       ddx   ddx \           /     d2dy       ddy   ddy \           /     d2dz       ddz   ddz \ \       d2e 
467          -------  =  3 | dx**2 | dx ------- + 3 --- . --- |  +  dy**2 | dy ------- + 3 --- . --- |  +  dz**2 | dz ------- + 3 --- . --- | |  +  ------- , 
468          dOi.dOj       \       \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj / /     dOi.dOj 
469   
470      where:: 
471   
472            d2e       1 /           /       /     d2dx       ddx   ddx \           /     d2dz       ddz   ddz \  
473          -------  =  - | (1 + 3Dr) | dx**2 | dx ------- + 3 --- . --- |  +  dy**2 | dz -------  +  --- . --- | 
474          dOi.dOj     R \           \       \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj /  
475   
476                                            /     d2dy       ddy   ddy \            / ddy   ddz     ddz   ddy \ \  
477                                    + dz**2 | dy -------  +  --- . --- |  +  2dy.dz | --- . ---  +  --- . --- | | 
478                                            \    dOi.dOj     dOi   dOj /            \ dOi   dOj     dOi   dOj / / 
479   
480                                    /       /     d2dy       ddy   ddy \           /     d2dz       ddz   ddz \  
481                        + (1 - 3Dr) | dy**2 | dy ------- + 3 --- . --- |  +  dx**2 | dz -------  +  --- . --- | 
482                                    \       \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj /  
483   
484                                            /     d2dx       ddx   ddx \            / ddx   ddz     ddz   ddx \ \  
485                                    + dz**2 | dx -------  +  --- . --- |  +  2dx.dz | --- . ---  +  --- . --- | | 
486                                            \    dOi.dOj     dOi   dOj /            \ dOi   dOj     dOi   dOj / / 
487   
488                            /       /     d2dz       ddz   ddz \           /     d2dy       ddy   ddy \  
489                        - 2 | dz**2 | dz ------- + 3 --- . --- |  +  dx**2 | dy -------  +  --- . --- | 
490                            \       \    dOi.dOj     dOi   dOj /           \    dOi.dOj     dOi   dOj /  
491   
492                                    /     d2dx       ddx   ddx \            / ddx   ddy     ddy   ddx \ \ \  
493                            + dy**2 | dx -------  +  --- . --- |  +  2dx.dy | --- . ---  +  --- . --- | | | 
494                                    \    dOi.dOj     dOi   dOj /            \ dOi   dOj     dOi   dOj / / / 
495   
496   
497      Oi-tm partial derivatives 
498      ========================= 
499   
500      The equations are:: 
501   
502           d2c-2 
503          -------  =  0, 
504          dOi.dtm 
505   
506   
507           d2c-1 
508          -------  =  0, 
509          dOi.dtm 
510   
511   
512           d2c0 
513          -------  =  0, 
514          dOi.dtm 
515   
516   
517           d2c1 
518          -------  =  0, 
519          dOi.dtm 
520   
521   
522           d2c2 
523          -------  =  0. 
524          dOi.dtm 
525   
526   
527      Oi-Da partial derivatives 
528      ========================= 
529   
530      The equations are:: 
531   
532           d2c-2 
533          -------  =  0, 
534          dOi.dDa 
535   
536   
537           d2c-1 
538          -------  =  0, 
539          dOi.dDa 
540   
541   
542           d2c0 
543          -------  =  0, 
544          dOi.dDa 
545   
546   
547           d2c1 
548          -------  =  0, 
549          dOi.dDa 
550   
551   
552           d2c2 
553          -------  =  0. 
554          dOi.dDa 
555   
556   
557      Oi-Dr partial derivatives 
558      ========================= 
559   
560      The equations are:: 
561   
562           d2c-2            d2e 
563          -------  =  - 3 -------, 
564          dOi.dDr         dOi.dDr 
565   
566   
567           d2c-1 
568          -------  =  0, 
569          dOi.dDr 
570   
571   
572           d2c0 
573          -------  =  0, 
574          dOi.dDr 
575   
576   
577           d2c1 
578          -------  =  0, 
579          dOi.dDr 
580   
581   
582           d2c2           d2e 
583          -------  =  3 -------, 
584          dOi.dDr       dOi.dDr 
585   
586      where:: 
587   
588           d2e         1   /          /       ddx           /    ddz        ddy \ \  
589          -------  =  ---- | (1 - Dr) | dx**3 ---  +  dy.dz | dy ---  +  dz --- | | 
590          dOi.dDr     R**3 \          \       dOi           \    dOi        dOi / / 
591   
592                                      /       ddy           /    ddz        ddx \ \  
593                           - (1 + Dr) | dy**3 ---  +  dx.dz | dx ---  +  dz --- | | 
594                                      \       dOi           \    dOi        dOi / / 
595   
596                                 /       ddz           /    ddy        ddx \ \ \  
597                           + 2Dr | dz**3 ---  +  dx.dy | dx ---  +  dy --- | | | 
598                                 \       dOi           \    dOi        dOi / / / 
599   
600   
601      tm-tm partial derivatives 
602      ========================= 
603   
604      The equations are:: 
605   
606          d2c-2 
607          -----  =  0, 
608          dtm2 
609   
610   
611          d2c-1 
612          -----  =  0, 
613          dtm2 
614   
615   
616          d2c0 
617          ----   =  0, 
618          dtm2 
619   
620   
621          d2c1 
622          ----   =  0, 
623          dtm2 
624   
625   
626          d2c2 
627          ----   =  0. 
628          dtm2 
629   
630   
631      tm-Da partial derivatives 
632      ========================= 
633   
634      The equations are:: 
635   
636           d2c-2 
637          -------  =  0, 
638          dtm.dDa 
639   
640   
641           d2c-1 
642          -------  =  0, 
643          dtm.dDa 
644   
645   
646           d2c0 
647          -------  =  0, 
648          dtm.dDa 
649   
650   
651           d2c1 
652          -------  =  0, 
653          dtm.dDa 
654   
655   
656           d2c2 
657          -------  =  0. 
658          dtm.dDa 
659   
660   
661      tm-Dr partial derivatives 
662      ========================= 
663   
664      The equations are:: 
665   
666           d2c-2 
667          -------  =  0, 
668          dtm.dDr 
669   
670   
671           d2c-1 
672          -------  =  0, 
673          dtm.dDr 
674   
675   
676           d2c0 
677          -------  =  0, 
678          dtm.dDr 
679   
680   
681           d2c1 
682          -------  =  0, 
683          dtm.dDr 
684   
685   
686           d2c2 
687          -------  =  0. 
688          dtm.dDr 
689   
690   
691      Da-Da partial derivatives 
692      ========================= 
693   
694      The equations are:: 
695   
696          d2c-2 
697          ------  =  0, 
698          dDa**2 
699   
700   
701          d2c-1 
702          ------  =  0, 
703          dDa**2 
704   
705   
706           d2c0 
707          ------  =  0, 
708          dDa**2 
709   
710   
711           d2c1 
712          ------  =  0, 
713          dDa**2 
714   
715   
716           d2c2 
717          ------  =  0. 
718          dDa**2 
719   
720   
721      Da-Dr partial derivatives 
722      ========================= 
723   
724      The equations are:: 
725   
726           d2c-2 
727          -------  =  0, 
728          dDa.dDr 
729   
730   
731           d2c-1 
732          -------  =  0, 
733          dDa.dDr 
734   
735   
736           d2c0 
737          -------  =  0, 
738          dDa.dDr 
739   
740   
741           d2c1 
742          -------  =  0, 
743          dDa.dDr 
744   
745   
746           d2c2 
747          -------  =  0. 
748          dDa.dDr 
749   
750   
751      Dr-Dr partial derivatives 
752      ========================= 
753   
754      The equations are:: 
755   
756          d2c-2        3  d2e 
757          ------  =  - - ------, 
758          dDr**2       4 dDr**2 
759   
760   
761          d2c-1 
762          ------  =  0, 
763          dDr**2 
764   
765   
766           d2c0 
767          ------  =  0, 
768          dDr**2 
769   
770   
771           d2c1 
772          ------  =  0, 
773          dDr**2 
774   
775   
776           d2c2      3  d2e 
777          ------  =  - ------, 
778          dDr**2     4 dDr**2 
779   
780      where:: 
781   
782           d2e        1   /                                                                                                                           \  
783          ------  =  ---- | (6Dr**2 - 9Dr - 1)(dx**4 + 2dy**2.dz**2) + (6Dr**2 + 9Dr - 1)(dy**4 + 2dx**2.dz**2) - 2(6Dr**2 - 1)(ddz*4 + 2dx**2.dy**2) | 
784          dDr**2     R**5 \                                                                                                                           / 
785       """ 
786   
787       
788       
789   
790       
791      op_xx = outer(data.ddx_dO, data.ddx_dO) 
792      op_yy = outer(data.ddy_dO, data.ddy_dO) 
793      op_zz = outer(data.ddz_dO, data.ddz_dO) 
794   
795      op_xy = outer(data.ddx_dO, data.ddy_dO) 
796      op_yx = outer(data.ddy_dO, data.ddx_dO) 
797   
798      op_xz = outer(data.ddx_dO, data.ddz_dO) 
799      op_zx = outer(data.ddz_dO, data.ddx_dO) 
800   
801      op_yz = outer(data.ddy_dO, data.ddz_dO) 
802      op_zy = outer(data.ddz_dO, data.ddy_dO) 
803   
804       
805      x_comp = data.dx * data.d2dx_dO2 + op_xx 
806      y_comp = data.dy * data.d2dy_dO2 + op_yy 
807      z_comp = data.dz * data.d2dz_dO2 + op_zz 
808   
809      x3_comp = data.dx_sqrd * (data.dx * data.d2dx_dO2 + 3.0 * op_xx) 
810      y3_comp = data.dy_sqrd * (data.dy * data.d2dy_dO2 + 3.0 * op_yy) 
811      z3_comp = data.dz_sqrd * (data.dz * data.d2dz_dO2 + 3.0 * op_zz) 
812   
813      xy_comp = data.dx_sqrd * y_comp  +  2.0 * data.dx * data.dy * (op_xy + op_yx)  +  data.dy_sqrd * x_comp 
814      xz_comp = data.dx_sqrd * z_comp  +  2.0 * data.dx * data.dz * (op_xz + op_zx)  +  data.dz_sqrd * x_comp 
815      yz_comp = data.dy_sqrd * z_comp  +  2.0 * data.dy * data.dz * (op_yz + op_zy)  +  data.dz_sqrd * y_comp 
816   
817       
818      d2d_dOidOj = 3.0 * (x3_comp + y3_comp + z3_comp) 
819   
820       
821      d2e_dOidOj = data.inv_R * (data.one_3Dr * (x3_comp + yz_comp) + data.one_m3Dr * (y3_comp + xz_comp) - 2.0 * (z3_comp + xy_comp)) 
822   
823       
824      data.d2ci[3:, 3:, 0] = d2d_dOidOj - d2e_dOidOj 
825   
826       
827      data.d2ci[3:, 3:, 1] = 6.0 * yz_comp 
828   
829       
830      data.d2ci[3:, 3:, 2] = 6.0 * xz_comp 
831   
832       
833      data.d2ci[3:, 3:, 3] = 6.0 * xy_comp 
834   
835       
836      data.d2ci[3:, 3:, 4] = d2d_dOidOj + d2e_dOidOj 
837   
838   
839       
840       
841   
842       
843      d2e_dOidDr = data.inv_R_cubed * (data.one_mDr * data.ci_X  -  data.one_Dr * data.ci_Y  +  2.0 * diff_data.params[2] * data.ci_Z) 
844   
845       
846      data.d2ci[3:, 2, 0] = data.d2ci[2, 3:, 0] = -3.0 * d2e_dOidDr 
847   
848       
849      data.d2ci[3:, 2, 4] = data.d2ci[2, 3:, 4] = 3.0 * d2e_dOidDr 
850   
851   
852       
853       
854   
855       
856      d2e1_dDr2 = 6.0 * diff_data.params[2]**2  -  9.0 * diff_data.params[2]  -  1.0 
857      d2e2_dDr2 = 6.0 * diff_data.params[2]**2  +  9.0 * diff_data.params[2]  -  1.0 
858      d2e3_dDr2 = 6.0 * diff_data.params[2]**2  -  1.0 
859   
860       
861      d2e_dDr2 = data.inv_R**5 * (d2e1_dDr2 * data.ex  +  d2e2_dDr2 * data.ey  -  2.0 * d2e3_dDr2 * data.ez) 
862   
863       
864      data.d2ci[2, 2, 0] = -0.75 * d2e_dDr2 
865   
866       
867      data.d2ci[2, 2, 4] = 0.75 * d2e_dDr2 
 868