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