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 || (esvTerm && !extendedSystem.useTotalChargeCorrection)) {
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 if (esvTerm){
245 for(int i=0; i<atoms.length;i++){
246 if(extendedSystem.isTitrating(i)){
247 double[] in = globalMultipole[0][i];
248 double[] indot = titrationMultipole[0][i];
249 double chargeCorrectiondUdL = -0.5 * electric * PI * indot[t000] * (2*totalCharge) / denom;
250 extendedSystem.addPermElecDeriv(i,chargeCorrectiondUdL,0.0);
251 }
252 }
253 }
254 chargeCorrectionEnergy = chargeCorrectionEnergy * permanentScale;
255 }
256 }
257
258 public double getInducedDipoleReciprocalEnergy() {
259 return inducedDipoleRecipEnergy.get();
260 }
261
262 public double getInducedDipoleSelfEnergy() {
263 return inducedDipoleSelfEnergy.get();
264 }
265
266 public double getPermanentReciprocalEnergy() {
267 return permanentReciprocalEnergy;
268 }
269
270 public double getPermanentSelfEnergy() {
271 return permanentSelfEnergy;
272 }
273
274 public double getPermanentChargeCorrectionEnergy() {
275 return chargeCorrectionEnergy;
276 }
277
278 public void init(
279 Atom[] atoms,
280 Crystal crystal,
281 boolean gradient,
282 boolean lambdaTerm,
283 boolean esvTerm,
284 boolean[] use,
285 double[][][] globalMultipole,
286 double[][][] globalFracMultipole,
287 double[][][] titrationMultipole,
288 double[][][] tautomerMultipole,
289 double[][] cartMultipolePhi,
290 double[][] fracMultipolePhi,
291 Polarization polarization,
292 double[][][] inducedDipole,
293 double[][][] inducedDipoleCR,
294 double[][] cartesianDipolePhi,
295 double[][] cartesianDipolePhiCR,
296 double[][] fracInducedDipolePhi,
297 double[][] fracInducedDipolePhiCR,
298 ReciprocalSpace reciprocalSpace,
299 AlchemicalParameters alchemicalParameters,
300 ExtendedSystem extendedSystem,
301 AtomicDoubleArray3D grad,
302 AtomicDoubleArray3D torque,
303 AtomicDoubleArray3D lambdaGrad,
304 AtomicDoubleArray3D lambdaTorque,
305 SharedDouble shareddEdLambda,
306 SharedDouble sharedd2EdLambda2
307 ) {
308 this.atoms = atoms;
309 this.crystal = crystal;
310 this.gradient = gradient;
311 this.lambdaTerm = lambdaTerm;
312 this.esvTerm = esvTerm;
313 this.use = use;
314
315 this.globalMultipole = globalMultipole;
316 this.globalFracMultipole = globalFracMultipole;
317 this.titrationMultipole = titrationMultipole;
318 this.tautomerMultipole = tautomerMultipole;
319 this.cartMultipolePhi = cartMultipolePhi;
320 this.fracMultipolePhi = fracMultipolePhi;
321
322 this.polarization = polarization;
323 this.inducedDipole = inducedDipole;
324 this.inducedDipoleCR = inducedDipoleCR;
325 this.cartesianDipolePhi = cartesianDipolePhi;
326 this.cartesianDipolePhiCR = cartesianDipolePhiCR;
327 this.fracInducedDipolePhi = fracInducedDipolePhi;
328 this.fracInducedDipolePhiCR = fracInducedDipolePhiCR;
329 this.reciprocalSpace = reciprocalSpace;
330 this.extendedSystem = extendedSystem;
331
332 this.permanentScale = alchemicalParameters.permanentScale;
333 this.dlPowPerm = alchemicalParameters.dlPowPerm;
334 this.d2lPowPerm = alchemicalParameters.d2lPowPerm;
335 this.polarizationScale = alchemicalParameters.polarizationScale;
336 this.dlPowPol = alchemicalParameters.dlPowPol;
337 this.d2lPowPol = alchemicalParameters.d2lPowPol;
338 this.dEdLSign = alchemicalParameters.dEdLSign;
339
340 this.grad = grad;
341 this.torque = torque;
342 this.lambdaGrad = lambdaGrad;
343 this.lambdaTorque = lambdaTorque;
344 this.shareddEdLambda = shareddEdLambda;
345 this.sharedd2EdLambda2 = sharedd2EdLambda2;
346 }
347
348 @Override
349 public void run() throws Exception {
350 int threadIndex = getThreadIndex();
351 if (permanentReciprocalEnergyLoop[threadIndex] == null) {
352 permanentReciprocalEnergyLoop[threadIndex] = new PermanentReciprocalEnergyLoop();
353 inducedDipoleReciprocalEnergyLoop[threadIndex] = new InducedDipoleReciprocalEnergyLoop();
354 }
355 try {
356 int nAtoms = atoms.length;
357 execute(0, nAtoms - 1, permanentReciprocalEnergyLoop[threadIndex]);
358 if (polarization != Polarization.NONE) {
359 execute(0, nAtoms - 1, inducedDipoleReciprocalEnergyLoop[threadIndex]);
360 }
361 } catch (Exception e) {
362 String message =
363 "Fatal exception computing the real space field in thread " + threadIndex + "\n";
364 logger.log(Level.SEVERE, message, e);
365 }
366 }
367
368 @Override
369 public void start() {
370 multipole = globalMultipole[0];
371 fracMultipole = globalFracMultipole[0];
372 ind = inducedDipole[0];
373 indCR = inducedDipoleCR[0];
374 inducedDipoleSelfEnergy.set(0.0);
375 inducedDipoleRecipEnergy.set(0.0);
376 nfftX = reciprocalSpace.getXDim();
377 nfftY = reciprocalSpace.getYDim();
378 nfftZ = reciprocalSpace.getZDim();
379 }
380
381 private class PermanentReciprocalEnergyLoop extends IntegerForLoop {
382
383 int threadID;
384 double totalCharge;
385 double eSelf;
386 double eRecip;
387
388 @Override
389 public void finish() {
390 eSelf *= permanentScale;
391 eRecip *= permanentScale * 0.5 * electric;
392 }
393
394 @Override
395 public void run(int lb, int ub) {
396
397
398 for (int i = lb; i <= ub; i++) {
399 if (use[i]) {
400 double[] in = globalMultipole[0][i];
401 totalCharge += in[t000];
402 double cii = in[t000] * in[t000];
403 double dii = in[t100] * in[t100] + in[t010] * in[t010] + in[t001] * in[t001];
404 double qii = in[t200] * in[t200] + in[t020] * in[t020] + in[t002] * in[t002]
405 + 2.0 * (in[t110] * in[t110] + in[t101] * in[t101] + in[t011] * in[t011]);
406 eSelf += aewald1 * (cii + aewald2 * (dii / 3.0 + 2.0 * aewald2 * qii / 45.0));
407
408 if (esvTerm && extendedSystem.isTitrating(i)) {
409 double[] indot = titrationMultipole[0][i];
410 double ciidot = indot[t000] * in[t000];
411 double diidot = indot[t100] * in[t100] + indot[t010] * in[t010] + indot[t001] * in[t001];
412 double qiidot = indot[t200] * in[t200] + indot[t020] * in[t020] + indot[t002] * in[t002]
413 + 2.0 * (indot[t110] * in[t110] + indot[t101] * in[t101] + indot[t011] * in[t011]);
414 double eSelfTitrDot =
415 2.0 * aewald1 * (ciidot + aewald2 * (diidot / 3.0 + 2.0 * aewald2 * qiidot / 45.0));
416 double eSelfTautDot = 0.0;
417 if (extendedSystem.isTautomerizing(i)) {
418 indot = tautomerMultipole[0][i];
419 ciidot = indot[t000] * in[t000];
420 diidot = indot[t100] * in[t100] + indot[t010] * in[t010] + indot[t001] * in[t001];
421 qiidot = indot[t200] * in[t200] + indot[t020] * in[t020] + indot[t002] * in[t002]
422 + 2.0 * (indot[t110] * in[t110] + indot[t101] * in[t101] + indot[t011] * in[t011]);
423 eSelfTautDot = 2.0 * aewald1 * (ciidot + aewald2 * (diidot / 3.0
424 + 2.0 * aewald2 * qiidot / 45.0));
425 }
426 double factor = permanentScale;
427 extendedSystem.addPermElecDeriv(i, eSelfTitrDot * factor, eSelfTautDot * factor);
428 }
429 }
430 }
431 if (lambdaTerm) {
432 shareddEdLambda.addAndGet(eSelf * dlPowPerm * dEdLSign);
433 sharedd2EdLambda2.addAndGet(eSelf * d2lPowPerm * dEdLSign);
434 }
435
436
437 final double[][] recip = crystal.getUnitCell().A;
438
439 double dUdL = 0.0;
440 double d2UdL2 = 0.0;
441 for (int i = lb; i <= ub; i++) {
442 if (use[i]) {
443 final double[] phi = cartMultipolePhi[i];
444 final double[] mpole = multipole[i];
445 final double[] fmpole = fracMultipole[i];
446 double e = mpole[t000] * phi[t000]
447 + mpole[t100] * phi[t100] + mpole[t010] * phi[t010] + mpole[t001] * phi[t001]
448 + oneThird * (mpole[t200] * phi[t200] + mpole[t020] * phi[t020]
449 + mpole[t002] * phi[t002]
450 + 2.0 * (mpole[t110] * phi[t110] + mpole[t101] * phi[t101] + mpole[t011] * phi[t011]));
451 eRecip += e;
452 if (esvTerm && extendedSystem.isTitrating(i)) {
453 double[] mpoleDot = titrationMultipole[0][i];
454 double edotTitr = 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 double edotTaut = 0.0;
462 if (extendedSystem.isTautomerizing(i)) {
463 mpoleDot = tautomerMultipole[0][i];
464 edotTaut = mpoleDot[t000] * phi[t000]
465 + mpoleDot[t100] * phi[t100] + mpoleDot[t010] * phi[t010]
466 + mpoleDot[t001] * phi[t001]
467 + oneThird * (mpoleDot[t200] * phi[t200] + mpoleDot[t020] * phi[t020]
468 + mpoleDot[t002] * phi[t002]
469 + 2.0 * (mpoleDot[t110] * phi[t110] + mpoleDot[t101] * phi[t101]
470 + mpoleDot[t011] * phi[t011]));
471 }
472 double factor = permanentScale * electric;
473 extendedSystem.addPermElecDeriv(i, edotTitr * factor, edotTaut * factor);
474 }
475 if (gradient || lambdaTerm) {
476 final double[] fPhi = fracMultipolePhi[i];
477 double gx = fmpole[t000] * fPhi[t100]
478 + fmpole[t100] * fPhi[t200] + fmpole[t010] * fPhi[t110] + fmpole[t001] * fPhi[t101]
479 + fmpole[t200] * fPhi[t300] + fmpole[t020] * fPhi[t120] + fmpole[t002] * fPhi[t102]
480 + fmpole[t110] * fPhi[t210] + fmpole[t101] * fPhi[t201] + fmpole[t011] * fPhi[t111];
481 double gy = fmpole[t000] * fPhi[t010]
482 + fmpole[t100] * fPhi[t110] + fmpole[t010] * fPhi[t020] + fmpole[t001] * fPhi[t011]
483 + fmpole[t200] * fPhi[t210] + fmpole[t020] * fPhi[t030] + fmpole[t002] * fPhi[t012]
484 + fmpole[t110] * fPhi[t120] + fmpole[t101] * fPhi[t111] + fmpole[t011] * fPhi[t021];
485 double gz = fmpole[t000] * fPhi[t001]
486 + fmpole[t100] * fPhi[t101] + fmpole[t010] * fPhi[t011] + fmpole[t001] * fPhi[t002]
487 + fmpole[t200] * fPhi[t201] + fmpole[t020] * fPhi[t021] + fmpole[t002] * fPhi[t003]
488 + fmpole[t110] * fPhi[t111] + fmpole[t101] * fPhi[t102] + fmpole[t011] * fPhi[t012];
489 gx *= nfftX;
490 gy *= nfftY;
491 gz *= nfftZ;
492 final double dfx = recip[0][0] * gx + recip[0][1] * gy + recip[0][2] * gz;
493 final double dfy = recip[1][0] * gx + recip[1][1] * gy + recip[1][2] * gz;
494 final double dfz = recip[2][0] * gx + recip[2][1] * gy + recip[2][2] * gz;
495
496 double tqx = -mpole[t010] * phi[t001] + mpole[t001] * phi[t010];
497 double tqy = -mpole[t001] * phi[t100] + mpole[t100] * phi[t001];
498 double tqz = -mpole[t100] * phi[t010] + mpole[t010] * phi[t100];
499
500 tqx -= twoThirds * (
501 mpole[t110] * phi[t101] + mpole[t020] * phi[t011] + mpole[t011] * phi[t002]
502 - mpole[t101] * phi[t110] - mpole[t011] * phi[t020] - mpole[t002] * phi[t011]);
503 tqy -= twoThirds * (
504 mpole[t101] * phi[t200] + mpole[t011] * phi[t110] + mpole[t002] * phi[t101]
505 - mpole[t200] * phi[t101] - mpole[t110] * phi[t011] - mpole[t101] * phi[t002]);
506 tqz -= twoThirds * (
507 mpole[t200] * phi[t110] + mpole[t110] * phi[t020] + mpole[t101] * phi[t011]
508 - mpole[t110] * phi[t200] - mpole[t020] * phi[t110] - mpole[t011] * phi[t101]);
509 if (gradient) {
510 double factor = permanentScale * electric;
511 grad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
512 torque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
513 }
514 if (lambdaTerm) {
515 dUdL += dEdLSign * dlPowPerm * e;
516 d2UdL2 += dEdLSign * d2lPowPerm * e;
517 double factor = dEdLSign * dlPowPerm * electric;
518 lambdaGrad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
519 lambdaTorque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
520 }
521 }
522 }
523 }
524
525 if (lambdaTerm) {
526 shareddEdLambda.addAndGet(0.5 * dUdL * electric);
527 sharedd2EdLambda2.addAndGet(0.5 * d2UdL2 * electric);
528 }
529 }
530
531 @Override
532 public IntegerSchedule schedule() {
533 return IntegerSchedule.fixed();
534 }
535
536 @Override
537 public void start() {
538 totalCharge = 0.0;
539 eSelf = 0.0;
540 eRecip = 0.0;
541 threadID = getThreadIndex();
542 }
543 }
544
545 private class InducedDipoleReciprocalEnergyLoop extends IntegerForLoop {
546
547 private final double[] sfPhi = new double[tensorCount];
548 private final double[] sPhi = new double[tensorCount];
549 private double eSelf;
550 private double eRecip;
551 private int threadID;
552
553 @Override
554 public void finish() {
555 inducedDipoleSelfEnergy.addAndGet(polarizationScale * eSelf);
556 inducedDipoleRecipEnergy.addAndGet(polarizationScale * eRecip);
557 }
558
559 @Override
560 public void run(int lb, int ub) throws Exception {
561
562 for (int i = lb; i <= ub; i++) {
563 if (use[i]) {
564 final double[] indi = ind[i];
565 final double[] multipolei = multipole[i];
566 final double dix = multipolei[t100];
567 final double diy = multipolei[t010];
568 final double diz = multipolei[t001];
569 final double dii = indi[0] * dix + indi[1] * diy + indi[2] * diz;
570 eSelf += aewald3 * dii;
571 }
572 }
573 if (lambdaTerm) {
574 shareddEdLambda.addAndGet(dEdLSign * dlPowPol * eSelf);
575 sharedd2EdLambda2.addAndGet(dEdLSign * d2lPowPol * eSelf);
576 }
577 if (gradient || esvTerm) {
578 for (int i = lb; i <= ub; i++) {
579 if (use[i]) {
580 final double[] indi = ind[i];
581 final double[] indpi = indCR[i];
582 final double[] multipolei = multipole[i];
583 final double dix = multipolei[t100];
584 final double diy = multipolei[t010];
585 final double diz = multipolei[t001];
586 final double uix = 0.5 * (indi[0] + indpi[0]);
587 final double uiy = 0.5 * (indi[1] + indpi[1]);
588 final double uiz = 0.5 * (indi[2] + indpi[2]);
589 if (gradient) {
590 final double tix = aewald4 * (diy * uiz - diz * uiy);
591 final double tiy = aewald4 * (diz * uix - dix * uiz);
592 final double tiz = aewald4 * (dix * uiy - diy * uix);
593 torque.add(threadID, i, polarizationScale * tix, polarizationScale * tiy, polarizationScale * tiz);
594 if (lambdaTerm) {
595 double factor = dEdLSign * dlPowPol;
596 lambdaTorque.add(threadID, i, factor * tix, factor * tiy, factor * tiz);
597 }
598 }
599 if (esvTerm && extendedSystem.isTitrating(i)) {
600 double[] titrDot = titrationMultipole[0][i];
601 double indiDotMdot = uix * titrDot[t100] + uiy * titrDot[t010] + uiz * titrDot[t001];
602 double dIndSelfTitrdLi = polarizationScale * 2.0 * aewald3 * indiDotMdot;
603 double dIndSelfTautdLi = 0.0;
604 if (extendedSystem.isTautomerizing(i)) {
605 double[] tautDot = tautomerMultipole[0][i];
606 indiDotMdot = uix * tautDot[t100] + uiy * tautDot[t010] + uiz * tautDot[t001];
607 dIndSelfTautdLi = polarizationScale * 2.0 * aewald3 * indiDotMdot;
608 }
609 extendedSystem.addIndElecDeriv(i, dIndSelfTitrdLi, dIndSelfTautdLi);
610 }
611 }
612 }
613 }
614
615
616 for (int i = lb; i <= ub; i++) {
617 double[] fracInd = new double[3];
618 double[] fracIndCR = new double[3];
619 if (use[i]) {
620 final double[] fPhi = fracMultipolePhi[i];
621 reciprocalSpace.cartToFracInducedDipole(inducedDipole[0][i], inducedDipoleCR[0][i], fracInd, fracIndCR);
622 final double indx = fracInd[0];
623 final double indy = fracInd[1];
624 final double indz = fracInd[2];
625 eRecip += indx * fPhi[t100] + indy * fPhi[t010] + indz * fPhi[t001];
626 if (gradient || esvTerm) {
627 final double[] iPhi = cartesianDipolePhi[i];
628 final double[] iCRPhi = cartesianDipolePhiCR[i];
629 final double[] fiPhi = fracInducedDipolePhi[i];
630 final double[] fiCRPhi = fracInducedDipolePhiCR[i];
631 final double[] mpolei = multipole[i];
632 final double[] fmpolei = fracMultipole[i];
633 final double inpx = fracIndCR[0];
634 final double inpy = fracIndCR[1];
635 final double inpz = fracIndCR[2];
636 final double insx = indx + inpx;
637 final double insy = indy + inpy;
638 final double insz = indz + inpz;
639 for (int t = 0; t < tensorCount; t++) {
640 sPhi[t] = 0.5 * (iPhi[t] + iCRPhi[t]);
641 sfPhi[t] = fiPhi[t] + fiCRPhi[t];
642 }
643 double gx = insx * fPhi[t200] + insy * fPhi[t110] + insz * fPhi[t101];
644 double gy = insx * fPhi[t110] + insy * fPhi[t020] + insz * fPhi[t011];
645 double gz = insx * fPhi[t101] + insy * fPhi[t011] + insz * fPhi[t002];
646 if (polarization == Polarization.MUTUAL) {
647 gx += indx * fiCRPhi[t200] + inpx * fiPhi[t200] + indy * fiCRPhi[t110]
648 + inpy * fiPhi[t110] + indz * fiCRPhi[t101] + inpz * fiPhi[t101];
649 gy += indx * fiCRPhi[t110] + inpx * fiPhi[t110] + indy * fiCRPhi[t020]
650 + inpy * fiPhi[t020] + indz * fiCRPhi[t011] + inpz * fiPhi[t011];
651 gz += indx * fiCRPhi[t101] + inpx * fiPhi[t101] + indy * fiCRPhi[t011]
652 + inpy * fiPhi[t011] + indz * fiCRPhi[t002] + inpz * fiPhi[t002];
653 }
654 gx += fmpolei[t000] * sfPhi[t100] + fmpolei[t100] * sfPhi[t200]
655 + fmpolei[t010] * sfPhi[t110] + fmpolei[t001] * sfPhi[t101]
656 + fmpolei[t200] * sfPhi[t300] + fmpolei[t020] * sfPhi[t120]
657 + fmpolei[t002] * sfPhi[t102] + fmpolei[t110] * sfPhi[t210]
658 + fmpolei[t101] * sfPhi[t201] + fmpolei[t011] * sfPhi[t111];
659 gy += fmpolei[t000] * sfPhi[t010] + fmpolei[t100] * sfPhi[t110]
660 + fmpolei[t010] * sfPhi[t020] + fmpolei[t001] * sfPhi[t011]
661 + fmpolei[t200] * sfPhi[t210] + fmpolei[t020] * sfPhi[t030]
662 + fmpolei[t002] * sfPhi[t012] + fmpolei[t110] * sfPhi[t120]
663 + fmpolei[t101] * sfPhi[t111] + fmpolei[t011] * sfPhi[t021];
664 gz += fmpolei[t000] * sfPhi[t001] + fmpolei[t100] * sfPhi[t101]
665 + fmpolei[t010] * sfPhi[t011] + fmpolei[t001] * sfPhi[t002]
666 + fmpolei[t200] * sfPhi[t201] + fmpolei[t020] * sfPhi[t021]
667 + fmpolei[t002] * sfPhi[t003] + fmpolei[t110] * sfPhi[t111]
668 + fmpolei[t101] * sfPhi[t102] + fmpolei[t011] * sfPhi[t012];
669 gx *= nfftX;
670 gy *= nfftY;
671 gz *= nfftZ;
672 double[][] recip = crystal.getUnitCell().A;
673 double dfx = recip[0][0] * gx + recip[0][1] * gy + recip[0][2] * gz;
674 double dfy = recip[1][0] * gx + recip[1][1] * gy + recip[1][2] * gz;
675 double dfz = recip[2][0] * gx + recip[2][1] * gy + recip[2][2] * gz;
676 dfx *= 0.5 * electric;
677 dfy *= 0.5 * electric;
678 dfz *= 0.5 * electric;
679
680 double tqx = -mpolei[t010] * sPhi[t001] + mpolei[t001] * sPhi[t010];
681 double tqy = -mpolei[t001] * sPhi[t100] + mpolei[t100] * sPhi[t001];
682 double tqz = -mpolei[t100] * sPhi[t010] + mpolei[t010] * sPhi[t100];
683
684 tqx -= twoThirds * (
685 mpolei[t110] * sPhi[t101] + mpolei[t020] * sPhi[t011] + mpolei[t011] * sPhi[t002]
686 - mpolei[t101] * sPhi[t110] - mpolei[t011] * sPhi[t020]
687 - mpolei[t002] * sPhi[t011]);
688 tqy -= twoThirds * (
689 mpolei[t101] * sPhi[t200] + mpolei[t011] * sPhi[t110] + mpolei[t002] * sPhi[t101]
690 - mpolei[t200] * sPhi[t101] - mpolei[t110] * sPhi[t011]
691 - mpolei[t101] * sPhi[t002]);
692 tqz -= twoThirds * (
693 mpolei[t200] * sPhi[t110] + mpolei[t110] * sPhi[t020] + mpolei[t101] * sPhi[t011]
694 - mpolei[t110] * sPhi[t200] - mpolei[t020] * sPhi[t110]
695 - mpolei[t011] * sPhi[t101]);
696 tqx *= electric;
697 tqy *= electric;
698 tqz *= electric;
699 grad.add(threadID, i, polarizationScale * dfx, polarizationScale * dfy,
700 polarizationScale * dfz);
701 torque.add(threadID, i, polarizationScale * tqx, polarizationScale * tqy,
702 polarizationScale * tqz);
703 if (lambdaTerm) {
704 double factor = dEdLSign * dlPowPol;
705 lambdaGrad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
706 lambdaTorque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
707 }
708 if (esvTerm && extendedSystem.isTitrating(i)) {
709 final double[] mTitrDot = titrationMultipole[0][i];
710 double eTitrDot = mTitrDot[t000] * sPhi[t000]
711 + mTitrDot[t100] * sPhi[t100] + mTitrDot[t010] * sPhi[t010] + mTitrDot[t001] * sPhi[t001]
712 + oneThird * (mTitrDot[t200] * sPhi[t200] + mTitrDot[t020] * sPhi[t020] + mTitrDot[t002] * sPhi[t002]
713 + 2.0 * (mTitrDot[t110] * sPhi[t110] + mTitrDot[t101] * sPhi[t101] + mTitrDot[t011] * sPhi[t011]));
714 double eTautDot = 0.0;
715 if (extendedSystem.isTautomerizing(i)) {
716 final double[] mTautDot = tautomerMultipole[0][i];
717 eTautDot = mTautDot[t000] * sPhi[t000]
718 + mTautDot[t100] * sPhi[t100] + mTautDot[t010] * sPhi[t010] + mTautDot[t001] * sPhi[t001]
719 + oneThird * (mTautDot[t200] * sPhi[t200] + mTautDot[t020] * sPhi[t020] + mTautDot[t002] * sPhi[t002]
720 + 2.0 * (mTautDot[t110] * sPhi[t110] + mTautDot[t101] * sPhi[t101] + mTautDot[t011] * sPhi[t011]));
721 }
722 double factor = polarizationScale * electric;
723 extendedSystem.addIndElecDeriv(i, eTitrDot * factor, eTautDot * factor);
724 }
725 }
726 }
727 }
728 eRecip *= 0.5 * electric;
729 if (lambdaTerm) {
730 shareddEdLambda.addAndGet(dEdLSign * dlPowPol * eRecip);
731 sharedd2EdLambda2.addAndGet(dEdLSign * d2lPowPol * eRecip);
732 }
733 }
734
735 @Override
736 public IntegerSchedule schedule() {
737 return IntegerSchedule.fixed();
738 }
739
740 @Override
741 public void start() {
742 eSelf = 0.0;
743 eRecip = 0.0;
744 threadID = getThreadIndex();
745 }
746 }
747 }