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