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
41
42
43
44
45
46
47 public class GKTensorQI extends CoulombTensorQI {
48
49
50
51
52 protected final GKMultipoleOrder multipoleOrder;
53
54
55
56
57 private final double c;
58
59 private final GKSource gkSource;
60
61
62
63
64
65
66
67
68
69
70 public GKTensorQI(GKMultipoleOrder 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 return multipoleEnergy(mI, mK);
292 }
293
294
295
296
297
298
299
300
301
302
303 @Override
304 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
305 double scaleEnergy) {
306 return polarizationEnergy(mI, mK);
307 }
308
309
310
311
312
313
314
315
316 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
317 return switch (multipoleOrder) {
318 default -> {
319
320 chargeIPotentialAtK(mI, 1);
321
322 double eK = polarizationEnergy(mK);
323
324 chargeKPotentialAtI(mK, 1);
325
326 double eI = polarizationEnergy(mI);
327 yield c * 0.5 * (eK + eI);
328 }
329 case DIPOLE -> {
330
331 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
332
333 double eK = polarizationEnergy(mK);
334
335 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
336
337 eK += 0.5 * multipoleEnergy(mK);
338
339 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
340
341 double eI = polarizationEnergy(mI);
342
343 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
344
345 eI += 0.5 * multipoleEnergy(mI);
346 yield c * 0.5 * (eK + eI);
347 }
348 case QUADRUPOLE -> {
349
350 quadrupoleIPotentialAtK(mI, 1);
351
352 double eK = polarizationEnergy(mK);
353
354 quadrupoleKPotentialAtI(mK, 1);
355
356 double eI = polarizationEnergy(mI);
357 yield c * 0.5 * (eK + eI);
358 }
359 };
360 }
361
362
363
364
365
366
367
368
369 public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
370 return switch (multipoleOrder) {
371 default -> {
372
373 chargeIPotentialAtK(mI, 1);
374
375 double eK = polarizationEnergyS(mK);
376
377 chargeKPotentialAtI(mK, 1);
378
379 double eI = polarizationEnergyS(mI);
380 yield c * 0.5 * (eK + eI);
381 }
382 case DIPOLE -> {
383
384 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
385
386 double eK = polarizationEnergyS(mK);
387
388 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
389
390 eK += 0.5 * multipoleEnergy(mK);
391
392 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
393
394 double eI = polarizationEnergyS(mI);
395
396 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
397
398 eI += 0.5 * multipoleEnergy(mI);
399 yield c * 0.5 * (eK + eI);
400 }
401 case QUADRUPOLE -> {
402
403 quadrupoleIPotentialAtK(mI, 1);
404
405 double eK = polarizationEnergyS(mK);
406
407 quadrupoleKPotentialAtI(mK, 1);
408
409 double eI = polarizationEnergyS(mI);
410 yield c * 0.5 * (eK + eI);
411 }
412 };
413 }
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430 @Override
431 public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
432 double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
433 double[] Tk) {
434 return switch (multipoleOrder) {
435 default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
436 case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
437 case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
438 };
439 }
440
441
442
443
444
445
446
447
448
449 public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
450 PolarizableMultipole mK, double[] Gi) {
451
452
453 chargeIPotentialAtK(mI, 2);
454
455 double eK = polarizationEnergy(mK);
456
457 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
458 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
459 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
460
461
462 chargeKPotentialAtI(mK, 2);
463
464 double eI = polarizationEnergy(mI);
465
466 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
467 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
468 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
469
470 double scale = c * 0.5;
471 Gi[0] *= scale;
472 Gi[1] *= scale;
473 Gi[2] *= scale;
474
475
476 return scale * (eI + eK);
477 }
478
479
480
481
482
483
484
485
486
487
488
489
490 public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
491 double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
492
493
494 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
495
496 double eK = polarizationEnergy(mK);
497
498 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
499 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
500 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
501
502 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
503 dipoleTorque(mI, Ti);
504
505
506 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
507
508 eK += 0.5 * multipoleEnergy(mK);
509
510 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 3);
511 double[] G = new double[3];
512 multipoleGradient(mK, G);
513 Gi[0] -= G[0];
514 Gi[1] -= G[1];
515 Gi[2] -= G[2];
516 multipoleTorque(mK, Tk);
517
518
519 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
520
521 double eI = polarizationEnergy(mI);
522
523 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
524 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
525 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
526
527 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
528 dipoleTorque(mK, Tk);
529
530
531 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
532
533 eI += 0.5 * multipoleEnergy(mI);
534
535 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 3);
536 G = new double[3];
537 multipoleGradient(mI, G);
538 Gi[0] += G[0];
539 Gi[1] += G[1];
540 Gi[2] += G[2];
541 multipoleTorque(mI, Ti);
542
543
544
545 if (mutualMask != 0.0) {
546
547 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
548 Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
549 Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
550 Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
551
552
553 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
554 Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
555 Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
556 Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
557 }
558
559
560 double scale = c * 0.5;
561 double energy = scale * (eI + eK);
562 Gi[0] *= scale;
563 Gi[1] *= scale;
564 Gi[2] *= scale;
565 Ti[0] *= scale;
566 Ti[1] *= scale;
567 Ti[2] *= scale;
568 Tk[0] *= scale;
569 Tk[1] *= scale;
570 Tk[2] *= scale;
571
572 return energy;
573 }
574
575
576
577
578
579
580
581
582
583
584
585 public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI,
586 PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {
587
588
589 quadrupoleIPotentialAtK(mI, 2);
590
591 double eK = polarizationEnergy(mK);
592
593 Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
594 Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
595 Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
596
597
598 quadrupoleKPotentialAtI(mK, 2);
599
600 double eI = polarizationEnergy(mI);
601
602 Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
603 Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
604 Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
605
606 double scale = c * 0.5;
607 Gi[0] *= scale;
608 Gi[1] *= scale;
609 Gi[2] *= scale;
610
611
612 dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
613 quadrupoleTorque(mK, Tk);
614
615
616 dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
617 quadrupoleTorque(mI, Ti);
618
619
620 return scale * (eI + eK);
621 }
622
623
624
625
626
627
628
629
630 public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
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 == GKMultipoleOrder.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 }