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