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.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 GKTensorGlobal extends CoulombTensorGlobal {
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 number of derivatives to complete.
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 GKTensorGlobal(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 a double.
86     */
87    @Override
88    public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
89      switch (multipoleOrder) {
90        default:
91        case MONOPOLE:
92          chargeIPotentialAtK(mI, 2);
93          double eK = multipoleEnergy(mK);
94          chargeKPotentialAtI(mK, 2);
95          double eI = multipoleEnergy(mI);
96          return c * 0.5 * (eK + eI);
97        case DIPOLE:
98          dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
99          eK = multipoleEnergy(mK);
100         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
101         eI = multipoleEnergy(mI);
102         return c * 0.5 * (eK + eI);
103       case QUADRUPOLE:
104         quadrupoleIPotentialAtK(mI, 2);
105         eK = multipoleEnergy(mK);
106         quadrupoleKPotentialAtI(mK, 2);
107         eI = multipoleEnergy(mI);
108         return c * 0.5 * (eK + eI);
109     }
110   }
111 
112   /**
113    * GK Permanent multipole energy and gradient.
114    *
115    * @param mI PolarizableMultipole at site I.
116    * @param mK PolarizableMultipole at site K.
117    * @param Gi Coordinate gradient at site I.
118    * @param Gk Coordinate gradient at site K.
119    * @param Ti Torque at site I.
120    * @param Tk Torque at site K.
121    * @return the permanent multipole GK energy.
122    */
123   @Override
124   public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
125       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
126     switch (multipoleOrder) {
127       default:
128       case MONOPOLE:
129         return monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
130       case DIPOLE:
131         return dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
132       case QUADRUPOLE:
133         return quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
134     }
135   }
136 
137   /**
138    * Permanent multipole energy and gradient using the GK monopole tensor.
139    *
140    * @param mI PolarizableMultipole at site I.
141    * @param mK PolarizableMultipole at site K.
142    * @param Gi Coordinate gradient at site I.
143    * @param Gk Coordinate gradient at site K.
144    * @param Ti Torque at site I.
145    * @param Tk Torque at site K.
146    * @return the permanent multipole GK energy.
147    */
148   protected double monopoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
149       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
150 
151     // Compute the potential due to a multipole component at site I.
152     chargeIPotentialAtK(mI, 3);
153     double eK = multipoleEnergy(mK);
154     multipoleGradient(mK, Gk);
155     multipoleTorque(mK, Tk);
156 
157     // Compute the potential due to a multipole component at site K.
158     chargeKPotentialAtI(mK, 3);
159     double eI = multipoleEnergy(mI);
160     multipoleGradient(mI, Gi);
161     multipoleTorque(mI, Ti);
162 
163     double scale = c * 0.5;
164     Gi[0] = scale * (Gi[0] - Gk[0]);
165     Gi[1] = scale * (Gi[1] - Gk[1]);
166     Gi[2] = scale * (Gi[2] - Gk[2]);
167     Gk[0] = -Gi[0];
168     Gk[1] = -Gi[1];
169     Gk[2] = -Gi[2];
170 
171     Ti[0] = scale * Ti[0];
172     Ti[1] = scale * Ti[1];
173     Ti[2] = scale * Ti[2];
174     Tk[0] = scale * Tk[0];
175     Tk[1] = scale * Tk[1];
176     Tk[2] = scale * Tk[2];
177 
178     return scale * (eK + eI);
179   }
180 
181   /**
182    * Permanent multipole energy and gradient using the GK dipole tensor.
183    *
184    * @param mI PolarizableMultipole at site I.
185    * @param mK PolarizableMultipole at site K.
186    * @param Gi Coordinate gradient at site I.
187    * @param Gk Coordinate gradient at site K.
188    * @param Ti Torque at site I.
189    * @param Tk Torque at site K.
190    * @return the permanent multipole GK energy.
191    */
192   protected double dipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
193       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
194 
195     // Compute the potential due to a multipole component at site I.
196     dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
197     double eK = multipoleEnergy(mK);
198     multipoleGradient(mK, Gk);
199     multipoleTorque(mK, Tk);
200 
201     // Need the torque on site I pole due to site K multipole.
202     // Only torque on the site I dipole.
203     multipoleKPotentialAtI(mK, 1);
204     dipoleTorque(mI, Ti);
205 
206     // Compute the potential due to a multipole component at site K.
207     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
208     double eI = multipoleEnergy(mI);
209     multipoleGradient(mI, Gi);
210     multipoleTorque(mI, Ti);
211 
212     // Need the torque on site K pole due to site I multipole.
213     // Only torque on the site K dipole.
214     multipoleIPotentialAtK(mI, 1);
215     dipoleTorque(mK, Tk);
216 
217     double scale = c * 0.5;
218     Gi[0] = scale * (Gi[0] - Gk[0]);
219     Gi[1] = scale * (Gi[1] - Gk[1]);
220     Gi[2] = scale * (Gi[2] - Gk[2]);
221     Gk[0] = -Gi[0];
222     Gk[1] = -Gi[1];
223     Gk[2] = -Gi[2];
224 
225     Ti[0] = scale * Ti[0];
226     Ti[1] = scale * Ti[1];
227     Ti[2] = scale * Ti[2];
228     Tk[0] = scale * Tk[0];
229     Tk[1] = scale * Tk[1];
230     Tk[2] = scale * Tk[2];
231 
232     return scale * (eK + eI);
233   }
234 
235   /**
236    * Permanent multipole energy and gradient using the GK quadrupole tensor.
237    *
238    * @param mI PolarizableMultipole at site I.
239    * @param mK PolarizableMultipole at site K.
240    * @param Gi Coordinate gradient at site I.
241    * @param Gk Coordinate gradient at site K.
242    * @param Ti Torque at site I.
243    * @param Tk Torque at site K.
244    * @return the permanent multipole GK energy.
245    */
246   protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
247       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
248 
249     // Compute the potential due to a multipole component at site I.
250     quadrupoleIPotentialAtK(mI, 3);
251     double eK = multipoleEnergy(mK);
252     multipoleGradient(mK, Gk);
253     multipoleTorque(mK, Tk);
254 
255     // Need the torque on site I quadrupole due to site K multipole.
256     multipoleKPotentialAtI(mK, 2);
257     quadrupoleTorque(mI, Ti);
258 
259     // Compute the potential due to a multipole component at site K.
260     quadrupoleKPotentialAtI(mK, 3);
261     double eI = multipoleEnergy(mI);
262     multipoleGradient(mI, Gi);
263     multipoleTorque(mI, Ti);
264 
265     // Need the torque on site K quadrupole due to site I multipole.
266     multipoleIPotentialAtK(mI, 2);
267     quadrupoleTorque(mK, Tk);
268 
269     double scale = c * 0.5;
270     Gi[0] = scale * (Gi[0] - Gk[0]);
271     Gi[1] = scale * (Gi[1] - Gk[1]);
272     Gi[2] = scale * (Gi[2] - Gk[2]);
273     Gk[0] = -Gi[0];
274     Gk[1] = -Gi[1];
275     Gk[2] = -Gi[2];
276 
277     Ti[0] = scale * Ti[0];
278     Ti[1] = scale * Ti[1];
279     Ti[2] = scale * Ti[2];
280     Tk[0] = scale * Tk[0];
281     Tk[1] = scale * Tk[1];
282     Tk[2] = scale * Tk[2];
283 
284     return scale * (eK + eI);
285   }
286 
287   /**
288    * GK Permanent multipole Born grad.
289    *
290    * @param mI PolarizableMultipole at site I.
291    * @param mK PolarizableMultipole at site K.
292    * @return a double.
293    */
294   public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
295     generateTensor();
296     return multipoleEnergy(mI, mK);
297   }
298 
299   /**
300    * GK Polarization Energy.
301    *
302    * @param mI PolarizableMultipole at site I.
303    * @param mK PolarizableMultipole at site K.
304    * @param scaleEnergy This is ignored, since masking/scaling is not applied to GK interactions
305    *     (everything is intermolecular).
306    * @return a double.
307    */
308   @Override
309   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
310       double scaleEnergy) {
311     return polarizationEnergy(mI, mK);
312   }
313 
314   /**
315    * GK Polarization Energy.
316    *
317    * @param mI PolarizableMultipole at site I.
318    * @param mK PolarizableMultipole at site K.
319    * @return a double.
320    */
321   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
322     switch (multipoleOrder) {
323       default:
324       case MONOPOLE:
325         // Find the GK charge potential of site I at site K.
326         chargeIPotentialAtK(mI, 1);
327         // Energy of induced dipole K in the field of permanent charge I.
328         double eK = polarizationEnergy(mK);
329         // Find the GK charge potential of site K at site I.
330         chargeKPotentialAtI(mK, 1);
331         // Energy of induced dipole I in the field of permanent charge K.
332         double eI = polarizationEnergy(mI);
333         return c * 0.5 * (eK + eI);
334       case DIPOLE:
335         // Find the GK dipole potential of site I at site K.
336         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
337         // Energy of induced dipole K in the field of permanent dipole I.
338         eK = polarizationEnergy(mK);
339         // Find the GK induced dipole potential of site I at site K.
340         dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
341         // Energy of permanent multipole K in the field of induced dipole I.
342         eK += 0.5 * multipoleEnergy(mK);
343         // Find the GK dipole potential of site K at site I.
344         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
345         // Energy of induced dipole I in the field of permanent dipole K.
346         eI = polarizationEnergy(mI);
347         // Find the GK induced dipole potential of site K at site I.
348         dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
349         // Energy of permanent multipole I in the field of induced dipole K.
350         eI += 0.5 * multipoleEnergy(mI);
351         return c * 0.5 * (eK + eI);
352       case QUADRUPOLE:
353         // Find the GK quadrupole potential of site I at site K.
354         quadrupoleIPotentialAtK(mI, 1);
355         // Energy of induced dipole K in the field of permanent quadrupole I.
356         eK = polarizationEnergy(mK);
357         // Find the GK quadrupole potential of site K at site I.
358         quadrupoleKPotentialAtI(mK, 1);
359         // Energy of induced dipole I in the field of permanent quadrupole K.
360         eI = polarizationEnergy(mI);
361         return c * 0.5 * (eK + eI);
362     }
363   }
364 
365   /**
366    * GK Polarization Energy.
367    *
368    * @param mI PolarizableMultipole at site I.
369    * @param mK PolarizableMultipole at site K.
370    * @return a double.
371    */
372   public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
373     switch (multipoleOrder) {
374       default:
375       case MONOPOLE:
376         // Find the GK charge potential of site I at site K.
377         chargeIPotentialAtK(mI, 1);
378         // Energy of induced dipole K in the field of permanent charge I.
379         double eK = polarizationEnergyS(mK);
380         // Find the GK charge potential of site K at site I.
381         chargeKPotentialAtI(mK, 1);
382         // Energy of induced dipole I in the field of permanent charge K.
383         double eI = polarizationEnergyS(mI);
384         return c * 0.5 * (eK + eI);
385       case DIPOLE:
386         // Find the GK dipole potential of site I at site K.
387         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
388         // Energy of induced dipole K in the field of permanent dipole I.
389         eK = polarizationEnergyS(mK);
390         // Find the GK induced dipole potential of site I at site K.
391         dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
392         // Energy of permanent multipole K in the field of induced dipole I.
393         eK += 0.5 * multipoleEnergy(mK);
394         // Find the GK dipole potential of site K at site I.
395         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
396         // Energy of induced dipole I in the field of permanent dipole K.
397         eI = polarizationEnergyS(mI);
398         // Find the GK induced dipole potential of site K at site I.
399         dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
400         // Energy of permanent multipole I in the field of induced dipole K.
401         eI += 0.5 * multipoleEnergy(mI);
402         return c * 0.5 * (eK + eI);
403       case QUADRUPOLE:
404         // Find the GK quadrupole potential of site I at site K.
405         quadrupoleIPotentialAtK(mI, 1);
406         // Energy of induced dipole K in the field of permanent quadrupole I.
407         eK = polarizationEnergyS(mK);
408         // Find the GK quadrupole potential of site K at site I.
409         quadrupoleKPotentialAtI(mK, 1);
410         // Energy of induced dipole I in the field of permanent quadrupole K.
411         eI = polarizationEnergyS(mI);
412         return c * 0.5 * (eK + eI);
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,
434       double[] Gi, double[] Ti, double[] Tk) {
435     switch (multipoleOrder) {
436       default:
437       case MONOPOLE:
438         return monopolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
439       case DIPOLE:
440         return dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
441       case QUADRUPOLE:
442         return quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
443     }
444   }
445 
446   /**
447    * Monopole Polarization Energy and Gradient.
448    *
449    * @param mI PolarizableMultipole at site I.
450    * @param mK PolarizableMultipole at site K.
451    * @param Gi an array of {@link double} objects.
452    * @param Ti an array of {@link double} objects.
453    * @param Tk an array of {@link double} objects.
454    * @return a double.
455    */
456   public double monopolePolarizationEnergyAndGradient(
457       PolarizableMultipole mI, PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {
458 
459     // Find the permanent multipole potential at site k.
460     chargeIPotentialAtK(mI, 2);
461     // Energy of induced dipole k in the field of multipole i.
462     double eK = polarizationEnergy(mK);
463     // Derivative with respect to moving atom k.
464     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
465     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
466     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
467 
468     // Find the permanent multipole potential and derivatives at site i.
469     chargeKPotentialAtI(mK, 2);
470     // Energy of induced dipole i in the field of multipole k.
471     double eI = polarizationEnergy(mI);
472     // Derivative with respect to moving atom i.
473     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
474     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
475     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
476 
477     double scale = c * 0.5;
478     Gi[0] *= scale;
479     Gi[1] *= scale;
480     Gi[2] *= scale;
481 
482     // Total polarization energy.
483     return scale * (eI + eK);
484   }
485 
486   /**
487    * Dipole Polarization Energy and Gradient.
488    *
489    * @param mI PolarizableMultipole at site I.
490    * @param mK PolarizableMultipole at site K.
491    * @param mutualMask This should be set to zero for direction polarization.
492    * @param Gi an array of {@link double} objects.
493    * @param Ti an array of {@link double} objects.
494    * @param Tk an array of {@link double} objects.
495    * @return a double.
496    */
497   public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
498       double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
499 
500     // Find the permanent multipole potential and derivatives at site k.
501     dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
502     // Energy of induced dipole k in the field of multipole i.
503     double eK = polarizationEnergy(mK);
504     // Derivative with respect to moving atom k.
505     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
506     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
507     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
508     // Find the potential at K due to the averaged induced dipole at i.
509     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
510     dipoleTorque(mK, Tk);
511 
512     // Find the GK induced dipole potential of site I at site K.
513     dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 3);
514     // Energy of permanent multipole K in the field of induced dipole I.
515     eK += 0.5 * multipoleEnergy(mK);
516     double[] G = new double[3];
517     multipoleGradient(mK, G);
518     Gi[0] -= G[0];
519     Gi[1] -= G[1];
520     Gi[2] -= G[2];
521     multipoleTorque(mK, Tk);
522 
523     // Find the permanent multipole potential and derivatives at site i.
524     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
525     // Energy of induced dipole i in the field of multipole k.
526     double eI = polarizationEnergy(mI);
527     // Derivative with respect to moving atom i.
528     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
529     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
530     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
531     // Find the potential at I due to the averaged induced dipole at k.
532     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
533     dipoleTorque(mI, Ti);
534 
535     // Find the GK induced dipole potential of site K at site I.
536     dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 3);
537     // Energy of permanent multipole I in the field of induced dipole K.
538     eI += 0.5 * multipoleEnergy(mI);
539     G = new double[3];
540     multipoleGradient(mI, G);
541     Gi[0] += G[0];
542     Gi[1] += G[1];
543     Gi[2] += G[2];
544     multipoleTorque(mI, Ti);
545 
546     // Get the induced-induced portion of the force (Ud . dC/dX . Up).
547     // This contribution does not exist for direct polarization (mutualMask == 0.0).
548     if (mutualMask != 0.0) {
549       // Find the potential and its derivatives at k due to induced dipole i.
550       dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
551       Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
552       Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
553       Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
554 
555       // Find the potential and its derivatives at i due to induced dipole k.
556       dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
557       Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
558       Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
559       Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
560     }
561 
562     // Total polarization energy.
563     double scale = c * 0.5;
564     double energy = scale * (eI + eK);
565     Gi[0] *= scale;
566     Gi[1] *= scale;
567     Gi[2] *= scale;
568     Ti[0] *= scale;
569     Ti[1] *= scale;
570     Ti[2] *= scale;
571     Tk[0] *= scale;
572     Tk[1] *= scale;
573     Tk[2] *= scale;
574 
575     return energy;
576   }
577 
578   /**
579    * Quadrupole Polarization Energy and Gradient.
580    *
581    * @param mI PolarizableMultipole at site I.
582    * @param mK PolarizableMultipole at site K.
583    * @param Gi an array of {@link double} objects.
584    * @param Ti an array of {@link double} objects.
585    * @param Tk an array of {@link double} objects.
586    * @return a double.
587    */
588   public double quadrupolePolarizationEnergyAndGradient(
589       PolarizableMultipole mI, PolarizableMultipole mK,
590       double[] Gi, double[] Ti, double[] Tk) {
591 
592     // Find the permanent multipole potential and derivatives at site k.
593     quadrupoleIPotentialAtK(mI, 2);
594     // Energy of induced dipole k in the field of multipole i.
595     double eK = polarizationEnergy(mK);
596     // Derivative with respect to moving atom k.
597     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
598     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
599     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
600 
601     // Find the permanent multipole potential and derivatives at site i.
602     quadrupoleKPotentialAtI(mK, 2);
603     // Energy of induced dipole i in the field of multipole k.
604     double eI = polarizationEnergy(mI);
605     // Derivative with respect to moving atom i.
606     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
607     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
608     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
609 
610     double scale = c * 0.5;
611     Gi[0] *= scale;
612     Gi[1] *= scale;
613     Gi[2] *= scale;
614 
615     // Find the potential and its derivatives at K due to the averaged induced dipole at i.
616     dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
617     quadrupoleTorque(mK, Tk);
618 
619     // Find the potential and its derivatives at I due to the averaged induced dipole at k.
620     dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
621     quadrupoleTorque(mI, Ti);
622 
623     // Total polarization energy.
624     return scale * (eI + eK);
625   }
626 
627   /**
628    * GK Direct Polarization Born grad.
629    *
630    * @param mI PolarizableMultipole at site I.
631    * @param mK PolarizableMultipole at site K.
632    * @return Partial derivative of the Polarization energy with respect to a Born grad.
633    */
634   public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
635     generateTensor();
636     return 2.0 * polarizationEnergyBorn(mI, mK);
637   }
638 
639   /**
640    * GK Mutual Polarization Contribution to the Born grad.
641    *
642    * @param mI PolarizableMultipole at site I.
643    * @param mK PolarizableMultipole at site K.
644    * @return Mutual Polarization contribution to the partial derivative with respect to a Born grad.
645    */
646   public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
647     double db = 0.0;
648     if (multipoleOrder == GK_MULTIPOLE_ORDER.DIPOLE) {
649       // Find the potential and its derivatives at k due to induced dipole i.
650       dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
651       db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);
652 
653       // Find the potential and its derivatives at i due to induced dipole k.
654       dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
655       db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
656     }
657     return c * db;
658   }
659 
660   /**
661    * Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
662    */
663   @Override
664   protected void source(double[] work) {
665     gkSource.source(work, multipoleOrder);
666   }
667 
668 }