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