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.rlmn;
41 import static ffx.numerics.multipole.MultipoleUtilities.term;
42 import static java.lang.Math.fma;
43 import static java.lang.String.format;
44 import static org.apache.commons.math3.util.FastMath.sqrt;
45
46
47
48
49
50
51
52
53
54
55
56
57 public class CoulombTensorQI extends MultipoleTensor {
58
59
60
61
62
63
64 public CoulombTensorQI(int order) {
65 super(CoordinateSystem.QI, order);
66 operator = Operator.COULOMB;
67 }
68
69
70
71
72 @Override
73 public void setR(double dx, double dy, double dz) {
74 x = 0.0;
75 y = 0.0;
76 r2 = dx * dx + dy * dy + dz * dz;
77 z = sqrt(r2);
78 R = z;
79 }
80
81
82
83
84
85
86 protected void source(double[] T000) {
87
88
89 double ir = 1.0 / R;
90 double ir2 = ir * ir;
91 for (int n = 0; n < o1; n++) {
92 T000[n] = coulombSource[n] * ir;
93 ir *= ir2;
94 }
95 }
96
97
98
99
100 @Override
101 protected void order1() {
102 source(work);
103 double term0000 = work[0];
104 double term0001 = work[1];
105 R000 = term0000;
106 R001 = z * term0001;
107 }
108
109
110
111
112 @Override
113 protected void order2() {
114 source(work);
115 double term0000 = work[0];
116 double term0001 = work[1];
117 double term0002 = work[2];
118 R000 = term0000;
119 R200 = term0001;
120 R020 = term0001;
121 R001 = z * term0001;
122 double term0011 = z * term0002;
123 R002 = fma(z, term0011, term0001);
124 }
125
126
127
128
129 @Override
130 protected void order3() {
131 source(work);
132 double term0000 = work[0];
133 double term0001 = work[1];
134 double term0002 = work[2];
135 double term0003 = work[3];
136 R000 = term0000;
137 R200 = term0001;
138 double term2001 = term0002;
139 R020 = term0001;
140 double term0201 = term0002;
141 R001 = z * term0001;
142 double term0011 = z * term0002;
143 R002 = fma(z, term0011, term0001);
144 double term0012 = z * term0003;
145 double term0021 = fma(z, term0012, term0002);
146 R003 = fma(z, term0021, 2 * term0011);
147 R021 = z * term0201;
148 R201 = z * term2001;
149 }
150
151
152
153
154 @Override
155 protected void order4() {
156 source(work);
157 double term0000 = work[0];
158 double term0001 = work[1];
159 double term0002 = work[2];
160 double term0003 = work[3];
161 double term0004 = work[4];
162 R000 = term0000;
163 R200 = term0001;
164 double term2001 = term0002;
165 double term2002 = term0003;
166 R400 = 3 * term2001;
167 R020 = term0001;
168 double term0201 = term0002;
169 double term0202 = term0003;
170 R040 = 3 * term0201;
171 R220 = term2001;
172 R001 = z * term0001;
173 double term0011 = z * term0002;
174 R002 = fma(z, term0011, term0001);
175 double term0012 = z * term0003;
176 double term0021 = fma(z, term0012, term0002);
177 R003 = fma(z, term0021, 2 * term0011);
178 double term0013 = z * term0004;
179 double term0022 = fma(z, term0013, term0003);
180 double term0031 = fma(z, term0022, 2 * term0012);
181 R004 = fma(z, term0031, 3 * term0021);
182 R021 = z * term0201;
183 double term0211 = z * term0202;
184 R022 = fma(z, term0211, term0201);
185 R201 = z * term2001;
186 double term2011 = z * term2002;
187 R202 = fma(z, term2011, term2001);
188 }
189
190
191
192
193 @Override
194 protected void order5() {
195 source(work);
196 double term0000 = work[0];
197 double term0001 = work[1];
198 double term0002 = work[2];
199 double term0003 = work[3];
200 double term0004 = work[4];
201 double term0005 = work[5];
202 R000 = term0000;
203 R200 = term0001;
204 double term2001 = term0002;
205 double term2002 = term0003;
206 R400 = 3 * term2001;
207 double term2003 = term0004;
208 double term4001 = 3 * term2002;
209 R020 = term0001;
210 double term0201 = term0002;
211 double term0202 = term0003;
212 R040 = 3 * term0201;
213 double term0203 = term0004;
214 double term0401 = 3 * term0202;
215 R220 = term2001;
216 double term2201 = term2002;
217 R001 = z * term0001;
218 double term0011 = z * term0002;
219 R002 = fma(z, term0011, term0001);
220 double term0012 = z * term0003;
221 double term0021 = fma(z, term0012, term0002);
222 R003 = fma(z, term0021, 2 * term0011);
223 double term0013 = z * term0004;
224 double term0022 = fma(z, term0013, term0003);
225 double term0031 = fma(z, term0022, 2 * term0012);
226 R004 = fma(z, term0031, 3 * term0021);
227 double term0014 = z * term0005;
228 double term0023 = fma(z, term0014, term0004);
229 double term0032 = fma(z, term0023, 2 * term0013);
230 double term0041 = fma(z, term0032, 3 * term0022);
231 R005 = fma(z, term0041, 4 * term0031);
232 R021 = z * term0201;
233 double term0211 = z * term0202;
234 R022 = fma(z, term0211, term0201);
235 double term0212 = z * term0203;
236 double term0221 = fma(z, term0212, term0202);
237 R023 = fma(z, term0221, 2 * term0211);
238 R041 = z * term0401;
239 R201 = z * term2001;
240 double term2011 = z * term2002;
241 R202 = fma(z, term2011, term2001);
242 double term2012 = z * term2003;
243 double term2021 = fma(z, term2012, term2002);
244 R203 = fma(z, term2021, 2 * term2011);
245 R221 = z * term2201;
246 R401 = z * term4001;
247 }
248
249
250
251
252 @Override
253 protected void order6() {
254 source(work);
255 double term0000 = work[0];
256 double term0001 = work[1];
257 double term0002 = work[2];
258 double term0003 = work[3];
259 double term0004 = work[4];
260 double term0005 = work[5];
261 double term0006 = work[6];
262 R000 = term0000;
263 R200 = term0001;
264 double term2001 = term0002;
265 double term2002 = term0003;
266 R400 = 3 * term2001;
267 double term2003 = term0004;
268 double term4001 = 3 * term2002;
269 double term2004 = term0005;
270 double term4002 = 3 * term2003;
271 R600 = 5 * term4001;
272 R020 = term0001;
273 double term0201 = term0002;
274 double term0202 = term0003;
275 R040 = 3 * term0201;
276 double term0203 = term0004;
277 double term0401 = 3 * term0202;
278 double term0204 = term0005;
279 double term0402 = 3 * term0203;
280 R060 = 5 * term0401;
281 R220 = term2001;
282 double term2201 = term2002;
283 double term2202 = term2003;
284 R240 = 3 * term2201;
285 R420 = term4001;
286 R001 = z * term0001;
287 double term0011 = z * term0002;
288 R002 = fma(z, term0011, term0001);
289 double term0012 = z * term0003;
290 double term0021 = fma(z, term0012, term0002);
291 R003 = fma(z, term0021, 2 * term0011);
292 double term0013 = z * term0004;
293 double term0022 = fma(z, term0013, term0003);
294 double term0031 = fma(z, term0022, 2 * term0012);
295 R004 = fma(z, term0031, 3 * term0021);
296 double term0014 = z * term0005;
297 double term0023 = fma(z, term0014, term0004);
298 double term0032 = fma(z, term0023, 2 * term0013);
299 double term0041 = fma(z, term0032, 3 * term0022);
300 R005 = fma(z, term0041, 4 * term0031);
301 double term0015 = z * term0006;
302 double term0024 = fma(z, term0015, term0005);
303 double term0033 = fma(z, term0024, 2 * term0014);
304 double term0042 = fma(z, term0033, 3 * term0023);
305 double term0051 = fma(z, term0042, 4 * term0032);
306 R006 = fma(z, term0051, 5 * term0041);
307 R021 = z * term0201;
308 double term0211 = z * term0202;
309 R022 = fma(z, term0211, term0201);
310 double term0212 = z * term0203;
311 double term0221 = fma(z, term0212, term0202);
312 R023 = fma(z, term0221, 2 * term0211);
313 double term0213 = z * term0204;
314 double term0222 = fma(z, term0213, term0203);
315 double term0231 = fma(z, term0222, 2 * term0212);
316 R024 = fma(z, term0231, 3 * term0221);
317 R041 = z * term0401;
318 double term0411 = z * term0402;
319 R042 = fma(z, term0411, term0401);
320 R201 = z * term2001;
321 double term2011 = z * term2002;
322 R202 = fma(z, term2011, term2001);
323 double term2012 = z * term2003;
324 double term2021 = fma(z, term2012, term2002);
325 R203 = fma(z, term2021, 2 * term2011);
326 double term2013 = z * term2004;
327 double term2022 = fma(z, term2013, term2003);
328 double term2031 = fma(z, term2022, 2 * term2012);
329 R204 = fma(z, term2031, 3 * term2021);
330 R221 = z * term2201;
331 double term2211 = z * term2202;
332 R222 = fma(z, term2211, term2201);
333 R401 = z * term4001;
334 double term4011 = z * term4002;
335 R402 = fma(z, term4011, term4001);
336 }
337
338
339
340
341 @SuppressWarnings("fallthrough")
342 @Override
343 protected void multipoleIPotentialAtK(PolarizableMultipole mI, int order) {
344 switch (order) {
345 default:
346 case 3:
347
348 double term300 = 0.0;
349 term300 = fma(mI.dx, -R400, term300);
350 term300 = fma(mI.qxz, R401, term300);
351 E300 = term300;
352 double term030 = 0.0;
353 term030 = fma(mI.dy, -R040, term030);
354 term030 = fma(mI.qyz, R041, term030);
355 E030 = term030;
356 double term003 = 0.0;
357 term003 = fma(mI.q, R003, term003);
358 term003 = fma(mI.dz, -R004, term003);
359 term003 = fma(mI.qxx, R203, term003);
360 term003 = fma(mI.qyy, R023, term003);
361 term003 = fma(mI.qzz, R005, term003);
362 E003 = term003;
363 double term210 = 0.0;
364 term210 = fma(mI.dy, -R220, term210);
365 term210 = fma(mI.qyz, R221, term210);
366 E210 = term210;
367 double term201 = 0.0;
368 term201 = fma(mI.q, R201, term201);
369 term201 = fma(mI.dz, -R202, term201);
370 term201 = fma(mI.qxx, R401, term201);
371 term201 = fma(mI.qyy, R221, term201);
372 term201 = fma(mI.qzz, R203, term201);
373 E201 = term201;
374 double term120 = 0.0;
375 term120 = fma(mI.dx, -R220, term120);
376 term120 = fma(mI.qxz, R221, term120);
377 E120 = term120;
378 double term021 = 0.0;
379 term021 = fma(mI.q, R021, term021);
380 term021 = fma(mI.dz, -R022, term021);
381 term021 = fma(mI.qxx, R221, term021);
382 term021 = fma(mI.qyy, R041, term021);
383 term021 = fma(mI.qzz, R023, term021);
384 E021 = term021;
385 double term102 = 0.0;
386 term102 = fma(mI.dx, -R202, term102);
387 term102 = fma(mI.qxz, R203, term102);
388 E102 = term102;
389 double term012 = 0.0;
390 term012 = fma(mI.dy, -R022, term012);
391 term012 = fma(mI.qyz, R023, term012);
392 E012 = term012;
393 double term111 = 0.0;
394 term111 = fma(mI.qxy, R221, term111);
395 E111 = term111;
396
397 case 2:
398
399 double term200 = 0.0;
400 term200 = fma(mI.q, R200, term200);
401 term200 = fma(mI.dz, -R201, term200);
402 term200 = fma(mI.qxx, R400, term200);
403 term200 = fma(mI.qyy, R220, term200);
404 term200 = fma(mI.qzz, R202, term200);
405 E200 = term200;
406 double term020 = 0.0;
407 term020 = fma(mI.q, R020, term020);
408 term020 = fma(mI.dz, -R021, term020);
409 term020 = fma(mI.qxx, R220, term020);
410 term020 = fma(mI.qyy, R040, term020);
411 term020 = fma(mI.qzz, R022, term020);
412 E020 = term020;
413 double term002 = 0.0;
414 term002 = fma(mI.q, R002, term002);
415 term002 = fma(mI.dz, -R003, term002);
416 term002 = fma(mI.qxx, R202, term002);
417 term002 = fma(mI.qyy, R022, term002);
418 term002 = fma(mI.qzz, R004, term002);
419 E002 = term002;
420 double term110 = 0.0;
421 term110 = fma(mI.qxy, R220, term110);
422 E110 = term110;
423 double term101 = 0.0;
424 term101 = fma(mI.dx, -R201, term101);
425 term101 = fma(mI.qxz, R202, term101);
426 E101 = term101;
427 double term011 = 0.0;
428 term011 = fma(mI.dy, -R021, term011);
429 term011 = fma(mI.qyz, R022, term011);
430 E011 = term011;
431
432 case 1:
433
434 double term100 = 0.0;
435 term100 = fma(mI.dx, -R200, term100);
436 term100 = fma(mI.qxz, R201, term100);
437 E100 = term100;
438 double term010 = 0.0;
439 term010 = fma(mI.dy, -R020, term010);
440 term010 = fma(mI.qyz, R021, term010);
441 E010 = term010;
442 double term001 = 0.0;
443 term001 = fma(mI.q, R001, term001);
444 term001 = fma(mI.dz, -R002, term001);
445 term001 = fma(mI.qxx, R201, term001);
446 term001 = fma(mI.qyy, R021, term001);
447 term001 = fma(mI.qzz, R003, term001);
448 E001 = term001;
449
450 case 0:
451 double term000 = 0.0;
452 term000 = fma(mI.q, R000, term000);
453 term000 = fma(mI.dz, -R001, term000);
454 term000 = fma(mI.qxx, R200, term000);
455 term000 = fma(mI.qyy, R020, term000);
456 term000 = fma(mI.qzz, R002, term000);
457 E000 = term000;
458 }
459 }
460
461
462
463
464 @SuppressWarnings("fallthrough")
465 @Override
466 protected void chargeIPotentialAtK(PolarizableMultipole mI, int order) {
467 switch (order) {
468 default:
469 case 3:
470
471 E300 = 0.0;
472 E030 = 0.0;
473 E003 = mI.q * R003;
474 E210 = 0.0;
475 E201 = mI.q * R201;
476 E120 = 0.0;
477 E021 = mI.q * R021;
478 E102 = 0.0;
479 E012 = 0.0;
480 E111 = 0.0;
481
482 case 2:
483
484 E200 = mI.q * R200;
485 E020 = mI.q * R020;
486 E002 = mI.q * R002;
487 E110 = 0.0;
488 E101 = 0.0;
489 E011 = 0.0;
490
491 case 1:
492
493 E100 = 0.0;
494 E010 = 0.0;
495 E001 = mI.q * R001;
496
497 case 0:
498 E000 = mI.q * R000;
499 }
500 }
501
502
503
504
505 @SuppressWarnings("fallthrough")
506 @Override
507 protected void dipoleIPotentialAtK(double uxi, double uyi, double uzi, int order) {
508 switch (order) {
509 case 3:
510 default:
511
512 E300 = -uxi * R400;
513 E030 = -uyi * R040;
514 E003 = -uzi * R004;
515 E210 = -uyi * R220;
516 E201 = -uzi * R202;
517 E120 = -uxi * R220;
518 E021 = -uzi * R022;
519 E102 = -uxi * R202;
520 E012 = -uyi * R022;
521 E111 = 0.0;
522
523 case 2:
524
525 E200 = -uzi * R201;
526 E020 = -uzi * R021;
527 E002 = -uzi * R003;
528 E110 = 0.0;
529 E101 = -uxi * R201;
530 E011 = -uyi * R021;
531
532 case 1:
533
534 E100 = -uxi * R200;
535 E010 = -uyi * R020;
536 E001 = -uzi * R002;
537
538 case 0:
539 E000 = -uzi * R001;
540 }
541 }
542
543
544
545
546 @SuppressWarnings("fallthrough")
547 @Override
548 protected void quadrupoleIPotentialAtK(PolarizableMultipole mI, int order) {
549 switch (order) {
550 default:
551 case 3:
552
553 E300 = mI.qxz * R401;
554 E030 = mI.qyz * R041;
555 double term003 = 0.0;
556 term003 = fma(mI.qxx, R203, term003);
557 term003 = fma(mI.qyy, R023, term003);
558 term003 = fma(mI.qzz, R005, term003);
559 E003 = term003;
560 E210 = mI.qyz * R221;
561 double term201 = 0.0;
562 term201 = fma(mI.qxx, R401, term201);
563 term201 = fma(mI.qyy, R221, term201);
564 term201 = fma(mI.qzz, R203, term201);
565 E201 = term201;
566 E120 = mI.qxz * R221;
567 double term021 = 0.0;
568 term021 = fma(mI.qxx, R221, term021);
569 term021 = fma(mI.qyy, R041, term021);
570 term021 = fma(mI.qzz, R023, term021);
571 E021 = term021;
572 E102 = mI.qxz * R203;
573 E012 = mI.qyz * R023;
574 E111 = mI.qxy * R221;
575
576 case 2:
577
578 double term200 = 0.0;
579 term200 = fma(mI.qxx, R400, term200);
580 term200 = fma(mI.qyy, R220, term200);
581 term200 = fma(mI.qzz, R202, term200);
582 E200 = term200;
583 double term020 = 0.0;
584 term020 = fma(mI.qxx, R220, term020);
585 term020 = fma(mI.qyy, R040, term020);
586 term020 = fma(mI.qzz, R022, term020);
587 E020 = term020;
588 double term002 = 0.0;
589 term002 = fma(mI.qxx, R202, term002);
590 term002 = fma(mI.qyy, R022, term002);
591 term002 = fma(mI.qzz, R004, term002);
592 E002 = term002;
593 E110 = mI.qxy * R220;
594 E101 = mI.qxz * R202;
595 E011 = mI.qyz * R022;
596
597 case 1:
598
599 E100 = mI.qxz * R201;
600 E010 = mI.qyz * R021;
601 double term001 = 0.0;
602 term001 = fma(mI.qxx, R201, term001);
603 term001 = fma(mI.qyy, R021, term001);
604 term001 = fma(mI.qzz, R003, term001);
605 E001 = term001;
606
607 case 0:
608 double term000 = 0.0;
609 term000 = fma(mI.qxx, R200, term000);
610 term000 = fma(mI.qyy, R020, term000);
611 term000 = fma(mI.qzz, R002, term000);
612 E000 = term000;
613 }
614 }
615
616
617
618
619 @SuppressWarnings("fallthrough")
620 @Override
621 protected void multipoleKPotentialAtI(PolarizableMultipole mK, int order) {
622 switch (order) {
623 default:
624 case 3:
625
626 double term300 = 0.0;
627 term300 = fma(mK.dx, R400, term300);
628 term300 = fma(mK.qxz, R401, term300);
629 E300 = -term300;
630 double term030 = 0.0;
631 term030 = fma(mK.dy, R040, term030);
632 term030 = fma(mK.qyz, R041, term030);
633 E030 = -term030;
634 double term003 = 0.0;
635 term003 = fma(mK.q, R003, term003);
636 term003 = fma(mK.dz, R004, term003);
637 term003 = fma(mK.qxx, R203, term003);
638 term003 = fma(mK.qyy, R023, term003);
639 term003 = fma(mK.qzz, R005, term003);
640 E003 = -term003;
641 double term210 = 0.0;
642 term210 = fma(mK.dy, R220, term210);
643 term210 = fma(mK.qyz, R221, term210);
644 E210 = -term210;
645 double term201 = 0.0;
646 term201 = fma(mK.q, R201, term201);
647 term201 = fma(mK.dz, R202, term201);
648 term201 = fma(mK.qxx, R401, term201);
649 term201 = fma(mK.qyy, R221, term201);
650 term201 = fma(mK.qzz, R203, term201);
651 E201 = -term201;
652 double term120 = 0.0;
653 term120 = fma(mK.dx, R220, term120);
654 term120 = fma(mK.qxz, R221, term120);
655 E120 = -term120;
656 double term021 = 0.0;
657 term021 = fma(mK.q, R021, term021);
658 term021 = fma(mK.dz, R022, term021);
659 term021 = fma(mK.qxx, R221, term021);
660 term021 = fma(mK.qyy, R041, term021);
661 term021 = fma(mK.qzz, R023, term021);
662 E021 = -term021;
663 double term102 = 0.0;
664 term102 = fma(mK.dx, R202, term102);
665 term102 = fma(mK.qxz, R203, term102);
666 E102 = -term102;
667 double term012 = 0.0;
668 term012 = fma(mK.dy, R022, term012);
669 term012 = fma(mK.qyz, R023, term012);
670 E012 = -term012;
671 double term111 = 0.0;
672 term111 = fma(mK.qxy, R221, term111);
673 E111 = -term111;
674
675 case 2:
676
677 double term200 = 0.0;
678 term200 = fma(mK.q, R200, term200);
679 term200 = fma(mK.dz, R201, term200);
680 term200 = fma(mK.qxx, R400, term200);
681 term200 = fma(mK.qyy, R220, term200);
682 term200 = fma(mK.qzz, R202, term200);
683 E200 = term200;
684 double term020 = 0.0;
685 term020 = fma(mK.q, R020, term020);
686 term020 = fma(mK.dz, R021, term020);
687 term020 = fma(mK.qxx, R220, term020);
688 term020 = fma(mK.qyy, R040, term020);
689 term020 = fma(mK.qzz, R022, term020);
690 E020 = term020;
691 double term002 = 0.0;
692 term002 = fma(mK.q, R002, term002);
693 term002 = fma(mK.dz, R003, term002);
694 term002 = fma(mK.qxx, R202, term002);
695 term002 = fma(mK.qyy, R022, term002);
696 term002 = fma(mK.qzz, R004, term002);
697 E002 = term002;
698 double term110 = 0.0;
699 term110 = fma(mK.qxy, R220, term110);
700 E110 = term110;
701 double term101 = 0.0;
702 term101 = fma(mK.dx, R201, term101);
703 term101 = fma(mK.qxz, R202, term101);
704 E101 = term101;
705 double term011 = 0.0;
706 term011 = fma(mK.dy, R021, term011);
707 term011 = fma(mK.qyz, R022, term011);
708 E011 = term011;
709
710 case 1:
711
712 double term100 = 0.0;
713 term100 = fma(mK.dx, R200, term100);
714 term100 = fma(mK.qxz, R201, term100);
715 E100 = -term100;
716 double term010 = 0.0;
717 term010 = fma(mK.dy, R020, term010);
718 term010 = fma(mK.qyz, R021, term010);
719 E010 = -term010;
720 double term001 = 0.0;
721 term001 = fma(mK.q, R001, term001);
722 term001 = fma(mK.dz, R002, term001);
723 term001 = fma(mK.qxx, R201, term001);
724 term001 = fma(mK.qyy, R021, term001);
725 term001 = fma(mK.qzz, R003, term001);
726 E001 = -term001;
727 case 0:
728 double term000 = 0.0;
729 term000 = fma(mK.q, R000, term000);
730 term000 = fma(mK.dz, R001, term000);
731 term000 = fma(mK.qxx, R200, term000);
732 term000 = fma(mK.qyy, R020, term000);
733 term000 = fma(mK.qzz, R002, term000);
734 E000 = term000;
735 }
736 }
737
738
739
740
741 @SuppressWarnings("fallthrough")
742 @Override
743 protected void chargeKPotentialAtI(PolarizableMultipole mK, int order) {
744 switch (order) {
745 default:
746 case 3:
747
748 E300 = 0.0;
749 E030 = 0.0;
750 E003 = -mK.q * R003;
751 E210 = 0.0;
752 E201 = -mK.q * R201;
753 E120 = 0.0;
754 E021 = -mK.q * R021;
755 E102 = 0.0;
756 E012 = 0.0;
757 E111 = 0.0;
758
759 case 2:
760
761 E200 = mK.q * R200;
762 E020 = mK.q * R020;
763 E002 = mK.q * R002;
764 E110 = 0.0;
765 E101 = 0.0;
766 E011 = 0.0;
767
768 case 1:
769
770 E100 = 0.0;
771 E010 = 0.0;
772 E001 = -mK.q * R001;
773 case 0:
774 E000 = mK.q * R000;
775 }
776 }
777
778
779
780
781 @SuppressWarnings("fallthrough")
782 @Override
783 protected void dipoleKPotentialAtI(double uxk, double uyk, double uzk, int order) {
784 switch (order) {
785 case 3:
786 default:
787
788 E300 = -uxk * R400;
789 E030 = -uyk * R040;
790 E003 = -uzk * R004;
791 E210 = -uyk * R220;
792 E201 = -uzk * R202;
793 E120 = -uxk * R220;
794 E021 = -uzk * R022;
795 E102 = -uxk * R202;
796 E012 = -uyk * R022;
797 E111 = 0.0;
798
799 case 2:
800 E200 = uzk * R201;
801 E020 = uzk * R021;
802 E002 = uzk * R003;
803 E110 = 0.0;
804 E101 = uxk * R201;
805 E011 = uyk * R021;
806 case 1:
807 E100 = -uxk * R200;
808 E010 = -uyk * R020;
809 E001 = -uzk * R002;
810 case 0:
811 E000 = uzk * R001;
812 }
813 }
814
815
816
817
818 @SuppressWarnings("fallthrough")
819 @Override
820 protected void quadrupoleKPotentialAtI(PolarizableMultipole mK, int order) {
821 switch (order) {
822 default:
823 case 3:
824
825 E300 = -mK.qxz * R401;
826 E030 = -mK.qyz * R041;
827 double term003 = 0.0;
828 term003 = fma(mK.qxx, R203, term003);
829 term003 = fma(mK.qyy, R023, term003);
830 term003 = fma(mK.qzz, R005, term003);
831 E003 = -term003;
832 E210 = -mK.qyz * R221;
833 double term201 = 0.0;
834 term201 = fma(mK.qxx, R401, term201);
835 term201 = fma(mK.qyy, R221, term201);
836 term201 = fma(mK.qzz, R203, term201);
837 E201 = -term201;
838 E120 = -mK.qxz * R221;
839 double term021 = 0.0;
840 term021 = fma(mK.qxx, R221, term021);
841 term021 = fma(mK.qyy, R041, term021);
842 term021 = fma(mK.qzz, R023, term021);
843 E021 = -term021;
844 E102 = -mK.qxz * R203;
845 E012 = -mK.qyz * R023;
846 E111 = -mK.qxy * R221;
847
848 case 2:
849
850 double term200 = 0.0;
851 term200 = fma(mK.qxx, R400, term200);
852 term200 = fma(mK.qyy, R220, term200);
853 term200 = fma(mK.qzz, R202, term200);
854 E200 = term200;
855 double term020 = 0.0;
856 term020 = fma(mK.qxx, R220, term020);
857 term020 = fma(mK.qyy, R040, term020);
858 term020 = fma(mK.qzz, R022, term020);
859 E020 = term020;
860 double term002 = 0.0;
861 term002 = fma(mK.qxx, R202, term002);
862 term002 = fma(mK.qyy, R022, term002);
863 term002 = fma(mK.qzz, R004, term002);
864 E002 = term002;
865 E110 = mK.qxy * R220;
866 E101 = mK.qxz * R202;
867 E011 = mK.qyz * R022;
868
869 case 1:
870
871 E100 = -mK.qxz * R201;
872 E010 = -mK.qyz * R021;
873 double term001 = 0.0;
874 term001 = fma(mK.qxx, R201, term001);
875 term001 = fma(mK.qyy, R021, term001);
876 term001 = fma(mK.qzz, R003, term001);
877 E001 = -term001;
878 case 0:
879 double term000 = 0.0;
880 term000 = fma(mK.qxx, R200, term000);
881 term000 = fma(mK.qyy, R020, term000);
882 term000 = fma(mK.qzz, R002, term000);
883 E000 = term000;
884 }
885 }
886
887
888
889
890 @Override
891 protected double Tlmnj(final int l, final int m, final int n, final int j, final double[] r, final double[] T000) {
892 double z = r[2];
893 assert (r[0] == 0.0 && r[1] == 0.0);
894
895 if (m == 0 && n == 0) {
896 if (l > 1) {
897 return (l - 1) * Tlmnj(l - 2, 0, 0, j + 1, r, T000);
898 } else if (l == 1) {
899 return 0.0;
900 } else {
901 return T000[j];
902 }
903 } else if (n == 0) {
904 if (m > 1) {
905 return (m - 1) * Tlmnj(l, m - 2, 0, j + 1, r, T000);
906 }
907 return 0.0;
908 } else {
909 if (n > 1) {
910 return z * Tlmnj(l, m, n - 1, j + 1, r, T000)
911 + (n - 1) * Tlmnj(l, m, n - 2, j + 1, r, T000);
912 }
913 return z * Tlmnj(l, m, 0, j + 1, r, T000);
914 }
915 }
916
917
918
919
920
921
922
923
924
925
926
927 @Override
928 protected void noStorageRecursion(double[] r, double[] tensor) {
929 setR(r);
930 noStorageRecursion(tensor);
931 }
932
933
934
935
936
937
938
939
940
941
942 @Override
943 protected void noStorageRecursion(double[] tensor) {
944 assert (x == 0.0 && y == 0.0);
945 double[] r = {x, y, z};
946 source(T000);
947
948 tensor[0] = T000[0];
949
950 for (int l = 1; l <= order; l++) {
951 tensor[ti(l, 0, 0)] = Tlmnj(l, 0, 0, 0, r, T000);
952 }
953
954 for (int l = 0; l <= o1; l++) {
955 for (int m = 1; m <= order - l; m++) {
956 tensor[ti(l, m, 0)] = Tlmnj(l, m, 0, 0, r, T000);
957 }
958 }
959
960 for (int l = 0; l <= o1; l++) {
961 for (int m = 0; m <= o1 - l; m++) {
962 for (int n = 1; n <= order - l - m; n++) {
963 tensor[ti(l, m, n)] = Tlmnj(l, m, n, 0, r, T000);
964 }
965 }
966 }
967 }
968
969
970
971
972 @Override
973 protected void recursion(final double[] r, final double[] tensor) {
974 setR(r);
975 recursion(tensor);
976 }
977
978
979
980
981 @Override
982 protected void recursion(final double[] tensor) {
983 assert (x == 0.0 && y == 0.0);
984 source(work);
985 tensor[0] = work[0];
986
987
988
989
990
991
992 double current;
993 double previous = work[1];
994 tensor[ti(1, 0, 0)] = 0.0;
995 for (int l = 2; l < o1; l++) {
996
997
998
999 current = 0.0;
1000 int iw = il + l - 1;
1001 work[iw] = current;
1002 for (int a = 1; a < l - 1; a++) {
1003
1004
1005
1006
1007 current = a * work[iw - il];
1008 iw += il - 1;
1009 work[iw] = current;
1010 }
1011
1012
1013 tensor[ti(l, 0, 0)] = (l - 1) * previous;
1014 previous = current;
1015 }
1016
1017
1018
1019 for (int l = 0; l < order; l++) {
1020
1021
1022 previous = work[l * il + 1];
1023 tensor[ti(l, 1, 0)] = 0.0;
1024 for (int m = 2; m + l < o1; m++) {
1025
1026 int iw = l * il + m;
1027 current = 0.0;
1028 iw += im - 1;
1029 work[iw] = current;
1030 for (int a = 1; a < m - 1; a++) {
1031
1032
1033
1034
1035 current = a * work[iw - im];
1036 iw += im - 1;
1037 work[iw] = current;
1038 }
1039
1040
1041 tensor[ti(l, m, 0)] = (m - 1) * previous;
1042 previous = current;
1043 }
1044 }
1045
1046
1047
1048 for (int l = 0; l < order; l++) {
1049 for (int m = 0; m + l < order; m++) {
1050
1051
1052 final int lm = m + l;
1053 final int lilmim = l * il + m * im;
1054 previous = work[lilmim + 1];
1055 tensor[ti(l, m, 1)] = z * previous;
1056 for (int n = 2; lm + n < o1; n++) {
1057
1058 int iw = lilmim + n;
1059 current = z * work[iw];
1060 iw += in - 1;
1061 work[iw] = current;
1062 final int n1 = n - 1;
1063 for (int a = 1; a < n1; a++) {
1064
1065
1066
1067
1068 current = z * current + a * work[iw - in];
1069 iw += in - 1;
1070 work[iw] = current;
1071 }
1072
1073
1074 tensor[ti(l, m, n)] = z * current + n1 * previous;
1075 previous = current;
1076 }
1077 }
1078 }
1079 }
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098 @Override
1099 protected String codeTensorRecursion(final double[] r, final double[] tensor) {
1100 setR(r);
1101 assert (x == 0.0 && y == 0.0);
1102 source(work);
1103 StringBuilder sb = new StringBuilder();
1104 tensor[0] = work[0];
1105 if (work[0] > 0) {
1106 sb.append(format("%s = %s;\n", rlmn(0, 0, 0), term(0, 0, 0, 0)));
1107 }
1108
1109
1110
1111
1112
1113
1114 double current;
1115 double previous = work[1];
1116 tensor[ti(1, 0, 0)] = 0.0;
1117 for (int l = 2; l < o1; l++) {
1118
1119
1120
1121 current = 0.0;
1122 int iw = il + (l - 1);
1123 work[iw] = current;
1124
1125 for (int a = 1; a < l - 1; a++) {
1126
1127
1128
1129
1130
1131 current = a * work[iw - il];
1132 iw += il - 1;
1133
1134 work[iw] = current;
1135 if (current != 0) {
1136 if (a > 2) {
1137 sb.append(
1138 format(
1139 "double %s = %d * %s;\n",
1140 term(a + 1, 0, 0, l - a - 1), a, term(a - 1, 0, 0, l - a)));
1141 } else {
1142 sb.append(
1143 format(
1144 "double %s = %s;\n", term(a + 1, 0, 0, l - a - 1), term(a - 1, 0, 0, l - a)));
1145 }
1146 }
1147 }
1148
1149
1150 int index = ti(l, 0, 0);
1151 tensor[index] = (l - 1) * previous;
1152 previous = current;
1153 if (tensor[index] != 0) {
1154 if (l > 2) {
1155 sb.append(format("%s = %d * %s;\n", rlmn(l, 0, 0), (l - 1), term(l - 2, 0, 0, 1)));
1156 } else {
1157 sb.append(format("%s = %s;\n", rlmn(l, 0, 0), term(l - 2, 0, 0, 1)));
1158 }
1159 }
1160 }
1161
1162
1163
1164 for (int l = 0; l < order; l++) {
1165
1166
1167 previous = work[l * il + 1];
1168 tensor[ti(l, 1, 0)] = 0.0;
1169 for (int m = 2; m + l < o1; m++) {
1170
1171 int iw = l * il + m;
1172 current = 0.0;
1173 iw += im - 1;
1174 work[iw] = current;
1175 for (int a = 1; a < m - 1; a++) {
1176
1177
1178
1179
1180 current = a * work[iw - im];
1181 iw += im - 1;
1182 work[iw] = current;
1183 if (current != 0) {
1184 if (a > 1) {
1185 sb.append(
1186 format(
1187 "double %s = %d * %s;\n",
1188 term(l, a + 1, 0, m - a - 1), a, term(l, a - 1, 0, m - a)));
1189 } else {
1190 sb.append(
1191 format(
1192 "double %s = %s;\n", term(l, a + 1, 0, m - a - 1), term(l, a - 1, 0, m - a)));
1193 }
1194 }
1195 }
1196
1197
1198 int index = ti(l, m, 0);
1199 tensor[index] = (m - 1) * previous;
1200 previous = current;
1201 if (tensor[index] != 0) {
1202 if (m > 2) {
1203 sb.append(format("%s = %d * %s;\n", rlmn(l, m, 0), (m - 1), term(l, m - 2, 0, 1)));
1204 } else {
1205 sb.append(format("%s = %s;\n", rlmn(l, m, 0), term(l, m - 2, 0, 1)));
1206 }
1207 }
1208 }
1209 }
1210
1211
1212
1213 for (int l = 0; l < order; l++) {
1214 for (int m = 0; m + l < order; m++) {
1215
1216
1217 final int lm = m + l;
1218 final int lilmim = l * il + m * im;
1219 previous = work[lilmim + 1];
1220 int index = ti(l, m, 1);
1221 tensor[index] = z * previous;
1222 if (tensor[index] != 0) {
1223 sb.append(format("%s = z * %s;\n", rlmn(l, m, 1), term(l, m, 0, 1)));
1224 }
1225 for (int n = 2; lm + n < o1; n++) {
1226
1227 int iw = lilmim + n;
1228 current = z * work[iw];
1229 iw += in - 1;
1230 work[iw] = current;
1231 if (current != 0) {
1232 sb.append(format("double %s = z * %s;\n", term(l, m, 1, n - 1), term(l, m, 0, n)));
1233 }
1234 final int n1 = n - 1;
1235 for (int a = 1; a < n1; a++) {
1236
1237
1238
1239
1240 current = z * current + a * work[iw - in];
1241 iw += in - 1;
1242 work[iw] = current;
1243 if (current != 0) {
1244 if (a > 1) {
1245 sb.append(
1246 format(
1247 "double %s = fma(z, %s, %d * %s);\n",
1248 term(l, m, a + 1, n - a - 1),
1249 term(l, m, a, n - a),
1250 a,
1251 term(l, m, a - 1, n - a)));
1252 } else {
1253 sb.append(
1254 format(
1255 "double %s = fma(z, %s, %s);\n",
1256 term(l, m, a + 1, n - a - 1),
1257 term(l, m, a, n - a),
1258 term(l, m, a - 1, n - a)));
1259 }
1260 }
1261 }
1262
1263
1264 index = ti(l, m, n);
1265 tensor[index] = z * current + n1 * previous;
1266 previous = current;
1267 if (tensor[index] != 0) {
1268 if (n > 2) {
1269 sb.append(
1270 format(
1271 "%s = fma(z, %s, %d * %s);\n",
1272 rlmn(l, m, n), term(l, m, n - 1, 1), (n - 1), term(l, m, n - 2, 1)));
1273 } else {
1274 sb.append(
1275 format(
1276 "%s = fma(z, %s, %s);\n",
1277 rlmn(l, m, n), term(l, m, n - 1, 1), term(l, m, n - 2, 1)));
1278 }
1279 }
1280 }
1281 }
1282 }
1283 return sb.toString();
1284 }
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303 protected String codeVectorTensorRecursion(final double[] r, final double[] tensor) {
1304 setR(r);
1305 assert (x == 0.0 && y == 0.0);
1306 source(work);
1307 StringBuilder sb = new StringBuilder();
1308 tensor[0] = work[0];
1309 if (work[0] > 0) {
1310 sb.append(format("%s = %s;\n", rlmn(0, 0, 0), term(0, 0, 0, 0)));
1311 }
1312
1313
1314
1315
1316
1317
1318 double current;
1319 double previous = work[1];
1320 tensor[ti(1, 0, 0)] = 0.0;
1321 for (int l = 2; l < o1; l++) {
1322
1323
1324
1325 current = 0.0;
1326 int iw = il + (l - 1);
1327 work[iw] = current;
1328
1329 for (int a = 1; a < l - 1; a++) {
1330
1331
1332
1333
1334
1335 current = a * work[iw - il];
1336 iw += il - 1;
1337
1338 work[iw] = current;
1339 if (current != 0) {
1340 if (a > 2) {
1341 sb.append(format("DoubleVector %s = %s.mul(%d);\n",
1342 term(a + 1, 0, 0, l - a - 1), term(a - 1, 0, 0, l - a), a));
1343 } else {
1344 sb.append(format("DoubleVector %s = %s;\n", term(a + 1, 0, 0, l - a - 1), term(a - 1, 0, 0, l - a)));
1345 }
1346 }
1347 }
1348
1349
1350 int index = ti(l, 0, 0);
1351 tensor[index] = (l - 1) * previous;
1352 previous = current;
1353 if (tensor[index] != 0) {
1354 if (l > 2) {
1355 sb.append(format("%s = %s.mul(%d);\n", rlmn(l, 0, 0), term(l - 2, 0, 0, 1), (l - 1)));
1356 } else {
1357 sb.append(format("%s = %s;\n", rlmn(l, 0, 0), term(0, 0, 0, 1)));
1358 }
1359 }
1360 }
1361
1362
1363
1364 for (int l = 0; l < order; l++) {
1365
1366
1367 previous = work[l * il + 1];
1368 tensor[ti(l, 1, 0)] = 0.0;
1369 for (int m = 2; m + l < o1; m++) {
1370
1371 int iw = l * il + m;
1372 current = 0.0;
1373 iw += im - 1;
1374 work[iw] = current;
1375 for (int a = 1; a < m - 1; a++) {
1376
1377
1378
1379
1380 current = a * work[iw - im];
1381 iw += im - 1;
1382 work[iw] = current;
1383 if (current != 0) {
1384 if (a > 1) {
1385 sb.append(format("DoubleVector %s = %s.mul(%d);\n",
1386 term(l, a + 1, 0, m - a - 1), term(l, a - 1, 0, m - a), a));
1387 } else {
1388 sb.append(format("DoubleVector %s = %s;\n",
1389 term(l, a + 1, 0, m - a - 1), term(l, 0, 0, m - a)));
1390 }
1391 }
1392 }
1393
1394
1395 int index = ti(l, m, 0);
1396 tensor[index] = (m - 1) * previous;
1397 previous = current;
1398 if (tensor[index] != 0) {
1399 if (m > 2) {
1400 sb.append(format("%s = %s.mul(%d);\n", rlmn(l, m, 0), term(l, m - 2, 0, 1), (m - 1)));
1401 } else {
1402 sb.append(format("%s = %s;\n", rlmn(l, m, 0), term(l, m - 2, 0, 1)));
1403 }
1404 }
1405 }
1406 }
1407
1408
1409
1410 for (int l = 0; l < order; l++) {
1411 for (int m = 0; m + l < order; m++) {
1412
1413
1414 final int lm = m + l;
1415 final int lilmim = l * il + m * im;
1416 previous = work[lilmim + 1];
1417 int index = ti(l, m, 1);
1418 tensor[index] = z * previous;
1419 if (tensor[index] != 0) {
1420 sb.append(format("%s = z.mul(%s);\n", rlmn(l, m, 1), term(l, m, 0, 1)));
1421 }
1422 for (int n = 2; lm + n < o1; n++) {
1423
1424 int iw = lilmim + n;
1425 current = z * work[iw];
1426 iw += in - 1;
1427 work[iw] = current;
1428 if (current != 0) {
1429 sb.append(format("DoubleVector %s = z.mul(%s);\n", term(l, m, 1, n - 1), term(l, m, 0, n)));
1430 }
1431 final int n1 = n - 1;
1432 for (int a = 1; a < n1; a++) {
1433
1434
1435
1436
1437 current = z * current + a * work[iw - in];
1438 iw += in - 1;
1439 work[iw] = current;
1440 if (current != 0) {
1441 if (a > 1) {
1442 sb.append(format("DoubleVector %s = z.fma(%s, %s.mul(%d));\n",
1443 term(l, m, a + 1, n - a - 1), term(l, m, a, n - a),
1444 term(l, m, a - 1, n - a), a));
1445 } else {
1446 sb.append(format("DoubleVector %s = z.fma(%s, %s);\n",
1447 term(l, m, a + 1, n - a - 1), term(l, m, a, n - a),
1448 term(l, m, 0, n - a)));
1449 }
1450 }
1451 }
1452
1453
1454 index = ti(l, m, n);
1455 tensor[index] = z * current + n1 * previous;
1456 previous = current;
1457 if (tensor[index] != 0) {
1458 if (n > 2) {
1459 sb.append(format("%s = z.fma(%s, %s.mul(%d));\n",
1460 rlmn(l, m, n), term(l, m, n - 1, 1), term(l, m, n - 2, 1), (n - 1)));
1461 } else {
1462 sb.append(format("%s = z.fma(%s, %s);\n",
1463 rlmn(l, m, n), term(l, m, n - 1, 1), term(l, m, n - 2, 1)));
1464 }
1465 }
1466 }
1467 }
1468 }
1469 return sb.toString();
1470 }
1471
1472 }