View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * Test the GK tensor evaluated in the global coordinate frame.
54   * <p>
55   * There is no quadrupole finite-difference test because the trace of the quadrupole potential is
56   * neglected; thus the derivatives of the quadrupole potential are correct when summed over the
57   * trace, but not on a per-element basis.
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    // Mutual polarization.
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    // Direct polarization.
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     // Test the "bn" method.
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     // Monopole potential and derivatives.
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     // Monopole potential Born chain-rule derivatives.
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     // Monpole potential Born chain-rule derivatives.
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     // Dipole potential and derivatives.
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     // Dipole potential Born chain-rule derivatives.
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     // Quadrupole potential and derivatives.
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     // Quadrupole potential Born chain-rule derivatives.
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     // Monopole potential and derivatives.
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     // Check Born radii chain rule terms.
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     // Dipole potential and derivatives.
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     // No charge potential.
418     assertEquals(" R000", 0.0, gkDipoleTensor.R000, tolerance);
419     // Ux Dipole potential
420     assertEquals(" R100", x * A10, gkDipoleTensor.R100, tolerance);
421     // Ux Dipole potential gradient
422     assertEquals(" R200", x * x * A11 + A10, gkDipoleTensor.R200, tolerance);
423     // Ux Dipole 2nd potential gradient
424     assertEquals(" R300", x * x * x * A12 + 3.0 * x * A11, gkDipoleTensor.R300, tolerance);
425     // Ux Dipole 3rd potential gradient
426     assertEquals(" R400", x * x * x * x * A13 + 6.0 * x * x * A12 + 3.0 * A11, gkDipoleTensor.R400,
427         tolerance);
428 
429     // Check Born radii chain rule terms.
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     // Quadrupole potential and derivatives.
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     // No charge potential.
463     assertEquals(" R000", 0.0, gkQuadrupoleTensor.R000, tolerance);
464     // No dipole potential.
465     assertEquals(" R100", 0.0, gkQuadrupoleTensor.R100, tolerance);
466 
467     // Potential for the quadrupole trace
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     // Partial derivative of the trace with respect to X (note that the x*A20 term sums to 2*x*A20 over the trace).
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     // Higher order Qxx terms.
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     // Qxy
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     // Check Born radii chain rule terms.
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     // No charge or dipole potential.
508     assertEquals(" B000", 0.0, gkQuadrupoleTensor.R000, tolerance);
509     assertEquals(" B100", 0.0, gkQuadrupoleTensor.R100, tolerance);
510     // Born chain rule for the potential for the quadrupole trace
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     // Born chain rule for the partial derivative of the trace with respect to X.
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     // Higher order term for Qxx.
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     // Order(L^4) recursion.
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     // Machine generated code.
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     // Order(L^4) recursion.
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     // Machine generated code.
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     // Test the partial derivatives for all tensor components.
694     for (int l = start; l < order; l++) {
695       // Test X derivative
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       // Test Y derivative
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       // Test Z derivative
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         // Test X derivative
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         // Test Y derivative
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         // Test Z derivative
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     // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l + m + n = 1..o (n >= 1)
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           // Test X derivative
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           // Test Y derivative
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           // Test Z derivative
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 }