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 GKTensorQI extends CoulombTensorQI {
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 GKTensorQI(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 multipoleKPotentialAtI(mK, 1);
201 dipoleTorque(mI, Ti);
202
203
204 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
205 double eI = multipoleEnergy(mI);
206 multipoleGradient(mI, Gi);
207 multipoleTorque(mI, Ti);
208
209
210 multipoleIPotentialAtK(mI, 1);
211 dipoleTorque(mK, Tk);
212
213 double scale = c * 0.5;
214 Gi[0] = scale * (Gi[0] - Gk[0]);
215 Gi[1] = scale * (Gi[1] - Gk[1]);
216 Gi[2] = scale * (Gi[2] - Gk[2]);
217 Gk[0] = -Gi[0];
218 Gk[1] = -Gi[1];
219 Gk[2] = -Gi[2];
220
221 Ti[0] = scale * Ti[0];
222 Ti[1] = scale * Ti[1];
223 Ti[2] = scale * Ti[2];
224 Tk[0] = scale * Tk[0];
225 Tk[1] = scale * Tk[1];
226 Tk[2] = scale * Tk[2];
227
228 return scale * (eK + eI);
229 }
230
231
232
233
234
235
236
237
238
239
240
241
242 protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
243 double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
244
245
246 quadrupoleIPotentialAtK(mI, 3);
247 double eK = multipoleEnergy(mK);
248 multipoleGradient(mK, Gk);
249 multipoleTorque(mK, Tk);
250
251
252 multipoleKPotentialAtI(mK, 2);
253 quadrupoleTorque(mI, Ti);
254
255
256 quadrupoleKPotentialAtI(mK, 3);
257 double eI = multipoleEnergy(mI);
258 multipoleGradient(mI, Gi);
259 multipoleTorque(mI, Ti);
260
261
262 multipoleIPotentialAtK(mI, 2);
263 quadrupoleTorque(mK, Tk);
264
265 double scale = c * 0.5;
266 Gi[0] = scale * (Gi[0] - Gk[0]);
267 Gi[1] = scale * (Gi[1] - Gk[1]);
268 Gi[2] = scale * (Gi[2] - Gk[2]);
269 Gk[0] = -Gi[0];
270 Gk[1] = -Gi[1];
271 Gk[2] = -Gi[2];
272
273 Ti[0] = scale * Ti[0];
274 Ti[1] = scale * Ti[1];
275 Ti[2] = scale * Ti[2];
276 Tk[0] = scale * Tk[0];
277 Tk[1] = scale * Tk[1];
278 Tk[2] = scale * Tk[2];
279
280 return scale * (eK + eI);
281 }
282
283
284
285
286
287
288
289
290 public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
291 generateTensor();
292 return multipoleEnergy(mI, mK);
293 }
294
295
296
297
298
299
300
301
302
303
304 @Override
305 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
306 double scaleEnergy) {
307 return polarizationEnergy(mI, mK);
308 }
309
310
311
312
313
314
315
316
317 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
318 return switch (multipoleOrder) {
319 default -> {
320
321 chargeIPotentialAtK(mI, 1);
322
323 double eK = polarizationEnergy(mK);
324
325 chargeKPotentialAtI(mK, 1);
326
327 double eI = polarizationEnergy(mI);
328 yield c * 0.5 * (eK + eI);
329 }
330 case DIPOLE -> {
331
332 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
333
334 double eK = polarizationEnergy(mK);
335
336 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
337
338 eK += 0.5 * multipoleEnergy(mK);
339
340 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
341
342 double eI = polarizationEnergy(mI);
343
344 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
345
346 eI += 0.5 * multipoleEnergy(mI);
347 yield c * 0.5 * (eK + eI);
348 }
349 case QUADRUPOLE -> {
350
351 quadrupoleIPotentialAtK(mI, 1);
352
353 double eK = polarizationEnergy(mK);
354
355 quadrupoleKPotentialAtI(mK, 1);
356
357 double eI = polarizationEnergy(mI);
358 yield c * 0.5 * (eK + eI);
359 }
360 };
361 }
362
363
364
365
366
367
368
369
370 public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
371 return switch (multipoleOrder) {
372 default -> {
373
374 chargeIPotentialAtK(mI, 1);
375
376 double eK = polarizationEnergyS(mK);
377
378 chargeKPotentialAtI(mK, 1);
379
380 double eI = polarizationEnergyS(mI);
381 yield c * 0.5 * (eK + eI);
382 }
383 case DIPOLE -> {
384
385 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
386
387 double eK = polarizationEnergyS(mK);
388
389 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
390
391 eK += 0.5 * multipoleEnergy(mK);
392
393 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
394
395 double eI = polarizationEnergyS(mI);
396
397 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
398
399 eI += 0.5 * multipoleEnergy(mI);
400 yield c * 0.5 * (eK + eI);
401 }
402 case QUADRUPOLE -> {
403
404 quadrupoleIPotentialAtK(mI, 1);
405
406 double eK = polarizationEnergyS(mK);
407
408 quadrupoleKPotentialAtI(mK, 1);
409
410 double eI = polarizationEnergyS(mI);
411 yield c * 0.5 * (eK + eI);
412 }
413 };
414 }
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431 @Override
432 public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
433 double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
434 double[] Tk) {
435 return switch (multipoleOrder) {
436 default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
437 case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
438 case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
439 };
440 }
441
442
443
444
445
446
447
448
449
450 public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
451 PolarizableMultipole mK, double[] Gi) {
452
453
454 chargeIPotentialAtK(mI, 2);
455
456 double eK = polarizationEnergy(mK);
457
458 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
459 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
460 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
461
462
463 chargeKPotentialAtI(mK, 2);
464
465 double eI = polarizationEnergy(mI);
466
467 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
468 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
469 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
470
471 double scale = c * 0.5;
472 Gi[0] *= scale;
473 Gi[1] *= scale;
474 Gi[2] *= scale;
475
476
477 return scale * (eI + eK);
478 }
479
480
481
482
483
484
485
486
487
488
489
490
491 public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
492 double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
493
494
495 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
496
497 double eK = polarizationEnergy(mK);
498
499 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
500 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
501 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
502
503 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
504 dipoleTorque(mI, Ti);
505
506
507 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
508
509 eK += 0.5 * multipoleEnergy(mK);
510
511 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 3);
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 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
529 dipoleTorque(mK, Tk);
530
531
532 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
533
534 eI += 0.5 * multipoleEnergy(mI);
535
536 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 3);
537 G = new double[3];
538 multipoleGradient(mI, G);
539 Gi[0] += G[0];
540 Gi[1] += G[1];
541 Gi[2] += G[2];
542 multipoleTorque(mI, Ti);
543
544
545
546 if (mutualMask != 0.0) {
547
548 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
549 Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
550 Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
551 Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
552
553
554 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
555 Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
556 Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
557 Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
558 }
559
560
561 double scale = c * 0.5;
562 double energy = scale * (eI + eK);
563 Gi[0] *= scale;
564 Gi[1] *= scale;
565 Gi[2] *= scale;
566 Ti[0] *= scale;
567 Ti[1] *= scale;
568 Ti[2] *= scale;
569 Tk[0] *= scale;
570 Tk[1] *= scale;
571 Tk[2] *= scale;
572
573 return energy;
574 }
575
576
577
578
579
580
581
582
583
584
585
586 public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI,
587 PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {
588
589
590 quadrupoleIPotentialAtK(mI, 2);
591
592 double eK = polarizationEnergy(mK);
593
594 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
595 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
596 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
597
598
599 quadrupoleKPotentialAtI(mK, 2);
600
601 double eI = polarizationEnergy(mI);
602
603 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
604 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
605 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
606
607 double scale = c * 0.5;
608 Gi[0] *= scale;
609 Gi[1] *= scale;
610 Gi[2] *= scale;
611
612
613 dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
614 quadrupoleTorque(mK, Tk);
615
616
617 dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
618 quadrupoleTorque(mI, Ti);
619
620
621 return scale * (eI + eK);
622 }
623
624
625
626
627
628
629
630
631 public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
632 generateTensor();
633 return 2.0 * polarizationEnergyBorn(mI, mK);
634 }
635
636
637
638
639
640
641
642
643 public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
644 double db = 0.0;
645 if (multipoleOrder == GK_MULTIPOLE_ORDER.DIPOLE) {
646
647 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
648 db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);
649
650
651 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
652 db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
653 }
654 return c * db;
655 }
656
657
658
659
660 @Override
661 protected void source(double[] work) {
662 gkSource.source(work, multipoleOrder);
663 }
664 }