1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.nonbonded.pme;
39
40 import static ffx.potential.parameters.MultipoleType.t000;
41 import static ffx.potential.parameters.MultipoleType.t001;
42 import static ffx.potential.parameters.MultipoleType.t002;
43 import static ffx.potential.parameters.MultipoleType.t003;
44 import static ffx.potential.parameters.MultipoleType.t010;
45 import static ffx.potential.parameters.MultipoleType.t011;
46 import static ffx.potential.parameters.MultipoleType.t012;
47 import static ffx.potential.parameters.MultipoleType.t020;
48 import static ffx.potential.parameters.MultipoleType.t021;
49 import static ffx.potential.parameters.MultipoleType.t030;
50 import static ffx.potential.parameters.MultipoleType.t100;
51 import static ffx.potential.parameters.MultipoleType.t101;
52 import static ffx.potential.parameters.MultipoleType.t102;
53 import static ffx.potential.parameters.MultipoleType.t110;
54 import static ffx.potential.parameters.MultipoleType.t111;
55 import static ffx.potential.parameters.MultipoleType.t120;
56 import static ffx.potential.parameters.MultipoleType.t200;
57 import static ffx.potential.parameters.MultipoleType.t201;
58 import static ffx.potential.parameters.MultipoleType.t210;
59 import static ffx.potential.parameters.MultipoleType.t300;
60 import static org.apache.commons.math3.util.FastMath.sqrt;
61
62 import edu.rit.pj.IntegerForLoop;
63 import edu.rit.pj.IntegerSchedule;
64 import edu.rit.pj.ParallelRegion;
65 import edu.rit.pj.ParallelTeam;
66 import edu.rit.pj.reduction.SharedDouble;
67 import ffx.crystal.Crystal;
68 import ffx.numerics.atomic.AtomicDoubleArray3D;
69 import ffx.numerics.multipole.MultipoleTensor;
70 import ffx.potential.bonded.Atom;
71 import ffx.potential.extended.ExtendedSystem;
72 import ffx.potential.nonbonded.ReciprocalSpace;
73 import java.util.logging.Level;
74 import java.util.logging.Logger;
75
76
77
78
79
80
81
82 public class ReciprocalEnergyRegion extends ParallelRegion {
83
84 private static final Logger logger = Logger.getLogger(ReciprocalEnergyRegion.class.getName());
85
86 private static final double SQRT_PI = sqrt(Math.PI);
87 private static final int tensorCount = MultipoleTensor.tensorCount(3);
88 private final double electric;
89 private final double aewald1;
90 private final double aewald2;
91 private final double aewald3;
92 private final double aewald4;
93 private final double oneThird;
94 private final double twoThirds;
95 private final int maxThreads;
96 private final SharedDouble inducedDipoleSelfEnergy;
97 private final SharedDouble inducedDipoleRecipEnergy;
98 private final PermanentReciprocalEnergyLoop[] permanentReciprocalEnergyLoop;
99 private final InducedDipoleReciprocalEnergyLoop[] inducedDipoleReciprocalEnergyLoop;
100
101 public double[][][] inducedDipole;
102 public double[][][] inducedDipoleCR;
103
104 private double[][] ind;
105
106 private double[][] indCR;
107 private double nfftX, nfftY, nfftZ;
108 private double[][] multipole;
109 private double[][] fracMultipole;
110 private double[][] fracMultipolePhi;
111 private double[][] fracInducedDipolePhi;
112 private double[][] fracInducedDipolePhiCR;
113 private double permanentSelfEnergy;
114 private double permanentReciprocalEnergy;
115
116 private Atom[] atoms;
117
118 private Crystal crystal;
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137 private boolean[] use;
138
139 private double[][][] globalMultipole;
140 private double[][][] globalFracMultipole;
141 private double[][][] titrationMultipole;
142 private double[][][] tautomerMultipole;
143 private ExtendedSystem extendedSystem = null;
144 private boolean esvTerm = false;
145 private double[][] cartMultipolePhi;
146 private double[][] cartesianDipolePhi;
147 private double[][] cartesianDipolePhiCR;
148
149 private ReciprocalSpace reciprocalSpace;
150 private Polarization polarization;
151
152 private AtomicDoubleArray3D grad;
153
154 private AtomicDoubleArray3D torque;
155
156 private AtomicDoubleArray3D lambdaGrad;
157
158 private AtomicDoubleArray3D lambdaTorque;
159
160 private boolean gradient;
161
162
163
164
165 private boolean lambdaTerm;
166
167 private SharedDouble shareddEdLambda;
168
169 private SharedDouble sharedd2EdLambda2;
170
171 private double permanentScale;
172 private double dlPowPerm;
173 private double d2lPowPerm;
174 private double polarizationScale;
175 private double dlPowPol;
176 private double d2lPowPol;
177 private double dEdLSign;
178
179 public ReciprocalEnergyRegion(int nt, double aewald, double electric) {
180 permanentReciprocalEnergyLoop = new PermanentReciprocalEnergyLoop[nt];
181 inducedDipoleReciprocalEnergyLoop = new InducedDipoleReciprocalEnergyLoop[nt];
182 inducedDipoleSelfEnergy = new SharedDouble();
183 inducedDipoleRecipEnergy = new SharedDouble();
184 maxThreads = nt;
185
186 this.electric = electric;
187 aewald1 = -electric * aewald / SQRT_PI;
188 aewald2 = 2.0 * aewald * aewald;
189 aewald3 = -2.0 / 3.0 * electric * aewald * aewald * aewald / SQRT_PI;
190 aewald4 = -2.0 * aewald3;
191 oneThird = 1.0 / 3.0;
192 twoThirds = 2.0 / 3.0;
193 }
194
195
196
197
198
199
200 public void executeWith(ParallelTeam parallelTeam) {
201 try {
202 parallelTeam.execute(this);
203 } catch (Exception e) {
204 String message = " Exception computing the electrostatic energy.\n";
205 logger.log(Level.SEVERE, message, e);
206 }
207 }
208
209 @Override
210 public void finish() {
211
212
213
214
215
216 permanentSelfEnergy = 0.0;
217 permanentReciprocalEnergy = 0.0;
218 for (int i = 0; i < maxThreads; i++) {
219 permanentSelfEnergy += permanentReciprocalEnergyLoop[i].eSelf;
220 permanentReciprocalEnergy += permanentReciprocalEnergyLoop[i].eRecip;
221 }
222 }
223
224 public double getInducedDipoleReciprocalEnergy() {
225 return inducedDipoleRecipEnergy.get();
226 }
227
228 public double getInducedDipoleSelfEnergy() {
229 return inducedDipoleSelfEnergy.get();
230 }
231
232 public double getPermanentReciprocalEnergy() {
233 return permanentReciprocalEnergy;
234 }
235
236 public double getPermanentSelfEnergy() {
237 return permanentSelfEnergy;
238 }
239
240 public void init(
241 Atom[] atoms,
242 Crystal crystal,
243 boolean gradient,
244 boolean lambdaTerm,
245 boolean esvTerm,
246 boolean[] use,
247 double[][][] globalMultipole,
248 double[][][] globalFracMultipole,
249 double[][][] titrationMultipole,
250 double[][][] tautomerMultipole,
251 double[][] cartMultipolePhi,
252 double[][] fracMultipolePhi,
253 Polarization polarization,
254 double[][][] inducedDipole,
255 double[][][] inducedDipoleCR,
256 double[][] cartesianDipolePhi,
257 double[][] cartesianDipolePhiCR,
258 double[][] fracInducedDipolePhi,
259 double[][] fracInducedDipolePhiCR,
260 ReciprocalSpace reciprocalSpace,
261 AlchemicalParameters alchemicalParameters,
262 ExtendedSystem extendedSystem,
263 AtomicDoubleArray3D grad,
264 AtomicDoubleArray3D torque,
265 AtomicDoubleArray3D lambdaGrad,
266 AtomicDoubleArray3D lambdaTorque,
267 SharedDouble shareddEdLambda,
268 SharedDouble sharedd2EdLambda2
269 ) {
270 this.atoms = atoms;
271 this.crystal = crystal;
272 this.gradient = gradient;
273 this.lambdaTerm = lambdaTerm;
274 this.esvTerm = esvTerm;
275 this.use = use;
276
277 this.globalMultipole = globalMultipole;
278 this.globalFracMultipole = globalFracMultipole;
279 this.titrationMultipole = titrationMultipole;
280 this.tautomerMultipole = tautomerMultipole;
281 this.cartMultipolePhi = cartMultipolePhi;
282 this.fracMultipolePhi = fracMultipolePhi;
283
284 this.polarization = polarization;
285 this.inducedDipole = inducedDipole;
286 this.inducedDipoleCR = inducedDipoleCR;
287 this.cartesianDipolePhi = cartesianDipolePhi;
288 this.cartesianDipolePhiCR = cartesianDipolePhiCR;
289 this.fracInducedDipolePhi = fracInducedDipolePhi;
290 this.fracInducedDipolePhiCR = fracInducedDipolePhiCR;
291 this.reciprocalSpace = reciprocalSpace;
292 this.extendedSystem = extendedSystem;
293
294 this.permanentScale = alchemicalParameters.permanentScale;
295 this.dlPowPerm = alchemicalParameters.dlPowPerm;
296 this.d2lPowPerm = alchemicalParameters.d2lPowPerm;
297 this.polarizationScale = alchemicalParameters.polarizationScale;
298 this.dlPowPol = alchemicalParameters.dlPowPol;
299 this.d2lPowPol = alchemicalParameters.d2lPowPol;
300 this.dEdLSign = alchemicalParameters.dEdLSign;
301
302 this.grad = grad;
303 this.torque = torque;
304 this.lambdaGrad = lambdaGrad;
305 this.lambdaTorque = lambdaTorque;
306 this.shareddEdLambda = shareddEdLambda;
307 this.sharedd2EdLambda2 = sharedd2EdLambda2;
308 }
309
310 @Override
311 public void run() throws Exception {
312 int threadIndex = getThreadIndex();
313 if (permanentReciprocalEnergyLoop[threadIndex] == null) {
314 permanentReciprocalEnergyLoop[threadIndex] = new PermanentReciprocalEnergyLoop();
315 inducedDipoleReciprocalEnergyLoop[threadIndex] = new InducedDipoleReciprocalEnergyLoop();
316 }
317 try {
318 int nAtoms = atoms.length;
319 execute(0, nAtoms - 1, permanentReciprocalEnergyLoop[threadIndex]);
320 if (polarization != Polarization.NONE) {
321 execute(0, nAtoms - 1, inducedDipoleReciprocalEnergyLoop[threadIndex]);
322 }
323 } catch (Exception e) {
324 String message =
325 "Fatal exception computing the real space field in thread " + threadIndex + "\n";
326 logger.log(Level.SEVERE, message, e);
327 }
328 }
329
330 @Override
331 public void start() {
332 multipole = globalMultipole[0];
333 fracMultipole = globalFracMultipole[0];
334 ind = inducedDipole[0];
335 indCR = inducedDipoleCR[0];
336 inducedDipoleSelfEnergy.set(0.0);
337 inducedDipoleRecipEnergy.set(0.0);
338 nfftX = reciprocalSpace.getXDim();
339 nfftY = reciprocalSpace.getYDim();
340 nfftZ = reciprocalSpace.getZDim();
341 }
342
343 private class PermanentReciprocalEnergyLoop extends IntegerForLoop {
344
345 int threadID;
346 double eSelf;
347 double eRecip;
348
349 @Override
350 public void finish() {
351 eSelf *= permanentScale;
352 eRecip *= permanentScale * 0.5 * electric;
353 }
354
355 @Override
356 public void run(int lb, int ub) throws Exception {
357
358
359 for (int i = lb; i <= ub; i++) {
360 if (use[i]) {
361 double[] in = globalMultipole[0][i];
362 double cii = in[t000] * in[t000];
363 double dii = in[t100] * in[t100] + in[t010] * in[t010] + in[t001] * in[t001];
364 double qii = in[t200] * in[t200] + in[t020] * in[t020] + in[t002] * in[t002]
365 + 2.0 * (in[t110] * in[t110] + in[t101] * in[t101] + in[t011] * in[t011]);
366 eSelf += aewald1 * (cii + aewald2 * (dii / 3.0 + 2.0 * aewald2 * qii / 45.0));
367
368 if (esvTerm && extendedSystem.isTitrating(i)) {
369 double[] indot = titrationMultipole[0][i];
370 double ciidot = indot[t000] * in[t000];
371 double diidot = indot[t100] * in[t100] + indot[t010] * in[t010] + indot[t001] * in[t001];
372 double qiidot = indot[t200] * in[t200] + indot[t020] * in[t020] + indot[t002] * in[t002]
373 + 2.0 * (indot[t110] * in[t110] + indot[t101] * in[t101] + indot[t011] * in[t011]);
374 double eSelfTitrDot =
375 2.0 * aewald1 * (ciidot + aewald2 * (diidot / 3.0 + 2.0 * aewald2 * qiidot / 45.0));
376 double eSelfTautDot = 0.0;
377 if (extendedSystem.isTautomerizing(i)) {
378 indot = tautomerMultipole[0][i];
379 ciidot = indot[t000] * in[t000];
380 diidot = indot[t100] * in[t100] + indot[t010] * in[t010] + indot[t001] * in[t001];
381 qiidot = indot[t200] * in[t200] + indot[t020] * in[t020] + indot[t002] * in[t002]
382 + 2.0 * (indot[t110] * in[t110] + indot[t101] * in[t101] + indot[t011] * in[t011]);
383 eSelfTautDot = 2.0 * aewald1 * (ciidot + aewald2 * (diidot / 3.0
384 + 2.0 * aewald2 * qiidot / 45.0));
385 }
386 double factor = permanentScale;
387 extendedSystem.addPermElecDeriv(i, eSelfTitrDot * factor, eSelfTautDot * factor);
388 }
389 }
390 }
391 if (lambdaTerm) {
392 shareddEdLambda.addAndGet(eSelf * dlPowPerm * dEdLSign);
393 sharedd2EdLambda2.addAndGet(eSelf * d2lPowPerm * dEdLSign);
394 }
395
396
397 final double[][] recip = crystal.getUnitCell().A;
398
399 double dUdL = 0.0;
400 double d2UdL2 = 0.0;
401 for (int i = lb; i <= ub; i++) {
402 if (use[i]) {
403 final double[] phi = cartMultipolePhi[i];
404 final double[] mpole = multipole[i];
405 final double[] fmpole = fracMultipole[i];
406 double e = mpole[t000] * phi[t000]
407 + mpole[t100] * phi[t100] + mpole[t010] * phi[t010] + mpole[t001] * phi[t001]
408 + oneThird * (mpole[t200] * phi[t200] + mpole[t020] * phi[t020]
409 + mpole[t002] * phi[t002]
410 + 2.0 * (mpole[t110] * phi[t110] + mpole[t101] * phi[t101] + mpole[t011] * phi[t011]));
411 eRecip += e;
412 if (esvTerm && extendedSystem.isTitrating(i)) {
413 double[] mpoleDot = titrationMultipole[0][i];
414 double edotTitr = mpoleDot[t000] * phi[t000]
415 + mpoleDot[t100] * phi[t100] + mpoleDot[t010] * phi[t010]
416 + mpoleDot[t001] * phi[t001]
417 + oneThird * (mpoleDot[t200] * phi[t200] + mpoleDot[t020] * phi[t020]
418 + mpoleDot[t002] * phi[t002]
419 + 2.0 * (mpoleDot[t110] * phi[t110] + mpoleDot[t101] * phi[t101]
420 + mpoleDot[t011] * phi[t011]));
421 double edotTaut = 0.0;
422 if (extendedSystem.isTautomerizing(i)) {
423 mpoleDot = tautomerMultipole[0][i];
424 edotTaut = mpoleDot[t000] * phi[t000]
425 + mpoleDot[t100] * phi[t100] + mpoleDot[t010] * phi[t010]
426 + mpoleDot[t001] * phi[t001]
427 + oneThird * (mpoleDot[t200] * phi[t200] + mpoleDot[t020] * phi[t020]
428 + mpoleDot[t002] * phi[t002]
429 + 2.0 * (mpoleDot[t110] * phi[t110] + mpoleDot[t101] * phi[t101]
430 + mpoleDot[t011] * phi[t011]));
431 }
432 double factor = permanentScale * electric;
433 extendedSystem.addPermElecDeriv(i, edotTitr * factor, edotTaut * factor);
434 }
435 if (gradient || lambdaTerm) {
436 final double[] fPhi = fracMultipolePhi[i];
437 double gx = fmpole[t000] * fPhi[t100]
438 + fmpole[t100] * fPhi[t200] + fmpole[t010] * fPhi[t110] + fmpole[t001] * fPhi[t101]
439 + fmpole[t200] * fPhi[t300] + fmpole[t020] * fPhi[t120] + fmpole[t002] * fPhi[t102]
440 + fmpole[t110] * fPhi[t210] + fmpole[t101] * fPhi[t201] + fmpole[t011] * fPhi[t111];
441 double gy = fmpole[t000] * fPhi[t010]
442 + fmpole[t100] * fPhi[t110] + fmpole[t010] * fPhi[t020] + fmpole[t001] * fPhi[t011]
443 + fmpole[t200] * fPhi[t210] + fmpole[t020] * fPhi[t030] + fmpole[t002] * fPhi[t012]
444 + fmpole[t110] * fPhi[t120] + fmpole[t101] * fPhi[t111] + fmpole[t011] * fPhi[t021];
445 double gz = fmpole[t000] * fPhi[t001]
446 + fmpole[t100] * fPhi[t101] + fmpole[t010] * fPhi[t011] + fmpole[t001] * fPhi[t002]
447 + fmpole[t200] * fPhi[t201] + fmpole[t020] * fPhi[t021] + fmpole[t002] * fPhi[t003]
448 + fmpole[t110] * fPhi[t111] + fmpole[t101] * fPhi[t102] + fmpole[t011] * fPhi[t012];
449 gx *= nfftX;
450 gy *= nfftY;
451 gz *= nfftZ;
452 final double dfx = recip[0][0] * gx + recip[0][1] * gy + recip[0][2] * gz;
453 final double dfy = recip[1][0] * gx + recip[1][1] * gy + recip[1][2] * gz;
454 final double dfz = recip[2][0] * gx + recip[2][1] * gy + recip[2][2] * gz;
455
456 double tqx = -mpole[t010] * phi[t001] + mpole[t001] * phi[t010];
457 double tqy = -mpole[t001] * phi[t100] + mpole[t100] * phi[t001];
458 double tqz = -mpole[t100] * phi[t010] + mpole[t010] * phi[t100];
459
460 tqx -= twoThirds * (
461 mpole[t110] * phi[t101] + mpole[t020] * phi[t011] + mpole[t011] * phi[t002]
462 - mpole[t101] * phi[t110] - mpole[t011] * phi[t020] - mpole[t002] * phi[t011]);
463 tqy -= twoThirds * (
464 mpole[t101] * phi[t200] + mpole[t011] * phi[t110] + mpole[t002] * phi[t101]
465 - mpole[t200] * phi[t101] - mpole[t110] * phi[t011] - mpole[t101] * phi[t002]);
466 tqz -= twoThirds * (
467 mpole[t200] * phi[t110] + mpole[t110] * phi[t020] + mpole[t101] * phi[t011]
468 - mpole[t110] * phi[t200] - mpole[t020] * phi[t110] - mpole[t011] * phi[t101]);
469 if (gradient) {
470 double factor = permanentScale * electric;
471 grad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
472 torque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
473 }
474 if (lambdaTerm) {
475 dUdL += dEdLSign * dlPowPerm * e;
476 d2UdL2 += dEdLSign * d2lPowPerm * e;
477 double factor = dEdLSign * dlPowPerm * electric;
478 lambdaGrad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
479 lambdaTorque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
480 }
481 }
482 }
483 }
484
485 if (lambdaTerm) {
486 shareddEdLambda.addAndGet(0.5 * dUdL * electric);
487 sharedd2EdLambda2.addAndGet(0.5 * d2UdL2 * electric);
488 }
489 }
490
491 @Override
492 public IntegerSchedule schedule() {
493 return IntegerSchedule.fixed();
494 }
495
496 @Override
497 public void start() {
498 eSelf = 0.0;
499 eRecip = 0.0;
500 threadID = getThreadIndex();
501 }
502 }
503
504 private class InducedDipoleReciprocalEnergyLoop extends IntegerForLoop {
505
506 private final double[] sfPhi = new double[tensorCount];
507 private final double[] sPhi = new double[tensorCount];
508 private double eSelf;
509 private double eRecip;
510 private int threadID;
511
512 @Override
513 public void finish() {
514 inducedDipoleSelfEnergy.addAndGet(polarizationScale * eSelf);
515 inducedDipoleRecipEnergy.addAndGet(polarizationScale * eRecip);
516 }
517
518 @Override
519 public void run(int lb, int ub) throws Exception {
520
521 for (int i = lb; i <= ub; i++) {
522 if (use[i]) {
523 final double[] indi = ind[i];
524 final double[] multipolei = multipole[i];
525 final double dix = multipolei[t100];
526 final double diy = multipolei[t010];
527 final double diz = multipolei[t001];
528 final double dii = indi[0] * dix + indi[1] * diy + indi[2] * diz;
529 eSelf += aewald3 * dii;
530 }
531 }
532 if (lambdaTerm) {
533 shareddEdLambda.addAndGet(dEdLSign * dlPowPol * eSelf);
534 sharedd2EdLambda2.addAndGet(dEdLSign * d2lPowPol * eSelf);
535 }
536 if (gradient || esvTerm) {
537 for (int i = lb; i <= ub; i++) {
538 if (use[i]) {
539 final double[] indi = ind[i];
540 final double[] indpi = indCR[i];
541 final double[] multipolei = multipole[i];
542 final double dix = multipolei[t100];
543 final double diy = multipolei[t010];
544 final double diz = multipolei[t001];
545 final double uix = 0.5 * (indi[0] + indpi[0]);
546 final double uiy = 0.5 * (indi[1] + indpi[1]);
547 final double uiz = 0.5 * (indi[2] + indpi[2]);
548 if (gradient) {
549 final double tix = aewald4 * (diy * uiz - diz * uiy);
550 final double tiy = aewald4 * (diz * uix - dix * uiz);
551 final double tiz = aewald4 * (dix * uiy - diy * uix);
552 torque.add(threadID, i, polarizationScale * tix, polarizationScale * tiy, polarizationScale * tiz);
553 if (lambdaTerm) {
554 double factor = dEdLSign * dlPowPol;
555 lambdaTorque.add(threadID, i, factor * tix, factor * tiy, factor * tiz);
556 }
557 }
558 if (esvTerm && extendedSystem.isTitrating(i)) {
559 double[] titrDot = titrationMultipole[0][i];
560 double indiDotMdot = uix * titrDot[t100] + uiy * titrDot[t010] + uiz * titrDot[t001];
561 double dIndSelfTitrdLi = polarizationScale * 2.0 * aewald3 * indiDotMdot;
562 double dIndSelfTautdLi = 0.0;
563 if (extendedSystem.isTautomerizing(i)) {
564 double[] tautDot = tautomerMultipole[0][i];
565 indiDotMdot = uix * tautDot[t100] + uiy * tautDot[t010] + uiz * tautDot[t001];
566 dIndSelfTautdLi = polarizationScale * 2.0 * aewald3 * indiDotMdot;
567 }
568 extendedSystem.addIndElecDeriv(i, dIndSelfTitrdLi, dIndSelfTautdLi);
569 }
570 }
571 }
572 }
573
574
575 for (int i = lb; i <= ub; i++) {
576 double[] fracInd = new double[3];
577 double[] fracIndCR = new double[3];
578 if (use[i]) {
579 final double[] fPhi = fracMultipolePhi[i];
580 reciprocalSpace.cartToFracInducedDipole(inducedDipole[0][i], inducedDipoleCR[0][i], fracInd, fracIndCR);
581 final double indx = fracInd[0];
582 final double indy = fracInd[1];
583 final double indz = fracInd[2];
584 eRecip += indx * fPhi[t100] + indy * fPhi[t010] + indz * fPhi[t001];
585 if (gradient || esvTerm) {
586 final double[] iPhi = cartesianDipolePhi[i];
587 final double[] iCRPhi = cartesianDipolePhiCR[i];
588 final double[] fiPhi = fracInducedDipolePhi[i];
589 final double[] fiCRPhi = fracInducedDipolePhiCR[i];
590 final double[] mpolei = multipole[i];
591 final double[] fmpolei = fracMultipole[i];
592 final double inpx = fracIndCR[0];
593 final double inpy = fracIndCR[1];
594 final double inpz = fracIndCR[2];
595 final double insx = indx + inpx;
596 final double insy = indy + inpy;
597 final double insz = indz + inpz;
598 for (int t = 0; t < tensorCount; t++) {
599 sPhi[t] = 0.5 * (iPhi[t] + iCRPhi[t]);
600 sfPhi[t] = fiPhi[t] + fiCRPhi[t];
601 }
602 double gx = insx * fPhi[t200] + insy * fPhi[t110] + insz * fPhi[t101];
603 double gy = insx * fPhi[t110] + insy * fPhi[t020] + insz * fPhi[t011];
604 double gz = insx * fPhi[t101] + insy * fPhi[t011] + insz * fPhi[t002];
605 if (polarization == Polarization.MUTUAL) {
606 gx += indx * fiCRPhi[t200] + inpx * fiPhi[t200] + indy * fiCRPhi[t110]
607 + inpy * fiPhi[t110] + indz * fiCRPhi[t101] + inpz * fiPhi[t101];
608 gy += indx * fiCRPhi[t110] + inpx * fiPhi[t110] + indy * fiCRPhi[t020]
609 + inpy * fiPhi[t020] + indz * fiCRPhi[t011] + inpz * fiPhi[t011];
610 gz += indx * fiCRPhi[t101] + inpx * fiPhi[t101] + indy * fiCRPhi[t011]
611 + inpy * fiPhi[t011] + indz * fiCRPhi[t002] + inpz * fiPhi[t002];
612 }
613 gx += fmpolei[t000] * sfPhi[t100] + fmpolei[t100] * sfPhi[t200]
614 + fmpolei[t010] * sfPhi[t110] + fmpolei[t001] * sfPhi[t101]
615 + fmpolei[t200] * sfPhi[t300] + fmpolei[t020] * sfPhi[t120]
616 + fmpolei[t002] * sfPhi[t102] + fmpolei[t110] * sfPhi[t210]
617 + fmpolei[t101] * sfPhi[t201] + fmpolei[t011] * sfPhi[t111];
618 gy += fmpolei[t000] * sfPhi[t010] + fmpolei[t100] * sfPhi[t110]
619 + fmpolei[t010] * sfPhi[t020] + fmpolei[t001] * sfPhi[t011]
620 + fmpolei[t200] * sfPhi[t210] + fmpolei[t020] * sfPhi[t030]
621 + fmpolei[t002] * sfPhi[t012] + fmpolei[t110] * sfPhi[t120]
622 + fmpolei[t101] * sfPhi[t111] + fmpolei[t011] * sfPhi[t021];
623 gz += fmpolei[t000] * sfPhi[t001] + fmpolei[t100] * sfPhi[t101]
624 + fmpolei[t010] * sfPhi[t011] + fmpolei[t001] * sfPhi[t002]
625 + fmpolei[t200] * sfPhi[t201] + fmpolei[t020] * sfPhi[t021]
626 + fmpolei[t002] * sfPhi[t003] + fmpolei[t110] * sfPhi[t111]
627 + fmpolei[t101] * sfPhi[t102] + fmpolei[t011] * sfPhi[t012];
628 gx *= nfftX;
629 gy *= nfftY;
630 gz *= nfftZ;
631 double[][] recip = crystal.getUnitCell().A;
632 double dfx = recip[0][0] * gx + recip[0][1] * gy + recip[0][2] * gz;
633 double dfy = recip[1][0] * gx + recip[1][1] * gy + recip[1][2] * gz;
634 double dfz = recip[2][0] * gx + recip[2][1] * gy + recip[2][2] * gz;
635 dfx *= 0.5 * electric;
636 dfy *= 0.5 * electric;
637 dfz *= 0.5 * electric;
638
639 double tqx = -mpolei[t010] * sPhi[t001] + mpolei[t001] * sPhi[t010];
640 double tqy = -mpolei[t001] * sPhi[t100] + mpolei[t100] * sPhi[t001];
641 double tqz = -mpolei[t100] * sPhi[t010] + mpolei[t010] * sPhi[t100];
642
643 tqx -= twoThirds * (
644 mpolei[t110] * sPhi[t101] + mpolei[t020] * sPhi[t011] + mpolei[t011] * sPhi[t002]
645 - mpolei[t101] * sPhi[t110] - mpolei[t011] * sPhi[t020]
646 - mpolei[t002] * sPhi[t011]);
647 tqy -= twoThirds * (
648 mpolei[t101] * sPhi[t200] + mpolei[t011] * sPhi[t110] + mpolei[t002] * sPhi[t101]
649 - mpolei[t200] * sPhi[t101] - mpolei[t110] * sPhi[t011]
650 - mpolei[t101] * sPhi[t002]);
651 tqz -= twoThirds * (
652 mpolei[t200] * sPhi[t110] + mpolei[t110] * sPhi[t020] + mpolei[t101] * sPhi[t011]
653 - mpolei[t110] * sPhi[t200] - mpolei[t020] * sPhi[t110]
654 - mpolei[t011] * sPhi[t101]);
655 tqx *= electric;
656 tqy *= electric;
657 tqz *= electric;
658 grad.add(threadID, i, polarizationScale * dfx, polarizationScale * dfy,
659 polarizationScale * dfz);
660 torque.add(threadID, i, polarizationScale * tqx, polarizationScale * tqy,
661 polarizationScale * tqz);
662 if (lambdaTerm) {
663 double factor = dEdLSign * dlPowPol;
664 lambdaGrad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
665 lambdaTorque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
666 }
667 if (esvTerm && extendedSystem.isTitrating(i)) {
668 final double[] mTitrDot = titrationMultipole[0][i];
669 double eTitrDot = mTitrDot[t000] * sPhi[t000]
670 + mTitrDot[t100] * sPhi[t100] + mTitrDot[t010] * sPhi[t010] + mTitrDot[t001] * sPhi[t001]
671 + oneThird * (mTitrDot[t200] * sPhi[t200] + mTitrDot[t020] * sPhi[t020] + mTitrDot[t002] * sPhi[t002]
672 + 2.0 * (mTitrDot[t110] * sPhi[t110] + mTitrDot[t101] * sPhi[t101] + mTitrDot[t011] * sPhi[t011]));
673 double eTautDot = 0.0;
674 if (extendedSystem.isTautomerizing(i)) {
675 final double[] mTautDot = tautomerMultipole[0][i];
676 eTautDot = mTautDot[t000] * sPhi[t000]
677 + mTautDot[t100] * sPhi[t100] + mTautDot[t010] * sPhi[t010] + mTautDot[t001] * sPhi[t001]
678 + oneThird * (mTautDot[t200] * sPhi[t200] + mTautDot[t020] * sPhi[t020] + mTautDot[t002] * sPhi[t002]
679 + 2.0 * (mTautDot[t110] * sPhi[t110] + mTautDot[t101] * sPhi[t101] + mTautDot[t011] * sPhi[t011]));
680 }
681 double factor = polarizationScale * electric;
682 extendedSystem.addIndElecDeriv(i, eTitrDot * factor, eTautDot * factor);
683 }
684 }
685 }
686 }
687 eRecip *= 0.5 * electric;
688 if (lambdaTerm) {
689 shareddEdLambda.addAndGet(dEdLSign * dlPowPol * eRecip);
690 sharedd2EdLambda2.addAndGet(dEdLSign * d2lPowPol * eRecip);
691 }
692 }
693
694 @Override
695 public IntegerSchedule schedule() {
696 return IntegerSchedule.fixed();
697 }
698
699 @Override
700 public void start() {
701 eSelf = 0.0;
702 eRecip = 0.0;
703 threadID = getThreadIndex();
704 }
705 }
706 }