1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.numerics.multipole;
39
40 import static java.lang.Math.fma;
41 import static java.lang.String.format;
42 import static org.apache.commons.math3.util.FastMath.sqrt;
43
44
45
46
47
48
49
50
51
52
53
54
55 public class CoulombTensorGlobal extends MultipoleTensor {
56
57
58
59
60
61
62 public CoulombTensorGlobal(int order) {
63 super(COORDINATES.GLOBAL, order);
64 operator = OPERATOR.COULOMB;
65 }
66
67
68
69
70
71
72 @Override
73 public double getd2EdZ2() {
74 return 0.0;
75 }
76
77
78
79
80
81
82 @Override
83 public double getdEdZ() {
84 return 0.0;
85 }
86
87
88
89
90 @Override
91 public void setR(double dx, double dy, double dz) {
92 x = dx;
93 y = dy;
94 z = dz;
95 r2 = (x * x + y * y + z * z);
96 R = sqrt(r2);
97 }
98
99
100
101
102
103
104 protected void source(double[] T000) {
105
106
107 double ir = 1.0 / R;
108 double ir2 = ir * ir;
109 for (int n = 0; n < o1; n++) {
110 T000[n] = coulombSource[n] * ir;
111 ir *= ir2;
112 }
113 }
114
115
116
117
118 @Override
119 protected void order1() {
120 source(work);
121 double term0000 = work[0];
122 double term0001 = work[1];
123 R000 = term0000;
124 R100 = x * term0001;
125 R010 = y * term0001;
126 R001 = z * term0001;
127 }
128
129
130
131
132 @Override
133 protected void order2() {
134 source(work);
135 double term0000 = work[0];
136 double term0001 = work[1];
137 double term0002 = work[2];
138 R000 = term0000;
139 R100 = x * term0001;
140 double term1001 = x * term0002;
141 R200 = fma(x, term1001, term0001);
142 R010 = y * term0001;
143 double term0101 = y * term0002;
144 R020 = fma(y, term0101, term0001);
145 R110 = y * term1001;
146 R001 = z * term0001;
147 double term0011 = z * term0002;
148 R002 = fma(z, term0011, term0001);
149 R011 = z * term0101;
150 R101 = z * term1001;
151 }
152
153
154
155
156 @Override
157 protected void order3() {
158 source(work);
159 double term0000 = work[0];
160 double term0001 = work[1];
161 double term0002 = work[2];
162 double term0003 = work[3];
163 R000 = term0000;
164 R100 = x * term0001;
165 double term1001 = x * term0002;
166 R200 = fma(x, term1001, term0001);
167 double term1002 = x * term0003;
168 double term2001 = fma(x, term1002, term0002);
169 R300 = fma(x, term2001, 2 * term1001);
170 R010 = y * term0001;
171 double term0101 = y * term0002;
172 R020 = fma(y, term0101, term0001);
173 double term0102 = y * term0003;
174 double term0201 = fma(y, term0102, term0002);
175 R030 = fma(y, term0201, 2 * term0101);
176 R110 = y * term1001;
177 double term1101 = y * term1002;
178 R120 = fma(y, term1101, term1001);
179 R210 = y * term2001;
180 R001 = z * term0001;
181 double term0011 = z * term0002;
182 R002 = fma(z, term0011, term0001);
183 double term0012 = z * term0003;
184 double term0021 = fma(z, term0012, term0002);
185 R003 = fma(z, term0021, 2 * term0011);
186 R011 = z * term0101;
187 double term0111 = z * term0102;
188 R012 = fma(z, term0111, term0101);
189 R021 = z * term0201;
190 R101 = z * term1001;
191 double term1011 = z * term1002;
192 R102 = fma(z, term1011, term1001);
193 R111 = z * term1101;
194 R201 = z * term2001;
195 }
196
197
198
199
200 @Override
201 protected void order4() {
202 source(work);
203 double term0000 = work[0];
204 double term0001 = work[1];
205 double term0002 = work[2];
206 double term0003 = work[3];
207 double term0004 = work[4];
208 R000 = term0000;
209 R100 = x * term0001;
210 double term1001 = x * term0002;
211 R200 = fma(x, term1001, term0001);
212 double term1002 = x * term0003;
213 double term2001 = fma(x, term1002, term0002);
214 R300 = fma(x, term2001, 2 * term1001);
215 double term1003 = x * term0004;
216 double term2002 = fma(x, term1003, term0003);
217 double term3001 = fma(x, term2002, 2 * term1002);
218 R400 = fma(x, term3001, 3 * term2001);
219 R010 = y * term0001;
220 double term0101 = y * term0002;
221 R020 = fma(y, term0101, term0001);
222 double term0102 = y * term0003;
223 double term0201 = fma(y, term0102, term0002);
224 R030 = fma(y, term0201, 2 * term0101);
225 double term0103 = y * term0004;
226 double term0202 = fma(y, term0103, term0003);
227 double term0301 = fma(y, term0202, 2 * term0102);
228 R040 = fma(y, term0301, 3 * term0201);
229 R110 = y * term1001;
230 double term1101 = y * term1002;
231 R120 = fma(y, term1101, term1001);
232 double term1102 = y * term1003;
233 double term1201 = fma(y, term1102, term1002);
234 R130 = fma(y, term1201, 2 * term1101);
235 R210 = y * term2001;
236 double term2101 = y * term2002;
237 R220 = fma(y, term2101, term2001);
238 R310 = y * term3001;
239 R001 = z * term0001;
240 double term0011 = z * term0002;
241 R002 = fma(z, term0011, term0001);
242 double term0012 = z * term0003;
243 double term0021 = fma(z, term0012, term0002);
244 R003 = fma(z, term0021, 2 * term0011);
245 double term0013 = z * term0004;
246 double term0022 = fma(z, term0013, term0003);
247 double term0031 = fma(z, term0022, 2 * term0012);
248 R004 = fma(z, term0031, 3 * term0021);
249 R011 = z * term0101;
250 double term0111 = z * term0102;
251 R012 = fma(z, term0111, term0101);
252 double term0112 = z * term0103;
253 double term0121 = fma(z, term0112, term0102);
254 R013 = fma(z, term0121, 2 * term0111);
255 R021 = z * term0201;
256 double term0211 = z * term0202;
257 R022 = fma(z, term0211, term0201);
258 R031 = z * term0301;
259 R101 = z * term1001;
260 double term1011 = z * term1002;
261 R102 = fma(z, term1011, term1001);
262 double term1012 = z * term1003;
263 double term1021 = fma(z, term1012, term1002);
264 R103 = fma(z, term1021, 2 * term1011);
265 R111 = z * term1101;
266 double term1111 = z * term1102;
267 R112 = fma(z, term1111, term1101);
268 R121 = z * term1201;
269 R201 = z * term2001;
270 double term2011 = z * term2002;
271 R202 = fma(z, term2011, term2001);
272 R211 = z * term2101;
273 R301 = z * term3001;
274 }
275
276
277
278
279 @Override
280 protected void order5() {
281 source(work);
282 double term0000 = work[0];
283 double term0001 = work[1];
284 double term0002 = work[2];
285 double term0003 = work[3];
286 double term0004 = work[4];
287 double term0005 = work[5];
288 R000 = term0000;
289 R100 = x * term0001;
290 double term1001 = x * term0002;
291 R200 = fma(x, term1001, term0001);
292 double term1002 = x * term0003;
293 double term2001 = fma(x, term1002, term0002);
294 R300 = fma(x, term2001, 2 * term1001);
295 double term1003 = x * term0004;
296 double term2002 = fma(x, term1003, term0003);
297 double term3001 = fma(x, term2002, 2 * term1002);
298 R400 = fma(x, term3001, 3 * term2001);
299 double term1004 = x * term0005;
300 double term2003 = fma(x, term1004, term0004);
301 double term3002 = fma(x, term2003, 2 * term1003);
302 double term4001 = fma(x, term3002, 3 * term2002);
303 R500 = fma(x, term4001, 4 * term3001);
304 R010 = y * term0001;
305 double term0101 = y * term0002;
306 R020 = fma(y, term0101, term0001);
307 double term0102 = y * term0003;
308 double term0201 = fma(y, term0102, term0002);
309 R030 = fma(y, term0201, 2 * term0101);
310 double term0103 = y * term0004;
311 double term0202 = fma(y, term0103, term0003);
312 double term0301 = fma(y, term0202, 2 * term0102);
313 R040 = fma(y, term0301, 3 * term0201);
314 double term0104 = y * term0005;
315 double term0203 = fma(y, term0104, term0004);
316 double term0302 = fma(y, term0203, 2 * term0103);
317 double term0401 = fma(y, term0302, 3 * term0202);
318 R050 = fma(y, term0401, 4 * term0301);
319 R110 = y * term1001;
320 double term1101 = y * term1002;
321 R120 = fma(y, term1101, term1001);
322 double term1102 = y * term1003;
323 double term1201 = fma(y, term1102, term1002);
324 R130 = fma(y, term1201, 2 * term1101);
325 double term1103 = y * term1004;
326 double term1202 = fma(y, term1103, term1003);
327 double term1301 = fma(y, term1202, 2 * term1102);
328 R140 = fma(y, term1301, 3 * term1201);
329 R210 = y * term2001;
330 double term2101 = y * term2002;
331 R220 = fma(y, term2101, term2001);
332 double term2102 = y * term2003;
333 double term2201 = fma(y, term2102, term2002);
334 R230 = fma(y, term2201, 2 * term2101);
335 R310 = y * term3001;
336 double term3101 = y * term3002;
337 R320 = fma(y, term3101, term3001);
338 R410 = y * term4001;
339 R001 = z * term0001;
340 double term0011 = z * term0002;
341 R002 = fma(z, term0011, term0001);
342 double term0012 = z * term0003;
343 double term0021 = fma(z, term0012, term0002);
344 R003 = fma(z, term0021, 2 * term0011);
345 double term0013 = z * term0004;
346 double term0022 = fma(z, term0013, term0003);
347 double term0031 = fma(z, term0022, 2 * term0012);
348 R004 = fma(z, term0031, 3 * term0021);
349 double term0014 = z * term0005;
350 double term0023 = fma(z, term0014, term0004);
351 double term0032 = fma(z, term0023, 2 * term0013);
352 double term0041 = fma(z, term0032, 3 * term0022);
353 R005 = fma(z, term0041, 4 * term0031);
354 R011 = z * term0101;
355 double term0111 = z * term0102;
356 R012 = fma(z, term0111, term0101);
357 double term0112 = z * term0103;
358 double term0121 = fma(z, term0112, term0102);
359 R013 = fma(z, term0121, 2 * term0111);
360 double term0113 = z * term0104;
361 double term0122 = fma(z, term0113, term0103);
362 double term0131 = fma(z, term0122, 2 * term0112);
363 R014 = fma(z, term0131, 3 * term0121);
364 R021 = z * term0201;
365 double term0211 = z * term0202;
366 R022 = fma(z, term0211, term0201);
367 double term0212 = z * term0203;
368 double term0221 = fma(z, term0212, term0202);
369 R023 = fma(z, term0221, 2 * term0211);
370 R031 = z * term0301;
371 double term0311 = z * term0302;
372 R032 = fma(z, term0311, term0301);
373 R041 = z * term0401;
374 R101 = z * term1001;
375 double term1011 = z * term1002;
376 R102 = fma(z, term1011, term1001);
377 double term1012 = z * term1003;
378 double term1021 = fma(z, term1012, term1002);
379 R103 = fma(z, term1021, 2 * term1011);
380 double term1013 = z * term1004;
381 double term1022 = fma(z, term1013, term1003);
382 double term1031 = fma(z, term1022, 2 * term1012);
383 R104 = fma(z, term1031, 3 * term1021);
384 R111 = z * term1101;
385 double term1111 = z * term1102;
386 R112 = fma(z, term1111, term1101);
387 double term1112 = z * term1103;
388 double term1121 = fma(z, term1112, term1102);
389 R113 = fma(z, term1121, 2 * term1111);
390 R121 = z * term1201;
391 double term1211 = z * term1202;
392 R122 = fma(z, term1211, term1201);
393 R131 = z * term1301;
394 R201 = z * term2001;
395 double term2011 = z * term2002;
396 R202 = fma(z, term2011, term2001);
397 double term2012 = z * term2003;
398 double term2021 = fma(z, term2012, term2002);
399 R203 = fma(z, term2021, 2 * term2011);
400 R211 = z * term2101;
401 double term2111 = z * term2102;
402 R212 = fma(z, term2111, term2101);
403 R221 = z * term2201;
404 R301 = z * term3001;
405 double term3011 = z * term3002;
406 R302 = fma(z, term3011, term3001);
407 R311 = z * term3101;
408 R401 = z * term4001;
409 }
410
411
412
413
414 @Override
415 protected void order6() {
416 source(work);
417 double term0000 = work[0];
418 double term0001 = work[1];
419 double term0002 = work[2];
420 double term0003 = work[3];
421 double term0004 = work[4];
422 double term0005 = work[5];
423 double term0006 = work[6];
424 R000 = term0000;
425 R100 = x * term0001;
426 double term1001 = x * term0002;
427 R200 = fma(x, term1001, term0001);
428 double term1002 = x * term0003;
429 double term2001 = fma(x, term1002, term0002);
430 R300 = fma(x, term2001, 2 * term1001);
431 double term1003 = x * term0004;
432 double term2002 = fma(x, term1003, term0003);
433 double term3001 = fma(x, term2002, 2 * term1002);
434 R400 = fma(x, term3001, 3 * term2001);
435 double term1004 = x * term0005;
436 double term2003 = fma(x, term1004, term0004);
437 double term3002 = fma(x, term2003, 2 * term1003);
438 double term4001 = fma(x, term3002, 3 * term2002);
439 R500 = fma(x, term4001, 4 * term3001);
440 double term1005 = x * term0006;
441 double term2004 = fma(x, term1005, term0005);
442 double term3003 = fma(x, term2004, 2 * term1004);
443 double term4002 = fma(x, term3003, 3 * term2003);
444 double term5001 = fma(x, term4002, 4 * term3002);
445 R600 = fma(x, term5001, 5 * term4001);
446 R010 = y * term0001;
447 double term0101 = y * term0002;
448 R020 = fma(y, term0101, term0001);
449 double term0102 = y * term0003;
450 double term0201 = fma(y, term0102, term0002);
451 R030 = fma(y, term0201, 2 * term0101);
452 double term0103 = y * term0004;
453 double term0202 = fma(y, term0103, term0003);
454 double term0301 = fma(y, term0202, 2 * term0102);
455 R040 = fma(y, term0301, 3 * term0201);
456 double term0104 = y * term0005;
457 double term0203 = fma(y, term0104, term0004);
458 double term0302 = fma(y, term0203, 2 * term0103);
459 double term0401 = fma(y, term0302, 3 * term0202);
460 R050 = fma(y, term0401, 4 * term0301);
461 double term0105 = y * term0006;
462 double term0204 = fma(y, term0105, term0005);
463 double term0303 = fma(y, term0204, 2 * term0104);
464 double term0402 = fma(y, term0303, 3 * term0203);
465 double term0501 = fma(y, term0402, 4 * term0302);
466 R060 = fma(y, term0501, 5 * term0401);
467 R110 = y * term1001;
468 double term1101 = y * term1002;
469 R120 = fma(y, term1101, term1001);
470 double term1102 = y * term1003;
471 double term1201 = fma(y, term1102, term1002);
472 R130 = fma(y, term1201, 2 * term1101);
473 double term1103 = y * term1004;
474 double term1202 = fma(y, term1103, term1003);
475 double term1301 = fma(y, term1202, 2 * term1102);
476 R140 = fma(y, term1301, 3 * term1201);
477 double term1104 = y * term1005;
478 double term1203 = fma(y, term1104, term1004);
479 double term1302 = fma(y, term1203, 2 * term1103);
480 double term1401 = fma(y, term1302, 3 * term1202);
481 R150 = fma(y, term1401, 4 * term1301);
482 R210 = y * term2001;
483 double term2101 = y * term2002;
484 R220 = fma(y, term2101, term2001);
485 double term2102 = y * term2003;
486 double term2201 = fma(y, term2102, term2002);
487 R230 = fma(y, term2201, 2 * term2101);
488 double term2103 = y * term2004;
489 double term2202 = fma(y, term2103, term2003);
490 double term2301 = fma(y, term2202, 2 * term2102);
491 R240 = fma(y, term2301, 3 * term2201);
492 R310 = y * term3001;
493 double term3101 = y * term3002;
494 R320 = fma(y, term3101, term3001);
495 double term3102 = y * term3003;
496 double term3201 = fma(y, term3102, term3002);
497 R330 = fma(y, term3201, 2 * term3101);
498 R410 = y * term4001;
499 double term4101 = y * term4002;
500 R420 = fma(y, term4101, term4001);
501 R510 = y * term5001;
502 R001 = z * term0001;
503 double term0011 = z * term0002;
504 R002 = fma(z, term0011, term0001);
505 double term0012 = z * term0003;
506 double term0021 = fma(z, term0012, term0002);
507 R003 = fma(z, term0021, 2 * term0011);
508 double term0013 = z * term0004;
509 double term0022 = fma(z, term0013, term0003);
510 double term0031 = fma(z, term0022, 2 * term0012);
511 R004 = fma(z, term0031, 3 * term0021);
512 double term0014 = z * term0005;
513 double term0023 = fma(z, term0014, term0004);
514 double term0032 = fma(z, term0023, 2 * term0013);
515 double term0041 = fma(z, term0032, 3 * term0022);
516 R005 = fma(z, term0041, 4 * term0031);
517 double term0015 = z * term0006;
518 double term0024 = fma(z, term0015, term0005);
519 double term0033 = fma(z, term0024, 2 * term0014);
520 double term0042 = fma(z, term0033, 3 * term0023);
521 double term0051 = fma(z, term0042, 4 * term0032);
522 R006 = fma(z, term0051, 5 * term0041);
523 R011 = z * term0101;
524 double term0111 = z * term0102;
525 R012 = fma(z, term0111, term0101);
526 double term0112 = z * term0103;
527 double term0121 = fma(z, term0112, term0102);
528 R013 = fma(z, term0121, 2 * term0111);
529 double term0113 = z * term0104;
530 double term0122 = fma(z, term0113, term0103);
531 double term0131 = fma(z, term0122, 2 * term0112);
532 R014 = fma(z, term0131, 3 * term0121);
533 double term0114 = z * term0105;
534 double term0123 = fma(z, term0114, term0104);
535 double term0132 = fma(z, term0123, 2 * term0113);
536 double term0141 = fma(z, term0132, 3 * term0122);
537 R015 = fma(z, term0141, 4 * term0131);
538 R021 = z * term0201;
539 double term0211 = z * term0202;
540 R022 = fma(z, term0211, term0201);
541 double term0212 = z * term0203;
542 double term0221 = fma(z, term0212, term0202);
543 R023 = fma(z, term0221, 2 * term0211);
544 double term0213 = z * term0204;
545 double term0222 = fma(z, term0213, term0203);
546 double term0231 = fma(z, term0222, 2 * term0212);
547 R024 = fma(z, term0231, 3 * term0221);
548 R031 = z * term0301;
549 double term0311 = z * term0302;
550 R032 = fma(z, term0311, term0301);
551 double term0312 = z * term0303;
552 double term0321 = fma(z, term0312, term0302);
553 R033 = fma(z, term0321, 2 * term0311);
554 R041 = z * term0401;
555 double term0411 = z * term0402;
556 R042 = fma(z, term0411, term0401);
557 R051 = z * term0501;
558 R101 = z * term1001;
559 double term1011 = z * term1002;
560 R102 = fma(z, term1011, term1001);
561 double term1012 = z * term1003;
562 double term1021 = fma(z, term1012, term1002);
563 R103 = fma(z, term1021, 2 * term1011);
564 double term1013 = z * term1004;
565 double term1022 = fma(z, term1013, term1003);
566 double term1031 = fma(z, term1022, 2 * term1012);
567 R104 = fma(z, term1031, 3 * term1021);
568 double term1014 = z * term1005;
569 double term1023 = fma(z, term1014, term1004);
570 double term1032 = fma(z, term1023, 2 * term1013);
571 double term1041 = fma(z, term1032, 3 * term1022);
572 R105 = fma(z, term1041, 4 * term1031);
573 R111 = z * term1101;
574 double term1111 = z * term1102;
575 R112 = fma(z, term1111, term1101);
576 double term1112 = z * term1103;
577 double term1121 = fma(z, term1112, term1102);
578 R113 = fma(z, term1121, 2 * term1111);
579 double term1113 = z * term1104;
580 double term1122 = fma(z, term1113, term1103);
581 double term1131 = fma(z, term1122, 2 * term1112);
582 R114 = fma(z, term1131, 3 * term1121);
583 R121 = z * term1201;
584 double term1211 = z * term1202;
585 R122 = fma(z, term1211, term1201);
586 double term1212 = z * term1203;
587 double term1221 = fma(z, term1212, term1202);
588 R123 = fma(z, term1221, 2 * term1211);
589 R131 = z * term1301;
590 double term1311 = z * term1302;
591 R132 = fma(z, term1311, term1301);
592 R141 = z * term1401;
593 R201 = z * term2001;
594 double term2011 = z * term2002;
595 R202 = fma(z, term2011, term2001);
596 double term2012 = z * term2003;
597 double term2021 = fma(z, term2012, term2002);
598 R203 = fma(z, term2021, 2 * term2011);
599 double term2013 = z * term2004;
600 double term2022 = fma(z, term2013, term2003);
601 double term2031 = fma(z, term2022, 2 * term2012);
602 R204 = fma(z, term2031, 3 * term2021);
603 R211 = z * term2101;
604 double term2111 = z * term2102;
605 R212 = fma(z, term2111, term2101);
606 double term2112 = z * term2103;
607 double term2121 = fma(z, term2112, term2102);
608 R213 = fma(z, term2121, 2 * term2111);
609 R221 = z * term2201;
610 double term2211 = z * term2202;
611 R222 = fma(z, term2211, term2201);
612 R231 = z * term2301;
613 R301 = z * term3001;
614 double term3011 = z * term3002;
615 R302 = fma(z, term3011, term3001);
616 double term3012 = z * term3003;
617 double term3021 = fma(z, term3012, term3002);
618 R303 = fma(z, term3021, 2 * term3011);
619 R311 = z * term3101;
620 double term3111 = z * term3102;
621 R312 = fma(z, term3111, term3101);
622 R321 = z * term3201;
623 R401 = z * term4001;
624 double term4011 = z * term4002;
625 R402 = fma(z, term4011, term4001);
626 R411 = z * term4101;
627 R501 = z * term5001;
628 }
629
630
631
632
633 @SuppressWarnings("fallthrough")
634 @Override
635 protected void multipoleIPotentialAtK(PolarizableMultipole mI, int order) {
636 switch (order) {
637 default:
638 case 3:
639
640 double term300 = 0.0;
641 term300 = fma(mI.q, R300, term300);
642 term300 = fma(mI.dx, -R400, term300);
643 term300 = fma(mI.dy, -R310, term300);
644 term300 = fma(mI.dz, -R301, term300);
645 term300 = fma(mI.qxx, R500, term300);
646 term300 = fma(mI.qyy, R320, term300);
647 term300 = fma(mI.qzz, R302, term300);
648 term300 = fma(mI.qxy, R410, term300);
649 term300 = fma(mI.qxz, R401, term300);
650 term300 = fma(mI.qyz, R311, term300);
651 E300 = term300;
652 double term030 = 0.0;
653 term030 = fma(mI.q, R030, term030);
654 term030 = fma(mI.dx, -R130, term030);
655 term030 = fma(mI.dy, -R040, term030);
656 term030 = fma(mI.dz, -R031, term030);
657 term030 = fma(mI.qxx, R230, term030);
658 term030 = fma(mI.qyy, R050, term030);
659 term030 = fma(mI.qzz, R032, term030);
660 term030 = fma(mI.qxy, R140, term030);
661 term030 = fma(mI.qxz, R131, term030);
662 term030 = fma(mI.qyz, R041, term030);
663 E030 = term030;
664 double term003 = 0.0;
665 term003 = fma(mI.q, R003, term003);
666 term003 = fma(mI.dx, -R103, term003);
667 term003 = fma(mI.dy, -R013, term003);
668 term003 = fma(mI.dz, -R004, term003);
669 term003 = fma(mI.qxx, R203, term003);
670 term003 = fma(mI.qyy, R023, term003);
671 term003 = fma(mI.qzz, R005, term003);
672 term003 = fma(mI.qxy, R113, term003);
673 term003 = fma(mI.qxz, R104, term003);
674 term003 = fma(mI.qyz, R014, term003);
675 E003 = term003;
676 double term210 = 0.0;
677 term210 = fma(mI.q, R210, term210);
678 term210 = fma(mI.dx, -R310, term210);
679 term210 = fma(mI.dy, -R220, term210);
680 term210 = fma(mI.dz, -R211, term210);
681 term210 = fma(mI.qxx, R410, term210);
682 term210 = fma(mI.qyy, R230, term210);
683 term210 = fma(mI.qzz, R212, term210);
684 term210 = fma(mI.qxy, R320, term210);
685 term210 = fma(mI.qxz, R311, term210);
686 term210 = fma(mI.qyz, R221, term210);
687 E210 = term210;
688 double term201 = 0.0;
689 term201 = fma(mI.q, R201, term201);
690 term201 = fma(mI.dx, -R301, term201);
691 term201 = fma(mI.dy, -R211, term201);
692 term201 = fma(mI.dz, -R202, term201);
693 term201 = fma(mI.qxx, R401, term201);
694 term201 = fma(mI.qyy, R221, term201);
695 term201 = fma(mI.qzz, R203, term201);
696 term201 = fma(mI.qxy, R311, term201);
697 term201 = fma(mI.qxz, R302, term201);
698 term201 = fma(mI.qyz, R212, term201);
699 E201 = term201;
700 double term120 = 0.0;
701 term120 = fma(mI.q, R120, term120);
702 term120 = fma(mI.dx, -R220, term120);
703 term120 = fma(mI.dy, -R130, term120);
704 term120 = fma(mI.dz, -R121, term120);
705 term120 = fma(mI.qxx, R320, term120);
706 term120 = fma(mI.qyy, R140, term120);
707 term120 = fma(mI.qzz, R122, term120);
708 term120 = fma(mI.qxy, R230, term120);
709 term120 = fma(mI.qxz, R221, term120);
710 term120 = fma(mI.qyz, R131, term120);
711 E120 = term120;
712 double term021 = 0.0;
713 term021 = fma(mI.q, R021, term021);
714 term021 = fma(mI.dx, -R121, term021);
715 term021 = fma(mI.dy, -R031, term021);
716 term021 = fma(mI.dz, -R022, term021);
717 term021 = fma(mI.qxx, R221, term021);
718 term021 = fma(mI.qyy, R041, term021);
719 term021 = fma(mI.qzz, R023, term021);
720 term021 = fma(mI.qxy, R131, term021);
721 term021 = fma(mI.qxz, R122, term021);
722 term021 = fma(mI.qyz, R032, term021);
723 E021 = term021;
724 double term102 = 0.0;
725 term102 = fma(mI.q, R102, term102);
726 term102 = fma(mI.dx, -R202, term102);
727 term102 = fma(mI.dy, -R112, term102);
728 term102 = fma(mI.dz, -R103, term102);
729 term102 = fma(mI.qxx, R302, term102);
730 term102 = fma(mI.qyy, R122, term102);
731 term102 = fma(mI.qzz, R104, term102);
732 term102 = fma(mI.qxy, R212, term102);
733 term102 = fma(mI.qxz, R203, term102);
734 term102 = fma(mI.qyz, R113, term102);
735 E102 = term102;
736 double term012 = 0.0;
737 term012 = fma(mI.q, R012, term012);
738 term012 = fma(mI.dx, -R112, term012);
739 term012 = fma(mI.dy, -R022, term012);
740 term012 = fma(mI.dz, -R013, term012);
741 term012 = fma(mI.qxx, R212, term012);
742 term012 = fma(mI.qyy, R032, term012);
743 term012 = fma(mI.qzz, R014, term012);
744 term012 = fma(mI.qxy, R122, term012);
745 term012 = fma(mI.qxz, R113, term012);
746 term012 = fma(mI.qyz, R023, term012);
747 E012 = term012;
748 double term111 = 0.0;
749 term111 = fma(mI.q, R111, term111);
750 term111 = fma(mI.dx, -R211, term111);
751 term111 = fma(mI.dy, -R121, term111);
752 term111 = fma(mI.dz, -R112, term111);
753 term111 = fma(mI.qxx, R311, term111);
754 term111 = fma(mI.qyy, R131, term111);
755 term111 = fma(mI.qzz, R113, term111);
756 term111 = fma(mI.qxy, R221, term111);
757 term111 = fma(mI.qxz, R212, term111);
758 term111 = fma(mI.qyz, R122, term111);
759 E111 = term111;
760
761 case 2:
762
763 double term200 = 0.0;
764 term200 = fma(mI.q, R200, term200);
765 term200 = fma(mI.dx, -R300, term200);
766 term200 = fma(mI.dy, -R210, term200);
767 term200 = fma(mI.dz, -R201, term200);
768 term200 = fma(mI.qxx, R400, term200);
769 term200 = fma(mI.qyy, R220, term200);
770 term200 = fma(mI.qzz, R202, term200);
771 term200 = fma(mI.qxy, R310, term200);
772 term200 = fma(mI.qxz, R301, term200);
773 term200 = fma(mI.qyz, R211, term200);
774 E200 = term200;
775 double term020 = 0.0;
776 term020 = fma(mI.q, R020, term020);
777 term020 = fma(mI.dx, -R120, term020);
778 term020 = fma(mI.dy, -R030, term020);
779 term020 = fma(mI.dz, -R021, term020);
780 term020 = fma(mI.qxx, R220, term020);
781 term020 = fma(mI.qyy, R040, term020);
782 term020 = fma(mI.qzz, R022, term020);
783 term020 = fma(mI.qxy, R130, term020);
784 term020 = fma(mI.qxz, R121, term020);
785 term020 = fma(mI.qyz, R031, term020);
786 E020 = term020;
787 double term002 = 0.0;
788 term002 = fma(mI.q, R002, term002);
789 term002 = fma(mI.dx, -R102, term002);
790 term002 = fma(mI.dy, -R012, term002);
791 term002 = fma(mI.dz, -R003, term002);
792 term002 = fma(mI.qxx, R202, term002);
793 term002 = fma(mI.qyy, R022, term002);
794 term002 = fma(mI.qzz, R004, term002);
795 term002 = fma(mI.qxy, R112, term002);
796 term002 = fma(mI.qxz, R103, term002);
797 term002 = fma(mI.qyz, R013, term002);
798 E002 = term002;
799 double term110 = 0.0;
800 term110 = fma(mI.q, R110, term110);
801 term110 = fma(mI.dx, -R210, term110);
802 term110 = fma(mI.dy, -R120, term110);
803 term110 = fma(mI.dz, -R111, term110);
804 term110 = fma(mI.qxx, R310, term110);
805 term110 = fma(mI.qyy, R130, term110);
806 term110 = fma(mI.qzz, R112, term110);
807 term110 = fma(mI.qxy, R220, term110);
808 term110 = fma(mI.qxz, R211, term110);
809 term110 = fma(mI.qyz, R121, term110);
810 E110 = term110;
811 double term101 = 0.0;
812 term101 = fma(mI.q, R101, term101);
813 term101 = fma(mI.dx, -R201, term101);
814 term101 = fma(mI.dy, -R111, term101);
815 term101 = fma(mI.dz, -R102, term101);
816 term101 = fma(mI.qxx, R301, term101);
817 term101 = fma(mI.qyy, R121, term101);
818 term101 = fma(mI.qzz, R103, term101);
819 term101 = fma(mI.qxy, R211, term101);
820 term101 = fma(mI.qxz, R202, term101);
821 term101 = fma(mI.qyz, R112, term101);
822 E101 = term101;
823 double term011 = 0.0;
824 term011 = fma(mI.q, R011, term011);
825 term011 = fma(mI.dx, -R111, term011);
826 term011 = fma(mI.dy, -R021, term011);
827 term011 = fma(mI.dz, -R012, term011);
828 term011 = fma(mI.qxx, R211, term011);
829 term011 = fma(mI.qyy, R031, term011);
830 term011 = fma(mI.qzz, R013, term011);
831 term011 = fma(mI.qxy, R121, term011);
832 term011 = fma(mI.qxz, R112, term011);
833 term011 = fma(mI.qyz, R022, term011);
834 E011 = term011;
835
836 case 1:
837
838
839 double term100 = 0.0;
840 term100 = fma(mI.q, R100, term100);
841 term100 = fma(mI.dx, -R200, term100);
842 term100 = fma(mI.dy, -R110, term100);
843 term100 = fma(mI.dz, -R101, term100);
844 term100 = fma(mI.qxx, R300, term100);
845 term100 = fma(mI.qyy, R120, term100);
846 term100 = fma(mI.qzz, R102, term100);
847 term100 = fma(mI.qxy, R210, term100);
848 term100 = fma(mI.qxz, R201, term100);
849 term100 = fma(mI.qyz, R111, term100);
850 E100 = term100;
851
852 double term010 = 0.0;
853 term010 = fma(mI.q, R010, term010);
854 term010 = fma(mI.dx, -R110, term010);
855 term010 = fma(mI.dy, -R020, term010);
856 term010 = fma(mI.dz, -R011, term010);
857 term010 = fma(mI.qxx, R210, term010);
858 term010 = fma(mI.qyy, R030, term010);
859 term010 = fma(mI.qzz, R012, term010);
860 term010 = fma(mI.qxy, R120, term010);
861 term010 = fma(mI.qxz, R111, term010);
862 term010 = fma(mI.qyz, R021, term010);
863 E010 = term010;
864 double term001 = 0.0;
865 term001 = fma(mI.q, R001, term001);
866 term001 = fma(mI.dx, -R101, term001);
867 term001 = fma(mI.dy, -R011, term001);
868 term001 = fma(mI.dz, -R002, term001);
869 term001 = fma(mI.qxx, R201, term001);
870 term001 = fma(mI.qyy, R021, term001);
871 term001 = fma(mI.qzz, R003, term001);
872 term001 = fma(mI.qxy, R111, term001);
873 term001 = fma(mI.qxz, R102, term001);
874 term001 = fma(mI.qyz, R012, term001);
875 E001 = term001;
876
877 case 0:
878
879 double term000 = 0.0;
880 term000 = fma(mI.q, R000, term000);
881 term000 = fma(mI.dx, -R100, term000);
882 term000 = fma(mI.dy, -R010, term000);
883 term000 = fma(mI.dz, -R001, term000);
884 term000 = fma(mI.qxx, R200, term000);
885 term000 = fma(mI.qyy, R020, term000);
886 term000 = fma(mI.qzz, R002, term000);
887 term000 = fma(mI.qxy, R110, term000);
888 term000 = fma(mI.qxz, R101, term000);
889 term000 = fma(mI.qyz, R011, term000);
890 E000 = term000;
891 }
892 }
893
894
895
896
897 @SuppressWarnings("fallthrough")
898 @Override
899 protected void chargeIPotentialAtK(PolarizableMultipole mI, int order) {
900 switch (order) {
901 default:
902 case 3:
903 E300 = mI.q * R300;
904 E030 = mI.q * R030;
905 E003 = mI.q * R003;
906 E210 = mI.q * R210;
907 E201 = mI.q * R201;
908 E120 = mI.q * R120;
909 E021 = mI.q * R021;
910 E102 = mI.q * R102;
911 E012 = mI.q * R012;
912 E111 = mI.q * R111;
913
914 case 2:
915 E200 = mI.q * R200;
916 E020 = mI.q * R020;
917 E002 = mI.q * R002;
918 E110 = mI.q * R110;
919 E101 = mI.q * R101;
920 E011 = mI.q * R011;
921
922 case 1:
923
924
925 E100 = mI.q * R100;
926 E010 = mI.q * R010;
927 E001 = mI.q * R001;
928
929 case 0:
930
931 E000 = mI.q * R000;
932 }
933 }
934
935
936
937
938 @SuppressWarnings("fallthrough")
939 @Override
940 protected void dipoleIPotentialAtK(double uxi, double uyi, double uzi, int order) {
941 switch (order) {
942 case 3:
943
944 double term300 = 0.0;
945 term300 = fma(uxi, -R400, term300);
946 term300 = fma(uyi, -R310, term300);
947 term300 = fma(uzi, -R301, term300);
948 E300 = term300;
949 double term030 = 0.0;
950 term030 = fma(uxi, -R130, term030);
951 term030 = fma(uyi, -R040, term030);
952 term030 = fma(uzi, -R031, term030);
953 E030 = term030;
954 double term003 = 0.0;
955 term003 = fma(uxi, -R103, term003);
956 term003 = fma(uyi, -R013, term003);
957 term003 = fma(uzi, -R004, term003);
958 E003 = term003;
959 double term210 = 0.0;
960 term210 = fma(uxi, -R310, term210);
961 term210 = fma(uyi, -R220, term210);
962 term210 = fma(uzi, -R211, term210);
963 E210 = term210;
964 double term201 = 0.0;
965 term201 = fma(uxi, -R301, term201);
966 term201 = fma(uyi, -R211, term201);
967 term201 = fma(uzi, -R202, term201);
968 E201 = term201;
969 double term120 = 0.0;
970 term120 = fma(uxi, -R220, term120);
971 term120 = fma(uyi, -R130, term120);
972 term120 = fma(uzi, -R121, term120);
973 E120 = term120;
974 double term021 = 0.0;
975 term021 = fma(uxi, -R121, term021);
976 term021 = fma(uyi, -R031, term021);
977 term021 = fma(uzi, -R022, term021);
978 E021 = term021;
979 double term102 = 0.0;
980 term102 = fma(uxi, -R202, term102);
981 term102 = fma(uyi, -R112, term102);
982 term102 = fma(uzi, -R103, term102);
983 E102 = term102;
984 double term012 = 0.0;
985 term012 = fma(uxi, -R112, term012);
986 term012 = fma(uyi, -R022, term012);
987 term012 = fma(uzi, -R013, term012);
988 E012 = term012;
989 double term111 = 0.0;
990 term111 = fma(uxi, -R211, term111);
991 term111 = fma(uyi, -R121, term111);
992 term111 = fma(uzi, -R112, term111);
993 E111 = term111;
994
995 case 2:
996
997 double term200 = -uxi * R300;
998 term200 -= uyi * R210;
999 term200 -= uzi * R201;
1000 E200 = term200;
1001 double term020 = -uxi * R120;
1002 term020 -= uyi * R030;
1003 term020 -= uzi * R021;
1004 E020 = term020;
1005 double term002 = -uxi * R102;
1006 term002 -= uyi * R012;
1007 term002 -= uzi * R003;
1008 E002 = term002;
1009 double term110 = -uxi * R210;
1010 term110 -= uyi * R120;
1011 term110 -= uzi * R111;
1012 E110 = term110;
1013 double term101 = -uxi * R201;
1014 term101 -= uyi * R111;
1015 term101 -= uzi * R102;
1016 E101 = term101;
1017 double term011 = -uxi * R111;
1018 term011 -= uyi * R021;
1019 term011 -= uzi * R012;
1020 E011 = term011;
1021
1022 case 1:
1023
1024 double term100 = -uxi * R200;
1025 term100 -= uyi * R110;
1026 term100 -= uzi * R101;
1027 E100 = term100;
1028 double term010 = -uxi * R110;
1029 term010 -= uyi * R020;
1030 term010 -= uzi * R011;
1031 E010 = term010;
1032 double term001 = -uxi * R101;
1033 term001 -= uyi * R011;
1034 term001 -= uzi * R002;
1035 E001 = term001;
1036
1037 case 0:
1038 double term000 = -uxi * R100;
1039 term000 -= uyi * R010;
1040 term000 -= uzi * R001;
1041 E000 = term000;
1042 }
1043 }
1044
1045
1046
1047
1048 @SuppressWarnings("fallthrough")
1049 @Override
1050 protected void quadrupoleIPotentialAtK(PolarizableMultipole mI, int order) {
1051 switch (order) {
1052 default:
1053 case 3:
1054
1055 double term300 = 0.0;
1056 term300 = fma(mI.qxx, R500, term300);
1057 term300 = fma(mI.qyy, R320, term300);
1058 term300 = fma(mI.qzz, R302, term300);
1059 term300 = fma(mI.qxy, R410, term300);
1060 term300 = fma(mI.qxz, R401, term300);
1061 term300 = fma(mI.qyz, R311, term300);
1062 E300 = term300;
1063 double term030 = 0.0;
1064 term030 = fma(mI.qxx, R230, term030);
1065 term030 = fma(mI.qyy, R050, term030);
1066 term030 = fma(mI.qzz, R032, term030);
1067 term030 = fma(mI.qxy, R140, term030);
1068 term030 = fma(mI.qxz, R131, term030);
1069 term030 = fma(mI.qyz, R041, term030);
1070 E030 = term030;
1071 double term003 = 0.0;
1072 term003 = fma(mI.qxx, R203, term003);
1073 term003 = fma(mI.qyy, R023, term003);
1074 term003 = fma(mI.qzz, R005, term003);
1075 term003 = fma(mI.qxy, R113, term003);
1076 term003 = fma(mI.qxz, R104, term003);
1077 term003 = fma(mI.qyz, R014, term003);
1078 E003 = term003;
1079 double term210 = 0.0;
1080 term210 = fma(mI.qxx, R410, term210);
1081 term210 = fma(mI.qyy, R230, term210);
1082 term210 = fma(mI.qzz, R212, term210);
1083 term210 = fma(mI.qxy, R320, term210);
1084 term210 = fma(mI.qxz, R311, term210);
1085 term210 = fma(mI.qyz, R221, term210);
1086 E210 = term210;
1087 double term201 = 0.0;
1088 term201 = fma(mI.qxx, R401, term201);
1089 term201 = fma(mI.qyy, R221, term201);
1090 term201 = fma(mI.qzz, R203, term201);
1091 term201 = fma(mI.qxy, R311, term201);
1092 term201 = fma(mI.qxz, R302, term201);
1093 term201 = fma(mI.qyz, R212, term201);
1094 E201 = term201;
1095 double term120 = 0.0;
1096 term120 = fma(mI.qxx, R320, term120);
1097 term120 = fma(mI.qyy, R140, term120);
1098 term120 = fma(mI.qzz, R122, term120);
1099 term120 = fma(mI.qxy, R230, term120);
1100 term120 = fma(mI.qxz, R221, term120);
1101 term120 = fma(mI.qyz, R131, term120);
1102 E120 = term120;
1103 double term021 = 0.0;
1104 term021 = fma(mI.qxx, R221, term021);
1105 term021 = fma(mI.qyy, R041, term021);
1106 term021 = fma(mI.qzz, R023, term021);
1107 term021 = fma(mI.qxy, R131, term021);
1108 term021 = fma(mI.qxz, R122, term021);
1109 term021 = fma(mI.qyz, R032, term021);
1110 E021 = term021;
1111 double term102 = 0.0;
1112 term102 = fma(mI.qxx, R302, term102);
1113 term102 = fma(mI.qyy, R122, term102);
1114 term102 = fma(mI.qzz, R104, term102);
1115 term102 = fma(mI.qxy, R212, term102);
1116 term102 = fma(mI.qxz, R203, term102);
1117 term102 = fma(mI.qyz, R113, term102);
1118 E102 = term102;
1119 double term012 = 0.0;
1120 term012 = fma(mI.qxx, R212, term012);
1121 term012 = fma(mI.qyy, R032, term012);
1122 term012 = fma(mI.qzz, R014, term012);
1123 term012 = fma(mI.qxy, R122, term012);
1124 term012 = fma(mI.qxz, R113, term012);
1125 term012 = fma(mI.qyz, R023, term012);
1126 E012 = term012;
1127 double term111 = 0.0;
1128 term111 = fma(mI.qxx, R311, term111);
1129 term111 = fma(mI.qyy, R131, term111);
1130 term111 = fma(mI.qzz, R113, term111);
1131 term111 = fma(mI.qxy, R221, term111);
1132 term111 = fma(mI.qxz, R212, term111);
1133 term111 = fma(mI.qyz, R122, term111);
1134 E111 = term111;
1135
1136 case 2:
1137
1138 double term200 = 0.0;
1139 term200 = fma(mI.qxx, R400, term200);
1140 term200 = fma(mI.qyy, R220, term200);
1141 term200 = fma(mI.qzz, R202, term200);
1142 term200 = fma(mI.qxy, R310, term200);
1143 term200 = fma(mI.qxz, R301, term200);
1144 term200 = fma(mI.qyz, R211, term200);
1145 E200 = term200;
1146 double term020 = 0.0;
1147 term020 = fma(mI.qxx, R220, term020);
1148 term020 = fma(mI.qyy, R040, term020);
1149 term020 = fma(mI.qzz, R022, term020);
1150 term020 = fma(mI.qxy, R130, term020);
1151 term020 = fma(mI.qxz, R121, term020);
1152 term020 = fma(mI.qyz, R031, term020);
1153 E020 = term020;
1154 double term002 = 0.0;
1155 term002 = fma(mI.qxx, R202, term002);
1156 term002 = fma(mI.qyy, R022, term002);
1157 term002 = fma(mI.qzz, R004, term002);
1158 term002 = fma(mI.qxy, R112, term002);
1159 term002 = fma(mI.qxz, R103, term002);
1160 term002 = fma(mI.qyz, R013, term002);
1161 E002 = term002;
1162 double term110 = 0.0;
1163 term110 = fma(mI.qxx, R310, term110);
1164 term110 = fma(mI.qyy, R130, term110);
1165 term110 = fma(mI.qzz, R112, term110);
1166 term110 = fma(mI.qxy, R220, term110);
1167 term110 = fma(mI.qxz, R211, term110);
1168 term110 = fma(mI.qyz, R121, term110);
1169 E110 = term110;
1170 double term101 = 0.0;
1171 term101 = fma(mI.qxx, R301, term101);
1172 term101 = fma(mI.qyy, R121, term101);
1173 term101 = fma(mI.qzz, R103, term101);
1174 term101 = fma(mI.qxy, R211, term101);
1175 term101 = fma(mI.qxz, R202, term101);
1176 term101 = fma(mI.qyz, R112, term101);
1177 E101 = term101;
1178 double term011 = 0.0;
1179 term011 = fma(mI.qxx, R211, term011);
1180 term011 = fma(mI.qyy, R031, term011);
1181 term011 = fma(mI.qzz, R013, term011);
1182 term011 = fma(mI.qxy, R121, term011);
1183 term011 = fma(mI.qxz, R112, term011);
1184 term011 = fma(mI.qyz, R022, term011);
1185 E011 = term011;
1186
1187 case 1:
1188
1189
1190 double term100 = 0.0;
1191 term100 = fma(mI.qxx, R300, term100);
1192 term100 = fma(mI.qyy, R120, term100);
1193 term100 = fma(mI.qzz, R102, term100);
1194 term100 = fma(mI.qxy, R210, term100);
1195 term100 = fma(mI.qxz, R201, term100);
1196 term100 = fma(mI.qyz, R111, term100);
1197 E100 = term100;
1198
1199 double term010 = 0.0;
1200 term010 = fma(mI.qxx, R210, term010);
1201 term010 = fma(mI.qyy, R030, term010);
1202 term010 = fma(mI.qzz, R012, term010);
1203 term010 = fma(mI.qxy, R120, term010);
1204 term010 = fma(mI.qxz, R111, term010);
1205 term010 = fma(mI.qyz, R021, term010);
1206 E010 = term010;
1207 double term001 = 0.0;
1208 term001 = fma(mI.qxx, R201, term001);
1209 term001 = fma(mI.qyy, R021, term001);
1210 term001 = fma(mI.qzz, R003, term001);
1211 term001 = fma(mI.qxy, R111, term001);
1212 term001 = fma(mI.qxz, R102, term001);
1213 term001 = fma(mI.qyz, R012, term001);
1214 E001 = term001;
1215
1216 case 0:
1217
1218 double term000 = 0.0;
1219 term000 = fma(mI.qxx, R200, term000);
1220 term000 = fma(mI.qyy, R020, term000);
1221 term000 = fma(mI.qzz, R002, term000);
1222 term000 = fma(mI.qxy, R110, term000);
1223 term000 = fma(mI.qxz, R101, term000);
1224 term000 = fma(mI.qyz, R011, term000);
1225 E000 = term000;
1226 }
1227 }
1228
1229
1230
1231
1232 @SuppressWarnings("fallthrough")
1233 @Override
1234 protected void multipoleKPotentialAtI(PolarizableMultipole mK, int order) {
1235 switch (order) {
1236 default:
1237 case 3:
1238
1239
1240
1241 double term300 = 0.0;
1242 term300 = fma(mK.q, R300, term300);
1243 term300 = fma(mK.dx, R400, term300);
1244 term300 = fma(mK.dy, R310, term300);
1245 term300 = fma(mK.dz, R301, term300);
1246 term300 = fma(mK.qxx, R500, term300);
1247 term300 = fma(mK.qyy, R320, term300);
1248 term300 = fma(mK.qzz, R302, term300);
1249 term300 = fma(mK.qxy, R410, term300);
1250 term300 = fma(mK.qxz, R401, term300);
1251 term300 = fma(mK.qyz, R311, term300);
1252 E300 = -term300;
1253 double term030 = 0.0;
1254 term030 = fma(mK.q, R030, term030);
1255 term030 = fma(mK.dx, R130, term030);
1256 term030 = fma(mK.dy, R040, term030);
1257 term030 = fma(mK.dz, R031, term030);
1258 term030 = fma(mK.qxx, R230, term030);
1259 term030 = fma(mK.qyy, R050, term030);
1260 term030 = fma(mK.qzz, R032, term030);
1261 term030 = fma(mK.qxy, R140, term030);
1262 term030 = fma(mK.qxz, R131, term030);
1263 term030 = fma(mK.qyz, R041, term030);
1264 E030 = -term030;
1265 double term003 = 0.0;
1266 term003 = fma(mK.q, R003, term003);
1267 term003 = fma(mK.dx, R103, term003);
1268 term003 = fma(mK.dy, R013, term003);
1269 term003 = fma(mK.dz, R004, term003);
1270 term003 = fma(mK.qxx, R203, term003);
1271 term003 = fma(mK.qyy, R023, term003);
1272 term003 = fma(mK.qzz, R005, term003);
1273 term003 = fma(mK.qxy, R113, term003);
1274 term003 = fma(mK.qxz, R104, term003);
1275 term003 = fma(mK.qyz, R014, term003);
1276 E003 = -term003;
1277 double term210 = 0.0;
1278 term210 = fma(mK.q, R210, term210);
1279 term210 = fma(mK.dx, R310, term210);
1280 term210 = fma(mK.dy, R220, term210);
1281 term210 = fma(mK.dz, R211, term210);
1282 term210 = fma(mK.qxx, R410, term210);
1283 term210 = fma(mK.qyy, R230, term210);
1284 term210 = fma(mK.qzz, R212, term210);
1285 term210 = fma(mK.qxy, R320, term210);
1286 term210 = fma(mK.qxz, R311, term210);
1287 term210 = fma(mK.qyz, R221, term210);
1288 E210 = -term210;
1289 double term201 = 0.0;
1290 term201 = fma(mK.q, R201, term201);
1291 term201 = fma(mK.dx, R301, term201);
1292 term201 = fma(mK.dy, R211, term201);
1293 term201 = fma(mK.dz, R202, term201);
1294 term201 = fma(mK.qxx, R401, term201);
1295 term201 = fma(mK.qyy, R221, term201);
1296 term201 = fma(mK.qzz, R203, term201);
1297 term201 = fma(mK.qxy, R311, term201);
1298 term201 = fma(mK.qxz, R302, term201);
1299 term201 = fma(mK.qyz, R212, term201);
1300 E201 = -term201;
1301 double term120 = 0.0;
1302 term120 = fma(mK.q, R120, term120);
1303 term120 = fma(mK.dx, R220, term120);
1304 term120 = fma(mK.dy, R130, term120);
1305 term120 = fma(mK.dz, R121, term120);
1306 term120 = fma(mK.qxx, R320, term120);
1307 term120 = fma(mK.qyy, R140, term120);
1308 term120 = fma(mK.qzz, R122, term120);
1309 term120 = fma(mK.qxy, R230, term120);
1310 term120 = fma(mK.qxz, R221, term120);
1311 term120 = fma(mK.qyz, R131, term120);
1312 E120 = -term120;
1313 double term021 = 0.0;
1314 term021 = fma(mK.q, R021, term021);
1315 term021 = fma(mK.dx, R121, term021);
1316 term021 = fma(mK.dy, R031, term021);
1317 term021 = fma(mK.dz, R022, term021);
1318 term021 = fma(mK.qxx, R221, term021);
1319 term021 = fma(mK.qyy, R041, term021);
1320 term021 = fma(mK.qzz, R023, term021);
1321 term021 = fma(mK.qxy, R131, term021);
1322 term021 = fma(mK.qxz, R122, term021);
1323 term021 = fma(mK.qyz, R032, term021);
1324 E021 = -term021;
1325 double term102 = 0.0;
1326 term102 = fma(mK.q, R102, term102);
1327 term102 = fma(mK.dx, R202, term102);
1328 term102 = fma(mK.dy, R112, term102);
1329 term102 = fma(mK.dz, R103, term102);
1330 term102 = fma(mK.qxx, R302, term102);
1331 term102 = fma(mK.qyy, R122, term102);
1332 term102 = fma(mK.qzz, R104, term102);
1333 term102 = fma(mK.qxy, R212, term102);
1334 term102 = fma(mK.qxz, R203, term102);
1335 term102 = fma(mK.qyz, R113, term102);
1336 E102 = -term102;
1337 double term012 = 0.0;
1338 term012 = fma(mK.q, R012, term012);
1339 term012 = fma(mK.dx, R112, term012);
1340 term012 = fma(mK.dy, R022, term012);
1341 term012 = fma(mK.dz, R013, term012);
1342 term012 = fma(mK.qxx, R212, term012);
1343 term012 = fma(mK.qyy, R032, term012);
1344 term012 = fma(mK.qzz, R014, term012);
1345 term012 = fma(mK.qxy, R122, term012);
1346 term012 = fma(mK.qxz, R113, term012);
1347 term012 = fma(mK.qyz, R023, term012);
1348 E012 = -term012;
1349 double term111 = 0.0;
1350 term111 = fma(mK.q, R111, term111);
1351 term111 = fma(mK.dx, R211, term111);
1352 term111 = fma(mK.dy, R121, term111);
1353 term111 = fma(mK.dz, R112, term111);
1354 term111 = fma(mK.qxx, R311, term111);
1355 term111 = fma(mK.qyy, R131, term111);
1356 term111 = fma(mK.qzz, R113, term111);
1357 term111 = fma(mK.qxy, R221, term111);
1358 term111 = fma(mK.qxz, R212, term111);
1359 term111 = fma(mK.qyz, R122, term111);
1360 E111 = -term111;
1361
1362 case 2:
1363
1364 double term200 = 0.0;
1365 term200 = fma(mK.q, R200, term200);
1366 term200 = fma(mK.dx, R300, term200);
1367 term200 = fma(mK.dy, R210, term200);
1368 term200 = fma(mK.dz, R201, term200);
1369 term200 = fma(mK.qxx, R400, term200);
1370 term200 = fma(mK.qyy, R220, term200);
1371 term200 = fma(mK.qzz, R202, term200);
1372 term200 = fma(mK.qxy, R310, term200);
1373 term200 = fma(mK.qxz, R301, term200);
1374 term200 = fma(mK.qyz, R211, term200);
1375 E200 = term200;
1376 double term020 = 0.0;
1377 term020 = fma(mK.q, R020, term020);
1378 term020 = fma(mK.dx, R120, term020);
1379 term020 = fma(mK.dy, R030, term020);
1380 term020 = fma(mK.dz, R021, term020);
1381 term020 = fma(mK.qxx, R220, term020);
1382 term020 = fma(mK.qyy, R040, term020);
1383 term020 = fma(mK.qzz, R022, term020);
1384 term020 = fma(mK.qxy, R130, term020);
1385 term020 = fma(mK.qxz, R121, term020);
1386 term020 = fma(mK.qyz, R031, term020);
1387 E020 = term020;
1388 double term002 = 0.0;
1389 term002 = fma(mK.q, R002, term002);
1390 term002 = fma(mK.dx, R102, term002);
1391 term002 = fma(mK.dy, R012, term002);
1392 term002 = fma(mK.dz, R003, term002);
1393 term002 = fma(mK.qxx, R202, term002);
1394 term002 = fma(mK.qyy, R022, term002);
1395 term002 = fma(mK.qzz, R004, term002);
1396 term002 = fma(mK.qxy, R112, term002);
1397 term002 = fma(mK.qxz, R103, term002);
1398 term002 = fma(mK.qyz, R013, term002);
1399 E002 = term002;
1400 double term110 = 0.0;
1401 term110 = fma(mK.q, R110, term110);
1402 term110 = fma(mK.dx, R210, term110);
1403 term110 = fma(mK.dy, R120, term110);
1404 term110 = fma(mK.dz, R111, term110);
1405 term110 = fma(mK.qxx, R310, term110);
1406 term110 = fma(mK.qyy, R130, term110);
1407 term110 = fma(mK.qzz, R112, term110);
1408 term110 = fma(mK.qxy, R220, term110);
1409 term110 = fma(mK.qxz, R211, term110);
1410 term110 = fma(mK.qyz, R121, term110);
1411 E110 = term110;
1412 double term101 = 0.0;
1413 term101 = fma(mK.q, R101, term101);
1414 term101 = fma(mK.dx, R201, term101);
1415 term101 = fma(mK.dy, R111, term101);
1416 term101 = fma(mK.dz, R102, term101);
1417 term101 = fma(mK.qxx, R301, term101);
1418 term101 = fma(mK.qyy, R121, term101);
1419 term101 = fma(mK.qzz, R103, term101);
1420 term101 = fma(mK.qxy, R211, term101);
1421 term101 = fma(mK.qxz, R202, term101);
1422 term101 = fma(mK.qyz, R112, term101);
1423 E101 = term101;
1424 double term011 = 0.0;
1425 term011 = fma(mK.q, R011, term011);
1426 term011 = fma(mK.dx, R111, term011);
1427 term011 = fma(mK.dy, R021, term011);
1428 term011 = fma(mK.dz, R012, term011);
1429 term011 = fma(mK.qxx, R211, term011);
1430 term011 = fma(mK.qyy, R031, term011);
1431 term011 = fma(mK.qzz, R013, term011);
1432 term011 = fma(mK.qxy, R121, term011);
1433 term011 = fma(mK.qxz, R112, term011);
1434 term011 = fma(mK.qyz, R022, term011);
1435 E011 = term011;
1436
1437 case 1:
1438
1439
1440 double term100 = 0.0;
1441 term100 = fma(mK.q, R100, term100);
1442 term100 = fma(mK.dx, R200, term100);
1443 term100 = fma(mK.dy, R110, term100);
1444 term100 = fma(mK.dz, R101, term100);
1445 term100 = fma(mK.qxx, R300, term100);
1446 term100 = fma(mK.qyy, R120, term100);
1447 term100 = fma(mK.qzz, R102, term100);
1448 term100 = fma(mK.qxy, R210, term100);
1449 term100 = fma(mK.qxz, R201, term100);
1450 term100 = fma(mK.qyz, R111, term100);
1451 E100 = -term100;
1452 double term010 = 0.0;
1453 term010 = fma(mK.q, R010, term010);
1454 term010 = fma(mK.dx, R110, term010);
1455 term010 = fma(mK.dy, R020, term010);
1456 term010 = fma(mK.dz, R011, term010);
1457 term010 = fma(mK.qxx, R210, term010);
1458 term010 = fma(mK.qyy, R030, term010);
1459 term010 = fma(mK.qzz, R012, term010);
1460 term010 = fma(mK.qxy, R120, term010);
1461 term010 = fma(mK.qxz, R111, term010);
1462 term010 = fma(mK.qyz, R021, term010);
1463 E010 = -term010;
1464 double term001 = 0.0;
1465 term001 = fma(mK.q, R001, term001);
1466 term001 = fma(mK.dx, R101, term001);
1467 term001 = fma(mK.dy, R011, term001);
1468 term001 = fma(mK.dz, R002, term001);
1469 term001 = fma(mK.qxx, R201, term001);
1470 term001 = fma(mK.qyy, R021, term001);
1471 term001 = fma(mK.qzz, R003, term001);
1472 term001 = fma(mK.qxy, R111, term001);
1473 term001 = fma(mK.qxz, R102, term001);
1474 term001 = fma(mK.qyz, R012, term001);
1475 E001 = -term001;
1476
1477 case 0:
1478
1479
1480 double term000 = 0.0;
1481 term000 = fma(mK.q, R000, term000);
1482 term000 = fma(mK.dx, R100, term000);
1483 term000 = fma(mK.dy, R010, term000);
1484 term000 = fma(mK.dz, R001, term000);
1485 term000 = fma(mK.qxx, R200, term000);
1486 term000 = fma(mK.qyy, R020, term000);
1487 term000 = fma(mK.qzz, R002, term000);
1488 term000 = fma(mK.qxy, R110, term000);
1489 term000 = fma(mK.qxz, R101, term000);
1490 term000 = fma(mK.qyz, R011, term000);
1491 E000 = term000;
1492 }
1493 }
1494
1495
1496
1497
1498 @SuppressWarnings("fallthrough")
1499 @Override
1500 protected void chargeKPotentialAtI(PolarizableMultipole mK, int order) {
1501 switch (order) {
1502 default:
1503 case 3:
1504
1505
1506
1507 E300 = -mK.q * R300;
1508 E030 = -mK.q * R030;
1509 E003 = -mK.q * R003;
1510 E210 = -mK.q * R210;
1511 E201 = -mK.q * R201;
1512 E120 = -mK.q * R120;
1513 E021 = -mK.q * R021;
1514 E102 = -mK.q * R102;
1515 E012 = -mK.q * R012;
1516 E111 = -mK.q * R111;
1517
1518 case 2:
1519
1520 E200 = mK.q * R200;
1521 E020 = mK.q * R020;
1522 E002 = mK.q * R002;
1523 E110 = mK.q * R110;
1524 E101 = mK.q * R101;
1525 E011 = mK.q * R011;
1526
1527 case 1:
1528
1529
1530 E100 = -mK.q * R100;
1531 E010 = -mK.q * R010;
1532 E001 = -mK.q * R001;
1533
1534 case 0:
1535
1536
1537 E000 = mK.q * R000;
1538 }
1539 }
1540
1541
1542
1543
1544 @SuppressWarnings("fallthrough")
1545 @Override
1546 protected void dipoleKPotentialAtI(double uxk, double uyk, double uzk, int order) {
1547 switch (order) {
1548 case 3:
1549
1550 double term300 = 0.0;
1551 term300 = fma(uxk, R400, term300);
1552 term300 = fma(uyk, R310, term300);
1553 term300 = fma(uzk, R301, term300);
1554 E300 = -term300;
1555 double term030 = 0.0;
1556 term030 = fma(uxk, R130, term030);
1557 term030 = fma(uyk, R040, term030);
1558 term030 = fma(uzk, R031, term030);
1559 E030 = -term030;
1560 double term003 = 0.0;
1561 term003 = fma(uxk, R103, term003);
1562 term003 = fma(uyk, R013, term003);
1563 term003 = fma(uzk, R004, term003);
1564 E003 = -term003;
1565 double term210 = 0.0;
1566 term210 = fma(uxk, R310, term210);
1567 term210 = fma(uyk, R220, term210);
1568 term210 = fma(uzk, R211, term210);
1569 E210 = -term210;
1570 double term201 = 0.0;
1571 term201 = fma(uxk, R301, term201);
1572 term201 = fma(uyk, R211, term201);
1573 term201 = fma(uzk, R202, term201);
1574 E201 = -term201;
1575 double term120 = 0.0;
1576 term120 = fma(uxk, R220, term120);
1577 term120 = fma(uyk, R130, term120);
1578 term120 = fma(uzk, R121, term120);
1579 E120 = -term120;
1580 double term021 = 0.0;
1581 term021 = fma(uxk, R121, term021);
1582 term021 = fma(uyk, R031, term021);
1583 term021 = fma(uzk, R022, term021);
1584 E021 = -term021;
1585 double term102 = 0.0;
1586 term102 = fma(uxk, R202, term102);
1587 term102 = fma(uyk, R112, term102);
1588 term102 = fma(uzk, R103, term102);
1589 E102 = -term102;
1590 double term012 = 0.0;
1591 term012 = fma(uxk, R112, term012);
1592 term012 = fma(uyk, R022, term012);
1593 term012 = fma(uzk, R013, term012);
1594 E012 = -term012;
1595 double term111 = 0.0;
1596 term111 = fma(uxk, R211, term111);
1597 term111 = fma(uyk, R121, term111);
1598 term111 = fma(uzk, R112, term111);
1599 E111 = -term111;
1600
1601 case 2:
1602
1603 double term200 = uxk * R300;
1604 term200 += uyk * R210;
1605 term200 += uzk * R201;
1606 E200 = term200;
1607 double term020 = uxk * R120;
1608 term020 += uyk * R030;
1609 term020 += uzk * R021;
1610 E020 = term020;
1611 double term002 = uxk * R102;
1612 term002 += uyk * R012;
1613 term002 += uzk * R003;
1614 E002 = term002;
1615 double term110 = uxk * R210;
1616 term110 += uyk * R120;
1617 term110 += uzk * R111;
1618 E110 = term110;
1619 double term101 = uxk * R201;
1620 term101 += uyk * R111;
1621 term101 += uzk * R102;
1622 E101 = term101;
1623 double term011 = uxk * R111;
1624 term011 += uyk * R021;
1625 term011 += uzk * R012;
1626 E011 = term011;
1627
1628 case 1:
1629
1630 double term100 = uxk * R200;
1631 term100 += uyk * R110;
1632 term100 += uzk * R101;
1633 E100 = -term100;
1634 double term010 = uxk * R110;
1635 term010 += uyk * R020;
1636 term010 += uzk * R011;
1637 E010 = -term010;
1638 double term001 = uxk * R101;
1639 term001 += uyk * R011;
1640 term001 += uzk * R002;
1641 E001 = -term001;
1642
1643 case 0:
1644 double term000 = uxk * R100;
1645 term000 += uyk * R010;
1646 term000 += uzk * R001;
1647 E000 = term000;
1648 }
1649 }
1650
1651
1652
1653
1654 @SuppressWarnings("fallthrough")
1655 @Override
1656 protected void quadrupoleKPotentialAtI(PolarizableMultipole mK, int order) {
1657 switch (order) {
1658 default:
1659 case 3:
1660
1661
1662
1663 double term300 = 0.0;
1664 term300 = fma(mK.qxx, R500, term300);
1665 term300 = fma(mK.qyy, R320, term300);
1666 term300 = fma(mK.qzz, R302, term300);
1667 term300 = fma(mK.qxy, R410, term300);
1668 term300 = fma(mK.qxz, R401, term300);
1669 term300 = fma(mK.qyz, R311, term300);
1670 E300 = -term300;
1671 double term030 = 0.0;
1672 term030 = fma(mK.qxx, R230, term030);
1673 term030 = fma(mK.qyy, R050, term030);
1674 term030 = fma(mK.qzz, R032, term030);
1675 term030 = fma(mK.qxy, R140, term030);
1676 term030 = fma(mK.qxz, R131, term030);
1677 term030 = fma(mK.qyz, R041, term030);
1678 E030 = -term030;
1679 double term003 = 0.0;
1680 term003 = fma(mK.qxx, R203, term003);
1681 term003 = fma(mK.qyy, R023, term003);
1682 term003 = fma(mK.qzz, R005, term003);
1683 term003 = fma(mK.qxy, R113, term003);
1684 term003 = fma(mK.qxz, R104, term003);
1685 term003 = fma(mK.qyz, R014, term003);
1686 E003 = -term003;
1687 double term210 = 0.0;
1688 term210 = fma(mK.qxx, R410, term210);
1689 term210 = fma(mK.qyy, R230, term210);
1690 term210 = fma(mK.qzz, R212, term210);
1691 term210 = fma(mK.qxy, R320, term210);
1692 term210 = fma(mK.qxz, R311, term210);
1693 term210 = fma(mK.qyz, R221, term210);
1694 E210 = -term210;
1695 double term201 = 0.0;
1696 term201 = fma(mK.qxx, R401, term201);
1697 term201 = fma(mK.qyy, R221, term201);
1698 term201 = fma(mK.qzz, R203, term201);
1699 term201 = fma(mK.qxy, R311, term201);
1700 term201 = fma(mK.qxz, R302, term201);
1701 term201 = fma(mK.qyz, R212, term201);
1702 E201 = -term201;
1703 double term120 = 0.0;
1704 term120 = fma(mK.qxx, R320, term120);
1705 term120 = fma(mK.qyy, R140, term120);
1706 term120 = fma(mK.qzz, R122, term120);
1707 term120 = fma(mK.qxy, R230, term120);
1708 term120 = fma(mK.qxz, R221, term120);
1709 term120 = fma(mK.qyz, R131, term120);
1710 E120 = -term120;
1711 double term021 = 0.0;
1712 term021 = fma(mK.qxx, R221, term021);
1713 term021 = fma(mK.qyy, R041, term021);
1714 term021 = fma(mK.qzz, R023, term021);
1715 term021 = fma(mK.qxy, R131, term021);
1716 term021 = fma(mK.qxz, R122, term021);
1717 term021 = fma(mK.qyz, R032, term021);
1718 E021 = -term021;
1719 double term102 = 0.0;
1720 term102 = fma(mK.qxx, R302, term102);
1721 term102 = fma(mK.qyy, R122, term102);
1722 term102 = fma(mK.qzz, R104, term102);
1723 term102 = fma(mK.qxy, R212, term102);
1724 term102 = fma(mK.qxz, R203, term102);
1725 term102 = fma(mK.qyz, R113, term102);
1726 E102 = -term102;
1727 double term012 = 0.0;
1728 term012 = fma(mK.qxx, R212, term012);
1729 term012 = fma(mK.qyy, R032, term012);
1730 term012 = fma(mK.qzz, R014, term012);
1731 term012 = fma(mK.qxy, R122, term012);
1732 term012 = fma(mK.qxz, R113, term012);
1733 term012 = fma(mK.qyz, R023, term012);
1734 E012 = -term012;
1735 double term111 = 0.0;
1736 term111 = fma(mK.qxx, R311, term111);
1737 term111 = fma(mK.qyy, R131, term111);
1738 term111 = fma(mK.qzz, R113, term111);
1739 term111 = fma(mK.qxy, R221, term111);
1740 term111 = fma(mK.qxz, R212, term111);
1741 term111 = fma(mK.qyz, R122, term111);
1742 E111 = -term111;
1743
1744 case 2:
1745
1746 double term200 = 0.0;
1747 term200 = fma(mK.qxx, R400, term200);
1748 term200 = fma(mK.qyy, R220, term200);
1749 term200 = fma(mK.qzz, R202, term200);
1750 term200 = fma(mK.qxy, R310, term200);
1751 term200 = fma(mK.qxz, R301, term200);
1752 term200 = fma(mK.qyz, R211, term200);
1753 E200 = term200;
1754 double term020 = 0.0;
1755 term020 = fma(mK.qxx, R220, term020);
1756 term020 = fma(mK.qyy, R040, term020);
1757 term020 = fma(mK.qzz, R022, term020);
1758 term020 = fma(mK.qxy, R130, term020);
1759 term020 = fma(mK.qxz, R121, term020);
1760 term020 = fma(mK.qyz, R031, term020);
1761 E020 = term020;
1762 double term002 = 0.0;
1763 term002 = fma(mK.qxx, R202, term002);
1764 term002 = fma(mK.qyy, R022, term002);
1765 term002 = fma(mK.qzz, R004, term002);
1766 term002 = fma(mK.qxy, R112, term002);
1767 term002 = fma(mK.qxz, R103, term002);
1768 term002 = fma(mK.qyz, R013, term002);
1769 E002 = term002;
1770 double term110 = 0.0;
1771 term110 = fma(mK.qxx, R310, term110);
1772 term110 = fma(mK.qyy, R130, term110);
1773 term110 = fma(mK.qzz, R112, term110);
1774 term110 = fma(mK.qxy, R220, term110);
1775 term110 = fma(mK.qxz, R211, term110);
1776 term110 = fma(mK.qyz, R121, term110);
1777 E110 = term110;
1778 double term101 = 0.0;
1779 term101 = fma(mK.qxx, R301, term101);
1780 term101 = fma(mK.qyy, R121, term101);
1781 term101 = fma(mK.qzz, R103, term101);
1782 term101 = fma(mK.qxy, R211, term101);
1783 term101 = fma(mK.qxz, R202, term101);
1784 term101 = fma(mK.qyz, R112, term101);
1785 E101 = term101;
1786 double term011 = 0.0;
1787 term011 = fma(mK.qxx, R211, term011);
1788 term011 = fma(mK.qyy, R031, term011);
1789 term011 = fma(mK.qzz, R013, term011);
1790 term011 = fma(mK.qxy, R121, term011);
1791 term011 = fma(mK.qxz, R112, term011);
1792 term011 = fma(mK.qyz, R022, term011);
1793 E011 = term011;
1794
1795 case 1:
1796
1797
1798 double term100 = 0.0;
1799 term100 = fma(mK.qxx, R300, term100);
1800 term100 = fma(mK.qyy, R120, term100);
1801 term100 = fma(mK.qzz, R102, term100);
1802 term100 = fma(mK.qxy, R210, term100);
1803 term100 = fma(mK.qxz, R201, term100);
1804 term100 = fma(mK.qyz, R111, term100);
1805 E100 = -term100;
1806 double term010 = 0.0;
1807 term010 = fma(mK.qxx, R210, term010);
1808 term010 = fma(mK.qyy, R030, term010);
1809 term010 = fma(mK.qzz, R012, term010);
1810 term010 = fma(mK.qxy, R120, term010);
1811 term010 = fma(mK.qxz, R111, term010);
1812 term010 = fma(mK.qyz, R021, term010);
1813 E010 = -term010;
1814 double term001 = 0.0;
1815 term001 = fma(mK.qxx, R201, term001);
1816 term001 = fma(mK.qyy, R021, term001);
1817 term001 = fma(mK.qzz, R003, term001);
1818 term001 = fma(mK.qxy, R111, term001);
1819 term001 = fma(mK.qxz, R102, term001);
1820 term001 = fma(mK.qyz, R012, term001);
1821 E001 = -term001;
1822
1823 case 0:
1824
1825
1826 double term000 = 0.0;
1827 term000 = fma(mK.qxx, R200, term000);
1828 term000 = fma(mK.qyy, R020, term000);
1829 term000 = fma(mK.qzz, R002, term000);
1830 term000 = fma(mK.qxy, R110, term000);
1831 term000 = fma(mK.qxz, R101, term000);
1832 term000 = fma(mK.qyz, R011, term000);
1833 E000 = term000;
1834 }
1835 }
1836
1837
1838
1839
1840 @Override
1841 protected double Tlmnj(
1842 final int l, final int m, final int n, final int j, final double[] r, final double[] T000) {
1843 if (m == 0 && n == 0) {
1844 if (l > 1) {
1845 return r[0] * Tlmnj(l - 1, 0, 0, j + 1, r, T000)
1846 + (l - 1) * Tlmnj(l - 2, 0, 0, j + 1, r, T000);
1847 } else if (l == 1) {
1848 return r[0] * Tlmnj(0, 0, 0, j + 1, r, T000);
1849 } else {
1850
1851 return T000[j];
1852 }
1853 } else if (n == 0) {
1854
1855 if (m > 1) {
1856 return r[1] * Tlmnj(l, m - 1, 0, j + 1, r, T000)
1857 + (m - 1) * Tlmnj(l, m - 2, 0, j + 1, r, T000);
1858 }
1859
1860 return r[1] * Tlmnj(l, 0, 0, j + 1, r, T000);
1861 } else {
1862
1863 if (n > 1) {
1864 return r[2] * Tlmnj(l, m, n - 1, j + 1, r, T000)
1865 + (n - 1) * Tlmnj(l, m, n - 2, j + 1, r, T000);
1866 }
1867
1868 return r[2] * Tlmnj(l, m, 0, j + 1, r, T000);
1869 }
1870 }
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880 @Override
1881 protected void noStorageRecursion(double[] r, double[] tensor) {
1882 setR(r);
1883 noStorageRecursion(tensor);
1884 }
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894 @Override
1895 protected void noStorageRecursion(double[] tensor) {
1896 source(T000);
1897 double[] r = {x, y, z};
1898
1899 tensor[0] = T000[0];
1900
1901 for (int l = 1; l <= order; l++) {
1902 tensor[ti(l, 0, 0)] = Tlmnj(l, 0, 0, 0, r, T000);
1903 }
1904
1905 for (int l = 0; l <= o1; l++) {
1906 for (int m = 1; m <= order - l; m++) {
1907 tensor[ti(l, m, 0)] = Tlmnj(l, m, 0, 0, r, T000);
1908 }
1909 }
1910
1911 for (int l = 0; l <= o1; l++) {
1912 for (int m = 0; m <= o1 - l; m++) {
1913 for (int n = 1; n <= order - l - m; n++) {
1914 tensor[ti(l, m, n)] = Tlmnj(l, m, n, 0, r, T000);
1915 }
1916 }
1917 }
1918 }
1919
1920
1921
1922
1923 @Override
1924 protected void recursion(double[] r, double[] tensor) {
1925 setR(r);
1926 recursion(tensor);
1927 }
1928
1929
1930
1931
1932 @Override
1933 protected void recursion(double[] tensor) {
1934 source(work);
1935 tensor[0] = work[0];
1936
1937
1938
1939
1940 double current;
1941 double previous = work[1];
1942
1943 tensor[ti(1, 0, 0)] = x * previous;
1944
1945 for (int l = 2; l < o1; l++) {
1946
1947
1948
1949 current = x * work[l];
1950 int iw = il + l - 1;
1951 work[iw] = current;
1952 for (int a = 1; a < l - 1; a++) {
1953
1954
1955
1956
1957 current = x * current + a * work[iw - il];
1958 iw += il - 1;
1959 work[iw] = current;
1960 }
1961
1962
1963 tensor[ti(l, 0, 0)] = x * current + (l - 1) * previous;
1964 previous = current;
1965 }
1966
1967
1968
1969 for (int l = 0; l < order; l++) {
1970
1971
1972 previous = work[l * il + 1];
1973 tensor[ti(l, 1, 0)] = y * previous;
1974 for (int m = 2; m + l < o1; m++) {
1975
1976 int iw = l * il + m;
1977 current = y * work[iw];
1978 iw += im - 1;
1979 work[iw] = current;
1980 for (int a = 1; a < m - 1; a++) {
1981
1982
1983
1984
1985 current = y * current + a * work[iw - im];
1986 iw += im - 1;
1987 work[iw] = current;
1988 }
1989
1990
1991 tensor[ti(l, m, 0)] = y * current + (m - 1) * previous;
1992 previous = current;
1993 }
1994 }
1995
1996
1997
1998 for (int l = 0; l < order; l++) {
1999 for (int m = 0; m + l < order; m++) {
2000
2001
2002 final int lm = m + l;
2003 final int lilmim = l * il + m * im;
2004 previous = work[lilmim + 1];
2005 tensor[ti(l, m, 1)] = z * previous;
2006 for (int n = 2; lm + n < o1; n++) {
2007
2008 int iw = lilmim + n;
2009 current = z * work[iw];
2010 iw += in - 1;
2011 work[iw] = current;
2012 final int n1 = n - 1;
2013 for (int a = 1; a < n1; a++) {
2014
2015
2016
2017
2018 current = z * current + a * work[iw - in];
2019 iw += in - 1;
2020 work[iw] = current;
2021 }
2022
2023
2024 tensor[ti(l, m, n)] = z * current + n1 * previous;
2025 previous = current;
2026 }
2027 }
2028 }
2029 }
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049 @Override
2050 protected String codeTensorRecursion(final double[] r, final double[] tensor) {
2051 setR(r);
2052 source(work);
2053 StringBuilder sb = new StringBuilder();
2054 tensor[0] = work[0];
2055 if (work[0] > 0) {
2056 sb.append(format("%s = %s;\n", rlmn(0, 0, 0), term(0, 0, 0, 0)));
2057 }
2058
2059
2060
2061
2062 double current;
2063 double previous = work[1];
2064
2065 tensor[ti(1, 0, 0)] = x * previous;
2066 sb.append(format("%s = x * %s;\n", rlmn(1, 0, 0), term(0, 0, 0, 1)));
2067
2068 for (int l = 2; l < o1; l++) {
2069
2070
2071
2072 current = x * work[l];
2073 int iw = il + l - 1;
2074 work[iw] = current;
2075 sb.append(format("double %s = x * %s;\n", term(1, 0, 0, l - 1), term(0, 0, 0, l)));
2076 for (int a = 1; a < l - 1; a++) {
2077
2078
2079
2080
2081 current = x * current + a * work[iw - il];
2082 iw += il - 1;
2083 work[iw] = current;
2084 if (a > 1) {
2085 sb.append(
2086 format(
2087 "double %s = fma(x, %s, %d * %s);\n",
2088 term(a + 1, 0, 0, l - a - 1), term(a, 0, 0, l - a), a, term(a - 1, 0, 0, l - a)));
2089 } else {
2090 sb.append(
2091 format(
2092 "double %s = fma(x, %s, %s);\n",
2093 term(a + 1, 0, 0, l - a - 1), term(a, 0, 0, l - a), term(a - 1, 0, 0, l - a)));
2094 }
2095 }
2096
2097
2098 tensor[ti(l, 0, 0)] = x * current + (l - 1) * previous;
2099 previous = current;
2100 if (l > 2) {
2101 sb.append(
2102 format(
2103 "%s = fma(x, %s, %d * %s);\n",
2104 rlmn(l, 0, 0), term(l - 1, 0, 0, 1), (l - 1), term(l - 2, 0, 0, 1)));
2105 } else {
2106 sb.append(
2107 format(
2108 "%s = fma(x, %s, %s);\n", rlmn(l, 0, 0), term(l - 1, 0, 0, 1),
2109 term(l - 2, 0, 0, 1)));
2110 }
2111 }
2112
2113
2114
2115 for (int l = 0; l < order; l++) {
2116
2117
2118 previous = work[l * il + 1];
2119 tensor[ti(l, 1, 0)] = y * previous;
2120 sb.append(format("%s = y * %s;\n", rlmn(l, 1, 0), term(l, 0, 0, 1)));
2121 for (int m = 2; m + l < o1; m++) {
2122
2123 int iw = l * il + m;
2124 current = y * work[iw];
2125 iw += im - 1;
2126 work[iw] = current;
2127 sb.append(format("double %s = y * %s;\n", term(l, 1, 0, m - 1), term(l, 0, 0, m)));
2128 for (int a = 1; a < m - 1; a++) {
2129
2130
2131
2132
2133 current = y * current + a * work[iw - im];
2134 iw += im - 1;
2135 work[iw] = current;
2136 if (a > 1) {
2137 sb.append(
2138 format(
2139 "double %s = fma(y, %s, %d * %s);\n",
2140 term(l, a + 1, 0, m - a - 1),
2141 term(l, a, 0, m - a),
2142 a,
2143 term(l, a - 1, 0, m - a)));
2144 } else {
2145 sb.append(
2146 format(
2147 "double %s = fma(y, %s, %s);\n",
2148 term(l, a + 1, 0, m - a - 1), term(l, a, 0, m - a), term(l, a - 1, 0, m - a)));
2149 }
2150 }
2151
2152
2153 tensor[ti(l, m, 0)] = y * current + (m - 1) * previous;
2154 previous = current;
2155 if (m > 2) {
2156 sb.append(
2157 format(
2158 "%s = fma(y, %s, %d * %s);\n",
2159 rlmn(l, m, 0), term(l, m - 1, 0, 1), (m - 1), term(l, m - 2, 0, 1)));
2160 } else {
2161 sb.append(
2162 format(
2163 "%s = fma(y, %s, %s);\n",
2164 rlmn(l, m, 0), term(l, m - 1, 0, 1), term(l, m - 2, 0, 1)));
2165 }
2166 }
2167 }
2168
2169
2170
2171 for (int l = 0; l < order; l++) {
2172 for (int m = 0; m + l < order; m++) {
2173
2174
2175 final int lm = m + l;
2176 final int lilmim = l * il + m * im;
2177 previous = work[lilmim + 1];
2178 tensor[ti(l, m, 1)] = z * previous;
2179 sb.append(format("%s = z * %s;\n", rlmn(l, m, 1), term(l, m, 0, 1)));
2180 for (int n = 2; lm + n < o1; n++) {
2181
2182 int iw = lilmim + n;
2183 current = z * work[iw];
2184 iw += in - 1;
2185 work[iw] = current;
2186 sb.append(format("double %s = z * %s;\n", term(l, m, 1, n - 1), term(l, m, 0, n)));
2187 final int n1 = n - 1;
2188 for (int a = 1; a < n1; a++) {
2189
2190
2191
2192
2193 current = z * current + a * work[iw - in];
2194 iw += in - 1;
2195 work[iw] = current;
2196 if (a > 1) {
2197 sb.append(
2198 format(
2199 "double %s = fma(z, %s, %d * %s);\n",
2200 term(l, m, a + 1, n - a - 1),
2201 term(l, m, a, n - a),
2202 a,
2203 term(l, m, a - 1, n - a)));
2204 } else {
2205 sb.append(
2206 format(
2207 "double %s = fma(z, %s, %s);\n",
2208 term(l, m, a + 1, n - a - 1),
2209 term(l, m, a, n - a),
2210 term(l, m, a - 1, n - a)));
2211 }
2212 }
2213
2214
2215 tensor[ti(l, m, n)] = z * current + n1 * previous;
2216 previous = current;
2217 if (n > 2) {
2218 sb.append(
2219 format(
2220 "%s = fma(z, %s, %d * %s);\n",
2221 rlmn(l, m, n), term(l, m, n - 1, 1), (n - 1), term(l, m, n - 2, 1)));
2222 } else {
2223 sb.append(
2224 format(
2225 "%s = fma(z, %s, %s);\n",
2226 rlmn(l, m, n), term(l, m, n - 1, 1), term(l, m, n - 2, 1)));
2227 }
2228 }
2229 }
2230 }
2231 return sb.toString();
2232 }
2233 }