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   public DoubleVector[] getSource() {
423     source(work);
424     return work;
425   }
426 
427   /**
428    * Contract a multipole with the potential and its derivatives.
429    *
430    * @param m PolarizableMultipole at the site of the potential.
431    * @return The permanent multipole energy.
432    */
433   protected final DoubleVector multipoleEnergy(PolarizableMultipoleSIMD m) {
434     DoubleVector total = m.q.mul(E000);
435     total = m.dx.fma(E100, total);
436     total = m.dy.fma(E010, total);
437     total = m.dz.fma(E001, total);
438     total = m.qxx.fma(E200, total);
439     total = m.qyy.fma(E020, total);
440     total = m.qzz.fma(E002, total);
441     total = m.qxy.fma(E110, total);
442     total = m.qxz.fma(E101, total);
443     total = m.qyz.fma(E011, total);
444     return total;
445   }
446 
447   public DoubleVector multipoleEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
448     multipoleIPotentialAtK(mI, 2);
449     return multipoleEnergy(mK);
450   }
451 
452   /**
453    * Compute the permanent multipole gradient.
454    *
455    * @param m PolarizableMultipole at the site of the potential.
456    * @param g The atomic gradient.
457    */
458   protected final void multipoleGradient(PolarizableMultipoleSIMD m, DoubleVector[] g) {
459     // dEnergy/dX
460     DoubleVector total = m.q.mul(E100);
461     total = m.dx.fma(E200, total);
462     total = m.dy.fma(E110, total);
463     total = m.dz.fma(E101, total);
464     total = m.qxx.fma(E300, total);
465     total = m.qyy.fma(E120, total);
466     total = m.qzz.fma(E102, total);
467     total = m.qxy.fma(E210, total);
468     total = m.qxz.fma(E201, total);
469     total = m.qyz.fma(E111, total);
470     g[0] = total;
471 
472     // dEnergy/dY
473     total = m.q.mul(E010);
474     total = m.dx.fma(E110, total);
475     total = m.dy.fma(E020, total);
476     total = m.dz.fma(E011, total);
477     total = m.qxx.fma(E210, total);
478     total = m.qyy.fma(E030, total);
479     total = m.qzz.fma(E012, total);
480     total = m.qxy.fma(E120, total);
481     total = m.qxz.fma(E111, total);
482     total = m.qyz.fma(E021, total);
483     g[1] = total;
484 
485     // dEnergy/dZ
486     total = m.q.mul(E001);
487     total = m.dx.fma(E101, total);
488     total = m.dy.fma(E011, total);
489     total = m.dz.fma(E002, total);
490     total = m.qxx.fma(E201, total);
491     total = m.qyy.fma(E021, total);
492     total = m.qzz.fma(E003, total);
493     total = m.qxy.fma(E111, total);
494     total = m.qxz.fma(E102, total);
495     total = m.qyz.fma(E012, total);
496     g[2] = total;
497   }
498 
499   /**
500    * Compute the torque on a permanent multipole.
501    *
502    * @param m      PolarizableMultipole at the site of the potential.
503    * @param torque an array of {@link double} objects.
504    */
505   protected final void multipoleTorque(PolarizableMultipoleSIMD m, DoubleVector[] torque) {
506     // Torque on the permanent dipole due to the field.
507     DoubleVector dx = m.dy.mul(E001).sub(m.dz.mul(E010));
508     DoubleVector dy = m.dz.mul(E100).sub(m.dx.mul(E001));
509     DoubleVector dz = m.dx.mul(E010).sub(m.dy.mul(E100));
510 
511     // Torque on the permanent quadrupole due to the gradient of the field.
512     DoubleVector qx = m.qxy.mul(E101).add(m.qyy.mul(E011).mul(2.0)).add(m.qyz.mul(E002))
513         .sub(m.qxz.mul(E110).add(m.qyz.mul(E020)).add(m.qzz.mul(E011).mul(2.0)));
514     DoubleVector qy = m.qxz.mul(E200).add(m.qyz.mul(E110)).add(m.qzz.mul(E101).mul(2.0))
515         .sub(m.qxx.mul(E101).mul(2.0).add(m.qxy.mul(E011)).add(m.qxz.mul(E002)));
516     DoubleVector qz = m.qxx.mul(E110).mul(2.0).add(m.qxy.mul(E020)).add(m.qxz.mul(E011))
517         .sub(m.qxy.mul(E200).add(m.qyy.mul(E110).mul(2.0)).add(m.qyz.mul(E101)));
518 
519     // The field along X is -E001, so we need a negative sign.
520     torque[0] = torque[0].sub(dx.add(qx));
521     torque[1] = torque[1].sub(dy.add(qy));
522     torque[2] = torque[2].sub(dz.add(qz));
523   }
524 
525   /**
526    * Compute the torque on a permanent dipole.
527    *
528    * @param m      PolarizableMultipole at the site of the potential.
529    * @param torque an array of {@link double} objects.
530    */
531   protected final void dipoleTorque(PolarizableMultipoleSIMD m, DoubleVector[] torque) {
532     // Torque on the permanent dipole due to the field.
533     DoubleVector dx = m.dy.mul(E001).sub(m.dz.mul(E010));
534     DoubleVector dy = m.dz.mul(E100).sub(m.dx.mul(E001));
535     DoubleVector dz = m.dx.mul(E010).sub(m.dy.mul(E100));
536 
537     // The field along X is -E001, so we need a negative sign.
538     torque[0] = torque[0].sub(dx);
539     torque[1] = torque[1].sub(dy);
540     torque[2] = torque[2].sub(dz);
541   }
542 
543   /**
544    * Compute the torque on a permanent quadrupole.
545    *
546    * @param m      PolarizableMultipole at the site of the potential.
547    * @param torque an array of {@link double} objects.
548    */
549   protected final void quadrupoleTorque(PolarizableMultipoleSIMD m, DoubleVector[] torque) {
550     // Torque on the permanent quadrupole due to the gradient of the field.
551     DoubleVector qx = m.qxy.mul(E101).add(m.qyy.mul(E011).mul(2.0)).add(m.qyz.mul(E002))
552         .sub(m.qxz.mul(E110).add(m.qyz.mul(E020)).add(m.qzz.mul(E011).mul(2.0)));
553     DoubleVector qy = m.qxz.mul(E200).add(m.qyz.mul(E110)).add(m.qzz.mul(E101).mul(2.0))
554         .sub(m.qxx.mul(E101).mul(2.0).add(m.qxy.mul(E011)).add(m.qxz.mul(E002)));
555     DoubleVector qz = m.qxx.mul(E110).mul(2.0).add(m.qxy.mul(E020)).add(m.qxz.mul(E011))
556         .sub(m.qxy.mul(E200).add(m.qyy.mul(E110).mul(2.0)).add(m.qyz.mul(E101)));
557 
558     // The field along X is -E001, so we need a negative sign.
559     torque[0] = torque[0].sub(qx);
560     torque[1] = torque[1].sub(qy);
561     torque[2] = torque[2].sub(qz);
562   }
563 
564   /**
565    * Contract an induced dipole with the potential and its derivatives.
566    *
567    * @param m PolarizableMultipole at the site of the potential.
568    * @return The polarization energy.
569    */
570   protected final DoubleVector polarizationEnergy(PolarizableMultipoleSIMD m) {
571     // E = -1/2 * u.E
572     // No negative sign because the field E = [-E100, -E010, -E001].
573     return (m.ux.mul(E100).add(m.uy.mul(E010)).add(m.uz.mul(E001))).mul(.5);
574   }
575 
576   /**
577    * Contract an induced dipole with the potential and its derivatives.
578    *
579    * @param m PolarizableMultipole at the site of the potential.
580    * @return The polarization energy.
581    */
582   protected final DoubleVector polarizationEnergyS(PolarizableMultipoleSIMD m) {
583     // E = -1/2 * u.E
584     // No negative sign because the field E = [-E100, -E010, -E001].
585     return (m.sx.mul(E100).add(m.sy.mul(E010)).add(m.sz.mul(E001))).mul(.5);
586   }
587 
588   /**
589    * Polarization Energy and Gradient.
590    *
591    * @param mI            PolarizableMultipole at site I.
592    * @param mK            PolarizableMultipole at site K.
593    * @param inductionMask a double.
594    * @param energyMask    a double.
595    * @param mutualMask    a double.
596    * @param Gi            an array of {@link double} objects.
597    * @param Ti            an array of {@link double} objects.
598    * @param Tk            an array of {@link double} objects.
599    * @return a double.
600    */
601   public DoubleVector polarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
602                                                     DoubleVector inductionMask, DoubleVector energyMask, DoubleVector mutualMask,
603                                                     DoubleVector[] Gi, DoubleVector[] Ti, DoubleVector[] Tk) {
604 
605     // Add the induction and energy masks to create an "averaged" induced dipole (sx, sy, sz).
606     mI.applyMasks(inductionMask, energyMask);
607     mK.applyMasks(inductionMask, energyMask);
608 
609     // Find the permanent multipole potential and derivatives at site k.
610     multipoleIPotentialAtK(mI, 2);
611     // Energy of induced dipole k in the field of multipole i.
612     // The field E_x = -E100.
613     DoubleVector eK = polarizationEnergy(mK);
614     // Derivative with respect to moving atom k.
615     Gi[0] = mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101)).neg();
616     Gi[1] = mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011)).neg();
617     Gi[2] = mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002)).neg();
618 
619     // Find the permanent multipole potential and derivatives at site i.
620     multipoleKPotentialAtI(mK, 2);
621     // Energy of induced dipole i in the field of multipole k.
622     DoubleVector eI = polarizationEnergy(mI);
623     // Derivative with respect to moving atom i.
624     Gi[0] = Gi[0].add(mI.sx.mul(E200).add(mI.sy.mul(E110)).add(mI.sz.mul(E101)));
625     Gi[1] = Gi[1].add(mI.sx.mul(E110).add(mI.sy.mul(E020)).add(mI.sz.mul(E011)));
626     Gi[2] = Gi[2].add(mI.sx.mul(E101).add(mI.sy.mul(E011)).add(mI.sz.mul(E002)));
627 
628     // Total polarization energy.
629     DoubleVector energy = eI.add(eK).mul(energyMask);
630 
631     // Get the induced-induced portion of the force (Ud . dC/dX . Up).
632     // This contribution does not exist for direct polarization (mutualMask == 0.0).
633     // For SIMD code, we are unable to hide this in an if statement, calculation is still correct with it
634     // Find the potential and its derivatives at k due to induced dipole i.
635     dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
636     Gi[0] = Gi[0].sub(mK.px.mul(E200).add(mK.py.mul(E110)).add(mK.pz.mul(E101)).mul(0.5).mul(mutualMask));
637     Gi[1] = Gi[1].sub(mK.px.mul(E110).add(mK.py.mul(E020)).add(mK.pz.mul(E011)).mul(0.5).mul(mutualMask));
638     Gi[2] = Gi[2].sub(mK.px.mul(E101).add(mK.py.mul(E011)).add(mK.pz.mul(E002)).mul(0.5).mul(mutualMask));
639 
640     // Find the potential and its derivatives at i due to induced dipole k.
641     dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
642     Gi[0] = Gi[0].add(mI.px.mul(E200).add(mI.py.mul(E110)).add(mI.pz.mul(E101)).mul(0.5).mul(mutualMask));
643     Gi[1] = Gi[1].add(mI.px.mul(E110).add(mI.py.mul(E020)).add(mI.pz.mul(E011)).mul(0.5).mul(mutualMask));
644     Gi[2] = Gi[2].add(mI.px.mul(E101).add(mI.py.mul(E011)).add(mI.pz.mul(E002)).mul(0.5).mul(mutualMask));
645 
646     // Find the potential and its derivatives at K due to the averaged induced dipole at site i.
647     dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
648     multipoleTorque(mK, Tk);
649 
650     // Find the potential and its derivatives at I due to the averaged induced dipole at site k.
651     dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
652     multipoleTorque(mI, Ti);
653 
654     return energy;
655   }
656 
657   /**
658    * Permanent multipole energy and gradient.
659    *
660    * @param mI PolarizableMultipole at site I.
661    * @param mK PolarizableMultipole at site K.
662    * @param Gi Coordinate gradient at site I.
663    * @param Gk Coordinate gradient at site K.
664    * @param Ti Torque at site I.
665    * @param Tk Torque at site K.
666    * @return the permanent multipole energy.
667    */
668   public DoubleVector multipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
669                                                  DoubleVector[] Gi, DoubleVector[] Gk, DoubleVector[] Ti, DoubleVector[] Tk) {
670     multipoleIPotentialAtK(mI, 3);
671     DoubleVector energy = multipoleEnergy(mK);
672     multipoleGradient(mK, Gk);
673     Gi[0] = Gi[0].sub(Gk[0]);
674     Gi[1] = Gi[1].sub(Gk[1]);
675     Gi[2] = Gi[2].sub(Gk[2]);
676 
677     // Torques
678     multipoleTorque(mK, Tk);
679     multipoleKPotentialAtI(mK, 2);
680     multipoleTorque(mI, Ti);
681 
682     // dEdZ = -Gi[2];
683     // if (order >= 6) {
684     //  multipoleIdZ2(mI);
685     //  d2EdZ2 = dotMultipole(mK);
686     // }
687 
688     return energy;
689   }
690 
691   /**
692    * Generate source terms for the Challacombe et al. recursion.
693    *
694    * @param T000 Location to store the source terms.
695    */
696   protected abstract void source(DoubleVector[] T000);
697 
698   /**
699    * Hard coded tensor computation up to 1st order. This code is auto-generated for both Global and QI frames.
700    */
701   protected abstract void order1();
702 
703   /**
704    * Hard coded tensor computation up to 2nd order. This code is auto-generated for both Global and QI frames.
705    */
706   protected abstract void order2();
707 
708   /**
709    * Hard coded tensor computation up to 3rd order. This code is auto-generated for both Global and QI frames.
710    */
711   protected abstract void order3();
712 
713   /**
714    * Hard coded tensor computation up to 4th order. This code is auto-generated for both Global and QI frames.
715    */
716   protected abstract void order4();
717 
718   /**
719    * Hard coded tensor computation up to 5th order. This code is auto-generated for both Global and QI frames.
720    * <br>
721    * The 5th order recursion is needed for quadrupole-quadrupole forces.
722    */
723   protected abstract void order5();
724 
725   /**
726    * Hard coded tensor computation up to 6th order. This code is auto-generated for both Global and QI frames.
727    * <br>
728    * This is needed for quadrupole-quadrupole forces and orthogonal space sampling.
729    */
730   protected abstract void order6();
731 
732   @SuppressWarnings("fallthrough")
733   protected abstract void multipoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order);
734 
735   @SuppressWarnings("fallthrough")
736   protected abstract void multipoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order);
737 
738   /**
739    * Compute the field components due to site K charge at site I.
740    *
741    * @param mK    MultipoleTensorSIMD at site K.
742    * @param order Potential order.
743    */
744   protected abstract void chargeKPotentialAtI(PolarizableMultipoleSIMD mK, int order);
745 
746   /**
747    * Compute the induced dipole field components due to site K at site I.
748    *
749    * @param uxk   X-dipole component.
750    * @param uyk   Y-dipole component.
751    * @param uzk   Z-dipole component.
752    * @param order Potential order.
753    */
754   protected abstract void dipoleKPotentialAtI(DoubleVector uxk, DoubleVector uyk, DoubleVector uzk, int order);
755 
756   /**
757    * Compute the field components due to site K quadrupole at site I.
758    *
759    * @param mK    MultipoleTensorSIMD at site K.
760    * @param order Potential order.
761    */
762   protected abstract void quadrupoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order);
763 
764   /**
765    * Compute the field components due to site I charge at site K.
766    *
767    * @param mI    MultipoleTensorSIMD at site I.
768    * @param order Potential order.
769    */
770   protected abstract void chargeIPotentialAtK(PolarizableMultipoleSIMD mI, int order);
771 
772   /**
773    * Compute the induced dipole field components due to site I at site K.
774    *
775    * @param uxi   X-dipole component.
776    * @param uyi   Y-dipole component.
777    * @param uzi   Z-dipole component.
778    * @param order Potential order.
779    */
780   protected abstract void dipoleIPotentialAtK(DoubleVector uxi, DoubleVector uyi, DoubleVector uzi, int order);
781 
782   /**
783    * Compute the field components due to site I quadrupole at site K.
784    *
785    * @param mI    MultipoleTensorSIMD at site I.
786    * @param order Potential order.
787    */
788   protected abstract void quadrupoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order);
789 
790   /**
791    * The index is based on the idea of filling tetrahedron.
792    * <p>
793    * 1/r has an index of 0.
794    * <br>
795    * derivatives of x are first; indices from 1..o for d/dx..(d/dx)^o
796    * <br>
797    * derivatives of x and y are second; base triangle of size (o+1)(o+2)/2
798    * <br>
799    * derivatives of x, y and z are last; total size (o+1)*(o+2)*(o+3)/6
800    * <br>
801    * <p>
802    * This function is useful to set up masking constants:
803    * <br>
804    * static int Tlmn = ti(l,m,n,order)
805    * <br>
806    * For example the (d/dy)^2 (1/R) storage location:
807    * <br>
808    * static int T020 = ti(0,2,0,order)
809    *
810    * @param dx int The number of d/dx operations.
811    * @param dy int The number of d/dy operations.
812    * @param dz int The number of d/dz operations.
813    * @return int in the range (0..binomial(order + 3, 3) - 1)
814    */
815   protected final int ti(int dx, int dy, int dz) {
816     return MultipoleUtilities.ti(dx, dy, dz, order);
817   }
818 
819   // l + m + n = 0 (1)
820   /**
821    * No derivatives.
822    */
823   protected final int t000;
824   // l + m + n = 1 (3)   4
825   /**
826    * First derivative with respect to x.
827    */
828   protected final int t100;
829   /**
830    * First derivative with respect to y.
831    */
832   protected final int t010;
833   /**
834    * First derivative with respect to z.
835    */
836   protected final int t001;
837   // l + m + n = 2 (6)  10
838   /**
839    * Second derivative with respect to x.
840    */
841   protected final int t200;
842   /**
843    * Second derivative with respect to y.
844    */
845   protected final int t020;
846   /**
847    * Second derivative with respect to z.
848    */
849   protected final int t002;
850   /**
851    * Derivatives with respect to x and y.
852    */
853   protected final int t110;
854   /**
855    * Derivatives with respect to x and z.
856    */
857   protected final int t101;
858   /**
859    * Derivatives with respect to y and z.
860    */
861   protected final int t011;
862   // l + m + n = 3 (10) 20
863   /**
864    * Third derivative with respect to x.
865    */
866   protected final int t300;
867   /**
868    * Third derivative with respect to y.
869    */
870   protected final int t030;
871   /**
872    * Third derivative with respect to z.
873    */
874   protected final int t003;
875   /**
876    * Derivatives with respect to x and y.
877    */
878   protected final int t210;
879   /**
880    * Derivatives with respect to x and z.
881    */
882   protected final int t201;
883   /**
884    * Derivatives with respect to x and y.
885    */
886   protected final int t120;
887   /**
888    * Derivatives with respect to y and z.
889    */
890   protected final int t021;
891   /**
892    * Derivatives with respect to x and z.
893    */
894   protected final int t102;
895   /**
896    * Derivatives with respect to y and z.
897    */
898   protected final int t012;
899   /**
900    * Derivatives with respect to x, y and z.
901    */
902   protected final int t111;
903   // l + m + n = 4 (15) 35
904   /**
905    * Fourth derivative with respect to x.
906    */
907   protected final int t400;
908   /**
909    * Fourth derivative with respect to y.
910    */
911   protected final int t040;
912   /**
913    * Fourth derivative with respect to z.
914    */
915   protected final int t004;
916   /**
917    * Derivatives with respect to x and y.
918    */
919   protected final int t310;
920   /**
921    * Derivatives with respect to x and z.
922    */
923   protected final int t301;
924   /**
925    * Derivatives with respect to x and y.
926    */
927   protected final int t130;
928   /**
929    * Derivatives with respect to y and z.
930    */
931   protected final int t031;
932   /**
933    * Derivatives with respect to x and z.
934    */
935   protected final int t103;
936   /**
937    * Derivatives with respect to y and z.
938    */
939   protected final int t013;
940   /**
941    * Derivatives with respect to x and y.
942    */
943   protected final int t220;
944   /**
945    * Derivatives with respect to x and z.
946    */
947   protected final int t202;
948   /**
949    * Derivatives with respect to y and z.
950    */
951   protected final int t022;
952   /**
953    * Derivatives with respect to x, y and z.
954    */
955   protected final int t211;
956   /**
957    * Derivatives with respect to x, y and z.
958    */
959   protected final int t121;
960   /**
961    * Derivatives with respect to x, y and z.
962    */
963   protected final int t112;
964 
965   // l + m + n = 5 (21) 56
966   /**
967    * Fifth derivative with respect to x.
968    */
969   protected final int t500;
970   /**
971    * Fifth derivative with respect to y.
972    */
973   protected final int t050;
974   /**
975    * Fifth derivative with respect to z.
976    */
977   protected final int t005;
978   /**
979    * Derivatives with respect to x and y.
980    */
981   protected final int t410;
982   /**
983    * Derivatives with respect to x and z.
984    */
985   protected final int t401;
986   /**
987    * Derivatives with respect to x and y.
988    */
989   protected final int t140;
990   /**
991    * Derivatives with respect to y and z.
992    */
993   protected final int t041;
994   /**
995    * Derivatives with respect to x and z.
996    */
997   protected final int t104;
998   /**
999    * Derivatives with respect to y and z.
1000    */
1001   protected final int t014;
1002   /**
1003    * Derivatives with respect to x and y.
1004    */
1005   protected final int t320;
1006   /**
1007    * Derivatives with respect to x and z.
1008    */
1009   protected final int t302;
1010   /**
1011    * Derivatives with respect to x and y.
1012    */
1013   protected final int t230;
1014   /**
1015    * Derivatives with respect to y and z.
1016    */
1017   protected final int t032;
1018   /**
1019    * Derivatives with respect to x and z.
1020    */
1021   protected final int t203;
1022   /**
1023    * Derivatives with respect to y and z.
1024    */
1025   protected final int t023;
1026   /**
1027    * Derivatives with respect to x, y and z.
1028    */
1029   protected final int t311;
1030   /**
1031    * Derivatives with respect to x, y and z.
1032    */
1033   protected final int t131;
1034   /**
1035    * Derivatives with respect to x, y and z.
1036    */
1037   protected final int t113;
1038   /**
1039    * Derivatives with respect to x, y and z.
1040    */
1041   protected final int t221;
1042   /**
1043    * Derivatives with respect to x, y and z.
1044    */
1045   protected final int t212;
1046   /**
1047    * Derivatives with respect to x, y and z.
1048    */
1049   protected final int t122;
1050 
1051   // l + m + n = 6 (28) 84
1052   /**
1053    * Sixth derivative with respect to x.
1054    */
1055   protected final int t600;
1056   /**
1057    * Sixth derivative with respect to y.
1058    */
1059   protected final int t060;
1060   /**
1061    * Sixth derivative with respect to z.
1062    */
1063   protected final int t006;
1064   /**
1065    * Derivatives with respect to x and y.
1066    */
1067   protected final int t510;
1068   /**
1069    * Derivatives with respect to x and z.
1070    */
1071   protected final int t501;
1072   /**
1073    * Derivatives with respect to x and y.
1074    */
1075   protected final int t150;
1076   /**
1077    * Derivatives with respect to y and z.
1078    */
1079   protected final int t051;
1080   /**
1081    * Derivatives with respect to x and z.
1082    */
1083   protected final int t105;
1084   /**
1085    * Derivatives with respect to y and z.
1086    */
1087   protected final int t015;
1088   /**
1089    * Derivatives with respect to x and y.
1090    */
1091   protected final int t420;
1092   /**
1093    * Derivatives with respect to x and z.
1094    */
1095   protected final int t402;
1096   /**
1097    * Derivatives with respect to x and y.
1098    */
1099   protected final int t240;
1100   /**
1101    * Derivatives with respect to y and z.
1102    */
1103   protected final int t042;
1104   /**
1105    * Derivatives with respect to x and z.
1106    */
1107   protected final int t204;
1108   /**
1109    * Derivatives with respect to y and z.
1110    */
1111   protected final int t024;
1112   /**
1113    * Derivatives with respect to x, y and z.
1114    */
1115   protected final int t411;
1116   /**
1117    * Derivatives with respect to x, y and z.
1118    */
1119   protected final int t141;
1120   /**
1121    * Derivatives with respect to x, y and z.
1122    */
1123   protected final int t114;
1124   /**
1125    * Derivatives with respect to x and y.
1126    */
1127   protected final int t330;
1128   /**
1129    * Derivatives with respect to x and z.
1130    */
1131   protected final int t303;
1132   /**
1133    * Derivatives with respect to y and z.
1134    */
1135   protected final int t033;
1136   /**
1137    * Derivatives with respect to x, y and z.
1138    */
1139   protected final int t321;
1140   /**
1141    * Derivatives with respect to x, y and z.
1142    */
1143   protected final int t231;
1144   /**
1145    * Derivatives with respect to x, y and z.
1146    */
1147   protected final int t213;
1148   /**
1149    * Derivatives with respect to x, y and z.
1150    */
1151   protected final int t312;
1152   /**
1153    * Derivatives with respect to x, y and z.
1154    */
1155   protected final int t132;
1156   /**
1157    * Derivatives with respect to x, y and z.
1158    */
1159   protected final int t123;
1160   /**
1161    * Derivatives with respect to x, y and z.
1162    */
1163   protected final int t222;
1164 
1165 }