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