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  /**
43   * The CoulombTensorQISIMD class computes derivatives of 1/|<b>r</b>| via recursion to arbitrary order
44   * for Cartesian multipoles in a quasi-internal frame using SIMD instructions.
45   *
46   * @author Michael J. Schnieders
47   * @see <a href="http://doi.org/10.1142/9789812830364_0002" target="_blank"> Matt Challacombe, Eric
48   * Schwegler and Jan Almlof, Modern developments in Hartree-Fock theory: Fast methods for
49   * computing the Coulomb matrix. Computational Chemistry: Review of Current Trends. pp. 53-107,
50   * Ed. J. Leczszynski, World Scientifc, 1996. </a>
51   * @since 1.0
52   */
53  public class CoulombTensorQISIMD extends MultipoleTensorSIMD {
54  
55    private static final DoubleVector zero = DoubleVector.zero(DoubleVector.SPECIES_PREFERRED);
56  
57    /**
58     * Create a new CoulombTensorQI object.
59     *
60     * @param order The tensor order.
61     */
62    public CoulombTensorQISIMD(int order) {
63      super(CoordinateSystem.QI, order);
64      operator = Operator.COULOMB;
65    }
66  
67    /**
68     * {@inheritDoc}
69     */
70    @Override
71    public void setR(DoubleVector dx, DoubleVector dy, DoubleVector dz) {
72      // QI separation along x is zero.
73      x = DoubleVector.broadcast(dz.species(), 0.0);
74      // QI separation along y is zero.
75      y = DoubleVector.broadcast(dz.species(), 0.0);
76      // Compute the separation distance squared.
77      DoubleVector dx2 = dx.mul(dx);
78      DoubleVector dy2 = dy.mul(dy);
79      DoubleVector dz2 = dz.mul(dz);
80      r2 = dx2.add(dy2).add(dz2);
81      // QI separation vector is along z.
82      z = r2.sqrt();
83      R = z;
84    }
85  
86    /**
87     * Generate source terms for the Coulomb Challacombe et al. recursion.
88     *
89     * @param T000 Location to store the source terms.
90     */
91    protected void source(DoubleVector[] T000) {
92      // Challacombe et al. Equation 21, last factor.
93      // == (1/r) * (1/r^3) * (1/r^5) * (1/r^7) * ...
94      DoubleVector ONE = DoubleVector.broadcast(R.species(), 1.0);
95      DoubleVector ir = ONE.div(R);
96      DoubleVector ir2 = ir.mul(ir);
97      for (int n = 0; n < o1; n++) {
98        T000[n] = ir.mul(coulombSource[n]);
99        ir = ir.mul(ir2);
100     }
101   }
102 
103   /**
104    * {@inheritDoc}
105    */
106   @Override
107   protected void order1() {
108     source(work);
109     DoubleVector term0000 = work[0];
110     DoubleVector term0001 = work[1];
111     R000 = term0000;
112     R001 = z.mul(term0001);
113   }
114 
115   /**
116    * {@inheritDoc}
117    */
118   @Override
119   protected void order2() {
120     source(work);
121     DoubleVector term0000 = work[0];
122     DoubleVector term0001 = work[1];
123     DoubleVector term0002 = work[2];
124     R000 = term0000;
125     R200 = term0001;
126     R020 = term0001;
127     R001 = z.mul(term0001);
128     DoubleVector term0011 = z.mul(term0002);
129     R002 = z.fma(term0011, term0001);
130   }
131 
132   /**
133    * {@inheritDoc}
134    */
135   @Override
136   protected void order3() {
137     source(work);
138     DoubleVector term0000 = work[0];
139     DoubleVector term0001 = work[1];
140     DoubleVector term0002 = work[2];
141     DoubleVector term0003 = work[3];
142     R000 = term0000;
143     R200 = term0001;
144     R020 = term0001;
145     R001 = z.mul(term0001);
146     DoubleVector term0011 = z.mul(term0002);
147     R002 = z.fma(term0011, term0001);
148     DoubleVector term0012 = z.mul(term0003);
149     DoubleVector term0021 = z.fma(term0012, term0002);
150     R003 = z.fma(term0021, term0011.mul(2));
151     R021 = z.mul(term0002);
152     R201 = z.mul(term0002);
153   }
154 
155   /**
156    * {@inheritDoc}
157    */
158   @Override
159   protected void order4() {
160     source(work);
161     DoubleVector term0000 = work[0];
162     DoubleVector term0001 = work[1];
163     DoubleVector term0002 = work[2];
164     DoubleVector term0003 = work[3];
165     DoubleVector term0004 = work[4];
166     R000 = term0000;
167     R200 = term0001;
168     R400 = term0002.mul(3);
169     R020 = term0001;
170     R040 = term0002.mul(3);
171     R220 = term0002;
172     R001 = z.mul(term0001);
173     DoubleVector term0011 = z.mul(term0002);
174     R002 = z.fma(term0011, term0001);
175     DoubleVector term0012 = z.mul(term0003);
176     DoubleVector term0021 = z.fma(term0012, term0002);
177     R003 = z.fma(term0021, term0011.mul(2));
178     DoubleVector term0013 = z.mul(term0004);
179     DoubleVector term0022 = z.fma(term0013, term0003);
180     DoubleVector term0031 = z.fma(term0022, term0012.mul(2));
181     R004 = z.fma(term0031, term0021.mul(3));
182     R021 = z.mul(term0002);
183     DoubleVector term0211 = z.mul(term0003);
184     R022 = z.fma(term0211, term0002);
185     R201 = z.mul(term0002);
186     DoubleVector term2011 = z.mul(term0003);
187     R202 = z.fma(term2011, term0002);
188   }
189 
190   /**
191    * {@inheritDoc}
192    */
193   @Override
194   protected void order5() {
195     source(work);
196     DoubleVector term0000 = work[0];
197     DoubleVector term0001 = work[1];
198     DoubleVector term0002 = work[2];
199     DoubleVector term0003 = work[3];
200     DoubleVector term0004 = work[4];
201     DoubleVector term0005 = work[5];
202     R000 = term0000;
203     R200 = term0001;
204     R400 = term0002.mul(3);
205     DoubleVector term4001 = term0003.mul(3);
206     R020 = term0001;
207     R040 = term0002.mul(3);
208     DoubleVector term0401 = term0003.mul(3);
209     R220 = term0002;
210     R001 = z.mul(term0001);
211     DoubleVector term0011 = z.mul(term0002);
212     R002 = z.fma(term0011, term0001);
213     DoubleVector term0012 = z.mul(term0003);
214     DoubleVector term0021 = z.fma(term0012, term0002);
215     R003 = z.fma(term0021, term0011.mul(2));
216     DoubleVector term0013 = z.mul(term0004);
217     DoubleVector term0022 = z.fma(term0013, term0003);
218     DoubleVector term0031 = z.fma(term0022, term0012.mul(2));
219     R004 = z.fma(term0031, term0021.mul(3));
220     DoubleVector term0014 = z.mul(term0005);
221     DoubleVector term0023 = z.fma(term0014, term0004);
222     DoubleVector term0032 = z.fma(term0023, term0013.mul(2));
223     DoubleVector term0041 = z.fma(term0032, term0022.mul(3));
224     R005 = z.fma(term0041, term0031.mul(4));
225     R021 = z.mul(term0002);
226     DoubleVector term0211 = z.mul(term0003);
227     R022 = z.fma(term0211, term0002);
228     DoubleVector term0212 = z.mul(term0004);
229     DoubleVector term0221 = z.fma(term0212, term0003);
230     R023 = z.fma(term0221, term0211.mul(2));
231     R041 = z.mul(term0401);
232     R201 = z.mul(term0002);
233     DoubleVector term2011 = z.mul(term0003);
234     R202 = z.fma(term2011, term0002);
235     DoubleVector term2012 = z.mul(term0004);
236     DoubleVector term2021 = z.fma(term2012, term0003);
237     R203 = z.fma(term2021, term2011.mul(2));
238     R221 = z.mul(term0003);
239     R401 = z.mul(term4001);
240   }
241 
242   /**
243    * {@inheritDoc}
244    */
245   @Override
246   protected void order6() {
247     source(work);
248     DoubleVector term0000 = work[0];
249     DoubleVector term0001 = work[1];
250     DoubleVector term0002 = work[2];
251     DoubleVector term0003 = work[3];
252     DoubleVector term0004 = work[4];
253     DoubleVector term0005 = work[5];
254     DoubleVector term0006 = work[6];
255     R000 = term0000;
256     R200 = term0001;
257     R400 = term0002.mul(3);
258     DoubleVector term4001 = term0003.mul(3);
259     DoubleVector term4002 = term0004.mul(3);
260     R600 = term4001.mul(5);
261     R020 = term0001;
262     R040 = term0002.mul(3);
263     DoubleVector term0401 = term0003.mul(3);
264     DoubleVector term0402 = term0004.mul(3);
265     R060 = term0401.mul(5);
266     R220 = term0002;
267     R240 = term0003.mul(3);
268     R420 = term4001;
269     R001 = z.mul(term0001);
270     DoubleVector term0011 = z.mul(term0002);
271     R002 = z.fma(term0011, term0001);
272     DoubleVector term0012 = z.mul(term0003);
273     DoubleVector term0021 = z.fma(term0012, term0002);
274     R003 = z.fma(term0021, term0011.mul(2));
275     DoubleVector term0013 = z.mul(term0004);
276     DoubleVector term0022 = z.fma(term0013, term0003);
277     DoubleVector term0031 = z.fma(term0022, term0012.mul(2));
278     R004 = z.fma(term0031, term0021.mul(3));
279     DoubleVector term0014 = z.mul(term0005);
280     DoubleVector term0023 = z.fma(term0014, term0004);
281     DoubleVector term0032 = z.fma(term0023, term0013.mul(2));
282     DoubleVector term0041 = z.fma(term0032, term0022.mul(3));
283     R005 = z.fma(term0041, term0031.mul(4));
284     DoubleVector term0015 = z.mul(term0006);
285     DoubleVector term0024 = z.fma(term0015, term0005);
286     DoubleVector term0033 = z.fma(term0024, term0014.mul(2));
287     DoubleVector term0042 = z.fma(term0033, term0023.mul(3));
288     DoubleVector term0051 = z.fma(term0042, term0032.mul(4));
289     R006 = z.fma(term0051, term0041.mul(5));
290     R021 = z.mul(term0002);
291     DoubleVector term0211 = z.mul(term0003);
292     R022 = z.fma(term0211, term0002);
293     DoubleVector term0212 = z.mul(term0004);
294     DoubleVector term0221 = z.fma(term0212, term0003);
295     R023 = z.fma(term0221, term0211.mul(2));
296     DoubleVector term0213 = z.mul(term0005);
297     DoubleVector term0222 = z.fma(term0213, term0004);
298     DoubleVector term0231 = z.fma(term0222, term0212.mul(2));
299     R024 = z.fma(term0231, term0221.mul(3));
300     R041 = z.mul(term0401);
301     DoubleVector term0411 = z.mul(term0402);
302     R042 = z.fma(term0411, term0401);
303     R201 = z.mul(term0002);
304     DoubleVector term2011 = z.mul(term0003);
305     R202 = z.fma(term2011, term0002);
306     DoubleVector term2012 = z.mul(term0004);
307     DoubleVector term2021 = z.fma(term2012, term0003);
308     R203 = z.fma(term2021, term2011.mul(2));
309     DoubleVector term2013 = z.mul(term0005);
310     DoubleVector term2022 = z.fma(term2013, term0004);
311     DoubleVector term2031 = z.fma(term2022, term2012.mul(2));
312     R204 = z.fma(term2031, term2021.mul(3));
313     R221 = z.mul(term0003);
314     DoubleVector term2211 = z.mul(term0004);
315     R222 = z.fma(term2211, term0003);
316     R401 = z.mul(term4001);
317     DoubleVector term4011 = z.mul(term4002);
318     R402 = z.fma(term4011, term4001);
319   }
320 
321   /**
322    * {@inheritDoc}
323    */
324   @SuppressWarnings("fallthrough")
325   @Override
326   protected void multipoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order) {
327     switch (order) {
328       default:
329       case 3:
330         DoubleVector term300 = mK.dx.mul(R400);
331         term300 = mK.qxz.fma(R401, term300);
332         E300 = term300.neg();
333         DoubleVector term030 = mK.dy.mul(R040);
334         term030 = mK.qyz.fma(R041, term030);
335         E030 = term030.neg();
336         DoubleVector term003 = mK.q.mul(R003);
337         term003 = mK.dz.fma(R004, term003);
338         term003 = mK.qxx.fma(R203, term003);
339         term003 = mK.qyy.fma(R023, term003);
340         term003 = mK.qzz.fma(R005, term003);
341         E003 = term003.neg();
342         DoubleVector term210 = mK.dy.mul(R220);
343         term210 = mK.qyz.fma(R221, term210);
344         E210 = term210.neg();
345         DoubleVector term201 = mK.q.mul(R201);
346         term201 = mK.dz.fma(R202, term201);
347         term201 = mK.qxx.fma(R401, term201);
348         term201 = mK.qyy.fma(R221, term201);
349         term201 = mK.qzz.fma(R203, term201);
350         E201 = term201.neg();
351         DoubleVector term120 = mK.dx.mul(R220);
352         term120 = mK.qxz.fma(R221, term120);
353         E120 = term120.neg();
354         DoubleVector term021 = mK.q.mul(R021);
355         term021 = mK.dz.fma(R022, term021);
356         term021 = mK.qxx.fma(R221, term021);
357         term021 = mK.qyy.fma(R041, term021);
358         term021 = mK.qzz.fma(R023, term021);
359         E021 = term021.neg();
360         DoubleVector term102 = mK.dx.mul(R202);
361         term102 = mK.qxz.fma(R203, term102);
362         E102 = term102.neg();
363         DoubleVector term012 = mK.dy.mul(R022);
364         term012 = mK.qyz.fma(R023, term012);
365         E012 = term012.neg();
366         DoubleVector term111 = mK.qxy.mul(R221);
367         E111 = term111.neg();
368       case 2:
369         DoubleVector term200 = mK.q.mul(R200);
370         term200 = mK.dz.fma(R201, term200);
371         term200 = mK.qxx.fma(R400, term200);
372         term200 = mK.qyy.fma(R220, term200);
373         term200 = mK.qzz.fma(R202, term200);
374         E200 = term200;
375         DoubleVector term020 = mK.q.mul(R020);
376         term020 = mK.dz.fma(R021, term020);
377         term020 = mK.qxx.fma(R220, term020);
378         term020 = mK.qyy.fma(R040, term020);
379         term020 = mK.qzz.fma(R022, term020);
380         E020 = term020;
381         DoubleVector term002 = mK.q.mul(R002);
382         term002 = mK.dz.fma(R003, term002);
383         term002 = mK.qxx.fma(R202, term002);
384         term002 = mK.qyy.fma(R022, term002);
385         term002 = mK.qzz.fma(R004, term002);
386         E002 = term002;
387         E110 = mK.qxy.mul(R220);
388         DoubleVector term101 = mK.dx.mul(R201);
389         term101 = mK.qxz.fma(R202, term101);
390         E101 = term101;
391         DoubleVector term011 = mK.dy.mul(R021);
392         term011 = mK.qyz.fma(R022, term011);
393         E011 = term011;
394       case 1:
395         DoubleVector term100 = mK.dx.mul(R200);
396         term100 = mK.qxz.fma(R201, term100);
397         E100 = term100.neg();
398         DoubleVector term010 = mK.dy.mul(R020);
399         term010 = mK.qyz.fma(R021, term010);
400         E010 = term010.neg();
401         DoubleVector term001 = mK.q.mul(R001);
402         term001 = mK.dz.fma(R002, term001);
403         term001 = mK.qxx.fma(R201, term001);
404         term001 = mK.qyy.fma(R021, term001);
405         term001 = mK.qzz.fma(R003, term001);
406         E001 = term001.neg();
407       case 0:
408         DoubleVector term000 = mK.q.mul(R000);
409         term000 = mK.dz.fma(R001, term000);
410         term000 = mK.qxx.fma(R200, term000);
411         term000 = mK.qyy.fma(R020, term000);
412         term000 = mK.qzz.fma(R002, term000);
413         E000 = term000;
414     }
415   }
416 
417   /**
418    * {@inheritDoc}
419    */
420   @SuppressWarnings("fallthrough")
421   @Override
422   protected void chargeKPotentialAtI(PolarizableMultipoleSIMD mK, int order) {
423     switch (order) {
424       default:
425       case 3:
426         // Order 3
427         E300 = zero;
428         E030 = zero;
429         E003 = mK.q.mul(R003).neg();
430         E210 = zero;
431         E201 = mK.q.mul(R201).neg();
432         E120 = zero;
433         E021 = mK.q.mul(R021).neg();
434         E102 = zero;
435         E012 = zero;
436         E111 = zero;
437         // Fall through to 2nd order.
438       case 2:
439         // Order 2
440         E200 = mK.q.mul(R200);
441         E020 = mK.q.mul(R020);
442         E002 = mK.q.mul(R002);
443         E110 = zero;
444         E101 = zero;
445         E011 = zero;
446         // Fall through to 1st order.
447       case 1:
448         // Order 1
449         E100 = zero;
450         E010 = zero;
451         E001 = mK.q.mul(R001).neg();
452       case 0:
453         E000 = mK.q.mul(R000);
454     }
455   }
456 
457   /**
458    * {@inheritDoc}
459    */
460   @SuppressWarnings("fallthrough")
461   @Override
462   protected void dipoleKPotentialAtI(DoubleVector uxk, DoubleVector uyk, DoubleVector uzk, int order) {
463     switch (order) {
464       default:
465       case 3:
466         // Order 3
467         E300 = uxk.mul(R400).neg();
468         E030 = uyk.mul(R040).neg();
469         E003 = uzk.mul(R004).neg();
470         E210 = uyk.mul(R220).neg();
471         E201 = uzk.mul(R202).neg();
472         E120 = uxk.mul(R220).neg();
473         E021 = uzk.mul(R022).neg();
474         E102 = uxk.mul(R202).neg();
475         E012 = uyk.mul(R022).neg();
476         E111 = zero;
477         // Fall through to 2nd order.
478       case 2:
479         // Order 2
480         E200 = uzk.mul(R201);
481         E020 = uzk.mul(R021);
482         E002 = uzk.mul(R003);
483         E110 = zero;
484         E101 = uxk.mul(R201);
485         E011 = uyk.mul(R021);
486         // Fall through to 1st order.
487       case 1:
488         // Order 1
489         E100 = uxk.mul(R200).neg();
490         E010 = uyk.mul(R020).neg();
491         E001 = uzk.mul(R002).neg();
492         // Fall through to the potential.
493       case 0:
494         E000 = uzk.mul(R001);
495     }
496   }
497 
498   /**
499    * {@inheritDoc}
500    */
501   @SuppressWarnings("fallthrough")
502   @Override
503   protected void quadrupoleKPotentialAtI(PolarizableMultipoleSIMD mK, int order) {
504     switch (order) {
505       default:
506       case 3:
507         // Order 3
508         E300 = mK.qxz.mul(R401).neg();
509         E030 = mK.qyz.mul(R041).neg();
510         DoubleVector term003 = mK.qxx.mul(R203);
511         term003 = mK.qyy.fma(R023, term003);
512         term003 = mK.qzz.fma(R005, term003);
513         E003 = term003;
514         E210 = mK.qyz.mul(R221).neg();
515         DoubleVector term201 = mK.qxx.mul(R401);
516         term201 = mK.qyy.fma(R221, term201);
517         term201 = mK.qzz.fma(R203, term201);
518         E201 = term201.neg();
519         E120 = mK.qxz.mul(R221).neg();
520         DoubleVector term021 = mK.qxx.mul(R221);
521         term021 = mK.qyy.fma(R041, term021);
522         term021 = mK.qzz.fma(R023, term021);
523         E021 = term021.neg();
524         E102 = mK.qxz.mul(R203).neg();
525         E012 = mK.qyz.mul(R023).neg();
526         E111 = mK.qxy.mul(R221).neg();
527         // Fall through to 2nd order.
528       case 2:
529         // Order 2
530         DoubleVector term200 = mK.qxx.mul(R400);
531         term200 = mK.qyy.fma(R220, term200);
532         term200 = mK.qzz.fma(R202, term200);
533         E200 = term200;
534         DoubleVector term020 = mK.qxx.mul(R220);
535         term020 = mK.qyy.fma(R040, term020);
536         term020 = mK.qzz.fma(R022, term020);
537         E020 = term020;
538         DoubleVector term002 = mK.qxx.mul(R202);
539         term002 = mK.qyy.fma(R022, term002);
540         term002 = mK.qzz.fma(R004, term002);
541         E002 = term002;
542         E110 = mK.qxy.mul(R220);
543         E101 = mK.qxz.mul(R202);
544         E011 = mK.qyz.mul(R022);
545         // Fall through to 1st order.
546       case 1:
547         // Order 1
548         E100 = mK.qxz.mul(R201).neg();
549         E010 = mK.qyz.mul(R021).neg();
550         DoubleVector term001 = mK.qxx.mul(R201);
551         term001 = mK.qyy.fma(R021, term001);
552         term001 = mK.qzz.fma(R003, term001);
553         E001 = term001.neg();
554       case 0:
555         DoubleVector term000 = mK.qxx.mul(R200);
556         term000 = mK.qyy.fma(R020, term000);
557         term000 = mK.qzz.fma(R002, term000);
558         E000 = term000;
559     }
560   }
561 
562   /**
563    * {@inheritDoc}
564    */
565   @SuppressWarnings("fallthrough")
566   @Override
567   protected void multipoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order) {
568     switch (order) {
569       default:
570       case 3:
571         DoubleVector term300 = mI.dx.mul(R400.neg());
572         term300 = mI.qxz.fma(R401, term300);
573         E300 = term300;
574         DoubleVector term030 = mI.dy.mul(R040.neg());
575         term030 = mI.qyz.fma(R041, term030);
576         E030 = term030;
577         DoubleVector term003 = mI.q.mul(R003);
578         term003 = mI.dz.fma(R004.neg(), term003);
579         term003 = mI.qxx.fma(R203, term003);
580         term003 = mI.qyy.fma(R023, term003);
581         term003 = mI.qzz.fma(R005, term003);
582         E003 = term003;
583         DoubleVector term210 = mI.dy.mul(R220.neg());
584         term210 = mI.qyz.fma(R221, term210);
585         E210 = term210;
586         DoubleVector term201 = mI.q.mul(R201);
587         term201 = mI.dz.fma(R202.neg(), term201);
588         term201 = mI.qxx.fma(R401, term201);
589         term201 = mI.qyy.fma(R221, term201);
590         term201 = mI.qzz.fma(R203, term201);
591         E201 = term201;
592         DoubleVector term120 = mI.dx.mul(R220.neg());
593         term120 = mI.qxz.fma(R221, term120);
594         E120 = term120;
595         DoubleVector term021 = mI.q.mul(R021);
596         term021 = mI.dz.fma(R022.neg(), term021);
597         term021 = mI.qxx.fma(R221, term021);
598         term021 = mI.qyy.fma(R041, term021);
599         term021 = mI.qzz.fma(R023, term021);
600         E021 = term021;
601         DoubleVector term102 = mI.dx.mul(R202.neg());
602         term102 = mI.qxz.fma(R203, term102);
603         E102 = term102;
604         DoubleVector term012 = mI.dy.mul(R022.neg());
605         term012 = mI.qyz.fma(R023, term012);
606         E012 = term012;
607         E111 = mI.qxy.mul(R221);
608       case 2:
609         DoubleVector term200 = mI.q.mul(R200);
610         term200 = mI.dz.fma(R201.neg(), term200);
611         term200 = mI.qxx.fma(R400, term200);
612         term200 = mI.qyy.fma(R220, term200);
613         term200 = mI.qzz.fma(R202, term200);
614         E200 = term200;
615         DoubleVector term020 = mI.q.mul(R020);
616         term020 = mI.dz.fma(R021.neg(), term020);
617         term020 = mI.qxx.fma(R220, term020);
618         term020 = mI.qyy.fma(R040, term020);
619         term020 = mI.qzz.fma(R022, term020);
620         E020 = term020;
621         DoubleVector term002 = mI.q.mul(R002);
622         term002 = mI.dz.fma(R003.neg(), term002);
623         term002 = mI.qxx.fma(R202, term002);
624         term002 = mI.qyy.fma(R022, term002);
625         term002 = mI.qzz.fma(R004, term002);
626         E002 = term002;
627         E110 = mI.qxy.mul(R220);
628         DoubleVector term101 = mI.dx.mul(R201.neg());
629         term101 = mI.qxz.fma(R202, term101);
630         E101 = term101;
631         DoubleVector term011 = mI.dy.mul(R021.neg());
632         term011 = mI.qyz.fma(R022, term011);
633         E011 = term011;
634       case 1:
635         DoubleVector term100 = mI.dx.mul(R200.neg());
636         term100 = mI.qxz.fma(R201, term100);
637         E100 = term100;
638         DoubleVector term010 = mI.dy.mul(R020.neg());
639         term010 = mI.qyz.fma(R021, term010);
640         E010 = term010;
641         DoubleVector term001 = mI.q.mul(R001);
642         term001 = mI.dz.fma(R002.neg(), term001);
643         term001 = mI.qxx.fma(R201, term001);
644         term001 = mI.qyy.fma(R021, term001);
645         term001 = mI.qzz.fma(R003, term001);
646         E001 = term001;
647       case 0:
648         DoubleVector term000 = mI.q.mul(R000);
649         term000 = mI.dz.fma(R001.neg(), term000);
650         term000 = mI.qxx.fma(R200, term000);
651         term000 = mI.qyy.fma(R020, term000);
652         term000 = mI.qzz.fma(R002, term000);
653         E000 = term000;
654     }
655   }
656 
657   /**
658    * {@inheritDoc}
659    */
660   @SuppressWarnings("fallthrough")
661   @Override
662   protected void chargeIPotentialAtK(PolarizableMultipoleSIMD mI, int order) {
663     switch (order) {
664       default:
665       case 3:
666         // Order 3.
667         E300 = zero;
668         E030 = zero;
669         E003 = mI.q.mul(R003);
670         E210 = zero;
671         E201 = mI.q.mul(R201);
672         E120 = zero;
673         E021 = mI.q.mul(R021);
674         E102 = zero;
675         E012 = zero;
676         E111 = zero;
677         // Fall through to 2nd order.
678       case 2:
679         // Order 2.
680         E200 = mI.q.mul(R200);
681         E020 = mI.q.mul(R020);
682         E002 = mI.q.mul(R002);
683         E110 = zero;
684         E101 = zero;
685         E011 = zero;
686         // Fall through to 1st order.
687       case 1:
688         // Order 1.
689         E100 = zero;
690         E010 = zero;
691         E001 = mI.q.mul(R001);
692         // Fall through to the potential.
693       case 0:
694         E000 = mI.q.mul(R000);
695     }
696   }
697 
698   /**
699    * {@inheritDoc}
700    */
701   @SuppressWarnings("fallthrough")
702   @Override
703   protected void dipoleIPotentialAtK(DoubleVector uxi, DoubleVector uyi, DoubleVector uzi, int order) {
704     switch (order) {
705       default:
706       case 3:
707         // Order 3
708         E300 = uxi.mul(R400).neg();
709         E030 = uyi.mul(R040).neg();
710         E003 = uzi.mul(R004).neg();
711         E210 = uyi.mul(R220).neg();
712         E201 = uzi.mul(R202).neg();
713         E120 = uxi.mul(R220).neg();
714         E021 = uzi.mul(R022).neg();
715         E102 = uxi.mul(R202).neg();
716         E012 = uyi.mul(R022).neg();
717         E111 = zero;
718         // Fall through to 2nd order.
719       case 2:
720         E200 = uzi.mul(R201.neg());
721         E020 = uzi.mul(R021.neg());
722         E002 = uzi.mul(R003.neg());
723         E110 = zero;
724         E101 = uxi.mul(R201.neg());
725         E011 = uyi.mul(R021.neg());
726         // Fall through to 1st order.
727       case 1:
728         E100 = uxi.mul(R200.neg());
729         E010 = uyi.mul(R020.neg());
730         E001 = uzi.mul(R002.neg());
731         // Fall through to the potential.
732       case 0:
733         E000 = uzi.mul(R001.neg());
734     }
735   }
736 
737   /**
738    * {@inheritDoc}
739    */
740   @SuppressWarnings("fallthrough")
741   @Override
742   protected void quadrupoleIPotentialAtK(PolarizableMultipoleSIMD mI, int order) {
743     switch (order) {
744       default:
745       case 3:
746         // Order 3.
747         E300 = mI.qxz.mul(R401);
748         E030 = mI.qyz.mul(R041);
749         DoubleVector term003 = mI.qxx.mul(R203);
750         term003 = mI.qyy.fma(R023, term003);
751         term003 = mI.qzz.fma(R005, term003);
752         E003 = term003;
753         E210 = mI.qyz.mul(R221);
754         DoubleVector term201 = mI.qxx.mul(R401);
755         term201 = mI.qyy.fma(R221, term201);
756         term201 = mI.qzz.fma(R203, term201);
757         E201 = term201;
758         E120 = mI.qxz.mul(R221);
759         DoubleVector term021 = mI.qxx.mul(R221);
760         term021 = mI.qyy.fma(R041, term021);
761         term021 = mI.qzz.fma(R023, term021);
762         E021 = term021;
763         E102 = mI.qxz.mul(R203);
764         E012 = mI.qyz.mul(R023);
765         E111 = mI.qxy.mul(R221);
766         // Fall through to 2nd order.
767       case 2:
768         // Order 2.
769         DoubleVector term200 = mI.qxx.mul(R400);
770         term200 = mI.qyy.fma(R220, term200);
771         term200 = mI.qzz.fma(R202, term200);
772         E200 = term200;
773         DoubleVector term020 = mI.qxx.mul(R220);
774         term020 = mI.qyy.fma(R040, term020);
775         term020 = mI.qzz.fma(R022, term020);
776         E020 = term020;
777         DoubleVector term002 = mI.qxx.mul(R202);
778         term002 = mI.qyy.fma(R022, term002);
779         term002 = mI.qzz.fma(R004, term002);
780         E002 = term002;
781         E110 = mI.qxy.mul(R220);
782         E101 = mI.qxz.mul(R202);
783         E011 = mI.qyz.mul(R022);
784         // Fall through to 1st order.
785       case 1:
786         // Order 1.
787         E100 = mI.qxz.mul(R201);
788         E010 = mI.qyz.mul(R021);
789         DoubleVector term001 = mI.qxx.mul(R201);
790         term001 = mI.qyy.fma(R021, term001);
791         term001 = mI.qzz.fma(R003, term001);
792         E001 = term001;
793         // Fall through to the potential.
794       case 0:
795         DoubleVector term000 = mI.qxx.mul(R200);
796         term000 = mI.qyy.fma(R020, term000);
797         term000 = mI.qzz.fma(R002, term000);
798         E000 = term000;
799     }
800   }
801 }