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 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      switch (multipoleOrder.getOrder()) {
90        default:
91        case 0:
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 1:
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 2:
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 dipole due to site K multipole.
202     multipoleKPotentialAtI(mK, 1);
203     dipoleTorque(mI, Ti);
204 
205     // Compute the potential due to a multipole component at site K.
206     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
207     double eI = multipoleEnergy(mI);
208     multipoleGradient(mI, Gi);
209     multipoleTorque(mI, Ti);
210 
211     // Need the torque on site K dipole due to site I multipole.
212     multipoleIPotentialAtK(mI, 1);
213     dipoleTorque(mK, Tk);
214 
215     double scale = c * 0.5;
216     Gi[0] = scale * (Gi[0] - Gk[0]);
217     Gi[1] = scale * (Gi[1] - Gk[1]);
218     Gi[2] = scale * (Gi[2] - Gk[2]);
219     Gk[0] = -Gi[0];
220     Gk[1] = -Gi[1];
221     Gk[2] = -Gi[2];
222 
223     Ti[0] = scale * Ti[0];
224     Ti[1] = scale * Ti[1];
225     Ti[2] = scale * Ti[2];
226     Tk[0] = scale * Tk[0];
227     Tk[1] = scale * Tk[1];
228     Tk[2] = scale * Tk[2];
229 
230     return scale * (eK + eI);
231   }
232 
233   /**
234    * Permanent multipole energy and gradient using the GK quadrupole tensor.
235    *
236    * @param mI PolarizableMultipole at site I.
237    * @param mK PolarizableMultipole at site K.
238    * @param Gi Coordinate gradient at site I.
239    * @param Gk Coordinate gradient at site K.
240    * @param Ti Torque at site I.
241    * @param Tk Torque at site K.
242    * @return the permanent multipole GK energy.
243    */
244   protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
245       double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
246 
247     // Compute the potential due to a multipole component at site I.
248     quadrupoleIPotentialAtK(mI, 3);
249     double eK = multipoleEnergy(mK);
250     multipoleGradient(mK, Gk);
251     multipoleTorque(mK, Tk);
252 
253     // Need the torque on site I quadrupole due to site K multipole.
254     multipoleKPotentialAtI(mK, 2);
255     quadrupoleTorque(mI, Ti);
256 
257     // Compute the potential due to a multipole component at site K.
258     quadrupoleKPotentialAtI(mK, 3);
259     double eI = multipoleEnergy(mI);
260     multipoleGradient(mI, Gi);
261     multipoleTorque(mI, Ti);
262 
263     // Need the torque on site K quadrupole due to site I multipole.
264     multipoleIPotentialAtK(mI, 2);
265     quadrupoleTorque(mK, Tk);
266 
267     double scale = c * 0.5;
268     Gi[0] = scale * (Gi[0] - Gk[0]);
269     Gi[1] = scale * (Gi[1] - Gk[1]);
270     Gi[2] = scale * (Gi[2] - Gk[2]);
271     Gk[0] = -Gi[0];
272     Gk[1] = -Gi[1];
273     Gk[2] = -Gi[2];
274 
275     Ti[0] = scale * Ti[0];
276     Ti[1] = scale * Ti[1];
277     Ti[2] = scale * Ti[2];
278     Tk[0] = scale * Tk[0];
279     Tk[1] = scale * Tk[1];
280     Tk[2] = scale * Tk[2];
281 
282     return scale * (eK + eI);
283   }
284 
285   /**
286    * GK Permanent multipole Born grad.
287    *
288    * @param mI PolarizableMultipole at site I.
289    * @param mK PolarizableMultipole at site K.
290    * @return a double.
291    */
292   public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
293     generateTensor();
294     return multipoleEnergy(mI, mK);
295   }
296 
297   /**
298    * GK Polarization Energy.
299    *
300    * @param mI PolarizableMultipole at site I.
301    * @param mK PolarizableMultipole at site K.
302    * @param scaleEnergy This is ignored, since masking/scaling is not applied to GK
303    *     interactions.
304    * @return a double.
305    */
306   @Override
307   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
308       double scaleEnergy) {
309     return polarizationEnergy(mI, mK);
310   }
311 
312   /**
313    * GK Polarization Energy.
314    *
315    * @param mI PolarizableMultipole at site I.
316    * @param mK PolarizableMultipole at site K.
317    * @return a double.
318    */
319   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
320     switch (multipoleOrder) {
321       default:
322       case MONOPOLE:
323         // Find the GK charge potential of site I at site K.
324         chargeIPotentialAtK(mI, 1);
325         // Energy of induced dipole K in the field of permanent charge I.
326         double eK = polarizationEnergy(mK);
327         // Find the GK charge potential of site K at site I.
328         chargeKPotentialAtI(mK, 1);
329         // Energy of induced dipole I in the field of permanent charge K.
330         double eI = polarizationEnergy(mI);
331         return c * 0.5 * (eK + eI);
332       case DIPOLE:
333         // Find the GK dipole potential of site I at site K.
334         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
335         // Energy of induced dipole K in the field of permanent dipole I.
336         eK = polarizationEnergy(mK);
337         // Find the GK induced dipole potential of site I at site K.
338         dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
339         // Energy of permanent multipole K in the field of induced dipole I.
340         eK += 0.5 * multipoleEnergy(mK);
341         // Find the GK dipole potential of site K at site I.
342         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
343         // Energy of induced dipole I in the field of permanent dipole K.
344         eI = polarizationEnergy(mI);
345         // Find the GK induced dipole potential of site K at site I.
346         dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
347         // Energy of permanent multipole I in the field of induced dipole K.
348         eI += 0.5 * multipoleEnergy(mI);
349         return c * 0.5 * (eK + eI);
350       case QUADRUPOLE:
351         // Find the GK quadrupole potential of site I at site K.
352         quadrupoleIPotentialAtK(mI, 1);
353         // Energy of induced dipole K in the field of permanent quadrupole I.
354         eK = polarizationEnergy(mK);
355         // Find the GK quadrupole potential of site K at site I.
356         quadrupoleKPotentialAtI(mK, 1);
357         // Energy of induced dipole I in the field of permanent quadrupole K.
358         eI = polarizationEnergy(mI);
359         return c * 0.5 * (eK + eI);
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     switch (multipoleOrder) {
372       default:
373       case MONOPOLE:
374         // Find the GK charge potential of site I at site K.
375         chargeIPotentialAtK(mI, 1);
376         // Energy of induced dipole K in the field of permanent charge I.
377         double eK = polarizationEnergyS(mK);
378         // Find the GK charge potential of site K at site I.
379         chargeKPotentialAtI(mK, 1);
380         // Energy of induced dipole I in the field of permanent charge K.
381         double eI = polarizationEnergyS(mI);
382         return c * 0.5 * (eK + eI);
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         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         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         return c * 0.5 * (eK + eI);
401       case QUADRUPOLE:
402         // Find the GK quadrupole potential of site I at site K.
403         quadrupoleIPotentialAtK(mI, 1);
404         // Energy of induced dipole K in the field of permanent quadrupole I.
405         eK = polarizationEnergyS(mK);
406         // Find the GK quadrupole potential of site K at site I.
407         quadrupoleKPotentialAtI(mK, 1);
408         // Energy of induced dipole I in the field of permanent quadrupole K.
409         eI = polarizationEnergyS(mI);
410         return c * 0.5 * (eK + eI);
411     }
412   }
413 
414   /**
415    * Polarization Energy and Gradient.
416    *
417    * @param mI PolarizableMultipole at site I.
418    * @param mK PolarizableMultipole at site K.
419    * @param inductionMask This is ignored, since masking/scaling is not applied to GK
420    *     interactions (everything is intermolecular).
421    * @param energyMask This is ignored, since masking/scaling is not applied to GK interactions
422    *     (everything is intermolecular).
423    * @param mutualMask This should be set to zero for direction polarization.
424    * @param Gi an array of {@link double} objects.
425    * @param Ti an array of {@link double} objects.
426    * @param Tk an array of {@link double} objects.
427    * @return a double.
428    */
429   @Override
430   public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
431       double inductionMask, double energyMask, double mutualMask,
432       double[] Gi, double[] Ti, double[] Tk) {
433     switch (multipoleOrder) {
434       default:
435       case MONOPOLE:
436         return monopolePolarizationEnergyAndGradient(mI, mK, Gi);
437       case DIPOLE:
438         return dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
439       case QUADRUPOLE:
440         return quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
441     }
442   }
443 
444   /**
445    * Monopole Polarization Energy and Gradient.
446    *
447    * @param mI PolarizableMultipole at site I.
448    * @param mK PolarizableMultipole at site K.
449    * @param Gi an array of {@link double} objects.
450    * @return a double.
451    */
452   public double monopolePolarizationEnergyAndGradient(
453       PolarizableMultipole mI, PolarizableMultipole mK, double[] Gi) {
454 
455     // Find the permanent multipole potential at site k.
456     chargeIPotentialAtK(mI, 2);
457     // Energy of induced dipole k in the field of multipole i.
458     double eK = polarizationEnergy(mK);
459     // Derivative with respect to moving atom k.
460     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
461     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
462     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
463 
464     // Find the permanent multipole potential and derivatives at site i.
465     chargeKPotentialAtI(mK, 2);
466     // Energy of induced dipole i in the field of multipole k.
467     double eI = polarizationEnergy(mI);
468     // Derivative with respect to moving atom i.
469     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
470     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
471     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
472 
473     double scale = c * 0.5;
474     Gi[0] *= scale;
475     Gi[1] *= scale;
476     Gi[2] *= scale;
477 
478     // Total polarization energy.
479     return scale * (eI + eK);
480   }
481 
482   /**
483    * Dipole Polarization Energy and Gradient.
484    *
485    * @param mI PolarizableMultipole at site I.
486    * @param mK PolarizableMultipole at site K.
487    * @param mutualMask This should be set to zero for direction polarization.
488    * @param Gi an array of {@link double} objects.
489    * @param Ti an array of {@link double} objects.
490    * @param Tk an array of {@link double} objects.
491    * @return a double.
492    */
493   public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
494       double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
495 
496     // Find the permanent dipole potential and derivatives at site K.
497     dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
498     // Energy of induced dipole k in the field of dipole I.
499     double eK = polarizationEnergy(mK);
500     // Derivative with respect to moving atom K.
501     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
502     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
503     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
504     // Find the potential at K due to the averaged induced dipole at I.
505     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
506     dipoleTorque(mI, Ti);
507 
508     // Find the GK induced dipole potential of site I at site K.
509     dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
510     // Energy of permanent multipole K in the field of induced dipole I.
511     eK += 0.5 * multipoleEnergy(mK);
512     // Find the GK induced dipole potential of site I at site K.
513     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 3);
514     double[] G = new double[3];
515     multipoleGradient(mK, G);
516     Gi[0] -= G[0];
517     Gi[1] -= G[1];
518     Gi[2] -= G[2];
519     multipoleTorque(mK, Tk);
520 
521     // Find the permanent multipole potential and derivatives at site i.
522     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
523     // Energy of induced dipole i in the field of multipole k.
524     double eI = polarizationEnergy(mI);
525     // Derivative with respect to moving atom i.
526     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
527     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
528     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
529     // Find the potential at I due to the averaged induced dipole at k.
530     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
531     dipoleTorque(mK, Tk);
532 
533     // Find the GK induced dipole potential of site K at site I.
534     dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
535     // Energy of permanent multipole I in the field of induced dipole K.
536     eI += 0.5 * multipoleEnergy(mI);
537     // Find the GK induced dipole potential of site K at site I.
538     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 3);
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 }