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