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 ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER;
41
42
43
44
45
46
47
48
49 public class GKTensorGlobal extends CoulombTensorGlobal {
50
51
52
53
54 protected final GK_MULTIPOLE_ORDER multipoleOrder;
55
56
57
58
59 private final double c;
60
61 private final GKSource gkSource;
62
63
64
65
66
67
68
69
70 public GKTensorGlobal(GK_MULTIPOLE_ORDER multipoleOrder, int order, GKSource gkSource, double Eh,
71 double Es) {
72 super(order);
73 this.multipoleOrder = multipoleOrder;
74 this.gkSource = gkSource;
75
76
77 c = GKSource.cn(multipoleOrder.getOrder(), Eh, Es);
78 }
79
80
81
82
83
84
85
86
87 @Override
88 public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
89 return switch (multipoleOrder) {
90 default -> {
91 chargeIPotentialAtK(mI, 2);
92 double eK = multipoleEnergy(mK);
93 chargeKPotentialAtI(mK, 2);
94 double eI = multipoleEnergy(mI);
95 yield c * 0.5 * (eK + eI);
96 }
97 case DIPOLE -> {
98 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
99 double eK = multipoleEnergy(mK);
100 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
101 double eI = multipoleEnergy(mI);
102 yield c * 0.5 * (eK + eI);
103 }
104 case QUADRUPOLE -> {
105 quadrupoleIPotentialAtK(mI, 2);
106 double eK = multipoleEnergy(mK);
107 quadrupoleKPotentialAtI(mK, 2);
108 double eI = multipoleEnergy(mI);
109 yield c * 0.5 * (eK + eI);
110 }
111 };
112 }
113
114
115
116
117
118
119
120
121
122
123
124
125 @Override
126 public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
127 double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
128 return switch (multipoleOrder) {
129 default -> monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
130 case DIPOLE -> dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
131 case QUADRUPOLE -> quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
132 };
133 }
134
135
136
137
138
139
140
141
142
143
144
145
146 protected double monopoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
147 double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
148
149
150 chargeIPotentialAtK(mI, 3);
151 double eK = multipoleEnergy(mK);
152 multipoleGradient(mK, Gk);
153 multipoleTorque(mK, Tk);
154
155
156 chargeKPotentialAtI(mK, 3);
157 double eI = multipoleEnergy(mI);
158 multipoleGradient(mI, Gi);
159 multipoleTorque(mI, Ti);
160
161 double scale = c * 0.5;
162 Gi[0] = scale * (Gi[0] - Gk[0]);
163 Gi[1] = scale * (Gi[1] - Gk[1]);
164 Gi[2] = scale * (Gi[2] - Gk[2]);
165 Gk[0] = -Gi[0];
166 Gk[1] = -Gi[1];
167 Gk[2] = -Gi[2];
168
169 Ti[0] = scale * Ti[0];
170 Ti[1] = scale * Ti[1];
171 Ti[2] = scale * Ti[2];
172 Tk[0] = scale * Tk[0];
173 Tk[1] = scale * Tk[1];
174 Tk[2] = scale * Tk[2];
175
176 return scale * (eK + eI);
177 }
178
179
180
181
182
183
184
185
186
187
188
189
190 protected double dipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
191 double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
192
193
194 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
195 double eK = multipoleEnergy(mK);
196 multipoleGradient(mK, Gk);
197 multipoleTorque(mK, Tk);
198
199
200
201 multipoleKPotentialAtI(mK, 1);
202 dipoleTorque(mI, Ti);
203
204
205 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
206 double eI = multipoleEnergy(mI);
207 multipoleGradient(mI, Gi);
208 multipoleTorque(mI, Ti);
209
210
211
212 multipoleIPotentialAtK(mI, 1);
213 dipoleTorque(mK, Tk);
214
215 double scale = c * 0.5;
216 Gi[0] = scale * (Gi[0] - Gk[0]);
217 Gi[1] = scale * (Gi[1] - Gk[1]);
218 Gi[2] = scale * (Gi[2] - Gk[2]);
219 Gk[0] = -Gi[0];
220 Gk[1] = -Gi[1];
221 Gk[2] = -Gi[2];
222
223 Ti[0] = scale * Ti[0];
224 Ti[1] = scale * Ti[1];
225 Ti[2] = scale * Ti[2];
226 Tk[0] = scale * Tk[0];
227 Tk[1] = scale * Tk[1];
228 Tk[2] = scale * Tk[2];
229
230 return scale * (eK + eI);
231 }
232
233
234
235
236
237
238
239
240
241
242
243
244 protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
245 double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
246
247
248 quadrupoleIPotentialAtK(mI, 3);
249 double eK = multipoleEnergy(mK);
250 multipoleGradient(mK, Gk);
251 multipoleTorque(mK, Tk);
252
253
254 multipoleKPotentialAtI(mK, 2);
255 quadrupoleTorque(mI, Ti);
256
257
258 quadrupoleKPotentialAtI(mK, 3);
259 double eI = multipoleEnergy(mI);
260 multipoleGradient(mI, Gi);
261 multipoleTorque(mI, Ti);
262
263
264 multipoleIPotentialAtK(mI, 2);
265 quadrupoleTorque(mK, Tk);
266
267 double scale = c * 0.5;
268 Gi[0] = scale * (Gi[0] - Gk[0]);
269 Gi[1] = scale * (Gi[1] - Gk[1]);
270 Gi[2] = scale * (Gi[2] - Gk[2]);
271 Gk[0] = -Gi[0];
272 Gk[1] = -Gi[1];
273 Gk[2] = -Gi[2];
274
275 Ti[0] = scale * Ti[0];
276 Ti[1] = scale * Ti[1];
277 Ti[2] = scale * Ti[2];
278 Tk[0] = scale * Tk[0];
279 Tk[1] = scale * Tk[1];
280 Tk[2] = scale * Tk[2];
281
282 return scale * (eK + eI);
283 }
284
285
286
287
288
289
290
291
292 public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
293 generateTensor();
294 return multipoleEnergy(mI, mK);
295 }
296
297
298
299
300
301
302
303
304
305
306 @Override
307 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
308 double scaleEnergy) {
309 return polarizationEnergy(mI, mK);
310 }
311
312
313
314
315
316
317
318
319 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
320 return switch (multipoleOrder) {
321 default -> {
322
323 chargeIPotentialAtK(mI, 1);
324
325 double eK = polarizationEnergy(mK);
326
327 chargeKPotentialAtI(mK, 1);
328
329 double eI = polarizationEnergy(mI);
330 yield c * 0.5 * (eK + eI);
331 }
332 case DIPOLE -> {
333
334 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
335
336 double eK = polarizationEnergy(mK);
337
338 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
339
340 eK += 0.5 * multipoleEnergy(mK);
341
342 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
343
344 double eI = polarizationEnergy(mI);
345
346 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
347
348 eI += 0.5 * multipoleEnergy(mI);
349 yield c * 0.5 * (eK + eI);
350 }
351 case QUADRUPOLE -> {
352
353 quadrupoleIPotentialAtK(mI, 1);
354
355 double eK = polarizationEnergy(mK);
356
357 quadrupoleKPotentialAtI(mK, 1);
358
359 double eI = polarizationEnergy(mI);
360 yield c * 0.5 * (eK + eI);
361 }
362 };
363 }
364
365
366
367
368
369
370
371
372 public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
373 return switch (multipoleOrder) {
374 default -> {
375
376 chargeIPotentialAtK(mI, 1);
377
378 double eK = polarizationEnergyS(mK);
379
380 chargeKPotentialAtI(mK, 1);
381
382 double eI = polarizationEnergyS(mI);
383 yield c * 0.5 * (eK + eI);
384 }
385 case DIPOLE -> {
386
387 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
388
389 double eK = polarizationEnergyS(mK);
390
391 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
392
393 eK += 0.5 * multipoleEnergy(mK);
394
395 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
396
397 double eI = polarizationEnergyS(mI);
398
399 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
400
401 eI += 0.5 * multipoleEnergy(mI);
402 yield c * 0.5 * (eK + eI);
403 }
404 case QUADRUPOLE -> {
405
406 quadrupoleIPotentialAtK(mI, 1);
407
408 double eK = polarizationEnergyS(mK);
409
410 quadrupoleKPotentialAtI(mK, 1);
411
412 double eI = polarizationEnergyS(mI);
413 yield c * 0.5 * (eK + eI);
414 }
415 };
416 }
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433 @Override
434 public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
435 double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
436 double[] Tk) {
437 return switch (multipoleOrder) {
438 default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
439 case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
440 case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
441 };
442 }
443
444
445
446
447
448
449
450
451
452 public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
453 PolarizableMultipole mK, double[] Gi) {
454
455
456 chargeIPotentialAtK(mI, 2);
457
458 double eK = polarizationEnergy(mK);
459
460 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
461 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
462 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
463
464
465 chargeKPotentialAtI(mK, 2);
466
467 double eI = polarizationEnergy(mI);
468
469 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
470 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
471 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
472
473 double scale = c * 0.5;
474 Gi[0] *= scale;
475 Gi[1] *= scale;
476 Gi[2] *= scale;
477
478
479 return scale * (eI + eK);
480 }
481
482
483
484
485
486
487
488
489
490
491
492
493 public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
494 double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
495
496
497 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
498
499 double eK = polarizationEnergy(mK);
500
501 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
502 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
503 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
504
505 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
506 dipoleTorque(mK, Tk);
507
508
509 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 3);
510
511 eK += 0.5 * multipoleEnergy(mK);
512 double[] G = new double[3];
513 multipoleGradient(mK, G);
514 Gi[0] -= G[0];
515 Gi[1] -= G[1];
516 Gi[2] -= G[2];
517 multipoleTorque(mK, Tk);
518
519
520 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
521
522 double eI = polarizationEnergy(mI);
523
524 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
525 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
526 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
527
528 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
529 dipoleTorque(mI, Ti);
530
531
532 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 3);
533
534 eI += 0.5 * multipoleEnergy(mI);
535 G = new double[3];
536 multipoleGradient(mI, G);
537 Gi[0] += G[0];
538 Gi[1] += G[1];
539 Gi[2] += G[2];
540 multipoleTorque(mI, Ti);
541
542
543
544 if (mutualMask != 0.0) {
545
546 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
547 Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
548 Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
549 Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
550
551
552 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
553 Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
554 Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
555 Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
556 }
557
558
559 double scale = c * 0.5;
560 double energy = scale * (eI + eK);
561 Gi[0] *= scale;
562 Gi[1] *= scale;
563 Gi[2] *= scale;
564 Ti[0] *= scale;
565 Ti[1] *= scale;
566 Ti[2] *= scale;
567 Tk[0] *= scale;
568 Tk[1] *= scale;
569 Tk[2] *= scale;
570
571 return energy;
572 }
573
574
575
576
577
578
579
580
581
582
583
584 public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI,
585 PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {
586
587
588 quadrupoleIPotentialAtK(mI, 2);
589
590 double eK = polarizationEnergy(mK);
591
592 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
593 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
594 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
595
596
597 quadrupoleKPotentialAtI(mK, 2);
598
599 double eI = polarizationEnergy(mI);
600
601 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
602 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
603 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
604
605 double scale = c * 0.5;
606 Gi[0] *= scale;
607 Gi[1] *= scale;
608 Gi[2] *= scale;
609
610
611 dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
612 quadrupoleTorque(mK, Tk);
613
614
615 dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
616 quadrupoleTorque(mI, Ti);
617
618
619 return scale * (eI + eK);
620 }
621
622
623
624
625
626
627
628
629 public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
630 generateTensor();
631 return 2.0 * polarizationEnergyBorn(mI, mK);
632 }
633
634
635
636
637
638
639
640
641 public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
642 double db = 0.0;
643 if (multipoleOrder == GK_MULTIPOLE_ORDER.DIPOLE) {
644
645 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
646 db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);
647
648
649 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
650 db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
651 }
652 return c * db;
653 }
654
655
656
657
658 @Override
659 protected void source(double[] work) {
660 gkSource.source(work, multipoleOrder);
661 }
662
663 }