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