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.utilities.FFXTest;
41 import org.junit.Test;
42
43 import static ffx.numerics.math.DoubleMath.length;
44 import static ffx.numerics.math.DoubleMath.length2;
45 import static ffx.numerics.multipole.GKMultipoleOrder.DIPOLE;
46 import static ffx.numerics.multipole.GKMultipoleOrder.MONOPOLE;
47 import static ffx.numerics.multipole.GKMultipoleOrder.QUADRUPOLE;
48 import static ffx.numerics.multipole.GKTensorGlobalTest.bornI;
49 import static ffx.numerics.multipole.GKTensorGlobalTest.bornK;
50 import static ffx.numerics.multipole.GKTensorGlobalTest.multI;
51 import static ffx.numerics.multipole.GKTensorGlobalTest.multK;
52 import static ffx.numerics.multipole.GKTensorGlobalTest.permGradBorn;
53 import static ffx.numerics.multipole.GKTensorGlobalTest.permGradI;
54 import static ffx.numerics.multipole.GKTensorGlobalTest.permTorqueI;
55 import static ffx.numerics.multipole.GKTensorGlobalTest.permTorqueK;
56 import static ffx.numerics.multipole.GKTensorGlobalTest.polGradBorn;
57 import static ffx.numerics.multipole.GKTensorGlobalTest.polGradBornDirect;
58 import static ffx.numerics.multipole.GKTensorGlobalTest.polGradI;
59 import static ffx.numerics.multipole.GKTensorGlobalTest.polGradIDirect;
60 import static ffx.numerics.multipole.GKTensorGlobalTest.polTorqueI;
61 import static ffx.numerics.multipole.GKTensorGlobalTest.polTorqueIDirect;
62 import static ffx.numerics.multipole.GKTensorGlobalTest.polTorqueK;
63 import static ffx.numerics.multipole.GKTensorGlobalTest.polTorqueKDirect;
64 import static ffx.numerics.multipole.GKTensorGlobalTest.rWater;
65 import static ffx.numerics.multipole.GKTensorGlobalTest.uI;
66 import static ffx.numerics.multipole.GKTensorGlobalTest.uIDirect;
67 import static ffx.numerics.multipole.GKTensorGlobalTest.uK;
68 import static ffx.numerics.multipole.GKTensorGlobalTest.uKDirect;
69 import static ffx.numerics.multipole.GKTensorGlobalTest.watPermEnergy;
70 import static ffx.numerics.multipole.GKTensorGlobalTest.watPolDirect;
71 import static ffx.numerics.multipole.GKTensorGlobalTest.watPolEnergy;
72 import static ffx.numerics.multipole.GKTensorMode.BORN;
73 import static ffx.numerics.multipole.GKTensorMode.POTENTIAL;
74 import static org.junit.Assert.assertEquals;
75
76
77
78
79
80
81
82
83 public class GKTensorQITest extends FFXTest {
84
85 private final double tolerance = 1.0e-9;
86 private final double fdTolerance = 1.0e-6;
87 private static final double gc = 2.455;
88 private static final double Eh = 1.0;
89 private static final double Es = 78.3;
90 private static final double Ai = 2.0;
91 private static final double Aj = 2.0;
92 private static final double[] r = {0.7, 0.8, 0.9};
93
94 @Test
95 public void permanentEnergyTest() {
96 PolarizableMultipole mI = new PolarizableMultipole(multI, uI, uI);
97 PolarizableMultipole mK = new PolarizableMultipole(multK, uK, uK);
98 QIFrame qiFrame = new QIFrame(rWater);
99 qiFrame.rotatePolarizableMultipole(mI);
100 qiFrame.rotatePolarizableMultipole(mK);
101
102 GKEnergyQI gkEnergyQI = new GKEnergyQI(Eh, Es, gc, false);
103 gkEnergyQI.initPotential(rWater, length2(rWater), bornI, bornK);
104 var e = gkEnergyQI.multipoleEnergy(mI, mK);
105 assertEquals("GK Permanent Energy", watPermEnergy, e, tolerance);
106 }
107
108 @Test
109 public void permanentEnergyAndGradientTest() {
110 PolarizableMultipole mI = new PolarizableMultipole(multI, uI, uI);
111 PolarizableMultipole mK = new PolarizableMultipole(multK, uK, uK);
112 QIFrame qiFrame = new QIFrame(rWater);
113 qiFrame.rotatePolarizableMultipole(mI);
114 qiFrame.rotatePolarizableMultipole(mK);
115
116 double[] gradI = new double[3];
117 double[] torqueI = new double[3];
118 double[] torqueK = new double[3];
119
120 double r2 = length2(rWater);
121 GKEnergyQI gkEnergyQI = new GKEnergyQI(Eh, Es, gc, true);
122 gkEnergyQI.initPotential(rWater, r2, bornI, bornK);
123 var e = gkEnergyQI.multipoleEnergyAndGradient(mI, mK, gradI, torqueI, torqueK);
124
125 gkEnergyQI.initBorn(rWater, r2, bornI, bornK);
126 var db = gkEnergyQI.multipoleEnergyBornGrad(mI, mK);
127
128 qiFrame.toGlobal(gradI);
129 qiFrame.toGlobal(torqueI);
130 qiFrame.toGlobal(torqueK);
131
132 assertEquals("GK Permanent Energy", watPermEnergy, e, tolerance);
133 assertEquals("GK Permanent Grad X", permGradI[0], gradI[0], tolerance);
134 assertEquals("GK Permanent Grad Y", permGradI[1], gradI[1], tolerance);
135 assertEquals("GK Permanent Grad Z", permGradI[2], gradI[2], tolerance);
136 assertEquals("GK Permanent Torque I X", permTorqueI[0], torqueI[0], tolerance);
137 assertEquals("GK Permanent Torque I Y", permTorqueI[1], torqueI[1], tolerance);
138 assertEquals("GK Permanent Torque I Z", permTorqueI[2], torqueI[2], tolerance);
139 assertEquals("GK Permanent Torque K X", permTorqueK[0], torqueK[0], tolerance);
140 assertEquals("GK Permanent Torque K Y", permTorqueK[1], torqueK[1], tolerance);
141 assertEquals("GK Permanent Torque K Z", permTorqueK[2], torqueK[2], tolerance);
142 assertEquals("GK Born Grad I", permGradBorn, db * bornI, tolerance);
143 }
144
145 @Test
146 public void polarizationEnergyTest() {
147 PolarizableMultipole mI = new PolarizableMultipole(multI, uI, uI);
148 PolarizableMultipole mK = new PolarizableMultipole(multK, uK, uK);
149
150 QIFrame qiFrame = new QIFrame(rWater);
151 qiFrame.rotatePolarizableMultipole(mI);
152 qiFrame.rotatePolarizableMultipole(mK);
153
154 GKEnergyQI gkEnergyQI = new GKEnergyQI(Eh, Es, gc, false);
155 gkEnergyQI.initPotential(rWater, length2(rWater), bornI, bornK);
156 var e = gkEnergyQI.polarizationEnergy(mI, mK);
157
158 assertEquals("GK Polarization Energy", watPolEnergy, e, tolerance);
159 }
160
161 @Test
162 public void polarizationEnergyDirectTest() {
163 PolarizableMultipole mI = new PolarizableMultipole(multI, uIDirect, uIDirect);
164 PolarizableMultipole mK = new PolarizableMultipole(multK, uKDirect, uKDirect);
165
166 QIFrame qiFrame = new QIFrame(rWater);
167 qiFrame.rotatePolarizableMultipole(mI);
168 qiFrame.rotatePolarizableMultipole(mK);
169
170 GKEnergyQI gkEnergyQI = new GKEnergyQI(Eh, Es, gc, false);
171 gkEnergyQI.initPotential(rWater, length2(rWater), bornI, bornK);
172 var e = gkEnergyQI.polarizationEnergy(mI, mK);
173
174 assertEquals("GK Direct Polarization Energy", watPolDirect, e, tolerance);
175 }
176
177 @Test
178 public void polarizationEnergyAndGradientTest() {
179 PolarizableMultipole mI = new PolarizableMultipole(multI, uI, uI);
180 PolarizableMultipole mK = new PolarizableMultipole(multK, uK, uK);
181
182 QIFrame qiFrame = new QIFrame(rWater);
183 qiFrame.rotatePolarizableMultipole(mI);
184 qiFrame.rotatePolarizableMultipole(mK);
185
186 double[] gradI = new double[3];
187 double[] torqueI = new double[3];
188 double[] torqueK = new double[3];
189
190 double r2 = length2(rWater);
191 GKEnergyQI gkEnergyQI = new GKEnergyQI(Eh, Es, gc, true);
192 gkEnergyQI.initPotential(rWater, r2, bornI, bornK);
193 var e = gkEnergyQI.polarizationEnergyAndGradient(mI, mK, 1.0, gradI, torqueI, torqueK);
194
195 gkEnergyQI.initBorn(rWater, r2, bornI, bornK);
196 var db = gkEnergyQI.polarizationEnergyBornGrad(mI, mK, true);
197
198 qiFrame.toGlobal(gradI);
199 qiFrame.toGlobal(torqueI);
200 qiFrame.toGlobal(torqueK);
201
202 assertEquals("GK Polarization Energy", watPolEnergy, e, tolerance);
203 assertEquals("GK Polarization Grad X", polGradI[0], gradI[0], tolerance);
204 assertEquals("GK Polarization Grad Y", polGradI[1], gradI[1], tolerance);
205 assertEquals("GK Polarization Grad Z", polGradI[2], gradI[2], tolerance);
206 assertEquals("GK Polarization Torque I X", polTorqueI[0], torqueI[0], tolerance);
207 assertEquals("GK Polarization Torque I Y", polTorqueI[1], torqueI[1], tolerance);
208 assertEquals("GK Polarization Torque I Z", polTorqueI[2], torqueI[2], tolerance);
209 assertEquals("GK Polarization Torque K X", polTorqueK[0], torqueK[0], tolerance);
210 assertEquals("GK Polarization Torque K Y", polTorqueK[1], torqueK[1], tolerance);
211 assertEquals("GK Polarization Torque K Z", polTorqueK[2], torqueK[2], tolerance);
212 assertEquals("GK Born Grad I", polGradBorn, db * bornI, tolerance);
213 }
214
215 @Test
216 public void polarizationEnergyAndGradientDirectTest() {
217 PolarizableMultipole mI = new PolarizableMultipole(multI, uIDirect, uIDirect);
218 PolarizableMultipole mK = new PolarizableMultipole(multK, uKDirect, uKDirect);
219
220 QIFrame qiFrame = new QIFrame(rWater);
221 qiFrame.rotatePolarizableMultipole(mI);
222 qiFrame.rotatePolarizableMultipole(mK);
223
224 double[] gradI = new double[3];
225 double[] torqueI = new double[3];
226 double[] torqueK = new double[3];
227
228 double r2 = length2(rWater);
229 GKEnergyQI gkEnergyQI = new GKEnergyQI(Eh, Es, gc, true);
230 gkEnergyQI.initPotential(rWater, r2, bornI, bornK);
231 var e = gkEnergyQI.polarizationEnergyAndGradient(mI, mK, 0.0, gradI, torqueI, torqueK);
232
233 gkEnergyQI.initBorn(rWater, r2, bornI, bornK);
234 var db = gkEnergyQI.polarizationEnergyBornGrad(mI, mK, false);
235
236 qiFrame.toGlobal(gradI);
237 qiFrame.toGlobal(torqueI);
238 qiFrame.toGlobal(torqueK);
239
240 assertEquals("GK Polarization Energy", watPolDirect, e, tolerance);
241 assertEquals("GK Polarization Grad X", polGradIDirect[0], gradI[0], tolerance);
242 assertEquals("GK Polarization Grad Y", polGradIDirect[1], gradI[1], tolerance);
243 assertEquals("GK Polarization Grad Z", polGradIDirect[2], gradI[2], tolerance);
244 assertEquals("GK Polarization Torque I X", polTorqueIDirect[0], torqueI[0], tolerance);
245 assertEquals("GK Polarization Torque I Y", polTorqueIDirect[1], torqueI[1], tolerance);
246 assertEquals("GK Polarization Torque I Z", polTorqueIDirect[2], torqueI[2], tolerance);
247 assertEquals("GK Polarization Torque K X", polTorqueKDirect[0], torqueK[0], tolerance);
248 assertEquals("GK Polarization Torque K Y", polTorqueKDirect[1], torqueK[1], tolerance);
249 assertEquals("GK Polarization Torque K Z", polTorqueKDirect[2], torqueK[2], tolerance);
250 assertEquals("GK Born Grad I", polGradBornDirect, db * bornI, tolerance);
251 }
252
253 @Test
254 public void tensorAuxiliaryTest() {
255 int order = 6;
256
257 double r2 = length2(r);
258 GKSource gkSource = new GKSource(order, gc);
259 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, Ai, Aj);
260
261 GKTensorQI tensorQI = new GKTensorQI(MONOPOLE, order, gkSource, Eh, Es);
262 tensorQI.setR(r);
263 double[] work = new double[order + 1];
264 tensorQI.source(work);
265
266
267 gkSource.bn(5);
268 double bn0 = gkSource.bn[0];
269 double bn1 = gkSource.bn[1];
270 double bn2 = gkSource.bn[2];
271 double bn3 = gkSource.bn[3];
272 double bn4 = gkSource.bn[4];
273 double mapleBN0 = 0.4914375691;
274 double mapleBN1 = -0.01651130007;
275 double mapleBN2 = -0.01365916860;
276 double mapleBN3 = 0.006248702126;
277 double mapleBN4 = -0.001978716128;
278 assertEquals("bn0", mapleBN0, bn0, tolerance);
279 assertEquals("bn1", mapleBN1, bn1, tolerance);
280 assertEquals("bn2", mapleBN2, bn2, tolerance);
281 assertEquals("bn3", mapleBN3, bn3, tolerance);
282 assertEquals("bn4", mapleBN4, bn4, tolerance);
283
284
285 double c = GKSource.cn(0, Eh, Es);
286 double A00 = c * work[0];
287 double A01 = c * work[1];
288 double A02 = c * work[2];
289 double A03 = c * work[3];
290 final double mapleA00 = -0.4319767286;
291 final double mapleA01 = 0.05505753880;
292 final double mapleA02 = -0.01542067105;
293 final double mapleA03 = 0.005809287830;
294 assertEquals("A00", mapleA00, A00, tolerance);
295 assertEquals("A01", mapleA01, A01, tolerance);
296 assertEquals("A02", mapleA02, A02, tolerance);
297 assertEquals("A03", mapleA03, A03, tolerance);
298
299
300 gkSource.generateSource(BORN, QUADRUPOLE, r2, Ai, Aj);
301 tensorQI.source(work);
302 double B00 = c * work[0];
303 double B01 = c * work[1];
304 double B02 = c * work[2];
305 final double mapleB00 = 0.04064563792;
306 final double mapleB01 = -0.01690706430;
307 final double mapleB02 = 0.008229167114;
308 assertEquals("B00", mapleB00, B00, tolerance);
309 assertEquals("B01", mapleB01, B01, tolerance);
310 assertEquals("B02", mapleB02, B02, tolerance);
311
312
313 work = new double[order + 1];
314 tensorQI = new GKTensorQI(DIPOLE, order, gkSource, Eh, Es);
315 tensorQI.setR(r);
316 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, Ai, Aj);
317 tensorQI.source(work);
318 c = GKSource.cn(1, Eh, Es);
319 double A10 = c * work[1];
320 double A11 = c * work[2];
321 double A12 = c * work[3];
322 double A13 = c * work[4];
323 final double mapleA10 = 0.08218283800;
324 final double mapleA11 = -0.03142380936;
325 final double mapleA12 = 0.01681150454;
326 final double mapleA13 = -0.01106715285;
327 assertEquals("A10", mapleA10, A10, tolerance);
328 assertEquals("A11", mapleA11, A11, tolerance);
329 assertEquals("A12", mapleA12, A12, tolerance);
330 assertEquals("A13", mapleA13, A13, tolerance);
331
332
333 gkSource.generateSource(BORN, QUADRUPOLE, r2, Ai, Aj);
334 tensorQI.source(work);
335 double B10 = c * work[1];
336 double B11 = c * work[2];
337 double B12 = c * work[3];
338 final double mapleB10 = -0.02319829046;
339 final double mapleB11 = 0.01556309100;
340 final double mapleB12 = -0.01202628190;
341 assertEquals("B10", mapleB10, B10, tolerance);
342 assertEquals("B11", mapleB11, B11, tolerance);
343 assertEquals("B12", mapleB12, B12, tolerance);
344
345
346 work = new double[order + 1];
347 tensorQI = new GKTensorQI(QUADRUPOLE, order, gkSource, Eh, Es);
348 tensorQI.setR(r);
349 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, Ai, Aj);
350 tensorQI.source(work);
351 c = GKSource.cn(2, Eh, Es);
352 double A20 = c * work[2];
353 double A21 = c * work[3];
354 double A22 = c * work[4];
355 double A23 = c * work[5];
356 final double mapleA20 = -0.04710532877;
357 final double mapleA21 = 0.03001901831;
358 final double mapleA22 = -0.02371209206;
359 final double mapleA23 = 0.02187861148;
360 assertEquals("A20", mapleA20, A20, tolerance);
361 assertEquals("A21", mapleA21, A21, tolerance);
362 assertEquals("A22", mapleA22, A22, tolerance);
363 assertEquals("A23", mapleA23, A23, tolerance);
364
365
366 gkSource.generateSource(BORN, QUADRUPOLE, r2, Ai, Aj);
367 tensorQI.source(work);
368 double B20 = c * work[2];
369 double B21 = c * work[3];
370 double B22 = c * work[4];
371 final double mapleB20 = 0.02216121853;
372 final double mapleB21 = -0.02051645869;
373 final double mapleB22 = 0.02137054028;
374 assertEquals("B20", mapleB20, B20, tolerance);
375 assertEquals("B21", mapleB21, B21, tolerance);
376 assertEquals("B22", mapleB22, B22, tolerance);
377 }
378
379 @Test
380 public void chargeTensorTest() {
381
382 double r2 = length2(r);
383 GKSource gkSource = new GKSource(3, gc);
384 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
385
386
387 int order = 3;
388 GKTensorQI gkMonopoleTensor = new GKTensorQI(MONOPOLE, order, gkSource, Eh, Es);
389 gkMonopoleTensor.setR(r);
390 double[] work = new double[order + 1];
391 gkMonopoleTensor.source(work);
392 double A00 = work[0];
393 double A01 = work[1];
394 double A02 = work[2];
395 double A03 = work[3];
396 gkMonopoleTensor.generateTensor();
397
398 double z = length(r);
399 assertEquals("R000", A00, gkMonopoleTensor.R000, tolerance);
400 assertEquals("R001", z * A01, gkMonopoleTensor.R001, tolerance);
401 assertEquals("R002", z * z * A02 + A01, gkMonopoleTensor.R002, tolerance);
402 assertEquals("R003", z * z * z * A03 + 3.0 * z * A02, gkMonopoleTensor.R003, tolerance);
403
404
405 gkSource.generateSource(BORN, QUADRUPOLE, r2, bornI, bornK);
406 gkMonopoleTensor.source(work);
407 double B00 = work[0];
408 double B01 = work[1];
409 double B02 = work[2];
410 gkMonopoleTensor.generateTensor();
411 assertEquals("B000", B00, gkMonopoleTensor.R000, tolerance);
412 assertEquals("B001", z * B01, gkMonopoleTensor.R001, tolerance);
413 assertEquals("B002", z * z * B02 + B01, gkMonopoleTensor.R002, tolerance);
414 }
415
416 @Test
417 public void dipoleTensorTest() {
418 double r2 = length2(r);
419 GKSource gkSource = new GKSource(4, gc);
420 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
421
422
423 int order = 4;
424 GKTensorQI gkDipoleTensor = new GKTensorQI(DIPOLE, order, gkSource, Eh, Es);
425 gkDipoleTensor.setR(r);
426 double[] work = new double[order + 1];
427 gkDipoleTensor.source(work);
428 double A10 = work[1];
429 double A11 = work[2];
430 double A12 = work[3];
431 double A13 = work[4];
432 gkDipoleTensor.generateTensor();
433
434
435 assertEquals("R000", 0.0, gkDipoleTensor.R000, tolerance);
436
437 double z = length(r);
438 assertEquals("R001", z * A10, gkDipoleTensor.R001, tolerance);
439 assertEquals("R002", z * z * A11 + A10, gkDipoleTensor.R002, tolerance);
440 assertEquals("R003", z * z * z * A12 + 3.0 * z * A11, gkDipoleTensor.R003, tolerance);
441 assertEquals("R004", z * z * z * z * A13 + 6.0 * z * z * A12 + 3.0 * A11, gkDipoleTensor.R004,
442 tolerance);
443
444
445 gkSource.generateSource(BORN, QUADRUPOLE, r2, bornI, bornK);
446 gkDipoleTensor.source(work);
447 double B10 = work[1];
448 double B11 = work[2];
449 double B12 = work[3];
450 gkDipoleTensor.generateTensor();
451 assertEquals("B001", z * B10, gkDipoleTensor.R001, tolerance);
452 assertEquals("B002", z * z * B11 + B10, gkDipoleTensor.R002, tolerance);
453 assertEquals("B003", z * z * z * B12 + 3.0 * z * B11, gkDipoleTensor.R003, tolerance);
454 }
455
456 @Test
457 public void quadrupoleTensorTest() {
458 double r2 = length2(r);
459 GKSource gkSource = new GKSource(4, gc);
460 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
461
462
463 int order = 5;
464 GKTensorQI gkQuadrupoleTensor = new GKTensorQI(QUADRUPOLE, order, gkSource, Eh, Es);
465 gkQuadrupoleTensor.setR(r);
466 double[] work = new double[order + 1];
467 gkQuadrupoleTensor.source(work);
468 double A20 = work[2];
469 double A21 = work[3];
470 double A22 = work[4];
471 double A23 = work[5];
472 gkQuadrupoleTensor.generateTensor();
473
474
475 assertEquals("R000", 0.0, gkQuadrupoleTensor.R000, tolerance);
476
477 assertEquals("R100", 0.0, gkQuadrupoleTensor.R100, tolerance);
478
479 double z = length(r);
480 double x = 0.0;
481 double y = 0.0;
482
483
484 assertEquals("R200", x * x * A20, gkQuadrupoleTensor.R200, tolerance);
485 assertEquals("R020", y * y * A20, gkQuadrupoleTensor.R020, tolerance);
486 assertEquals("R002", z * z * A20, gkQuadrupoleTensor.R002, tolerance);
487
488
489 assertEquals("R003", z * z * z * A21 + 3.0 * z * A20, gkQuadrupoleTensor.R003, tolerance);
490 assertEquals("R201", z * A20, gkQuadrupoleTensor.R201, tolerance);
491 assertEquals("R021", z * A20, gkQuadrupoleTensor.R021, tolerance);
492
493
494 assertEquals("R004", z * z * z * z * A22 + 6.0 * z * z * A21 + 3.0 * A20,
495 gkQuadrupoleTensor.R004, tolerance);
496 assertEquals("R005", z * z * z * z * z * A23 + 10.0 * z * z * z * A22 + 15.0 * z * A21,
497 gkQuadrupoleTensor.R005, tolerance);
498
499
500 assertEquals("R110", 0.0, gkQuadrupoleTensor.R110, tolerance);
501 assertEquals("R111", 0.0, gkQuadrupoleTensor.R111, tolerance);
502 assertEquals("R310", 0.0, gkQuadrupoleTensor.R310, tolerance);
503 assertEquals("R220", A20, gkQuadrupoleTensor.R220, tolerance);
504 assertEquals("R211", 0.0, gkQuadrupoleTensor.R211, tolerance);
505 assertEquals("R401", 3.0 * z * A21, gkQuadrupoleTensor.R401, tolerance);
506 assertEquals("R023", z * z * z * A22 + 3.0 * z * A21, gkQuadrupoleTensor.R023, tolerance);
507 assertEquals("R311", 0.0, gkQuadrupoleTensor.R311, tolerance);
508
509
510 gkSource.generateSource(BORN, QUADRUPOLE, r2, bornI, bornK);
511 gkQuadrupoleTensor.source(work);
512 double B20 = work[2];
513 double B21 = work[3];
514 double B22 = work[4];
515 gkQuadrupoleTensor.generateTensor();
516
517
518 assertEquals("B000", 0.0, gkQuadrupoleTensor.R000, tolerance);
519 assertEquals("B100", 0.0, gkQuadrupoleTensor.R100, tolerance);
520
521 assertEquals("B200", 0.0, gkQuadrupoleTensor.R200, tolerance);
522 assertEquals("B020", 0.0, gkQuadrupoleTensor.R020, tolerance);
523 assertEquals("B002", z * z * B20, gkQuadrupoleTensor.R002, tolerance);
524
525 assertEquals("B003", z * z * z * B21 + 3.0 * z * B20, gkQuadrupoleTensor.R003, tolerance);
526 assertEquals("B201", z * B20, gkQuadrupoleTensor.R201, tolerance);
527 assertEquals("B021", z * B20, gkQuadrupoleTensor.R021, tolerance);
528
529 assertEquals("B004", z * z * z * z * B22 + 6.0 * z * z * B21 + 3.0 * B20,
530 gkQuadrupoleTensor.R004, tolerance);
531 }
532
533 @Test
534 public void chargeFiniteDifferenceTest() {
535 double[] r = {0.7, 0.8, 0.9};
536 r[2] = length(r);
537 r[0] = 0.0;
538 r[1] = 0.0;
539 int order = 6;
540 GKSource gkSource = new GKSource(order, gc);
541 double r2 = length2(r);
542 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
543 GKTensorQI gkTensorQI = new GKTensorQI(MONOPOLE, order, gkSource, Eh, Es);
544
545 gkTensorQI.setR(r);
546 int tensorCount = MultipoleUtilities.tensorCount(order);
547 double[] tensor = new double[tensorCount];
548 gkTensorQI.noStorageRecursion(tensor);
549 double[] tensorsPz = new double[tensorCount];
550 double[] tensorsNz = new double[tensorCount];
551 double delta = 1.0e-5;
552 double delta2 = delta * 2;
553
554 r[2] += delta;
555 r2 = length2(r);
556 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
557 gkTensorQI.setR(r);
558 gkTensorQI.noStorageRecursion(tensorsPz);
559 r[2] -= delta2;
560 r2 = length2(r);
561 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
562 gkTensorQI.setR(r);
563 gkTensorQI.noStorageRecursion(tensorsNz);
564 r[2] += delta;
565
566 tensorFiniteDifference(gkTensorQI, delta2, 3, tensor, tensorsPz, tensorsNz);
567
568
569 r2 = length2(r);
570 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
571 gkTensorQI.setR(r);
572 gkTensorQI.recursion(tensor);
573 tensorFiniteDifference(gkTensorQI, delta2, 3, tensor, tensorsPz, tensorsNz);
574
575
576 gkTensorQI.setR(r);
577 gkTensorQI.generateTensor();
578 gkTensorQI.getTensor(tensor);
579 tensorFiniteDifference(gkTensorQI, delta2, 3, tensor, tensorsPz, tensorsNz);
580 }
581
582 @Test
583 public void dipoleFiniteDifferenceTest() {
584
585 double[] r = {0.7, 0.8, 0.9};
586 r[2] = length(r);
587 r[0] = 0.0;
588 r[1] = 0.0;
589 double r2 = length2(r);
590 int order = 6;
591 GKSource gkSource = new GKSource(order, gc);
592 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
593
594 GKTensorQI gkTensorQI = new GKTensorQI(DIPOLE, order, gkSource, Eh, Es);
595 gkTensorQI.setR(r);
596 int tensorCount = MultipoleUtilities.tensorCount(order);
597 double[] tensor = new double[tensorCount];
598 gkTensorQI.noStorageRecursion(tensor);
599 double[] tensorsPz = new double[tensorCount];
600 double[] tensorsNz = new double[tensorCount];
601 double delta = 1.0e-5;
602 double delta2 = delta * 2;
603
604 r[2] += delta;
605 r2 = length2(r);
606 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
607 gkTensorQI.setR(r);
608 gkTensorQI.noStorageRecursion(tensorsPz);
609 r[2] -= delta2;
610 r2 = length2(r);
611 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
612 gkTensorQI.setR(r);
613 gkTensorQI.noStorageRecursion(tensorsNz);
614 r[2] += delta;
615 tensorFiniteDifference(gkTensorQI, delta2, 3, tensor, tensorsPz, tensorsNz);
616
617
618 r2 = length2(r);
619 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, bornI, bornK);
620 gkTensorQI.setR(r);
621 gkTensorQI.recursion(tensor);
622 tensorFiniteDifference(gkTensorQI, delta2, 3, tensor, tensorsPz, tensorsNz);
623
624
625 gkTensorQI.setR(r);
626 gkTensorQI.generateTensor();
627 gkTensorQI.getTensor(tensor);
628 tensorFiniteDifference(gkTensorQI, delta2, 3, tensor, tensorsPz, tensorsNz);
629 }
630
631 private void tensorFiniteDifference(GKTensorQI gkTensorQI, double delta2, int order,
632 double[] tensor, double[] tensorsPz, double[] tensorsNz) {
633
634 int start = gkTensorQI.multipoleOrder.getOrder();
635
636 String info = "QK QI";
637
638
639 for (int l = start; l < order; l++) {
640
641 double expect = tensor[gkTensorQI.ti(l, 0, 1)];
642 double actual =
643 (tensorsPz[gkTensorQI.ti(l, 0, 0)] - tensorsNz[gkTensorQI.ti(l, 0, 0)]) / delta2;
644 assertEquals(info + " (d/dx): " + l, expect, actual, fdTolerance);
645 }
646 for (int l = 0; l < order; l++) {
647 for (int m = 1; m < order - l; m++) {
648
649 double expect = tensor[gkTensorQI.ti(l, m, 1)];
650 double actual =
651 (tensorsPz[gkTensorQI.ti(l, m, 0)] - tensorsNz[gkTensorQI.ti(l, m, 0)])
652 / delta2;
653 assertEquals(info + " (d/dx): " + l + " (d/dy): " + m, expect, actual, fdTolerance);
654 }
655 }
656
657 for (int l = 0; l < order; l++) {
658 for (int m = 0; m < order - l; m++) {
659 for (int n = 1; n < order - l - m; n++) {
660
661 double expect = tensor[gkTensorQI.ti(l, m, n + 1)];
662 double actual =
663 (tensorsPz[gkTensorQI.ti(l, m, n)] - tensorsNz[gkTensorQI.ti(l, m, n)])
664 / delta2;
665 assertEquals(info + " (d/dx): " + l + " (d/dy): " + m + " (d/dz):" + n, expect, actual,
666 fdTolerance);
667 }
668 }
669 }
670 }
671 }