View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.numerics.multipole;
39  
40  import ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER;
41  
42  /**
43   * The GeneralizedKirkwoodTensor class contains utilities for generated Generalized Kirkwood
44   * interaction tensors.
45   *
46   * @author Michael J. Schnieders
47   * @since 1.0
48   */
49  public class GKTensorQI extends CoulombTensorQI {
50  
51    /**
52     * The GK tensor can be constructed for monopoles (GB), dipoles or quadrupoles.
53     */
54    protected final GK_MULTIPOLE_ORDER multipoleOrder;
55  
56    /**
57     * The Kirkwood dielectric function for the given multipole order.
58     */
59    private final double c;
60  
61    private final GKSource gkSource;
62  
63    /**
64     * @param multipoleOrder The multipole order.
65     * @param order The tensor order.
66     * @param gkSource Generate the source terms for the GK recurrence.
67     * @param Eh Homogeneous dielectric constant.
68     * @param Es Solvent dielectric constant.
69     */
70    public GKTensorQI(GK_MULTIPOLE_ORDER multipoleOrder, int order, GKSource gkSource, double Eh,
71        double Es) {
72      super(order);
73      this.multipoleOrder = multipoleOrder;
74      this.gkSource = gkSource;
75  
76      // Load the dielectric function
77      c = GKSource.cn(multipoleOrder.getOrder(), Eh, Es);
78    }
79  
80    /**
81     * GK Permanent multipole energy.
82     *
83     * @param mI PolarizableMultipole at site I.
84     * @param mK PolarizableMultipole at site K.
85     * @return the GK permanent multipole energy.
86     */
87    @Override
88    public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
89      return switch (multipoleOrder) {
90        default -> {
91          chargeIPotentialAtK(mI, 2);
92          double eK = multipoleEnergy(mK);
93          chargeKPotentialAtI(mK, 2);
94          double eI = multipoleEnergy(mI);
95          yield c * 0.5 * (eK + eI);
96        }
97        case DIPOLE -> {
98          dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
99          double eK = multipoleEnergy(mK);
100         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
101         double eI = multipoleEnergy(mI);
102         yield c * 0.5 * (eK + eI);
103       }
104       case QUADRUPOLE -> {
105         quadrupoleIPotentialAtK(mI, 2);
106         double eK = multipoleEnergy(mK);
107         quadrupoleKPotentialAtI(mK, 2);
108         double eI = multipoleEnergy(mI);
109         yield c * 0.5 * (eK + eI);
110       }
111     };
112   }
113 
114   /**
115    * GK Permanent multipole energy and gradient.
116    *
117    * @param mI PolarizableMultipole at site I.
118    * @param mK PolarizableMultipole at site K.
119    * @param Gi Coordinate gradient at site I.
120    * @param Gk Coordinate gradient at site K.
121    * @param Ti Torque at site I.
122    * @param Tk Torque at site K.
123    * @return the permanent multipole GK energy.
124    */
125   @Override
126   public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
127       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
128     return switch (multipoleOrder) {
129       default -> monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
130       case DIPOLE -> dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
131       case QUADRUPOLE -> quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
132     };
133   }
134 
135   /**
136    * Permanent multipole energy and gradient using the GK monopole tensor.
137    *
138    * @param mI PolarizableMultipole at site I.
139    * @param mK PolarizableMultipole at site K.
140    * @param Gi Coordinate gradient at site I.
141    * @param Gk Coordinate gradient at site K.
142    * @param Ti Torque at site I.
143    * @param Tk Torque at site K.
144    * @return the permanent multipole GK energy.
145    */
146   protected double monopoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
147       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
148 
149     // Compute the potential due to a multipole component at site I.
150     chargeIPotentialAtK(mI, 3);
151     double eK = multipoleEnergy(mK);
152     multipoleGradient(mK, Gk);
153     multipoleTorque(mK, Tk);
154 
155     // Compute the potential due to a multipole component at site K.
156     chargeKPotentialAtI(mK, 3);
157     double eI = multipoleEnergy(mI);
158     multipoleGradient(mI, Gi);
159     multipoleTorque(mI, Ti);
160 
161     double scale = c * 0.5;
162     Gi[0] = scale * (Gi[0] - Gk[0]);
163     Gi[1] = scale * (Gi[1] - Gk[1]);
164     Gi[2] = scale * (Gi[2] - Gk[2]);
165     Gk[0] = -Gi[0];
166     Gk[1] = -Gi[1];
167     Gk[2] = -Gi[2];
168 
169     Ti[0] = scale * Ti[0];
170     Ti[1] = scale * Ti[1];
171     Ti[2] = scale * Ti[2];
172     Tk[0] = scale * Tk[0];
173     Tk[1] = scale * Tk[1];
174     Tk[2] = scale * Tk[2];
175 
176     return scale * (eK + eI);
177   }
178 
179   /**
180    * Permanent multipole energy and gradient using the GK dipole tensor.
181    *
182    * @param mI PolarizableMultipole at site I.
183    * @param mK PolarizableMultipole at site K.
184    * @param Gi Coordinate gradient at site I.
185    * @param Gk Coordinate gradient at site K.
186    * @param Ti Torque at site I.
187    * @param Tk Torque at site K.
188    * @return the permanent multipole GK energy.
189    */
190   protected double dipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
191       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
192 
193     // Compute the potential due to a multipole component at site I.
194     dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
195     double eK = multipoleEnergy(mK);
196     multipoleGradient(mK, Gk);
197     multipoleTorque(mK, Tk);
198 
199     // Need the torque on site I dipole due to site K multipole.
200     multipoleKPotentialAtI(mK, 1);
201     dipoleTorque(mI, Ti);
202 
203     // Compute the potential due to a multipole component at site K.
204     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
205     double eI = multipoleEnergy(mI);
206     multipoleGradient(mI, Gi);
207     multipoleTorque(mI, Ti);
208 
209     // Need the torque on site K dipole due to site I multipole.
210     multipoleIPotentialAtK(mI, 1);
211     dipoleTorque(mK, Tk);
212 
213     double scale = c * 0.5;
214     Gi[0] = scale * (Gi[0] - Gk[0]);
215     Gi[1] = scale * (Gi[1] - Gk[1]);
216     Gi[2] = scale * (Gi[2] - Gk[2]);
217     Gk[0] = -Gi[0];
218     Gk[1] = -Gi[1];
219     Gk[2] = -Gi[2];
220 
221     Ti[0] = scale * Ti[0];
222     Ti[1] = scale * Ti[1];
223     Ti[2] = scale * Ti[2];
224     Tk[0] = scale * Tk[0];
225     Tk[1] = scale * Tk[1];
226     Tk[2] = scale * Tk[2];
227 
228     return scale * (eK + eI);
229   }
230 
231   /**
232    * Permanent multipole energy and gradient using the GK quadrupole tensor.
233    *
234    * @param mI PolarizableMultipole at site I.
235    * @param mK PolarizableMultipole at site K.
236    * @param Gi Coordinate gradient at site I.
237    * @param Gk Coordinate gradient at site K.
238    * @param Ti Torque at site I.
239    * @param Tk Torque at site K.
240    * @return the permanent multipole GK energy.
241    */
242   protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
243       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
244 
245     // Compute the potential due to a multipole component at site I.
246     quadrupoleIPotentialAtK(mI, 3);
247     double eK = multipoleEnergy(mK);
248     multipoleGradient(mK, Gk);
249     multipoleTorque(mK, Tk);
250 
251     // Need the torque on site I quadrupole due to site K multipole.
252     multipoleKPotentialAtI(mK, 2);
253     quadrupoleTorque(mI, Ti);
254 
255     // Compute the potential due to a multipole component at site K.
256     quadrupoleKPotentialAtI(mK, 3);
257     double eI = multipoleEnergy(mI);
258     multipoleGradient(mI, Gi);
259     multipoleTorque(mI, Ti);
260 
261     // Need the torque on site K quadrupole due to site I multipole.
262     multipoleIPotentialAtK(mI, 2);
263     quadrupoleTorque(mK, Tk);
264 
265     double scale = c * 0.5;
266     Gi[0] = scale * (Gi[0] - Gk[0]);
267     Gi[1] = scale * (Gi[1] - Gk[1]);
268     Gi[2] = scale * (Gi[2] - Gk[2]);
269     Gk[0] = -Gi[0];
270     Gk[1] = -Gi[1];
271     Gk[2] = -Gi[2];
272 
273     Ti[0] = scale * Ti[0];
274     Ti[1] = scale * Ti[1];
275     Ti[2] = scale * Ti[2];
276     Tk[0] = scale * Tk[0];
277     Tk[1] = scale * Tk[1];
278     Tk[2] = scale * Tk[2];
279 
280     return scale * (eK + eI);
281   }
282 
283   /**
284    * GK Permanent multipole Born grad.
285    *
286    * @param mI PolarizableMultipole at site I.
287    * @param mK PolarizableMultipole at site K.
288    * @return a double.
289    */
290   public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
291     generateTensor();
292     return multipoleEnergy(mI, mK);
293   }
294 
295   /**
296    * GK Polarization Energy.
297    *
298    * @param mI PolarizableMultipole at site I.
299    * @param mK PolarizableMultipole at site K.
300    * @param scaleEnergy This is ignored, since masking/scaling is not applied to GK
301    *     interactions.
302    * @return a double.
303    */
304   @Override
305   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
306       double scaleEnergy) {
307     return polarizationEnergy(mI, mK);
308   }
309 
310   /**
311    * GK Polarization Energy.
312    *
313    * @param mI PolarizableMultipole at site I.
314    * @param mK PolarizableMultipole at site K.
315    * @return a double.
316    */
317   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
318     return switch (multipoleOrder) {
319       default -> {
320         // Find the GK charge potential of site I at site K.
321         chargeIPotentialAtK(mI, 1);
322         // Energy of induced dipole K in the field of permanent charge I.
323         double eK = polarizationEnergy(mK);
324         // Find the GK charge potential of site K at site I.
325         chargeKPotentialAtI(mK, 1);
326         // Energy of induced dipole I in the field of permanent charge K.
327         double eI = polarizationEnergy(mI);
328         yield c * 0.5 * (eK + eI);
329       }
330       case DIPOLE -> {
331         // Find the GK dipole potential of site I at site K.
332         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
333         // Energy of induced dipole K in the field of permanent dipole I.
334         double eK = polarizationEnergy(mK);
335         // Find the GK induced dipole potential of site I at site K.
336         dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
337         // Energy of permanent multipole K in the field of induced dipole I.
338         eK += 0.5 * multipoleEnergy(mK);
339         // Find the GK dipole potential of site K at site I.
340         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
341         // Energy of induced dipole I in the field of permanent dipole K.
342         double eI = polarizationEnergy(mI);
343         // Find the GK induced dipole potential of site K at site I.
344         dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
345         // Energy of permanent multipole I in the field of induced dipole K.
346         eI += 0.5 * multipoleEnergy(mI);
347         yield c * 0.5 * (eK + eI);
348       }
349       case QUADRUPOLE -> {
350         // Find the GK quadrupole potential of site I at site K.
351         quadrupoleIPotentialAtK(mI, 1);
352         // Energy of induced dipole K in the field of permanent quadrupole I.
353         double eK = polarizationEnergy(mK);
354         // Find the GK quadrupole potential of site K at site I.
355         quadrupoleKPotentialAtI(mK, 1);
356         // Energy of induced dipole I in the field of permanent quadrupole K.
357         double eI = polarizationEnergy(mI);
358         yield c * 0.5 * (eK + eI);
359       }
360     };
361   }
362 
363   /**
364    * GK Polarization Energy.
365    *
366    * @param mI PolarizableMultipole at site I.
367    * @param mK PolarizableMultipole at site K.
368    * @return a double.
369    */
370   public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
371     return switch (multipoleOrder) {
372       default -> {
373         // Find the GK charge potential of site I at site K.
374         chargeIPotentialAtK(mI, 1);
375         // Energy of induced dipole K in the field of permanent charge I.
376         double eK = polarizationEnergyS(mK);
377         // Find the GK charge potential of site K at site I.
378         chargeKPotentialAtI(mK, 1);
379         // Energy of induced dipole I in the field of permanent charge K.
380         double eI = polarizationEnergyS(mI);
381         yield c * 0.5 * (eK + eI);
382       }
383       case DIPOLE -> {
384         // Find the GK dipole potential of site I at site K.
385         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
386         // Energy of induced dipole K in the field of permanent dipole I.
387         double eK = polarizationEnergyS(mK);
388         // Find the GK induced dipole potential of site I at site K.
389         dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
390         // Energy of permanent multipole K in the field of induced dipole I.
391         eK += 0.5 * multipoleEnergy(mK);
392         // Find the GK dipole potential of site K at site I.
393         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
394         // Energy of induced dipole I in the field of permanent dipole K.
395         double eI = polarizationEnergyS(mI);
396         // Find the GK induced dipole potential of site K at site I.
397         dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
398         // Energy of permanent multipole I in the field of induced dipole K.
399         eI += 0.5 * multipoleEnergy(mI);
400         yield c * 0.5 * (eK + eI);
401       }
402       case QUADRUPOLE -> {
403         // Find the GK quadrupole potential of site I at site K.
404         quadrupoleIPotentialAtK(mI, 1);
405         // Energy of induced dipole K in the field of permanent quadrupole I.
406         double eK = polarizationEnergyS(mK);
407         // Find the GK quadrupole potential of site K at site I.
408         quadrupoleKPotentialAtI(mK, 1);
409         // Energy of induced dipole I in the field of permanent quadrupole K.
410         double eI = polarizationEnergyS(mI);
411         yield c * 0.5 * (eK + eI);
412       }
413     };
414   }
415 
416   /**
417    * Polarization Energy and Gradient.
418    *
419    * @param mI PolarizableMultipole at site I.
420    * @param mK PolarizableMultipole at site K.
421    * @param inductionMask This is ignored, since masking/scaling is not applied to GK
422    *     interactions (everything is intermolecular).
423    * @param energyMask This is ignored, since masking/scaling is not applied to GK interactions
424    *     (everything is intermolecular).
425    * @param mutualMask This should be set to zero for direction polarization.
426    * @param Gi an array of {@link double} objects.
427    * @param Ti an array of {@link double} objects.
428    * @param Tk an array of {@link double} objects.
429    * @return a double.
430    */
431   @Override
432   public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
433       double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
434       double[] Tk) {
435     return switch (multipoleOrder) {
436       default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
437       case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
438       case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
439     };
440   }
441 
442   /**
443    * Monopole Polarization Energy and Gradient.
444    *
445    * @param mI PolarizableMultipole at site I.
446    * @param mK PolarizableMultipole at site K.
447    * @param Gi an array of {@link double} objects.
448    * @return a double.
449    */
450   public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
451       PolarizableMultipole mK, double[] Gi) {
452 
453     // Find the permanent multipole potential at site k.
454     chargeIPotentialAtK(mI, 2);
455     // Energy of induced dipole k in the field of multipole i.
456     double eK = polarizationEnergy(mK);
457     // Derivative with respect to moving atom k.
458     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
459     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
460     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
461 
462     // Find the permanent multipole potential and derivatives at site i.
463     chargeKPotentialAtI(mK, 2);
464     // Energy of induced dipole i in the field of multipole k.
465     double eI = polarizationEnergy(mI);
466     // Derivative with respect to moving atom i.
467     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
468     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
469     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
470 
471     double scale = c * 0.5;
472     Gi[0] *= scale;
473     Gi[1] *= scale;
474     Gi[2] *= scale;
475 
476     // Total polarization energy.
477     return scale * (eI + eK);
478   }
479 
480   /**
481    * Dipole Polarization Energy and Gradient.
482    *
483    * @param mI PolarizableMultipole at site I.
484    * @param mK PolarizableMultipole at site K.
485    * @param mutualMask This should be set to zero for direction polarization.
486    * @param Gi an array of {@link double} objects.
487    * @param Ti an array of {@link double} objects.
488    * @param Tk an array of {@link double} objects.
489    * @return a double.
490    */
491   public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
492       double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
493 
494     // Find the permanent dipole potential and derivatives at site K.
495     dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
496     // Energy of induced dipole k in the field of dipole I.
497     double eK = polarizationEnergy(mK);
498     // Derivative with respect to moving atom K.
499     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
500     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
501     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
502     // Find the potential at K due to the averaged induced dipole at site I.
503     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
504     dipoleTorque(mI, Ti);
505 
506     // Find the GK induced dipole potential of site I at site K.
507     dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
508     // Energy of permanent multipole K in the field of induced dipole I.
509     eK += 0.5 * multipoleEnergy(mK);
510     // Find the GK induced dipole potential of site I at site K.
511     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 3);
512     double[] G = new double[3];
513     multipoleGradient(mK, G);
514     Gi[0] -= G[0];
515     Gi[1] -= G[1];
516     Gi[2] -= G[2];
517     multipoleTorque(mK, Tk);
518 
519     // Find the permanent multipole potential and derivatives at site i.
520     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
521     // Energy of induced dipole i in the field of multipole k.
522     double eI = polarizationEnergy(mI);
523     // Derivative with respect to moving atom i.
524     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
525     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
526     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
527     // Find the potential at I due to the averaged induced dipole at k.
528     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
529     dipoleTorque(mK, Tk);
530 
531     // Find the GK induced dipole potential of site K at site I.
532     dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
533     // Energy of permanent multipole I in the field of induced dipole K.
534     eI += 0.5 * multipoleEnergy(mI);
535     // Find the GK induced dipole potential of site K at site I.
536     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 3);
537     G = new double[3];
538     multipoleGradient(mI, G);
539     Gi[0] += G[0];
540     Gi[1] += G[1];
541     Gi[2] += G[2];
542     multipoleTorque(mI, Ti);
543 
544     // Get the induced-induced portion of the force (Ud . dC/dX . Up).
545     // This contribution does not exist for direct polarization (mutualMask == 0.0).
546     if (mutualMask != 0.0) {
547       // Find the potential and its derivatives at k due to induced dipole i.
548       dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
549       Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
550       Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
551       Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
552 
553       // Find the potential and its derivatives at i due to induced dipole k.
554       dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
555       Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
556       Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
557       Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
558     }
559 
560     // Total polarization energy.
561     double scale = c * 0.5;
562     double energy = scale * (eI + eK);
563     Gi[0] *= scale;
564     Gi[1] *= scale;
565     Gi[2] *= scale;
566     Ti[0] *= scale;
567     Ti[1] *= scale;
568     Ti[2] *= scale;
569     Tk[0] *= scale;
570     Tk[1] *= scale;
571     Tk[2] *= scale;
572 
573     return energy;
574   }
575 
576   /**
577    * Quadrupole Polarization Energy and Gradient.
578    *
579    * @param mI PolarizableMultipole at site I.
580    * @param mK PolarizableMultipole at site K.
581    * @param Gi an array of {@link double} objects.
582    * @param Ti an array of {@link double} objects.
583    * @param Tk an array of {@link double} objects.
584    * @return a double.
585    */
586   public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI,
587       PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {
588 
589     // Find the permanent multipole potential and derivatives at site k.
590     quadrupoleIPotentialAtK(mI, 2);
591     // Energy of induced dipole k in the field of multipole i.
592     double eK = polarizationEnergy(mK);
593     // Derivative with respect to moving atom k.
594     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
595     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
596     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
597 
598     // Find the permanent multipole potential and derivatives at site i.
599     quadrupoleKPotentialAtI(mK, 2);
600     // Energy of induced dipole i in the field of multipole k.
601     double eI = polarizationEnergy(mI);
602     // Derivative with respect to moving atom i.
603     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
604     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
605     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
606 
607     double scale = c * 0.5;
608     Gi[0] *= scale;
609     Gi[1] *= scale;
610     Gi[2] *= scale;
611 
612     // Find the potential and its derivatives at K due to the averaged induced dipole at site i.
613     dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
614     quadrupoleTorque(mK, Tk);
615 
616     // Find the potential and its derivatives at I due to the averaged induced dipole at k.
617     dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
618     quadrupoleTorque(mI, Ti);
619 
620     // Total polarization energy.
621     return scale * (eI + eK);
622   }
623 
624   /**
625    * GK Direct Polarization Born grad.
626    *
627    * @param mI PolarizableMultipole at site I.
628    * @param mK PolarizableMultipole at site K.
629    * @return Partial derivative of the Polarization energy with respect to a Born grad.
630    */
631   public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
632     generateTensor();
633     return 2.0 * polarizationEnergyBorn(mI, mK);
634   }
635 
636   /**
637    * GK Mutual Polarization Contribution to the Born grad.
638    *
639    * @param mI PolarizableMultipole at site I.
640    * @param mK PolarizableMultipole at site K.
641    * @return Mutual Polarization contribution to the partial derivative with respect to a Born grad.
642    */
643   public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
644     double db = 0.0;
645     if (multipoleOrder == GK_MULTIPOLE_ORDER.DIPOLE) {
646       // Find the potential and its derivatives at k due to induced dipole i.
647       dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
648       db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);
649 
650       // Find the potential and its derivatives at i due to induced dipole k.
651       dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
652       db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
653     }
654     return c * db;
655   }
656 
657   /**
658    * Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
659    */
660   @Override
661   protected void source(double[] work) {
662     gkSource.source(work, multipoleOrder);
663   }
664 }