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.math.ScalarMath.doubleFactorial;
41 import static java.lang.System.arraycopy;
42 import static java.util.Arrays.fill;
43 import static org.apache.commons.math3.util.FastMath.exp;
44 import static org.apache.commons.math3.util.FastMath.pow;
45 import static org.apache.commons.math3.util.FastMath.sqrt;
46
47 public class GKSource {
48
49
50
51
52 public enum GK_TENSOR_MODE {POTENTIAL, BORN}
53
54
55
56
57
58 public enum GK_MULTIPOLE_ORDER {
59 MONOPOLE(0), DIPOLE(1), QUADRUPOLE(2);
60
61 private final int order;
62
63 GK_MULTIPOLE_ORDER(int order) {
64 this.order = order;
65 }
66
67 public int getOrder() {
68 return order;
69 }
70 }
71
72 private GK_TENSOR_MODE mode = GK_TENSOR_MODE.POTENTIAL;
73
74
75
76
77 private double rb2;
78
79
80
81
82 private double expTerm;
83
84
85
86
87 private double f;
88
89
90
91
92 private double f1;
93
94
95
96
97 private double f2;
98
99
100
101
102 private final double gc;
103
104
105
106
107 private final double igc;
108
109
110
111
112 private double gcAiAj;
113
114
115
116
117 private double ratio;
118
119
120
121
122 private double fr;
123
124
125
126
127 private double r2;
128
129
130
131
132 private final int order;
133
134
135
136
137 private final double[] kirkwoodSource;
138
139
140
141
142 private final double[][] anmc;
143
144
145
146
147
148 private final double[] fn;
149
150
151
152
153
154 protected final double[] bn;
155
156
157
158
159 private final double[][] anm;
160
161
162
163
164
165 private final double[][] bnm;
166
167 public GKSource(int order, double gc) {
168 this.order = order;
169 this.gc = gc;
170 this.igc = 1.0 / gc;
171
172
173 kirkwoodSource = new double[order + 1];
174 for (short n = 0; n <= order; n++) {
175 kirkwoodSource[n] = pow(-1, n) * doubleFactorial(2 * n - 1);
176 }
177
178 anmc = new double[order + 1][];
179 for (int n = 0; n <= order; n++) {
180 anmc[n] = anmc(n);
181 }
182
183 fn = new double[order + 1];
184 anm = new double[order + 1][order + 1];
185 bn = new double[order + 1];
186 bnm = new double[order + 1][order + 1];
187 }
188
189
190
191
192 protected void source(double[] work, GK_MULTIPOLE_ORDER multipoleOrder) {
193 int mpoleOrder = multipoleOrder.getOrder();
194 fill(work, 0, mpoleOrder, 0.0);
195 if (mode == GK_TENSOR_MODE.POTENTIAL) {
196
197 int derivatives = order - mpoleOrder;
198 arraycopy(anm[mpoleOrder], 0, work, mpoleOrder, derivatives + 1);
199 } else {
200
201 int derivatives = (order - 1) - mpoleOrder;
202 arraycopy(bnm[mpoleOrder], 0, work, mpoleOrder, derivatives + 1);
203 }
204 }
205
206
207
208
209 public void generateSource(GK_TENSOR_MODE mode, GK_MULTIPOLE_ORDER multipole, double r2,
210 double ai, double aj) {
211 int multipoleOrder = multipole.getOrder();
212 this.mode = mode;
213
214 this.r2 = r2;
215 rb2 = ai * aj;
216 gcAiAj = gc * rb2;
217 ratio = -r2 / gcAiAj;
218 expTerm = exp(ratio);
219 fr = -2.0 / gcAiAj;
220
221 f = sqrt(r2 + rb2 * expTerm);
222 f1 = 1.0 - expTerm * igc;
223 f2 = 2.0 * expTerm / (gc * gcAiAj);
224
225 if (mode == GK_TENSOR_MODE.POTENTIAL) {
226
227 anm(order, order - multipoleOrder);
228 } else {
229
230 bnm(order - 1, (order - 1) - multipoleOrder);
231 }
232 }
233
234
235
236
237
238
239
240
241
242
243
244 private void anm(int n, int derivatives) {
245
246 fn(n);
247
248
249 an0(n);
250
251
252 for (int d = 1; d <= derivatives; d++) {
253
254 var coef = anmc[d];
255
256 int limit = n - d;
257
258 for (int order = 0; order <= limit; order++) {
259
260 var terms = anm[order + 1];
261 var sum = 0.0;
262 for (int i = 1; i <= d; i++) {
263 sum += coef[i - 1] * fn[i] * terms[d - i];
264 }
265 anm[order][d] = sum;
266 }
267 }
268 }
269
270
271
272
273
274
275
276
277
278
279 private void bnm(int n, int derivatives) {
280
281 bn(n);
282
283
284 bn0(n);
285
286
287 for (int d = 1; d <= derivatives; d++) {
288
289 var coef = anmc[d];
290
291 int limit = n - d;
292
293 for (int order = 0; order <= limit; order++) {
294
295 var terma = anm[order + 1];
296 var termb = bnm[order + 1];
297 var sum = 0.0;
298 for (int i = 1; i <= d; i++) {
299 sum += coef[i - 1] * bn[i] * terma[d - i];
300 sum += coef[i - 1] * fn[i] * termb[d - i];
301 }
302 bnm[order][d] = sum;
303 }
304 }
305 }
306
307
308
309
310
311
312
313 private void fn(int n) {
314 switch (n) {
315 case 0 -> fn[0] = f;
316 case 1 -> {
317 fn[0] = f;
318 fn[1] = f1;
319 }
320 case 2 -> {
321 fn[0] = f;
322 fn[1] = f1;
323 fn[2] = f2;
324 }
325 case 3 -> {
326 fn[0] = f;
327 fn[1] = f1;
328 fn[2] = f2;
329 fn[3] = fr * f2;
330 }
331 case 4 -> {
332 fn[0] = f;
333 fn[1] = f1;
334 fn[2] = f2;
335 fn[3] = fr * f2;
336 fn[4] = fr * fn[3];
337 }
338 case 5 -> {
339 fn[0] = f;
340 fn[1] = f1;
341 fn[2] = f2;
342 fn[3] = fr * f2;
343 fn[4] = fr * fn[3];
344 fn[5] = fr * fn[4];
345 }
346 case 6 -> {
347 fn[0] = f;
348 fn[1] = f1;
349 fn[2] = f2;
350 fn[3] = fr * f2;
351 fn[4] = fr * fn[3];
352 fn[5] = fr * fn[4];
353 fn[6] = fr * fn[5];
354 }
355 default -> {
356 fn[0] = f;
357 fn[1] = f1;
358 fn[2] = f2;
359 for (int i = 3; i <= n; i++) {
360 fn[i] = fr * fn[i - 1];
361 }
362 }
363 }
364 }
365
366
367
368
369
370
371
372 protected void bn(int n) {
373 var b2 = 2.0 * expTerm / (gcAiAj * gcAiAj) * (-ratio - 1.0);
374 switch (n) {
375 case 0 -> bn[0] = 0.5 * expTerm * (1.0 - ratio);
376 case 1 -> {
377 bn[0] = 0.5 * expTerm * (1.0 - ratio);
378 bn[1] = -r2 * expTerm / (gcAiAj * gcAiAj);
379 }
380 case 2 -> {
381 bn[0] = 0.5 * expTerm * (1.0 - ratio);
382 bn[1] = -r2 * expTerm / (gcAiAj * gcAiAj);
383 bn[2] = b2;
384 }
385 default -> {
386 bn[0] = 0.5 * expTerm * (1.0 - ratio);
387 bn[1] = -r2 * expTerm / (gcAiAj * gcAiAj);
388 bn[2] = b2;
389 var br = 2.0 / (gcAiAj * rb2);
390 var f2 = 2.0 / (gc * gcAiAj) * expTerm;
391 var frA = 1.0;
392 var frB = fr;
393 for (int i = 3; i < n; i++) {
394
395 bn[i] = (i - 2) * frA * br * f2 + frB * b2;
396 frA *= fr;
397 frB *= fr;
398 }
399 }
400 }
401 }
402
403
404
405
406
407
408 private void an0(int n) {
409 double inverseF = 1.0 / f;
410 double inverseF2 = inverseF * inverseF;
411 for (int i = 0; i <= n; i++) {
412 anm[i][0] = kirkwoodSource[i] * inverseF;
413 inverseF *= inverseF2;
414 }
415 }
416
417
418
419
420
421
422 private void bn0(int n) {
423 for (int i = 0; i <= n; i++) {
424 bnm[i][0] = bn[0] * anm[i + 1][0];
425 }
426 }
427
428
429
430
431
432
433
434 protected static double[] anmc(int n) {
435 double[] ret = new double[n + 1];
436 ret[0] = 1.0;
437 switch (n) {
438 case 0 -> {
439 return ret;
440 }
441 case 1 -> {
442 ret[1] = 1.0;
443 return ret;
444 }
445 default -> {
446 ret[1] = 1.0;
447 var prev = new double[n];
448 prev[0] = 1.0;
449 prev[1] = 1.0;
450 for (int i = 3; i <= n; i++) {
451 for (int j = 2; j <= i - 1; j++) {
452 ret[j - 1] = prev[j - 2] + prev[j - 1];
453 }
454 ret[i - 1] = 1.0;
455 arraycopy(ret, 0, prev, 0, i);
456 }
457 return ret;
458 }
459 }
460 }
461
462
463
464
465
466
467
468
469
470 public static double cn(int n, double Eh, double Es) {
471 var ret = (n + 1) * (Eh - Es) / ((n + 1) * Es + n * Eh);
472 return ret / Eh;
473 }
474
475 public static double selfEnergy(PolarizableMultipole polarizableMultipole,
476 double ai, double Eh, double Es) {
477 double q2 = polarizableMultipole.q * polarizableMultipole.q;
478 double dx = polarizableMultipole.dx;
479 double dy = polarizableMultipole.dy;
480 double dz = polarizableMultipole.dz;
481 double dx2 = dx * dx;
482 double dy2 = dy * dy;
483 double dz2 = dz * dz;
484 double ux = polarizableMultipole.ux;
485 double uy = polarizableMultipole.uy;
486 double uz = polarizableMultipole.uz;
487 double qxy2 = polarizableMultipole.qxy * polarizableMultipole.qxy;
488 double qxz2 = polarizableMultipole.qxz * polarizableMultipole.qxz;
489 double qyz2 = polarizableMultipole.qyz * polarizableMultipole.qyz;
490 double qxx2 = polarizableMultipole.qxx * polarizableMultipole.qxx;
491 double qyy2 = polarizableMultipole.qyy * polarizableMultipole.qyy;
492 double qzz2 = polarizableMultipole.qzz * polarizableMultipole.qzz;
493
494 double a2 = ai * ai;
495 double a3 = ai * a2;
496 double a5 = a2 * a3;
497
498
499 double e0 = cn(0, Eh, Es) * q2 / ai;
500
501 double e1 = cn(1, Eh, Es) * (dx2 + dy2 + dz2) / a3;
502
503 double e2 = cn(2, Eh, Es) * (3.0 * (qxy2 + qxz2 + qyz2) + 6.0 * (qxx2 + qyy2 + qzz2)) / a5;
504
505
506 double ei = cn(1, Eh, Es) * (dx * ux + dy * uy + dz * uz) / a3;
507
508 return 0.5 * (e0 + e1 + e2 + ei);
509 }
510
511 }