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.potential.nonbonded.implicit;
39  
40  import edu.rit.pj.IntegerForLoop;
41  import edu.rit.pj.ParallelRegion;
42  import edu.rit.pj.ParallelTeam;
43  import edu.rit.pj.reduction.SharedDouble;
44  import edu.rit.pj.reduction.SharedInteger;
45  import ffx.crystal.Crystal;
46  import ffx.crystal.SymOp;
47  import ffx.numerics.atomic.AtomicDoubleArray;
48  import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
49  import ffx.numerics.atomic.AtomicDoubleArray3D;
50  import ffx.numerics.multipole.GKEnergyQI;
51  import ffx.numerics.multipole.PolarizableMultipole;
52  import ffx.numerics.multipole.QIFrame;
53  import ffx.potential.bonded.Atom;
54  import ffx.potential.nonbonded.GeneralizedKirkwood.NonPolarModel;
55  import ffx.potential.nonbonded.pme.Polarization;
56  
57  import java.util.List;
58  import java.util.logging.Level;
59  import java.util.logging.Logger;
60  
61  import static ffx.numerics.multipole.GKSource.cn;
62  import static ffx.potential.parameters.MultipoleType.t000;
63  import static ffx.potential.parameters.MultipoleType.t001;
64  import static ffx.potential.parameters.MultipoleType.t002;
65  import static ffx.potential.parameters.MultipoleType.t010;
66  import static ffx.potential.parameters.MultipoleType.t011;
67  import static ffx.potential.parameters.MultipoleType.t020;
68  import static ffx.potential.parameters.MultipoleType.t100;
69  import static ffx.potential.parameters.MultipoleType.t101;
70  import static ffx.potential.parameters.MultipoleType.t110;
71  import static ffx.potential.parameters.MultipoleType.t200;
72  import static java.lang.String.format;
73  import static org.apache.commons.math3.util.FastMath.exp;
74  import static org.apache.commons.math3.util.FastMath.sqrt;
75  
76  /**
77   * Parallel calculation of the Generalized Kirkwood reaction field energy.
78   *
79   * @author Michael J. Schnieders
80   * @since 1.0
81   */
82  public class GKEnergyRegion extends ParallelRegion {
83  
84    private static final Logger logger = Logger.getLogger(GKEnergyRegion.class.getName());
85    /** Constant factor used with quadrupoles. */
86    protected static final double oneThird = 1.0 / 3.0;
87    /** Conversion from electron**2/Ang to kcal/mole. */
88    protected final double electric;
89    /** Treatment of polarization. */
90    protected final Polarization polarization;
91    protected final NonPolarModel nonPolarModel;
92  
93    private final double soluteDielectric;
94    private final double solventDielectric;
95  
96    /**
97     * Dielectric offset from:
98     *
99     * <p>W. C. Still, A. Tempczyk, R. C. Hawley and T. Hendrickson, "A Semianalytical Treatment of
100    * Solvation for Molecular Mechanics and Dynamics", J. Amer. Chem. Soc., 112, 6127-6129 (1990)
101    */
102   protected final double dOffset = 0.09;
103   /** Cavitation surface tension coefficient (kcal/mol/A^2). */
104   protected final double surfaceTension;
105   /** Empirical constant that controls the GK cross-term. */
106   protected final double gkc;
107   /** Kirkwood monopole reaction field constant. */
108   private final double fc;
109   /** Kirkwood dipole reaction field constant. */
110   private final double fd;
111   /** Kirkwood quadrupole reaction field constant. */
112   private final double fq;
113   /** Water probe radius. */
114   private final double probe;
115   private final SharedDouble sharedPermanentGKEnergy;
116   private final SharedDouble sharedPolarizationGKEnergy;
117   private final SharedDouble sharedGKEnergy;
118   private final SharedInteger sharedInteractions;
119   private final IntegerForLoop[] gkEnergyLoop;
120   /** An ordered array of atoms in the system. */
121   private Atom[] atoms;
122   /** Induced dipoles for each symmetry operator. */
123   private double[][][] inducedDipole;
124   /** Induced dipole chain rule terms for each symmetry operator. */
125   private double[][][] inducedDipoleCR;
126   /** Multipole moments for each symmetry operator. */
127   private double[][][] globalMultipole;
128   /** Periodic boundary conditions and symmetry. */
129   private Crystal crystal;
130   /** Atomic coordinates for each symmetry operator. */
131   private double[][][] sXYZ;
132   /** Neighbor lists for each atom and symmetry operator. */
133   private int[][][] neighborLists;
134   /** Flag to indicate if an atom should be included. */
135   private boolean[] use = null;
136   /** GK cut-off distance squared. */
137   private double cut2;
138   /** Base radius of each atom (for Born radii based nonpolar energy). */
139   private double[] baseRadius;
140   /** Born radius of each atom. */
141   private double[] born;
142 
143   private boolean gradient = false;
144   /** Atomic 3D Gradient array. */
145   private AtomicDoubleArray3D grad;
146   /** Atomic 3D Torque array. */
147   private AtomicDoubleArray3D torque;
148   /** Shared array for computation of Born radii gradient. */
149   private AtomicDoubleArray sharedBornGrad;
150   /** Self-energy for each atom */
151   private AtomicDoubleArray selfEnergy;
152   /** Cross-term energy for each atom */
153   private AtomicDoubleArray crossEnergy;
154 
155 
156   public GKEnergyRegion(
157       int nt, Polarization polarization, NonPolarModel nonPolarModel,
158       double surfaceTension, double probe, double electric,
159       double soluteDieletric, double solventDieletric, double gkc, boolean gkQI) {
160 
161     // Set the Kirkwood multipolar reaction field constants.
162     this.soluteDielectric = soluteDieletric;
163     this.solventDielectric = solventDieletric;
164     fc = cn(0, soluteDieletric, solventDieletric);
165     fd = cn(1, soluteDieletric, solventDieletric);
166     fq = cn(2, soluteDieletric, solventDieletric);
167     this.gkc = gkc;
168 
169     this.polarization = polarization;
170     this.nonPolarModel = nonPolarModel;
171     this.surfaceTension = surfaceTension;
172     this.probe = probe;
173     // Set the conversion from electron**2/Ang to kcal/mole
174     this.electric = electric;
175 
176     gkEnergyLoop = new IntegerForLoop[nt];
177     for (int i = 0; i < nt; i++) {
178       if (gkQI) {
179         gkEnergyLoop[i] = new GKEnergyLoopQI();
180       } else {
181         gkEnergyLoop[i] = new GKEnergyLoop();
182       }
183     }
184     sharedPermanentGKEnergy = new SharedDouble();
185     sharedPolarizationGKEnergy = new SharedDouble();
186     sharedGKEnergy = new SharedDouble();
187     sharedInteractions = new SharedInteger();
188   }
189 
190   public double getEnergy() {
191     return sharedGKEnergy.get();
192   }
193 
194   public double getPermanentEnergy() {
195     return sharedPermanentGKEnergy.get();
196   }
197 
198   public double getPolarizationEnergy() {
199     return sharedPolarizationGKEnergy.get();
200   }
201 
202   public AtomicDoubleArray getSelfEnergy() {
203     int nAtoms = atoms.length;
204     selfEnergy.reduce(0, nAtoms - 1);
205     return selfEnergy;
206   }
207 
208   public int getInteractions() {
209     return sharedInteractions.get();
210   }
211 
212   public void init(
213       Atom[] atoms,
214       double[][][] globalMultipole,
215       double[][][] inducedDipole,
216       double[][][] inducedDipoleCR,
217       Crystal crystal,
218       double[][][] sXYZ,
219       int[][][] neighborLists,
220       boolean[] use,
221       double cut2,
222       double[] baseRadius,
223       double[] born,
224       boolean gradient,
225       ParallelTeam parallelTeam,
226       AtomicDoubleArray3D grad,
227       AtomicDoubleArray3D torque,
228       AtomicDoubleArray sharedBornGrad) {
229     // Input
230     this.atoms = atoms;
231     this.globalMultipole = globalMultipole;
232     this.inducedDipole = inducedDipole;
233     this.inducedDipoleCR = inducedDipoleCR;
234     this.crystal = crystal;
235     this.sXYZ = sXYZ;
236     this.neighborLists = neighborLists;
237     this.use = use;
238     this.cut2 = cut2;
239     this.baseRadius = baseRadius;
240     this.born = born;
241     this.gradient = gradient;
242     // Output
243     this.grad = grad;
244     this.torque = torque;
245     this.sharedBornGrad = sharedBornGrad;
246 
247     int nAtoms = atoms.length;
248     int nThreads = gkEnergyLoop.length;
249     if (selfEnergy == null || selfEnergy.size() != atoms.length) {
250       selfEnergy = AtomicDoubleArray.atomicDoubleArrayFactory(
251           AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
252       crossEnergy = AtomicDoubleArray.atomicDoubleArrayFactory(
253           AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
254     } else {
255       selfEnergy.reset(parallelTeam, 0, nAtoms - 1);
256       crossEnergy.reset(parallelTeam, 0, nAtoms - 1);
257     }
258   }
259 
260   @Override
261   public void run() {
262     try {
263       int nAtoms = atoms.length;
264       int threadIndex = getThreadIndex();
265       execute(0, nAtoms - 1, gkEnergyLoop[threadIndex]);
266     } catch (Exception e) {
267       String message = "Fatal exception computing GK Energy in thread " + getThreadIndex() + "\n";
268       logger.log(Level.SEVERE, message, e);
269     }
270   }
271 
272   @Override
273   public void start() {
274     sharedPermanentGKEnergy.set(0.0);
275     sharedPolarizationGKEnergy.set(0.0);
276     sharedGKEnergy.set(0.0);
277     sharedInteractions.set(0);
278   }
279 
280   @Override
281   public void finish() {
282     if (logger.isLoggable(Level.FINER)) {
283       int nAtoms = atoms.length;
284       selfEnergy.reduce(0, nAtoms - 1);
285       crossEnergy.reduce(0, nAtoms - 1);
286       logger.finer(" Generalized Kirkwood Self-Energies and Cross-Energies\n");
287       double selfSum = 0.0;
288       double crossSum = 0.0;
289       for (int i = 0; i < nAtoms; i++) {
290         double self = selfEnergy.get(i);
291         double cross = crossEnergy.get(i);
292         logger.finer(format("GKSELF   %5d %16.8f %16.8f", i, self, cross));
293         selfSum += self;
294         crossSum += cross;
295       }
296       logger.finer(format("       %16.8f %16.8f %16.8f\n",
297           selfSum, crossSum, selfSum + crossSum));
298     }
299   }
300 
301   /**
302    * Compute the GK Energy.
303    *
304    * @since 1.0
305    */
306   private class GKEnergyLoop extends IntegerForLoop {
307 
308     private final double[][] a;
309     private final double[][] b;
310     private final double[] gc;
311     private final double[] gux;
312     private final double[] guy;
313     private final double[] guz;
314     private final double[] gqxx;
315     private final double[] gqyy;
316     private final double[] gqzz;
317     private final double[] gqxy;
318     private final double[] gqxz;
319     private final double[] gqyz;
320     private final double[] dx_local;
321     private double ci, uxi, uyi, uzi, qxxi, qxyi, qxzi, qyyi, qyzi, qzzi;
322     private double ck, uxk, uyk, uzk, qxxk, qxyk, qxzk, qyyk, qyzk, qzzk;
323     private double dxi, dyi, dzi, pxi, pyi, pzi, sxi, syi, szi;
324     private double dxk, dyk, dzk, pxk, pyk, pzk, sxk, syk, szk;
325     private double xr, yr, zr, xr2, yr2, zr2, rbi, rbk;
326     private double xi, yi, zi;
327     private double dedxi, dedyi, dedzi;
328     private double dborni;
329     private double trqxi, trqyi, trqzi;
330     private int count;
331     private int iSymm;
332     private int threadID;
333     private final double[][] transOp;
334     private double gkPermanentEnergy;
335     private double gkPolarizationEnergy;
336     private double gkEnergy;
337 
338     GKEnergyLoop() {
339       a = new double[6][4];
340       b = new double[5][3];
341       gc = new double[31];
342       gux = new double[31];
343       guy = new double[31];
344       guz = new double[31];
345       gqxx = new double[31];
346       gqyy = new double[31];
347       gqzz = new double[31];
348       gqxy = new double[31];
349       gqxz = new double[31];
350       gqyz = new double[31];
351       dx_local = new double[3];
352       transOp = new double[3][3];
353     }
354 
355     @Override
356     public void finish() {
357       sharedInteractions.addAndGet(count);
358       sharedGKEnergy.addAndGet(gkEnergy);
359       sharedPermanentGKEnergy.addAndGet(gkPermanentEnergy);
360       sharedPolarizationGKEnergy.addAndGet(gkPolarizationEnergy);
361     }
362 
363     @Override
364     public void run(int lb, int ub) {
365 
366       double[] x = sXYZ[0][0];
367       double[] y = sXYZ[0][1];
368       double[] z = sXYZ[0][2];
369 
370       int nSymm = crystal.spaceGroup.symOps.size();
371       List<SymOp> symOps = crystal.spaceGroup.symOps;
372 
373       for (iSymm = 0; iSymm < nSymm; iSymm++) {
374         SymOp symOp = symOps.get(iSymm);
375         crystal.getTransformationOperator(symOp, transOp);
376         for (int i = lb; i <= ub; i++) {
377           if (!use[i]) {
378             continue;
379           }
380 
381           // Zero out force accumulation for atom i.
382           dedxi = 0.0;
383           dedyi = 0.0;
384           dedzi = 0.0;
385           dborni = 0.0;
386           trqxi = 0.0;
387           trqyi = 0.0;
388           trqzi = 0.0;
389 
390           xi = x[i];
391           yi = y[i];
392           zi = z[i];
393           final double[] multipolei = globalMultipole[0][i];
394           ci = multipolei[t000];
395           uxi = multipolei[t100];
396           uyi = multipolei[t010];
397           uzi = multipolei[t001];
398           qxxi = multipolei[t200] * oneThird;
399           qyyi = multipolei[t020] * oneThird;
400           qzzi = multipolei[t002] * oneThird;
401           qxyi = multipolei[t110] * oneThird;
402           qxzi = multipolei[t101] * oneThird;
403           qyzi = multipolei[t011] * oneThird;
404           dxi = inducedDipole[0][i][0];
405           dyi = inducedDipole[0][i][1];
406           dzi = inducedDipole[0][i][2];
407           pxi = inducedDipoleCR[0][i][0];
408           pyi = inducedDipoleCR[0][i][1];
409           pzi = inducedDipoleCR[0][i][2];
410           sxi = dxi + pxi;
411           syi = dyi + pyi;
412           szi = dzi + pzi;
413           rbi = born[i];
414           int[] list = neighborLists[iSymm][i];
415           for (int k : list) {
416             if (!use[k]) {
417               continue;
418             }
419             interaction(i, k);
420           }
421           if (iSymm == 0) {
422             // Include self-interactions for the asymmetric unit atoms.
423             interaction(i, i);
424             /*
425              Formula for Born energy approximation for cavitation energy is:
426              e = surfaceTension / 6 * (ri + probe)^2 * (ri/rb)^6.
427              ri is the base atomic radius the atom.
428              rb is Born radius of the atom.
429             */
430             switch (nonPolarModel) {
431               case BORN_SOLV:
432               case BORN_CAV_DISP:
433                 double r = baseRadius[i] + dOffset + probe;
434                 double ratio = (baseRadius[i] + dOffset) / born[i];
435                 ratio *= ratio;
436                 ratio *= (ratio * ratio);
437                 double saTerm = surfaceTension * r * r * ratio / 6.0;
438                 gkEnergy += saTerm;
439                 sharedBornGrad.sub(threadID, i, 6.0 * saTerm / born[i]);
440                 break;
441               default:
442                 break;
443             }
444           }
445           if (gradient) {
446             grad.add(threadID, i, dedxi, dedyi, dedzi);
447             torque.add(threadID, i, trqxi, trqyi, trqzi);
448             sharedBornGrad.add(threadID, i, dborni);
449           }
450         }
451       }
452     }
453 
454     @Override
455     public void start() {
456       gkEnergy = 0.0;
457       gkPermanentEnergy = 0.0;
458       gkPolarizationEnergy = 0.0;
459       count = 0;
460       threadID = getThreadIndex();
461     }
462 
463     private void interaction(int i, int k) {
464       dx_local[0] = sXYZ[iSymm][0][k] - xi;
465       dx_local[1] = sXYZ[iSymm][1][k] - yi;
466       dx_local[2] = sXYZ[iSymm][2][k] - zi;
467       double r2 = crystal.image(dx_local);
468       if (r2 > cut2) {
469         return;
470       }
471       xr = dx_local[0];
472       yr = dx_local[1];
473       zr = dx_local[2];
474       xr2 = xr * xr;
475       yr2 = yr * yr;
476       zr2 = zr * zr;
477       rbk = born[k];
478 
479       final double[] multipolek = globalMultipole[iSymm][k];
480       ck = multipolek[t000];
481       uxk = multipolek[t100];
482       uyk = multipolek[t010];
483       uzk = multipolek[t001];
484       qxxk = multipolek[t200] * oneThird;
485       qyyk = multipolek[t020] * oneThird;
486       qzzk = multipolek[t002] * oneThird;
487       qxyk = multipolek[t110] * oneThird;
488       qxzk = multipolek[t101] * oneThird;
489       qyzk = multipolek[t011] * oneThird;
490       dxk = inducedDipole[iSymm][k][0];
491       dyk = inducedDipole[iSymm][k][1];
492       dzk = inducedDipole[iSymm][k][2];
493       pxk = inducedDipoleCR[iSymm][k][0];
494       pyk = inducedDipoleCR[iSymm][k][1];
495       pzk = inducedDipoleCR[iSymm][k][2];
496       sxk = dxk + pxk;
497       syk = dyk + pyk;
498       szk = dzk + pzk;
499       final double rb2 = rbi * rbk;
500       final double expterm = exp(-r2 / (gkc * rb2));
501       final double expc = expterm / gkc;
502       final double expc1 = 1.0 - expc;
503       final double expcr = r2 * expterm / (gkc * gkc * rb2 * rb2);
504       final double dexpc = -2.0 / (gkc * rb2);
505       double expcdexpc = -expc * dexpc;
506       final double dexpcr = 2.0 / (gkc * rb2 * rb2);
507       final double dgfdr = 0.5 * expterm * (1.0 + r2 / (rb2 * gkc));
508       final double gf2 = 1.0 / (r2 + rb2 * expterm);
509       final double gf = sqrt(gf2);
510       final double gf3 = gf2 * gf;
511       final double gf5 = gf3 * gf2;
512       final double gf7 = gf5 * gf2;
513       final double gf9 = gf7 * gf2;
514       final double gf11 = gf9 * gf2;
515 
516       // Reaction potential auxiliary terms.
517       a[0][0] = gf;
518       a[1][0] = -gf3;
519       a[2][0] = 3.0 * gf5;
520       a[3][0] = -15.0 * gf7;
521       a[4][0] = 105.0 * gf9;
522       a[5][0] = -945.0 * gf11;
523 
524       // Reaction potential gradient auxiliary terms.
525       a[0][1] = expc1 * a[1][0];
526       a[1][1] = expc1 * a[2][0];
527       a[2][1] = expc1 * a[3][0];
528       a[3][1] = expc1 * a[4][0];
529       a[4][1] = expc1 * a[5][0];
530 
531       // 2nd reaction potential gradient auxiliary terms.
532       a[0][2] = expc1 * a[1][1] + expcdexpc * a[1][0];
533       a[1][2] = expc1 * a[2][1] + expcdexpc * a[2][0];
534       a[2][2] = expc1 * a[3][1] + expcdexpc * a[3][0];
535       a[3][2] = expc1 * a[4][1] + expcdexpc * a[4][0];
536 
537       if (gradient) {
538 
539         // 3rd reaction potential gradient auxiliary terms.
540         expcdexpc = 2.0 * expcdexpc;
541         a[0][3] = expc1 * a[1][2] + expcdexpc * a[1][1];
542         a[1][3] = expc1 * a[2][2] + expcdexpc * a[2][1];
543         a[2][3] = expc1 * a[3][2] + expcdexpc * a[3][1];
544         expcdexpc = -expc * dexpc * dexpc;
545         a[0][3] = a[0][3] + expcdexpc * a[1][0];
546         a[1][3] = a[1][3] + expcdexpc * a[2][0];
547         a[2][3] = a[2][3] + expcdexpc * a[3][0];
548 
549         // Born radii derivatives of reaction potential auxiliary terms.
550         b[0][0] = dgfdr * a[1][0];
551         b[1][0] = dgfdr * a[2][0];
552         b[2][0] = dgfdr * a[3][0];
553         b[3][0] = dgfdr * a[4][0];
554         b[4][0] = dgfdr * a[5][0];
555 
556         // Born radii gradients of reaction potential gradient auxiliary terms.
557         b[0][1] = b[1][0] - expcr * a[1][0] - expc * b[1][0];
558         b[1][1] = b[2][0] - expcr * a[2][0] - expc * b[2][0];
559         b[2][1] = b[3][0] - expcr * a[3][0] - expc * b[3][0];
560         b[3][1] = b[4][0] - expcr * a[4][0] - expc * b[4][0];
561 
562         // Born radii derivatives of the 2nd reaction potential gradient auxiliary terms.
563         b[0][2] = b[1][1] - (expcr * (a[1][1] + dexpc * a[1][0])
564             + expc * (b[1][1] + dexpcr * a[1][0] + dexpc * b[1][0]));
565         b[1][2] = b[2][1] - (expcr * (a[2][1] + dexpc * a[2][0])
566             + expc * (b[2][1] + dexpcr * a[2][0] + dexpc * b[2][0]));
567         b[2][2] = b[3][1] - (expcr * (a[3][1] + dexpc * a[3][0])
568             + expc * (b[3][1] + dexpcr * a[3][0] + dexpc * b[3][0]));
569 
570         // Multiply the Born radii auxiliary terms by their dielectric functions.
571         b[0][0] = electric * fc * b[0][0];
572         b[0][1] = electric * fc * b[0][1];
573         b[0][2] = electric * fc * b[0][2];
574         b[1][0] = electric * fd * b[1][0];
575         b[1][1] = electric * fd * b[1][1];
576         b[1][2] = electric * fd * b[1][2];
577         b[2][0] = electric * fq * b[2][0];
578         b[2][1] = electric * fq * b[2][1];
579         b[2][2] = electric * fq * b[2][2];
580       }
581 
582       // Multiply the potential auxiliary terms by their dielectric functions.
583       a[0][0] = electric * fc * a[0][0];
584       a[0][1] = electric * fc * a[0][1];
585       a[0][2] = electric * fc * a[0][2];
586       a[0][3] = electric * fc * a[0][3];
587       a[1][0] = electric * fd * a[1][0];
588       a[1][1] = electric * fd * a[1][1];
589       a[1][2] = electric * fd * a[1][2];
590       a[1][3] = electric * fd * a[1][3];
591       a[2][0] = electric * fq * a[2][0];
592       a[2][1] = electric * fq * a[2][1];
593       a[2][2] = electric * fq * a[2][2];
594       a[2][3] = electric * fq * a[2][3];
595 
596       // Compute the GK tensors required to compute the energy.
597       energyTensors();
598 
599       // Compute the GK interaction energy.
600       double eik = energy(i, k);
601       gkEnergy += eik;
602       count++;
603 
604       // List<Atom> i12 = atoms[i].get12List();
605       // List<Atom> i13 = atoms[i].get13List();
606       // List<Atom> i14 = atoms[i].get14List();
607       // List<Atom> i15 = atoms[i].get15List();
608       // Atom aK = atoms[k];
609       // if (i != k && !i12.contains(aK) && !i13.contains(aK) && !i14.contains(aK) && !i15.contains(aK)) {
610       //   logger.info(format(" GK %s %7.4f %s %7.4f %7.4f %7.4f", atoms[i], rbi, atoms[k], rbk, sqrt(r2), eik));
611       // }
612 
613       if (gradient) {
614         // Compute the additional GK tensors required to compute the energy gradient.
615         gradientTensors();
616 
617         // Compute the permanent GK energy gradient.
618         permanentEnergyGradient(i, k);
619 
620         if (polarization != Polarization.NONE) {
621           // Compute the induced GK energy gradient.
622           polarizationEnergyGradient(i, k);
623         }
624       }
625     }
626 
627     private void energyTensors() {
628       // Unweighted reaction potential tensor.
629       gc[1] = a[0][0];
630       gux[1] = xr * a[1][0];
631       guy[1] = yr * a[1][0];
632       guz[1] = zr * a[1][0];
633       gqxx[1] = xr2 * a[2][0];
634       gqyy[1] = yr2 * a[2][0];
635       gqzz[1] = zr2 * a[2][0];
636       gqxy[1] = xr * yr * a[2][0];
637       gqxz[1] = xr * zr * a[2][0];
638       gqyz[1] = yr * zr * a[2][0];
639 
640       // Unweighted reaction potential gradient tensor.
641       gc[2] = xr * a[0][1];
642       gc[3] = yr * a[0][1];
643       gc[4] = zr * a[0][1];
644       gux[2] = a[1][0] + xr2 * a[1][1];
645       gux[3] = xr * yr * a[1][1];
646       gux[4] = xr * zr * a[1][1];
647       guy[2] = gux[3];
648       guy[3] = a[1][0] + yr2 * a[1][1];
649       guy[4] = yr * zr * a[1][1];
650       guz[2] = gux[4];
651       guz[3] = guy[4];
652       guz[4] = a[1][0] + zr2 * a[1][1];
653       gqxx[2] = xr * (2.0 * a[2][0] + xr2 * a[2][1]);
654       gqxx[3] = yr * xr2 * a[2][1];
655       gqxx[4] = zr * xr2 * a[2][1];
656       gqyy[2] = xr * yr2 * a[2][1];
657       gqyy[3] = yr * (2.0 * a[2][0] + yr2 * a[2][1]);
658       gqyy[4] = zr * yr2 * a[2][1];
659       gqzz[2] = xr * zr2 * a[2][1];
660       gqzz[3] = yr * zr2 * a[2][1];
661       gqzz[4] = zr * (2.0 * a[2][0] + zr2 * a[2][1]);
662       gqxy[2] = yr * (a[2][0] + xr2 * a[2][1]);
663       gqxy[3] = xr * (a[2][0] + yr2 * a[2][1]);
664       gqxy[4] = zr * xr * yr * a[2][1];
665       gqxz[2] = zr * (a[2][0] + xr2 * a[2][1]);
666       gqxz[3] = gqxy[4];
667       gqxz[4] = xr * (a[2][0] + zr2 * a[2][1]);
668       gqyz[2] = gqxy[4];
669       gqyz[3] = zr * (a[2][0] + yr2 * a[2][1]);
670       gqyz[4] = yr * (a[2][0] + zr2 * a[2][1]);
671 
672       // Unweighted 2nd reaction potential gradient tensor.
673       gc[5] = a[0][1] + xr2 * a[0][2];
674       gc[6] = xr * yr * a[0][2];
675       gc[7] = xr * zr * a[0][2];
676       gc[8] = a[0][1] + yr2 * a[0][2];
677       gc[9] = yr * zr * a[0][2];
678       gc[10] = a[0][1] + zr2 * a[0][2];
679       gux[5] = xr * (3.0 * a[1][1] + xr2 * a[1][2]);
680       gux[6] = yr * (a[1][1] + xr2 * a[1][2]);
681       gux[7] = zr * (a[1][1] + xr2 * a[1][2]);
682       gux[8] = xr * (a[1][1] + yr2 * a[1][2]);
683       gux[9] = zr * xr * yr * a[1][2];
684       gux[10] = xr * (a[1][1] + zr2 * a[1][2]);
685       guy[5] = yr * (a[1][1] + xr2 * a[1][2]);
686       guy[6] = xr * (a[1][1] + yr2 * a[1][2]);
687       guy[7] = gux[9];
688       guy[8] = yr * (3.0 * a[1][1] + yr2 * a[1][2]);
689       guy[9] = zr * (a[1][1] + yr2 * a[1][2]);
690       guy[10] = yr * (a[1][1] + zr2 * a[1][2]);
691       guz[5] = zr * (a[1][1] + xr2 * a[1][2]);
692       guz[6] = gux[9];
693       guz[7] = xr * (a[1][1] + zr2 * a[1][2]);
694       guz[8] = zr * (a[1][1] + yr2 * a[1][2]);
695       guz[9] = yr * (a[1][1] + zr2 * a[1][2]);
696       guz[10] = zr * (3.0 * a[1][1] + zr2 * a[1][2]);
697       gqxx[5] = 2.0 * a[2][0] + xr2 * (5.0 * a[2][1] + xr2 * a[2][2]);
698       gqxx[6] = yr * xr * (2.0 * a[2][1] + xr2 * a[2][2]);
699       gqxx[7] = zr * xr * (2.0 * a[2][1] + xr2 * a[2][2]);
700       gqxx[8] = xr2 * (a[2][1] + yr2 * a[2][2]);
701       gqxx[9] = zr * yr * xr2 * a[2][2];
702       gqxx[10] = xr2 * (a[2][1] + zr2 * a[2][2]);
703       gqyy[5] = yr2 * (a[2][1] + xr2 * a[2][2]);
704       gqyy[6] = xr * yr * (2.0 * a[2][1] + yr2 * a[2][2]);
705       gqyy[7] = xr * zr * yr2 * a[2][2];
706       gqyy[8] = 2.0 * a[2][0] + yr2 * (5.0 * a[2][1] + yr2 * a[2][2]);
707       gqyy[9] = yr * zr * (2.0 * a[2][1] + yr2 * a[2][2]);
708       gqyy[10] = yr2 * (a[2][1] + zr2 * a[2][2]);
709       gqzz[5] = zr2 * (a[2][1] + xr2 * a[2][2]);
710       gqzz[6] = xr * yr * zr2 * a[2][2];
711       gqzz[7] = xr * zr * (2.0 * a[2][1] + zr2 * a[2][2]);
712       gqzz[8] = zr2 * (a[2][1] + yr2 * a[2][2]);
713       gqzz[9] = yr * zr * (2.0 * a[2][1] + zr2 * a[2][2]);
714       gqzz[10] = 2.0 * a[2][0] + zr2 * (5.0 * a[2][1] + zr2 * a[2][2]);
715       gqxy[5] = xr * yr * (3.0 * a[2][1] + xr2 * a[2][2]);
716       gqxy[6] = a[2][0] + (xr2 + yr2) * a[2][1] + xr2 * yr2 * a[2][2];
717       gqxy[7] = zr * yr * (a[2][1] + xr2 * a[2][2]);
718       gqxy[8] = xr * yr * (3.0 * a[2][1] + yr2 * a[2][2]);
719       gqxy[9] = zr * xr * (a[2][1] + yr2 * a[2][2]);
720       gqxy[10] = xr * yr * (a[2][1] + zr2 * a[2][2]);
721       gqxz[5] = xr * zr * (3.0 * a[2][1] + xr2 * a[2][2]);
722       gqxz[6] = yr * zr * (a[2][1] + xr2 * a[2][2]);
723       gqxz[7] = a[2][0] + (xr2 + zr2) * a[2][1] + xr2 * zr2 * a[2][2];
724       gqxz[8] = xr * zr * (a[2][1] + yr2 * a[2][2]);
725       gqxz[9] = xr * yr * (a[2][1] + zr2 * a[2][2]);
726       gqxz[10] = xr * zr * (3.0 * a[2][1] + zr2 * a[2][2]);
727       gqyz[5] = zr * yr * (a[2][1] + xr2 * a[2][2]);
728       gqyz[6] = xr * zr * (a[2][1] + yr2 * a[2][2]);
729       gqyz[7] = xr * yr * (a[2][1] + zr2 * a[2][2]);
730       gqyz[8] = yr * zr * (3.0 * a[2][1] + yr2 * a[2][2]);
731       gqyz[9] = a[2][0] + (yr2 + zr2) * a[2][1] + yr2 * zr2 * a[2][2];
732       gqyz[10] = yr * zr * (3.0 * a[2][1] + zr2 * a[2][2]);
733     }
734 
735     private double energy(int i, int k) {
736       // Electrostatic solvation energy of the permanent multipoles in their own GK reaction
737       // potential.
738       double esym =
739           ci * ck * gc[1]
740               - (uxi * (uxk * gux[2] + uyk * guy[2] + uzk * guz[2])
741               + uyi * (uxk * gux[3] + uyk * guy[3] + uzk * guz[3])
742               + uzi * (uxk * gux[4] + uyk * guy[4] + uzk * guz[4]));
743       double ewi =
744           ci * (uxk * gc[2] + uyk * gc[3] + uzk * gc[4])
745               - ck * (uxi * gux[1] + uyi * guy[1] + uzi * guz[1])
746               + ci
747               * (qxxk * gc[5]
748               + qyyk * gc[8]
749               + qzzk * gc[10]
750               + 2.0 * (qxyk * gc[6] + qxzk * gc[7] + qyzk * gc[9]))
751               + ck
752               * (qxxi * gqxx[1]
753               + qyyi * gqyy[1]
754               + qzzi * gqzz[1]
755               + 2.0 * (qxyi * gqxy[1] + qxzi * gqxz[1] + qyzi * gqyz[1]))
756               - uxi
757               * (qxxk * gux[5]
758               + qyyk * gux[8]
759               + qzzk * gux[10]
760               + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
761               - uyi
762               * (qxxk * guy[5]
763               + qyyk * guy[8]
764               + qzzk * guy[10]
765               + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
766               - uzi
767               * (qxxk * guz[5]
768               + qyyk * guz[8]
769               + qzzk * guz[10]
770               + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
771               + uxk
772               * (qxxi * gqxx[2]
773               + qyyi * gqyy[2]
774               + qzzi * gqzz[2]
775               + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
776               + uyk
777               * (qxxi * gqxx[3]
778               + qyyi * gqyy[3]
779               + qzzi * gqzz[3]
780               + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
781               + uzk
782               * (qxxi * gqxx[4]
783               + qyyi * gqyy[4]
784               + qzzi * gqzz[4]
785               + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]))
786               + qxxi
787               * (qxxk * gqxx[5]
788               + qyyk * gqxx[8]
789               + qzzk * gqxx[10]
790               + 2.0 * (qxyk * gqxx[6] + qxzk * gqxx[7] + qyzk * gqxx[9]))
791               + qyyi
792               * (qxxk * gqyy[5]
793               + qyyk * gqyy[8]
794               + qzzk * gqyy[10]
795               + 2.0 * (qxyk * gqyy[6] + qxzk * gqyy[7] + qyzk * gqyy[9]))
796               + qzzi
797               * (qxxk * gqzz[5]
798               + qyyk * gqzz[8]
799               + qzzk * gqzz[10]
800               + 2.0 * (qxyk * gqzz[6] + qxzk * gqzz[7] + qyzk * gqzz[9]))
801               + 2.0
802               * (qxyi
803               * (qxxk * gqxy[5]
804               + qyyk * gqxy[8]
805               + qzzk * gqxy[10]
806               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxy[7] + qyzk * gqxy[9]))
807               + qxzi
808               * (qxxk * gqxz[5]
809               + qyyk * gqxz[8]
810               + qzzk * gqxz[10]
811               + 2.0 * (qxyk * gqxz[6] + qxzk * gqxz[7] + qyzk * gqxz[9]))
812               + qyzi
813               * (qxxk * gqyz[5]
814               + qyyk * gqyz[8]
815               + qzzk * gqyz[10]
816               + 2.0 * (qxyk * gqyz[6] + qxzk * gqyz[7] + qyzk * gqyz[9])));
817       double ewk =
818           ci * (uxk * gux[1] + uyk * guy[1] + uzk * guz[1])
819               - ck * (uxi * gc[2] + uyi * gc[3] + uzi * gc[4])
820               + ci
821               * (qxxk * gqxx[1]
822               + qyyk * gqyy[1]
823               + qzzk * gqzz[1]
824               + 2.0 * (qxyk * gqxy[1] + qxzk * gqxz[1] + qyzk * gqyz[1]))
825               + ck
826               * (qxxi * gc[5]
827               + qyyi * gc[8]
828               + qzzi * gc[10]
829               + 2.0 * (qxyi * gc[6] + qxzi * gc[7] + qyzi * gc[9]))
830               - uxi
831               * (qxxk * gqxx[2]
832               + qyyk * gqyy[2]
833               + qzzk * gqzz[2]
834               + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
835               - uyi
836               * (qxxk * gqxx[3]
837               + qyyk * gqyy[3]
838               + qzzk * gqzz[3]
839               + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
840               - uzi
841               * (qxxk * gqxx[4]
842               + qyyk * gqyy[4]
843               + qzzk * gqzz[4]
844               + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
845               + uxk
846               * (qxxi * gux[5]
847               + qyyi * gux[8]
848               + qzzi * gux[10]
849               + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
850               + uyk
851               * (qxxi * guy[5]
852               + qyyi * guy[8]
853               + qzzi * guy[10]
854               + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
855               + uzk
856               * (qxxi * guz[5]
857               + qyyi * guz[8]
858               + qzzi * guz[10]
859               + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]))
860               + qxxi
861               * (qxxk * gqxx[5]
862               + qyyk * gqyy[5]
863               + qzzk * gqzz[5]
864               + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
865               + qyyi
866               * (qxxk * gqxx[8]
867               + qyyk * gqyy[8]
868               + qzzk * gqzz[8]
869               + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
870               + qzzi
871               * (qxxk * gqxx[10]
872               + qyyk * gqyy[10]
873               + qzzk * gqzz[10]
874               + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
875               + 2.0
876               * (qxyi
877               * (qxxk * gqxx[6]
878               + qyyk * gqyy[6]
879               + qzzk * gqzz[6]
880               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
881               + qxzi
882               * (qxxk * gqxx[7]
883               + qyyk * gqyy[7]
884               + qzzk * gqzz[7]
885               + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
886               + qyzi
887               * (qxxk * gqxx[9]
888               + qyyk * gqyy[9]
889               + qzzk * gqzz[9]
890               + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9])));
891       double e = esym + 0.5 * (ewi + ewk);
892       double ei = 0.0;
893 
894       // Electrostatic solvation energy of the permanent multipoles in the
895       // GK reaction potential of the induced dipoles.
896       if (polarization != Polarization.NONE) {
897         double esymi =
898             -uxi * (dxk * gux[2] + dyk * guy[2] + dzk * guz[2])
899                 - uyi * (dxk * gux[3] + dyk * guy[3] + dzk * guz[3])
900                 - uzi * (dxk * gux[4] + dyk * guy[4] + dzk * guz[4])
901                 - uxk * (dxi * gux[2] + dyi * guy[2] + dzi * guz[2])
902                 - uyk * (dxi * gux[3] + dyi * guy[3] + dzi * guz[3])
903                 - uzk * (dxi * gux[4] + dyi * guy[4] + dzi * guz[4]);
904         double ewii =
905             ci * (dxk * gc[2] + dyk * gc[3] + dzk * gc[4])
906                 - ck * (dxi * gux[1] + dyi * guy[1] + dzi * guz[1])
907                 - dxi
908                 * (qxxk * gux[5]
909                 + qyyk * gux[8]
910                 + qzzk * gux[10]
911                 + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
912                 - dyi
913                 * (qxxk * guy[5]
914                 + qyyk * guy[8]
915                 + qzzk * guy[10]
916                 + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
917                 - dzi
918                 * (qxxk * guz[5]
919                 + qyyk * guz[8]
920                 + qzzk * guz[10]
921                 + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
922                 + dxk
923                 * (qxxi * gqxx[2]
924                 + qyyi * gqyy[2]
925                 + qzzi * gqzz[2]
926                 + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
927                 + dyk
928                 * (qxxi * gqxx[3]
929                 + qyyi * gqyy[3]
930                 + qzzi * gqzz[3]
931                 + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
932                 + dzk
933                 * (qxxi * gqxx[4]
934                 + qyyi * gqyy[4]
935                 + qzzi * gqzz[4]
936                 + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
937         double ewki =
938             ci * (dxk * gux[1] + dyk * guy[1] + dzk * guz[1])
939                 - ck * (dxi * gc[2] + dyi * gc[3] + dzi * gc[4])
940                 - dxi
941                 * (qxxk * gqxx[2]
942                 + qyyk * gqyy[2]
943                 + qzzk * gqzz[2]
944                 + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
945                 - dyi
946                 * (qxxk * gqxx[3]
947                 + qyyk * gqyy[3]
948                 + qzzk * gqzz[3]
949                 + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
950                 - dzi
951                 * (qxxk * gqxx[4]
952                 + qyyk * gqyy[4]
953                 + qzzk * gqzz[4]
954                 + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
955                 + dxk
956                 * (qxxi * gux[5]
957                 + qyyi * gux[8]
958                 + qzzi * gux[10]
959                 + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
960                 + dyk
961                 * (qxxi * guy[5]
962                 + qyyi * guy[8]
963                 + qzzi * guy[10]
964                 + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
965                 + dzk
966                 * (qxxi * guz[5]
967                 + qyyi * guz[8]
968                 + qzzi * guz[10]
969                 + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]));
970         ei = 0.5 * (esymi + 0.5 * (ewii + ewki));
971       }
972 
973       if (i == k) {
974         e *= 0.5;
975         ei *= 0.5;
976         selfEnergy.add(threadID, i, e + ei);
977       } else {
978         double half = 0.5 * (e + ei);
979         crossEnergy.add(threadID, i, half);
980         crossEnergy.add(threadID, k, half);
981       }
982 
983       gkPermanentEnergy += e;
984       gkPolarizationEnergy += ei;
985 
986       return e + ei;
987     }
988 
989     private void gradientTensors() {
990       // Born radii gradients of unweighted reaction potential tensor.
991       gc[21] = b[0][0];
992       gux[21] = xr * b[1][0];
993       guy[21] = yr * b[1][0];
994       guz[21] = zr * b[1][0];
995       gqxx[21] = xr2 * b[2][0];
996       gqyy[21] = yr2 * b[2][0];
997       gqzz[21] = zr2 * b[2][0];
998       gqxy[21] = xr * yr * b[2][0];
999       gqxz[21] = xr * zr * b[2][0];
1000       gqyz[21] = yr * zr * b[2][0];
1001 
1002       // Born gradients of the unweighted reaction potential gradient tensor
1003       gc[22] = xr * b[0][1];
1004       gc[23] = yr * b[0][1];
1005       gc[24] = zr * b[0][1];
1006       gux[22] = b[1][0] + xr2 * b[1][1];
1007       gux[23] = xr * yr * b[1][1];
1008       gux[24] = xr * zr * b[1][1];
1009       guy[22] = gux[23];
1010       guy[23] = b[1][0] + yr2 * b[1][1];
1011       guy[24] = yr * zr * b[1][1];
1012       guz[22] = gux[24];
1013       guz[23] = guy[24];
1014       guz[24] = b[1][0] + zr2 * b[1][1];
1015       gqxx[22] = xr * (2.0 * b[2][0] + xr2 * b[2][1]);
1016       gqxx[23] = yr * xr2 * b[2][1];
1017       gqxx[24] = zr * xr2 * b[2][1];
1018       gqyy[22] = xr * yr2 * b[2][1];
1019       gqyy[23] = yr * (2.0 * b[2][0] + yr2 * b[2][1]);
1020       gqyy[24] = zr * yr2 * b[2][1];
1021       gqzz[22] = xr * zr2 * b[2][1];
1022       gqzz[23] = yr * zr2 * b[2][1];
1023       gqzz[24] = zr * (2.0 * b[2][0] + zr2 * b[2][1]);
1024       gqxy[22] = yr * (b[2][0] + xr2 * b[2][1]);
1025       gqxy[23] = xr * (b[2][0] + yr2 * b[2][1]);
1026       gqxy[24] = zr * xr * yr * b[2][1];
1027       gqxz[22] = zr * (b[2][0] + xr2 * b[2][1]);
1028       gqxz[23] = gqxy[24];
1029       gqxz[24] = xr * (b[2][0] + zr2 * b[2][1]);
1030       gqyz[22] = gqxy[24];
1031       gqyz[23] = zr * (b[2][0] + yr2 * b[2][1]);
1032       gqyz[24] = yr * (b[2][0] + zr2 * b[2][1]);
1033 
1034       // Born radii derivatives of the unweighted 2nd reaction potential gradient tensor.
1035       gc[25] = b[0][1] + xr2 * b[0][2];
1036       gc[26] = xr * yr * b[0][2];
1037       gc[27] = xr * zr * b[0][2];
1038       gc[28] = b[0][1] + yr2 * b[0][2];
1039       gc[29] = yr * zr * b[0][2];
1040       gc[30] = b[0][1] + zr2 * b[0][2];
1041       gux[25] = xr * (3.0 * b[1][1] + xr2 * b[1][2]);
1042       gux[26] = yr * (b[1][1] + xr2 * b[1][2]);
1043       gux[27] = zr * (b[1][1] + xr2 * b[1][2]);
1044       gux[28] = xr * (b[1][1] + yr2 * b[1][2]);
1045       gux[29] = zr * xr * yr * b[1][2];
1046       gux[30] = xr * (b[1][1] + zr2 * b[1][2]);
1047       guy[25] = yr * (b[1][1] + xr2 * b[1][2]);
1048       guy[26] = xr * (b[1][1] + yr2 * b[1][2]);
1049       guy[27] = gux[29];
1050       guy[28] = yr * (3.0 * b[1][1] + yr2 * b[1][2]);
1051       guy[29] = zr * (b[1][1] + yr2 * b[1][2]);
1052       guy[30] = yr * (b[1][1] + zr2 * b[1][2]);
1053       guz[25] = zr * (b[1][1] + xr2 * b[1][2]);
1054       guz[26] = gux[29];
1055       guz[27] = xr * (b[1][1] + zr2 * b[1][2]);
1056       guz[28] = zr * (b[1][1] + yr2 * b[1][2]);
1057       guz[29] = yr * (b[1][1] + zr2 * b[1][2]);
1058       guz[30] = zr * (3.0 * b[1][1] + zr2 * b[1][2]);
1059       gqxx[25] = 2.0 * b[2][0] + xr2 * (5.0 * b[2][1] + xr2 * b[2][2]);
1060       gqxx[26] = yr * xr * (2.0 * b[2][1] + xr2 * b[2][2]);
1061       gqxx[27] = zr * xr * (2.0 * b[2][1] + xr2 * b[2][2]);
1062       gqxx[28] = xr2 * (b[2][1] + yr2 * b[2][2]);
1063       gqxx[29] = zr * yr * xr2 * b[2][2];
1064       gqxx[30] = xr2 * (b[2][1] + zr2 * b[2][2]);
1065       gqyy[25] = yr2 * (b[2][1] + xr2 * b[2][2]);
1066       gqyy[26] = xr * yr * (2.0 * b[2][1] + yr2 * b[2][2]);
1067       gqyy[27] = xr * zr * yr2 * b[2][2];
1068       gqyy[28] = 2.0 * b[2][0] + yr2 * (5.0 * b[2][1] + yr2 * b[2][2]);
1069       gqyy[29] = yr * zr * (2.0 * b[2][1] + yr2 * b[2][2]);
1070       gqyy[30] = yr2 * (b[2][1] + zr2 * b[2][2]);
1071       gqzz[25] = zr2 * (b[2][1] + xr2 * b[2][2]);
1072       gqzz[26] = xr * yr * zr2 * b[2][2];
1073       gqzz[27] = xr * zr * (2.0 * b[2][1] + zr2 * b[2][2]);
1074       gqzz[28] = zr2 * (b[2][1] + yr2 * b[2][2]);
1075       gqzz[29] = yr * zr * (2.0 * b[2][1] + zr2 * b[2][2]);
1076       gqzz[30] = 2.0 * b[2][0] + zr2 * (5.0 * b[2][1] + zr2 * b[2][2]);
1077       gqxy[25] = xr * yr * (3.0 * b[2][1] + xr2 * b[2][2]);
1078       gqxy[26] = b[2][0] + (xr2 + yr2) * b[2][1] + xr2 * yr2 * b[2][2];
1079       gqxy[27] = zr * yr * (b[2][1] + xr2 * b[2][2]);
1080       gqxy[28] = xr * yr * (3.0 * b[2][1] + yr2 * b[2][2]);
1081       gqxy[29] = zr * xr * (b[2][1] + yr2 * b[2][2]);
1082       gqxy[30] = xr * yr * (b[2][1] + zr2 * b[2][2]);
1083       gqxz[25] = xr * zr * (3.0 * b[2][1] + xr2 * b[2][2]);
1084       gqxz[26] = yr * zr * (b[2][1] + xr2 * b[2][2]);
1085       gqxz[27] = b[2][0] + (xr2 + zr2) * b[2][1] + xr2 * zr2 * b[2][2];
1086       gqxz[28] = xr * zr * (b[2][1] + yr2 * b[2][2]);
1087       gqxz[29] = xr * yr * (b[2][1] + zr2 * b[2][2]);
1088       gqxz[30] = xr * zr * (3.0 * b[2][1] + zr2 * b[2][2]);
1089       gqyz[25] = zr * yr * (b[2][1] + xr2 * b[2][2]);
1090       gqyz[26] = xr * zr * (b[2][1] + yr2 * b[2][2]);
1091       gqyz[27] = xr * yr * (b[2][1] + zr2 * b[2][2]);
1092       gqyz[28] = yr * zr * (3.0 * b[2][1] + yr2 * b[2][2]);
1093       gqyz[29] = b[2][0] + (yr2 + zr2) * b[2][1] + yr2 * zr2 * b[2][2];
1094       gqyz[30] = yr * zr * (3.0 * b[2][1] + zr2 * b[2][2]);
1095 
1096       // Unweighted 3rd reaction potential gradient tensor.
1097       gc[11] = xr * (3.0 * a[0][2] + xr2 * a[0][3]);
1098       gc[12] = yr * (a[0][2] + xr2 * a[0][3]);
1099       gc[13] = zr * (a[0][2] + xr2 * a[0][3]);
1100       gc[14] = xr * (a[0][2] + yr2 * a[0][3]);
1101       gc[15] = xr * yr * zr * a[0][3];
1102       gc[16] = xr * (a[0][2] + zr2 * a[0][3]);
1103       gc[17] = yr * (3.0 * a[0][2] + yr2 * a[0][3]);
1104       gc[18] = zr * (a[0][2] + yr2 * a[0][3]);
1105       gc[19] = yr * (a[0][2] + zr2 * a[0][3]);
1106       gc[20] = zr * (3.0 * a[0][2] + zr2 * a[0][3]);
1107 
1108       gux[11] = 3.0 * a[1][1] + xr2 * (6.0 * a[1][2] + xr2 * a[1][3]);
1109       gux[12] = xr * yr * (3.0 * a[1][2] + xr2 * a[1][3]);
1110       gux[13] = xr * zr * (3.0 * a[1][2] + xr2 * a[1][3]);
1111       gux[14] = a[1][1] + (xr2 + yr2) * a[1][2] + xr2 * yr2 * a[1][3];
1112       gux[15] = yr * zr * (a[1][2] + xr2 * a[1][3]);
1113       gux[16] = a[1][1] + (xr2 + zr2) * a[1][2] + xr2 * zr2 * a[1][3];
1114       gux[17] = xr * yr * (3.0 * a[1][2] + yr2 * a[1][3]);
1115       gux[18] = xr * zr * (a[1][2] + yr2 * a[1][3]);
1116       gux[19] = xr * yr * (a[1][2] + zr2 * a[1][3]);
1117       gux[20] = xr * zr * (3.0 * a[1][2] + zr2 * a[1][3]);
1118 
1119       guy[11] = gux[12];
1120       guy[12] = gux[14];
1121       guy[13] = gux[15];
1122       guy[14] = gux[17];
1123       guy[15] = gux[18];
1124       guy[16] = gux[19];
1125       guy[17] = 3.0 * a[1][1] + yr2 * (6.0 * a[1][2] + yr2 * a[1][3]);
1126       guy[18] = yr * zr * (3.0 * a[1][2] + yr2 * a[1][3]);
1127       guy[19] = a[1][1] + (yr2 + zr2) * a[1][2] + yr2 * zr2 * a[1][3];
1128       guy[20] = yr * zr * (3.0 * a[1][2] + zr2 * a[1][3]);
1129 
1130       guz[11] = gux[13];
1131       guz[12] = gux[15];
1132       guz[13] = gux[16];
1133       guz[14] = gux[18];
1134       guz[15] = gux[19];
1135       guz[16] = gux[20];
1136       guz[17] = guy[18];
1137       guz[18] = guy[19];
1138       guz[19] = guy[20];
1139       guz[20] = 3.0 * a[1][1] + zr2 * (6.0 * a[1][2] + zr2 * a[1][3]);
1140 
1141       gqxx[11] = xr * (12.0 * a[2][1] + xr2 * (9.0 * a[2][2] + xr2 * a[2][3]));
1142       gqxx[12] = yr * (2.0 * a[2][1] + xr2 * (5.0 * a[2][2] + xr2 * a[2][3]));
1143       gqxx[13] = zr * (2.0 * a[2][1] + xr2 * (5.0 * a[2][2] + xr2 * a[2][3]));
1144       gqxx[14] = xr * (2.0 * a[2][1] + yr2 * 2.0 * a[2][2] + xr2 * (a[2][2] + yr2 * a[2][3]));
1145       gqxx[15] = xr * yr * zr * (2.0 * a[2][2] + xr2 * a[2][3]);
1146       gqxx[16] = xr * (2.0 * a[2][1] + zr2 * 2.0 * a[2][2] + xr2 * (a[2][2] + zr2 * a[2][3]));
1147       gqxx[17] = yr * xr2 * (3.0 * a[2][2] + yr2 * a[2][3]);
1148       gqxx[18] = zr * xr2 * (a[2][2] + yr2 * a[2][3]);
1149       gqxx[19] = yr * xr2 * (a[2][2] + zr2 * a[2][3]);
1150       gqxx[20] = zr * xr2 * (3.0 * a[2][2] + zr2 * a[2][3]);
1151 
1152       gqxy[11] = yr * (3.0 * a[2][1] + xr2 * (6.0 * a[2][2] + xr2 * a[2][3]));
1153       gqxy[12] = xr * (3.0 * (a[2][1] + yr2 * a[2][2]) + xr2 * (a[2][2] + yr2 * a[2][3]));
1154       gqxy[13] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
1155       gqxy[14] = yr * (3.0 * (a[2][1] + xr2 * a[2][2]) + yr2 * (a[2][2] + xr2 * a[2][3]));
1156       gqxy[15] = zr * (a[2][1] + (yr2 + xr2) * a[2][2] + yr2 * xr2 * a[2][3]);
1157       gqxy[16] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + xr2 * zr2 * a[2][3]);
1158       gqxy[17] = xr * (3.0 * (a[2][1] + yr2 * a[2][2]) + yr2 * (3.0 * a[2][2] + yr2 * a[2][3]));
1159       gqxy[18] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
1160       gqxy[19] = xr * (a[2][1] + (yr2 + zr2) * a[2][2] + yr2 * zr2 * a[2][3]);
1161       gqxy[20] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
1162 
1163       gqxz[11] = zr * (3.0 * a[2][1] + xr2 * (6.0 * a[2][2] + xr2 * a[2][3]));
1164       gqxz[12] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
1165       gqxz[13] = xr * (3.0 * (a[2][1] + zr2 * a[2][2]) + xr2 * (a[2][2] + zr2 * a[2][3]));
1166       gqxz[14] = zr * (a[2][1] + (xr2 + yr2) * a[2][2] + xr2 * yr2 * a[2][3]);
1167       gqxz[15] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + zr2 * xr2 * a[2][3]);
1168       gqxz[16] = zr * (3.0 * (a[2][1] + xr2 * a[2][2]) + zr2 * (a[2][2] + xr2 * a[2][3]));
1169       gqxz[17] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
1170       gqxz[18] = xr * (a[2][1] + (zr2 + yr2) * a[2][2] + zr2 * yr2 * a[2][3]);
1171       gqxz[19] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
1172       gqxz[20] = xr * (3.0 * a[2][1] + zr2 * (6.0 * a[2][2] + zr2 * a[2][3]));
1173 
1174       gqyy[11] = xr * yr2 * (3.0 * a[2][2] + xr2 * a[2][3]);
1175       gqyy[12] = yr * (2.0 * a[2][1] + xr2 * 2.0 * a[2][2] + yr2 * (a[2][2] + xr2 * a[2][3]));
1176       gqyy[13] = zr * yr2 * (a[2][2] + xr2 * a[2][3]);
1177       gqyy[14] = xr * (2.0 * a[2][1] + yr2 * (5.0 * a[2][2] + yr2 * a[2][3]));
1178       gqyy[15] = xr * yr * zr * (2.0 * a[2][2] + yr2 * a[2][3]);
1179       gqyy[16] = xr * yr2 * (a[2][2] + zr2 * a[2][3]);
1180       gqyy[17] = yr * (12.0 * a[2][1] + yr2 * (9.0 * a[2][2] + yr2 * a[2][3]));
1181       gqyy[18] = zr * (2.0 * a[2][1] + yr2 * (5.0 * a[2][2] + yr2 * a[2][3]));
1182       gqyy[19] = yr * (2.0 * a[2][1] + zr2 * 2.0 * a[2][2] + yr2 * (a[2][2] + zr2 * a[2][3]));
1183       gqyy[20] = zr * yr2 * (3.0 * a[2][2] + zr2 * a[2][3]);
1184 
1185       gqyz[11] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
1186       gqyz[12] = zr * (a[2][1] + (xr2 + yr2) * a[2][2] + xr2 * yr2 * a[2][3]);
1187       gqyz[13] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + xr2 * zr2 * a[2][3]);
1188       gqyz[14] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
1189       gqyz[15] = xr * (a[2][1] + (yr2 + zr2) * a[2][2] + yr2 * zr2 * a[2][3]);
1190       gqyz[16] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
1191       gqyz[17] = zr * (3.0 * a[2][1] + yr2 * (6.0 * a[2][2] + yr2 * a[2][3]));
1192       gqyz[18] = yr * (3.0 * (a[2][1] + zr2 * a[2][2]) + yr2 * (a[2][2] + zr2 * a[2][3]));
1193       gqyz[19] = zr * (3.0 * (a[2][1] + yr2 * a[2][2]) + zr2 * (a[2][2] + yr2 * a[2][3]));
1194       gqyz[20] = yr * (3.0 * a[2][1] + zr2 * (6.0 * a[2][2] + zr2 * a[2][3]));
1195 
1196       gqzz[11] = xr * zr2 * (3.0 * a[2][2] + xr2 * a[2][3]);
1197       gqzz[12] = yr * (zr2 * a[2][2] + xr2 * (zr2 * a[2][3]));
1198       gqzz[13] = zr * (2.0 * a[2][1] + xr2 * 2.0 * a[2][2] + zr2 * (a[2][2] + xr2 * a[2][3]));
1199       gqzz[14] = xr * zr2 * (a[2][2] + yr2 * a[2][3]);
1200       gqzz[15] = xr * yr * zr * (2.0 * a[2][2] + zr2 * a[2][3]);
1201       gqzz[16] = xr * (2.0 * a[2][1] + zr2 * (5.0 * a[2][2] + zr2 * a[2][3]));
1202       gqzz[17] = yr * zr2 * (3.0 * a[2][2] + yr2 * a[2][3]);
1203       gqzz[18] = zr * (2.0 * a[2][1] + yr2 * 2.0 * a[2][2] + zr2 * (a[2][2] + yr2 * a[2][3]));
1204       gqzz[19] = yr * (2.0 * a[2][1] + zr2 * (5.0 * a[2][2] + zr2 * a[2][3]));
1205       gqzz[20] = zr * (12.0 * a[2][1] + zr2 * (9.0 * a[2][2] + zr2 * a[2][3]));
1206     }
1207 
1208     private void permanentEnergyGradient(int i, int k) {
1209       final double desymdr =
1210           ci * ck * gc[21]
1211               - (uxi * (uxk * gux[22] + uyk * guy[22] + uzk * guz[22])
1212               + uyi * (uxk * gux[23] + uyk * guy[23] + uzk * guz[23])
1213               + uzi * (uxk * gux[24] + uyk * guy[24] + uzk * guz[24]));
1214       final double dewidr =
1215           ci * (uxk * gc[22] + uyk * gc[23] + uzk * gc[24])
1216               - ck * (uxi * gux[21] + uyi * guy[21] + uzi * guz[21])
1217               + ci
1218               * (qxxk * gc[25]
1219               + qyyk * gc[28]
1220               + qzzk * gc[30]
1221               + 2.0 * (qxyk * gc[26] + qxzk * gc[27] + qyzk * gc[29]))
1222               + ck
1223               * (qxxi * gqxx[21]
1224               + qyyi * gqyy[21]
1225               + qzzi * gqzz[21]
1226               + 2.0 * (qxyi * gqxy[21] + qxzi * gqxz[21] + qyzi * gqyz[21]))
1227               - uxi
1228               * (qxxk * gux[25]
1229               + qyyk * gux[28]
1230               + qzzk * gux[30]
1231               + 2.0 * (qxyk * gux[26] + qxzk * gux[27] + qyzk * gux[29]))
1232               - uyi
1233               * (qxxk * guy[25]
1234               + qyyk * guy[28]
1235               + qzzk * guy[30]
1236               + 2.0 * (qxyk * guy[26] + qxzk * guy[27] + qyzk * guy[29]))
1237               - uzi
1238               * (qxxk * guz[25]
1239               + qyyk * guz[28]
1240               + qzzk * guz[30]
1241               + 2.0 * (qxyk * guz[26] + qxzk * guz[27] + qyzk * guz[29]))
1242               + uxk
1243               * (qxxi * gqxx[22]
1244               + qyyi * gqyy[22]
1245               + qzzi * gqzz[22]
1246               + 2.0 * (qxyi * gqxy[22] + qxzi * gqxz[22] + qyzi * gqyz[22]))
1247               + uyk
1248               * (qxxi * gqxx[23]
1249               + qyyi * gqyy[23]
1250               + qzzi * gqzz[23]
1251               + 2.0 * (qxyi * gqxy[23] + qxzi * gqxz[23] + qyzi * gqyz[23]))
1252               + uzk
1253               * (qxxi * gqxx[24]
1254               + qyyi * gqyy[24]
1255               + qzzi * gqzz[24]
1256               + 2.0 * (qxyi * gqxy[24] + qxzi * gqxz[24] + qyzi * gqyz[24]))
1257               + qxxi
1258               * (qxxk * gqxx[25]
1259               + qyyk * gqxx[28]
1260               + qzzk * gqxx[30]
1261               + 2.0 * (qxyk * gqxx[26] + qxzk * gqxx[27] + qyzk * gqxx[29]))
1262               + qyyi
1263               * (qxxk * gqyy[25]
1264               + qyyk * gqyy[28]
1265               + qzzk * gqyy[30]
1266               + 2.0 * (qxyk * gqyy[26] + qxzk * gqyy[27] + qyzk * gqyy[29]))
1267               + qzzi
1268               * (qxxk * gqzz[25]
1269               + qyyk * gqzz[28]
1270               + qzzk * gqzz[30]
1271               + 2.0 * (qxyk * gqzz[26] + qxzk * gqzz[27] + qyzk * gqzz[29]))
1272               + 2.0
1273               * (qxyi
1274               * (qxxk * gqxy[25]
1275               + qyyk * gqxy[28]
1276               + qzzk * gqxy[30]
1277               + 2.0 * (qxyk * gqxy[26] + qxzk * gqxy[27] + qyzk * gqxy[29]))
1278               + qxzi
1279               * (qxxk * gqxz[25]
1280               + qyyk * gqxz[28]
1281               + qzzk * gqxz[30]
1282               + 2.0 * (qxyk * gqxz[26] + qxzk * gqxz[27] + qyzk * gqxz[29]))
1283               + qyzi
1284               * (qxxk * gqyz[25]
1285               + qyyk * gqyz[28]
1286               + qzzk * gqyz[30]
1287               + 2.0 * (qxyk * gqyz[26] + qxzk * gqyz[27] + qyzk * gqyz[29])));
1288       final double dewkdr =
1289           ci * (uxk * gux[21] + uyk * guy[21] + uzk * guz[21])
1290               - ck * (uxi * gc[22] + uyi * gc[23] + uzi * gc[24])
1291               + ci
1292               * (qxxk * gqxx[21]
1293               + qyyk * gqyy[21]
1294               + qzzk * gqzz[21]
1295               + 2.0 * (qxyk * gqxy[21] + qxzk * gqxz[21] + qyzk * gqyz[21]))
1296               + ck
1297               * (qxxi * gc[25]
1298               + qyyi * gc[28]
1299               + qzzi * gc[30]
1300               + 2.0 * (qxyi * gc[26] + qxzi * gc[27] + qyzi * gc[29]))
1301               - uxi
1302               * (qxxk * gqxx[22]
1303               + qyyk * gqyy[22]
1304               + qzzk * gqzz[22]
1305               + 2.0 * (qxyk * gqxy[22] + qxzk * gqxz[22] + qyzk * gqyz[22]))
1306               - uyi
1307               * (qxxk * gqxx[23]
1308               + qyyk * gqyy[23]
1309               + qzzk * gqzz[23]
1310               + 2.0 * (qxyk * gqxy[23] + qxzk * gqxz[23] + qyzk * gqyz[23]))
1311               - uzi
1312               * (qxxk * gqxx[24]
1313               + qyyk * gqyy[24]
1314               + qzzk * gqzz[24]
1315               + 2.0 * (qxyk * gqxy[24] + qxzk * gqxz[24] + qyzk * gqyz[24]))
1316               + uxk
1317               * (qxxi * gux[25]
1318               + qyyi * gux[28]
1319               + qzzi * gux[30]
1320               + 2.0 * (qxyi * gux[26] + qxzi * gux[27] + qyzi * gux[29]))
1321               + uyk
1322               * (qxxi * guy[25]
1323               + qyyi * guy[28]
1324               + qzzi * guy[30]
1325               + 2.0 * (qxyi * guy[26] + qxzi * guy[27] + qyzi * guy[29]))
1326               + uzk
1327               * (qxxi * guz[25]
1328               + qyyi * guz[28]
1329               + qzzi * guz[30]
1330               + 2.0 * (qxyi * guz[26] + qxzi * guz[27] + qyzi * guz[29]))
1331               + qxxi
1332               * (qxxk * gqxx[25]
1333               + qyyk * gqyy[25]
1334               + qzzk * gqzz[25]
1335               + 2.0 * (qxyk * gqxy[25] + qxzk * gqxz[25] + qyzk * gqyz[25]))
1336               + qyyi
1337               * (qxxk * gqxx[28]
1338               + qyyk * gqyy[28]
1339               + qzzk * gqzz[28]
1340               + 2.0 * (qxyk * gqxy[28] + qxzk * gqxz[28] + qyzk * gqyz[28]))
1341               + qzzi
1342               * (qxxk * gqxx[30]
1343               + qyyk * gqyy[30]
1344               + qzzk * gqzz[30]
1345               + 2.0 * (qxyk * gqxy[30] + qxzk * gqxz[30] + qyzk * gqyz[30]))
1346               + 2.0
1347               * (qxyi
1348               * (qxxk * gqxx[26]
1349               + qyyk * gqyy[26]
1350               + qzzk * gqzz[26]
1351               + 2.0 * (qxyk * gqxy[26] + qxzk * gqxz[26] + qyzk * gqyz[26]))
1352               + qxzi
1353               * (qxxk * gqxx[27]
1354               + qyyk * gqyy[27]
1355               + qzzk * gqzz[27]
1356               + 2.0 * (qxyk * gqxy[27] + qxzk * gqxz[27] + qyzk * gqyz[27]))
1357               + qyzi
1358               * (qxxk * gqxx[29]
1359               + qyyk * gqyy[29]
1360               + qzzk * gqzz[29]
1361               + 2.0 * (qxyk * gqxy[29] + qxzk * gqxz[29] + qyzk * gqyz[29])));
1362       final double dsumdr = desymdr + 0.5 * (dewidr + dewkdr);
1363       final double drbi = rbk * dsumdr;
1364       final double drbk = rbi * dsumdr;
1365 
1366       double selfScale = 1.0;
1367       if (i == k) {
1368         if (iSymm == 0) {
1369           sharedBornGrad.add(threadID, i, drbi);
1370           return;
1371         } else {
1372           selfScale = 0.5;
1373         }
1374       }
1375 
1376       // Increment the gradients and Born chain rule term.
1377       final double dedx = selfScale * dEdX();
1378       final double dedy = selfScale * dEdY();
1379       final double dedz = selfScale * dEdZ();
1380       dedxi -= dedx;
1381       dedyi -= dedy;
1382       dedzi -= dedz;
1383       dborni += selfScale * drbi;
1384 
1385       final double dedxk = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
1386       final double dedyk = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
1387       final double dedzk = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
1388       grad.add(threadID, k, dedxk, dedyk, dedzk);
1389       sharedBornGrad.add(threadID, k, selfScale * drbk);
1390       permanentEnergyTorque(i, k);
1391     }
1392 
1393     private double dEdZ() {
1394       final double desymdz =
1395           ci * ck * gc[4]
1396               - (uxi * (uxk * gux[7] + uyk * guy[7] + uzk * guz[7])
1397               + uyi * (uxk * gux[9] + uyk * guy[9] + uzk * guz[9])
1398               + uzi * (uxk * gux[10] + uyk * guy[10] + uzk * guz[10]));
1399       final double dewidz =
1400           ci * (uxk * gc[7] + uyk * gc[9] + uzk * gc[10])
1401               - ck * (uxi * gux[4] + uyi * guy[4] + uzi * guz[4])
1402               + ci
1403               * (qxxk * gc[13]
1404               + qyyk * gc[18]
1405               + qzzk * gc[20]
1406               + 2.0 * (qxyk * gc[15] + qxzk * gc[16] + qyzk * gc[19]))
1407               + ck
1408               * (qxxi * gqxx[4]
1409               + qyyi * gqyy[4]
1410               + qzzi * gqzz[4]
1411               + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]))
1412               - uxi
1413               * (qxxk * gux[13]
1414               + qyyk * gux[18]
1415               + qzzk * gux[20]
1416               + 2.0 * (qxyk * gux[15] + qxzk * gux[16] + qyzk * gux[19]))
1417               - uyi
1418               * (qxxk * guy[13]
1419               + qyyk * guy[18]
1420               + qzzk * guy[20]
1421               + 2.0 * (qxyk * guy[15] + qxzk * guy[16] + qyzk * guy[19]))
1422               - uzi
1423               * (qxxk * guz[13]
1424               + qyyk * guz[18]
1425               + qzzk * guz[20]
1426               + 2.0 * (qxyk * guz[15] + qxzk * guz[16] + qyzk * guz[19]))
1427               + uxk
1428               * (qxxi * gqxx[7]
1429               + qyyi * gqyy[7]
1430               + qzzi * gqzz[7]
1431               + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
1432               + uyk
1433               * (qxxi * gqxx[9]
1434               + qyyi * gqyy[9]
1435               + qzzi * gqzz[9]
1436               + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
1437               + uzk
1438               * (qxxi * gqxx[10]
1439               + qyyi * gqyy[10]
1440               + qzzi * gqzz[10]
1441               + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]))
1442               + qxxi
1443               * (qxxk * gqxx[13]
1444               + qyyk * gqxx[18]
1445               + qzzk * gqxx[20]
1446               + 2.0 * (qxyk * gqxx[15] + qxzk * gqxx[16] + qyzk * gqxx[19]))
1447               + qyyi
1448               * (qxxk * gqyy[13]
1449               + qyyk * gqyy[18]
1450               + qzzk * gqyy[20]
1451               + 2.0 * (qxyk * gqyy[15] + qxzk * gqyy[16] + qyzk * gqyy[19]))
1452               + qzzi
1453               * (qxxk * gqzz[13]
1454               + qyyk * gqzz[18]
1455               + qzzk * gqzz[20]
1456               + 2.0 * (qxyk * gqzz[15] + qxzk * gqzz[16] + qyzk * gqzz[19]))
1457               + 2.0
1458               * (qxyi
1459               * (qxxk * gqxy[13]
1460               + qyyk * gqxy[18]
1461               + qzzk * gqxy[20]
1462               + 2.0 * (qxyk * gqxy[15] + qxzk * gqxy[16] + qyzk * gqxy[19]))
1463               + qxzi
1464               * (qxxk * gqxz[13]
1465               + qyyk * gqxz[18]
1466               + qzzk * gqxz[20]
1467               + 2.0 * (qxyk * gqxz[15] + qxzk * gqxz[16] + qyzk * gqxz[19]))
1468               + qyzi
1469               * (qxxk * gqyz[13]
1470               + qyyk * gqyz[18]
1471               + qzzk * gqyz[20]
1472               + 2.0 * (qxyk * gqyz[15] + qxzk * gqyz[16] + qyzk * gqyz[19])));
1473       final double dewkdz =
1474           ci * (uxk * gux[4] + uyk * guy[4] + uzk * guz[4])
1475               - ck * (uxi * gc[7] + uyi * gc[9] + uzi * gc[10])
1476               + ci
1477               * (qxxk * gqxx[4]
1478               + qyyk * gqyy[4]
1479               + qzzk * gqzz[4]
1480               + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
1481               + ck
1482               * (qxxi * gc[13]
1483               + qyyi * gc[18]
1484               + qzzi * gc[20]
1485               + 2.0 * (qxyi * gc[15] + qxzi * gc[16] + qyzi * gc[19]))
1486               - uxi
1487               * (qxxk * gqxx[7]
1488               + qyyk * gqyy[7]
1489               + qzzk * gqzz[7]
1490               + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
1491               - uyi
1492               * (qxxk * gqxx[9]
1493               + qyyk * gqyy[9]
1494               + qzzk * gqzz[9]
1495               + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
1496               - uzi
1497               * (qxxk * gqxx[10]
1498               + qyyk * gqyy[10]
1499               + qzzk * gqzz[10]
1500               + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
1501               + uxk
1502               * (qxxi * gux[13]
1503               + qyyi * gux[18]
1504               + qzzi * gux[20]
1505               + 2.0 * (qxyi * gux[15] + qxzi * gux[16] + qyzi * gux[19]))
1506               + uyk
1507               * (qxxi * guy[13]
1508               + qyyi * guy[18]
1509               + qzzi * guy[20]
1510               + 2.0 * (qxyi * guy[15] + qxzi * guy[16] + qyzi * guy[19]))
1511               + uzk
1512               * (qxxi * guz[13]
1513               + qyyi * guz[18]
1514               + qzzi * guz[20]
1515               + 2.0 * (qxyi * guz[15] + qxzi * guz[16] + qyzi * guz[19]))
1516               + qxxi
1517               * (qxxk * gqxx[13]
1518               + qyyk * gqyy[13]
1519               + qzzk * gqzz[13]
1520               + 2.0 * (qxyk * gqxy[13] + qxzk * gqxz[13] + qyzk * gqyz[13]))
1521               + qyyi
1522               * (qxxk * gqxx[18]
1523               + qyyk * gqyy[18]
1524               + qzzk * gqzz[18]
1525               + 2.0 * (qxyk * gqxy[18] + qxzk * gqxz[18] + qyzk * gqyz[18]))
1526               + qzzi
1527               * (qxxk * gqxx[20]
1528               + qyyk * gqyy[20]
1529               + qzzk * gqzz[20]
1530               + 2.0 * (qxyk * gqxy[20] + qxzk * gqxz[20] + qyzk * gqyz[20]))
1531               + 2.0
1532               * (qxyi
1533               * (qxxk * gqxx[15]
1534               + qyyk * gqyy[15]
1535               + qzzk * gqzz[15]
1536               + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15]))
1537               + qxzi
1538               * (qxxk * gqxx[16]
1539               + qyyk * gqyy[16]
1540               + qzzk * gqzz[16]
1541               + 2.0 * (qxyk * gqxy[16] + qxzk * gqxz[16] + qyzk * gqyz[16]))
1542               + qyzi
1543               * (qxxk * gqxx[19]
1544               + qyyk * gqyy[19]
1545               + qzzk * gqzz[19]
1546               + 2.0 * (qxyk * gqxy[19] + qxzk * gqxz[19] + qyzk * gqyz[19])));
1547       return desymdz + 0.5 * (dewidz + dewkdz);
1548     }
1549 
1550     private double dEdY() {
1551       final double desymdy =
1552           ci * ck * gc[3]
1553               - (uxi * (uxk * gux[6] + uyk * guy[6] + uzk * guz[6])
1554               + uyi * (uxk * gux[8] + uyk * guy[8] + uzk * guz[8])
1555               + uzi * (uxk * gux[9] + uyk * guy[9] + uzk * guz[9]));
1556       final double dewidy =
1557           ci * (uxk * gc[6] + uyk * gc[8] + uzk * gc[9])
1558               - ck * (uxi * gux[3] + uyi * guy[3] + uzi * guz[3])
1559               + ci
1560               * (qxxk * gc[12]
1561               + qyyk * gc[17]
1562               + qzzk * gc[19]
1563               + 2.0 * (qxyk * gc[14] + qxzk * gc[15] + qyzk * gc[18]))
1564               + ck
1565               * (qxxi * gqxx[3]
1566               + qyyi * gqyy[3]
1567               + qzzi * gqzz[3]
1568               + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
1569               - uxi
1570               * (qxxk * gux[12]
1571               + qyyk * gux[17]
1572               + qzzk * gux[19]
1573               + 2.0 * (qxyk * gux[14] + qxzk * gux[15] + qyzk * gux[18]))
1574               - uyi
1575               * (qxxk * guy[12]
1576               + qyyk * guy[17]
1577               + qzzk * guy[19]
1578               + 2.0 * (qxyk * guy[14] + qxzk * guy[15] + qyzk * guy[18]))
1579               - uzi
1580               * (qxxk * guz[12]
1581               + qyyk * guz[17]
1582               + qzzk * guz[19]
1583               + 2.0 * (qxyk * guz[14] + qxzk * guz[15] + qyzk * guz[18]))
1584               + uxk
1585               * (qxxi * gqxx[6]
1586               + qyyi * gqyy[6]
1587               + qzzi * gqzz[6]
1588               + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
1589               + uyk
1590               * (qxxi * gqxx[8]
1591               + qyyi * gqyy[8]
1592               + qzzi * gqzz[8]
1593               + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]))
1594               + uzk
1595               * (qxxi * gqxx[9]
1596               + qyyi * gqyy[9]
1597               + qzzi * gqzz[9]
1598               + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
1599               + qxxi
1600               * (qxxk * gqxx[12]
1601               + qyyk * gqxx[17]
1602               + qzzk * gqxx[19]
1603               + 2.0 * (qxyk * gqxx[14] + qxzk * gqxx[15] + qyzk * gqxx[18]))
1604               + qyyi
1605               * (qxxk * gqyy[12]
1606               + qyyk * gqyy[17]
1607               + qzzk * gqyy[19]
1608               + 2.0 * (qxyk * gqyy[14] + qxzk * gqyy[15] + qyzk * gqyy[18]))
1609               + qzzi
1610               * (qxxk * gqzz[12]
1611               + qyyk * gqzz[17]
1612               + qzzk * gqzz[19]
1613               + 2.0 * (qxyk * gqzz[14] + qxzk * gqzz[15] + qyzk * gqzz[18]))
1614               + 2.0
1615               * (qxyi
1616               * (qxxk * gqxy[12]
1617               + qyyk * gqxy[17]
1618               + qzzk * gqxy[19]
1619               + 2.0 * (qxyk * gqxy[14] + qxzk * gqxy[15] + qyzk * gqxy[18]))
1620               + qxzi
1621               * (qxxk * gqxz[12]
1622               + qyyk * gqxz[17]
1623               + qzzk * gqxz[19]
1624               + 2.0 * (qxyk * gqxz[14] + qxzk * gqxz[15] + qyzk * gqxz[18]))
1625               + qyzi
1626               * (qxxk * gqyz[12]
1627               + qyyk * gqyz[17]
1628               + qzzk * gqyz[19]
1629               + 2.0 * (qxyk * gqyz[14] + qxzk * gqyz[15] + qyzk * gqyz[18])));
1630       final double dewkdy =
1631           ci * (uxk * gux[3] + uyk * guy[3] + uzk * guz[3])
1632               - ck * (uxi * gc[6] + uyi * gc[8] + uzi * gc[9])
1633               + ci
1634               * (qxxk * gqxx[3]
1635               + qyyk * gqyy[3]
1636               + qzzk * gqzz[3]
1637               + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
1638               + ck
1639               * (qxxi * gc[12]
1640               + qyyi * gc[17]
1641               + qzzi * gc[19]
1642               + 2.0 * (qxyi * gc[14] + qxzi * gc[15] + qyzi * gc[18]))
1643               - uxi
1644               * (qxxk * gqxx[6]
1645               + qyyk * gqyy[6]
1646               + qzzk * gqzz[6]
1647               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
1648               - uyi
1649               * (qxxk * gqxx[8]
1650               + qyyk * gqyy[8]
1651               + qzzk * gqzz[8]
1652               + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
1653               - uzi
1654               * (qxxk * gqxx[9]
1655               + qyyk * gqyy[9]
1656               + qzzk * gqzz[9]
1657               + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
1658               + uxk
1659               * (qxxi * gux[12]
1660               + qyyi * gux[17]
1661               + qzzi * gux[19]
1662               + 2.0 * (qxyi * gux[14] + qxzi * gux[15] + qyzi * gux[18]))
1663               + uyk
1664               * (qxxi * guy[12]
1665               + qyyi * guy[17]
1666               + qzzi * guy[19]
1667               + 2.0 * (qxyi * guy[14] + qxzi * guy[15] + qyzi * guy[18]))
1668               + uzk
1669               * (qxxi * guz[12]
1670               + qyyi * guz[17]
1671               + qzzi * guz[19]
1672               + 2.0 * (qxyi * guz[14] + qxzi * guz[15] + qyzi * guz[18]))
1673               + qxxi
1674               * (qxxk * gqxx[12]
1675               + qyyk * gqyy[12]
1676               + qzzk * gqzz[12]
1677               + 2.0 * (qxyk * gqxy[12] + qxzk * gqxz[12] + qyzk * gqyz[12]))
1678               + qyyi
1679               * (qxxk * gqxx[17]
1680               + qyyk * gqyy[17]
1681               + qzzk * gqzz[17]
1682               + 2.0 * (qxyk * gqxy[17] + qxzk * gqxz[17] + qyzk * gqyz[17]))
1683               + qzzi
1684               * (qxxk * gqxx[19]
1685               + qyyk * gqyy[19]
1686               + qzzk * gqzz[19]
1687               + 2.0 * (qxyk * gqxy[19] + qxzk * gqxz[19] + qyzk * gqyz[19]))
1688               + 2.0
1689               * (qxyi
1690               * (qxxk * gqxx[14]
1691               + qyyk * gqyy[14]
1692               + qzzk * gqzz[14]
1693               + 2.0 * (qxyk * gqxy[14] + qxzk * gqxz[14] + qyzk * gqyz[14]))
1694               + qxzi
1695               * (qxxk * gqxx[15]
1696               + qyyk * gqyy[15]
1697               + qzzk * gqzz[15]
1698               + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15]))
1699               + qyzi
1700               * (qxxk * gqxx[18]
1701               + qyyk * gqyy[18]
1702               + qzzk * gqzz[18]
1703               + 2.0 * (qxyk * gqxy[18] + qxzk * gqxz[18] + qyzk * gqyz[18])));
1704 
1705       return desymdy + 0.5 * (dewidy + dewkdy);
1706     }
1707 
1708     private double dEdX() {
1709       final double desymdx =
1710           ci * ck * gc[2]
1711               - (uxi * (uxk * gux[5] + uyk * guy[5] + uzk * guz[5])
1712               + uyi * (uxk * gux[6] + uyk * guy[6] + uzk * guz[6])
1713               + uzi * (uxk * gux[7] + uyk * guy[7] + uzk * guz[7]));
1714       final double dewidx =
1715           ci * (uxk * gc[5] + uyk * gc[6] + uzk * gc[7])
1716               - ck * (uxi * gux[2] + uyi * guy[2] + uzi * guz[2])
1717               + ci
1718               * (qxxk * gc[11]
1719               + qyyk * gc[14]
1720               + qzzk * gc[16]
1721               + 2.0 * (qxyk * gc[12] + qxzk * gc[13] + qyzk * gc[15]))
1722               + ck
1723               * (qxxi * gqxx[2]
1724               + qyyi * gqyy[2]
1725               + qzzi * gqzz[2]
1726               + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
1727               - uxi
1728               * (qxxk * gux[11]
1729               + qyyk * gux[14]
1730               + qzzk * gux[16]
1731               + 2.0 * (qxyk * gux[12] + qxzk * gux[13] + qyzk * gux[15]))
1732               - uyi
1733               * (qxxk * guy[11]
1734               + qyyk * guy[14]
1735               + qzzk * guy[16]
1736               + 2.0 * (qxyk * guy[12] + qxzk * guy[13] + qyzk * guy[15]))
1737               - uzi
1738               * (qxxk * guz[11]
1739               + qyyk * guz[14]
1740               + qzzk * guz[16]
1741               + 2.0 * (qxyk * guz[12] + qxzk * guz[13] + qyzk * guz[15]))
1742               + uxk
1743               * (qxxi * gqxx[5]
1744               + qyyi * gqyy[5]
1745               + qzzi * gqzz[5]
1746               + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]))
1747               + uyk
1748               * (qxxi * gqxx[6]
1749               + qyyi * gqyy[6]
1750               + qzzi * gqzz[6]
1751               + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
1752               + uzk
1753               * (qxxi * gqxx[7]
1754               + qyyi * gqyy[7]
1755               + qzzi * gqzz[7]
1756               + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
1757               + qxxi
1758               * (qxxk * gqxx[11]
1759               + qyyk * gqxx[14]
1760               + qzzk * gqxx[16]
1761               + 2.0 * (qxyk * gqxx[12] + qxzk * gqxx[13] + qyzk * gqxx[15]))
1762               + qyyi
1763               * (qxxk * gqyy[11]
1764               + qyyk * gqyy[14]
1765               + qzzk * gqyy[16]
1766               + 2.0 * (qxyk * gqyy[12] + qxzk * gqyy[13] + qyzk * gqyy[15]))
1767               + qzzi
1768               * (qxxk * gqzz[11]
1769               + qyyk * gqzz[14]
1770               + qzzk * gqzz[16]
1771               + 2.0 * (qxyk * gqzz[12] + qxzk * gqzz[13] + qyzk * gqzz[15]))
1772               + 2.0
1773               * (qxyi
1774               * (qxxk * gqxy[11]
1775               + qyyk * gqxy[14]
1776               + qzzk * gqxy[16]
1777               + 2.0 * (qxyk * gqxy[12] + qxzk * gqxy[13] + qyzk * gqxy[15]))
1778               + qxzi
1779               * (qxxk * gqxz[11]
1780               + qyyk * gqxz[14]
1781               + qzzk * gqxz[16]
1782               + 2.0 * (qxyk * gqxz[12] + qxzk * gqxz[13] + qyzk * gqxz[15]))
1783               + qyzi
1784               * (qxxk * gqyz[11]
1785               + qyyk * gqyz[14]
1786               + qzzk * gqyz[16]
1787               + 2.0 * (qxyk * gqyz[12] + qxzk * gqyz[13] + qyzk * gqyz[15])));
1788       final double dewkdx =
1789           ci * (uxk * gux[2] + uyk * guy[2] + uzk * guz[2])
1790               - ck * (uxi * gc[5] + uyi * gc[6] + uzi * gc[7])
1791               + ci
1792               * (qxxk * gqxx[2]
1793               + qyyk * gqyy[2]
1794               + qzzk * gqzz[2]
1795               + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
1796               + ck
1797               * (qxxi * gc[11]
1798               + qyyi * gc[14]
1799               + qzzi * gc[16]
1800               + 2.0 * (qxyi * gc[12] + qxzi * gc[13] + qyzi * gc[15]))
1801               - uxi
1802               * (qxxk * gqxx[5]
1803               + qyyk * gqyy[5]
1804               + qzzk * gqzz[5]
1805               + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
1806               - uyi
1807               * (qxxk * gqxx[6]
1808               + qyyk * gqyy[6]
1809               + qzzk * gqzz[6]
1810               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
1811               - uzi
1812               * (qxxk * gqxx[7]
1813               + qyyk * gqyy[7]
1814               + qzzk * gqzz[7]
1815               + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
1816               + uxk
1817               * (qxxi * gux[11]
1818               + qyyi * gux[14]
1819               + qzzi * gux[16]
1820               + 2.0 * (qxyi * gux[12] + qxzi * gux[13] + qyzi * gux[15]))
1821               + uyk
1822               * (qxxi * guy[11]
1823               + qyyi * guy[14]
1824               + qzzi * guy[16]
1825               + 2.0 * (qxyi * guy[12] + qxzi * guy[13] + qyzi * guy[15]))
1826               + uzk
1827               * (qxxi * guz[11]
1828               + qyyi * guz[14]
1829               + qzzi * guz[16]
1830               + 2.0 * (qxyi * guz[12] + qxzi * guz[13] + qyzi * guz[15]))
1831               + qxxi
1832               * (qxxk * gqxx[11]
1833               + qyyk * gqyy[11]
1834               + qzzk * gqzz[11]
1835               + 2.0 * (qxyk * gqxy[11] + qxzk * gqxz[11] + qyzk * gqyz[11]))
1836               + qyyi
1837               * (qxxk * gqxx[14]
1838               + qyyk * gqyy[14]
1839               + qzzk * gqzz[14]
1840               + 2.0 * (qxyk * gqxy[14] + qxzk * gqxz[14] + qyzk * gqyz[14]))
1841               + qzzi
1842               * (qxxk * gqxx[16]
1843               + qyyk * gqyy[16]
1844               + qzzk * gqzz[16]
1845               + 2.0 * (qxyk * gqxy[16] + qxzk * gqxz[16] + qyzk * gqyz[16]))
1846               + 2.0
1847               * (qxyi
1848               * (qxxk * gqxx[12]
1849               + qyyk * gqyy[12]
1850               + qzzk * gqzz[12]
1851               + 2.0 * (qxyk * gqxy[12] + qxzk * gqxz[12] + qyzk * gqyz[12]))
1852               + qxzi
1853               * (qxxk * gqxx[13]
1854               + qyyk * gqyy[13]
1855               + qzzk * gqzz[13]
1856               + 2.0 * (qxyk * gqxy[13] + qxzk * gqxz[13] + qyzk * gqyz[13]))
1857               + qyzi
1858               * (qxxk * gqxx[15]
1859               + qyyk * gqyy[15]
1860               + qzzk * gqzz[15]
1861               + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15])));
1862 
1863       return desymdx + 0.5 * (dewidx + dewkdx);
1864     }
1865 
1866     private void permanentEnergyTorque(int i, int k) {
1867       // Torque on permanent dipoles due to permanent reaction field.
1868       final double ix =
1869           uxk * gux[2]
1870               + uyk * gux[3]
1871               + uzk * gux[4]
1872               + 0.5
1873               * (ck * gux[1]
1874               + qxxk * gux[5]
1875               + qyyk * gux[8]
1876               + qzzk * gux[10]
1877               + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9])
1878               + ck * gc[2]
1879               + qxxk * gqxx[2]
1880               + qyyk * gqyy[2]
1881               + qzzk * gqzz[2]
1882               + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]));
1883       final double iy =
1884           uxk * guy[2]
1885               + uyk * guy[3]
1886               + uzk * guy[4]
1887               + 0.5
1888               * (ck * guy[1]
1889               + qxxk * guy[5]
1890               + qyyk * guy[8]
1891               + qzzk * guy[10]
1892               + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9])
1893               + ck * gc[3]
1894               + qxxk * gqxx[3]
1895               + qyyk * gqyy[3]
1896               + qzzk * gqzz[3]
1897               + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]));
1898       final double iz =
1899           uxk * guz[2]
1900               + uyk * guz[3]
1901               + uzk * guz[4]
1902               + 0.5
1903               * (ck * guz[1]
1904               + qxxk * guz[5]
1905               + qyyk * guz[8]
1906               + qzzk * guz[10]
1907               + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9])
1908               + ck * gc[4]
1909               + qxxk * gqxx[4]
1910               + qyyk * gqyy[4]
1911               + qzzk * gqzz[4]
1912               + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]));
1913       final double kx =
1914           uxi * gux[2]
1915               + uyi * gux[3]
1916               + uzi * gux[4]
1917               - 0.5
1918               * (ci * gux[1]
1919               + qxxi * gux[5]
1920               + qyyi * gux[8]
1921               + qzzi * gux[10]
1922               + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9])
1923               + ci * gc[2]
1924               + qxxi * gqxx[2]
1925               + qyyi * gqyy[2]
1926               + qzzi * gqzz[2]
1927               + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]));
1928       final double ky =
1929           uxi * guy[2]
1930               + uyi * guy[3]
1931               + uzi * guy[4]
1932               - 0.5
1933               * (ci * guy[1]
1934               + qxxi * guy[5]
1935               + qyyi * guy[8]
1936               + qzzi * guy[10]
1937               + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9])
1938               + ci * gc[3]
1939               + qxxi * gqxx[3]
1940               + qyyi * gqyy[3]
1941               + qzzi * gqzz[3]
1942               + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]));
1943       final double kz =
1944           uxi * guz[2]
1945               + uyi * guz[3]
1946               + uzi * guz[4]
1947               - 0.5
1948               * (ci * guz[1]
1949               + qxxi * guz[5]
1950               + qyyi * guz[8]
1951               + qzzi * guz[10]
1952               + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9])
1953               + ci * gc[4]
1954               + qxxi * gqxx[4]
1955               + qyyi * gqyy[4]
1956               + qzzi * gqzz[4]
1957               + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
1958       double tix = uyi * iz - uzi * iy;
1959       double tiy = uzi * ix - uxi * iz;
1960       double tiz = uxi * iy - uyi * ix;
1961       double tkx = uyk * kz - uzk * ky;
1962       double tky = uzk * kx - uxk * kz;
1963       double tkz = uxk * ky - uyk * kx;
1964 
1965       // Torque on quadrupoles due to permanent reaction field gradient.
1966       final double ixx =
1967           -0.5
1968               * (ck * gqxx[1]
1969               + uxk * gqxx[2]
1970               + uyk * gqxx[3]
1971               + uzk * gqxx[4]
1972               + qxxk * gqxx[5]
1973               + qyyk * gqxx[8]
1974               + qzzk * gqxx[10]
1975               + 2.0 * (qxyk * gqxx[6] + qxzk * gqxx[7] + qyzk * gqxx[9])
1976               + ck * gc[5]
1977               + uxk * gux[5]
1978               + uyk * guy[5]
1979               + uzk * guz[5]
1980               + qxxk * gqxx[5]
1981               + qyyk * gqyy[5]
1982               + qzzk * gqzz[5]
1983               + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]));
1984       final double ixy =
1985           -0.5
1986               * (ck * gqxy[1]
1987               + uxk * gqxy[2]
1988               + uyk * gqxy[3]
1989               + uzk * gqxy[4]
1990               + qxxk * gqxy[5]
1991               + qyyk * gqxy[8]
1992               + qzzk * gqxy[10]
1993               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxy[7] + qyzk * gqxy[9])
1994               + ck * gc[6]
1995               + uxk * gux[6]
1996               + uyk * guy[6]
1997               + uzk * guz[6]
1998               + qxxk * gqxx[6]
1999               + qyyk * gqyy[6]
2000               + qzzk * gqzz[6]
2001               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]));
2002       final double ixz =
2003           -0.5
2004               * (ck * gqxz[1]
2005               + uxk * gqxz[2]
2006               + uyk * gqxz[3]
2007               + uzk * gqxz[4]
2008               + qxxk * gqxz[5]
2009               + qyyk * gqxz[8]
2010               + qzzk * gqxz[10]
2011               + 2.0 * (qxyk * gqxz[6] + qxzk * gqxz[7] + qyzk * gqxz[9])
2012               + ck * gc[7]
2013               + uxk * gux[7]
2014               + uyk * guy[7]
2015               + uzk * guz[7]
2016               + qxxk * gqxx[7]
2017               + qyyk * gqyy[7]
2018               + qzzk * gqzz[7]
2019               + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]));
2020       final double iyy =
2021           -0.5
2022               * (ck * gqyy[1]
2023               + uxk * gqyy[2]
2024               + uyk * gqyy[3]
2025               + uzk * gqyy[4]
2026               + qxxk * gqyy[5]
2027               + qyyk * gqyy[8]
2028               + qzzk * gqyy[10]
2029               + 2.0 * (qxyk * gqyy[6] + qxzk * gqyy[7] + qyzk * gqyy[9])
2030               + ck * gc[8]
2031               + uxk * gux[8]
2032               + uyk * guy[8]
2033               + uzk * guz[8]
2034               + qxxk * gqxx[8]
2035               + qyyk * gqyy[8]
2036               + qzzk * gqzz[8]
2037               + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]));
2038       final double iyz =
2039           -0.5
2040               * (ck * gqyz[1]
2041               + uxk * gqyz[2]
2042               + uyk * gqyz[3]
2043               + uzk * gqyz[4]
2044               + qxxk * gqyz[5]
2045               + qyyk * gqyz[8]
2046               + qzzk * gqyz[10]
2047               + 2.0 * (qxyk * gqyz[6] + qxzk * gqyz[7] + qyzk * gqyz[9])
2048               + ck * gc[9]
2049               + uxk * gux[9]
2050               + uyk * guy[9]
2051               + uzk * guz[9]
2052               + qxxk * gqxx[9]
2053               + qyyk * gqyy[9]
2054               + qzzk * gqzz[9]
2055               + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]));
2056       final double izz =
2057           -0.5
2058               * (ck * gqzz[1]
2059               + uxk * gqzz[2]
2060               + uyk * gqzz[3]
2061               + uzk * gqzz[4]
2062               + qxxk * gqzz[5]
2063               + qyyk * gqzz[8]
2064               + qzzk * gqzz[10]
2065               + 2.0 * (qxyk * gqzz[6] + qxzk * gqzz[7] + qyzk * gqzz[9])
2066               + ck * gc[10]
2067               + uxk * gux[10]
2068               + uyk * guy[10]
2069               + uzk * guz[10]
2070               + qxxk * gqxx[10]
2071               + qyyk * gqyy[10]
2072               + qzzk * gqzz[10]
2073               + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]));
2074       final double iyx = ixy;
2075       final double izx = ixz;
2076       final double izy = iyz;
2077       final double kxx =
2078           -0.5
2079               * (ci * gqxx[1]
2080               - uxi * gqxx[2]
2081               - uyi * gqxx[3]
2082               - uzi * gqxx[4]
2083               + qxxi * gqxx[5]
2084               + qyyi * gqxx[8]
2085               + qzzi * gqxx[10]
2086               + 2.0 * (qxyi * gqxx[6] + qxzi * gqxx[7] + qyzi * gqxx[9])
2087               + ci * gc[5]
2088               - uxi * gux[5]
2089               - uyi * guy[5]
2090               - uzi * guz[5]
2091               + qxxi * gqxx[5]
2092               + qyyi * gqyy[5]
2093               + qzzi * gqzz[5]
2094               + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]));
2095       double kxy =
2096           -0.5
2097               * (ci * gqxy[1]
2098               - uxi * gqxy[2]
2099               - uyi * gqxy[3]
2100               - uzi * gqxy[4]
2101               + qxxi * gqxy[5]
2102               + qyyi * gqxy[8]
2103               + qzzi * gqxy[10]
2104               + 2.0 * (qxyi * gqxy[6] + qxzi * gqxy[7] + qyzi * gqxy[9])
2105               + ci * gc[6]
2106               - uxi * gux[6]
2107               - uyi * guy[6]
2108               - uzi * guz[6]
2109               + qxxi * gqxx[6]
2110               + qyyi * gqyy[6]
2111               + qzzi * gqzz[6]
2112               + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]));
2113       double kxz =
2114           -0.5
2115               * (ci * gqxz[1]
2116               - uxi * gqxz[2]
2117               - uyi * gqxz[3]
2118               - uzi * gqxz[4]
2119               + qxxi * gqxz[5]
2120               + qyyi * gqxz[8]
2121               + qzzi * gqxz[10]
2122               + 2.0 * (qxyi * gqxz[6] + qxzi * gqxz[7] + qyzi * gqxz[9])
2123               + ci * gc[7]
2124               - uxi * gux[7]
2125               - uyi * guy[7]
2126               - uzi * guz[7]
2127               + qxxi * gqxx[7]
2128               + qyyi * gqyy[7]
2129               + qzzi * gqzz[7]
2130               + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]));
2131       double kyy =
2132           -0.5
2133               * (ci * gqyy[1]
2134               - uxi * gqyy[2]
2135               - uyi * gqyy[3]
2136               - uzi * gqyy[4]
2137               + qxxi * gqyy[5]
2138               + qyyi * gqyy[8]
2139               + qzzi * gqyy[10]
2140               + 2.0 * (qxyi * gqyy[6] + qxzi * gqyy[7] + qyzi * gqyy[9])
2141               + ci * gc[8]
2142               - uxi * gux[8]
2143               - uyi * guy[8]
2144               - uzi * guz[8]
2145               + qxxi * gqxx[8]
2146               + qyyi * gqyy[8]
2147               + qzzi * gqzz[8]
2148               + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]));
2149       double kyz =
2150           -0.5
2151               * (ci * gqyz[1]
2152               - uxi * gqyz[2]
2153               - uyi * gqyz[3]
2154               - uzi * gqyz[4]
2155               + qxxi * gqyz[5]
2156               + qyyi * gqyz[8]
2157               + qzzi * gqyz[10]
2158               + 2.0 * (qxyi * gqyz[6] + qxzi * gqyz[7] + qyzi * gqyz[9])
2159               + ci * gc[9]
2160               - uxi * gux[9]
2161               - uyi * guy[9]
2162               - uzi * guz[9]
2163               + qxxi * gqxx[9]
2164               + qyyi * gqyy[9]
2165               + qzzi * gqzz[9]
2166               + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]));
2167       double kzz =
2168           -0.5
2169               * (ci * gqzz[1]
2170               - uxi * gqzz[2]
2171               - uyi * gqzz[3]
2172               - uzi * gqzz[4]
2173               + qxxi * gqzz[5]
2174               + qyyi * gqzz[8]
2175               + qzzi * gqzz[10]
2176               + 2.0 * (qxyi * gqzz[6] + qxzi * gqzz[7] + qyzi * gqzz[9])
2177               + ci * gc[10]
2178               - uxi * gux[10]
2179               - uyi * guy[10]
2180               - uzi * guz[10]
2181               + qxxi * gqxx[10]
2182               + qyyi * gqyy[10]
2183               + qzzi * gqzz[10]
2184               + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]));
2185       double kyx = kxy;
2186       double kzx = kxz;
2187       double kzy = kyz;
2188       tix += 2.0 * (qxyi * ixz + qyyi * iyz + qyzi * izz - qxzi * ixy - qyzi * iyy - qzzi * izy);
2189       tiy += 2.0 * (qxzi * ixx + qyzi * iyx + qzzi * izx - qxxi * ixz - qxyi * iyz - qxzi * izz);
2190       tiz += 2.0 * (qxxi * ixy + qxyi * iyy + qxzi * izy - qxyi * ixx - qyyi * iyx - qyzi * izx);
2191       tkx += 2.0 * (qxyk * kxz + qyyk * kyz + qyzk * kzz - qxzk * kxy - qyzk * kyy - qzzk * kzy);
2192       tky += 2.0 * (qxzk * kxx + qyzk * kyx + qzzk * kzx - qxxk * kxz - qxyk * kyz - qxzk * kzz);
2193       tkz += 2.0 * (qxxk * kxy + qxyk * kyy + qxzk * kzy - qxyk * kxx - qyyk * kyx - qyzk * kzx);
2194       if (i == k) {
2195         double selfScale = 0.5;
2196         tix *= selfScale;
2197         tiy *= selfScale;
2198         tiz *= selfScale;
2199         tkx *= selfScale;
2200         tky *= selfScale;
2201         tkz *= selfScale;
2202       }
2203       trqxi += tix;
2204       trqyi += tiy;
2205       trqzi += tiz;
2206 
2207       final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
2208       final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
2209       final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
2210       torque.add(threadID, k, rtkx, rtky, rtkz);
2211     }
2212 
2213     private void polarizationEnergyGradient(int i, int k) {
2214       // Electrostatic solvation free energy gradient of the permanent
2215       // multipoles in the reaction potential of the induced dipoles.
2216       final double dpsymdx =
2217           -uxi * (sxk * gux[5] + syk * guy[5] + szk * guz[5])
2218               - uyi * (sxk * gux[6] + syk * guy[6] + szk * guz[6])
2219               - uzi * (sxk * gux[7] + syk * guy[7] + szk * guz[7])
2220               - uxk * (sxi * gux[5] + syi * guy[5] + szi * guz[5])
2221               - uyk * (sxi * gux[6] + syi * guy[6] + szi * guz[6])
2222               - uzk * (sxi * gux[7] + syi * guy[7] + szi * guz[7]);
2223       final double dpwidx =
2224           ci * (sxk * gc[5] + syk * gc[6] + szk * gc[7])
2225               - ck * (sxi * gux[2] + syi * guy[2] + szi * guz[2])
2226               - sxi
2227               * (qxxk * gux[11]
2228               + qyyk * gux[14]
2229               + qzzk * gux[16]
2230               + 2.0 * (qxyk * gux[12] + qxzk * gux[13] + qyzk * gux[15]))
2231               - syi
2232               * (qxxk * guy[11]
2233               + qyyk * guy[14]
2234               + qzzk * guy[16]
2235               + 2.0 * (qxyk * guy[12] + qxzk * guy[13] + qyzk * guy[15]))
2236               - szi
2237               * (qxxk * guz[11]
2238               + qyyk * guz[14]
2239               + qzzk * guz[16]
2240               + 2.0 * (qxyk * guz[12] + qxzk * guz[13] + qyzk * guz[15]))
2241               + sxk
2242               * (qxxi * gqxx[5]
2243               + qyyi * gqyy[5]
2244               + qzzi * gqzz[5]
2245               + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]))
2246               + syk
2247               * (qxxi * gqxx[6]
2248               + qyyi * gqyy[6]
2249               + qzzi * gqzz[6]
2250               + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
2251               + szk
2252               * (qxxi * gqxx[7]
2253               + qyyi * gqyy[7]
2254               + qzzi * gqzz[7]
2255               + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]));
2256       final double dpwkdx =
2257           ci * (sxk * gux[2] + syk * guy[2] + szk * guz[2])
2258               - ck * (sxi * gc[5] + syi * gc[6] + szi * gc[7])
2259               - sxi
2260               * (qxxk * gqxx[5]
2261               + qyyk * gqyy[5]
2262               + qzzk * gqzz[5]
2263               + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
2264               - syi
2265               * (qxxk * gqxx[6]
2266               + qyyk * gqyy[6]
2267               + qzzk * gqzz[6]
2268               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
2269               - szi
2270               * (qxxk * gqxx[7]
2271               + qyyk * gqyy[7]
2272               + qzzk * gqzz[7]
2273               + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
2274               + sxk
2275               * (qxxi * gux[11]
2276               + qyyi * gux[14]
2277               + qzzi * gux[16]
2278               + 2.0 * (qxyi * gux[12] + qxzi * gux[13] + qyzi * gux[15]))
2279               + syk
2280               * (qxxi * guy[11]
2281               + qyyi * guy[14]
2282               + qzzi * guy[16]
2283               + 2.0 * (qxyi * guy[12] + qxzi * guy[13] + qyzi * guy[15]))
2284               + szk
2285               * (qxxi * guz[11]
2286               + qyyi * guz[14]
2287               + qzzi * guz[16]
2288               + 2.0 * (qxyi * guz[12] + qxzi * guz[13] + qyzi * guz[15]));
2289       final double dpsymdy =
2290           -uxi * (sxk * gux[6] + syk * guy[6] + szk * guz[6])
2291               - uyi * (sxk * gux[8] + syk * guy[8] + szk * guz[8])
2292               - uzi * (sxk * gux[9] + syk * guy[9] + szk * guz[9])
2293               - uxk * (sxi * gux[6] + syi * guy[6] + szi * guz[6])
2294               - uyk * (sxi * gux[8] + syi * guy[8] + szi * guz[8])
2295               - uzk * (sxi * gux[9] + syi * guy[9] + szi * guz[9]);
2296       final double dpwidy =
2297           ci * (sxk * gc[6] + syk * gc[8] + szk * gc[9])
2298               - ck * (sxi * gux[3] + syi * guy[3] + szi * guz[3])
2299               - sxi
2300               * (qxxk * gux[12]
2301               + qyyk * gux[17]
2302               + qzzk * gux[19]
2303               + 2.0 * (qxyk * gux[14] + qxzk * gux[15] + qyzk * gux[18]))
2304               - syi
2305               * (qxxk * guy[12]
2306               + qyyk * guy[17]
2307               + qzzk * guy[19]
2308               + 2.0 * (qxyk * guy[14] + qxzk * guy[15] + qyzk * guy[18]))
2309               - szi
2310               * (qxxk * guz[12]
2311               + qyyk * guz[17]
2312               + qzzk * guz[19]
2313               + 2.0 * (qxyk * guz[14] + qxzk * guz[15] + qyzk * guz[18]))
2314               + sxk
2315               * (qxxi * gqxx[6]
2316               + qyyi * gqyy[6]
2317               + qzzi * gqzz[6]
2318               + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
2319               + syk
2320               * (qxxi * gqxx[8]
2321               + qyyi * gqyy[8]
2322               + qzzi * gqzz[8]
2323               + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]))
2324               + szk
2325               * (qxxi * gqxx[9]
2326               + qyyi * gqyy[9]
2327               + qzzi * gqzz[9]
2328               + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]));
2329       final double dpwkdy =
2330           ci * (sxk * gux[3] + syk * guy[3] + szk * guz[3])
2331               - ck * (sxi * gc[6] + syi * gc[8] + szi * gc[9])
2332               - sxi
2333               * (qxxk * gqxx[6]
2334               + qyyk * gqyy[6]
2335               + qzzk * gqzz[6]
2336               + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
2337               - syi
2338               * (qxxk * gqxx[8]
2339               + qyyk * gqyy[8]
2340               + qzzk * gqzz[8]
2341               + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
2342               - szi
2343               * (qxxk * gqxx[9]
2344               + qyyk * gqyy[9]
2345               + qzzk * gqzz[9]
2346               + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
2347               + sxk
2348               * (qxxi * gux[12]
2349               + qyyi * gux[17]
2350               + qzzi * gux[19]
2351               + 2.0 * (qxyi * gux[14] + qxzi * gux[15] + qyzi * gux[18]))
2352               + syk
2353               * (qxxi * guy[12]
2354               + qyyi * guy[17]
2355               + qzzi * guy[19]
2356               + 2.0 * (qxyi * guy[14] + qxzi * guy[15] + qyzi * guy[18]))
2357               + szk
2358               * (qxxi * guz[12]
2359               + qyyi * guz[17]
2360               + qzzi * guz[19]
2361               + 2.0 * (qxyi * guz[14] + qxzi * guz[15] + qyzi * guz[18]));
2362       final double dpsymdz =
2363           -uxi * (sxk * gux[7] + syk * guy[7] + szk * guz[7])
2364               - uyi * (sxk * gux[9] + syk * guy[9] + szk * guz[9])
2365               - uzi * (sxk * gux[10] + syk * guy[10] + szk * guz[10])
2366               - uxk * (sxi * gux[7] + syi * guy[7] + szi * guz[7])
2367               - uyk * (sxi * gux[9] + syi * guy[9] + szi * guz[9])
2368               - uzk * (sxi * gux[10] + syi * guy[10] + szi * guz[10]);
2369       final double dpwidz =
2370           ci * (sxk * gc[7] + syk * gc[9] + szk * gc[10])
2371               - ck * (sxi * gux[4] + syi * guy[4] + szi * guz[4])
2372               - sxi
2373               * (qxxk * gux[13]
2374               + qyyk * gux[18]
2375               + qzzk * gux[20]
2376               + 2.0 * (qxyk * gux[15] + qxzk * gux[16] + qyzk * gux[19]))
2377               - syi
2378               * (qxxk * guy[13]
2379               + qyyk * guy[18]
2380               + qzzk * guy[20]
2381               + 2.0 * (qxyk * guy[15] + qxzk * guy[16] + qyzk * guy[19]))
2382               - szi
2383               * (qxxk * guz[13]
2384               + qyyk * guz[18]
2385               + qzzk * guz[20]
2386               + 2.0 * (qxyk * guz[15] + qxzk * guz[16] + qyzk * guz[19]))
2387               + sxk
2388               * (qxxi * gqxx[7]
2389               + qyyi * gqyy[7]
2390               + qzzi * gqzz[7]
2391               + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
2392               + syk
2393               * (qxxi * gqxx[9]
2394               + qyyi * gqyy[9]
2395               + qzzi * gqzz[9]
2396               + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
2397               + szk
2398               * (qxxi * gqxx[10]
2399               + qyyi * gqyy[10]
2400               + qzzi * gqzz[10]
2401               + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]));
2402       final double dpwkdz =
2403           ci * (sxk * gux[4] + syk * guy[4] + szk * guz[4])
2404               - ck * (sxi * gc[7] + syi * gc[9] + szi * gc[10])
2405               - sxi
2406               * (qxxk * gqxx[7]
2407               + qyyk * gqyy[7]
2408               + qzzk * gqzz[7]
2409               + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
2410               - syi
2411               * (qxxk * gqxx[9]
2412               + qyyk * gqyy[9]
2413               + qzzk * gqzz[9]
2414               + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
2415               - szi
2416               * (qxxk * gqxx[10]
2417               + qyyk * gqyy[10]
2418               + qzzk * gqzz[10]
2419               + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
2420               + sxk
2421               * (qxxi * gux[13]
2422               + qyyi * gux[18]
2423               + qzzi * gux[20]
2424               + 2.0 * (qxyi * gux[15] + qxzi * gux[16] + qyzi * gux[19]))
2425               + syk
2426               * (qxxi * guy[13]
2427               + qyyi * guy[18]
2428               + qzzi * guy[20]
2429               + 2.0 * (qxyi * guy[15] + qxzi * guy[16] + qyzi * guy[19]))
2430               + szk
2431               * (qxxi * guz[13]
2432               + qyyi * guz[18]
2433               + qzzi * guz[20]
2434               + 2.0 * (qxyi * guz[15] + qxzi * guz[16] + qyzi * guz[19]));
2435 
2436       // Effective radii chain rule terms for the electrostatic solvation free energy
2437       // gradient of the permanent multipoles in the reaction potential of the induced dipoles.
2438       final double dsymdr =
2439           -uxi * (sxk * gux[22] + syk * guy[22] + szk * guz[22])
2440               - uyi * (sxk * gux[23] + syk * guy[23] + szk * guz[23])
2441               - uzi * (sxk * gux[24] + syk * guy[24] + szk * guz[24])
2442               - uxk * (sxi * gux[22] + syi * guy[22] + szi * guz[22])
2443               - uyk * (sxi * gux[23] + syi * guy[23] + szi * guz[23])
2444               - uzk * (sxi * gux[24] + syi * guy[24] + szi * guz[24]);
2445       final double dwipdr =
2446           ci * (sxk * gc[22] + syk * gc[23] + szk * gc[24])
2447               - ck * (sxi * gux[21] + syi * guy[21] + szi * guz[21])
2448               - sxi
2449               * (qxxk * gux[25]
2450               + qyyk * gux[28]
2451               + qzzk * gux[30]
2452               + 2.0 * (qxyk * gux[26] + qxzk * gux[27] + qyzk * gux[29]))
2453               - syi
2454               * (qxxk * guy[25]
2455               + qyyk * guy[28]
2456               + qzzk * guy[30]
2457               + 2.0 * (qxyk * guy[26] + qxzk * guy[27] + qyzk * guy[29]))
2458               - szi
2459               * (qxxk * guz[25]
2460               + qyyk * guz[28]
2461               + qzzk * guz[30]
2462               + 2.0 * (qxyk * guz[26] + qxzk * guz[27] + qyzk * guz[29]))
2463               + sxk
2464               * (qxxi * gqxx[22]
2465               + qyyi * gqyy[22]
2466               + qzzi * gqzz[22]
2467               + 2.0 * (qxyi * gqxy[22] + qxzi * gqxz[22] + qyzi * gqyz[22]))
2468               + syk
2469               * (qxxi * gqxx[23]
2470               + qyyi * gqyy[23]
2471               + qzzi * gqzz[23]
2472               + 2.0 * (qxyi * gqxy[23] + qxzi * gqxz[23] + qyzi * gqyz[23]))
2473               + szk
2474               * (qxxi * gqxx[24]
2475               + qyyi * gqyy[24]
2476               + qzzi * gqzz[24]
2477               + 2.0 * (qxyi * gqxy[24] + qxzi * gqxz[24] + qyzi * gqyz[24]));
2478       final double dwkpdr =
2479           ci * (sxk * gux[21] + syk * guy[21] + szk * guz[21])
2480               - ck * (sxi * gc[22] + syi * gc[23] + szi * gc[24])
2481               - sxi
2482               * (qxxk * gqxx[22]
2483               + qyyk * gqyy[22]
2484               + qzzk * gqzz[22]
2485               + 2.0 * (qxyk * gqxy[22] + qxzk * gqxz[22] + qyzk * gqyz[22]))
2486               - syi
2487               * (qxxk * gqxx[23]
2488               + qyyk * gqyy[23]
2489               + qzzk * gqzz[23]
2490               + 2.0 * (qxyk * gqxy[23] + qxzk * gqxz[23] + qyzk * gqyz[23]))
2491               - szi
2492               * (qxxk * gqxx[24]
2493               + qyyk * gqyy[24]
2494               + qzzk * gqzz[24]
2495               + 2.0 * (qxyk * gqxy[24] + qxzk * gqxz[24] + qyzk * gqyz[24]))
2496               + sxk
2497               * (qxxi * gux[25]
2498               + qyyi * gux[28]
2499               + qzzi * gux[30]
2500               + 2.0 * (qxyi * gux[26] + qxzi * gux[27] + qyzi * gux[29]))
2501               + syk
2502               * (qxxi * guy[25]
2503               + qyyi * guy[28]
2504               + qzzi * guy[30]
2505               + 2.0 * (qxyi * guy[26] + qxzi * guy[27] + qyzi * guy[29]))
2506               + szk
2507               * (qxxi * guz[25]
2508               + qyyi * guz[28]
2509               + qzzi * guz[30]
2510               + 2.0 * (qxyi * guz[26] + qxzi * guz[27] + qyzi * guz[29]));
2511       double dpdx = 0.5 * (dpsymdx + 0.5 * (dpwidx + dpwkdx));
2512       double dpdy = 0.5 * (dpsymdy + 0.5 * (dpwidy + dpwkdy));
2513       double dpdz = 0.5 * (dpsymdz + 0.5 * (dpwidz + dpwkdz));
2514       double dsumdri = dsymdr + 0.5 * (dwipdr + dwkpdr);
2515       double dbi = 0.5 * rbk * dsumdri;
2516       double dbk = 0.5 * rbi * dsumdri;
2517       if (polarization == Polarization.MUTUAL) {
2518         dpdx -=
2519             0.5
2520                 * (dxi * (pxk * gux[5] + pyk * gux[6] + pzk * gux[7])
2521                 + dyi * (pxk * guy[5] + pyk * guy[6] + pzk * guy[7])
2522                 + dzi * (pxk * guz[5] + pyk * guz[6] + pzk * guz[7])
2523                 + dxk * (pxi * gux[5] + pyi * gux[6] + pzi * gux[7])
2524                 + dyk * (pxi * guy[5] + pyi * guy[6] + pzi * guy[7])
2525                 + dzk * (pxi * guz[5] + pyi * guz[6] + pzi * guz[7]));
2526         dpdy -=
2527             0.5
2528                 * (dxi * (pxk * gux[6] + pyk * gux[8] + pzk * gux[9])
2529                 + dyi * (pxk * guy[6] + pyk * guy[8] + pzk * guy[9])
2530                 + dzi * (pxk * guz[6] + pyk * guz[8] + pzk * guz[9])
2531                 + dxk * (pxi * gux[6] + pyi * gux[8] + pzi * gux[9])
2532                 + dyk * (pxi * guy[6] + pyi * guy[8] + pzi * guy[9])
2533                 + dzk * (pxi * guz[6] + pyi * guz[8] + pzi * guz[9]));
2534         dpdz -=
2535             0.5
2536                 * (dxi * (pxk * gux[7] + pyk * gux[9] + pzk * gux[10])
2537                 + dyi * (pxk * guy[7] + pyk * guy[9] + pzk * guy[10])
2538                 + dzi * (pxk * guz[7] + pyk * guz[9] + pzk * guz[10])
2539                 + dxk * (pxi * gux[7] + pyi * gux[9] + pzi * gux[10])
2540                 + dyk * (pxi * guy[7] + pyi * guy[9] + pzi * guy[10])
2541                 + dzk * (pxi * guz[7] + pyi * guz[9] + pzi * guz[10]));
2542         final double duvdr =
2543             dxi * (pxk * gux[22] + pyk * gux[23] + pzk * gux[24])
2544                 + dyi * (pxk * guy[22] + pyk * guy[23] + pzk * guy[24])
2545                 + dzi * (pxk * guz[22] + pyk * guz[23] + pzk * guz[24])
2546                 + dxk * (pxi * gux[22] + pyi * gux[23] + pzi * gux[24])
2547                 + dyk * (pxi * guy[22] + pyi * guy[23] + pzi * guy[24])
2548                 + dzk * (pxi * guz[22] + pyi * guz[23] + pzi * guz[24]);
2549         dbi -= 0.5 * rbk * duvdr;
2550         dbk -= 0.5 * rbi * duvdr;
2551       }
2552 
2553       // Increment the gradients and Born chain rule term.
2554       if (i == k && iSymm == 0) {
2555         dborni += dbi;
2556       } else {
2557         if (i == k) {
2558           dpdx *= 0.5;
2559           dpdy *= 0.5;
2560           dpdz *= 0.5;
2561           dbi *= 0.5;
2562           dbk *= 0.5;
2563         }
2564         dedxi -= dpdx;
2565         dedyi -= dpdy;
2566         dedzi -= dpdz;
2567         dborni += dbi;
2568         final double rdpdx = dpdx * transOp[0][0] + dpdy * transOp[1][0] + dpdz * transOp[2][0];
2569         final double rdpdy = dpdx * transOp[0][1] + dpdy * transOp[1][1] + dpdz * transOp[2][1];
2570         final double rdpdz = dpdx * transOp[0][2] + dpdy * transOp[1][2] + dpdz * transOp[2][2];
2571         grad.add(threadID, k, rdpdx, rdpdy, rdpdz);
2572         sharedBornGrad.add(threadID, k, dbk);
2573       }
2574       polarizationEnergyTorque(i, k);
2575     }
2576 
2577     private void polarizationEnergyTorque(int i, int k) {
2578       double fix = 0.5 * (sxk * gux[2] + syk * guy[2] + szk * guz[2]);
2579       double fiy = 0.5 * (sxk * gux[3] + syk * guy[3] + szk * guz[3]);
2580       double fiz = 0.5 * (sxk * gux[4] + syk * guy[4] + szk * guz[4]);
2581       double fkx = 0.5 * (sxi * gux[2] + syi * guy[2] + szi * guz[2]);
2582       double fky = 0.5 * (sxi * gux[3] + syi * guy[3] + szi * guz[3]);
2583       double fkz = 0.5 * (sxi * gux[4] + syi * guy[4] + szi * guz[4]);
2584       double fixx =
2585           -0.25
2586               * ((sxk * gqxx[2] + syk * gqxx[3] + szk * gqxx[4])
2587               + (sxk * gux[5] + syk * guy[5] + szk * guz[5]));
2588       double fixy =
2589           -0.25
2590               * ((sxk * gqxy[2] + syk * gqxy[3] + szk * gqxy[4])
2591               + (sxk * gux[6] + syk * guy[6] + szk * guz[6]));
2592       double fixz =
2593           -0.25
2594               * ((sxk * gqxz[2] + syk * gqxz[3] + szk * gqxz[4])
2595               + (sxk * gux[7] + syk * guy[7] + szk * guz[7]));
2596       double fiyy =
2597           -0.25
2598               * ((sxk * gqyy[2] + syk * gqyy[3] + szk * gqyy[4])
2599               + (sxk * gux[8] + syk * guy[8] + szk * guz[8]));
2600       double fiyz =
2601           -0.25
2602               * ((sxk * gqyz[2] + syk * gqyz[3] + szk * gqyz[4])
2603               + (sxk * gux[9] + syk * guy[9] + szk * guz[9]));
2604       double fizz =
2605           -0.25
2606               * ((sxk * gqzz[2] + syk * gqzz[3] + szk * gqzz[4])
2607               + (sxk * gux[10] + syk * guy[10] + szk * guz[10]));
2608       double fiyx = fixy;
2609       double fizx = fixz;
2610       double fizy = fiyz;
2611       double fkxx =
2612           0.25
2613               * ((sxi * gqxx[2] + syi * gqxx[3] + szi * gqxx[4])
2614               + (sxi * gux[5] + syi * guy[5] + szi * guz[5]));
2615       double fkxy =
2616           0.25
2617               * ((sxi * gqxy[2] + syi * gqxy[3] + szi * gqxy[4])
2618               + (sxi * gux[6] + syi * guy[6] + szi * guz[6]));
2619       double fkxz =
2620           0.25
2621               * ((sxi * gqxz[2] + syi * gqxz[3] + szi * gqxz[4])
2622               + (sxi * gux[7] + syi * guy[7] + szi * guz[7]));
2623       double fkyy =
2624           0.25
2625               * ((sxi * gqyy[2] + syi * gqyy[3] + szi * gqyy[4])
2626               + (sxi * gux[8] + syi * guy[8] + szi * guz[8]));
2627       double fkyz =
2628           0.25
2629               * ((sxi * gqyz[2] + syi * gqyz[3] + szi * gqyz[4])
2630               + (sxi * gux[9] + syi * guy[9] + szi * guz[9]));
2631       double fkzz =
2632           0.25
2633               * ((sxi * gqzz[2] + syi * gqzz[3] + szi * gqzz[4])
2634               + (sxi * gux[10] + syi * guy[10] + szi * guz[10]));
2635       double fkyx = fkxy;
2636       double fkzx = fkxz;
2637       double fkzy = fkyz;
2638       if (i == k) {
2639         fix *= 0.5;
2640         fiy *= 0.5;
2641         fiz *= 0.5;
2642         fkx *= 0.5;
2643         fky *= 0.5;
2644         fkz *= 0.5;
2645         fixx *= 0.5;
2646         fixy *= 0.5;
2647         fixz *= 0.5;
2648         fiyx *= 0.5;
2649         fiyy *= 0.5;
2650         fiyz *= 0.5;
2651         fizx *= 0.5;
2652         fizy *= 0.5;
2653         fizz *= 0.5;
2654         fkxx *= 0.5;
2655         fkxy *= 0.5;
2656         fkxz *= 0.5;
2657         fkyx *= 0.5;
2658         fkyy *= 0.5;
2659         fkyz *= 0.5;
2660         fkzx *= 0.5;
2661         fkzy *= 0.5;
2662         fkzz *= 0.5;
2663       }
2664 
2665       // Torque due to induced reaction field on permanent dipoles.
2666       double tix = uyi * fiz - uzi * fiy;
2667       double tiy = uzi * fix - uxi * fiz;
2668       double tiz = uxi * fiy - uyi * fix;
2669       double tkx = uyk * fkz - uzk * fky;
2670       double tky = uzk * fkx - uxk * fkz;
2671       double tkz = uxk * fky - uyk * fkx;
2672 
2673       // Torque due to induced reaction field gradient on quadrupoles.
2674       tix +=
2675           2.0 * (qxyi * fixz + qyyi * fiyz + qyzi * fizz - qxzi * fixy - qyzi * fiyy - qzzi * fizy);
2676       tiy +=
2677           2.0 * (qxzi * fixx + qyzi * fiyx + qzzi * fizx - qxxi * fixz - qxyi * fiyz - qxzi * fizz);
2678       tiz +=
2679           2.0 * (qxxi * fixy + qxyi * fiyy + qxzi * fizy - qxyi * fixx - qyyi * fiyx - qyzi * fizx);
2680       tkx +=
2681           2.0 * (qxyk * fkxz + qyyk * fkyz + qyzk * fkzz - qxzk * fkxy - qyzk * fkyy - qzzk * fkzy);
2682       tky +=
2683           2.0 * (qxzk * fkxx + qyzk * fkyx + qzzk * fkzx - qxxk * fkxz - qxyk * fkyz - qxzk * fkzz);
2684       tkz +=
2685           2.0 * (qxxk * fkxy + qxyk * fkyy + qxzk * fkzy - qxyk * fkxx - qyyk * fkyx - qyzk * fkzx);
2686       trqxi += tix;
2687       trqyi += tiy;
2688       trqzi += tiz;
2689 
2690       final double rx = tkx;
2691       final double ry = tky;
2692       final double rz = tkz;
2693       tkx = rx * transOp[0][0] + ry * transOp[1][0] + rz * transOp[2][0];
2694       tky = rx * transOp[0][1] + ry * transOp[1][1] + rz * transOp[2][1];
2695       tkz = rx * transOp[0][2] + ry * transOp[1][2] + rz * transOp[2][2];
2696       torque.add(threadID, k, tkx, tky, tkz);
2697     }
2698   }
2699 
2700   /**
2701    * Compute the GK Energy using a QI frame.
2702    *
2703    * @since 1.0
2704    */
2705   private class GKEnergyLoopQI extends IntegerForLoop {
2706 
2707     private final double[] dx_local;
2708     private final double[] gradI = new double[3];
2709     private final double[] torqueI = new double[3];
2710     private final double[] torqueK = new double[3];
2711     private final double[] gI = new double[3];
2712     private final double[] tI = new double[3];
2713     private final double[] tK = new double[3];
2714 
2715     private double rbi;
2716     PolarizableMultipole mI;
2717     PolarizableMultipole mK;
2718     QIFrame qiFrame;
2719     GKEnergyQI gkEnergyQI;
2720     private double xi, yi, zi;
2721     private double dedxi, dedyi, dedzi;
2722     private double dborni;
2723     private double trqxi, trqyi, trqzi;
2724     private int count;
2725     private int iSymm;
2726     private int threadID;
2727     private final double[][] transOp;
2728     private double gkEnergy;
2729     private double gkPermanentEnergy;
2730     private double gkPolarizationEnergy;
2731 
2732     GKEnergyLoopQI() {
2733       mI = new PolarizableMultipole();
2734       mK = new PolarizableMultipole();
2735       qiFrame = new QIFrame();
2736       dx_local = new double[3];
2737       transOp = new double[3][3];
2738     }
2739 
2740     @Override
2741     public void start() {
2742       gkEnergy = 0.0;
2743       gkPermanentEnergy = 0.0;
2744       gkPolarizationEnergy = 0.0;
2745       count = 0;
2746       threadID = getThreadIndex();
2747       gkEnergyQI = new GKEnergyQI(soluteDielectric, solventDielectric, gkc, gradient);
2748     }
2749 
2750     @Override
2751     public void finish() {
2752       sharedInteractions.addAndGet(count);
2753       sharedGKEnergy.addAndGet(gkEnergy);
2754       sharedPermanentGKEnergy.addAndGet(gkPermanentEnergy);
2755       sharedPolarizationGKEnergy.addAndGet(gkPolarizationEnergy);
2756     }
2757 
2758     @Override
2759     public void run(int lb, int ub) {
2760 
2761       double[] x = sXYZ[0][0];
2762       double[] y = sXYZ[0][1];
2763       double[] z = sXYZ[0][2];
2764 
2765       int nSymm = crystal.spaceGroup.symOps.size();
2766       List<SymOp> symOps = crystal.spaceGroup.symOps;
2767 
2768       for (iSymm = 0; iSymm < nSymm; iSymm++) {
2769         SymOp symOp = symOps.get(iSymm);
2770         crystal.getTransformationOperator(symOp, transOp);
2771         for (int i = lb; i <= ub; i++) {
2772           if (!use[i]) {
2773             continue;
2774           }
2775 
2776           // Zero out force accumulation for atom i.
2777           dedxi = 0.0;
2778           dedyi = 0.0;
2779           dedzi = 0.0;
2780           trqxi = 0.0;
2781           trqyi = 0.0;
2782           trqzi = 0.0;
2783           dborni = 0.0;
2784 
2785           // Collect the coordinates of atom i.
2786           xi = x[i];
2787           yi = y[i];
2788           zi = z[i];
2789           rbi = born[i];
2790           int[] list = neighborLists[iSymm][i];
2791           for (int k : list) {
2792             if (!use[k]) {
2793               continue;
2794             }
2795             interaction(i, k);
2796           }
2797           if (iSymm == 0) {
2798             // Include self-interactions for the asymmetric unit atoms.
2799             interaction(i, i);
2800             /*
2801              Formula for Born energy approximation for cavitation energy is:
2802              e = surfaceTension / 6 * (ri + probe)^2 * (ri/rb)^6.
2803              ri is the base atomic radius the atom.
2804              rb is Born radius of the atom.
2805             */
2806             switch (nonPolarModel) {
2807               case BORN_SOLV:
2808               case BORN_CAV_DISP:
2809                 double r = baseRadius[i] + dOffset + probe;
2810                 double ratio = (baseRadius[i] + dOffset) / born[i];
2811                 ratio *= ratio;
2812                 ratio *= (ratio * ratio);
2813                 double saTerm = surfaceTension * r * r * ratio / 6.0;
2814                 gkEnergy += saTerm;
2815                 sharedBornGrad.sub(threadID, i, 6.0 * saTerm / born[i]);
2816                 break;
2817               default:
2818                 break;
2819             }
2820           }
2821           if (gradient) {
2822             grad.add(threadID, i, dedxi, dedyi, dedzi);
2823             torque.add(threadID, i, trqxi, trqyi, trqzi);
2824             sharedBornGrad.add(threadID, i, dborni);
2825           }
2826         }
2827       }
2828     }
2829 
2830     private void interaction(int i, int k) {
2831       dx_local[0] = sXYZ[iSymm][0][k] - xi;
2832       dx_local[1] = sXYZ[iSymm][1][k] - yi;
2833       dx_local[2] = sXYZ[iSymm][2][k] - zi;
2834       double r2 = crystal.image(dx_local);
2835       if (r2 > cut2) {
2836         return;
2837       }
2838 
2839       // Set the multipole moments for site I.
2840       mI.setPermanentMultipole(globalMultipole[0][i]);
2841       mI.setInducedDipole(inducedDipole[0][i], inducedDipoleCR[0][i]);
2842       // Set the multipole moments for site K.
2843       mK.setPermanentMultipole(globalMultipole[iSymm][k]);
2844       mK.setInducedDipole(inducedDipole[iSymm][k], inducedDipoleCR[iSymm][k]);
2845 
2846       double selfScale = 1.0;
2847       if (i == k) {
2848         selfScale = 0.5;
2849       } else {
2850         qiFrame.setAndRotate(dx_local, mI, mK);
2851       }
2852 
2853       // Set the GK energy tensor.
2854       double rbk = born[k];
2855       gkEnergyQI.initPotential(dx_local, r2, rbi, rbk);
2856 
2857       // Compute the GK permanent multipole interaction.
2858       double eik;
2859       if (!gradient) {
2860         // Compute the GK permanent interaction energy.
2861         eik = electric * selfScale * gkEnergyQI.multipoleEnergy(mI, mK);
2862         gkPermanentEnergy += eik;
2863         if (polarization != Polarization.NONE) {
2864           // Compute the GK polarization interaction energy.
2865           double ep = electric * selfScale * gkEnergyQI.polarizationEnergy(mI, mK);
2866           gkPolarizationEnergy += ep;
2867           eik += ep;
2868         }
2869       } else {
2870         if (i == k) {
2871           // Compute the GK permanent interaction energy.
2872           eik = electric * selfScale * gkEnergyQI.multipoleEnergy(mI, mK);
2873           gkPermanentEnergy += eik;
2874           // There is no dE/dR or torque for the i == k permanent self-energy.
2875           if (polarization != Polarization.NONE) {
2876             // Compute the GK polarization interaction energy.
2877             // There is no dE/dR, but there can be a torque for the i == k polarization self-energy
2878             // if the permanent dipole and induced dipole are not aligned.
2879             double mutualMask = 1.0;
2880             if (polarization == Polarization.DIRECT) {
2881               mutualMask = 0.0;
2882             }
2883             double ep =
2884                 electric * selfScale * gkEnergyQI.polarizationEnergyAndGradient(mI, mK, mutualMask,
2885                     gradI, torqueI, torqueK);
2886             gkPolarizationEnergy += ep;
2887             eik += ep;
2888             // The torque for the i == k polarization self-energy.
2889             trqxi += electric * selfScale * torqueI[0];
2890             trqyi += electric * selfScale * torqueI[1];
2891             trqzi += electric * selfScale * torqueI[2];
2892             double tkx = electric * selfScale * torqueK[0];
2893             double tky = electric * selfScale * torqueK[1];
2894             double tkz = electric * selfScale * torqueK[2];
2895             final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
2896             final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
2897             final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
2898             torque.add(threadID, k, rtkx, rtky, rtkz);
2899           }
2900           // Compute the GK permanent Born chain-rule term.
2901           gkEnergyQI.initBorn(dx_local, r2, rbi, rbk);
2902           double db = gkEnergyQI.multipoleEnergyBornGrad(mI, mK);
2903           if (polarization != Polarization.NONE) {
2904             db += gkEnergyQI.polarizationEnergyBornGrad(mI, mK, polarization == Polarization.MUTUAL);
2905           }
2906           dborni += electric * rbi * db;
2907         } else {
2908           // Sum the GK permanent interaction energy.
2909           eik = electric * gkEnergyQI.multipoleEnergyAndGradient(mI, mK, gradI, torqueI, torqueK);
2910           gkPermanentEnergy += eik;
2911           if (polarization != Polarization.NONE) {
2912             double mutualMask = 1.0;
2913             if (polarization == Polarization.DIRECT) {
2914               mutualMask = 0.0;
2915             }
2916             // Sum the GK polarization interaction energy.
2917             double ep =
2918                 electric * gkEnergyQI.polarizationEnergyAndGradient(mI, mK, mutualMask, gI, tI, tK);
2919             eik += ep;
2920             gkPolarizationEnergy += ep;
2921             for (int j = 0; j < 3; j++) {
2922               gradI[j] += gI[j];
2923               torqueI[j] += tI[j];
2924               torqueK[j] += tK[j];
2925             }
2926           }
2927 
2928           // Convert to units of kcal/mol
2929           for (int j = 0; j < 3; j++) {
2930             gradI[j] *= electric;
2931             torqueI[j] *= electric;
2932             torqueK[j] *= electric;
2933           }
2934 
2935           // Rotate gradient and torques from the QI frame into the Global frame.
2936           qiFrame.toGlobal(gradI);
2937           qiFrame.toGlobal(torqueI);
2938           qiFrame.toGlobal(torqueK);
2939 
2940           // Accumulate partial derivatives on site I.
2941           dedxi += gradI[0];
2942           dedyi += gradI[1];
2943           dedzi += gradI[2];
2944           trqxi += torqueI[0];
2945           trqyi += torqueI[1];
2946           trqzi += torqueI[2];
2947 
2948           // Accumulate partial derivatives on site K.
2949           double dedx = -gradI[0];
2950           double dedy = -gradI[1];
2951           double dedz = -gradI[2];
2952           // Handle symmetry mate rotation.
2953           final double dedxk = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
2954           final double dedyk = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
2955           final double dedzk = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
2956           double tkx = torqueK[0];
2957           double tky = torqueK[1];
2958           double tkz = torqueK[2];
2959           final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
2960           final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
2961           final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
2962           grad.add(threadID, k, dedxk, dedyk, dedzk);
2963           torque.add(threadID, k, rtkx, rtky, rtkz);
2964 
2965           // Compute the Born-chain rule term.
2966           gkEnergyQI.initBorn(dx_local, r2, rbi, rbk);
2967           double db = gkEnergyQI.multipoleEnergyBornGrad(mI, mK);
2968           if (polarization != Polarization.NONE) {
2969             db += gkEnergyQI.polarizationEnergyBornGrad(mI, mK, polarization == Polarization.MUTUAL);
2970           }
2971           db *= electric;
2972           dborni += rbk * db;
2973           sharedBornGrad.add(threadID, k, rbi * db);
2974         }
2975       }
2976 
2977       if (i == k) {
2978         selfEnergy.add(threadID, i, eik);
2979       } else {
2980         double half = 0.5 * eik;
2981         crossEnergy.add(threadID, i, half);
2982         crossEnergy.add(threadID, k, half);
2983       }
2984 
2985       gkEnergy += eik;
2986       count++;
2987     }
2988   }
2989 }