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