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 static ffx.numerics.multipole.CoordinateSystem.GLOBAL;
43  
44  /**
45   * The CoulombTensorGlobal class computes derivatives of 1/|<b>r</b>| via recursion to arbitrary
46   * order for Cartesian multipoles in the global frame using SIMD instructions.
47   *
48   * @author Michael J. Schnieders
49   * @see <a href="http://doi.org/10.1142/9789812830364_0002" target="_blank"> Matt Challacombe, Eric
50   * Schwegler and Jan Almlof, Modern developments in Hartree-Fock theory: Fast methods for
51   * computing the Coulomb matrix. Computational Chemistry: Review of Current Trends. pp. 53-107,
52   * Ed. J. Leczszynski, World Scientifc, 1996. </a>
53   * @since 1.0
54   */
55  public class CoulombTensorGlobalSIMD extends MultipoleTensorSIMD {
56  
57    /**
58     * Constructor for CoulombTensorGlobalSIMD.
59     *
60     * @param order The order of the tensor.
61     */
62    public CoulombTensorGlobalSIMD(int order) {
63      super(GLOBAL, order);
64      operator = Operator.COULOMB;
65    }
66  
67    @Override
68    public void setR(DoubleVector dx, DoubleVector dy, DoubleVector dz) {
69      x = dx;
70      y = dy;
71      z = dz;
72      DoubleVector x2 = x.mul(x);
73      DoubleVector y2 = y.mul(y);
74      DoubleVector z2 = z.mul(z);
75      r2 = x2.add(y2).add(z2);
76      R = r2.sqrt();
77    }
78  
79    @Override
80    protected void source(DoubleVector[] T000) {
81      DoubleVector ONE = DoubleVector.broadcast(R.species(), 1.0);
82      DoubleVector ir = ONE.div(R);
83      DoubleVector ir2 = ir.mul(ir);
84      for (int n = 0; n < o1; n++) {
85        T000[n] = ir.mul(coulombSource[n]);
86        ir = ir.mul(ir2);
87      }
88    }
89  
90    @Override
91    protected void order1() {
92      source(work);
93      DoubleVector term0000 = work[0];
94      DoubleVector term0001 = work[1];
95      R000 = term0000;
96      R100 = x.mul(term0001);
97      R010 = y.mul(term0001);
98      R001 = z.mul(term0001);
99    }
100 
101   @Override
102   protected void order2() {
103     source(work);
104     DoubleVector term0000 = work[0];
105     DoubleVector term0001 = work[1];
106     DoubleVector term0002 = work[2];
107     R000 = term0000;
108     R100 = x.mul(term0001);
109     DoubleVector term1001 = x.mul(term0002);
110     R200 = x.fma(term1001, term0001);
111     R010 = y.mul(term0001);
112     DoubleVector term0101 = y.mul(term0002);
113     R020 = y.fma(term0101, term0001);
114     R110 = y.mul(term1001);
115     R001 = z.mul(term0001);
116     DoubleVector term0011 = z.mul(term0002);
117     R002 = z.fma(term0011, term0001);
118     R011 = z.mul(term0101);
119     R101 = z.mul(term1001);
120   }
121 
122   @Override
123   protected void order3() {
124     source(work);
125     DoubleVector term0000 = work[0];
126     DoubleVector term0001 = work[1];
127     DoubleVector term0002 = work[2];
128     DoubleVector term0003 = work[3];
129     R000 = term0000;
130     R100 = x.mul(term0001);
131     DoubleVector term1001 = x.mul(term0002);
132     R200 = x.fma(term1001, term0001);
133     DoubleVector term1002 = x.mul(term0003);
134     DoubleVector term2001 = x.fma(term1002, term0002);
135     R300 = x.fma(term2001, term1001.mul(2));
136     R010 = y.mul(term0001);
137     DoubleVector term0101 = y.mul(term0002);
138     R020 = y.fma(term0101, term0001);
139     DoubleVector term0102 = y.mul(term0003);
140     DoubleVector term0201 = y.fma(term0102, term0002);
141     R030 = y.fma(term0201, term0101.mul(2));
142     R110 = y.mul(term1001);
143     DoubleVector term1101 = y.mul(term1002);
144     R120 = y.fma(term1101, term1001);
145     R210 = y.mul(term2001);
146     R001 = z.mul(term0001);
147     DoubleVector term0011 = z.mul(term0002);
148     R002 = z.fma(term0011, term0001);
149     DoubleVector term0012 = z.mul(term0003);
150     DoubleVector term0021 = z.fma(term0012, term0002);
151     R003 = z.fma(term0021, term0011.mul(2));
152     R011 = z.mul(term0101);
153     DoubleVector term0111 = z.mul(term0102);
154     R012 = z.fma(term0111, term0101);
155     R021 = z.mul(term0201);
156     R101 = z.mul(term1001);
157     DoubleVector term1011 = z.mul(term1002);
158     R102 = z.fma(term1011, term1001);
159     R111 = z.mul(term1101);
160     R201 = z.mul(term2001);
161   }
162 
163   @Override
164   protected void order4() {
165     source(work);
166     DoubleVector term0000 = work[0];
167     DoubleVector term0001 = work[1];
168     DoubleVector term0002 = work[2];
169     DoubleVector term0003 = work[3];
170     DoubleVector term0004 = work[4];
171     R000 = term0000;
172     R100 = x.mul(term0001);
173     DoubleVector term1001 = x.mul(term0002);
174     R200 = x.fma(term1001, term0001);
175     DoubleVector term1002 = x.mul(term0003);
176     DoubleVector term2001 = x.fma(term1002, term0002);
177     R300 = x.fma(term2001, term1001.mul(2));
178     DoubleVector term1003 = x.mul(term0004);
179     DoubleVector term2002 = x.fma(term1003, term0003);
180     DoubleVector term3001 = x.fma(term2002, term1002.mul(2));
181     R400 = x.fma(term3001, term2001.mul(3));
182     R010 = y.mul(term0001);
183     DoubleVector term0101 = y.mul(term0002);
184     R020 = y.fma(term0101, term0001);
185     DoubleVector term0102 = y.mul(term0003);
186     DoubleVector term0201 = y.fma(term0102, term0002);
187     R030 = y.fma(term0201, term0101.mul(2));
188     DoubleVector term0103 = y.mul(term0004);
189     DoubleVector term0202 = y.fma(term0103, term0003);
190     DoubleVector term0301 = y.fma(term0202, term0102.mul(2));
191     R040 = y.fma(term0301, term0201.mul(3));
192     R110 = y.mul(term1001);
193     DoubleVector term1101 = y.mul(term1002);
194     R120 = y.fma(term1101, term1001);
195     DoubleVector term1102 = y.mul(term1003);
196     DoubleVector term1201 = y.fma(term1102, term1002);
197     R130 = y.fma(term1201, term1101.mul(2));
198     R210 = y.mul(term2001);
199     DoubleVector term2101 = y.mul(term2002);
200     R220 = y.fma(term2101, term2001);
201     R310 = y.mul(term3001);
202     R001 = z.mul(term0001);
203     DoubleVector term0011 = z.mul(term0002);
204     R002 = z.fma(term0011, term0001);
205     DoubleVector term0012 = z.mul(term0003);
206     DoubleVector term0021 = z.fma(term0012, term0002);
207     R003 = z.fma(term0021, term0011.mul(2));
208     DoubleVector term0013 = z.mul(term0004);
209     DoubleVector term0022 = z.fma(term0013, term0003);
210     DoubleVector term0031 = z.fma(term0022, term0012.mul(2));
211     R004 = z.fma(term0031, term0021.mul(3));
212     R011 = z.mul(term0101);
213     DoubleVector term0111 = z.mul(term0102);
214     R012 = z.fma(term0111, term0101);
215     DoubleVector term0112 = z.mul(term0103);
216     DoubleVector term0121 = z.fma(term0112, term0102);
217     R013 = z.fma(term0121, term0111.mul(2));
218     R021 = z.mul(term0201);
219     DoubleVector term0211 = z.mul(term0202);
220     R022 = z.fma(term0211, term0201);
221     R031 = z.mul(term0301);
222     R101 = z.mul(term1001);
223     DoubleVector term1011 = z.mul(term1002);
224     R102 = z.fma(term1011, term1001);
225     DoubleVector term1012 = z.mul(term1003);
226     DoubleVector term1021 = z.fma(term1012, term1002);
227     R103 = z.fma(term1021, term1011.mul(2));
228     R111 = z.mul(term1101);
229     DoubleVector term1111 = z.mul(term1102);
230     R112 = z.fma(term1111, term1101);
231     R121 = z.mul(term1201);
232     R201 = z.mul(term2001);
233     DoubleVector term2011 = z.mul(term2002);
234     R202 = z.fma(term2011, term2001);
235     R211 = z.mul(term2101);
236     R301 = z.mul(term3001);
237   }
238 
239   @Override
240   protected void order5() {
241     source(work);
242     DoubleVector term0000 = work[0];
243     DoubleVector term0001 = work[1];
244     DoubleVector term0002 = work[2];
245     DoubleVector term0003 = work[3];
246     DoubleVector term0004 = work[4];
247     DoubleVector term0005 = work[5];
248     R000 = term0000;
249     R100 = x.mul(term0001);
250     DoubleVector term1001 = x.mul(term0002);
251     R200 = x.fma(term1001, term0001);
252     DoubleVector term1002 = x.mul(term0003);
253     DoubleVector term2001 = x.fma(term1002, term0002);
254     R300 = x.fma(term2001, term1001.mul(2));
255     DoubleVector term1003 = x.mul(term0004);
256     DoubleVector term2002 = x.fma(term1003, term0003);
257     DoubleVector term3001 = x.fma(term2002, term1002.mul(2));
258     R400 = x.fma(term3001, term2001.mul(3));
259     DoubleVector term1004 = x.mul(term0005);
260     DoubleVector term2003 = x.fma(term1004, term0004);
261     DoubleVector term3002 = x.fma(term2003, term1003.mul(2));
262     DoubleVector term4001 = x.fma(term3002, term2002.mul(3));
263     R500 = x.fma(term4001, term3001.mul(4));
264     R010 = y.mul(term0001);
265     DoubleVector term0101 = y.mul(term0002);
266     R020 = y.fma(term0101, term0001);
267     DoubleVector term0102 = y.mul(term0003);
268     DoubleVector term0201 = y.fma(term0102, term0002);
269     R030 = y.fma(term0201, term0101.mul(2));
270     DoubleVector term0103 = y.mul(term0004);
271     DoubleVector term0202 = y.fma(term0103, term0003);
272     DoubleVector term0301 = y.fma(term0202, term0102.mul(2));
273     R040 = y.fma(term0301, term0201.mul(3));
274     DoubleVector term0104 = y.mul(term0005);
275     DoubleVector term0203 = y.fma(term0104, term0004);
276     DoubleVector term0302 = y.fma(term0203, term0103.mul(2));
277     DoubleVector term0401 = y.fma(term0302, term0202.mul(3));
278     R050 = y.fma(term0401, term0301.mul(4));
279     R110 = y.mul(term1001);
280     DoubleVector term1101 = y.mul(term1002);
281     R120 = y.fma(term1101, term1001);
282     DoubleVector term1102 = y.mul(term1003);
283     DoubleVector term1201 = y.fma(term1102, term1002);
284     R130 = y.fma(term1201, term1101.mul(2));
285     DoubleVector term1103 = y.mul(term1004);
286     DoubleVector term1202 = y.fma(term1103, term1003);
287     DoubleVector term1301 = y.fma(term1202, term1102.mul(2));
288     R140 = y.fma(term1301, term1201.mul(3));
289     R210 = y.mul(term2001);
290     DoubleVector term2101 = y.mul(term2002);
291     R220 = y.fma(term2101, term2001);
292     DoubleVector term2102 = y.mul(term2003);
293     DoubleVector term2201 = y.fma(term2102, term2002);
294     R230 = y.fma(term2201, term2101.mul(2));
295     R310 = y.mul(term3001);
296     DoubleVector term3101 = y.mul(term3002);
297     R320 = y.fma(term3101, term3001);
298     R410 = y.mul(term4001);
299     R001 = z.mul(term0001);
300     DoubleVector term0011 = z.mul(term0002);
301     R002 = z.fma(term0011, term0001);
302     DoubleVector term0012 = z.mul(term0003);
303     DoubleVector term0021 = z.fma(term0012, term0002);
304     R003 = z.fma(term0021, term0011.mul(2));
305     DoubleVector term0013 = z.mul(term0004);
306     DoubleVector term0022 = z.fma(term0013, term0003);
307     DoubleVector term0031 = z.fma(term0022, term0012.mul(2));
308     R004 = z.fma(term0031, term0021.mul(3));
309     DoubleVector term0014 = z.mul(term0005);
310     DoubleVector term0023 = z.fma(term0014, term0004);
311     DoubleVector term0032 = z.fma(term0023, term0013.mul(2));
312     DoubleVector term0041 = z.fma(term0032, term0022.mul(3));
313     R005 = z.fma(term0041, term0031.mul(4));
314     R011 = z.mul(term0101);
315     DoubleVector term0111 = z.mul(term0102);
316     R012 = z.fma(term0111, term0101);
317     DoubleVector term0112 = z.mul(term0103);
318     DoubleVector term0121 = z.fma(term0112, term0102);
319     R013 = z.fma(term0121, term0111.mul(2));
320     DoubleVector term0113 = z.mul(term0104);
321     DoubleVector term0122 = z.fma(term0113, term0103);
322     DoubleVector term0131 = z.fma(term0122, term0112.mul(2));
323     R014 = z.fma(term0131, term0121.mul(3));
324     R021 = z.mul(term0201);
325     DoubleVector term0211 = z.mul(term0202);
326     R022 = z.fma(term0211, term0201);
327     DoubleVector term0212 = z.mul(term0203);
328     DoubleVector term0221 = z.fma(term0212, term0202);
329     R023 = z.fma(term0221, term0211.mul(2));
330     R031 = z.mul(term0301);
331     DoubleVector term0311 = z.mul(term0302);
332     R032 = z.fma(term0311, term0301);
333     R041 = z.mul(term0401);
334     R101 = z.mul(term1001);
335     DoubleVector term1011 = z.mul(term1002);
336     R102 = z.fma(term1011, term1001);
337     DoubleVector term1012 = z.mul(term1003);
338     DoubleVector term1021 = z.fma(term1012, term1002);
339     R103 = z.fma(term1021, term1011.mul(2));
340     DoubleVector term1013 = z.mul(term1004);
341     DoubleVector term1022 = z.fma(term1013, term1003);
342     DoubleVector term1031 = z.fma(term1022, term1012.mul(2));
343     R104 = z.fma(term1031, term1021.mul(3));
344     R111 = z.mul(term1101);
345     DoubleVector term1111 = z.mul(term1102);
346     R112 = z.fma(term1111, term1101);
347     DoubleVector term1112 = z.mul(term1103);
348     DoubleVector term1121 = z.fma(term1112, term1102);
349     R113 = z.fma(term1121, term1111.mul(2));
350     R121 = z.mul(term1201);
351     DoubleVector term1211 = z.mul(term1202);
352     R122 = z.fma(term1211, term1201);
353     R131 = z.mul(term1301);
354     R201 = z.mul(term2001);
355     DoubleVector term2011 = z.mul(term2002);
356     R202 = z.fma(term2011, term2001);
357     DoubleVector term2012 = z.mul(term2003);
358     DoubleVector term2021 = z.fma(term2012, term2002);
359     R203 = z.fma(term2021, term2011.mul(2));
360     R211 = z.mul(term2101);
361     DoubleVector term2111 = z.mul(term2102);
362     R212 = z.fma(term2111, term2101);
363     R221 = z.mul(term2201);
364     R301 = z.mul(term3001);
365     DoubleVector term3011 = z.mul(term3002);
366     R302 = z.fma(term3011, term3001);
367     R311 = z.mul(term3101);
368     R401 = z.mul(term4001);
369   }
370 
371   @Override
372   protected void order6() {
373     source(work);
374     DoubleVector term0000 = work[0];
375     DoubleVector term0001 = work[1];
376     DoubleVector term0002 = work[2];
377     DoubleVector term0003 = work[3];
378     DoubleVector term0004 = work[4];
379     DoubleVector term0005 = work[5];
380     DoubleVector term0006 = work[6];
381     R000 = term0000;
382     R100 = x.mul(term0001);
383     DoubleVector term1001 = x.mul(term0002);
384     R200 = x.fma(term1001, term0001);
385     DoubleVector term1002 = x.mul(term0003);
386     DoubleVector term2001 = x.fma(term1002, term0002);
387     R300 = x.fma(term2001, term1001.mul(2));
388     DoubleVector term1003 = x.mul(term0004);
389     DoubleVector term2002 = x.fma(term1003, term0003);
390     DoubleVector term3001 = x.fma(term2002, term1002.mul(2));
391     R400 = x.fma(term3001, term2001.mul(3));
392     DoubleVector term1004 = x.mul(term0005);
393     DoubleVector term2003 = x.fma(term1004, term0004);
394     DoubleVector term3002 = x.fma(term2003, term1003.mul(2));
395     DoubleVector term4001 = x.fma(term3002, term2002.mul(3));
396     R500 = x.fma(term4001, term3001.mul(4));
397     DoubleVector term1005 = x.mul(term0006);
398     DoubleVector term2004 = x.fma(term1005, term0005);
399     DoubleVector term3003 = x.fma(term2004, term1004.mul(2));
400     DoubleVector term4002 = x.fma(term3003, term2003.mul(3));
401     DoubleVector term5001 = x.fma(term4002, term3002.mul(4));
402     R600 = x.fma(term5001, term4001.mul(5));
403     R010 = y.mul(term0001);
404     DoubleVector term0101 = y.mul(term0002);
405     R020 = y.fma(term0101, term0001);
406     DoubleVector term0102 = y.mul(term0003);
407     DoubleVector term0201 = y.fma(term0102, term0002);
408     R030 = y.fma(term0201, term0101.mul(2));
409     DoubleVector term0103 = y.mul(term0004);
410     DoubleVector term0202 = y.fma(term0103, term0003);
411     DoubleVector term0301 = y.fma(term0202, term0102.mul(2));
412     R040 = y.fma(term0301, term0201.mul(3));
413     DoubleVector term0104 = y.mul(term0005);
414     DoubleVector term0203 = y.fma(term0104, term0004);
415     DoubleVector term0302 = y.fma(term0203, term0103.mul(2));
416     DoubleVector term0401 = y.fma(term0302, term0202.mul(3));
417     R050 = y.fma(term0401, term0301.mul(4));
418     DoubleVector term0105 = y.mul(term0006);
419     DoubleVector term0204 = y.fma(term0105, term0005);
420     DoubleVector term0303 = y.fma(term0204, term0104.mul(2));
421     DoubleVector term0402 = y.fma(term0303, term0203.mul(3));
422     DoubleVector term0501 = y.fma(term0402, term0302.mul(4));
423     R060 = y.fma(term0501, term0401.mul(5));
424     R110 = y.mul(term1001);
425     DoubleVector term1101 = y.mul(term1002);
426     R120 = y.fma(term1101, term1001);
427     DoubleVector term1102 = y.mul(term1003);
428     DoubleVector term1201 = y.fma(term1102, term1002);
429     R130 = y.fma(term1201, term1101.mul(2));
430     DoubleVector term1103 = y.mul(term1004);
431     DoubleVector term1202 = y.fma(term1103, term1003);
432     DoubleVector term1301 = y.fma(term1202, term1102.mul(2));
433     R140 = y.fma(term1301, term1201.mul(3));
434     DoubleVector term1104 = y.mul(term1005);
435     DoubleVector term1203 = y.fma(term1104, term1004);
436     DoubleVector term1302 = y.fma(term1203, term1103.mul(2));
437     DoubleVector term1401 = y.fma(term1302, term1202.mul(3));
438     R150 = y.fma(term1401, term1301.mul(4));
439     R210 = y.mul(term2001);
440     DoubleVector term2101 = y.mul(term2002);
441     R220 = y.fma(term2101, term2001);
442     DoubleVector term2102 = y.mul(term2003);
443     DoubleVector term2201 = y.fma(term2102, term2002);
444     R230 = y.fma(term2201, term2101.mul(2));
445     DoubleVector term2103 = y.mul(term2004);
446     DoubleVector term2202 = y.fma(term2103, term2003);
447     DoubleVector term2301 = y.fma(term2202, term2102.mul(2));
448     R240 = y.fma(term2301, term2201.mul(3));
449     R310 = y.mul(term3001);
450     DoubleVector term3101 = y.mul(term3002);
451     R320 = y.fma(term3101, term3001);
452     DoubleVector term3102 = y.mul(term3003);
453     DoubleVector term3201 = y.fma(term3102, term3002);
454     R330 = y.fma(term3201, term3101.mul(2));
455     R410 = y.mul(term4001);
456     DoubleVector term4101 = y.mul(term4002);
457     R420 = y.fma(term4101, term4001);
458     R510 = y.mul(term5001);
459     R001 = z.mul(term0001);
460     DoubleVector term0011 = z.mul(term0002);
461     R002 = z.fma(term0011, term0001);
462     DoubleVector term0012 = z.mul(term0003);
463     DoubleVector term0021 = z.fma(term0012, term0002);
464     R003 = z.fma(term0021, term0011.mul(2));
465     DoubleVector term0013 = z.mul(term0004);
466     DoubleVector term0022 = z.fma(term0013, term0003);
467     DoubleVector term0031 = z.fma(term0022, term0012.mul(2));
468     R004 = z.fma(term0031, term0021.mul(3));
469     DoubleVector term0014 = z.mul(term0005);
470     DoubleVector term0023 = z.fma(term0014, term0004);
471     DoubleVector term0032 = z.fma(term0023, term0013.mul(2));
472     DoubleVector term0041 = z.fma(term0032, term0022.mul(3));
473     R005 = z.fma(term0041, term0031.mul(4));
474     DoubleVector term0015 = z.mul(term0006);
475     DoubleVector term0024 = z.fma(term0015, term0005);
476     DoubleVector term0033 = z.fma(term0024, term0014.mul(2));
477     DoubleVector term0042 = z.fma(term0033, term0023.mul(3));
478     DoubleVector term0051 = z.fma(term0042, term0032.mul(4));
479     R006 = z.fma(term0051, term0041.mul(5));
480     R011 = z.mul(term0101);
481     DoubleVector term0111 = z.mul(term0102);
482     R012 = z.fma(term0111, term0101);
483     DoubleVector term0112 = z.mul(term0103);
484     DoubleVector term0121 = z.fma(term0112, term0102);
485     R013 = z.fma(term0121, term0111.mul(2));
486     DoubleVector term0113 = z.mul(term0104);
487     DoubleVector term0122 = z.fma(term0113, term0103);
488     DoubleVector term0131 = z.fma(term0122, term0112.mul(2));
489     R014 = z.fma(term0131, term0121.mul(3));
490     DoubleVector term0114 = z.mul(term0105);
491     DoubleVector term0123 = z.fma(term0114, term0104);
492     DoubleVector term0132 = z.fma(term0123, term0113.mul(2));
493     DoubleVector term0141 = z.fma(term0132, term0122.mul(3));
494     R015 = z.fma(term0141, term0131.mul(4));
495     R021 = z.mul(term0201);
496     DoubleVector term0211 = z.mul(term0202);
497     R022 = z.fma(term0211, term0201);
498     DoubleVector term0212 = z.mul(term0203);
499     DoubleVector term0221 = z.fma(term0212, term0202);
500     R023 = z.fma(term0221, term0211.mul(2));
501     DoubleVector term0213 = z.mul(term0204);
502     DoubleVector term0222 = z.fma(term0213, term0203);
503     DoubleVector term0231 = z.fma(term0222, term0212.mul(2));
504     R024 = z.fma(term0231, term0221.mul(3));
505     R031 = z.mul(term0301);
506     DoubleVector term0311 = z.mul(term0302);
507     R032 = z.fma(term0311, term0301);
508     DoubleVector term0312 = z.mul(term0303);
509     DoubleVector term0321 = z.fma(term0312, term0302);
510     R033 = z.fma(term0321, term0311.mul(2));
511     R041 = z.mul(term0401);
512     DoubleVector term0411 = z.mul(term0402);
513     R042 = z.fma(term0411, term0401);
514     R051 = z.mul(term0501);
515     R101 = z.mul(term1001);
516     DoubleVector term1011 = z.mul(term1002);
517     R102 = z.fma(term1011, term1001);
518     DoubleVector term1012 = z.mul(term1003);
519     DoubleVector term1021 = z.fma(term1012, term1002);
520     R103 = z.fma(term1021, term1011.mul(2));
521     DoubleVector term1013 = z.mul(term1004);
522     DoubleVector term1022 = z.fma(term1013, term1003);
523     DoubleVector term1031 = z.fma(term1022, term1012.mul(2));
524     R104 = z.fma(term1031, term1021.mul(3));
525     DoubleVector term1014 = z.mul(term1005);
526     DoubleVector term1023 = z.fma(term1014, term1004);
527     DoubleVector term1032 = z.fma(term1023, term1013.mul(2));
528     DoubleVector term1041 = z.fma(term1032, term1022.mul(3));
529     R105 = z.fma(term1041, term1031.mul(4));
530     R111 = z.mul(term1101);
531     DoubleVector term1111 = z.mul(term1102);
532     R112 = z.fma(term1111, term1101);
533     DoubleVector term1112 = z.mul(term1103);
534     DoubleVector term1121 = z.fma(term1112, term1102);
535     R113 = z.fma(term1121, term1111.mul(2));
536     DoubleVector term1113 = z.mul(term1104);
537     DoubleVector term1122 = z.fma(term1113, term1103);
538     DoubleVector term1131 = z.fma(term1122, term1112.mul(2));
539     R114 = z.fma(term1131, term1121.mul(3));
540     R121 = z.mul(term1201);
541     DoubleVector term1211 = z.mul(term1202);
542     R122 = z.fma(term1211, term1201);
543     DoubleVector term1212 = z.mul(term1203);
544     DoubleVector term1221 = z.fma(term1212, term1202);
545     R123 = z.fma(term1221, term1211.mul(2));
546     R131 = z.mul(term1301);
547     DoubleVector term1311 = z.mul(term1302);
548     R132 = z.fma(term1311, term1301);
549     R141 = z.mul(term1401);
550     R201 = z.mul(term2001);
551     DoubleVector term2011 = z.mul(term2002);
552     R202 = z.fma(term2011, term2001);
553     DoubleVector term2012 = z.mul(term2003);
554     DoubleVector term2021 = z.fma(term2012, term2002);
555     R203 = z.fma(term2021, term2011.mul(2));
556     DoubleVector term2013 = z.mul(term2004);
557     DoubleVector term2022 = z.fma(term2013, term2003);
558     DoubleVector term2031 = z.fma(term2022, term2012.mul(2));
559     R204 = z.fma(term2031, term2021.mul(3));
560     R211 = z.mul(term2101);
561     DoubleVector term2111 = z.mul(term2102);
562     R212 = z.fma(term2111, term2101);
563     DoubleVector term2112 = z.mul(term2103);
564     DoubleVector term2121 = z.fma(term2112, term2102);
565     R213 = z.fma(term2121, term2111.mul(2));
566     R221 = z.mul(term2201);
567     DoubleVector term2211 = z.mul(term2202);
568     R222 = z.fma(term2211, term2201);
569     R231 = z.mul(term2301);
570     R301 = z.mul(term3001);
571     DoubleVector term3011 = z.mul(term3002);
572     R302 = z.fma(term3011, term3001);
573     DoubleVector term3012 = z.mul(term3003);
574     DoubleVector term3021 = z.fma(term3012, term3002);
575     R303 = z.fma(term3021, term3011.mul(2));
576     R311 = z.mul(term3101);
577     DoubleVector term3111 = z.mul(term3102);
578     R312 = z.fma(term3111, term3101);
579     R321 = z.mul(term3201);
580     R401 = z.mul(term4001);
581     DoubleVector term4011 = z.mul(term4002);
582     R402 = z.fma(term4011, term4001);
583     R411 = z.mul(term4101);
584     R501 = z.mul(term5001);
585   }
586 
587   /**
588    * {@inheritDoc}
589    */
590   @SuppressWarnings("fallthrough")
591   @Override
592   protected void multipoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order) {
593     switch (order) {
594       default:
595       case 3:
596         // Order 3 (need a minus sign).
597         DoubleVector term300 = mK.q.mul(R300);
598         term300 = mK.dx.fma(R400, term300);
599         term300 = mK.dy.fma(R310, term300);
600         term300 = mK.dz.fma(R301, term300);
601         term300 = mK.qxx.fma(R500, term300);
602         term300 = mK.qyy.fma(R320, term300);
603         term300 = mK.qzz.fma(R302, term300);
604         term300 = mK.qxy.fma(R410, term300);
605         term300 = mK.qxz.fma(R401, term300);
606         term300 = mK.qyz.fma(R311, term300);
607         E300 = term300.neg();
608         DoubleVector term030 = mK.q.mul(R030);
609         term030 = mK.dx.fma(R130, term030);
610         term030 = mK.dy.fma(R040, term030);
611         term030 = mK.dz.fma(R031, term030);
612         term030 = mK.qxx.fma(R230, term030);
613         term030 = mK.qyy.fma(R050, term030);
614         term030 = mK.qzz.fma(R032, term030);
615         term030 = mK.qxy.fma(R140, term030);
616         term030 = mK.qxz.fma(R131, term030);
617         term030 = mK.qyz.fma(R041, term030);
618         E030 = term030.neg();
619         DoubleVector term003 = mK.q.mul(R003);
620         term003 = mK.dx.fma(R103, term003);
621         term003 = mK.dy.fma(R013, term003);
622         term003 = mK.dz.fma(R004, term003);
623         term003 = mK.qxx.fma(R203, term003);
624         term003 = mK.qyy.fma(R023, term003);
625         term003 = mK.qzz.fma(R005, term003);
626         term003 = mK.qxy.fma(R113, term003);
627         term003 = mK.qxz.fma(R104, term003);
628         term003 = mK.qyz.fma(R014, term003);
629         E003 = term003.neg();
630         DoubleVector term210 = mK.q.mul(R210);
631         term210 = mK.dx.fma(R310, term210);
632         term210 = mK.dy.fma(R220, term210);
633         term210 = mK.dz.fma(R211, term210);
634         term210 = mK.qxx.fma(R410, term210);
635         term210 = mK.qyy.fma(R230, term210);
636         term210 = mK.qzz.fma(R212, term210);
637         term210 = mK.qxy.fma(R320, term210);
638         term210 = mK.qxz.fma(R311, term210);
639         term210 = mK.qyz.fma(R221, term210);
640         E210 = term210.neg();
641         DoubleVector term201 = mK.q.mul(R201);
642         term201 = mK.dx.fma(R301, term201);
643         term201 = mK.dy.fma(R211, term201);
644         term201 = mK.dz.fma(R202, term201);
645         term201 = mK.qxx.fma(R401, term201);
646         term201 = mK.qyy.fma(R221, term201);
647         term201 = mK.qzz.fma(R203, term201);
648         term201 = mK.qxy.fma(R311, term201);
649         term201 = mK.qxz.fma(R302, term201);
650         term201 = mK.qyz.fma(R212, term201);
651         E201 = term201.neg();
652         DoubleVector term120 = mK.q.mul(R120);
653         term120 = mK.dx.fma(R220, term120);
654         term120 = mK.dy.fma(R130, term120);
655         term120 = mK.dz.fma(R121, term120);
656         term120 = mK.qxx.fma(R320, term120);
657         term120 = mK.qyy.fma(R140, term120);
658         term120 = mK.qzz.fma(R122, term120);
659         term120 = mK.qxy.fma(R230, term120);
660         term120 = mK.qxz.fma(R221, term120);
661         term120 = mK.qyz.fma(R131, term120);
662         E120 = term120.neg();
663         DoubleVector term021 = mK.q.mul(R021);
664         term021 = mK.dx.fma(R121, term021);
665         term021 = mK.dy.fma(R031, term021);
666         term021 = mK.dz.fma(R022, term021);
667         term021 = mK.qxx.fma(R221, term021);
668         term021 = mK.qyy.fma(R041, term021);
669         term021 = mK.qzz.fma(R023, term021);
670         term021 = mK.qxy.fma(R131, term021);
671         term021 = mK.qxz.fma(R122, term021);
672         term021 = mK.qyz.fma(R032, term021);
673         E021 = term021.neg();
674         DoubleVector term102 = mK.q.mul(R102);
675         term102 = mK.dx.fma(R202, term102);
676         term102 = mK.dy.fma(R112, term102);
677         term102 = mK.dz.fma(R103, term102);
678         term102 = mK.qxx.fma(R302, term102);
679         term102 = mK.qyy.fma(R122, term102);
680         term102 = mK.qzz.fma(R104, term102);
681         term102 = mK.qxy.fma(R212, term102);
682         term102 = mK.qxz.fma(R203, term102);
683         term102 = mK.qyz.fma(R113, term102);
684         E102 = term102.neg();
685         DoubleVector term012 = mK.q.mul(R012);
686         term012 = mK.dx.fma(R112, term012);
687         term012 = mK.dy.fma(R022, term012);
688         term012 = mK.dz.fma(R013, term012);
689         term012 = mK.qxx.fma(R212, term012);
690         term012 = mK.qyy.fma(R032, term012);
691         term012 = mK.qzz.fma(R014, term012);
692         term012 = mK.qxy.fma(R122, term012);
693         term012 = mK.qxz.fma(R113, term012);
694         term012 = mK.qyz.fma(R023, term012);
695         E012 = term012.neg();
696         DoubleVector term111 = mK.q.mul(R111);
697         term111 = mK.dx.fma(R211, term111);
698         term111 = mK.dy.fma(R121, term111);
699         term111 = mK.dz.fma(R112, term111);
700         term111 = mK.qxx.fma(R311, term111);
701         term111 = mK.qyy.fma(R131, term111);
702         term111 = mK.qzz.fma(R113, term111);
703         term111 = mK.qxy.fma(R221, term111);
704         term111 = mK.qxz.fma(R212, term111);
705         term111 = mK.qyz.fma(R122, term111);
706         E111 = term111.neg();
707       case 2:
708         // Order 2.
709         DoubleVector term200 = mK.q.mul(R200);
710         term200 = mK.dx.fma(R300, term200);
711         term200 = mK.dy.fma(R210, term200);
712         term200 = mK.dz.fma(R201, term200);
713         term200 = mK.qxx.fma(R400, term200);
714         term200 = mK.qyy.fma(R220, term200);
715         term200 = mK.qzz.fma(R202, term200);
716         term200 = mK.qxy.fma(R310, term200);
717         term200 = mK.qxz.fma(R301, term200);
718         term200 = mK.qyz.fma(R211, term200);
719         E200 = term200;
720         DoubleVector term020 = mK.q.mul(R020);
721         term020 = mK.dx.fma(R120, term020);
722         term020 = mK.dy.fma(R030, term020);
723         term020 = mK.dz.fma(R021, term020);
724         term020 = mK.qxx.fma(R220, term020);
725         term020 = mK.qyy.fma(R040, term020);
726         term020 = mK.qzz.fma(R022, term020);
727         term020 = mK.qxy.fma(R130, term020);
728         term020 = mK.qxz.fma(R121, term020);
729         term020 = mK.qyz.fma(R031, term020);
730         E020 = term020;
731         DoubleVector term002 = mK.q.mul(R002);
732         term002 = mK.dx.fma(R102, term002);
733         term002 = mK.dy.fma(R012, term002);
734         term002 = mK.dz.fma(R003, term002);
735         term002 = mK.qxx.fma(R202, term002);
736         term002 = mK.qyy.fma(R022, term002);
737         term002 = mK.qzz.fma(R004, term002);
738         term002 = mK.qxy.fma(R112, term002);
739         term002 = mK.qxz.fma(R103, term002);
740         term002 = mK.qyz.fma(R013, term002);
741         E002 = term002;
742         DoubleVector term110 = mK.q.mul(R110);
743         term110 = mK.dx.fma(R210, term110);
744         term110 = mK.dy.fma(R120, term110);
745         term110 = mK.dz.fma(R111, term110);
746         term110 = mK.qxx.fma(R310, term110);
747         term110 = mK.qyy.fma(R130, term110);
748         term110 = mK.qzz.fma(R112, term110);
749         term110 = mK.qxy.fma(R220, term110);
750         term110 = mK.qxz.fma(R211, term110);
751         term110 = mK.qyz.fma(R121, term110);
752         E110 = term110;
753         DoubleVector term101 = mK.q.mul(R101);
754         term101 = mK.dx.fma(R201, term101);
755         term101 = mK.dy.fma(R111, term101);
756         term101 = mK.dz.fma(R102, term101);
757         term101 = mK.qxx.fma(R301, term101);
758         term101 = mK.qyy.fma(R121, term101);
759         term101 = mK.qzz.fma(R103, term101);
760         term101 = mK.qxy.fma(R211, term101);
761         term101 = mK.qxz.fma(R202, term101);
762         term101 = mK.qyz.fma(R112, term101);
763         E101 = term101;
764         DoubleVector term011 = mK.q.mul(R011);
765         term011 = mK.dx.fma(R111, term011);
766         term011 = mK.dy.fma(R021, term011);
767         term011 = mK.dz.fma(R012, term011);
768         term011 = mK.qxx.fma(R211, term011);
769         term011 = mK.qyy.fma(R031, term011);
770         term011 = mK.qzz.fma(R013, term011);
771         term011 = mK.qxy.fma(R121, term011);
772         term011 = mK.qxz.fma(R112, term011);
773         term011 = mK.qyz.fma(R022, term011);
774         E011 = term011;
775       case 1:
776         // Order 1 (need a minus sign).
777         DoubleVector term100 = mK.q.mul(R100);
778         term100 = mK.dx.fma(R200, term100);
779         term100 = mK.dy.fma(R110, term100);
780         term100 = mK.dz.fma(R101, term100);
781         term100 = mK.qxx.fma(R300, term100);
782         term100 = mK.qyy.fma(R120, term100);
783         term100 = mK.qzz.fma(R102, term100);
784         term100 = mK.qxy.fma(R210, term100);
785         term100 = mK.qxz.fma(R201, term100);
786         term100 = mK.qyz.fma(R111, term100);
787         E100 = term100.neg();
788         DoubleVector term010 = mK.q.mul(R010);
789         term010 = mK.dx.fma(R110, term010);
790         term010 = mK.dy.fma(R020, term010);
791         term010 = mK.dz.fma(R011, term010);
792         term010 = mK.qxx.fma(R210, term010);
793         term010 = mK.qyy.fma(R030, term010);
794         term010 = mK.qzz.fma(R012, term010);
795         term010 = mK.qxy.fma(R120, term010);
796         term010 = mK.qxz.fma(R111, term010);
797         term010 = mK.qyz.fma(R021, term010);
798         E010 = term010.neg();
799         DoubleVector term001 = mK.q.mul(R001);
800         term001 = mK.dx.fma(R101, term001);
801         term001 = mK.dy.fma(R011, term001);
802         term001 = mK.dz.fma(R002, term001);
803         term001 = mK.qxx.fma(R201, term001);
804         term001 = mK.qyy.fma(R021, term001);
805         term001 = mK.qzz.fma(R003, term001);
806         term001 = mK.qxy.fma(R111, term001);
807         term001 = mK.qxz.fma(R102, term001);
808         term001 = mK.qyz.fma(R012, term001);
809         E001 = term001.neg();
810       case 0:
811         DoubleVector term000 = mK.q.mul(R000);
812         term000 = mK.dx.fma(R100, term000);
813         term000 = mK.dy.fma(R010, term000);
814         term000 = mK.dz.fma(R001, term000);
815         term000 = mK.qxx.fma(R200, term000);
816         term000 = mK.qyy.fma(R020, term000);
817         term000 = mK.qzz.fma(R002, term000);
818         term000 = mK.qxy.fma(R110, term000);
819         term000 = mK.qxz.fma(R101, term000);
820         term000 = mK.qyz.fma(R011, term000);
821         E000 = term000;
822     }
823   }
824 
825   /**
826    * {@inheritDoc}
827    */
828   @SuppressWarnings("fallthrough")
829   @Override
830   protected void chargeKPotentialAtI(PolarizableMultipoleSIMD mK, int order) {
831     switch (order) {
832       default:
833       case 3:
834         // Order 3 (need a minus sign).
835         E300 = mK.q.mul(R300).neg();
836         E030 = mK.q.mul(R030).neg();
837         E003 = mK.q.mul(R003).neg();
838         E210 = mK.q.mul(R210).neg();
839         E201 = mK.q.mul(R201).neg();
840         E120 = mK.q.mul(R120).neg();
841         E021 = mK.q.mul(R021).neg();
842         E102 = mK.q.mul(R102).neg();
843         E012 = mK.q.mul(R012).neg();
844         E111 = mK.q.mul(R111).neg();
845       case 2:
846         // Order 2.
847         E200 = mK.q.mul(R200);
848         E020 = mK.q.mul(R020);
849         E002 = mK.q.mul(R002);
850         E110 = mK.q.mul(R110);
851         E101 = mK.q.mul(R101);
852         E011 = mK.q.mul(R011);
853       case 1:
854         // Order 1 (need a minus sign).
855         E100 = mK.q.mul(R100).neg();
856         E010 = mK.q.mul(R010).neg();
857         E001 = mK.q.mul(R001).neg();
858       case 0:
859         E000 = mK.q.mul(R000);
860     }
861   }
862 
863   /**
864    * {@inheritDoc}
865    */
866   @SuppressWarnings("fallthrough")
867   @Override
868   protected void dipoleKPotentialAtI(DoubleVector uxk, DoubleVector uyk, DoubleVector uzk, int order) {
869     switch (order) {
870       default:
871       case 3:
872         // Order 3
873         DoubleVector term300 = uxk.mul(R400);
874         term300 = uyk.fma(R310, term300);
875         term300 = uzk.fma(R301, term300);
876         E300 = term300.neg();
877         DoubleVector term030 = uxk.mul(R130);
878         term030 = uyk.fma(R040, term030);
879         term030 = uzk.fma(R031, term030);
880         E030 = term030.neg();
881         DoubleVector term003 = uxk.mul(R103);
882         term003 = uyk.fma(R013, term003);
883         term003 = uzk.fma(R004, term003);
884         E003 = term003.neg();
885         DoubleVector term210 = uxk.mul(R310);
886         term210 = uyk.fma(R220, term210);
887         term210 = uzk.fma(R211, term210);
888         E210 = term210.neg();
889         DoubleVector term201 = uxk.mul(R301);
890         term201 = uyk.fma(R211, term201);
891         term201 = uzk.fma(R202, term201);
892         E201 = term201.neg();
893         DoubleVector term120 = uxk.mul(R220);
894         term120 = uyk.fma(R130, term120);
895         term120 = uzk.fma(R121, term120);
896         E120 = term120.neg();
897         DoubleVector term021 = uxk.mul(R121);
898         term021 = uyk.fma(R031, term021);
899         term021 = uzk.fma(R022, term021);
900         E021 = term021.neg();
901         DoubleVector term102 = uxk.mul(R202);
902         term102 = uyk.fma(R112, term102);
903         term102 = uzk.fma(R103, term102);
904         E102 = term102.neg();
905         DoubleVector term012 = uxk.mul(R112);
906         term012 = uyk.fma(R022, term012);
907         term012 = uzk.fma(R013, term012);
908         E012 = term012.neg();
909         DoubleVector term111 = uxk.mul(R211);
910         term111 = uyk.fma(R121, term111);
911         term111 = uzk.fma(R112, term111);
912         E111 = term111.neg();
913         // Foll through to 2nd order.
914       case 2:
915         // Order 2.
916         DoubleVector term200 = uxk.mul(R300);
917         term200 = uyk.fma(R210, term200);
918         term200 = uzk.fma(R201, term200);
919         E200 = term200;
920         DoubleVector term020 = uxk.mul(R120);
921         term020 = uyk.fma(R030, term020);
922         term020 = uzk.fma(R021, term020);
923         E020 = term020;
924         DoubleVector term002 = uxk.mul(R102);
925         term002 = uyk.fma(R012, term002);
926         term002 = uzk.fma(R003, term002);
927         E002 = term002;
928         DoubleVector term110 = uxk.mul(R210);
929         term110 = uyk.fma(R120, term110);
930         term110 = uzk.fma(R111, term110);
931         E110 = term110;
932         DoubleVector term101 = uxk.mul(R201);
933         term101 = uyk.fma(R111, term101);
934         term101 = uzk.fma(R102, term101);
935         E101 = term101;
936         DoubleVector term011 = uxk.mul(R111);
937         term011 = uyk.fma(R021, term011);
938         term011 = uzk.fma(R012, term011);
939         E011 = term011;
940       case 1:
941         // Order 1 (need a minus sign).
942         DoubleVector term100 = uxk.mul(R200);
943         term100 = uyk.fma(R110, term100);
944         term100 = uzk.fma(R101, term100);
945         E100 = term100.neg();
946         DoubleVector term010 = uxk.mul(R110);
947         term010 = uyk.fma(R020, term010);
948         term010 = uzk.fma(R011, term010);
949         E010 = term010.neg();
950         DoubleVector term001 = uxk.mul(R101);
951         term001 = uyk.fma(R011, term001);
952         term001 = uzk.fma(R002, term001);
953         E001 = term001.neg();
954       case 0:
955         DoubleVector term000 = uxk.mul(R100);
956         term000 = uyk.fma(R010, term000);
957         term000 = uzk.fma(R001, term000);
958         E000 = term000;
959     }
960   }
961 
962   /**
963    * {@inheritDoc}
964    */
965   @SuppressWarnings("fallthrough")
966   @Override
967   protected void quadrupoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order) {
968     switch (order) {
969       default:
970       case 3:
971         // Order 3 (need a minus sign).
972         DoubleVector term300 = mK.qxx.mul(R500);
973         term300 = mK.qyy.fma(R320, term300);
974         term300 = mK.qzz.fma(R302, term300);
975         term300 = mK.qxy.fma(R410, term300);
976         term300 = mK.qxz.fma(R401, term300);
977         term300 = mK.qyz.fma(R311, term300);
978         E300 = term300.neg();
979         DoubleVector term030 = mK.qxx.mul(R230);
980         term030 = mK.qyy.fma(R050, term030);
981         term030 = mK.qzz.fma(R032, term030);
982         term030 = mK.qxy.fma(R140, term030);
983         term030 = mK.qxz.fma(R131, term030);
984         term030 = mK.qyz.fma(R041, term030);
985         E030 = term030.neg();
986         DoubleVector term003 = mK.qxx.mul(R203);
987         term003 = mK.qyy.fma(R023, term003);
988         term003 = mK.qzz.fma(R005, term003);
989         term003 = mK.qxy.fma(R113, term003);
990         term003 = mK.qxz.fma(R104, term003);
991         term003 = mK.qyz.fma(R014, term003);
992         E003 = term003.neg();
993         DoubleVector term210 = mK.qxx.mul(R410);
994         term210 = mK.qyy.fma(R230, term210);
995         term210 = mK.qzz.fma(R212, term210);
996         term210 = mK.qxy.fma(R320, term210);
997         term210 = mK.qxz.fma(R311, term210);
998         term210 = mK.qyz.fma(R221, term210);
999         E210 = term210.neg();
1000         DoubleVector term201 = mK.qxx.mul(R401);
1001         term201 = mK.qyy.fma(R221, term201);
1002         term201 = mK.qzz.fma(R203, term201);
1003         term201 = mK.qxy.fma(R311, term201);
1004         term201 = mK.qxz.fma(R302, term201);
1005         term201 = mK.qyz.fma(R212, term201);
1006         E201 = term201.neg();
1007         DoubleVector term120 = mK.qxx.mul(R320);
1008         term120 = mK.qyy.fma(R140, term120);
1009         term120 = mK.qzz.fma(R122, term120);
1010         term120 = mK.qxy.fma(R230, term120);
1011         term120 = mK.qxz.fma(R221, term120);
1012         term120 = mK.qyz.fma(R131, term120);
1013         E120 = term120.neg();
1014         DoubleVector term021 = mK.qxx.mul(R221);
1015         term021 = mK.qyy.fma(R041, term021);
1016         term021 = mK.qzz.fma(R023, term021);
1017         term021 = mK.qxy.fma(R131, term021);
1018         term021 = mK.qxz.fma(R122, term021);
1019         term021 = mK.qyz.fma(R032, term021);
1020         E021 = term021.neg();
1021         DoubleVector term102 = mK.qxx.mul(R302);
1022         term102 = mK.qyy.fma(R122, term102);
1023         term102 = mK.qzz.fma(R104, term102);
1024         term102 = mK.qxy.fma(R212, term102);
1025         term102 = mK.qxz.fma(R203, term102);
1026         term102 = mK.qyz.fma(R113, term102);
1027         E102 = term102.neg();
1028         DoubleVector term012 = mK.qxx.mul(R212);
1029         term012 = mK.qyy.fma(R032, term012);
1030         term012 = mK.qzz.fma(R014, term012);
1031         term012 = mK.qxy.fma(R122, term012);
1032         term012 = mK.qxz.fma(R113, term012);
1033         term012 = mK.qyz.fma(R023, term012);
1034         E012 = term012.neg();
1035         DoubleVector term111 = mK.qxx.mul(R311);
1036         term111 = mK.qyy.fma(R131, term111);
1037         term111 = mK.qzz.fma(R113, term111);
1038         term111 = mK.qxy.fma(R221, term111);
1039         term111 = mK.qxz.fma(R212, term111);
1040         term111 = mK.qyz.fma(R122, term111);
1041         E111 = term111.neg();
1042       case 2:
1043         // Order 2.
1044         DoubleVector term200 = mK.qxx.mul(R400);
1045         term200 = mK.qyy.fma(R220, term200);
1046         term200 = mK.qzz.fma(R202, term200);
1047         term200 = mK.qxy.fma(R310, term200);
1048         term200 = mK.qxz.fma(R301, term200);
1049         term200 = mK.qyz.fma(R211, term200);
1050         E200 = term200;
1051         DoubleVector term020 = mK.qxx.mul(R220);
1052         term020 = mK.qyy.fma(R040, term020);
1053         term020 = mK.qzz.fma(R022, term020);
1054         term020 = mK.qxy.fma(R130, term020);
1055         term020 = mK.qxz.fma(R121, term020);
1056         term020 = mK.qyz.fma(R031, term020);
1057         E020 = term020;
1058         DoubleVector term002 = mK.qxx.mul(R202);
1059         term002 = mK.qyy.fma(R022, term002);
1060         term002 = mK.qzz.fma(R004, term002);
1061         term002 = mK.qxy.fma(R112, term002);
1062         term002 = mK.qxz.fma(R103, term002);
1063         term002 = mK.qyz.fma(R013, term002);
1064         E002 = term002;
1065         DoubleVector term110 = mK.qxx.mul(R310);
1066         term110 = mK.qyy.fma(R130, term110);
1067         term110 = mK.qzz.fma(R112, term110);
1068         term110 = mK.qxy.fma(R220, term110);
1069         term110 = mK.qxz.fma(R211, term110);
1070         term110 = mK.qyz.fma(R121, term110);
1071         E110 = term110;
1072         DoubleVector term101 = mK.qxx.mul(R301);
1073         term101 = mK.qyy.fma(R121, term101);
1074         term101 = mK.qzz.fma(R103, term101);
1075         term101 = mK.qxy.fma(R211, term101);
1076         term101 = mK.qxz.fma(R202, term101);
1077         term101 = mK.qyz.fma(R112, term101);
1078         E101 = term101;
1079         DoubleVector term011 = mK.qxx.mul(R211);
1080         term011 = mK.qyy.fma(R031, term011);
1081         term011 = mK.qzz.fma(R013, term011);
1082         term011 = mK.qxy.fma(R121, term011);
1083         term011 = mK.qxz.fma(R112, term011);
1084         term011 = mK.qyz.fma(R022, term011);
1085         E011 = term011;
1086       case 1:
1087         // Order 1 (need a minus sign).
1088         DoubleVector term100 = mK.qxx.mul(R300);
1089         term100 = mK.qyy.fma(R120, term100);
1090         term100 = mK.qzz.fma(R102, term100);
1091         term100 = mK.qxy.fma(R210, term100);
1092         term100 = mK.qxz.fma(R201, term100);
1093         term100 = mK.qyz.fma(R111, term100);
1094         E100 = term100.neg();
1095         DoubleVector term010 = mK.qxx.mul(R210);
1096         term010 = mK.qyy.fma(R030, term010);
1097         term010 = mK.qzz.fma(R012, term010);
1098         term010 = mK.qxy.fma(R120, term010);
1099         term010 = mK.qxz.fma(R111, term010);
1100         term010 = mK.qyz.fma(R021, term010);
1101         E010 = term010.neg();
1102         DoubleVector term001 = mK.qxx.mul(R201);
1103         term001 = mK.qyy.fma(R021, term001);
1104         term001 = mK.qzz.fma(R003, term001);
1105         term001 = mK.qxy.fma(R111, term001);
1106         term001 = mK.qxz.fma(R102, term001);
1107         term001 = mK.qyz.fma(R012, term001);
1108         E001 = term001.neg();
1109       case 0:
1110         DoubleVector term000 = mK.qxx.mul(R200);
1111         term000 = mK.qyy.fma(R020, term000);
1112         term000 = mK.qzz.fma(R002, term000);
1113         term000 = mK.qxy.fma(R110, term000);
1114         term000 = mK.qxz.fma(R101, term000);
1115         term000 = mK.qyz.fma(R011, term000);
1116         E000 = term000;
1117     }
1118   }
1119 
1120   /**
1121    * {@inheritDoc}
1122    */
1123   @SuppressWarnings("fallthrough")
1124   @Override
1125   protected void multipoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order) {
1126     switch (order) {
1127       default:
1128       case 3:
1129         DoubleVector term300 = mI.q.mul(R300);
1130         term300 = mI.dx.fma(R400.neg(), term300);
1131         term300 = mI.dy.fma(R310.neg(), term300);
1132         term300 = mI.dz.fma(R301.neg(), term300);
1133         term300 = mI.qxx.fma(R500, term300);
1134         term300 = mI.qyy.fma(R320, term300);
1135         term300 = mI.qzz.fma(R302, term300);
1136         term300 = mI.qxy.fma(R410, term300);
1137         term300 = mI.qxz.fma(R401, term300);
1138         term300 = mI.qyz.fma(R311, term300);
1139         E300 = term300;
1140         DoubleVector term030 = mI.q.mul(R030);
1141         term030 = mI.dx.fma(R130.neg(), term030);
1142         term030 = mI.dy.fma(R040.neg(), term030);
1143         term030 = mI.dz.fma(R031.neg(), term030);
1144         term030 = mI.qxx.fma(R230, term030);
1145         term030 = mI.qyy.fma(R050, term030);
1146         term030 = mI.qzz.fma(R032, term030);
1147         term030 = mI.qxy.fma(R140, term030);
1148         term030 = mI.qxz.fma(R131, term030);
1149         term030 = mI.qyz.fma(R041, term030);
1150         E030 = term030;
1151         DoubleVector term003 = mI.q.mul(R003);
1152         term003 = mI.dx.fma(R103.neg(), term003);
1153         term003 = mI.dy.fma(R013.neg(), term003);
1154         term003 = mI.dz.fma(R004.neg(), term003);
1155         term003 = mI.qxx.fma(R203, term003);
1156         term003 = mI.qyy.fma(R023, term003);
1157         term003 = mI.qzz.fma(R005, term003);
1158         term003 = mI.qxy.fma(R113, term003);
1159         term003 = mI.qxz.fma(R104, term003);
1160         term003 = mI.qyz.fma(R014, term003);
1161         E003 = term003;
1162         DoubleVector term210 = mI.q.mul(R210);
1163         term210 = mI.dx.fma(R310.neg(), term210);
1164         term210 = mI.dy.fma(R220.neg(), term210);
1165         term210 = mI.dz.fma(R211.neg(), term210);
1166         term210 = mI.qxx.fma(R410, term210);
1167         term210 = mI.qyy.fma(R230, term210);
1168         term210 = mI.qzz.fma(R212, term210);
1169         term210 = mI.qxy.fma(R320, term210);
1170         term210 = mI.qxz.fma(R311, term210);
1171         term210 = mI.qyz.fma(R221, term210);
1172         E210 = term210;
1173         DoubleVector term201 = mI.q.mul(R201);
1174         term201 = mI.dx.fma(R301.neg(), term201);
1175         term201 = mI.dy.fma(R211.neg(), term201);
1176         term201 = mI.dz.fma(R202.neg(), term201);
1177         term201 = mI.qxx.fma(R401, term201);
1178         term201 = mI.qyy.fma(R221, term201);
1179         term201 = mI.qzz.fma(R203, term201);
1180         term201 = mI.qxy.fma(R311, term201);
1181         term201 = mI.qxz.fma(R302, term201);
1182         term201 = mI.qyz.fma(R212, term201);
1183         E201 = term201;
1184         DoubleVector term120 = mI.q.mul(R120);
1185         term120 = mI.dx.fma(R220.neg(), term120);
1186         term120 = mI.dy.fma(R130.neg(), term120);
1187         term120 = mI.dz.fma(R121.neg(), term120);
1188         term120 = mI.qxx.fma(R320, term120);
1189         term120 = mI.qyy.fma(R140, term120);
1190         term120 = mI.qzz.fma(R122, term120);
1191         term120 = mI.qxy.fma(R230, term120);
1192         term120 = mI.qxz.fma(R221, term120);
1193         term120 = mI.qyz.fma(R131, term120);
1194         E120 = term120;
1195         DoubleVector term021 = mI.q.mul(R021);
1196         term021 = mI.dx.fma(R121.neg(), term021);
1197         term021 = mI.dy.fma(R031.neg(), term021);
1198         term021 = mI.dz.fma(R022.neg(), term021);
1199         term021 = mI.qxx.fma(R221, term021);
1200         term021 = mI.qyy.fma(R041, term021);
1201         term021 = mI.qzz.fma(R023, term021);
1202         term021 = mI.qxy.fma(R131, term021);
1203         term021 = mI.qxz.fma(R122, term021);
1204         term021 = mI.qyz.fma(R032, term021);
1205         E021 = term021;
1206         DoubleVector term102 = mI.q.mul(R102);
1207         term102 = mI.dx.fma(R202.neg(), term102);
1208         term102 = mI.dy.fma(R112.neg(), term102);
1209         term102 = mI.dz.fma(R103.neg(), term102);
1210         term102 = mI.qxx.fma(R302, term102);
1211         term102 = mI.qyy.fma(R122, term102);
1212         term102 = mI.qzz.fma(R104, term102);
1213         term102 = mI.qxy.fma(R212, term102);
1214         term102 = mI.qxz.fma(R203, term102);
1215         term102 = mI.qyz.fma(R113, term102);
1216         E102 = term102;
1217         DoubleVector term012 = mI.q.mul(R012);
1218         term012 = mI.dx.fma(R112.neg(), term012);
1219         term012 = mI.dy.fma(R022.neg(), term012);
1220         term012 = mI.dz.fma(R013.neg(), term012);
1221         term012 = mI.qxx.fma(R212, term012);
1222         term012 = mI.qyy.fma(R032, term012);
1223         term012 = mI.qzz.fma(R014, term012);
1224         term012 = mI.qxy.fma(R122, term012);
1225         term012 = mI.qxz.fma(R113, term012);
1226         term012 = mI.qyz.fma(R023, term012);
1227         E012 = term012;
1228         DoubleVector term111 = mI.q.mul(R111);
1229         term111 = mI.dx.fma(R211.neg(), term111);
1230         term111 = mI.dy.fma(R121.neg(), term111);
1231         term111 = mI.dz.fma(R112.neg(), term111);
1232         term111 = mI.qxx.fma(R311, term111);
1233         term111 = mI.qyy.fma(R131, term111);
1234         term111 = mI.qzz.fma(R113, term111);
1235         term111 = mI.qxy.fma(R221, term111);
1236         term111 = mI.qxz.fma(R212, term111);
1237         term111 = mI.qyz.fma(R122, term111);
1238         E111 = term111;
1239       case 2:
1240         DoubleVector term200 = mI.q.mul(R200);
1241         term200 = mI.dx.fma(R300.neg(), term200);
1242         term200 = mI.dy.fma(R210.neg(), term200);
1243         term200 = mI.dz.fma(R201.neg(), term200);
1244         term200 = mI.qxx.fma(R400, term200);
1245         term200 = mI.qyy.fma(R220, term200);
1246         term200 = mI.qzz.fma(R202, term200);
1247         term200 = mI.qxy.fma(R310, term200);
1248         term200 = mI.qxz.fma(R301, term200);
1249         term200 = mI.qyz.fma(R211, term200);
1250         E200 = term200;
1251         DoubleVector term020 = mI.q.mul(R020);
1252         term020 = mI.dx.fma(R120.neg(), term020);
1253         term020 = mI.dy.fma(R030.neg(), term020);
1254         term020 = mI.dz.fma(R021.neg(), term020);
1255         term020 = mI.qxx.fma(R220, term020);
1256         term020 = mI.qyy.fma(R040, term020);
1257         term020 = mI.qzz.fma(R022, term020);
1258         term020 = mI.qxy.fma(R130, term020);
1259         term020 = mI.qxz.fma(R121, term020);
1260         term020 = mI.qyz.fma(R031, term020);
1261         E020 = term020;
1262         DoubleVector term002 = mI.q.mul(R002);
1263         term002 = mI.dx.fma(R102.neg(), term002);
1264         term002 = mI.dy.fma(R012.neg(), term002);
1265         term002 = mI.dz.fma(R003.neg(), term002);
1266         term002 = mI.qxx.fma(R202, term002);
1267         term002 = mI.qyy.fma(R022, term002);
1268         term002 = mI.qzz.fma(R004, term002);
1269         term002 = mI.qxy.fma(R112, term002);
1270         term002 = mI.qxz.fma(R103, term002);
1271         term002 = mI.qyz.fma(R013, term002);
1272         E002 = term002;
1273         DoubleVector term110 = mI.q.mul(R110);
1274         term110 = mI.dx.fma(R210.neg(), term110);
1275         term110 = mI.dy.fma(R120.neg(), term110);
1276         term110 = mI.dz.fma(R111.neg(), term110);
1277         term110 = mI.qxx.fma(R310, term110);
1278         term110 = mI.qyy.fma(R130, term110);
1279         term110 = mI.qzz.fma(R112, term110);
1280         term110 = mI.qxy.fma(R220, term110);
1281         term110 = mI.qxz.fma(R211, term110);
1282         term110 = mI.qyz.fma(R121, term110);
1283         E110 = term110;
1284         DoubleVector term101 = mI.q.mul(R101);
1285         term101 = mI.dx.fma(R201.neg(), term101);
1286         term101 = mI.dy.fma(R111.neg(), term101);
1287         term101 = mI.dz.fma(R102.neg(), term101);
1288         term101 = mI.qxx.fma(R301, term101);
1289         term101 = mI.qyy.fma(R121, term101);
1290         term101 = mI.qzz.fma(R103, term101);
1291         term101 = mI.qxy.fma(R211, term101);
1292         term101 = mI.qxz.fma(R202, term101);
1293         term101 = mI.qyz.fma(R112, term101);
1294         E101 = term101;
1295         DoubleVector term011 = mI.q.mul(R011);
1296         term011 = mI.dx.fma(R111.neg(), term011);
1297         term011 = mI.dy.fma(R021.neg(), term011);
1298         term011 = mI.dz.fma(R012.neg(), term011);
1299         term011 = mI.qxx.fma(R211, term011);
1300         term011 = mI.qyy.fma(R031, term011);
1301         term011 = mI.qzz.fma(R013, term011);
1302         term011 = mI.qxy.fma(R121, term011);
1303         term011 = mI.qxz.fma(R112, term011);
1304         term011 = mI.qyz.fma(R022, term011);
1305         E011 = term011;
1306       case 1:
1307         DoubleVector term100 = mI.q.mul(R100);
1308         term100 = mI.dx.fma(R200.neg(), term100);
1309         term100 = mI.dy.fma(R110.neg(), term100);
1310         term100 = mI.dz.fma(R101.neg(), term100);
1311         term100 = mI.qxx.fma(R300, term100);
1312         term100 = mI.qyy.fma(R120, term100);
1313         term100 = mI.qzz.fma(R102, term100);
1314         term100 = mI.qxy.fma(R210, term100);
1315         term100 = mI.qxz.fma(R201, term100);
1316         term100 = mI.qyz.fma(R111, term100);
1317         E100 = term100;
1318         DoubleVector term010 = mI.q.mul(R010);
1319         term010 = mI.dx.fma(R110.neg(), term010);
1320         term010 = mI.dy.fma(R020.neg(), term010);
1321         term010 = mI.dz.fma(R011.neg(), term010);
1322         term010 = mI.qxx.fma(R210, term010);
1323         term010 = mI.qyy.fma(R030, term010);
1324         term010 = mI.qzz.fma(R012, term010);
1325         term010 = mI.qxy.fma(R120, term010);
1326         term010 = mI.qxz.fma(R111, term010);
1327         term010 = mI.qyz.fma(R021, term010);
1328         E010 = term010;
1329         DoubleVector term001 = mI.q.mul(R001);
1330         term001 = mI.dx.fma(R101.neg(), term001);
1331         term001 = mI.dy.fma(R011.neg(), term001);
1332         term001 = mI.dz.fma(R002.neg(), term001);
1333         term001 = mI.qxx.fma(R201, term001);
1334         term001 = mI.qyy.fma(R021, term001);
1335         term001 = mI.qzz.fma(R003, term001);
1336         term001 = mI.qxy.fma(R111, term001);
1337         term001 = mI.qxz.fma(R102, term001);
1338         term001 = mI.qyz.fma(R012, term001);
1339         E001 = term001;
1340       case 0:
1341         DoubleVector term000 = mI.q.mul(R000);
1342         term000 = mI.dx.fma(R100.neg(), term000);
1343         term000 = mI.dy.fma(R010.neg(), term000);
1344         term000 = mI.dz.fma(R001.neg(), term000);
1345         term000 = mI.qxx.fma(R200, term000);
1346         term000 = mI.qyy.fma(R020, term000);
1347         term000 = mI.qzz.fma(R002, term000);
1348         term000 = mI.qxy.fma(R110, term000);
1349         term000 = mI.qxz.fma(R101, term000);
1350         term000 = mI.qyz.fma(R011, term000);
1351         E000 = term000;
1352     }
1353   }
1354 
1355   /**
1356    * {@inheritDoc}
1357    */
1358   @SuppressWarnings("fallthrough")
1359   @Override
1360   protected void chargeIPotentialAtK(PolarizableMultipoleSIMD mI, int order) {
1361     switch (order) {
1362       default:
1363       case 3:
1364         E300 = mI.q.mul(R300);
1365         E030 = mI.q.mul(R030);
1366         E003 = mI.q.mul(R003);
1367         E210 = mI.q.mul(R210);
1368         E201 = mI.q.mul(R201);
1369         E120 = mI.q.mul(R120);
1370         E021 = mI.q.mul(R021);
1371         E102 = mI.q.mul(R102);
1372         E012 = mI.q.mul(R012);
1373         E111 = mI.q.mul(R111);
1374       case 2:
1375         E200 = mI.q.mul(R200);
1376         E020 = mI.q.mul(R020);
1377         E002 = mI.q.mul(R002);
1378         E110 = mI.q.mul(R110);
1379         E101 = mI.q.mul(R101);
1380         E011 = mI.q.mul(R011);
1381       case 1:
1382         E100 = mI.q.mul(R100);
1383         E010 = mI.q.mul(R010);
1384         E001 = mI.q.mul(R001);
1385       case 0:
1386         E000 = mI.q.mul(R000);
1387     }
1388   }
1389 
1390   /**
1391    * {@inheritDoc}
1392    */
1393   @SuppressWarnings("fallthrough")
1394   @Override
1395   protected void dipoleIPotentialAtK(DoubleVector uxi, DoubleVector uyi, DoubleVector uzi, int order) {
1396     switch (order) {
1397       default:
1398       case 3:
1399         // Order 3
1400         DoubleVector term300 = uxi.mul(R400.neg());
1401         term300 = uyi.fma(R310.neg(), term300);
1402         term300 = uzi.fma(R301.neg(), term300);
1403         E300 = term300;
1404         DoubleVector term030 = uxi.mul(R130.neg());
1405         term030 = uyi.fma(R040.neg(), term030);
1406         term030 = uzi.fma(R031.neg(), term030);
1407         E030 = term030;
1408         DoubleVector term003 = uxi.mul(R103.neg());
1409         term003 = uyi.fma(R013.neg(), term003);
1410         term003 = uzi.fma(R004.neg(), term003);
1411         E003 = term003;
1412         DoubleVector term210 = uxi.mul(R310.neg());
1413         term210 = uyi.fma(R220.neg(), term210);
1414         term210 = uzi.fma(R211.neg(), term210);
1415         E210 = term210;
1416         DoubleVector term201 = uxi.mul(R301.neg());
1417         term201 = uyi.fma(R211.neg(), term201);
1418         term201 = uzi.fma(R202.neg(), term201);
1419         E201 = term201;
1420         DoubleVector term120 = uxi.mul(R220.neg());
1421         term120 = uyi.fma(R130.neg(), term120);
1422         term120 = uzi.fma(R121.neg(), term120);
1423         E120 = term120;
1424         DoubleVector term021 = uxi.mul(R121.neg());
1425         term021 = uyi.fma(R031.neg(), term021);
1426         term021 = uzi.fma(R022.neg(), term021);
1427         E021 = term021;
1428         DoubleVector term102 = uxi.mul(R202.neg());
1429         term102 = uyi.fma(R112.neg(), term102);
1430         term102 = uzi.fma(R103.neg(), term102);
1431         E102 = term102;
1432         DoubleVector term012 = uxi.mul(R112.neg());
1433         term012 = uyi.fma(R022.neg(), term012);
1434         term012 = uzi.fma(R013.neg(), term012);
1435         E012 = term012;
1436         DoubleVector term111 = uxi.mul(R211.neg());
1437         term111 = uyi.fma(R121.neg(), term111);
1438         term111 = uzi.fma(R112.neg(), term111);
1439         E111 = term111;
1440         // Fall through to 2nd order.
1441       case 2:
1442         // Order 2
1443         DoubleVector term200 = uxi.mul(R300.neg());
1444         term200 = uyi.fma(R210.neg(), term200);
1445         term200 = uzi.fma(R201.neg(), term200);
1446         E200 = term200;
1447         DoubleVector term020 = uxi.mul(R120.neg());
1448         term020 = uyi.fma(R030.neg(), term020);
1449         term020 = uzi.fma(R021.neg(), term020);
1450         E020 = term020;
1451         DoubleVector term002 = uxi.mul(R102.neg());
1452         term002 = uyi.fma(R012.neg(), term002);
1453         term002 = uzi.fma(R003.neg(), term002);
1454         E002 = term002;
1455         DoubleVector term110 = uxi.mul(R210.neg());
1456         term110 = uyi.fma(R120.neg(), term110);
1457         term110 = uzi.fma(R111.neg(), term110);
1458         E110 = term110;
1459         DoubleVector term101 = uxi.mul(R201.neg());
1460         term101 = uyi.fma(R111.neg(), term101);
1461         term101 = uzi.fma(R102.neg(), term101);
1462         E101 = term101;
1463         DoubleVector term011 = uxi.mul(R111.neg());
1464         term011 = uyi.fma(R021.neg(), term011);
1465         term011 = uzi.fma(R012.neg(), term011);
1466         E011 = term011;
1467         // Fall through to 1st order.
1468       case 1:
1469         // Order 1
1470         // This is d/dX of equation 3.1.3 in the Stone book.
1471         DoubleVector term100 = uxi.mul(R200.neg());
1472         term100 = uyi.fma(R110.neg(), term100);
1473         term100 = uzi.fma(R101.neg(), term100);
1474         E100 = term100;
1475         // This is d/dY of equation 3.1.3 in the Stone book.
1476         DoubleVector term010 = uxi.mul(R110.neg());
1477         term010 = uyi.fma(R020.neg(), term010);
1478         term010 = uzi.fma(R011.neg(), term010);
1479         E010 = term010;
1480         DoubleVector term001 = uxi.mul(R101.neg());
1481         term001 = uyi.fma(R011.neg(), term001);
1482         term001 = uzi.fma(R002.neg(), term001);
1483         E001 = term001;
1484         // Fall through to the potential.
1485       case 0:
1486         // This is equation 3.1.3 in the Stone book.
1487         DoubleVector term000 = uxi.mul(R100.neg());
1488         term000 = uyi.fma(R010.neg(), term000);
1489         term000 = uzi.fma(R001.neg(), term000);
1490         E000 = term000;
1491     }
1492   }
1493 
1494   /**
1495    * {@inheritDoc}
1496    */
1497   @SuppressWarnings("fallthrough")
1498   @Override
1499   protected void quadrupoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order) {
1500     switch (order) {
1501       default:
1502       case 3:
1503         DoubleVector term300 = mI.qxx.mul(R500);
1504         term300 = mI.qyy.fma(R320, term300);
1505         term300 = mI.qzz.fma(R302, term300);
1506         term300 = mI.qxy.fma(R410, term300);
1507         term300 = mI.qxz.fma(R401, term300);
1508         term300 = mI.qyz.fma(R311, term300);
1509         E300 = term300;
1510         DoubleVector term030 = mI.qxx.mul(R230);
1511         term030 = mI.qyy.fma(R050, term030);
1512         term030 = mI.qzz.fma(R032, term030);
1513         term030 = mI.qxy.fma(R140, term030);
1514         term030 = mI.qxz.fma(R131, term030);
1515         term030 = mI.qyz.fma(R041, term030);
1516         E030 = term030;
1517         DoubleVector term003 = mI.qxx.mul(R203);
1518         term003 = mI.qyy.fma(R023, term003);
1519         term003 = mI.qzz.fma(R005, term003);
1520         term003 = mI.qxy.fma(R113, term003);
1521         term003 = mI.qxz.fma(R104, term003);
1522         term003 = mI.qyz.fma(R014, term003);
1523         E003 = term003;
1524         DoubleVector term210 = mI.qxx.mul(R410);
1525         term210 = mI.qyy.fma(R230, term210);
1526         term210 = mI.qzz.fma(R212, term210);
1527         term210 = mI.qxy.fma(R320, term210);
1528         term210 = mI.qxz.fma(R311, term210);
1529         term210 = mI.qyz.fma(R221, term210);
1530         E210 = term210;
1531         DoubleVector term201 = mI.qxx.mul(R401);
1532         term201 = mI.qyy.fma(R221, term201);
1533         term201 = mI.qzz.fma(R203, term201);
1534         term201 = mI.qxy.fma(R311, term201);
1535         term201 = mI.qxz.fma(R302, term201);
1536         term201 = mI.qyz.fma(R212, term201);
1537         E201 = term201;
1538         DoubleVector term120 = mI.qxx.mul(R320);
1539         term120 = mI.qyy.fma(R140, term120);
1540         term120 = mI.qzz.fma(R122, term120);
1541         term120 = mI.qxy.fma(R230, term120);
1542         term120 = mI.qxz.fma(R221, term120);
1543         term120 = mI.qyz.fma(R131, term120);
1544         E120 = term120;
1545         DoubleVector term021 = mI.qxx.mul(R221);
1546         term021 = mI.qyy.fma(R041, term021);
1547         term021 = mI.qzz.fma(R023, term021);
1548         term021 = mI.qxy.fma(R131, term021);
1549         term021 = mI.qxz.fma(R122, term021);
1550         term021 = mI.qyz.fma(R032, term021);
1551         E021 = term021;
1552         DoubleVector term102 = mI.qxx.mul(R302);
1553         term102 = mI.qyy.fma(R122, term102);
1554         term102 = mI.qzz.fma(R104, term102);
1555         term102 = mI.qxy.fma(R212, term102);
1556         term102 = mI.qxz.fma(R203, term102);
1557         term102 = mI.qyz.fma(R113, term102);
1558         E102 = term102;
1559         DoubleVector term012 = mI.qxx.mul(R212);
1560         term012 = mI.qyy.fma(R032, term012);
1561         term012 = mI.qzz.fma(R014, term012);
1562         term012 = mI.qxy.fma(R122, term012);
1563         term012 = mI.qxz.fma(R113, term012);
1564         term012 = mI.qyz.fma(R023, term012);
1565         E012 = term012;
1566         DoubleVector term111 = mI.qxx.mul(R311);
1567         term111 = mI.qyy.fma(R131, term111);
1568         term111 = mI.qzz.fma(R113, term111);
1569         term111 = mI.qxy.fma(R221, term111);
1570         term111 = mI.qxz.fma(R212, term111);
1571         term111 = mI.qyz.fma(R122, term111);
1572         E111 = term111;
1573       case 2:
1574         DoubleVector term200 = mI.qxx.mul(R400);
1575         term200 = mI.qyy.fma(R220, term200);
1576         term200 = mI.qzz.fma(R202, term200);
1577         term200 = mI.qxy.fma(R310, term200);
1578         term200 = mI.qxz.fma(R301, term200);
1579         term200 = mI.qyz.fma(R211, term200);
1580         E200 = term200;
1581         DoubleVector term020 = mI.qxx.mul(R220);
1582         term020 = mI.qyy.fma(R040, term020);
1583         term020 = mI.qzz.fma(R022, term020);
1584         term020 = mI.qxy.fma(R130, term020);
1585         term020 = mI.qxz.fma(R121, term020);
1586         term020 = mI.qyz.fma(R031, term020);
1587         E020 = term020;
1588         DoubleVector term002 = mI.qxx.mul(R202);
1589         term002 = mI.qyy.fma(R022, term002);
1590         term002 = mI.qzz.fma(R004, term002);
1591         term002 = mI.qxy.fma(R112, term002);
1592         term002 = mI.qxz.fma(R103, term002);
1593         term002 = mI.qyz.fma(R013, term002);
1594         E002 = term002;
1595         DoubleVector term110 = mI.qxx.mul(R310);
1596         term110 = mI.qyy.fma(R130, term110);
1597         term110 = mI.qzz.fma(R112, term110);
1598         term110 = mI.qxy.fma(R220, term110);
1599         term110 = mI.qxz.fma(R211, term110);
1600         term110 = mI.qyz.fma(R121, term110);
1601         E110 = term110;
1602         DoubleVector term101 = mI.qxx.mul(R301);
1603         term101 = mI.qyy.fma(R121, term101);
1604         term101 = mI.qzz.fma(R103, term101);
1605         term101 = mI.qxy.fma(R211, term101);
1606         term101 = mI.qxz.fma(R202, term101);
1607         term101 = mI.qyz.fma(R112, term101);
1608         E101 = term101;
1609         DoubleVector term011 = mI.qxx.mul(R211);
1610         term011 = mI.qyy.fma(R031, term011);
1611         term011 = mI.qzz.fma(R013, term011);
1612         term011 = mI.qxy.fma(R121, term011);
1613         term011 = mI.qxz.fma(R112, term011);
1614         term011 = mI.qyz.fma(R022, term011);
1615         E011 = term011;
1616       case 1:
1617         DoubleVector term100 = mI.qxx.mul(R300);
1618         term100 = mI.qyy.fma(R120, term100);
1619         term100 = mI.qzz.fma(R102, term100);
1620         term100 = mI.qxy.fma(R210, term100);
1621         term100 = mI.qxz.fma(R201, term100);
1622         term100 = mI.qyz.fma(R111, term100);
1623         E100 = term100;
1624         DoubleVector term010 = mI.qxx.mul(R210);
1625         term010 = mI.qyy.fma(R030, term010);
1626         term010 = mI.qzz.fma(R012, term010);
1627         term010 = mI.qxy.fma(R120, term010);
1628         term010 = mI.qxz.fma(R111, term010);
1629         term010 = mI.qyz.fma(R021, term010);
1630         E010 = term010;
1631         DoubleVector term001 = mI.qxx.mul(R201);
1632         term001 = mI.qyy.fma(R021, term001);
1633         term001 = mI.qzz.fma(R003, term001);
1634         term001 = mI.qxy.fma(R111, term001);
1635         term001 = mI.qxz.fma(R102, term001);
1636         term001 = mI.qyz.fma(R012, term001);
1637         E001 = term001;
1638       case 0:
1639         DoubleVector term000 = mI.qxx.mul(R200);
1640         term000 = mI.qyy.fma(R020, term000);
1641         term000 = mI.qzz.fma(R002, term000);
1642         term000 = mI.qxy.fma(R110, term000);
1643         term000 = mI.qxz.fma(R101, term000);
1644         term000 = mI.qyz.fma(R011, term000);
1645         E000 = term000;
1646     }
1647   }
1648 }