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  import jdk.incubator.vector.DoubleVector;
41  
42  import java.util.logging.Level;
43  import java.util.logging.Logger;
44  
45  import static ffx.numerics.math.ScalarMath.doubleFactorial;
46  import static org.apache.commons.math3.util.FastMath.pow;
47  
48  /**
49   * The abstract MultipoleTensorSIMD is extended by classes that compute derivatives of 1/|<b>r</b>| via
50   * recursion to arbitrary order using Cartesian multipoles in either a global frame or a
51   * quasi-internal frame.
52   * <br>
53   * This class serves as the abstract parent to both and defines all shared logic. Non-abstract methods
54   * are declared final to disallow unnecessary overrides.
55   *
56   * @author Michael J. Schnieders
57   * @since 1.0
58   */
59  public abstract class MultipoleTensorSIMD {
60  
61    /**
62     * Logger for the MultipoleTensor class.
63     */
64    private static final Logger logger = Logger.getLogger(MultipoleTensorSIMD.class.getName());
65  
66    /**
67     * Order of the tensor recursion (5th is needed for AMOEBA forces).
68     */
69    protected final int order;
70  
71    /**
72     * Order plus 1.
73     */
74    protected final int o1;
75  
76    /**
77     * Order plus one.
78     */
79    protected final int il;
80  
81    /**
82     * im = (Order plus one)^2.
83     */
84    protected final int im;
85  
86    /**
87     * in = (Order plus one)^3.
88     */
89    protected final int in;
90  
91    /**
92     * Size = (order + 1) * (order + 2) * (order + 3) / 6;
93     */
94    protected final int size;
95  
96    /**
97     * The OPERATOR in use.
98     */
99    protected Operator operator;
100 
101   /**
102    * These are the "source" terms for the recursion for the Coulomb operator (1/R).
103    */
104   protected final double[] coulombSource;
105 
106   /**
107    * The coordinate system in use (global or QI).
108    */
109   protected final CoordinateSystem coordinates;
110 
111   /**
112    * Separation distance.
113    */
114   protected DoubleVector R;
115 
116   /**
117    * Separation distance squared.
118    */
119   protected DoubleVector r2;
120 
121   /**
122    * Xk - Xi.
123    */
124   protected DoubleVector x;
125 
126   /**
127    * Yk - Yi.
128    */
129   protected DoubleVector y;
130 
131   /**
132    * Zk - Zi.
133    */
134   protected DoubleVector z;
135 
136   /**
137    * A work array with auxiliary terms for the recursion.
138    */
139   protected final DoubleVector[] work;
140 
141   // Cartesian tensor elements (for 1/R, erfc(Beta*R)/R or DoubleVectorhole damping.
142   protected DoubleVector R000;
143   // l + m + n = 1 (3)   4
144   protected DoubleVector R100;
145   protected DoubleVector R010;
146   protected DoubleVector R001;
147   // l + m + n = 2 (6)  10
148   protected DoubleVector R200;
149   protected DoubleVector R020;
150   protected DoubleVector R002;
151   protected DoubleVector R110;
152   protected DoubleVector R101;
153   protected DoubleVector R011;
154   // l + m + n = 3 (10) 20
155   protected DoubleVector R300;
156   protected DoubleVector R030;
157   protected DoubleVector R003;
158   protected DoubleVector R210;
159   protected DoubleVector R201;
160   protected DoubleVector R120;
161   protected DoubleVector R021;
162   protected DoubleVector R102;
163   protected DoubleVector R012;
164   protected DoubleVector R111;
165   // l + m + n = 4 (15) 35
166   protected DoubleVector R400;
167   protected DoubleVector R040;
168   protected DoubleVector R004;
169   protected DoubleVector R310;
170   protected DoubleVector R301;
171   protected DoubleVector R130;
172   protected DoubleVector R031;
173   protected DoubleVector R103;
174   protected DoubleVector R013;
175   protected DoubleVector R220;
176   protected DoubleVector R202;
177   protected DoubleVector R022;
178   protected DoubleVector R211;
179   protected DoubleVector R121;
180   protected DoubleVector R112;
181   // l + m + n = 5 (21) 56
182   protected DoubleVector R500;
183   protected DoubleVector R050;
184   protected DoubleVector R005;
185   protected DoubleVector R410;
186   protected DoubleVector R401;
187   protected DoubleVector R140;
188   protected DoubleVector R041;
189   protected DoubleVector R104;
190   protected DoubleVector R014;
191   protected DoubleVector R320;
192   protected DoubleVector R302;
193   protected DoubleVector R230;
194   protected DoubleVector R032;
195   protected DoubleVector R203;
196   protected DoubleVector R023;
197   protected DoubleVector R311;
198   protected DoubleVector R131;
199   protected DoubleVector R113;
200   protected DoubleVector R221;
201   protected DoubleVector R212;
202   protected DoubleVector R122;
203   // l + m + n = 6 (28) 84
204   protected DoubleVector R006;
205   protected DoubleVector R402;
206   protected DoubleVector R042;
207   protected DoubleVector R204;
208   protected DoubleVector R024;
209   protected DoubleVector R222;
210   protected DoubleVector R600;
211   protected DoubleVector R060;
212   protected DoubleVector R510;
213   protected DoubleVector R501;
214   protected DoubleVector R150;
215   protected DoubleVector R051;
216   protected DoubleVector R105;
217   protected DoubleVector R015;
218   protected DoubleVector R420;
219   protected DoubleVector R240;
220   protected DoubleVector R411;
221   protected DoubleVector R141;
222   protected DoubleVector R114;
223   protected DoubleVector R330;
224   protected DoubleVector R303;
225   protected DoubleVector R033;
226   protected DoubleVector R321;
227   protected DoubleVector R231;
228   protected DoubleVector R213;
229   protected DoubleVector R312;
230   protected DoubleVector R132;
231   protected DoubleVector R123;
232 
233   // Components of the potential, field and field gradient.
234   protected DoubleVector E000; // Potential
235   // l + m + n = 1 (3)   4
236   protected DoubleVector E100; // d/dX
237   protected DoubleVector E010; // d/dY
238   protected DoubleVector E001; // d/dz
239   // l + m + n = 2 (6)  10
240   protected DoubleVector E200; // d^2/dXdX
241   protected DoubleVector E020; // d^2/dYdY
242   protected DoubleVector E002; // d^2/dZdZ
243   protected DoubleVector E110; // d^2/dXdY
244   protected DoubleVector E101; // d^2/dXdZ
245   protected DoubleVector E011; // d^2/dYdZ
246   // l + m + n = 3 (10) 20
247   protected DoubleVector E300; // d^3/dXdXdX
248   protected DoubleVector E030; // d^3/dYdYdY
249   protected DoubleVector E003; // d^3/dZdZdZ
250   protected DoubleVector E210; // d^3/dXdXdY
251   protected DoubleVector E201; // d^3/dXdXdZ
252   protected DoubleVector E120; // d^3/dXdYdY
253   protected DoubleVector E021; // d^3/dYdYdZ
254   protected DoubleVector E102; // d^3/dXdZdZ
255   protected DoubleVector E012; // d^3/dYdZdZ
256   protected DoubleVector E111; // d^3/dXdYdZ
257 
258   /**
259    * Constructor for MultipoleTensor.
260    *
261    * @param order       The order of the tensor.
262    * @param coordinates a {@link CoordinateSystem} object.
263    */
264   public MultipoleTensorSIMD(CoordinateSystem coordinates, int order) {
265     assert (order > 0);
266     o1 = order + 1;
267     il = o1;
268     im = il * o1;
269     in = im * o1;
270     size = (order + 1) * (order + 2) * (order + 3) / 6;
271 
272     this.order = order;
273     this.coordinates = coordinates;
274     this.operator = Operator.COULOMB;
275 
276     // Auxiliary terms for Coulomb and Thole Screening.
277     coulombSource = new double[o1];
278     for (short n = 0; n <= order; n++) {
279       /*
280        Math.pow(-1.0, j) returns positive for all j, with -1.0 as the //
281        argument rather than -1. This is a bug?
282        Challacombe Eq. 21, first two factors.
283       */
284       coulombSource[n] = pow(-1, n) * doubleFactorial(2 * n - 1);
285     }
286 
287     work = new DoubleVector[in * o1];
288 
289     // l + m + n = 0 (1)
290     t000 = MultipoleUtilities.ti(0, 0, 0, order);
291     // l + m + n = 1 (3)   4
292     t100 = MultipoleUtilities.ti(1, 0, 0, order);
293     t010 = MultipoleUtilities.ti(0, 1, 0, order);
294     t001 = MultipoleUtilities.ti(0, 0, 1, order);
295     // l + m + n = 2 (6)  10
296     t200 = MultipoleUtilities.ti(2, 0, 0, order);
297     t020 = MultipoleUtilities.ti(0, 2, 0, order);
298     t002 = MultipoleUtilities.ti(0, 0, 2, order);
299     t110 = MultipoleUtilities.ti(1, 1, 0, order);
300     t101 = MultipoleUtilities.ti(1, 0, 1, order);
301     t011 = MultipoleUtilities.ti(0, 1, 1, order);
302     // l + m + n = 3 (10) 20
303     t300 = MultipoleUtilities.ti(3, 0, 0, order);
304     t030 = MultipoleUtilities.ti(0, 3, 0, order);
305     t003 = MultipoleUtilities.ti(0, 0, 3, order);
306     t210 = MultipoleUtilities.ti(2, 1, 0, order);
307     t201 = MultipoleUtilities.ti(2, 0, 1, order);
308     t120 = MultipoleUtilities.ti(1, 2, 0, order);
309     t021 = MultipoleUtilities.ti(0, 2, 1, order);
310     t102 = MultipoleUtilities.ti(1, 0, 2, order);
311     t012 = MultipoleUtilities.ti(0, 1, 2, order);
312     t111 = MultipoleUtilities.ti(1, 1, 1, order);
313     // l + m + n = 4 (15) 35
314     t400 = MultipoleUtilities.ti(4, 0, 0, order);
315     t040 = MultipoleUtilities.ti(0, 4, 0, order);
316     t004 = MultipoleUtilities.ti(0, 0, 4, order);
317     t310 = MultipoleUtilities.ti(3, 1, 0, order);
318     t301 = MultipoleUtilities.ti(3, 0, 1, order);
319     t130 = MultipoleUtilities.ti(1, 3, 0, order);
320     t031 = MultipoleUtilities.ti(0, 3, 1, order);
321     t103 = MultipoleUtilities.ti(1, 0, 3, order);
322     t013 = MultipoleUtilities.ti(0, 1, 3, order);
323     t220 = MultipoleUtilities.ti(2, 2, 0, order);
324     t202 = MultipoleUtilities.ti(2, 0, 2, order);
325     t022 = MultipoleUtilities.ti(0, 2, 2, order);
326     t211 = MultipoleUtilities.ti(2, 1, 1, order);
327     t121 = MultipoleUtilities.ti(1, 2, 1, order);
328     t112 = MultipoleUtilities.ti(1, 1, 2, order);
329     // l + m + n = 5 (21) 56
330     t500 = MultipoleUtilities.ti(5, 0, 0, order);
331     t050 = MultipoleUtilities.ti(0, 5, 0, order);
332     t005 = MultipoleUtilities.ti(0, 0, 5, order);
333     t410 = MultipoleUtilities.ti(4, 1, 0, order);
334     t401 = MultipoleUtilities.ti(4, 0, 1, order);
335     t140 = MultipoleUtilities.ti(1, 4, 0, order);
336     t041 = MultipoleUtilities.ti(0, 4, 1, order);
337     t104 = MultipoleUtilities.ti(1, 0, 4, order);
338     t014 = MultipoleUtilities.ti(0, 1, 4, order);
339     t320 = MultipoleUtilities.ti(3, 2, 0, order);
340     t302 = MultipoleUtilities.ti(3, 0, 2, order);
341     t230 = MultipoleUtilities.ti(2, 3, 0, order);
342     t032 = MultipoleUtilities.ti(0, 3, 2, order);
343     t203 = MultipoleUtilities.ti(2, 0, 3, order);
344     t023 = MultipoleUtilities.ti(0, 2, 3, order);
345     t311 = MultipoleUtilities.ti(3, 1, 1, order);
346     t131 = MultipoleUtilities.ti(1, 3, 1, order);
347     t113 = MultipoleUtilities.ti(1, 1, 3, order);
348     t221 = MultipoleUtilities.ti(2, 2, 1, order);
349     t212 = MultipoleUtilities.ti(2, 1, 2, order);
350     t122 = MultipoleUtilities.ti(1, 2, 2, order);
351     // l + m + n = 6 (28) 84
352     t600 = MultipoleUtilities.ti(6, 0, 0, order);
353     t060 = MultipoleUtilities.ti(0, 6, 0, order);
354     t006 = MultipoleUtilities.ti(0, 0, 6, order);
355     t510 = MultipoleUtilities.ti(5, 1, 0, order);
356     t501 = MultipoleUtilities.ti(5, 0, 1, order);
357     t150 = MultipoleUtilities.ti(1, 5, 0, order);
358     t051 = MultipoleUtilities.ti(0, 5, 1, order);
359     t105 = MultipoleUtilities.ti(1, 0, 5, order);
360     t015 = MultipoleUtilities.ti(0, 1, 5, order);
361     t420 = MultipoleUtilities.ti(4, 2, 0, order);
362     t402 = MultipoleUtilities.ti(4, 0, 2, order);
363     t240 = MultipoleUtilities.ti(2, 4, 0, order);
364     t042 = MultipoleUtilities.ti(0, 4, 2, order);
365     t204 = MultipoleUtilities.ti(2, 0, 4, order);
366     t024 = MultipoleUtilities.ti(0, 2, 4, order);
367     t411 = MultipoleUtilities.ti(4, 1, 1, order);
368     t141 = MultipoleUtilities.ti(1, 4, 1, order);
369     t114 = MultipoleUtilities.ti(1, 1, 4, order);
370     t330 = MultipoleUtilities.ti(3, 3, 0, order);
371     t303 = MultipoleUtilities.ti(3, 0, 3, order);
372     t033 = MultipoleUtilities.ti(0, 3, 3, order);
373     t321 = MultipoleUtilities.ti(3, 2, 1, order);
374     t231 = MultipoleUtilities.ti(2, 3, 1, order);
375     t213 = MultipoleUtilities.ti(2, 1, 3, order);
376     t312 = MultipoleUtilities.ti(3, 1, 2, order);
377     t132 = MultipoleUtilities.ti(1, 3, 2, order);
378     t123 = MultipoleUtilities.ti(1, 2, 3, order);
379     t222 = MultipoleUtilities.ti(2, 2, 2, order);
380   }
381 
382   /**
383    * Set the separation vector.
384    *
385    * @param r The separation vector.
386    */
387   public final void setR(DoubleVector[] r) {
388     setR(r[0], r[1], r[2]);
389   }
390 
391   /**
392    * Set the separation vector.
393    *
394    * @param dx Separation along the X-axis.
395    * @param dy Separation along the Y-axis.
396    * @param dz Separation along the Z-axis.
397    */
398   public abstract void setR(DoubleVector dx, DoubleVector dy, DoubleVector dz);
399 
400   /**
401    * Generate the tensor using hard-coded methods.
402    */
403   public void generateTensor() {
404     switch (order) {
405       case 1 -> order1();
406       case 2 -> order2();
407       case 3 -> order3();
408       case 4 -> order4();
409       case 5 -> order5();
410       case 6 -> order6();
411       default -> {
412         if (logger.isLoggable(Level.WARNING)) {
413           logger.severe("Order " + order + " not supported.");
414         }
415       }
416     }
417   }
418 
419   /**
420    * Return the source terms.
421    *
422    * @return A DoubleVector array of source terms.
423    */
424   public DoubleVector[] getSource() {
425     source(work);
426     return work;
427   }
428 
429   /**
430    * Contract a multipole with the potential and its derivatives.
431    *
432    * @param m PolarizableMultipole at the site of the potential.
433    * @return The permanent multipole energy.
434    */
435   protected final DoubleVector multipoleEnergy(PolarizableMultipoleSIMD m) {
436     DoubleVector total = m.q.mul(E000);
437     total = m.dx.fma(E100, total);
438     total = m.dy.fma(E010, total);
439     total = m.dz.fma(E001, total);
440     total = m.qxx.fma(E200, total);
441     total = m.qyy.fma(E020, total);
442     total = m.qzz.fma(E002, total);
443     total = m.qxy.fma(E110, total);
444     total = m.qxz.fma(E101, total);
445     total = m.qyz.fma(E011, total);
446     return total;
447   }
448 
449   /**
450    * Contract a multipole with the potential and its derivatives.
451    *
452    * @param mI PolarizableMultipole at site I.
453    * @param mK PolarizableMultipole at site K.
454    * @return The permanent multipole energy.
455    */
456   public DoubleVector multipoleEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
457     multipoleIPotentialAtK(mI, 2);
458     return multipoleEnergy(mK);
459   }
460 
461   /**
462    * Compute the permanent multipole gradient.
463    *
464    * @param m PolarizableMultipole at the site of the potential.
465    * @param g The atomic gradient.
466    */
467   protected final void multipoleGradient(PolarizableMultipoleSIMD m, DoubleVector[] g) {
468     // dEnergy/dX
469     DoubleVector total = m.q.mul(E100);
470     total = m.dx.fma(E200, total);
471     total = m.dy.fma(E110, total);
472     total = m.dz.fma(E101, total);
473     total = m.qxx.fma(E300, total);
474     total = m.qyy.fma(E120, total);
475     total = m.qzz.fma(E102, total);
476     total = m.qxy.fma(E210, total);
477     total = m.qxz.fma(E201, total);
478     total = m.qyz.fma(E111, total);
479     g[0] = total;
480 
481     // dEnergy/dY
482     total = m.q.mul(E010);
483     total = m.dx.fma(E110, total);
484     total = m.dy.fma(E020, total);
485     total = m.dz.fma(E011, total);
486     total = m.qxx.fma(E210, total);
487     total = m.qyy.fma(E030, total);
488     total = m.qzz.fma(E012, total);
489     total = m.qxy.fma(E120, total);
490     total = m.qxz.fma(E111, total);
491     total = m.qyz.fma(E021, total);
492     g[1] = total;
493 
494     // dEnergy/dZ
495     total = m.q.mul(E001);
496     total = m.dx.fma(E101, total);
497     total = m.dy.fma(E011, total);
498     total = m.dz.fma(E002, total);
499     total = m.qxx.fma(E201, total);
500     total = m.qyy.fma(E021, total);
501     total = m.qzz.fma(E003, total);
502     total = m.qxy.fma(E111, total);
503     total = m.qxz.fma(E102, total);
504     total = m.qyz.fma(E012, total);
505     g[2] = total;
506   }
507 
508   /**
509    * Compute the torque on a permanent multipole.
510    *
511    * @param m      PolarizableMultipole at the site of the potential.
512    * @param torque an array of double values.
513    */
514   protected final void multipoleTorque(PolarizableMultipoleSIMD m, DoubleVector[] torque) {
515     // Torque on the permanent dipole due to the field.
516     DoubleVector dx = m.dy.mul(E001).sub(m.dz.mul(E010));
517     DoubleVector dy = m.dz.mul(E100).sub(m.dx.mul(E001));
518     DoubleVector dz = m.dx.mul(E010).sub(m.dy.mul(E100));
519 
520     // Torque on the permanent quadrupole due to the gradient of the field.
521     DoubleVector qx = m.qxy.mul(E101).add(m.qyy.mul(E011).mul(2.0)).add(m.qyz.mul(E002))
522         .sub(m.qxz.mul(E110).add(m.qyz.mul(E020)).add(m.qzz.mul(E011).mul(2.0)));
523     DoubleVector qy = m.qxz.mul(E200).add(m.qyz.mul(E110)).add(m.qzz.mul(E101).mul(2.0))
524         .sub(m.qxx.mul(E101).mul(2.0).add(m.qxy.mul(E011)).add(m.qxz.mul(E002)));
525     DoubleVector qz = m.qxx.mul(E110).mul(2.0).add(m.qxy.mul(E020)).add(m.qxz.mul(E011))
526         .sub(m.qxy.mul(E200).add(m.qyy.mul(E110).mul(2.0)).add(m.qyz.mul(E101)));
527 
528     // The field along X is -E001, so we need a negative sign.
529     torque[0] = torque[0].sub(dx.add(qx));
530     torque[1] = torque[1].sub(dy.add(qy));
531     torque[2] = torque[2].sub(dz.add(qz));
532   }
533 
534   /**
535    * Compute the torque on a permanent dipole.
536    *
537    * @param m      PolarizableMultipole at the site of the potential.
538    * @param torque an array of double values.
539    */
540   protected final void dipoleTorque(PolarizableMultipoleSIMD m, DoubleVector[] torque) {
541     // Torque on the permanent dipole due to the field.
542     DoubleVector dx = m.dy.mul(E001).sub(m.dz.mul(E010));
543     DoubleVector dy = m.dz.mul(E100).sub(m.dx.mul(E001));
544     DoubleVector dz = m.dx.mul(E010).sub(m.dy.mul(E100));
545 
546     // The field along X is -E001, so we need a negative sign.
547     torque[0] = torque[0].sub(dx);
548     torque[1] = torque[1].sub(dy);
549     torque[2] = torque[2].sub(dz);
550   }
551 
552   /**
553    * Compute the torque on a permanent quadrupole.
554    *
555    * @param m      PolarizableMultipole at the site of the potential.
556    * @param torque an array of double values.
557    */
558   protected final void quadrupoleTorque(PolarizableMultipoleSIMD m, DoubleVector[] torque) {
559     // Torque on the permanent quadrupole due to the gradient of the field.
560     DoubleVector qx = m.qxy.mul(E101).add(m.qyy.mul(E011).mul(2.0)).add(m.qyz.mul(E002))
561         .sub(m.qxz.mul(E110).add(m.qyz.mul(E020)).add(m.qzz.mul(E011).mul(2.0)));
562     DoubleVector qy = m.qxz.mul(E200).add(m.qyz.mul(E110)).add(m.qzz.mul(E101).mul(2.0))
563         .sub(m.qxx.mul(E101).mul(2.0).add(m.qxy.mul(E011)).add(m.qxz.mul(E002)));
564     DoubleVector qz = m.qxx.mul(E110).mul(2.0).add(m.qxy.mul(E020)).add(m.qxz.mul(E011))
565         .sub(m.qxy.mul(E200).add(m.qyy.mul(E110).mul(2.0)).add(m.qyz.mul(E101)));
566 
567     // The field along X is -E001, so we need a negative sign.
568     torque[0] = torque[0].sub(qx);
569     torque[1] = torque[1].sub(qy);
570     torque[2] = torque[2].sub(qz);
571   }
572 
573   /**
574    * Contract an induced dipole with the potential and its derivatives.
575    *
576    * @param m PolarizableMultipole at the site of the potential.
577    * @return The polarization energy.
578    */
579   protected final DoubleVector polarizationEnergy(PolarizableMultipoleSIMD m) {
580     // E = -1/2 * u.E
581     // No negative sign because the field E = [-E100, -E010, -E001].
582     return (m.ux.mul(E100).add(m.uy.mul(E010)).add(m.uz.mul(E001))).mul(.5);
583   }
584 
585   /**
586    * Contract an induced dipole with the potential and its derivatives.
587    *
588    * @param m PolarizableMultipole at the site of the potential.
589    * @return The polarization energy.
590    */
591   protected final DoubleVector polarizationEnergyS(PolarizableMultipoleSIMD m) {
592     // E = -1/2 * u.E
593     // No negative sign because the field E = [-E100, -E010, -E001].
594     return (m.sx.mul(E100).add(m.sy.mul(E010)).add(m.sz.mul(E001))).mul(.5);
595   }
596 
597   /**
598    * Polarization Energy and Gradient.
599    *
600    * @param mI            PolarizableMultipole at site I.
601    * @param mK            PolarizableMultipole at site K.
602    * @param inductionMask a double.
603    * @param energyMask    a double.
604    * @param mutualMask    a double.
605    * @param Gi            an array of double values.
606    * @param Ti            an array of double values.
607    * @param Tk            an array of double values.
608    * @return a double.
609    */
610   public DoubleVector polarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
611                                                     DoubleVector inductionMask, DoubleVector energyMask, DoubleVector mutualMask,
612                                                     DoubleVector[] Gi, DoubleVector[] Ti, DoubleVector[] Tk) {
613 
614     // Add the induction and energy masks to create an "averaged" induced dipole (sx, sy, sz).
615     mI.applyMasks(inductionMask, energyMask);
616     mK.applyMasks(inductionMask, energyMask);
617 
618     // Find the permanent multipole potential and derivatives at site k.
619     multipoleIPotentialAtK(mI, 2);
620     // Energy of induced dipole k in the field of multipole i.
621     // The field E_x = -E100.
622     DoubleVector eK = polarizationEnergy(mK);
623     // Derivative with respect to moving atom k.
624     Gi[0] = mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101)).neg();
625     Gi[1] = mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011)).neg();
626     Gi[2] = mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002)).neg();
627 
628     // Find the permanent multipole potential and derivatives at site i.
629     multipoleKPotentialAtI(mK, 2);
630     // Energy of induced dipole i in the field of multipole k.
631     DoubleVector eI = polarizationEnergy(mI);
632     // Derivative with respect to moving atom i.
633     Gi[0] = Gi[0].add(mI.sx.mul(E200).add(mI.sy.mul(E110)).add(mI.sz.mul(E101)));
634     Gi[1] = Gi[1].add(mI.sx.mul(E110).add(mI.sy.mul(E020)).add(mI.sz.mul(E011)));
635     Gi[2] = Gi[2].add(mI.sx.mul(E101).add(mI.sy.mul(E011)).add(mI.sz.mul(E002)));
636 
637     // Total polarization energy.
638     DoubleVector energy = eI.add(eK).mul(energyMask);
639 
640     // Get the induced-induced portion of the force (Ud . dC/dX . Up).
641     // This contribution does not exist for direct polarization (mutualMask == 0.0).
642     // For SIMD code, we are unable to hide this in an if statement, calculation is still correct with it
643     // Find the potential and its derivatives at k due to induced dipole i.
644     dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
645     Gi[0] = Gi[0].sub(mK.px.mul(E200).add(mK.py.mul(E110)).add(mK.pz.mul(E101)).mul(0.5).mul(mutualMask));
646     Gi[1] = Gi[1].sub(mK.px.mul(E110).add(mK.py.mul(E020)).add(mK.pz.mul(E011)).mul(0.5).mul(mutualMask));
647     Gi[2] = Gi[2].sub(mK.px.mul(E101).add(mK.py.mul(E011)).add(mK.pz.mul(E002)).mul(0.5).mul(mutualMask));
648 
649     // Find the potential and its derivatives at i due to induced dipole k.
650     dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
651     Gi[0] = Gi[0].add(mI.px.mul(E200).add(mI.py.mul(E110)).add(mI.pz.mul(E101)).mul(0.5).mul(mutualMask));
652     Gi[1] = Gi[1].add(mI.px.mul(E110).add(mI.py.mul(E020)).add(mI.pz.mul(E011)).mul(0.5).mul(mutualMask));
653     Gi[2] = Gi[2].add(mI.px.mul(E101).add(mI.py.mul(E011)).add(mI.pz.mul(E002)).mul(0.5).mul(mutualMask));
654 
655     // Find the potential and its derivatives at K due to the averaged induced dipole at site i.
656     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
657     multipoleTorque(mK, Tk);
658 
659     // Find the potential and its derivatives at I due to the averaged induced dipole at site k.
660     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
661     multipoleTorque(mI, Ti);
662 
663     return energy;
664   }
665 
666   /**
667    * Permanent multipole energy and gradient.
668    *
669    * @param mI PolarizableMultipole at site I.
670    * @param mK PolarizableMultipole at site K.
671    * @param Gi Coordinate gradient at site I.
672    * @param Gk Coordinate gradient at site K.
673    * @param Ti Torque at site I.
674    * @param Tk Torque at site K.
675    * @return the permanent multipole energy.
676    */
677   public DoubleVector multipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
678                                                  DoubleVector[] Gi, DoubleVector[] Gk, DoubleVector[] Ti, DoubleVector[] Tk) {
679     multipoleIPotentialAtK(mI, 3);
680     DoubleVector energy = multipoleEnergy(mK);
681     multipoleGradient(mK, Gk);
682     Gi[0] = Gi[0].sub(Gk[0]);
683     Gi[1] = Gi[1].sub(Gk[1]);
684     Gi[2] = Gi[2].sub(Gk[2]);
685 
686     // Torques
687     multipoleTorque(mK, Tk);
688     multipoleKPotentialAtI(mK, 2);
689     multipoleTorque(mI, Ti);
690 
691     // dEdZ = -Gi[2];
692     // if (order >= 6) {
693     //  multipoleIdZ2(mI);
694     //  d2EdZ2 = dotMultipole(mK);
695     // }
696 
697     return energy;
698   }
699 
700   /**
701    * Generate source terms for the Challacombe et al. recursion.
702    *
703    * @param T000 Location to store the source terms.
704    */
705   protected abstract void source(DoubleVector[] T000);
706 
707   /**
708    * Hard coded tensor computation up to 1st order. This code is auto-generated for both Global and QI frames.
709    */
710   protected abstract void order1();
711 
712   /**
713    * Hard coded tensor computation up to 2nd order. This code is auto-generated for both Global and QI frames.
714    */
715   protected abstract void order2();
716 
717   /**
718    * Hard coded tensor computation up to 3rd order. This code is auto-generated for both Global and QI frames.
719    */
720   protected abstract void order3();
721 
722   /**
723    * Hard coded tensor computation up to 4th order. This code is auto-generated for both Global and QI frames.
724    */
725   protected abstract void order4();
726 
727   /**
728    * Hard coded tensor computation up to 5th order. This code is auto-generated for both Global and QI frames.
729    * <br>
730    * The 5th order recursion is needed for quadrupole-quadrupole forces.
731    */
732   protected abstract void order5();
733 
734   /**
735    * Hard coded tensor computation up to 6th order. This code is auto-generated for both Global and QI frames.
736    * <br>
737    * This is needed for quadrupole-quadrupole forces and orthogonal space sampling.
738    */
739   protected abstract void order6();
740 
741   /**
742    * Compute the field components due to site I multipole at site K.
743    *
744    * @param mI    PolarizableMultipoleSIMD at site I.
745    * @param order Potential order.
746    */
747   @SuppressWarnings("fallthrough")
748   protected abstract void multipoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order);
749 
750   /**
751    * Compute the field components due to site K multipole at site I.
752    *
753    * @param mK    PolarizableMultipoleSIMD at site I.
754    * @param order Potential order.
755    */
756   @SuppressWarnings("fallthrough")
757   protected abstract void multipoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order);
758 
759   /**
760    * Compute the field components due to site K charge at site I.
761    *
762    * @param mK    PolarizableMultipoleSIMD at site K.
763    * @param order Potential order.
764    */
765   protected abstract void chargeKPotentialAtI(PolarizableMultipoleSIMD mK, int order);
766 
767   /**
768    * Compute the induced dipole field components due to site K at site I.
769    *
770    * @param uxk   X-dipole component.
771    * @param uyk   Y-dipole component.
772    * @param uzk   Z-dipole component.
773    * @param order Potential order.
774    */
775   protected abstract void dipoleKPotentialAtI(DoubleVector uxk, DoubleVector uyk, DoubleVector uzk, int order);
776 
777   /**
778    * Compute the field components due to site K quadrupole at site I.
779    *
780    * @param mK    MultipoleTensorSIMD at site K.
781    * @param order Potential order.
782    */
783   protected abstract void quadrupoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order);
784 
785   /**
786    * Compute the field components due to site I charge at site K.
787    *
788    * @param mI    PolarizableMultipoleSIMD at site I.
789    * @param order Potential order.
790    */
791   protected abstract void chargeIPotentialAtK(PolarizableMultipoleSIMD mI, int order);
792 
793   /**
794    * Compute the induced dipole field components due to site I at site K.
795    *
796    * @param uxi   X-dipole component.
797    * @param uyi   Y-dipole component.
798    * @param uzi   Z-dipole component.
799    * @param order Potential order.
800    */
801   protected abstract void dipoleIPotentialAtK(DoubleVector uxi, DoubleVector uyi, DoubleVector uzi, int order);
802 
803   /**
804    * Compute the field components due to site I quadrupole at site K.
805    *
806    * @param mI    MultipoleTensorSIMD at site I.
807    * @param order Potential order.
808    */
809   protected abstract void quadrupoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order);
810 
811   /**
812    * The index is based on the idea of filling tetrahedron.
813    * <p>
814    * 1/r has an index of 0.
815    * <br>
816    * derivatives of x are first; indices from 1..o for d/dx..(d/dx)^o
817    * <br>
818    * derivatives of x and y are second; base triangle of size (o+1)(o+2)/2
819    * <br>
820    * derivatives of x, y and z are last; total size (o+1)*(o+2)*(o+3)/6
821    * <br>
822    * <p>
823    * This function is useful to set up masking constants:
824    * <br>
825    * static int Tlmn = ti(l,m,n,order)
826    * <br>
827    * For example the (d/dy)^2 (1/R) storage location:
828    * <br>
829    * static int T020 = ti(0,2,0,order)
830    *
831    * @param dx int The number of d/dx operations.
832    * @param dy int The number of d/dy operations.
833    * @param dz int The number of d/dz operations.
834    * @return int in the range (0..binomial(order + 3, 3) - 1)
835    */
836   protected final int ti(int dx, int dy, int dz) {
837     return MultipoleUtilities.ti(dx, dy, dz, order);
838   }
839 
840   // l + m + n = 0 (1)
841   /**
842    * No derivatives.
843    */
844   protected final int t000;
845   // l + m + n = 1 (3)   4
846   /**
847    * First derivative with respect to x.
848    */
849   protected final int t100;
850   /**
851    * First derivative with respect to y.
852    */
853   protected final int t010;
854   /**
855    * First derivative with respect to z.
856    */
857   protected final int t001;
858   // l + m + n = 2 (6)  10
859   /**
860    * Second derivative with respect to x.
861    */
862   protected final int t200;
863   /**
864    * Second derivative with respect to y.
865    */
866   protected final int t020;
867   /**
868    * Second derivative with respect to z.
869    */
870   protected final int t002;
871   /**
872    * Derivatives with respect to x and y.
873    */
874   protected final int t110;
875   /**
876    * Derivatives with respect to x and z.
877    */
878   protected final int t101;
879   /**
880    * Derivatives with respect to y and z.
881    */
882   protected final int t011;
883   // l + m + n = 3 (10) 20
884   /**
885    * Third derivative with respect to x.
886    */
887   protected final int t300;
888   /**
889    * Third derivative with respect to y.
890    */
891   protected final int t030;
892   /**
893    * Third derivative with respect to z.
894    */
895   protected final int t003;
896   /**
897    * Derivatives with respect to x and y.
898    */
899   protected final int t210;
900   /**
901    * Derivatives with respect to x and z.
902    */
903   protected final int t201;
904   /**
905    * Derivatives with respect to x and y.
906    */
907   protected final int t120;
908   /**
909    * Derivatives with respect to y and z.
910    */
911   protected final int t021;
912   /**
913    * Derivatives with respect to x and z.
914    */
915   protected final int t102;
916   /**
917    * Derivatives with respect to y and z.
918    */
919   protected final int t012;
920   /**
921    * Derivatives with respect to x, y and z.
922    */
923   protected final int t111;
924   // l + m + n = 4 (15) 35
925   /**
926    * Fourth derivative with respect to x.
927    */
928   protected final int t400;
929   /**
930    * Fourth derivative with respect to y.
931    */
932   protected final int t040;
933   /**
934    * Fourth derivative with respect to z.
935    */
936   protected final int t004;
937   /**
938    * Derivatives with respect to x and y.
939    */
940   protected final int t310;
941   /**
942    * Derivatives with respect to x and z.
943    */
944   protected final int t301;
945   /**
946    * Derivatives with respect to x and y.
947    */
948   protected final int t130;
949   /**
950    * Derivatives with respect to y and z.
951    */
952   protected final int t031;
953   /**
954    * Derivatives with respect to x and z.
955    */
956   protected final int t103;
957   /**
958    * Derivatives with respect to y and z.
959    */
960   protected final int t013;
961   /**
962    * Derivatives with respect to x and y.
963    */
964   protected final int t220;
965   /**
966    * Derivatives with respect to x and z.
967    */
968   protected final int t202;
969   /**
970    * Derivatives with respect to y and z.
971    */
972   protected final int t022;
973   /**
974    * Derivatives with respect to x, y and z.
975    */
976   protected final int t211;
977   /**
978    * Derivatives with respect to x, y and z.
979    */
980   protected final int t121;
981   /**
982    * Derivatives with respect to x, y and z.
983    */
984   protected final int t112;
985 
986   // l + m + n = 5 (21) 56
987   /**
988    * Fifth derivative with respect to x.
989    */
990   protected final int t500;
991   /**
992    * Fifth derivative with respect to y.
993    */
994   protected final int t050;
995   /**
996    * Fifth derivative with respect to z.
997    */
998   protected final int t005;
999   /**
1000    * Derivatives with respect to x and y.
1001    */
1002   protected final int t410;
1003   /**
1004    * Derivatives with respect to x and z.
1005    */
1006   protected final int t401;
1007   /**
1008    * Derivatives with respect to x and y.
1009    */
1010   protected final int t140;
1011   /**
1012    * Derivatives with respect to y and z.
1013    */
1014   protected final int t041;
1015   /**
1016    * Derivatives with respect to x and z.
1017    */
1018   protected final int t104;
1019   /**
1020    * Derivatives with respect to y and z.
1021    */
1022   protected final int t014;
1023   /**
1024    * Derivatives with respect to x and y.
1025    */
1026   protected final int t320;
1027   /**
1028    * Derivatives with respect to x and z.
1029    */
1030   protected final int t302;
1031   /**
1032    * Derivatives with respect to x and y.
1033    */
1034   protected final int t230;
1035   /**
1036    * Derivatives with respect to y and z.
1037    */
1038   protected final int t032;
1039   /**
1040    * Derivatives with respect to x and z.
1041    */
1042   protected final int t203;
1043   /**
1044    * Derivatives with respect to y and z.
1045    */
1046   protected final int t023;
1047   /**
1048    * Derivatives with respect to x, y and z.
1049    */
1050   protected final int t311;
1051   /**
1052    * Derivatives with respect to x, y and z.
1053    */
1054   protected final int t131;
1055   /**
1056    * Derivatives with respect to x, y and z.
1057    */
1058   protected final int t113;
1059   /**
1060    * Derivatives with respect to x, y and z.
1061    */
1062   protected final int t221;
1063   /**
1064    * Derivatives with respect to x, y and z.
1065    */
1066   protected final int t212;
1067   /**
1068    * Derivatives with respect to x, y and z.
1069    */
1070   protected final int t122;
1071 
1072   // l + m + n = 6 (28) 84
1073   /**
1074    * Sixth derivative with respect to x.
1075    */
1076   protected final int t600;
1077   /**
1078    * Sixth derivative with respect to y.
1079    */
1080   protected final int t060;
1081   /**
1082    * Sixth derivative with respect to z.
1083    */
1084   protected final int t006;
1085   /**
1086    * Derivatives with respect to x and y.
1087    */
1088   protected final int t510;
1089   /**
1090    * Derivatives with respect to x and z.
1091    */
1092   protected final int t501;
1093   /**
1094    * Derivatives with respect to x and y.
1095    */
1096   protected final int t150;
1097   /**
1098    * Derivatives with respect to y and z.
1099    */
1100   protected final int t051;
1101   /**
1102    * Derivatives with respect to x and z.
1103    */
1104   protected final int t105;
1105   /**
1106    * Derivatives with respect to y and z.
1107    */
1108   protected final int t015;
1109   /**
1110    * Derivatives with respect to x and y.
1111    */
1112   protected final int t420;
1113   /**
1114    * Derivatives with respect to x and z.
1115    */
1116   protected final int t402;
1117   /**
1118    * Derivatives with respect to x and y.
1119    */
1120   protected final int t240;
1121   /**
1122    * Derivatives with respect to y and z.
1123    */
1124   protected final int t042;
1125   /**
1126    * Derivatives with respect to x and z.
1127    */
1128   protected final int t204;
1129   /**
1130    * Derivatives with respect to y and z.
1131    */
1132   protected final int t024;
1133   /**
1134    * Derivatives with respect to x, y and z.
1135    */
1136   protected final int t411;
1137   /**
1138    * Derivatives with respect to x, y and z.
1139    */
1140   protected final int t141;
1141   /**
1142    * Derivatives with respect to x, y and z.
1143    */
1144   protected final int t114;
1145   /**
1146    * Derivatives with respect to x and y.
1147    */
1148   protected final int t330;
1149   /**
1150    * Derivatives with respect to x and z.
1151    */
1152   protected final int t303;
1153   /**
1154    * Derivatives with respect to y and z.
1155    */
1156   protected final int t033;
1157   /**
1158    * Derivatives with respect to x, y and z.
1159    */
1160   protected final int t321;
1161   /**
1162    * Derivatives with respect to x, y and z.
1163    */
1164   protected final int t231;
1165   /**
1166    * Derivatives with respect to x, y and z.
1167    */
1168   protected final int t213;
1169   /**
1170    * Derivatives with respect to x, y and z.
1171    */
1172   protected final int t312;
1173   /**
1174    * Derivatives with respect to x, y and z.
1175    */
1176   protected final int t132;
1177   /**
1178    * Derivatives with respect to x, y and z.
1179    */
1180   protected final int t123;
1181   /**
1182    * Derivatives with respect to x, y and z.
1183    */
1184   protected final int t222;
1185 
1186 }