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