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 GeneralizedKirkwoodTensor class contains utilities for generated Generalized Kirkwood
44 * interaction tensors.
45 *
46 * @author Michael J. Schnieders
47 * @since 1.0
48 */
49 public class GKTensorQISIMD extends CoulombTensorQISIMD {
50
51 /**
52 * The GK tensor can be constructed for monopoles (GB), dipoles or quadrupoles.
53 */
54 protected final GKMultipoleOrder multipoleOrder;
55
56 /**
57 * The Kirkwood dielectric function for the given multipole order.
58 */
59 private final double c;
60
61 private final GKSourceSIMD gkSource;
62
63 private final DoubleVector zero = DoubleVector.zero(DoubleVector.SPECIES_PREFERRED);
64
65 /**
66 * Create a new GKTensorQI object.
67 *
68 * @param multipoleOrder The multipole order.
69 * @param order The tensor order.
70 * @param gkSource Generate the source terms for the GK recurrence.
71 * @param Eh Homogeneous dielectric constant.
72 * @param Es Solvent dielectric constant.
73 */
74 public GKTensorQISIMD(GKMultipoleOrder multipoleOrder, int order, GKSourceSIMD gkSource, double Eh, double Es) {
75 super(order);
76 this.multipoleOrder = multipoleOrder;
77 this.gkSource = gkSource;
78
79 // Load the dielectric function
80 c = GKSource.cn(multipoleOrder.getOrder(), Eh, Es);
81 }
82
83 /**
84 * GK Permanent multipole energy.
85 *
86 * @param mI PolarizableMultipole at site I.
87 * @param mK PolarizableMultipole at site K.
88 * @return the GK permanent multipole energy.
89 */
90 @Override
91 public DoubleVector multipoleEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
92 return switch (multipoleOrder) {
93 default -> {
94 chargeIPotentialAtK(mI, 2);
95 DoubleVector eK = multipoleEnergy(mK);
96 chargeKPotentialAtI(mK, 2);
97 DoubleVector eI = multipoleEnergy(mI);
98 yield eK.add(eI).mul(c * 0.5);
99 }
100 case DIPOLE -> {
101 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
102 DoubleVector eK = multipoleEnergy(mK);
103 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
104 DoubleVector eI = multipoleEnergy(mI);
105 yield eK.add(eI).mul(c * 0.5);
106 }
107 case QUADRUPOLE -> {
108 quadrupoleIPotentialAtK(mI, 2);
109 DoubleVector eK = multipoleEnergy(mK);
110 quadrupoleKPotentialAtI(mK, 2);
111 DoubleVector eI = multipoleEnergy(mI);
112 yield eK.add(eI).mul(c * 0.5);
113 }
114 };
115 }
116
117 /**
118 * GK Permanent multipole energy and gradient.
119 *
120 * @param mI PolarizableMultipole at site I.
121 * @param mK PolarizableMultipole at site K.
122 * @param Gi Coordinate gradient at site I.
123 * @param Gk Coordinate gradient at site K.
124 * @param Ti Torque at site I.
125 * @param Tk Torque at site K.
126 * @return the permanent multipole GK energy.
127 */
128 @Override
129 public DoubleVector multipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
130 DoubleVector[] Gi, DoubleVector[] Gk, DoubleVector[] Ti, DoubleVector[] Tk) {
131 return switch (multipoleOrder) {
132 default -> monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
133 case DIPOLE -> dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
134 case QUADRUPOLE -> quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
135 };
136 }
137
138 /**
139 * Permanent multipole energy and gradient using the GK monopole tensor.
140 *
141 * @param mI PolarizableMultipole at site I.
142 * @param mK PolarizableMultipole at site K.
143 * @param Gi Coordinate gradient at site I.
144 * @param Gk Coordinate gradient at site K.
145 * @param Ti Torque at site I.
146 * @param Tk Torque at site K.
147 * @return the permanent multipole GK energy.
148 */
149 protected DoubleVector monopoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
150 DoubleVector[] Gi, DoubleVector[] Gk, DoubleVector[] Ti, DoubleVector[] Tk) {
151
152 // Compute the potential due to a multipole component at site I.
153 chargeIPotentialAtK(mI, 3);
154 DoubleVector eK = multipoleEnergy(mK);
155 multipoleGradient(mK, Gk);
156 multipoleTorque(mK, Tk);
157
158 // Compute the potential due to a multipole component at site K.
159 chargeKPotentialAtI(mK, 3);
160 DoubleVector eI = multipoleEnergy(mI);
161 multipoleGradient(mI, Gi);
162 multipoleTorque(mI, Ti);
163
164 double scale = c * 0.5;
165 Gi[0] = Gi[0].sub(Gk[0]).mul(scale);
166 Gi[1] = Gi[1].sub(Gk[1]).mul(scale);
167 Gi[2] = Gi[2].sub(Gk[2]).mul(scale);
168 Gk[0] = Gi[0].neg();
169 Gk[1] = Gi[1].neg();
170 Gk[2] = Gi[2].neg();
171
172 Ti[0] = Ti[0].mul(scale);
173 Ti[1] = Ti[1].mul(scale);
174 Ti[2] = Ti[2].mul(scale);
175 Tk[0] = Tk[0].mul(scale);
176 Tk[1] = Tk[1].mul(scale);
177 Tk[2] = Tk[2].mul(scale);
178
179 return eK.add(eI).mul(scale);
180 }
181
182 /**
183 * Permanent multipole energy and gradient using the GK dipole tensor.
184 *
185 * @param mI PolarizableMultipole at site I.
186 * @param mK PolarizableMultipole at site K.
187 * @param Gi Coordinate gradient at site I.
188 * @param Gk Coordinate gradient at site K.
189 * @param Ti Torque at site I.
190 * @param Tk Torque at site K.
191 * @return the permanent multipole GK energy.
192 */
193 protected DoubleVector dipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
194 DoubleVector[] Gi, DoubleVector[] Gk, DoubleVector[] Ti, DoubleVector[] Tk) {
195
196 // Compute the potential due to a multipole component at site I.
197 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
198 DoubleVector eK = multipoleEnergy(mK);
199 multipoleGradient(mK, Gk);
200 multipoleTorque(mK, Tk);
201
202 // Need the torque on site I pole due to site K multipole.
203 // Only torque on the site I dipole.
204 multipoleKPotentialAtI(mK, 1);
205 dipoleTorque(mI, Ti);
206
207 // Compute the potential due to a multipole component at site K.
208 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
209 DoubleVector eI = multipoleEnergy(mI);
210 multipoleGradient(mI, Gi);
211 multipoleTorque(mI, Ti);
212
213 // Need the torque on site K pole due to multipole on site I.
214 // Only torque on the site K dipole.
215 multipoleIPotentialAtK(mI, 1);
216 dipoleTorque(mK, Tk);
217
218 double scale = c * 0.5;
219 Gi[0] = Gi[0].sub(Gk[0]).mul(scale);
220 Gi[1] = Gi[1].sub(Gk[1]).mul(scale);
221 Gi[2] = Gi[2].sub(Gk[2]).mul(scale);
222 Gk[0] = Gi[0].neg();
223 Gk[1] = Gi[1].neg();
224 Gk[2] = Gi[2].neg();
225
226 Ti[0] = Ti[0].mul(scale);
227 Ti[1] = Ti[1].mul(scale);
228 Ti[2] = Ti[2].mul(scale);
229 Tk[0] = Tk[0].mul(scale);
230 Tk[1] = Tk[1].mul(scale);
231 Tk[2] = Tk[2].mul(scale);
232
233 return eK.add(eI).mul(scale);
234 }
235
236 /**
237 * Permanent multipole energy and gradient using the GK quadrupole tensor.
238 *
239 * @param mI PolarizableMultipole at site I.
240 * @param mK PolarizableMultipole at site K.
241 * @param Gi Coordinate gradient at site I.
242 * @param Gk Coordinate gradient at site K.
243 * @param Ti Torque at site I.
244 * @param Tk Torque at site K.
245 * @return the permanent multipole GK energy.
246 */
247 protected DoubleVector quadrupoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
248 DoubleVector[] Gi, DoubleVector[] Gk, DoubleVector[] Ti, DoubleVector[] Tk) {
249
250 // Compute the potential due to a multipole component at site I.
251 quadrupoleIPotentialAtK(mI, 3);
252 DoubleVector eK = multipoleEnergy(mK);
253 multipoleGradient(mK, Gk);
254 multipoleTorque(mK, Tk);
255
256 // Need the torque on site I quadrupole due to site K multipole.
257 multipoleKPotentialAtI(mK, 2);
258 quadrupoleTorque(mI, Ti);
259
260 // Compute the potential due to a multipole component at site K.
261 quadrupoleKPotentialAtI(mK, 3);
262 DoubleVector eI = multipoleEnergy(mI);
263 multipoleGradient(mI, Gi);
264 multipoleTorque(mI, Ti);
265
266 // Need the torque on site K quadrupole due to site I multipole.
267 multipoleIPotentialAtK(mI, 2);
268 quadrupoleTorque(mK, Tk);
269
270 double scale = c * 0.5;
271 Gi[0] = Gi[0].sub(Gk[0]).mul(scale);
272 Gi[1] = Gi[1].sub(Gk[1]).mul(scale);
273 Gi[2] = Gi[2].sub(Gk[2]).mul(scale);
274 Gk[0] = Gi[0].neg();
275 Gk[1] = Gi[1].neg();
276 Gk[2] = Gi[2].neg();
277
278 Ti[0] = Ti[0].mul(scale);
279 Ti[1] = Ti[1].mul(scale);
280 Ti[2] = Ti[2].mul(scale);
281 Tk[0] = Tk[0].mul(scale);
282 Tk[1] = Tk[1].mul(scale);
283 Tk[2] = Tk[2].mul(scale);
284
285 return eK.add(eI).mul(scale);
286 }
287
288 /**
289 * GK Permanent multipole Born grad.
290 *
291 * @param mI PolarizableMultipole at site I.
292 * @param mK PolarizableMultipole at site K.
293 * @return a double.
294 */
295 public DoubleVector multipoleEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
296 return multipoleEnergy(mI, mK);
297 }
298
299 /**
300 * GK Polarization Energy.
301 *
302 * @param mI PolarizableMultipole at site I.
303 * @param mK PolarizableMultipole at site K.
304 * @param scaleEnergy This is ignored, since masking/scaling is not applied to GK
305 * interactions.
306 * @return a double.
307 */
308 public DoubleVector polarizationEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK, DoubleVector scaleEnergy) {
309 return polarizationEnergy(mI, mK);
310 }
311
312 /**
313 * GK Polarization Energy.
314 *
315 * @param mI PolarizableMultipole at site I.
316 * @param mK PolarizableMultipole at site K.
317 * @return a double.
318 */
319 public DoubleVector polarizationEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
320 return switch (multipoleOrder) {
321 default -> {
322 // Find the GK charge potential of site I at site K.
323 chargeIPotentialAtK(mI, 1);
324 // Energy of induced dipole K in the field of permanent charge I.
325 DoubleVector eK = polarizationEnergy(mK);
326 // Find the GK charge potential of site K at site I.
327 chargeKPotentialAtI(mK, 1);
328 // Energy of induced dipole I in the field of permanent charge K.
329 DoubleVector eI = polarizationEnergy(mI);
330 yield eK.add(eI).mul(c * 0.5);
331 }
332 case DIPOLE -> {
333 // Find the GK dipole potential of site I at site K.
334 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
335 // Energy of induced dipole K in the field of permanent dipole I.
336 DoubleVector eK = polarizationEnergy(mK);
337 // Find the GK induced dipole potential of site I at site K.
338 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
339 // Energy of permanent multipole K in the field of induced dipole I.
340 eK = multipoleEnergy(mK).mul(0.5).add(eK);
341 // Find the GK dipole potential of site K at site I.
342 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
343 // Energy of induced dipole I in the field of permanent dipole K.
344 DoubleVector eI = polarizationEnergy(mI);
345 // Find the GK induced dipole potential of site K at site I.
346 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
347 // Energy of permanent multipole I in the field of induced dipole K.
348 eI = multipoleEnergy(mI).mul(0.5).add(eI);
349 yield eK.add(eI).mul(c * 0.5);
350 }
351 case QUADRUPOLE -> {
352 // Find the GK quadrupole potential of site I at site K.
353 quadrupoleIPotentialAtK(mI, 1);
354 // Energy of induced dipole K in the field of permanent quadrupole I.
355 DoubleVector eK = polarizationEnergy(mK);
356 // Find the GK quadrupole potential of site K at site I.
357 quadrupoleKPotentialAtI(mK, 1);
358 // Energy of induced dipole I in the field of permanent quadrupole K.
359 DoubleVector eI = polarizationEnergy(mI);
360 yield eK.add(eI).mul(c * 0.5);
361 }
362 };
363 }
364
365 /**
366 * GK Polarization Energy.
367 *
368 * @param mI PolarizableMultipole at site I.
369 * @param mK PolarizableMultipole at site K.
370 * @return a double.
371 */
372 public DoubleVector polarizationEnergyBorn(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
373 return switch (multipoleOrder) {
374 default -> {
375 // Find the GK charge potential of site I at site K.
376 chargeIPotentialAtK(mI, 1);
377 // Energy of induced dipole K in the field of permanent charge I.
378 DoubleVector eK = polarizationEnergyS(mK);
379 // Find the GK charge potential of site K at site I.
380 chargeKPotentialAtI(mK, 1);
381 // Energy of induced dipole I in the field of permanent charge K.
382 DoubleVector eI = polarizationEnergyS(mI);
383 yield eK.add(eI).mul(c * 0.5);
384 }
385 case DIPOLE -> {
386 // Find the GK dipole potential of site I at site K.
387 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
388 // Energy of induced dipole K in the field of permanent dipole I.
389 DoubleVector eK = polarizationEnergyS(mK);
390 // Find the GK induced dipole potential of site I at site K.
391 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
392 // Energy of permanent multipole K in the field of induced dipole I.
393 eK = multipoleEnergy(mK).mul(0.5).add(eK);
394 // Find the GK dipole potential of site K at site I.
395 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
396 // Energy of induced dipole I in the field of permanent dipole K.
397 DoubleVector eI = polarizationEnergyS(mI);
398 // Find the GK induced dipole potential of site K at site I.
399 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
400 // Energy of permanent multipole I in the field of induced dipole K.
401 eI = multipoleEnergy(mI).mul(0.5).add(eI);
402 yield eK.add(eI).mul(c * 0.5);
403 }
404 case QUADRUPOLE -> {
405 // Find the GK quadrupole potential of site I at site K.
406 quadrupoleIPotentialAtK(mI, 1);
407 // Energy of induced dipole K in the field of permanent quadrupole I.
408 DoubleVector eK = polarizationEnergyS(mK);
409 // Find the GK quadrupole potential of site K at site I.
410 quadrupoleKPotentialAtI(mK, 1);
411 // Energy of induced dipole I in the field of permanent quadrupole K.
412 DoubleVector eI = polarizationEnergyS(mI);
413 yield eK.add(eI).mul(c * 0.5);
414 }
415 };
416 }
417
418 /**
419 * Polarization Energy and Gradient.
420 *
421 * @param mI PolarizableMultipole at site I.
422 * @param mK PolarizableMultipole at site K.
423 * @param inductionMask This is ignored, since masking/scaling is not applied to GK
424 * interactions (everything is intermolecular).
425 * @param energyMask This is ignored, since masking/scaling is not applied to GK interactions
426 * (everything is intermolecular).
427 * @param mutualMask This should be set to zero for direction polarization.
428 * @param Gi an array of {@link double} objects.
429 * @param Ti an array of {@link double} objects.
430 * @param Tk an array of {@link double} objects.
431 * @return a double.
432 */
433 @Override
434 public DoubleVector polarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
435 DoubleVector inductionMask, DoubleVector energyMask, DoubleVector mutualMask,
436 DoubleVector[] Gi, DoubleVector[] Ti, DoubleVector[] Tk) {
437 return switch (multipoleOrder) {
438 default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
439 case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
440 case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
441 };
442 }
443
444 /**
445 * Monopole Polarization Energy and Gradient.
446 *
447 * @param mI PolarizableMultipole at site I.
448 * @param mK PolarizableMultipole at site K.
449 * @param Gi an array of {@link double} objects.
450 * @return a double.
451 */
452 public DoubleVector monopolePolarizationEnergyAndGradient(PolarizableMultipoleSIMD mI,
453 PolarizableMultipoleSIMD mK, DoubleVector[] Gi) {
454 // Find the permanent multipole potential at site k.
455 chargeIPotentialAtK(mI, 2);
456 // Energy of induced dipole k in the field of multipole i.
457 DoubleVector eK = polarizationEnergy(mK);
458 // Derivative with respect to moving atom k.
459 Gi[0] = mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101)).neg();
460 Gi[1] = mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011)).neg();
461 Gi[2] = mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002)).neg();
462
463 // Find the permanent multipole potential and derivatives at site i.
464 chargeKPotentialAtI(mK, 2);
465 // Energy of induced dipole i in the field of multipole k.
466 DoubleVector eI = polarizationEnergy(mI);
467 // Derivative with respect to moving atom i.
468 Gi[0] = Gi[0].add(mI.sx.mul(E200)).add(mI.sy.mul(E110)).add(mI.sz.mul(E101));
469 Gi[1] = Gi[1].add(mI.sx.mul(E110)).add(mI.sy.mul(E020)).add(mI.sz.mul(E011));
470 Gi[2] = Gi[2].add(mI.sx.mul(E101)).add(mI.sy.mul(E011)).add(mI.sz.mul(E002));
471
472 double scale = c * 0.5;
473 Gi[0] = Gi[0].mul(scale);
474 Gi[1] = Gi[1].mul(scale);
475 Gi[2] = Gi[2].mul(scale);
476
477 // Total polarization energy.
478 return eI.add(eK).mul(scale);
479 }
480
481 /**
482 * Dipole Polarization Energy and Gradient.
483 *
484 * @param mI PolarizableMultipole at site I.
485 * @param mK PolarizableMultipole at site K.
486 * @param mutualMask This should be set to zero for direction polarization.
487 * @param Gi an array of {@link double} objects.
488 * @param Ti an array of {@link double} objects.
489 * @param Tk an array of {@link double} objects.
490 * @return a double.
491 */
492 public DoubleVector dipolePolarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
493 DoubleVector mutualMask, DoubleVector[] Gi,
494 DoubleVector[] Ti, DoubleVector[] Tk) {
495
496 // Find the permanent multipole potential and derivatives at site k.
497 dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
498 // Energy of induced dipole k in the field of multipole i.
499 DoubleVector eK = polarizationEnergy(mK);
500 // Derivative with respect to moving atom k.
501 Gi[0] = mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101)).neg();
502 Gi[1] = mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011)).neg();
503 Gi[2] = mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002)).neg();
504
505 // Find the potential at K due to the averaged induced dipole at site i.
506 dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
507 dipoleTorque(mK, Tk);
508
509 // Find the GK induced dipole potential of site I at site K.
510 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 3);
511 // Energy of permanent multipole K in the field of induced dipole I.
512 eK = eK.add(multipoleEnergy(mK).mul(0.5));
513
514 DoubleVector[] G = new DoubleVector[3];
515 multipoleGradient(mK, G);
516 Gi[0] = Gi[0].sub(G[0]);
517 Gi[1] = Gi[1].sub(G[1]);
518 Gi[2] = Gi[2].sub(G[2]);
519 multipoleTorque(mK, Tk);
520
521 // Find the permanent multipole potential and derivatives at site i.
522 dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
523 // Energy of induced dipole i in the field of multipole k.
524 DoubleVector eI = polarizationEnergy(mI);
525 // Derivative with respect to moving atom i.
526 Gi[0] = Gi[0].add(mI.sx.mul(E200)).add(mI.sy.mul(E110)).add(mI.sz.mul(E101));
527 Gi[1] = Gi[1].add(mI.sx.mul(E110)).add(mI.sy.mul(E020)).add(mI.sz.mul(E011));
528 Gi[2] = Gi[2].add(mI.sx.mul(E101)).add(mI.sy.mul(E011)).add(mI.sz.mul(E002));
529
530 // Find the potential at I due to the averaged induced dipole at k.
531 dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
532 dipoleTorque(mI, Ti);
533
534 // Find the GK induced dipole potential of site K at site I.
535 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 3);
536 // Energy of permanent multipole I in the field of induced dipole K.
537 eI = eI.add(multipoleEnergy(mI).mul(0.5));
538
539 multipoleGradient(mI, G);
540 Gi[0] = Gi[0].add(G[0]);
541 Gi[1] = Gi[1].add(G[1]);
542 Gi[2] = Gi[2].add(G[2]);
543 multipoleTorque(mI, Ti);
544
545 // Get the induced-induced portion of the force (Ud . dC/dX . Up).
546 // This contribution does not exist for direct polarization (mutualMask == 0.0).
547 if (!mutualMask.eq(zero).allTrue()) {
548 // Find the potential and its derivatives at k due to induced dipole i.
549 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
550 Gi[0] = Gi[0].sub((mK.px.mul(E200).add(mK.py.mul(E110)).add(mK.pz.mul(E101)).mul(mutualMask)));
551 Gi[1] = Gi[1].sub((mK.px.mul(E110).add(mK.py.mul(E020)).add(mK.pz.mul(E011)).mul(mutualMask)));
552 Gi[2] = Gi[2].sub((mK.px.mul(E101).add(mK.py.mul(E011)).add(mK.pz.mul(E002)).mul(mutualMask)));
553
554 // Find the potential and its derivatives at i due to induced dipole k.
555 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
556 Gi[0] = Gi[0].sub((mI.px.mul(E200).add(mI.py.mul(E110)).add(mI.pz.mul(E101)).mul(mutualMask)));
557 Gi[1] = Gi[1].sub((mI.px.mul(E110).add(mI.py.mul(E020)).add(mI.pz.mul(E011)).mul(mutualMask)));
558 Gi[2] = Gi[2].sub((mI.px.mul(E101).add(mI.py.mul(E011)).add(mI.pz.mul(E002)).mul(mutualMask)));
559 }
560
561 // Total polarization energy.
562 double scale = c * 0.5;
563 DoubleVector energy = eI.add(eK).mul(scale);
564 Gi[0] = Gi[0].mul(scale);
565 Gi[1] = Gi[1].mul(scale);
566 Gi[2] = Gi[2].mul(scale);
567 Ti[0] = Ti[0].mul(scale);
568 Ti[1] = Ti[1].mul(scale);
569 Ti[2] = Ti[2].mul(scale);
570 Tk[0] = Tk[0].mul(scale);
571 Tk[1] = Tk[1].mul(scale);
572 Tk[2] = Tk[2].mul(scale);
573 return energy;
574 }
575
576 /**
577 * Quadrupole Polarization Energy and Gradient.
578 *
579 * @param mI PolarizableMultipole at site I.
580 * @param mK PolarizableMultipole at site K.
581 * @param Gi an array of {@link double} objects.
582 * @param Ti an array of {@link double} objects.
583 * @param Tk an array of {@link double} objects.
584 * @return a double.
585 */
586 public DoubleVector quadrupolePolarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
587 DoubleVector[] Gi, DoubleVector[] Ti, DoubleVector[] Tk) {
588
589 // Find the permanent multipole potential and derivatives at site k.
590 quadrupoleIPotentialAtK(mI, 2);
591 // Energy of induced dipole k in the field of multipole i.
592 DoubleVector eK = polarizationEnergy(mK);
593 // Derivative with respect to moving atom k.
594 Gi[0] = (mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101))).neg();
595 Gi[1] = (mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011))).neg();
596 Gi[2] = (mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002))).neg();
597
598 // Find the permanent multipole potential and derivatives at site i.
599 quadrupoleKPotentialAtI(mK, 2);
600 // Energy of induced dipole i in the field of multipole k.
601 DoubleVector eI = polarizationEnergy(mI);
602 // Derivative with respect to moving atom i.
603 Gi[0] = Gi[0].add(mI.sx.mul(E200)).add(mI.sy.mul(E110)).add(mI.sz.mul(E101));
604 Gi[1] = Gi[1].add(mI.sx.mul(E110)).add(mI.sy.mul(E020)).add(mI.sz.mul(E011));
605 Gi[2] = Gi[2].add(mI.sx.mul(E101)).add(mI.sy.mul(E011)).add(mI.sz.mul(E002));
606
607 double scale = c * 0.5;
608 Gi[0] = Gi[0].mul(scale);
609 Gi[1] = Gi[1].mul(scale);
610 Gi[2] = Gi[2].mul(scale);
611
612 // Find the potential and its derivatives at K due to the averaged induced dipole at site i.
613 dipoleIPotentialAtK(mI.sx.mul(scale), mI.sy.mul(scale), mI.sz.mul(scale), 2);
614 quadrupoleTorque(mK, Tk);
615
616 // Find the potential and its derivatives at I due to the averaged induced dipole at k.
617 dipoleKPotentialAtI(mK.sx.mul(scale), mK.sy.mul(scale), mK.sz.mul(scale), 2);
618 quadrupoleTorque(mI, Ti);
619
620 // Total polarization energy.
621 return eI.add(eK).mul(scale);
622 }
623
624 /**
625 * GK Direct Polarization Born grad.
626 *
627 * @param mI PolarizableMultipole at site I.
628 * @param mK PolarizableMultipole at site K.
629 * @return Partial derivative of the Polarization energy with respect to a Born grad.
630 */
631 public DoubleVector polarizationEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
632 return polarizationEnergyBorn(mI, mK).mul(2.0);
633 }
634
635 /**
636 * GK Mutual Polarization Contribution to the Born grad.
637 *
638 * @param mI PolarizableMultipole at site I.
639 * @param mK PolarizableMultipole at site K.
640 * @return Mutual Polarization contribution to the partial derivative with respect to a Born grad.
641 */
642 public DoubleVector mutualPolarizationEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
643 DoubleVector db = zero;
644 if (multipoleOrder == GKMultipoleOrder.DIPOLE) {
645 // Find the potential and its derivatives at k due to induced dipole i.
646 dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
647 db = mK.px.mul(E100).add(mK.py.mul(E010)).add(mK.pz.mul(E001)).mul(0.5);
648
649
650 // Find the potential and its derivatives at i due to induced dipole k.
651 dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
652 db = db.add(mI.px.mul(E100).add(mI.py.mul(E010)).add(mI.pz.mul(E001)).mul(0.5));
653
654 }
655 return db.mul(c);
656 }
657
658 /**
659 * Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
660 */
661 @Override
662 protected void source(DoubleVector[] work) {
663 gkSource.source(work, multipoleOrder);
664 }
665 }