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-2025.
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  /**
41   * The GeneralizedKirkwoodTensor class contains utilities for generated Generalized Kirkwood
42   * interaction tensors.
43   *
44   * @author Michael J. Schnieders
45   * @since 1.0
46   */
47  public class GKTensorQI extends CoulombTensorQI {
48  
49    /**
50     * The GK tensor can be constructed for monopoles (GB), dipoles or quadrupoles.
51     */
52    protected final GKMultipoleOrder multipoleOrder;
53  
54    /**
55     * The Kirkwood dielectric function for the given multipole order.
56     */
57    private final double c;
58  
59    private final GKSource gkSource;
60  
61    /**
62     * Create a new GKTensorQI object.
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(GKMultipoleOrder 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     return multipoleEnergy(mI, mK);
292   }
293 
294   /**
295    * GK Polarization Energy.
296    *
297    * @param mI          PolarizableMultipole at site I.
298    * @param mK          PolarizableMultipole at site K.
299    * @param scaleEnergy This is ignored, since masking/scaling is not applied to GK
300    *                    interactions.
301    * @return a double.
302    */
303   @Override
304   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
305                                    double scaleEnergy) {
306     return polarizationEnergy(mI, mK);
307   }
308 
309   /**
310    * GK Polarization Energy.
311    *
312    * @param mI PolarizableMultipole at site I.
313    * @param mK PolarizableMultipole at site K.
314    * @return a double.
315    */
316   public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
317     return switch (multipoleOrder) {
318       default -> {
319         // Find the GK charge potential of site I at site K.
320         chargeIPotentialAtK(mI, 1);
321         // Energy of induced dipole K in the field of permanent charge I.
322         double eK = polarizationEnergy(mK);
323         // Find the GK charge potential of site K at site I.
324         chargeKPotentialAtI(mK, 1);
325         // Energy of induced dipole I in the field of permanent charge K.
326         double eI = polarizationEnergy(mI);
327         yield c * 0.5 * (eK + eI);
328       }
329       case DIPOLE -> {
330         // Find the GK dipole potential of site I at site K.
331         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
332         // Energy of induced dipole K in the field of permanent dipole I.
333         double eK = polarizationEnergy(mK);
334         // Find the GK induced dipole potential of site I at site K.
335         dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
336         // Energy of permanent multipole K in the field of induced dipole I.
337         eK += 0.5 * multipoleEnergy(mK);
338         // Find the GK dipole potential of site K at site I.
339         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
340         // Energy of induced dipole I in the field of permanent dipole K.
341         double eI = polarizationEnergy(mI);
342         // Find the GK induced dipole potential of site K at site I.
343         dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
344         // Energy of permanent multipole I in the field of induced dipole K.
345         eI += 0.5 * multipoleEnergy(mI);
346         yield c * 0.5 * (eK + eI);
347       }
348       case QUADRUPOLE -> {
349         // Find the GK quadrupole potential of site I at site K.
350         quadrupoleIPotentialAtK(mI, 1);
351         // Energy of induced dipole K in the field of permanent quadrupole I.
352         double eK = polarizationEnergy(mK);
353         // Find the GK quadrupole potential of site K at site I.
354         quadrupoleKPotentialAtI(mK, 1);
355         // Energy of induced dipole I in the field of permanent quadrupole K.
356         double eI = polarizationEnergy(mI);
357         yield c * 0.5 * (eK + eI);
358       }
359     };
360   }
361 
362   /**
363    * GK Polarization Energy.
364    *
365    * @param mI PolarizableMultipole at site I.
366    * @param mK PolarizableMultipole at site K.
367    * @return a double.
368    */
369   public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
370     return switch (multipoleOrder) {
371       default -> {
372         // Find the GK charge potential of site I at site K.
373         chargeIPotentialAtK(mI, 1);
374         // Energy of induced dipole K in the field of permanent charge I.
375         double eK = polarizationEnergyS(mK);
376         // Find the GK charge potential of site K at site I.
377         chargeKPotentialAtI(mK, 1);
378         // Energy of induced dipole I in the field of permanent charge K.
379         double eI = polarizationEnergyS(mI);
380         yield c * 0.5 * (eK + eI);
381       }
382       case DIPOLE -> {
383         // Find the GK dipole potential of site I at site K.
384         dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
385         // Energy of induced dipole K in the field of permanent dipole I.
386         double eK = polarizationEnergyS(mK);
387         // Find the GK induced dipole potential of site I at site K.
388         dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
389         // Energy of permanent multipole K in the field of induced dipole I.
390         eK += 0.5 * multipoleEnergy(mK);
391         // Find the GK dipole potential of site K at site I.
392         dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
393         // Energy of induced dipole I in the field of permanent dipole K.
394         double eI = polarizationEnergyS(mI);
395         // Find the GK induced dipole potential of site K at site I.
396         dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
397         // Energy of permanent multipole I in the field of induced dipole K.
398         eI += 0.5 * multipoleEnergy(mI);
399         yield c * 0.5 * (eK + eI);
400       }
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         double 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         double eI = polarizationEnergyS(mI);
410         yield c * 0.5 * (eK + eI);
411       }
412     };
413   }
414 
415   /**
416    * Polarization Energy and Gradient.
417    *
418    * @param mI            PolarizableMultipole at site I.
419    * @param mK            PolarizableMultipole at site K.
420    * @param inductionMask This is ignored, since masking/scaling is not applied to GK
421    *                      interactions (everything is intermolecular).
422    * @param energyMask    This is ignored, since masking/scaling is not applied to GK interactions
423    *                      (everything is intermolecular).
424    * @param mutualMask    This should be set to zero for direction polarization.
425    * @param Gi            an array of {@link double} objects.
426    * @param Ti            an array of {@link double} objects.
427    * @param Tk            an array of {@link double} objects.
428    * @return a double.
429    */
430   @Override
431   public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
432                                               double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
433                                               double[] Tk) {
434     return switch (multipoleOrder) {
435       default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
436       case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
437       case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
438     };
439   }
440 
441   /**
442    * Monopole Polarization Energy and Gradient.
443    *
444    * @param mI PolarizableMultipole at site I.
445    * @param mK PolarizableMultipole at site K.
446    * @param Gi an array of {@link double} objects.
447    * @return a double.
448    */
449   public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
450                                                       PolarizableMultipole mK, double[] Gi) {
451 
452     // Find the permanent multipole potential at site k.
453     chargeIPotentialAtK(mI, 2);
454     // Energy of induced dipole k in the field of multipole i.
455     double eK = polarizationEnergy(mK);
456     // Derivative with respect to moving atom k.
457     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
458     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
459     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
460 
461     // Find the permanent multipole potential and derivatives at site i.
462     chargeKPotentialAtI(mK, 2);
463     // Energy of induced dipole i in the field of multipole k.
464     double eI = polarizationEnergy(mI);
465     // Derivative with respect to moving atom i.
466     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
467     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
468     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
469 
470     double scale = c * 0.5;
471     Gi[0] *= scale;
472     Gi[1] *= scale;
473     Gi[2] *= scale;
474 
475     // Total polarization energy.
476     return scale * (eI + eK);
477   }
478 
479   /**
480    * Dipole Polarization Energy and Gradient.
481    *
482    * @param mI         PolarizableMultipole at site I.
483    * @param mK         PolarizableMultipole at site K.
484    * @param mutualMask This should be set to zero for direction polarization.
485    * @param Gi         an array of {@link double} objects.
486    * @param Ti         an array of {@link double} objects.
487    * @param Tk         an array of {@link double} objects.
488    * @return a double.
489    */
490   public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
491                                                     double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
492 
493     // Find the permanent dipole potential and derivatives at site K.
494     dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
495     // Energy of induced dipole k in the field of dipole I.
496     double eK = polarizationEnergy(mK);
497     // Derivative with respect to moving atom K.
498     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
499     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
500     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
501     // Find the potential at K due to the averaged induced dipole at site I.
502     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
503     dipoleTorque(mI, Ti);
504 
505     // Find the GK induced dipole potential of site I at site K.
506     dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
507     // Energy of permanent multipole K in the field of induced dipole I.
508     eK += 0.5 * multipoleEnergy(mK);
509     // Find the GK induced dipole potential of site I at site K.
510     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 3);
511     double[] G = new double[3];
512     multipoleGradient(mK, G);
513     Gi[0] -= G[0];
514     Gi[1] -= G[1];
515     Gi[2] -= G[2];
516     multipoleTorque(mK, Tk);
517 
518     // Find the permanent multipole potential and derivatives at site i.
519     dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
520     // Energy of induced dipole i in the field of multipole k.
521     double eI = polarizationEnergy(mI);
522     // Derivative with respect to moving atom i.
523     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
524     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
525     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
526     // Find the potential at I due to the averaged induced dipole at k.
527     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
528     dipoleTorque(mK, Tk);
529 
530     // Find the GK induced dipole potential of site K at site I.
531     dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
532     // Energy of permanent multipole I in the field of induced dipole K.
533     eI += 0.5 * multipoleEnergy(mI);
534     // Find the GK induced dipole potential of site K at site I.
535     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 3);
536     G = new double[3];
537     multipoleGradient(mI, G);
538     Gi[0] += G[0];
539     Gi[1] += G[1];
540     Gi[2] += G[2];
541     multipoleTorque(mI, Ti);
542 
543     // Get the induced-induced portion of the force (Ud . dC/dX . Up).
544     // This contribution does not exist for direct polarization (mutualMask == 0.0).
545     if (mutualMask != 0.0) {
546       // Find the potential and its derivatives at k due to induced dipole i.
547       dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
548       Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
549       Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
550       Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
551 
552       // Find the potential and its derivatives at i due to induced dipole k.
553       dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
554       Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
555       Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
556       Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
557     }
558 
559     // Total polarization energy.
560     double scale = c * 0.5;
561     double energy = scale * (eI + eK);
562     Gi[0] *= scale;
563     Gi[1] *= scale;
564     Gi[2] *= scale;
565     Ti[0] *= scale;
566     Ti[1] *= scale;
567     Ti[2] *= scale;
568     Tk[0] *= scale;
569     Tk[1] *= scale;
570     Tk[2] *= scale;
571 
572     return energy;
573   }
574 
575   /**
576    * Quadrupole Polarization Energy and Gradient.
577    *
578    * @param mI PolarizableMultipole at site I.
579    * @param mK PolarizableMultipole at site K.
580    * @param Gi an array of {@link double} objects.
581    * @param Ti an array of {@link double} objects.
582    * @param Tk an array of {@link double} objects.
583    * @return a double.
584    */
585   public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI,
586                                                         PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {
587 
588     // Find the permanent multipole potential and derivatives at site k.
589     quadrupoleIPotentialAtK(mI, 2);
590     // Energy of induced dipole k in the field of multipole i.
591     double eK = polarizationEnergy(mK);
592     // Derivative with respect to moving atom k.
593     Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
594     Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
595     Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
596 
597     // Find the permanent multipole potential and derivatives at site i.
598     quadrupoleKPotentialAtI(mK, 2);
599     // Energy of induced dipole i in the field of multipole k.
600     double eI = polarizationEnergy(mI);
601     // Derivative with respect to moving atom i.
602     Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
603     Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
604     Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
605 
606     double scale = c * 0.5;
607     Gi[0] *= scale;
608     Gi[1] *= scale;
609     Gi[2] *= scale;
610 
611     // Find the potential and its derivatives at K due to the averaged induced dipole at site i.
612     dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
613     quadrupoleTorque(mK, Tk);
614 
615     // Find the potential and its derivatives at I due to the averaged induced dipole at k.
616     dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
617     quadrupoleTorque(mI, Ti);
618 
619     // Total polarization energy.
620     return scale * (eI + eK);
621   }
622 
623   /**
624    * GK Direct Polarization Born grad.
625    *
626    * @param mI PolarizableMultipole at site I.
627    * @param mK PolarizableMultipole at site K.
628    * @return Partial derivative of the Polarization energy with respect to a Born grad.
629    */
630   public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
631     return 2.0 * polarizationEnergyBorn(mI, mK);
632   }
633 
634   /**
635    * GK Mutual Polarization Contribution to the Born grad.
636    *
637    * @param mI PolarizableMultipole at site I.
638    * @param mK PolarizableMultipole at site K.
639    * @return Mutual Polarization contribution to the partial derivative with respect to a Born grad.
640    */
641   public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
642     double db = 0.0;
643     if (multipoleOrder == GKMultipoleOrder.DIPOLE) {
644       // Find the potential and its derivatives at k due to induced dipole i.
645       dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
646       db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);
647 
648       // Find the potential and its derivatives at i due to induced dipole k.
649       dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
650       db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
651     }
652     return c * db;
653   }
654 
655   /**
656    * Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
657    */
658   @Override
659   protected void source(double[] work) {
660     gkSource.source(work, multipoleOrder);
661   }
662 }