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.DoubleMath.length;
41 import static ffx.numerics.multipole.MultipoleTensorTest.Qi;
42 import static ffx.numerics.multipole.MultipoleTensorTest.Qk;
43 import static ffx.numerics.multipole.MultipoleTensorTest.Ui;
44 import static ffx.numerics.multipole.MultipoleTensorTest.UiEwald;
45 import static ffx.numerics.multipole.MultipoleTensorTest.Uk;
46 import static ffx.numerics.multipole.MultipoleTensorTest.UkEwald;
47 import static ffx.numerics.multipole.MultipoleTensorTest.permTorqueI;
48 import static ffx.numerics.multipole.MultipoleTensorTest.permTorqueIEwald;
49 import static ffx.numerics.multipole.MultipoleTensorTest.permTorqueK;
50 import static ffx.numerics.multipole.MultipoleTensorTest.permTorqueKEwald;
51 import static ffx.numerics.multipole.MultipoleTensorTest.permanentEnergy;
52 import static ffx.numerics.multipole.MultipoleTensorTest.permanentEnergyEwald;
53 import static ffx.numerics.multipole.MultipoleTensorTest.polarGradICoulomb;
54 import static ffx.numerics.multipole.MultipoleTensorTest.polarGradIEwald;
55 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueICoulomb;
56 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueIEwald;
57 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueKCoulomb;
58 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueKEwald;
59 import static ffx.numerics.multipole.MultipoleTensorTest.polarizationEnergyCoulomb;
60 import static ffx.numerics.multipole.MultipoleTensorTest.polarizationEnergyEwald;
61 import static ffx.numerics.multipole.MultipoleTensorTest.scaleMutual;
62 import static ffx.numerics.multipole.MultipoleTensorTest.polarizationEnergyThole;
63 import static ffx.numerics.multipole.MultipoleTensorTest.polarGradIThole;
64 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueIThole;
65 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueKThole;
66 import static ffx.numerics.multipole.MultipoleTensorTest.xyz;
67 import static java.lang.String.format;
68 import static java.util.Arrays.fill;
69 import static org.junit.Assert.assertEquals;
70
71 import ffx.numerics.multipole.MultipoleTensor.OPERATOR;
72 import java.util.Arrays;
73 import java.util.Collection;
74 import java.util.logging.Logger;
75
76 import ffx.utilities.FFXTest;
77 import org.junit.Test;
78 import org.junit.runner.RunWith;
79 import org.junit.runners.Parameterized;
80 import org.junit.runners.Parameterized.Parameters;
81
82
83
84
85
86
87
88 @RunWith(Parameterized.class)
89 public class QIMultipoleTensorTest extends FFXTest {
90
91
92 private static final Logger logger = Logger.getLogger(QIMultipoleTensorTest.class.getName());
93
94 private final double tolerance = 1.0e-13;
95 private final double fdTolerance = 1.0e-6;
96 private final double[] r = new double[3];
97 private final double[] tensor;
98 private final double[] noStorageTensor;
99 private final double[] fastTensor;
100 private final int order;
101 private final int tensorCount;
102 private final String info;
103 private final OPERATOR operator;
104 private final double beta = 0.545;
105 private final double thole = 0.39;
106 private final double AiAk = 1.061104559485911;
107
108 public QIMultipoleTensorTest(
109 String info,
110 int order,
111 OPERATOR operator) {
112 this.info = info;
113 this.order = order;
114 r[0] = xyz[0];
115 r[1] = xyz[1];
116 r[2] = xyz[2];
117 this.tensorCount = MultipoleTensor.tensorCount(order);
118 tensor = new double[tensorCount];
119 noStorageTensor = new double[tensorCount];
120 fastTensor = new double[tensorCount];
121 this.operator = operator;
122 }
123
124 @Parameters
125 public static Collection<Object[]> data() {
126 return Arrays.asList(
127 new Object[][] {
128 {
129 "Order 5 Coulomb",
130 5,
131 OPERATOR.COULOMB
132 },
133 {
134 "Order 5 Screened Coulomb",
135 5,
136 OPERATOR.SCREENED_COULOMB
137 },
138 {
139 "Order 4 Thole Field",
140 4,
141 OPERATOR.THOLE_FIELD
142 }
143 });
144 }
145
146 @Test
147 public void codeGenerationTest() {
148
149 if (operator != OPERATOR.COULOMB) {
150 return;
151 }
152
153 double[] Qi = {0.11, 0.22, 0.33, 0.44, 0.11, 0.22, -0.33, 0.12, 0.13, 0.14};
154 double[] r = {2.11, 2.12, 2.13};
155
156 double[] tensor = new double[tensorCount];
157 MultipoleTensor multipoleTensor = new CoulombTensorQI(order);
158 logger.info(format(" Writing QI Order %d tensor recursion code:", order));
159
160 r[2] = length(r);
161 r[0] = 0.0;
162 r[1] = 0.0;
163 String code = multipoleTensor.codeTensorRecursion(r, tensor);
164 logger.info(format("\n%s", code));
165
166 PolarizableMultipole polarizableMultipole = new PolarizableMultipole(Qi, Ui, Ui);
167 QIFrame qiFrame = new QIFrame(r);
168 qiFrame.rotatePolarizableMultipole(polarizableMultipole);
169 StringBuilder sb = new StringBuilder();
170 logger.info(" Writing QI potential code due to multipole I:");
171 multipoleTensor.codePotentialMultipoleI(polarizableMultipole, tensor, 0, 0, 0, sb);
172 logger.info("\n" + sb);
173
174 sb = new StringBuilder();
175 logger.info(" Writing QI potential code due to multipole K:");
176 multipoleTensor.codePotentialMultipoleK(polarizableMultipole, tensor, 0, 0, 0, sb);
177 logger.info("\n" + sb);
178 }
179
180 @Test
181 public void permanentMultipoleEnergyAndGradTest() {
182
183 if (operator == OPERATOR.THOLE_FIELD) {
184 return;
185 }
186
187 MultipoleTensor multipoleTensor = new CoulombTensorQI(order);
188 if (operator == OPERATOR.SCREENED_COULOMB) {
189 multipoleTensor = new EwaldTensorQI(order, beta);
190 }
191
192 double delta = 1.0e-5;
193 double delta2 = 2.0 * 1.0e-5;
194 double[] Gi = new double[3];
195 double[] Gk = new double[3];
196 double[] Ti = new double[3];
197 double[] Tk = new double[3];
198
199 PolarizableMultipole mI = new PolarizableMultipole(Qi, Ui, Ui);
200 PolarizableMultipole mK = new PolarizableMultipole(Qk, Uk, Uk);
201 QIFrame qiFrame = new QIFrame(r[0], r[1], r[2]);
202 qiFrame.rotatePolarizableMultipole(mI);
203 qiFrame.rotatePolarizableMultipole(mK);
204
205 multipoleTensor.generateTensor(r);
206 double e = multipoleTensor.multipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
207
208 qiFrame.toGlobal(Gk);
209 qiFrame.toGlobal(Ti);
210 qiFrame.toGlobal(Tk);
211
212 if (operator == OPERATOR.COULOMB) {
213 assertEquals(info + " QI Permanent Energy", permanentEnergy, e, tolerance);
214 assertEquals(info + " QI Multipole Torque I", permTorqueI[1], Ti[1], tolerance);
215 assertEquals(info + " QI Multipole Torque K", permTorqueK[1], Tk[1], tolerance);
216 } else if (operator == OPERATOR.SCREENED_COULOMB) {
217 assertEquals(info + " QI Ewald Permanent Energy", permanentEnergyEwald, e, tolerance);
218 assertEquals(info + " QI Ewald Multipole Torque I", permTorqueIEwald[1], Ti[1], tolerance);
219 assertEquals(info + " QI Ewald Multipole Torque K", permTorqueKEwald[1], Tk[1], tolerance);
220 }
221
222
223 double aX = Gk[0];
224 double aY = Gk[1];
225 double aZ = Gk[2];
226
227 r[0] += delta;
228 mI.set(Qi, Ui, Ui);
229 mK.set(Qk, Uk, Uk);
230 qiFrame.setAndRotate(r, mI, mK);
231 multipoleTensor.generateTensor(r);
232 double posX = multipoleTensor.multipoleEnergy(mI, mK);
233
234 r[0] -= delta2;
235 mI.set(Qi, Ui, Ui);
236 mK.set(Qk, Uk, Uk);
237 qiFrame.setAndRotate(r, mI, mK);
238 multipoleTensor.generateTensor(r);
239 double negX = multipoleTensor.multipoleEnergy(mI, mK);
240 r[0] += delta;
241
242 r[1] += delta;
243 mI.set(Qi, Ui, Ui);
244 mK.set(Qk, Uk, Uk);
245 qiFrame.setAndRotate(r, mI, mK);
246 multipoleTensor.generateTensor(r);
247 double posY = multipoleTensor.multipoleEnergy(mI, mK);
248 r[1] -= delta2;
249 mI.set(Qi, Ui, Ui);
250 mK.set(Qk, Uk, Uk);
251 qiFrame.setAndRotate(r, mI, mK);
252 multipoleTensor.generateTensor(r);
253 double negY = multipoleTensor.multipoleEnergy(mI, mK);
254 r[1] += delta;
255
256 r[2] += delta;
257 mI.set(Qi, Ui, Ui);
258 mK.set(Qk, Uk, Uk);
259 qiFrame.setAndRotate(r, mI, mK);
260 multipoleTensor.generateTensor(r);
261 double posZ = multipoleTensor.multipoleEnergy(mI, mK);
262 r[2] -= delta2;
263 mI.set(Qi, Ui, Ui);
264 mK.set(Qk, Uk, Uk);
265 qiFrame.setAndRotate(r, mI, mK);
266 multipoleTensor.generateTensor(r);
267 double negZ = multipoleTensor.multipoleEnergy(mI, mK);
268 r[2] += delta;
269
270 double expect = aX;
271 double actual = (posX - negX) / delta2;
272 assertEquals(info + " QI Force X", expect, actual, fdTolerance);
273
274 expect = aY;
275 actual = (posY - negY) / delta2;
276 assertEquals(info + " QI Force Y", expect, actual, fdTolerance);
277
278 expect = aZ;
279 actual = (posZ - negZ) / delta2;
280 assertEquals(info + " QI Force Z", expect, actual, fdTolerance);
281 }
282
283 @Test
284 public void polarizationEnergyAndGradTest() {
285
286 MultipoleTensor multipoleTensor = new CoulombTensorQI(order);
287 if (operator == OPERATOR.THOLE_FIELD) {
288 multipoleTensor = new TholeTensorQI(order, thole, AiAk);
289 } else if (operator == OPERATOR.SCREENED_COULOMB) {
290 multipoleTensor = new EwaldTensorQI(order, beta);
291 }
292
293 double delta = 1.0e-5;
294 double delta2 = 2.0 * 1.0e-5;
295 double[] Gi = new double[3];
296 double[] Ti = new double[3];
297 double[] Tk = new double[3];
298
299
300
301
302 PolarizableMultipole mI = new PolarizableMultipole(Qi, Ui, Ui);
303 PolarizableMultipole mK = new PolarizableMultipole(Qk, Uk, Uk);
304 if (operator == OPERATOR.SCREENED_COULOMB) {
305 mI.setInducedDipole(UiEwald, UiEwald);
306 mK.setInducedDipole(UkEwald, UkEwald);
307 }
308
309 QIFrame qiFrame = new QIFrame(r[0], r[1], r[2]);
310 qiFrame.rotatePolarizableMultipole(mI);
311 qiFrame.rotatePolarizableMultipole(mK);
312
313 multipoleTensor.generateTensor(r);
314 double e = multipoleTensor.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, scaleMutual, Gi, Ti,
315 Tk);
316
317 qiFrame.toGlobal(Gi);
318 qiFrame.toGlobal(Ti);
319 qiFrame.toGlobal(Tk);
320
321
322 double aX = Gi[0];
323 double aY = Gi[1];
324 double aZ = Gi[2];
325
326 double energy = polarizationEnergyCoulomb;
327 double[] gradI = polarGradICoulomb;
328 double[] torqueI = polarTorqueICoulomb;
329 double[] torqueK = polarTorqueKCoulomb;
330 if (operator == OPERATOR.THOLE_FIELD) {
331 energy = polarizationEnergyThole;
332 gradI = polarGradIThole;
333 torqueI = polarTorqueIThole;
334 torqueK = polarTorqueKThole;
335 } else if (operator == OPERATOR.SCREENED_COULOMB) {
336 energy = polarizationEnergyEwald;
337 gradI = polarGradIEwald;
338 torqueI = polarTorqueIEwald;
339 torqueK = polarTorqueKEwald;
340 }
341 assertEquals(info + " Polarization Energy", energy, e, tolerance);
342 assertEquals(info + " Polarization GradX", gradI[0], aX, tolerance);
343 assertEquals(info + " Polarization GradY", gradI[1], aY, tolerance);
344 assertEquals(info + " Polarization GradZ", gradI[2], aZ, tolerance);
345 assertEquals(info + " Polarization TorqueX I", torqueI[0], Ti[0], tolerance);
346 assertEquals(info + " Polarization TorqueY I", torqueI[1], Ti[1], tolerance);
347 assertEquals(info + " Polarization TorqueZ I", torqueI[2], Ti[2], tolerance);
348 assertEquals(info + " Polarization TorqueX K", torqueK[0], Tk[0], tolerance);
349 assertEquals(info + " Polarization TorqueY K", torqueK[1], Tk[1], tolerance);
350 assertEquals(info + " Polarization TorqueZ K", torqueK[2], Tk[2], tolerance);
351
352 if (scaleMutual == 0.0) {
353
354 r[0] += delta;
355 mI.set(Qi, Ui, Ui);
356 mK.set(Qk, Uk, Uk);
357 qiFrame.setAndRotate(r, mI, mK);
358 multipoleTensor.generateTensor(r);
359 double posX = multipoleTensor.polarizationEnergy(mI, mK, 1.0);
360
361 r[0] -= delta2;
362 mI.set(Qi, Ui, Ui);
363 mK.set(Qk, Uk, Uk);
364 qiFrame.setAndRotate(r, mI, mK);
365 multipoleTensor.generateTensor(r);
366 double negX = multipoleTensor.polarizationEnergy(mI, mK, 1.0);
367 r[0] += delta;
368
369 r[1] += delta;
370 mI.set(Qi, Ui, Ui);
371 mK.set(Qk, Uk, Uk);
372 qiFrame.setAndRotate(r, mI, mK);
373 multipoleTensor.generateTensor(r);
374 double posY = multipoleTensor.polarizationEnergy(mI, mK, 1.0);
375 r[1] -= delta2;
376 mI.set(Qi, Ui, Ui);
377 mK.set(Qk, Uk, Uk);
378 qiFrame.setAndRotate(r, mI, mK);
379 multipoleTensor.generateTensor(r);
380 double negY = multipoleTensor.polarizationEnergy(mI, mK, 1.0);
381 r[1] += delta;
382
383 r[2] += delta;
384 mI.set(Qi, Ui, Ui);
385 mK.set(Qk, Uk, Uk);
386 qiFrame.setAndRotate(r, mI, mK);
387 multipoleTensor.generateTensor(r);
388 double posZ = multipoleTensor.polarizationEnergy(mI, mK, 1.0);
389 r[2] -= delta2;
390 mI.set(Qi, Ui, Ui);
391 mK.set(Qk, Uk, Uk);
392 qiFrame.setAndRotate(r, mI, mK);
393 multipoleTensor.generateTensor(r);
394 double negZ = multipoleTensor.polarizationEnergy(mI, mK, 1.0);
395 r[2] += delta;
396
397
398
399 double expect = -aX;
400 double actual = 2.0 * (posX - negX) / delta2;
401 assertEquals(info + " QI Polarization Force X", expect, actual, fdTolerance);
402
403 expect = -aY;
404 actual = 2.0 * (posY - negY) / delta2;
405 assertEquals(info + " QI Polarization Force Y", expect, actual, fdTolerance);
406
407 expect = -aZ;
408 actual = 2.0 * (posZ - negZ) / delta2;
409 assertEquals(info + " QI Polarization Force Z", expect, actual, fdTolerance);
410 }
411 }
412
413 @Test
414 public void finiteDifferenceTest() {
415
416 MultipoleTensor multipoleTensor = new CoulombTensorQI(order);
417 if (operator == OPERATOR.THOLE_FIELD) {
418 multipoleTensor = new TholeTensorQI(order, thole, AiAk);
419 } else if (operator == OPERATOR.SCREENED_COULOMB) {
420 multipoleTensor = new EwaldTensorQI(order, beta);
421 }
422
423 r[0] = 0.0;
424 r[1] = 0.0;
425 r[2] = 2.1;
426
427 multipoleTensor.noStorageRecursion(r, tensor);
428 double[] tensorsPz = new double[tensorCount];
429 double[] tensorsNz = new double[tensorCount];
430 double delta = 1.0e-5;
431 double delta2 = delta * 2;
432
433 r[2] += delta;
434 multipoleTensor.noStorageRecursion(r, tensorsPz);
435 r[2] -= delta2;
436 multipoleTensor.noStorageRecursion(r, tensorsNz);
437 r[2] += delta;
438
439 tensorFiniteDifference(multipoleTensor, delta2, tensorsPz, tensorsNz);
440
441 multipoleTensor.recursion(r, tensor);
442 tensorFiniteDifference(multipoleTensor, delta2, tensorsPz, tensorsNz);
443
444 multipoleTensor.generateTensor();
445 multipoleTensor.getTensor(tensor);
446 tensorFiniteDifference(multipoleTensor, delta2, tensorsPz, tensorsNz);
447 }
448
449 @Test
450 public void multipoleTensorTest() {
451
452 MultipoleTensor multipoleTensor = new CoulombTensorQI(order);
453 if (operator == OPERATOR.THOLE_FIELD) {
454 multipoleTensor = new TholeTensorQI(order, thole, AiAk);
455 } else if (operator == OPERATOR.SCREENED_COULOMB) {
456 multipoleTensor = new EwaldTensorQI(order, beta);
457 }
458
459
460 multipoleTensor.noStorageRecursion(r, noStorageTensor);
461 multipoleTensor.recursion(r, tensor);
462 multipoleTensor.generateTensor();
463 multipoleTensor.getTensor(fastTensor);
464
465 for (int i = 0; i < tensorCount; i++) {
466 double expect = noStorageTensor[i];
467 double actual = tensor[i];
468 assertEquals(info + " @ " + i, expect, actual, tolerance);
469 if (order == 4 || order == 5) {
470 expect = noStorageTensor[i];
471 actual = fastTensor[i];
472 assertEquals(info + " @ " + i, expect, actual, tolerance);
473 }
474 }
475
476
477
478 r[2] = length(r);
479 r[0] = 0.0;
480 r[1] = 0.0;
481 fill(noStorageTensor, 0.0);
482 fill(tensor, 0.0);
483 fill(fastTensor, 0.0);
484 multipoleTensor.noStorageRecursion(r, noStorageTensor);
485 multipoleTensor.recursion(r, tensor);
486 multipoleTensor.setTensor(fastTensor);
487 multipoleTensor.generateTensor();
488 multipoleTensor.getTensor(fastTensor);
489
490 for (int i = 0; i < tensorCount; i++) {
491 double expect = noStorageTensor[i];
492 double actual = tensor[i];
493 assertEquals(info + " @ " + i, expect, actual, tolerance);
494 if (order == 4 || order == 5) {
495 expect = noStorageTensor[i];
496 actual = fastTensor[i];
497 assertEquals(info + " @ " + i, expect, actual, tolerance);
498 }
499 }
500 }
501
502 private void tensorFiniteDifference(
503 MultipoleTensor multipoleTensor, double delta2, double[] tensorsPz, double[] tensorsNz) {
504
505 int start = 0;
506
507
508 if (operator == OPERATOR.THOLE_FIELD) {
509 start = 1;
510 }
511
512
513 for (int l = start; l < order; l++) {
514
515 double expect = tensor[multipoleTensor.ti(l, 0, 1)];
516 double actual =
517 (tensorsPz[multipoleTensor.ti(l, 0, 0)] - tensorsNz[multipoleTensor.ti(l, 0, 0)]) / delta2;
518 assertEquals(info + " (d/dx): " + l, expect, actual, fdTolerance);
519 }
520 for (int l = 0; l < order; l++) {
521 for (int m = 1; m < order - l; m++) {
522
523 double expect = tensor[multipoleTensor.ti(l, m, 1)];
524 double actual =
525 (tensorsPz[multipoleTensor.ti(l, m, 0)] - tensorsNz[multipoleTensor.ti(l, m, 0)])
526 / delta2;
527 assertEquals(info + " (d/dx): " + l + " (d/dy): " + m, expect, actual, fdTolerance);
528 }
529 }
530
531 for (int l = 0; l < order; l++) {
532 for (int m = 0; m < order - l; m++) {
533 for (int n = 1; n < order - l - m; n++) {
534
535 double expect = tensor[multipoleTensor.ti(l, m, n + 1)];
536 double actual =
537 (tensorsPz[multipoleTensor.ti(l, m, n)] - tensorsNz[multipoleTensor.ti(l, m, n)])
538 / delta2;
539 assertEquals(info + " (d/dx): " + l + " (d/dy): " + m + " (d/dz):" + n, expect, actual,
540 fdTolerance);
541 }
542 }
543 }
544 }
545 }