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 edu.rit.pj.IntegerForLoop;
41 import edu.rit.pj.IntegerSchedule;
42 import edu.rit.pj.ParallelRegion;
43 import edu.rit.pj.ParallelTeam;
44 import edu.rit.pj.reduction.SharedDouble;
45 import edu.rit.pj.reduction.SharedInteger;
46 import ffx.crystal.Crystal;
47 import ffx.crystal.SymOp;
48 import ffx.numerics.atomic.AtomicDoubleArray3D;
49 import ffx.potential.bonded.Atom;
50 import ffx.potential.extended.ExtendedSystem;
51 import ffx.potential.nonbonded.MaskingInterface;
52 import ffx.potential.parameters.ForceField;
53 import ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition;
54 import ffx.potential.utils.EnergyException;
55
56 import java.util.List;
57 import java.util.logging.Level;
58 import java.util.logging.Logger;
59
60 import static ffx.numerics.special.Erf.erfc;
61 import static ffx.potential.parameters.MultipoleType.t000;
62 import static ffx.potential.parameters.MultipoleType.t001;
63 import static ffx.potential.parameters.MultipoleType.t002;
64 import static ffx.potential.parameters.MultipoleType.t010;
65 import static ffx.potential.parameters.MultipoleType.t011;
66 import static ffx.potential.parameters.MultipoleType.t020;
67 import static ffx.potential.parameters.MultipoleType.t100;
68 import static ffx.potential.parameters.MultipoleType.t101;
69 import static ffx.potential.parameters.MultipoleType.t110;
70 import static ffx.potential.parameters.MultipoleType.t200;
71 import static java.lang.Double.isInfinite;
72 import static java.lang.Double.isNaN;
73 import static java.lang.String.format;
74 import static java.util.Arrays.fill;
75 import static org.apache.commons.math3.util.FastMath.exp;
76 import static org.apache.commons.math3.util.FastMath.min;
77 import static org.apache.commons.math3.util.FastMath.sqrt;
78
79
80
81
82
83
84
85 public class RealSpaceEnergyRegion extends ParallelRegion implements MaskingInterface {
86 private static final Logger logger = Logger.getLogger(RealSpaceEnergyRegion.class.getName());
87
88
89
90 private static final double oneThird = 1.0 / 3.0;
91 private final int maxThreads;
92 private final double electric;
93 private final SharedInteger sharedInteractions;
94 private final ForceField.ELEC_FORM elecForm;
95 private final RealSpaceEnergyLoop[] realSpaceEnergyLoop;
96 private final FixedChargeEnergyLoop[] fixedChargeEnergyLoop;
97
98
99
100 private final boolean intermolecularSoftcore;
101
102
103
104 private final boolean intramolecularSoftcore;
105
106
107
108 public double[][][] inducedDipole;
109 public double[][][] inducedDipoleCR;
110
111
112
113 protected int[][] ip11;
114
115
116
117 private Atom[] atoms;
118
119
120
121 private Crystal crystal;
122
123
124
125 private double[][][] coordinates;
126
127
128
129 private MultipoleFrameDefinition[] frame;
130
131
132
133 private int[][] axisAtom;
134
135
136
137 private double[][][] globalMultipole;
138 private double[][][] titrationMultipole;
139 private double[][][] tautomerMultipole;
140
141 private ExtendedSystem extendedSystem = null;
142 private boolean esvTerm = false;
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161 private boolean[] use;
162
163
164
165 private int[] molecule;
166
167
168
169 private boolean[] isSoft;
170 private double[] ipdamp;
171 private double[] thole;
172
173
174
175 private int[][] mask12;
176 private int[][] mask13;
177 private int[][] mask14;
178 private int[][] mask15;
179
180
181
182 private int[][][] realSpaceLists;
183
184
185
186 private int[][] realSpaceCounts;
187
188
189
190 private IntegerSchedule realSpaceSchedule;
191 private long[] realSpaceEnergyTime;
192
193
194
195 private AtomicDoubleArray3D grad;
196
197
198
199 private AtomicDoubleArray3D torque;
200
201
202
203 private AtomicDoubleArray3D lambdaGrad;
204
205
206
207 private AtomicDoubleArray3D lambdaTorque;
208
209
210
211 private SharedDouble shareddEdLambda;
212
213
214
215 private SharedDouble sharedd2EdLambda2;
216
217
218
219 private boolean gradient;
220
221
222
223
224 private boolean lambdaTerm;
225
226
227
228 private boolean nnTerm;
229
230
231
232
233 private LambdaMode lambdaMode = LambdaMode.OFF;
234 private Polarization polarization;
235
236
237
238 private double lAlpha = 0.0;
239 private double dlAlpha = 0.0;
240 private double d2lAlpha = 0.0;
241 private double dEdLSign = 1.0;
242
243
244
245 private double dlPowPerm = 0.0;
246 private double d2lPowPerm = 0.0;
247 private boolean doPermanentRealSpace;
248 private double permanentScale = 1.0;
249
250
251
252 private double lPowPol = 1.0;
253 private double dlPowPol = 0.0;
254 private double d2lPowPol = 0.0;
255 private boolean doPolarization;
256
257
258
259
260
261
262
263
264
265
266
267
268
269 private double polarizationScale = 1.0;
270
271
272 private double aewald;
273 private double an0;
274 private double an1;
275 private double an2;
276 private double an3;
277 private double an4;
278 private double an5;
279 private double permanentEnergy;
280 private double polarizationEnergy;
281 private ScaleParameters scaleParameters;
282
283 public RealSpaceEnergyRegion(ForceField.ELEC_FORM elecForm, int nt, ForceField forceField, boolean lambdaTerm, double electric) {
284 maxThreads = nt;
285 this.electric = electric;
286 this.elecForm = elecForm;
287 sharedInteractions = new SharedInteger();
288 if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
289 fixedChargeEnergyLoop = new FixedChargeEnergyLoop[nt];
290 realSpaceEnergyLoop = null;
291 } else {
292 fixedChargeEnergyLoop = null;
293 realSpaceEnergyLoop = new RealSpaceEnergyLoop[nt];
294 }
295
296
297 if (lambdaTerm) {
298 intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
299 intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
300 } else {
301 intermolecularSoftcore = false;
302 intramolecularSoftcore = false;
303 }
304 }
305
306 @Override
307 public void applyMask(int i, boolean[] is14, double[]... masks) {
308 if (ip11[i] != null) {
309
310 double[] permanentEnergyMask = masks[0];
311 double[] permanentFieldMask = masks[1];
312 for (int value : mask12[i]) {
313 permanentFieldMask[value] = scaleParameters.p12scale;
314 permanentEnergyMask[value] = scaleParameters.m12scale;
315 }
316 for (int value : mask13[i]) {
317 permanentFieldMask[value] = scaleParameters.p13scale;
318 permanentEnergyMask[value] = scaleParameters.m13scale;
319 }
320 for (int value : mask14[i]) {
321 permanentFieldMask[value] = scaleParameters.p14scale;
322 permanentEnergyMask[value] = scaleParameters.m14scale;
323 for (int k : ip11[i]) {
324 if (k == value) {
325 permanentFieldMask[value] = scaleParameters.intra14Scale * scaleParameters.p14scale;
326 break;
327 }
328 }
329 }
330 for (int value : mask15[i]) {
331 permanentEnergyMask[value] = scaleParameters.m15scale;
332 }
333
334 double[] polarizationGroupMask = masks[2];
335 for (int j : ip11[i]) {
336 polarizationGroupMask[j] = scaleParameters.d11scale;
337 }
338 } else {
339
340 double[] permanentEnergyMask = masks[0];
341 for (int value : mask12[i]) {
342 permanentEnergyMask[value] = scaleParameters.m12scale;
343 }
344 for (int value : mask13[i]) {
345 permanentEnergyMask[value] = scaleParameters.m13scale;
346 }
347 for (int value : mask14[i]) {
348 permanentEnergyMask[value] = scaleParameters.m14scale;
349 }
350 }
351 }
352
353
354
355
356
357
358 public void executeWith(ParallelTeam parallelTeam) {
359 try {
360 parallelTeam.execute(this);
361 } catch (Exception e) {
362 String message = " Exception computing the electrostatic energy.\n";
363 logger.log(Level.SEVERE, message, e);
364 }
365 }
366
367 @Override
368 public void finish() {
369 permanentEnergy = 0.0;
370 polarizationEnergy = 0.0;
371 if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
372 for (int i = 0; i < maxThreads; i++) {
373 permanentEnergy += fixedChargeEnergyLoop[i].permanentEnergy;
374 }
375 } else {
376 for (int i = 0; i < maxThreads; i++) {
377 permanentEnergy += realSpaceEnergyLoop[i].permanentEnergy;
378 polarizationEnergy += realSpaceEnergyLoop[i].inducedEnergy;
379 }
380 }
381 permanentEnergy *= electric;
382 polarizationEnergy *= electric;
383 if (isNaN(permanentEnergy)) {
384 throw new EnergyException(
385 format(" The permanent multipole energy is %16.8f", permanentEnergy));
386 }
387 if (isNaN(polarizationEnergy)) {
388 throw new EnergyException(
389 format(" The polarization energy is %16.8f", polarizationEnergy));
390 }
391 }
392
393 public int getCount(int i) {
394 if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
395 return fixedChargeEnergyLoop[i].getCount();
396 } else {
397 return realSpaceEnergyLoop[i].getCount();
398 }
399 }
400
401 public int getInteractions() {
402 return sharedInteractions.get();
403 }
404
405 public double getPermanentEnergy() {
406 return permanentEnergy;
407 }
408
409 public double getPolarizationEnergy() {
410 return polarizationEnergy;
411 }
412
413 public void init(
414 Atom[] atoms,
415 Crystal crystal,
416 ExtendedSystem extendedSystem,
417 boolean esvTerm,
418 double[][][] coordinates,
419 MultipoleFrameDefinition[] frame,
420 int[][] axisAtom,
421 double[][][] globalMultipole,
422 double[][][] titrationMultipole,
423 double[][][] tautomerMultipole,
424 double[][][] inducedDipole,
425 double[][][] inducedDipoleCR,
426 boolean[] use,
427 int[] molecule,
428 int[][] ip11,
429 int[][] mask12,
430 int[][] mask13,
431 int[][] mask14,
432 int[][] mask15,
433 boolean[] isSoft,
434 double[] ipdamp,
435 double[] thole,
436 RealSpaceNeighborParameters realSpaceNeighborParameters,
437 boolean gradient,
438 boolean lambdaTerm,
439 boolean nnTerm,
440 LambdaMode lambdaMode,
441 Polarization polarization,
442 EwaldParameters ewaldParameters,
443 ScaleParameters scaleParameters,
444 AlchemicalParameters alchemicalParameters,
445 long[] realSpaceEnergyTime,
446 AtomicDoubleArray3D grad,
447 AtomicDoubleArray3D torque,
448 AtomicDoubleArray3D lambdaGrad,
449 AtomicDoubleArray3D lambdaTorque,
450 SharedDouble shareddEdLambda,
451 SharedDouble sharedd2EdLambda2) {
452
453 this.atoms = atoms;
454 this.crystal = crystal;
455 this.extendedSystem = extendedSystem;
456 this.esvTerm = esvTerm;
457 this.coordinates = coordinates;
458 this.frame = frame;
459 this.axisAtom = axisAtom;
460 this.globalMultipole = globalMultipole;
461 this.titrationMultipole = titrationMultipole;
462 this.tautomerMultipole = tautomerMultipole;
463 this.inducedDipole = inducedDipole;
464 this.inducedDipoleCR = inducedDipoleCR;
465 this.use = use;
466 this.molecule = molecule;
467 this.ip11 = ip11;
468 this.mask12 = mask12;
469 this.mask13 = mask13;
470 this.mask14 = mask14;
471 this.mask15 = mask15;
472 this.isSoft = isSoft;
473 this.ipdamp = ipdamp;
474 this.thole = thole;
475 this.realSpaceLists = realSpaceNeighborParameters.realSpaceLists;
476 this.realSpaceCounts = realSpaceNeighborParameters.realSpaceCounts;
477 this.realSpaceSchedule = realSpaceNeighborParameters.realSpaceSchedule;
478 this.gradient = gradient;
479 this.lambdaTerm = lambdaTerm;
480 this.nnTerm = nnTerm;
481 this.lambdaMode = lambdaMode;
482 this.polarization = polarization;
483 this.lAlpha = alchemicalParameters.lAlpha;
484 this.dlAlpha = alchemicalParameters.dlAlpha;
485 this.d2lAlpha = alchemicalParameters.d2lAlpha;
486 this.dEdLSign = alchemicalParameters.dEdLSign;
487 this.dlPowPerm = alchemicalParameters.dlPowPerm;
488 this.d2lPowPerm = alchemicalParameters.d2lPowPerm;
489 this.doPermanentRealSpace = alchemicalParameters.doPermanentRealSpace;
490 this.permanentScale = alchemicalParameters.permanentScale;
491 this.lPowPol = alchemicalParameters.lPowPol;
492 this.dlPowPol = alchemicalParameters.dlPowPol;
493 this.d2lPowPol = alchemicalParameters.d2lPowPol;
494 this.doPolarization = alchemicalParameters.doPolarization;
495 this.polarizationScale = alchemicalParameters.polarizationScale;
496 this.aewald = ewaldParameters.aewald;
497 this.an0 = ewaldParameters.an0;
498 this.an1 = ewaldParameters.an1;
499 this.an2 = ewaldParameters.an2;
500 this.an3 = ewaldParameters.an3;
501 this.an4 = ewaldParameters.an4;
502 this.an5 = ewaldParameters.an5;
503 this.scaleParameters = scaleParameters;
504
505 this.realSpaceEnergyTime = realSpaceEnergyTime;
506 this.grad = grad;
507 this.torque = torque;
508 this.lambdaGrad = lambdaGrad;
509 this.lambdaTorque = lambdaTorque;
510 this.shareddEdLambda = shareddEdLambda;
511 this.sharedd2EdLambda2 = sharedd2EdLambda2;
512 }
513
514 @Override
515 public void removeMask(int i, boolean[] is14, double[]... masks) {
516 if (ip11[i] != null) {
517
518 double[] permanentEnergyMask = masks[0];
519 double[] permanentFieldMask = masks[1];
520 for (int value : mask12[i]) {
521 permanentFieldMask[value] = 1.0;
522 permanentEnergyMask[value] = 1.0;
523 }
524 for (int value : mask13[i]) {
525 permanentFieldMask[value] = 1.0;
526 permanentEnergyMask[value] = 1.0;
527 }
528 for (int value : mask14[i]) {
529 permanentFieldMask[value] = 1.0;
530 permanentEnergyMask[value] = 1.0;
531 for (int k : ip11[i]) {
532 if (k == value) {
533 permanentFieldMask[value] = 1.0;
534 break;
535 }
536 }
537 }
538 for (int value : mask15[i]) {
539 permanentEnergyMask[value] = 1.0;
540 }
541
542 double[] polarizationGroupMask = masks[2];
543 for (int j : ip11[i]) {
544 polarizationGroupMask[j] = 1.0;
545 }
546 } else {
547
548 double[] permanentEnergyMask = masks[0];
549 for (int value : mask12[i]) {
550 permanentEnergyMask[value] = 1.0;
551 }
552 for (int value : mask13[i]) {
553 permanentEnergyMask[value] = 1.0;
554 }
555 for (int value : mask14[i]) {
556 permanentEnergyMask[value] = 1.0;
557 }
558 }
559 }
560
561 @Override
562 public void run() {
563 int threadIndex = getThreadIndex();
564 if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
565 if (fixedChargeEnergyLoop[threadIndex] == null) {
566 fixedChargeEnergyLoop[threadIndex] = new FixedChargeEnergyLoop();
567 }
568 try {
569 int nAtoms = atoms.length;
570 execute(0, nAtoms - 1, fixedChargeEnergyLoop[threadIndex]);
571 } catch (Exception e) {
572 String message =
573 "Fatal exception computing the real space energy in thread " + getThreadIndex() + "\n";
574 logger.log(Level.SEVERE, message, e);
575 }
576 } else {
577 if (realSpaceEnergyLoop[threadIndex] == null) {
578 realSpaceEnergyLoop[threadIndex] = new RealSpaceEnergyLoop();
579 }
580 try {
581 int nAtoms = atoms.length;
582 execute(0, nAtoms - 1, realSpaceEnergyLoop[threadIndex]);
583 } catch (Exception e) {
584 String message =
585 "Fatal exception computing the real space energy in thread " + getThreadIndex() + "\n";
586 logger.log(Level.SEVERE, message, e);
587 }
588 }
589 }
590
591 @Override
592 public void start() {
593 sharedInteractions.set(0);
594 }
595
596
597
598
599
600
601
602
603
604
605 private void log(int i, int k, double r, double eij, int count, double total) {
606 logger.info(format("%s %6d %6d %12.10f %12.10f %8d %16.10f",
607 "ELEC", atoms[i].getIndex(), atoms[k].getIndex(), r, eij, count, total));
608 }
609
610
611
612
613
614 private class RealSpaceEnergyLoop extends IntegerForLoop {
615
616 private final double[] dx_local;
617 private final double[][] rot_local;
618 private final Torque torques;
619 private double ci;
620 private double dix, diy, diz;
621 private double qixx, qiyy, qizz, qixy, qixz, qiyz;
622 private double ck;
623 private double dkx, dky, dkz;
624 private double qkxx, qkyy, qkzz, qkxy, qkxz, qkyz;
625 private double uix, uiy, uiz;
626 private double pix, piy, piz;
627 private double xr, yr, zr;
628 private double ukx, uky, ukz;
629 private double pkx, pky, pkz;
630 private double bn0, bn1, bn2, bn3, bn4, bn5, bn6;
631 private double rr1, rr3, rr5, rr7, rr9, rr11, rr13;
632 private double scale, scale3, scale5, scale7;
633 private double scalep, scaled;
634 private double ddsc3x, ddsc3y, ddsc3z;
635 private double ddsc5x, ddsc5y, ddsc5z;
636 private double ddsc7x, ddsc7y, ddsc7z;
637 private double l2;
638 private boolean soft;
639 private double selfScale;
640 private double permanentEnergy;
641 private double inducedEnergy;
642 private double dUdL, d2UdL2;
643 private int i, k, iSymm, count;
644 private double[] gxk_local, gyk_local, gzk_local;
645 private double[] txk_local, tyk_local, tzk_local;
646 private double[] lxk_local, lyk_local, lzk_local;
647 private double[] ltxk_local, ltyk_local, ltzk_local;
648 private double[] masking_local;
649 private double[] maskingp_local;
650 private double[] maskingd_local;
651 private int threadID;
652
653 RealSpaceEnergyLoop() {
654 super();
655 dx_local = new double[3];
656 rot_local = new double[3][3];
657 torques = new Torque();
658 }
659
660 @Override
661 public void finish() {
662 sharedInteractions.addAndGet(count);
663 if (lambdaTerm) {
664 shareddEdLambda.addAndGet(dUdL * electric);
665 sharedd2EdLambda2.addAndGet(d2UdL2 * electric);
666 }
667 realSpaceEnergyTime[getThreadIndex()] += System.nanoTime();
668 }
669
670 public int getCount() {
671 return count;
672 }
673
674 @Override
675 public void run(int lb, int ub) {
676 List<SymOp> symOps = crystal.spaceGroup.symOps;
677 int nSymm = symOps.size();
678 int nAtoms = atoms.length;
679 for (iSymm = 0; iSymm < nSymm; iSymm++) {
680 SymOp symOp = symOps.get(iSymm);
681 if (gradient) {
682 fill(gxk_local, 0.0);
683 fill(gyk_local, 0.0);
684 fill(gzk_local, 0.0);
685 fill(txk_local, 0.0);
686 fill(tyk_local, 0.0);
687 fill(tzk_local, 0.0);
688 }
689 if (lambdaTerm) {
690 fill(lxk_local, 0.0);
691 fill(lyk_local, 0.0);
692 fill(lzk_local, 0.0);
693 fill(ltxk_local, 0.0);
694 fill(ltyk_local, 0.0);
695 fill(ltzk_local, 0.0);
696 }
697 realSpaceChunk(lb, ub);
698 if (gradient) {
699
700 int[] frameIndex = {-1, -1, -1, -1};
701 double[][] g = new double[4][3];
702 double[] trq = new double[3];
703 for (int i = 0; i < nAtoms; i++) {
704 fill(frameIndex, -1);
705 for (int j = 0; j < 4; j++) {
706 fill(g[j], 0.0);
707 }
708 trq[0] = txk_local[i];
709 trq[1] = tyk_local[i];
710 trq[2] = tzk_local[i];
711 torques.torque(i, iSymm, trq, frameIndex, g);
712 for (int j = 0; j < 4; j++) {
713 int index = frameIndex[j];
714 if (index >= 0) {
715 double[] gj = g[j];
716 gxk_local[index] += gj[0];
717 gyk_local[index] += gj[1];
718 gzk_local[index] += gj[2];
719 }
720 }
721 }
722
723 if (iSymm != 0) {
724 crystal.applyTransSymRot(
725 nAtoms, gxk_local, gyk_local, gzk_local, gxk_local, gyk_local, gzk_local, symOp,
726 rot_local);
727 }
728
729 for (int j = 0; j < nAtoms; j++) {
730 grad.add(threadID, j, gxk_local[j], gyk_local[j], gzk_local[j]);
731 }
732 }
733 if (lambdaTerm) {
734
735 int[] frameIndex = {-1, -1, -1, -1};
736 double[][] g = new double[4][3];
737 double[] trq = new double[3];
738 for (int i = 0; i < nAtoms; i++) {
739 fill(frameIndex, -1);
740 for (int j = 0; j < 4; j++) {
741 fill(g[j], 0.0);
742 }
743 trq[0] = ltxk_local[i];
744 trq[1] = ltyk_local[i];
745 trq[2] = ltzk_local[i];
746 torques.torque(i, iSymm, trq, frameIndex, g);
747 for (int j = 0; j < 4; j++) {
748 int index = frameIndex[j];
749 if (index >= 0) {
750 double[] gj = g[j];
751 lxk_local[index] += gj[0];
752 lyk_local[index] += gj[1];
753 lzk_local[index] += gj[2];
754 }
755 }
756 }
757
758 if (iSymm != 0) {
759 crystal.applyTransSymRot(
760 nAtoms, lxk_local, lyk_local, lzk_local, lxk_local, lyk_local, lzk_local, symOp,
761 rot_local);
762 }
763
764 for (int j = 0; j < nAtoms; j++) {
765 lambdaGrad.add(threadID, j, lxk_local[j], lyk_local[j], lzk_local[j]);
766 }
767 }
768 }
769 }
770
771 @Override
772 public IntegerSchedule schedule() {
773 return realSpaceSchedule;
774 }
775
776 @Override
777 public void start() {
778 init();
779 threadID = getThreadIndex();
780 realSpaceEnergyTime[threadID] -= System.nanoTime();
781 permanentEnergy = 0.0;
782 inducedEnergy = 0.0;
783 count = 0;
784 if (lambdaTerm) {
785 dUdL = 0.0;
786 d2UdL2 = 0.0;
787 }
788 torques.init(axisAtom, frame, coordinates);
789 }
790
791 private void init() {
792 int nAtoms = atoms.length;
793 if (masking_local == null || masking_local.length < nAtoms) {
794 txk_local = new double[nAtoms];
795 tyk_local = new double[nAtoms];
796 tzk_local = new double[nAtoms];
797 gxk_local = new double[nAtoms];
798 gyk_local = new double[nAtoms];
799 gzk_local = new double[nAtoms];
800 lxk_local = new double[nAtoms];
801 lyk_local = new double[nAtoms];
802 lzk_local = new double[nAtoms];
803 ltxk_local = new double[nAtoms];
804 ltyk_local = new double[nAtoms];
805 ltzk_local = new double[nAtoms];
806 masking_local = new double[nAtoms];
807 maskingp_local = new double[nAtoms];
808 maskingd_local = new double[nAtoms];
809 fill(masking_local, 1.0);
810 fill(maskingp_local, 1.0);
811 fill(maskingd_local, 1.0);
812 }
813 }
814
815
816
817
818
819
820
821 private void realSpaceChunk(final int lb, final int ub) {
822 final double[] x = coordinates[0][0];
823 final double[] y = coordinates[0][1];
824 final double[] z = coordinates[0][2];
825 final double[][] mpole = globalMultipole[0];
826 final double[] zeropole = new double[10];
827 final double[][] ind = inducedDipole[0];
828 final double[][] indp = inducedDipoleCR[0];
829 final int[][] lists = realSpaceLists[iSymm];
830 final double[] neighborX = coordinates[iSymm][0];
831 final double[] neighborY = coordinates[iSymm][1];
832 final double[] neighborZ = coordinates[iSymm][2];
833 final double[][] neighborMultipole = globalMultipole[iSymm];
834 final double[][] neighborInducedDipole = inducedDipole[iSymm];
835 final double[][] neighborInducedDipolep = inducedDipoleCR[iSymm];
836 for (i = lb; i <= ub; i++) {
837 if (!use[i]) {
838 continue;
839 }
840 final int moleculei = molecule[i];
841 if (iSymm == 0) {
842 applyMask(i, null, masking_local, maskingp_local, maskingd_local);
843 }
844 final double xi = x[i];
845 final double yi = y[i];
846 final double zi = z[i];
847 final double[] globalMultipolei = mpole[i];
848 final double[] inducedDipolei = ind[i];
849 final double[] inducedDipolepi = indp[i];
850 setMultipoleI(globalMultipolei);
851 setInducedI(inducedDipolei);
852 setInducedpI(inducedDipolepi);
853 final boolean softi = isSoft[i];
854 final double pdi = ipdamp[i];
855 final double pti = thole[i];
856 final int[] list = lists[i];
857 final int npair = realSpaceCounts[iSymm][i];
858 for (int j = 0; j < npair; j++) {
859 k = list[j];
860 if (!use[k]) {
861 continue;
862 }
863 boolean sameMolecule = (moleculei == molecule[k]);
864 if (lambdaMode == LambdaMode.VAPOR) {
865 if ((intermolecularSoftcore && !sameMolecule)
866 || (intramolecularSoftcore && sameMolecule)) {
867 continue;
868 }
869 }
870 selfScale = 1.0;
871 if (i == k) {
872 selfScale = 0.5;
873 }
874 double beta = 0.0;
875 l2 = 1.0;
876 soft = (softi || isSoft[k]);
877 if (soft && doPermanentRealSpace) {
878 beta = lAlpha;
879 l2 = permanentScale;
880 } else if (nnTerm) {
881 l2 = permanentScale;
882 }
883 final double xk = neighborX[k];
884 final double yk = neighborY[k];
885 final double zk = neighborZ[k];
886 dx_local[0] = xk - xi;
887 dx_local[1] = yk - yi;
888 dx_local[2] = zk - zi;
889 double r2 = crystal.image(dx_local);
890 xr = dx_local[0];
891 yr = dx_local[1];
892 zr = dx_local[2];
893 final double[] globalMultipolek = neighborMultipole[k];
894 final double[] inducedDipolek = neighborInducedDipole[k];
895 final double[] inducedDipolepk = neighborInducedDipolep[k];
896 setMultipoleK(globalMultipolek);
897 setInducedK(inducedDipolek);
898 setInducedpK(inducedDipolepk);
899 final double pdk = ipdamp[k];
900 final double ptk = thole[k];
901 scale = masking_local[k];
902 scalep = maskingp_local[k];
903 scaled = maskingd_local[k];
904 scale3 = 1.0;
905 scale5 = 1.0;
906 scale7 = 1.0;
907 double r = sqrt(r2 + beta);
908 double ralpha = aewald * r;
909 double exp2a = exp(-ralpha * ralpha);
910 rr1 = 1.0 / r;
911 double rr2 = rr1 * rr1;
912 bn0 = erfc(ralpha) * rr1;
913 bn1 = (bn0 + an0 * exp2a) * rr2;
914 bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
915 bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
916 bn4 = (7.0 * bn3 + an3 * exp2a) * rr2;
917 bn5 = (9.0 * bn4 + an4 * exp2a) * rr2;
918 bn6 = (11.0 * bn5 + an5 * exp2a) * rr2;
919 rr3 = rr1 * rr2;
920 rr5 = 3.0 * rr3 * rr2;
921 rr7 = 5.0 * rr5 * rr2;
922 rr9 = 7.0 * rr7 * rr2;
923 rr11 = 9.0 * rr9 * rr2;
924 rr13 = 11.0 * rr11 * rr2;
925 ddsc3x = 0.0;
926 ddsc3y = 0.0;
927 ddsc3z = 0.0;
928 ddsc5x = 0.0;
929 ddsc5y = 0.0;
930 ddsc5z = 0.0;
931 ddsc7x = 0.0;
932 ddsc7y = 0.0;
933 ddsc7z = 0.0;
934 double damp = pdi * pdk;
935 double thole = min(pti, ptk);
936 double rdamp = r * damp;
937 damp = -thole * rdamp * rdamp * rdamp;
938 if (damp > -50.0) {
939 final double expdamp = exp(damp);
940 scale3 = 1.0 - expdamp;
941 scale5 = 1.0 - expdamp * (1.0 - damp);
942 scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
943 final double temp3 = -3.0 * damp * expdamp * rr2;
944 final double temp5 = -damp;
945 final double temp7 = -0.2 - 0.6 * damp;
946 ddsc3x = temp3 * xr;
947 ddsc3y = temp3 * yr;
948 ddsc3z = temp3 * zr;
949 ddsc5x = temp5 * ddsc3x;
950 ddsc5y = temp5 * ddsc3y;
951 ddsc5z = temp5 * ddsc3z;
952 ddsc7x = temp7 * ddsc5x;
953 ddsc7y = temp7 * ddsc5y;
954 ddsc7z = temp7 * ddsc5z;
955 }
956 if (doPermanentRealSpace) {
957 double ei = permanentPair(gradient, lambdaTerm);
958 if (isNaN(ei) || isInfinite(ei)) {
959 String message = format(" %s\n %s\n %s\n The permanent multipole energy between "
960 + "atoms %d and %d (%d) is %16.8f at %16.8f A.",
961 crystal.getUnitCell().toString(), atoms[i].toString(), atoms[k].toString(),
962 i, k, iSymm, ei, r);
963 throw new EnergyException(message);
964 }
965 if (ei != 0.0) {
966 permanentEnergy += ei;
967 count++;
968
969 }
970 if (esvTerm) {
971 if (extendedSystem.isTitrating(i)) {
972 double titrdUdL = 0.0;
973 double tautdUdL = 0.0;
974 final double[][] titrMpole = titrationMultipole[0];
975 final double[] titrMultipolei = titrMpole[i];
976 setMultipoleI(titrMultipolei);
977 titrdUdL = electric * permanentPair(false, false);
978 if (extendedSystem.isTautomerizing(i)) {
979 final double[][] tautMpole = tautomerMultipole[0];
980 final double[] tautMultipolei = tautMpole[i];
981 setMultipoleI(tautMultipolei);
982 tautdUdL = electric * permanentPair(false, false);
983 }
984 extendedSystem.addPermElecDeriv(i, titrdUdL, tautdUdL);
985
986 setMultipoleI(globalMultipolei);
987 }
988 if (extendedSystem.isTitrating(k)) {
989 double titrdUdL = 0.0;
990 double tautdUdL = 0.0;
991 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
992 final double[] titrMultipolek = titrNeighborMpole[k];
993 setMultipoleK(titrMultipolek);
994 titrdUdL = electric * permanentPair(false, false);
995 if (extendedSystem.isTautomerizing(k)) {
996 final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
997 final double[] tautMultipolek = tautNeighborMpole[k];
998 setMultipoleK(tautMultipolek);
999 tautdUdL = electric * permanentPair(false, false);
1000 }
1001 extendedSystem.addPermElecDeriv(k, titrdUdL, tautdUdL);
1002
1003 setMultipoleK(globalMultipolek);
1004 }
1005 }
1006 }
1007 if (polarization != Polarization.NONE && doPolarization) {
1008
1009 if (soft && doPermanentRealSpace) {
1010 scale3 = 1.0;
1011 scale5 = 1.0;
1012 scale7 = 1.0;
1013 r = sqrt(r2);
1014 ralpha = aewald * r;
1015 exp2a = exp(-ralpha * ralpha);
1016 rr1 = 1.0 / r;
1017 rr2 = rr1 * rr1;
1018 bn0 = erfc(ralpha) * rr1;
1019 bn1 = (bn0 + an0 * exp2a) * rr2;
1020 bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
1021 bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
1022 bn4 = (7.0 * bn3 + an3 * exp2a) * rr2;
1023 bn5 = (9.0 * bn4 + an4 * exp2a) * rr2;
1024 bn6 = (11.0 * bn5 + an5 * exp2a) * rr2;
1025 rr3 = rr1 * rr2;
1026 rr5 = 3.0 * rr3 * rr2;
1027 rr7 = 5.0 * rr5 * rr2;
1028 rr9 = 7.0 * rr7 * rr2;
1029 rr11 = 9.0 * rr9 * rr2;
1030 ddsc3x = 0.0;
1031 ddsc3y = 0.0;
1032 ddsc3z = 0.0;
1033 ddsc5x = 0.0;
1034 ddsc5y = 0.0;
1035 ddsc5z = 0.0;
1036 ddsc7x = 0.0;
1037 ddsc7y = 0.0;
1038 ddsc7z = 0.0;
1039 damp = pdi * pdk;
1040 thole = min(pti, ptk);
1041 rdamp = r * damp;
1042 damp = -thole * rdamp * rdamp * rdamp;
1043 if (damp > -50.0) {
1044 final double expdamp = exp(damp);
1045 scale3 = 1.0 - expdamp;
1046 scale5 = 1.0 - expdamp * (1.0 - damp);
1047 scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
1048 final double temp3 = -3.0 * damp * expdamp * rr2;
1049 final double temp5 = -damp;
1050 final double temp7 = -0.2 - 0.6 * damp;
1051 ddsc3x = temp3 * xr;
1052 ddsc3y = temp3 * yr;
1053 ddsc3z = temp3 * zr;
1054 ddsc5x = temp5 * ddsc3x;
1055 ddsc5y = temp5 * ddsc3y;
1056 ddsc5z = temp5 * ddsc3z;
1057 ddsc7x = temp7 * ddsc5x;
1058 ddsc7y = temp7 * ddsc5y;
1059 ddsc7z = temp7 * ddsc5z;
1060 }
1061 }
1062 double ei = polarizationPair(gradient, lambdaTerm);
1063 if (isNaN(ei) || isInfinite(ei)) {
1064 String message = format(" %s\n"
1065 + " %s\n with induced dipole: %8.3f %8.3f %8.3f\n"
1066 + " %s\n with induced dipole: %8.3f %8.3f %8.3f\n"
1067 + " The polarization energy due to atoms "
1068 + "%d and %d (%d) is %10.6f at %10.6f A.",
1069 crystal.getUnitCell(), atoms[i], uix, uiy, uiz, atoms[k], ukx, uky, ukz, i + 1, k + 1, iSymm, ei, r);
1070 throw new EnergyException(message);
1071 }
1072 inducedEnergy += ei;
1073 if (esvTerm) {
1074 int esvI = extendedSystem.getTitrationESVIndex(i);
1075 int esvK = extendedSystem.getTitrationESVIndex(k);
1076 if (extendedSystem.isTitrating(i)) {
1077 double titrdUdL = 0.0;
1078 double tautdUdL = 0.0;
1079 final double[][] titrMpole = titrationMultipole[0];
1080 final double[] titrMultipolei = titrMpole[i];
1081 setMultipoleI(titrMultipolei);
1082
1083 if (extendedSystem.isTitrating(k) && esvI == esvK) {
1084 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
1085 final double[] titrMultipolek = titrNeighborMpole[k];
1086 setMultipoleK(titrMultipolek);
1087 } else {
1088 setMultipoleK(zeropole);
1089 }
1090 titrdUdL = polarizationPair(false, false);
1091
1092 setInducedI(inducedDipolepi);
1093 setInducedK(inducedDipolepk);
1094 scaled = maskingp_local[k];
1095 scalep = maskingd_local[k];
1096 titrdUdL += polarizationPair(false, false);
1097
1098 setInducedI(inducedDipolei);
1099 setInducedK(inducedDipolek);
1100 scalep = maskingp_local[k];
1101 scaled = maskingd_local[k];
1102
1103 if (extendedSystem.isTautomerizing(i)) {
1104 final double[][] tautMpole = tautomerMultipole[0];
1105 final double[] tautMultipolei = tautMpole[i];
1106 setMultipoleI(tautMultipolei);
1107
1108 if (extendedSystem.isTautomerizing(k) && esvI == esvK) {
1109 final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
1110 final double[] tautMultipolek = tautNeighborMpole[k];
1111 setMultipoleK(tautMultipolek);
1112 } else {
1113 setMultipoleK(zeropole);
1114 }
1115 tautdUdL = polarizationPair(false, false);
1116
1117 setInducedI(inducedDipolepi);
1118 setInducedK(inducedDipolepk);
1119 scaled = maskingp_local[k];
1120 scalep = maskingd_local[k];
1121 tautdUdL += polarizationPair(false, false);
1122
1123 setInducedI(inducedDipolei);
1124 setInducedK(inducedDipolek);
1125 scalep = maskingp_local[k];
1126 scaled = maskingd_local[k];
1127 }
1128 extendedSystem.addIndElecDeriv(i, titrdUdL * electric, tautdUdL * electric);
1129
1130 setMultipoleI(globalMultipolei);
1131 setMultipoleK(globalMultipolek);
1132 }
1133
1134 if (extendedSystem.isTitrating(k) && esvI != esvK) {
1135 double titrdUdL = 0.0;
1136 double tautdUdL = 0.0;
1137 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
1138 final double[] titrMultipolek = titrNeighborMpole[k];
1139 setMultipoleK(titrMultipolek);
1140 setMultipoleI(zeropole);
1141 titrdUdL = polarizationPair(false, false);
1142
1143 setInducedI(inducedDipolepi);
1144 setInducedK(inducedDipolepk);
1145 scaled = maskingp_local[k];
1146 scalep = maskingd_local[k];
1147 titrdUdL += polarizationPair(false, false);
1148
1149 setInducedI(inducedDipolei);
1150 setInducedK(inducedDipolek);
1151 scalep = maskingp_local[k];
1152 scaled = maskingd_local[k];
1153 if (extendedSystem.isTautomerizing(k)) {
1154 final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
1155 final double[] tautMultipolek = tautNeighborMpole[k];
1156 setMultipoleK(tautMultipolek);
1157 tautdUdL = polarizationPair(false, false);
1158
1159 setInducedI(inducedDipolepi);
1160 setInducedK(inducedDipolepk);
1161 scaled = maskingp_local[k];
1162 scalep = maskingd_local[k];
1163 tautdUdL += polarizationPair(false, false);
1164
1165 setInducedI(inducedDipolei);
1166 setInducedK(inducedDipolek);
1167 scalep = maskingp_local[k];
1168 scaled = maskingd_local[k];
1169 }
1170 extendedSystem.addIndElecDeriv(k, titrdUdL * electric, tautdUdL * electric);
1171
1172 setMultipoleI(globalMultipolei);
1173 setMultipoleK(globalMultipolek);
1174 }
1175 }
1176 }
1177 }
1178 if (iSymm == 0) {
1179 removeMask(i, null, masking_local, maskingp_local, maskingd_local);
1180 }
1181 }
1182 }
1183
1184
1185
1186
1187
1188
1189 private double permanentPair(boolean gradientLocal, boolean lambdaTermLocal) {
1190 final double dixdkx = diy * dkz - diz * dky;
1191 final double dixdky = diz * dkx - dix * dkz;
1192 final double dixdkz = dix * dky - diy * dkx;
1193 final double dixrx = diy * zr - diz * yr;
1194 final double dixry = diz * xr - dix * zr;
1195 final double dixrz = dix * yr - diy * xr;
1196 final double dkxrx = dky * zr - dkz * yr;
1197 final double dkxry = dkz * xr - dkx * zr;
1198 final double dkxrz = dkx * yr - dky * xr;
1199 final double qirx = qixx * xr + qixy * yr + qixz * zr;
1200 final double qiry = qixy * xr + qiyy * yr + qiyz * zr;
1201 final double qirz = qixz * xr + qiyz * yr + qizz * zr;
1202 final double qkrx = qkxx * xr + qkxy * yr + qkxz * zr;
1203 final double qkry = qkxy * xr + qkyy * yr + qkyz * zr;
1204 final double qkrz = qkxz * xr + qkyz * yr + qkzz * zr;
1205 final double qiqkrx = qixx * qkrx + qixy * qkry + qixz * qkrz;
1206 final double qiqkry = qixy * qkrx + qiyy * qkry + qiyz * qkrz;
1207 final double qiqkrz = qixz * qkrx + qiyz * qkry + qizz * qkrz;
1208 final double qkqirx = qkxx * qirx + qkxy * qiry + qkxz * qirz;
1209 final double qkqiry = qkxy * qirx + qkyy * qiry + qkyz * qirz;
1210 final double qkqirz = qkxz * qirx + qkyz * qiry + qkzz * qirz;
1211 final double qixqkx =
1212 qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy - qiyz * qkyy - qizz * qkyz;
1213 final double qixqky =
1214 qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz - qixy * qkyz - qixz * qkzz;
1215 final double qixqkz =
1216 qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx - qiyy * qkxy - qiyz * qkxz;
1217 final double rxqirx = yr * qirz - zr * qiry;
1218 final double rxqiry = zr * qirx - xr * qirz;
1219 final double rxqirz = xr * qiry - yr * qirx;
1220 final double rxqkrx = yr * qkrz - zr * qkry;
1221 final double rxqkry = zr * qkrx - xr * qkrz;
1222 final double rxqkrz = xr * qkry - yr * qkrx;
1223 final double rxqikrx = yr * qiqkrz - zr * qiqkry;
1224 final double rxqikry = zr * qiqkrx - xr * qiqkrz;
1225 final double rxqikrz = xr * qiqkry - yr * qiqkrx;
1226 final double rxqkirx = yr * qkqirz - zr * qkqiry;
1227 final double rxqkiry = zr * qkqirx - xr * qkqirz;
1228 final double rxqkirz = xr * qkqiry - yr * qkqirx;
1229 final double qkrxqirx = qkry * qirz - qkrz * qiry;
1230 final double qkrxqiry = qkrz * qirx - qkrx * qirz;
1231 final double qkrxqirz = qkrx * qiry - qkry * qirx;
1232 final double qidkx = qixx * dkx + qixy * dky + qixz * dkz;
1233 final double qidky = qixy * dkx + qiyy * dky + qiyz * dkz;
1234 final double qidkz = qixz * dkx + qiyz * dky + qizz * dkz;
1235 final double qkdix = qkxx * dix + qkxy * diy + qkxz * diz;
1236 final double qkdiy = qkxy * dix + qkyy * diy + qkyz * diz;
1237 final double qkdiz = qkxz * dix + qkyz * diy + qkzz * diz;
1238 final double dixqkrx = diy * qkrz - diz * qkry;
1239 final double dixqkry = diz * qkrx - dix * qkrz;
1240 final double dixqkrz = dix * qkry - diy * qkrx;
1241 final double dkxqirx = dky * qirz - dkz * qiry;
1242 final double dkxqiry = dkz * qirx - dkx * qirz;
1243 final double dkxqirz = dkx * qiry - dky * qirx;
1244 final double rxqidkx = yr * qidkz - zr * qidky;
1245 final double rxqidky = zr * qidkx - xr * qidkz;
1246 final double rxqidkz = xr * qidky - yr * qidkx;
1247 final double rxqkdix = yr * qkdiz - zr * qkdiy;
1248 final double rxqkdiy = zr * qkdix - xr * qkdiz;
1249 final double rxqkdiz = xr * qkdiy - yr * qkdix;
1250
1251
1252 final double sc2 = dix * dkx + diy * dky + diz * dkz;
1253 final double sc3 = dix * xr + diy * yr + diz * zr;
1254 final double sc4 = dkx * xr + dky * yr + dkz * zr;
1255 final double sc5 = qirx * xr + qiry * yr + qirz * zr;
1256 final double sc6 = qkrx * xr + qkry * yr + qkrz * zr;
1257 final double sc7 = qirx * dkx + qiry * dky + qirz * dkz;
1258 final double sc8 = qkrx * dix + qkry * diy + qkrz * diz;
1259 final double sc9 = qirx * qkrx + qiry * qkry + qirz * qkrz;
1260 final double sc10 =
1261 2.0 * (qixy * qkxy + qixz * qkxz + qiyz * qkyz) + qixx * qkxx + qiyy * qkyy + qizz * qkzz;
1262
1263
1264 final double gl0 = ci * ck;
1265 final double gl1 = ck * sc3 - ci * sc4;
1266 final double gl2 = ci * sc6 + ck * sc5 - sc3 * sc4;
1267 final double gl3 = sc3 * sc6 - sc4 * sc5;
1268 final double gl4 = sc5 * sc6;
1269 final double gl5 = -4.0 * sc9;
1270 final double gl6 = sc2;
1271 final double gl7 = 2.0 * (sc7 - sc8);
1272 final double gl8 = 2.0 * sc10;
1273
1274
1275 final double scale1 = 1.0 - scale;
1276 final double ereal =
1277 gl0 * bn0 + (gl1 + gl6) * bn1 + (gl2 + gl7 + gl8) * bn2 + (gl3 + gl5) * bn3 + gl4 * bn4;
1278 final double efix =
1279 scale1
1280 * (gl0 * rr1
1281 + (gl1 + gl6) * rr3
1282 + (gl2 + gl7 + gl8) * rr5
1283 + (gl3 + gl5) * rr7
1284 + gl4 * rr9);
1285 final double e = selfScale * l2 * (ereal - efix);
1286 if (gradientLocal) {
1287 final double gf1 =
1288 bn1 * gl0 + bn2 * (gl1 + gl6) + bn3 * (gl2 + gl7 + gl8) + bn4 * (gl3 + gl5) + bn5 * gl4;
1289 final double gf2 = -ck * bn1 + sc4 * bn2 - sc6 * bn3;
1290 final double gf3 = ci * bn1 + sc3 * bn2 + sc5 * bn3;
1291 final double gf4 = 2.0 * bn2;
1292 final double gf5 = 2.0 * (-ck * bn2 + sc4 * bn3 - sc6 * bn4);
1293 final double gf6 = 2.0 * (-ci * bn2 - sc3 * bn3 - sc5 * bn4);
1294 final double gf7 = 4.0 * bn3;
1295
1296
1297 double ftm2x =
1298 gf1 * xr
1299 + gf2 * dix
1300 + gf3 * dkx
1301 + gf4 * (qkdix - qidkx)
1302 + gf5 * qirx
1303 + gf6 * qkrx
1304 + gf7 * (qiqkrx + qkqirx);
1305 double ftm2y =
1306 gf1 * yr
1307 + gf2 * diy
1308 + gf3 * dky
1309 + gf4 * (qkdiy - qidky)
1310 + gf5 * qiry
1311 + gf6 * qkry
1312 + gf7 * (qiqkry + qkqiry);
1313 double ftm2z =
1314 gf1 * zr
1315 + gf2 * diz
1316 + gf3 * dkz
1317 + gf4 * (qkdiz - qidkz)
1318 + gf5 * qirz
1319 + gf6 * qkrz
1320 + gf7 * (qiqkrz + qkqirz);
1321
1322
1323 double ttm2x =
1324 -bn1 * dixdkx
1325 + gf2 * dixrx
1326 + gf4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1327 - gf5 * rxqirx
1328 - gf7 * (rxqikrx + qkrxqirx);
1329 double ttm2y =
1330 -bn1 * dixdky
1331 + gf2 * dixry
1332 + gf4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1333 - gf5 * rxqiry
1334 - gf7 * (rxqikry + qkrxqiry);
1335 double ttm2z =
1336 -bn1 * dixdkz
1337 + gf2 * dixrz
1338 + gf4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1339 - gf5 * rxqirz
1340 - gf7 * (rxqikrz + qkrxqirz);
1341 double ttm3x =
1342 bn1 * dixdkx
1343 + gf3 * dkxrx
1344 - gf4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1345 - gf6 * rxqkrx
1346 - gf7 * (rxqkirx - qkrxqirx);
1347 double ttm3y =
1348 bn1 * dixdky
1349 + gf3 * dkxry
1350 - gf4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1351 - gf6 * rxqkry
1352 - gf7 * (rxqkiry - qkrxqiry);
1353 double ttm3z =
1354 bn1 * dixdkz
1355 + gf3 * dkxrz
1356 - gf4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1357 - gf6 * rxqkrz
1358 - gf7 * (rxqkirz - qkrxqirz);
1359
1360
1361 if (scale1 != 0.0) {
1362 final double gfr1 =
1363 rr3 * gl0
1364 + rr5 * (gl1 + gl6)
1365 + rr7 * (gl2 + gl7 + gl8)
1366 + rr9 * (gl3 + gl5)
1367 + rr11 * gl4;
1368 final double gfr2 = -ck * rr3 + sc4 * rr5 - sc6 * rr7;
1369 final double gfr3 = ci * rr3 + sc3 * rr5 + sc5 * rr7;
1370 final double gfr4 = 2.0 * rr5;
1371 final double gfr5 = 2.0 * (-ck * rr5 + sc4 * rr7 - sc6 * rr9);
1372 final double gfr6 = 2.0 * (-ci * rr5 - sc3 * rr7 - sc5 * rr9);
1373 final double gfr7 = 4.0 * rr7;
1374
1375
1376 final double ftm2rx =
1377 gfr1 * xr
1378 + gfr2 * dix
1379 + gfr3 * dkx
1380 + gfr4 * (qkdix - qidkx)
1381 + gfr5 * qirx
1382 + gfr6 * qkrx
1383 + gfr7 * (qiqkrx + qkqirx);
1384 final double ftm2ry =
1385 gfr1 * yr
1386 + gfr2 * diy
1387 + gfr3 * dky
1388 + gfr4 * (qkdiy - qidky)
1389 + gfr5 * qiry
1390 + gfr6 * qkry
1391 + gfr7 * (qiqkry + qkqiry);
1392 final double ftm2rz =
1393 gfr1 * zr
1394 + gfr2 * diz
1395 + gfr3 * dkz
1396 + gfr4 * (qkdiz - qidkz)
1397 + gfr5 * qirz
1398 + gfr6 * qkrz
1399 + gfr7 * (qiqkrz + qkqirz);
1400
1401
1402 final double ttm2rx =
1403 -rr3 * dixdkx
1404 + gfr2 * dixrx
1405 + gfr4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1406 - gfr5 * rxqirx
1407 - gfr7 * (rxqikrx + qkrxqirx);
1408 final double ttm2ry =
1409 -rr3 * dixdky
1410 + gfr2 * dixry
1411 + gfr4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1412 - gfr5 * rxqiry
1413 - gfr7 * (rxqikry + qkrxqiry);
1414 final double ttm2rz =
1415 -rr3 * dixdkz
1416 + gfr2 * dixrz
1417 + gfr4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1418 - gfr5 * rxqirz
1419 - gfr7 * (rxqikrz + qkrxqirz);
1420 final double ttm3rx =
1421 rr3 * dixdkx
1422 + gfr3 * dkxrx
1423 - gfr4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1424 - gfr6 * rxqkrx
1425 - gfr7 * (rxqkirx - qkrxqirx);
1426 final double ttm3ry =
1427 rr3 * dixdky
1428 + gfr3 * dkxry
1429 - gfr4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1430 - gfr6 * rxqkry
1431 - gfr7 * (rxqkiry - qkrxqiry);
1432 final double ttm3rz =
1433 rr3 * dixdkz
1434 + gfr3 * dkxrz
1435 - gfr4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1436 - gfr6 * rxqkrz
1437 - gfr7 * (rxqkirz - qkrxqirz);
1438 ftm2x -= scale1 * ftm2rx;
1439 ftm2y -= scale1 * ftm2ry;
1440 ftm2z -= scale1 * ftm2rz;
1441 ttm2x -= scale1 * ttm2rx;
1442 ttm2y -= scale1 * ttm2ry;
1443 ttm2z -= scale1 * ttm2rz;
1444 ttm3x -= scale1 * ttm3rx;
1445 ttm3y -= scale1 * ttm3ry;
1446 ttm3z -= scale1 * ttm3rz;
1447 }
1448 double prefactor = electric * selfScale * l2;
1449 grad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
1450 torque.add(threadID, i, prefactor * ttm2x, prefactor * ttm2y, prefactor * ttm2z);
1451
1452
1453
1454
1455
1456 gxk_local[k] -= prefactor * ftm2x;
1457 gyk_local[k] -= prefactor * ftm2y;
1458 gzk_local[k] -= prefactor * ftm2z;
1459 txk_local[k] += prefactor * ttm3x;
1460 tyk_local[k] += prefactor * ttm3y;
1461 tzk_local[k] += prefactor * ttm3z;
1462
1463
1464 if (lambdaTerm && soft) {
1465 prefactor = electric * selfScale * dEdLSign * dlPowPerm;
1466 lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
1467 lambdaTorque.add(threadID, i, prefactor * ttm2x, prefactor * ttm2y, prefactor * ttm2z);
1468 lxk_local[k] -= prefactor * ftm2x;
1469 lyk_local[k] -= prefactor * ftm2y;
1470 lzk_local[k] -= prefactor * ftm2z;
1471 ltxk_local[k] += prefactor * ttm3x;
1472 ltyk_local[k] += prefactor * ttm3y;
1473 ltzk_local[k] += prefactor * ttm3z;
1474 }
1475 }
1476 if (lambdaTermLocal && soft) {
1477 double dRealdL =
1478 gl0 * bn1 + (gl1 + gl6) * bn2 + (gl2 + gl7 + gl8) * bn3 + (gl3 + gl5) * bn4 + gl4 * bn5;
1479 double d2RealdL2 =
1480 gl0 * bn2 + (gl1 + gl6) * bn3 + (gl2 + gl7 + gl8) * bn4 + (gl3 + gl5) * bn5 + gl4 * bn6;
1481
1482 dUdL += selfScale * (dEdLSign * dlPowPerm * ereal + l2 * dlAlpha * dRealdL);
1483 d2UdL2 +=
1484 selfScale
1485 * (dEdLSign
1486 * (d2lPowPerm * ereal
1487 + dlPowPerm * dlAlpha * dRealdL
1488 + dlPowPerm * dlAlpha * dRealdL)
1489 + l2 * d2lAlpha * dRealdL
1490 + l2 * dlAlpha * dlAlpha * d2RealdL2);
1491
1492 double dFixdL =
1493 gl0 * rr3
1494 + (gl1 + gl6) * rr5
1495 + (gl2 + gl7 + gl8) * rr7
1496 + (gl3 + gl5) * rr9
1497 + gl4 * rr11;
1498 double d2FixdL2 =
1499 gl0 * rr5
1500 + (gl1 + gl6) * rr7
1501 + (gl2 + gl7 + gl8) * rr9
1502 + (gl3 + gl5) * rr11
1503 + gl4 * rr13;
1504 dFixdL *= scale1;
1505 d2FixdL2 *= scale1;
1506 dUdL -= selfScale * (dEdLSign * dlPowPerm * efix + l2 * dlAlpha * dFixdL);
1507 d2UdL2 -=
1508 selfScale
1509 * (dEdLSign
1510 * (d2lPowPerm * efix
1511 + dlPowPerm * dlAlpha * dFixdL
1512 + dlPowPerm * dlAlpha * dFixdL)
1513 + l2 * d2lAlpha * dFixdL
1514 + l2 * dlAlpha * dlAlpha * d2FixdL2);
1515
1516
1517 final double gf1 =
1518 bn2 * gl0 + bn3 * (gl1 + gl6) + bn4 * (gl2 + gl7 + gl8) + bn5 * (gl3 + gl5) + bn6 * gl4;
1519 final double gf2 = -ck * bn2 + sc4 * bn3 - sc6 * bn4;
1520 final double gf3 = ci * bn2 + sc3 * bn3 + sc5 * bn4;
1521 final double gf4 = 2.0 * bn3;
1522 final double gf5 = 2.0 * (-ck * bn3 + sc4 * bn4 - sc6 * bn5);
1523 final double gf6 = 2.0 * (-ci * bn3 - sc3 * bn4 - sc5 * bn5);
1524 final double gf7 = 4.0 * bn4;
1525
1526
1527 double ftm2x =
1528 gf1 * xr
1529 + gf2 * dix
1530 + gf3 * dkx
1531 + gf4 * (qkdix - qidkx)
1532 + gf5 * qirx
1533 + gf6 * qkrx
1534 + gf7 * (qiqkrx + qkqirx);
1535 double ftm2y =
1536 gf1 * yr
1537 + gf2 * diy
1538 + gf3 * dky
1539 + gf4 * (qkdiy - qidky)
1540 + gf5 * qiry
1541 + gf6 * qkry
1542 + gf7 * (qiqkry + qkqiry);
1543 double ftm2z =
1544 gf1 * zr
1545 + gf2 * diz
1546 + gf3 * dkz
1547 + gf4 * (qkdiz - qidkz)
1548 + gf5 * qirz
1549 + gf6 * qkrz
1550 + gf7 * (qiqkrz + qkqirz);
1551
1552
1553 double ttm2x =
1554 -bn2 * dixdkx
1555 + gf2 * dixrx
1556 + gf4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1557 - gf5 * rxqirx
1558 - gf7 * (rxqikrx + qkrxqirx);
1559 double ttm2y =
1560 -bn2 * dixdky
1561 + gf2 * dixry
1562 + gf4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1563 - gf5 * rxqiry
1564 - gf7 * (rxqikry + qkrxqiry);
1565 double ttm2z =
1566 -bn2 * dixdkz
1567 + gf2 * dixrz
1568 + gf4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1569 - gf5 * rxqirz
1570 - gf7 * (rxqikrz + qkrxqirz);
1571 double ttm3x =
1572 bn2 * dixdkx
1573 + gf3 * dkxrx
1574 - gf4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1575 - gf6 * rxqkrx
1576 - gf7 * (rxqkirx - qkrxqirx);
1577 double ttm3y =
1578 bn2 * dixdky
1579 + gf3 * dkxry
1580 - gf4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1581 - gf6 * rxqkry
1582 - gf7 * (rxqkiry - qkrxqiry);
1583 double ttm3z =
1584 bn2 * dixdkz
1585 + gf3 * dkxrz
1586 - gf4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1587 - gf6 * rxqkrz
1588 - gf7 * (rxqkirz - qkrxqirz);
1589
1590
1591 if (scale1 != 0.0) {
1592 final double gfr1 =
1593 rr5 * gl0
1594 + rr7 * (gl1 + gl6)
1595 + rr9 * (gl2 + gl7 + gl8)
1596 + rr11 * (gl3 + gl5)
1597 + rr13 * gl4;
1598 final double gfr2 = -ck * rr5 + sc4 * rr7 - sc6 * rr9;
1599 final double gfr3 = ci * rr5 + sc3 * rr7 + sc5 * rr9;
1600 final double gfr4 = 2.0 * rr7;
1601 final double gfr5 = 2.0 * (-ck * rr7 + sc4 * rr9 - sc6 * rr11);
1602 final double gfr6 = 2.0 * (-ci * rr7 - sc3 * rr9 - sc5 * rr11);
1603 final double gfr7 = 4.0 * rr9;
1604
1605
1606 final double ftm2rx =
1607 gfr1 * xr
1608 + gfr2 * dix
1609 + gfr3 * dkx
1610 + gfr4 * (qkdix - qidkx)
1611 + gfr5 * qirx
1612 + gfr6 * qkrx
1613 + gfr7 * (qiqkrx + qkqirx);
1614 final double ftm2ry =
1615 gfr1 * yr
1616 + gfr2 * diy
1617 + gfr3 * dky
1618 + gfr4 * (qkdiy - qidky)
1619 + gfr5 * qiry
1620 + gfr6 * qkry
1621 + gfr7 * (qiqkry + qkqiry);
1622 final double ftm2rz =
1623 gfr1 * zr
1624 + gfr2 * diz
1625 + gfr3 * dkz
1626 + gfr4 * (qkdiz - qidkz)
1627 + gfr5 * qirz
1628 + gfr6 * qkrz
1629 + gfr7 * (qiqkrz + qkqirz);
1630
1631
1632 final double ttm2rx =
1633 -rr5 * dixdkx
1634 + gfr2 * dixrx
1635 + gfr4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1636 - gfr5 * rxqirx
1637 - gfr7 * (rxqikrx + qkrxqirx);
1638 final double ttm2ry =
1639 -rr5 * dixdky
1640 + gfr2 * dixry
1641 + gfr4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1642 - gfr5 * rxqiry
1643 - gfr7 * (rxqikry + qkrxqiry);
1644 final double ttm2rz =
1645 -rr5 * dixdkz
1646 + gfr2 * dixrz
1647 + gfr4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1648 - gfr5 * rxqirz
1649 - gfr7 * (rxqikrz + qkrxqirz);
1650 final double ttm3rx =
1651 rr5 * dixdkx
1652 + gfr3 * dkxrx
1653 - gfr4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1654 - gfr6 * rxqkrx
1655 - gfr7 * (rxqkirx - qkrxqirx);
1656 final double ttm3ry =
1657 rr5 * dixdky
1658 + gfr3 * dkxry
1659 - gfr4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1660 - gfr6 * rxqkry
1661 - gfr7 * (rxqkiry - qkrxqiry);
1662 final double ttm3rz =
1663 rr5 * dixdkz
1664 + gfr3 * dkxrz
1665 - gfr4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1666 - gfr6 * rxqkrz
1667 - gfr7 * (rxqkirz - qkrxqirz);
1668 ftm2x -= scale1 * ftm2rx;
1669 ftm2y -= scale1 * ftm2ry;
1670 ftm2z -= scale1 * ftm2rz;
1671 ttm2x -= scale1 * ttm2rx;
1672 ttm2y -= scale1 * ttm2ry;
1673 ttm2z -= scale1 * ttm2rz;
1674 ttm3x -= scale1 * ttm3rx;
1675 ttm3y -= scale1 * ttm3ry;
1676 ttm3z -= scale1 * ttm3rz;
1677 }
1678
1679
1680 double prefactor = electric * selfScale * l2 * dlAlpha;
1681 lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
1682 lambdaTorque.add(threadID, i, prefactor * ttm2x, prefactor * ttm2y, prefactor * ttm2z);
1683 lxk_local[k] -= prefactor * ftm2x;
1684 lyk_local[k] -= prefactor * ftm2y;
1685 lzk_local[k] -= prefactor * ftm2z;
1686 ltxk_local[k] += prefactor * ttm3x;
1687 ltyk_local[k] += prefactor * ttm3y;
1688 ltzk_local[k] += prefactor * ttm3z;
1689 }
1690
1691 return e;
1692 }
1693
1694
1695
1696
1697
1698
1699 private double polarizationPair(boolean gradientLocal, boolean lambdaTermLocal) {
1700 final double dsc3 = 1.0 - scale3 * scaled;
1701 final double dsc5 = 1.0 - scale5 * scaled;
1702 final double dsc7 = 1.0 - scale7 * scaled;
1703 final double psc3 = 1.0 - scale3 * scalep;
1704 final double psc5 = 1.0 - scale5 * scalep;
1705 final double psc7 = 1.0 - scale7 * scalep;
1706 final double usc3 = 1.0 - scale3;
1707 final double usc5 = 1.0 - scale5;
1708
1709 final double dixukx = diy * ukz - diz * uky;
1710 final double dixuky = diz * ukx - dix * ukz;
1711 final double dixukz = dix * uky - diy * ukx;
1712 final double dkxuix = dky * uiz - dkz * uiy;
1713 final double dkxuiy = dkz * uix - dkx * uiz;
1714 final double dkxuiz = dkx * uiy - dky * uix;
1715 final double dixukpx = diy * pkz - diz * pky;
1716 final double dixukpy = diz * pkx - dix * pkz;
1717 final double dixukpz = dix * pky - diy * pkx;
1718 final double dkxuipx = dky * piz - dkz * piy;
1719 final double dkxuipy = dkz * pix - dkx * piz;
1720 final double dkxuipz = dkx * piy - dky * pix;
1721 final double dixrx = diy * zr - diz * yr;
1722 final double dixry = diz * xr - dix * zr;
1723 final double dixrz = dix * yr - diy * xr;
1724 final double dkxrx = dky * zr - dkz * yr;
1725 final double dkxry = dkz * xr - dkx * zr;
1726 final double dkxrz = dkx * yr - dky * xr;
1727 final double qirx = qixx * xr + qixy * yr + qixz * zr;
1728 final double qiry = qixy * xr + qiyy * yr + qiyz * zr;
1729 final double qirz = qixz * xr + qiyz * yr + qizz * zr;
1730 final double qkrx = qkxx * xr + qkxy * yr + qkxz * zr;
1731 final double qkry = qkxy * xr + qkyy * yr + qkyz * zr;
1732 final double qkrz = qkxz * xr + qkyz * yr + qkzz * zr;
1733 final double rxqirx = yr * qirz - zr * qiry;
1734 final double rxqiry = zr * qirx - xr * qirz;
1735 final double rxqirz = xr * qiry - yr * qirx;
1736 final double rxqkrx = yr * qkrz - zr * qkry;
1737 final double rxqkry = zr * qkrx - xr * qkrz;
1738 final double rxqkrz = xr * qkry - yr * qkrx;
1739 final double qiukx = qixx * ukx + qixy * uky + qixz * ukz;
1740 final double qiuky = qixy * ukx + qiyy * uky + qiyz * ukz;
1741 final double qiukz = qixz * ukx + qiyz * uky + qizz * ukz;
1742 final double qkuix = qkxx * uix + qkxy * uiy + qkxz * uiz;
1743 final double qkuiy = qkxy * uix + qkyy * uiy + qkyz * uiz;
1744 final double qkuiz = qkxz * uix + qkyz * uiy + qkzz * uiz;
1745 final double qiukpx = qixx * pkx + qixy * pky + qixz * pkz;
1746 final double qiukpy = qixy * pkx + qiyy * pky + qiyz * pkz;
1747 final double qiukpz = qixz * pkx + qiyz * pky + qizz * pkz;
1748 final double qkuipx = qkxx * pix + qkxy * piy + qkxz * piz;
1749 final double qkuipy = qkxy * pix + qkyy * piy + qkyz * piz;
1750 final double qkuipz = qkxz * pix + qkyz * piy + qkzz * piz;
1751 final double uixqkrx = uiy * qkrz - uiz * qkry;
1752 final double uixqkry = uiz * qkrx - uix * qkrz;
1753 final double uixqkrz = uix * qkry - uiy * qkrx;
1754 final double ukxqirx = uky * qirz - ukz * qiry;
1755 final double ukxqiry = ukz * qirx - ukx * qirz;
1756 final double ukxqirz = ukx * qiry - uky * qirx;
1757 final double uixqkrpx = piy * qkrz - piz * qkry;
1758 final double uixqkrpy = piz * qkrx - pix * qkrz;
1759 final double uixqkrpz = pix * qkry - piy * qkrx;
1760 final double ukxqirpx = pky * qirz - pkz * qiry;
1761 final double ukxqirpy = pkz * qirx - pkx * qirz;
1762 final double ukxqirpz = pkx * qiry - pky * qirx;
1763 final double rxqiukx = yr * qiukz - zr * qiuky;
1764 final double rxqiuky = zr * qiukx - xr * qiukz;
1765 final double rxqiukz = xr * qiuky - yr * qiukx;
1766 final double rxqkuix = yr * qkuiz - zr * qkuiy;
1767 final double rxqkuiy = zr * qkuix - xr * qkuiz;
1768 final double rxqkuiz = xr * qkuiy - yr * qkuix;
1769 final double rxqiukpx = yr * qiukpz - zr * qiukpy;
1770 final double rxqiukpy = zr * qiukpx - xr * qiukpz;
1771 final double rxqiukpz = xr * qiukpy - yr * qiukpx;
1772 final double rxqkuipx = yr * qkuipz - zr * qkuipy;
1773 final double rxqkuipy = zr * qkuipx - xr * qkuipz;
1774 final double rxqkuipz = xr * qkuipy - yr * qkuipx;
1775
1776
1777 final double sc3 = dix * xr + diy * yr + diz * zr;
1778 final double sc4 = dkx * xr + dky * yr + dkz * zr;
1779 final double sc5 = qirx * xr + qiry * yr + qirz * zr;
1780 final double sc6 = qkrx * xr + qkry * yr + qkrz * zr;
1781
1782
1783 final double sci1 = uix * dkx + uiy * dky + uiz * dkz + dix * ukx + diy * uky + diz * ukz;
1784 final double sci3 = uix * xr + uiy * yr + uiz * zr;
1785 final double sci4 = ukx * xr + uky * yr + ukz * zr;
1786 final double sci7 = qirx * ukx + qiry * uky + qirz * ukz;
1787 final double sci8 = qkrx * uix + qkry * uiy + qkrz * uiz;
1788 final double scip1 = pix * dkx + piy * dky + piz * dkz + dix * pkx + diy * pky + diz * pkz;
1789 final double scip2 = uix * pkx + uiy * pky + uiz * pkz + pix * ukx + piy * uky + piz * ukz;
1790 final double scip3 = pix * xr + piy * yr + piz * zr;
1791 final double scip4 = pkx * xr + pky * yr + pkz * zr;
1792 final double scip7 = qirx * pkx + qiry * pky + qirz * pkz;
1793 final double scip8 = qkrx * pix + qkry * piy + qkrz * piz;
1794
1795
1796 final double gli1 = ck * sci3 - ci * sci4;
1797 final double gli2 = -sc3 * sci4 - sci3 * sc4;
1798 final double gli3 = sci3 * sc6 - sci4 * sc5;
1799 final double gli6 = sci1;
1800 final double gli7 = 2.0 * (sci7 - sci8);
1801 final double glip1 = ck * scip3 - ci * scip4;
1802 final double glip2 = -sc3 * scip4 - scip3 * sc4;
1803 final double glip3 = scip3 * sc6 - scip4 * sc5;
1804 final double glip6 = scip1;
1805 final double glip7 = 2.0 * (scip7 - scip8);
1806
1807
1808 final double ereal = (gli1 + gli6) * bn1 + (gli2 + gli7) * bn2 + gli3 * bn3;
1809 final double efix =
1810 (gli1 + gli6) * rr3 * psc3 + (gli2 + gli7) * rr5 * psc5 + gli3 * rr7 * psc7;
1811 final double e = selfScale * 0.5 * (ereal - efix);
1812
1813
1814
1815
1816 if (!(gradientLocal || lambdaTermLocal)) {
1817 return polarizationScale * e;
1818 }
1819 boolean dorli = false;
1820 if (psc3 != 0.0 || dsc3 != 0.0 || usc3 != 0.0) {
1821 dorli = true;
1822 }
1823
1824
1825 final double gfi1 =
1826 0.5 * bn2 * (gli1 + glip1 + gli6 + glip6)
1827 + 0.5 * bn2 * scip2
1828 + 0.5 * bn3 * (gli2 + glip2 + gli7 + glip7)
1829 - 0.5 * bn3 * (sci3 * scip4 + scip3 * sci4)
1830 + 0.5 * bn4 * (gli3 + glip3);
1831 final double gfi2 = -ck * bn1 + sc4 * bn2 - sc6 * bn3;
1832 final double gfi3 = ci * bn1 + sc3 * bn2 + sc5 * bn3;
1833 final double gfi4 = 2.0 * bn2;
1834 final double gfi5 = bn3 * (sci4 + scip4);
1835 final double gfi6 = -bn3 * (sci3 + scip3);
1836 double ftm2ix =
1837 gfi1 * xr
1838 + 0.5
1839 * (gfi2 * (uix + pix)
1840 + bn2 * (sci4 * pix + scip4 * uix)
1841 + gfi3 * (ukx + pkx)
1842 + bn2 * (sci3 * pkx + scip3 * ukx)
1843 + (sci4 + scip4) * bn2 * dix
1844 + (sci3 + scip3) * bn2 * dkx
1845 + gfi4 * (qkuix + qkuipx - qiukx - qiukpx))
1846 + gfi5 * qirx
1847 + gfi6 * qkrx;
1848 double ftm2iy =
1849 gfi1 * yr
1850 + 0.5
1851 * (gfi2 * (uiy + piy)
1852 + bn2 * (sci4 * piy + scip4 * uiy)
1853 + gfi3 * (uky + pky)
1854 + bn2 * (sci3 * pky + scip3 * uky)
1855 + (sci4 + scip4) * bn2 * diy
1856 + (sci3 + scip3) * bn2 * dky
1857 + gfi4 * (qkuiy + qkuipy - qiuky - qiukpy))
1858 + gfi5 * qiry
1859 + gfi6 * qkry;
1860 double ftm2iz =
1861 gfi1 * zr
1862 + 0.5
1863 * (gfi2 * (uiz + piz)
1864 + bn2 * (sci4 * piz + scip4 * uiz)
1865 + gfi3 * (ukz + pkz)
1866 + bn2 * (sci3 * pkz + scip3 * ukz)
1867 + (sci4 + scip4) * bn2 * diz
1868 + (sci3 + scip3) * bn2 * dkz
1869 + gfi4 * (qkuiz + qkuipz - qiukz - qiukpz))
1870 + gfi5 * qirz
1871 + gfi6 * qkrz;
1872
1873
1874 final double gti2 = 0.5 * bn2 * (sci4 + scip4);
1875 final double gti3 = 0.5 * bn2 * (sci3 + scip3);
1876 final double gti4 = gfi4;
1877 final double gti5 = gfi5;
1878 final double gti6 = gfi6;
1879 double ttm2ix =
1880 -0.5 * bn1 * (dixukx + dixukpx)
1881 + gti2 * dixrx
1882 - gti5 * rxqirx
1883 + 0.5 * gti4 * (ukxqirx + rxqiukx + ukxqirpx + rxqiukpx);
1884 double ttm2iy =
1885 -0.5 * bn1 * (dixuky + dixukpy)
1886 + gti2 * dixry
1887 - gti5 * rxqiry
1888 + 0.5 * gti4 * (ukxqiry + rxqiuky + ukxqirpy + rxqiukpy);
1889 double ttm2iz =
1890 -0.5 * bn1 * (dixukz + dixukpz)
1891 + gti2 * dixrz
1892 - gti5 * rxqirz
1893 + 0.5 * gti4 * (ukxqirz + rxqiukz + ukxqirpz + rxqiukpz);
1894 double ttm3ix =
1895 -0.5 * bn1 * (dkxuix + dkxuipx)
1896 + gti3 * dkxrx
1897 - gti6 * rxqkrx
1898 - 0.5 * gti4 * (uixqkrx + rxqkuix + uixqkrpx + rxqkuipx);
1899 double ttm3iy =
1900 -0.5 * bn1 * (dkxuiy + dkxuipy)
1901 + gti3 * dkxry
1902 - gti6 * rxqkry
1903 - 0.5 * gti4 * (uixqkry + rxqkuiy + uixqkrpy + rxqkuipy);
1904 double ttm3iz =
1905 -0.5 * bn1 * (dkxuiz + dkxuipz)
1906 + gti3 * dkxrz
1907 - gti6 * rxqkrz
1908 - 0.5 * gti4 * (uixqkrz + rxqkuiz + uixqkrpz + rxqkuipz);
1909 double ftm2rix = 0.0;
1910 double ftm2riy = 0.0;
1911 double ftm2riz = 0.0;
1912 double ttm2rix = 0.0;
1913 double ttm2riy = 0.0;
1914 double ttm2riz = 0.0;
1915 double ttm3rix = 0.0;
1916 double ttm3riy = 0.0;
1917 double ttm3riz = 0.0;
1918 if (dorli) {
1919
1920 final double gfri1 =
1921 0.5 * rr5 * ((gli1 + gli6) * psc3 + (glip1 + glip6) * dsc3 + scip2 * usc3)
1922 + 0.5
1923 * rr7
1924 * ((gli7 + gli2) * psc5
1925 + (glip7 + glip2) * dsc5
1926 - (sci3 * scip4 + scip3 * sci4) * usc5)
1927 + 0.5 * rr9 * (gli3 * psc7 + glip3 * dsc7);
1928 final double gfri4 = 2.0 * rr5;
1929 final double gfri5 = rr7 * (sci4 * psc7 + scip4 * dsc7);
1930 final double gfri6 = -rr7 * (sci3 * psc7 + scip3 * dsc7);
1931 ftm2rix =
1932 gfri1 * xr
1933 + 0.5
1934 * (-rr3 * ck * (uix * psc3 + pix * dsc3)
1935 + rr5 * sc4 * (uix * psc5 + pix * dsc5)
1936 - rr7 * sc6 * (uix * psc7 + pix * dsc7))
1937 + (rr3 * ci * (ukx * psc3 + pkx * dsc3)
1938 + rr5 * sc3 * (ukx * psc5 + pkx * dsc5)
1939 + rr7 * sc5 * (ukx * psc7 + pkx * dsc7))
1940 * 0.5
1941 + rr5 * usc5 * (sci4 * pix + scip4 * uix + sci3 * pkx + scip3 * ukx) * 0.5
1942 + 0.5 * (sci4 * psc5 + scip4 * dsc5) * rr5 * dix
1943 + 0.5 * (sci3 * psc5 + scip3 * dsc5) * rr5 * dkx
1944 + 0.5 * gfri4 * ((qkuix - qiukx) * psc5 + (qkuipx - qiukpx) * dsc5)
1945 + gfri5 * qirx
1946 + gfri6 * qkrx;
1947 ftm2riy =
1948 gfri1 * yr
1949 + 0.5
1950 * (-rr3 * ck * (uiy * psc3 + piy * dsc3)
1951 + rr5 * sc4 * (uiy * psc5 + piy * dsc5)
1952 - rr7 * sc6 * (uiy * psc7 + piy * dsc7))
1953 + (rr3 * ci * (uky * psc3 + pky * dsc3)
1954 + rr5 * sc3 * (uky * psc5 + pky * dsc5)
1955 + rr7 * sc5 * (uky * psc7 + pky * dsc7))
1956 * 0.5
1957 + rr5 * usc5 * (sci4 * piy + scip4 * uiy + sci3 * pky + scip3 * uky) * 0.5
1958 + 0.5 * (sci4 * psc5 + scip4 * dsc5) * rr5 * diy
1959 + 0.5 * (sci3 * psc5 + scip3 * dsc5) * rr5 * dky
1960 + 0.5 * gfri4 * ((qkuiy - qiuky) * psc5 + (qkuipy - qiukpy) * dsc5)
1961 + gfri5 * qiry
1962 + gfri6 * qkry;
1963 ftm2riz =
1964 gfri1 * zr
1965 + 0.5
1966 * (-rr3 * ck * (uiz * psc3 + piz * dsc3)
1967 + rr5 * sc4 * (uiz * psc5 + piz * dsc5)
1968 - rr7 * sc6 * (uiz * psc7 + piz * dsc7))
1969 + (rr3 * ci * (ukz * psc3 + pkz * dsc3)
1970 + rr5 * sc3 * (ukz * psc5 + pkz * dsc5)
1971 + rr7 * sc5 * (ukz * psc7 + pkz * dsc7))
1972 * 0.5
1973 + rr5 * usc5 * (sci4 * piz + scip4 * uiz + sci3 * pkz + scip3 * ukz) * 0.5
1974 + 0.5 * (sci4 * psc5 + scip4 * dsc5) * rr5 * diz
1975 + 0.5 * (sci3 * psc5 + scip3 * dsc5) * rr5 * dkz
1976 + 0.5 * gfri4 * ((qkuiz - qiukz) * psc5 + (qkuipz - qiukpz) * dsc5)
1977 + gfri5 * qirz
1978 + gfri6 * qkrz;
1979
1980
1981 final double gtri2 = 0.5 * rr5 * (sci4 * psc5 + scip4 * dsc5);
1982 final double gtri3 = 0.5 * rr5 * (sci3 * psc5 + scip3 * dsc5);
1983 final double gtri4 = gfri4;
1984 final double gtri5 = gfri5;
1985 final double gtri6 = gfri6;
1986 ttm2rix =
1987 -rr3 * (dixukx * psc3 + dixukpx * dsc3) * 0.5
1988 + gtri2 * dixrx
1989 - gtri5 * rxqirx
1990 + gtri4 * ((ukxqirx + rxqiukx) * psc5 + (ukxqirpx + rxqiukpx) * dsc5) * 0.5;
1991 ttm2riy =
1992 -rr3 * (dixuky * psc3 + dixukpy * dsc3) * 0.5
1993 + gtri2 * dixry
1994 - gtri5 * rxqiry
1995 + gtri4 * ((ukxqiry + rxqiuky) * psc5 + (ukxqirpy + rxqiukpy) * dsc5) * 0.5;
1996 ttm2riz =
1997 -rr3 * (dixukz * psc3 + dixukpz * dsc3) * 0.5
1998 + gtri2 * dixrz
1999 - gtri5 * rxqirz
2000 + gtri4 * ((ukxqirz + rxqiukz) * psc5 + (ukxqirpz + rxqiukpz) * dsc5) * 0.5;
2001 ttm3rix =
2002 -rr3 * (dkxuix * psc3 + dkxuipx * dsc3) * 0.5
2003 + gtri3 * dkxrx
2004 - gtri6 * rxqkrx
2005 - gtri4 * ((uixqkrx + rxqkuix) * psc5 + (uixqkrpx + rxqkuipx) * dsc5) * 0.5;
2006 ttm3riy =
2007 -rr3 * (dkxuiy * psc3 + dkxuipy * dsc3) * 0.5
2008 + gtri3 * dkxry
2009 - gtri6 * rxqkry
2010 - gtri4 * ((uixqkry + rxqkuiy) * psc5 + (uixqkrpy + rxqkuipy) * dsc5) * 0.5;
2011 ttm3riz =
2012 -rr3 * (dkxuiz * psc3 + dkxuipz * dsc3) * 0.5
2013 + gtri3 * dkxrz
2014 - gtri6 * rxqkrz
2015 - gtri4 * ((uixqkrz + rxqkuiz) * psc5 + (uixqkrpz + rxqkuipz) * dsc5) * 0.5;
2016 }
2017
2018
2019 double temp3 = 0.5 * rr3 * ((gli1 + gli6) * scalep + (glip1 + glip6) * scaled);
2020 double temp5 = 0.5 * rr5 * ((gli2 + gli7) * scalep + (glip2 + glip7) * scaled);
2021 final double temp7 = 0.5 * rr7 * (gli3 * scalep + glip3 * scaled);
2022 final double fridmpx = temp3 * ddsc3x + temp5 * ddsc5x + temp7 * ddsc7x;
2023 final double fridmpy = temp3 * ddsc3y + temp5 * ddsc5y + temp7 * ddsc7y;
2024 final double fridmpz = temp3 * ddsc3z + temp5 * ddsc5z + temp7 * ddsc7z;
2025
2026
2027 temp3 = 0.5 * rr3 * scip2;
2028 temp5 = -0.5 * rr5 * (sci3 * scip4 + scip3 * sci4);
2029 final double findmpx = temp3 * ddsc3x + temp5 * ddsc5x;
2030 final double findmpy = temp3 * ddsc3y + temp5 * ddsc5y;
2031 final double findmpz = temp3 * ddsc3z + temp5 * ddsc5z;
2032
2033
2034
2035
2036
2037
2038 ftm2ix = ftm2ix - fridmpx - findmpx;
2039 ftm2iy = ftm2iy - fridmpy - findmpy;
2040 ftm2iz = ftm2iz - fridmpz - findmpz;
2041
2042
2043 if (polarization == Polarization.DIRECT) {
2044 final double gfd = 0.5 * (bn2 * scip2 - bn3 * (scip3 * sci4 + sci3 * scip4));
2045 final double gfdr = 0.5 * (rr5 * scip2 * usc3 - rr7 * (scip3 * sci4 + sci3 * scip4) * usc5);
2046 ftm2ix = ftm2ix - gfd * xr - 0.5 * bn2 * (sci4 * pix + scip4 * uix + sci3 * pkx + scip3 * ukx);
2047 ftm2iy = ftm2iy - gfd * yr - 0.5 * bn2 * (sci4 * piy + scip4 * uiy + sci3 * pky + scip3 * uky);
2048 ftm2iz = ftm2iz - gfd * zr - 0.5 * bn2 * (sci4 * piz + scip4 * uiz + sci3 * pkz + scip3 * ukz);
2049 final double fdirx = gfdr * xr + 0.5 * usc5 * rr5 * (sci4 * pix + scip4 * uix + sci3 * pkx + scip3 * ukx);
2050 final double fdiry = gfdr * yr + 0.5 * usc5 * rr5 * (sci4 * piy + scip4 * uiy + sci3 * pky + scip3 * uky);
2051 final double fdirz = gfdr * zr + 0.5 * usc5 * rr5 * (sci4 * piz + scip4 * uiz + sci3 * pkz + scip3 * ukz);
2052 ftm2ix = ftm2ix + fdirx + findmpx;
2053 ftm2iy = ftm2iy + fdiry + findmpy;
2054 ftm2iz = ftm2iz + fdirz + findmpz;
2055 }
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142 ftm2ix = ftm2ix - ftm2rix;
2143 ftm2iy = ftm2iy - ftm2riy;
2144 ftm2iz = ftm2iz - ftm2riz;
2145 ttm2ix = ttm2ix - ttm2rix;
2146 ttm2iy = ttm2iy - ttm2riy;
2147 ttm2iz = ttm2iz - ttm2riz;
2148 ttm3ix = ttm3ix - ttm3rix;
2149 ttm3iy = ttm3iy - ttm3riy;
2150 ttm3iz = ttm3iz - ttm3riz;
2151
2152 double scalar = electric * polarizationScale * selfScale;
2153
2154 grad.add(threadID, i, scalar * ftm2ix, scalar * ftm2iy, scalar * ftm2iz);
2155 torque.add(threadID, i, scalar * ttm2ix, scalar * ttm2iy, scalar * ttm2iz);
2156 gxk_local[k] -= scalar * ftm2ix;
2157 gyk_local[k] -= scalar * ftm2iy;
2158 gzk_local[k] -= scalar * ftm2iz;
2159 txk_local[k] += scalar * ttm3ix;
2160 tyk_local[k] += scalar * ttm3iy;
2161 tzk_local[k] += scalar * ttm3iz;
2162 if (lambdaTermLocal) {
2163 dUdL += dEdLSign * dlPowPol * e;
2164 d2UdL2 += dEdLSign * d2lPowPol * e;
2165 scalar = electric * dEdLSign * dlPowPol * selfScale;
2166 lambdaGrad.add(threadID, i, scalar * ftm2ix, scalar * ftm2iy, scalar * ftm2iz);
2167 lambdaTorque.add(threadID, i, scalar * ttm2ix, scalar * ttm2iy, scalar * ttm2iz);
2168 lxk_local[k] -= scalar * ftm2ix;
2169 lyk_local[k] -= scalar * ftm2iy;
2170 lzk_local[k] -= scalar * ftm2iz;
2171 ltxk_local[k] += scalar * ttm3ix;
2172 ltyk_local[k] += scalar * ttm3iy;
2173 ltzk_local[k] += scalar * ttm3iz;
2174 }
2175 return polarizationScale * e;
2176 }
2177
2178 private void setMultipoleI(double[] globalMultipolei) {
2179 ci = globalMultipolei[t000];
2180 dix = globalMultipolei[t100];
2181 diy = globalMultipolei[t010];
2182 diz = globalMultipolei[t001];
2183 qixx = globalMultipolei[t200] * oneThird;
2184 qiyy = globalMultipolei[t020] * oneThird;
2185 qizz = globalMultipolei[t002] * oneThird;
2186 qixy = globalMultipolei[t110] * oneThird;
2187 qixz = globalMultipolei[t101] * oneThird;
2188 qiyz = globalMultipolei[t011] * oneThird;
2189 }
2190
2191 private void setMultipoleK(double[] globalMultipolek) {
2192 ck = globalMultipolek[t000];
2193 dkx = globalMultipolek[t100];
2194 dky = globalMultipolek[t010];
2195 dkz = globalMultipolek[t001];
2196 qkxx = globalMultipolek[t200] * oneThird;
2197 qkyy = globalMultipolek[t020] * oneThird;
2198 qkzz = globalMultipolek[t002] * oneThird;
2199 qkxy = globalMultipolek[t110] * oneThird;
2200 qkxz = globalMultipolek[t101] * oneThird;
2201 qkyz = globalMultipolek[t011] * oneThird;
2202 }
2203
2204 private void setInducedI(double[] inducedDipolei) {
2205 uix = inducedDipolei[0];
2206 uiy = inducedDipolei[1];
2207 uiz = inducedDipolei[2];
2208 }
2209
2210 private void setInducedK(double[] inducedDipolek) {
2211 ukx = inducedDipolek[0];
2212 uky = inducedDipolek[1];
2213 ukz = inducedDipolek[2];
2214 }
2215
2216 private void setInducedpI(double[] inducedDipolepi) {
2217 pix = inducedDipolepi[0];
2218 piy = inducedDipolepi[1];
2219 piz = inducedDipolepi[2];
2220 }
2221
2222 private void setInducedpK(double[] inducedDipolepk) {
2223 pkx = inducedDipolepk[0];
2224 pky = inducedDipolepk[1];
2225 pkz = inducedDipolepk[2];
2226 }
2227
2228 }
2229
2230
2231
2232
2233
2234 private class FixedChargeEnergyLoop extends IntegerForLoop {
2235
2236 private final double[] dx_local;
2237 private final double[][] rot_local;
2238 private double ci;
2239 private double ck;
2240 private double xr, yr, zr;
2241 private double bn0, bn1, bn2;
2242 private double rr1, rr3, rr5;
2243 private double scale;
2244 private double l2;
2245 private boolean soft;
2246 private double selfScale;
2247 private double permanentEnergy;
2248 private double inducedEnergy;
2249 private double dUdL, d2UdL2;
2250 private int i, k, iSymm, count;
2251 private double[] gxk_local, gyk_local, gzk_local;
2252 private double[] lxk_local, lyk_local, lzk_local;
2253 private double[] masking_local;
2254 private int threadID;
2255
2256 FixedChargeEnergyLoop() {
2257 super();
2258 dx_local = new double[3];
2259 rot_local = new double[3][3];
2260 }
2261
2262 @Override
2263 public void finish() {
2264 sharedInteractions.addAndGet(count);
2265 if (lambdaTerm) {
2266 shareddEdLambda.addAndGet(dUdL * electric);
2267 sharedd2EdLambda2.addAndGet(d2UdL2 * electric);
2268 }
2269 realSpaceEnergyTime[getThreadIndex()] += System.nanoTime();
2270 }
2271
2272 public int getCount() {
2273 return count;
2274 }
2275
2276 @Override
2277 public void run(int lb, int ub) {
2278 List<SymOp> symOps = crystal.spaceGroup.symOps;
2279 int nSymm = symOps.size();
2280 int nAtoms = atoms.length;
2281 for (iSymm = 0; iSymm < nSymm; iSymm++) {
2282 SymOp symOp = symOps.get(iSymm);
2283 if (gradient) {
2284 fill(gxk_local, 0.0);
2285 fill(gyk_local, 0.0);
2286 fill(gzk_local, 0.0);
2287 }
2288 if (lambdaTerm) {
2289 fill(lxk_local, 0.0);
2290 fill(lyk_local, 0.0);
2291 fill(lzk_local, 0.0);
2292 }
2293 realSpaceChunk(lb, ub);
2294 if (gradient) {
2295
2296 if (iSymm != 0) {
2297 crystal.applyTransSymRot(nAtoms, gxk_local, gyk_local, gzk_local, gxk_local, gyk_local, gzk_local, symOp, rot_local);
2298 }
2299
2300 for (int j = 0; j < nAtoms; j++) {
2301 grad.add(threadID, j, gxk_local[j], gyk_local[j], gzk_local[j]);
2302 }
2303 }
2304 if (lambdaTerm) {
2305
2306 if (iSymm != 0) {
2307 crystal.applyTransSymRot(
2308 nAtoms, lxk_local, lyk_local, lzk_local, lxk_local, lyk_local, lzk_local, symOp,
2309 rot_local);
2310 }
2311
2312 for (int j = 0; j < nAtoms; j++) {
2313 lambdaGrad.add(threadID, j, lxk_local[j], lyk_local[j], lzk_local[j]);
2314 }
2315 }
2316 }
2317 }
2318
2319 @Override
2320 public IntegerSchedule schedule() {
2321 return realSpaceSchedule;
2322 }
2323
2324 @Override
2325 public void start() {
2326 init();
2327 threadID = getThreadIndex();
2328 realSpaceEnergyTime[threadID] -= System.nanoTime();
2329 permanentEnergy = 0.0;
2330 inducedEnergy = 0.0;
2331 count = 0;
2332 if (lambdaTerm) {
2333 dUdL = 0.0;
2334 d2UdL2 = 0.0;
2335 }
2336 }
2337
2338 private void init() {
2339 int nAtoms = atoms.length;
2340 if (masking_local == null || masking_local.length < nAtoms) {
2341 gxk_local = new double[nAtoms];
2342 gyk_local = new double[nAtoms];
2343 gzk_local = new double[nAtoms];
2344 lxk_local = new double[nAtoms];
2345 lyk_local = new double[nAtoms];
2346 lzk_local = new double[nAtoms];
2347 masking_local = new double[nAtoms];
2348 fill(masking_local, 1.0);
2349 }
2350 }
2351
2352
2353
2354
2355
2356
2357
2358 private void realSpaceChunk(final int lb, final int ub) {
2359 final double[] x = coordinates[0][0];
2360 final double[] y = coordinates[0][1];
2361 final double[] z = coordinates[0][2];
2362 final double[][] mpole = globalMultipole[0];
2363 final int[][] lists = realSpaceLists[iSymm];
2364 final double[] neighborX = coordinates[iSymm][0];
2365 final double[] neighborY = coordinates[iSymm][1];
2366 final double[] neighborZ = coordinates[iSymm][2];
2367 final double[][] neighborMultipole = globalMultipole[iSymm];
2368 for (i = lb; i <= ub; i++) {
2369 if (!use[i]) {
2370 continue;
2371 }
2372 final int moleculei = molecule[i];
2373 if (iSymm == 0) {
2374 applyMask(i, null, masking_local);
2375 }
2376 final double xi = x[i];
2377 final double yi = y[i];
2378 final double zi = z[i];
2379 final double[] globalMultipolei = mpole[i];
2380 setMultipoleI(globalMultipolei);
2381 final boolean softi = isSoft[i];
2382 final int[] list = lists[i];
2383 final int npair = realSpaceCounts[iSymm][i];
2384 for (int j = 0; j < npair; j++) {
2385 k = list[j];
2386 if (!use[k]) {
2387 continue;
2388 }
2389 boolean sameMolecule = (moleculei == molecule[k]);
2390 if (lambdaMode == LambdaMode.VAPOR) {
2391 if ((intermolecularSoftcore && !sameMolecule)
2392 || (intramolecularSoftcore && sameMolecule)) {
2393 continue;
2394 }
2395 }
2396 selfScale = 1.0;
2397 if (i == k) {
2398 selfScale = 0.5;
2399 }
2400 double beta = 0.0;
2401 l2 = 1.0;
2402 soft = (softi || isSoft[k]);
2403 if (soft && doPermanentRealSpace) {
2404 beta = lAlpha;
2405 l2 = permanentScale;
2406 } else if (nnTerm) {
2407 l2 = permanentScale;
2408 }
2409 final double xk = neighborX[k];
2410 final double yk = neighborY[k];
2411 final double zk = neighborZ[k];
2412 dx_local[0] = xk - xi;
2413 dx_local[1] = yk - yi;
2414 dx_local[2] = zk - zi;
2415 double r2 = crystal.image(dx_local);
2416 xr = dx_local[0];
2417 yr = dx_local[1];
2418 zr = dx_local[2];
2419 final double[] globalMultipolek = neighborMultipole[k];
2420 setMultipoleK(globalMultipolek);
2421 scale = masking_local[k];
2422 double r = sqrt(r2 + beta);
2423 double ralpha = aewald * r;
2424 double exp2a = exp(-ralpha * ralpha);
2425 rr1 = 1.0 / r;
2426 double rr2 = rr1 * rr1;
2427 bn0 = erfc(ralpha) * rr1;
2428 bn1 = (bn0 + an0 * exp2a) * rr2;
2429 bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
2430 rr3 = rr1 * rr2;
2431 rr5 = 3.0 * rr3 * rr2;
2432 if (doPermanentRealSpace) {
2433 double ei = permanentPair(gradient, lambdaTerm);
2434 if (isNaN(ei) || isInfinite(ei)) {
2435 String message = format(" %s\n %s\n %s\n The permanent multipole energy between "
2436 + "atoms %d and %d (%d) is %16.8f at %16.8f A.",
2437 crystal.getUnitCell().toString(), atoms[i].toString(), atoms[k].toString(),
2438 i, k, iSymm, ei, r);
2439 throw new EnergyException(message);
2440 }
2441 if (ei != 0.0) {
2442 permanentEnergy += ei;
2443 count++;
2444 }
2445 if (esvTerm) {
2446 if (extendedSystem.isTitrating(i)) {
2447 double titrdUdL = 0.0;
2448 double tautdUdL = 0.0;
2449 final double[][] titrMpole = titrationMultipole[0];
2450 final double[] titrMultipolei = titrMpole[i];
2451 setMultipoleI(titrMultipolei);
2452 titrdUdL = electric * permanentPair(false, false);
2453 if (extendedSystem.isTautomerizing(i)) {
2454 final double[][] tautMpole = tautomerMultipole[0];
2455 final double[] tautMultipolei = tautMpole[i];
2456 setMultipoleI(tautMultipolei);
2457 tautdUdL = electric * permanentPair(false, false);
2458 }
2459 extendedSystem.addPermElecDeriv(i, titrdUdL, tautdUdL);
2460
2461 setMultipoleI(globalMultipolei);
2462 }
2463 if (extendedSystem.isTitrating(k)) {
2464 double titrdUdL = 0.0;
2465 double tautdUdL = 0.0;
2466 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
2467 final double[] titrMultipolek = titrNeighborMpole[k];
2468 setMultipoleK(titrMultipolek);
2469 titrdUdL = electric * permanentPair(false, false);
2470 if (extendedSystem.isTautomerizing(k)) {
2471 final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
2472 final double[] tautMultipolek = tautNeighborMpole[k];
2473 setMultipoleK(tautMultipolek);
2474 tautdUdL = electric * permanentPair(false, false);
2475 }
2476 extendedSystem.addPermElecDeriv(k, titrdUdL, tautdUdL);
2477
2478 setMultipoleK(globalMultipolek);
2479 }
2480 }
2481 }
2482 }
2483 if (iSymm == 0) {
2484 removeMask(i, null, masking_local);
2485 }
2486 }
2487 }
2488
2489
2490
2491
2492
2493
2494 private double permanentPair(boolean gradientLocal, boolean lambdaTermLocal) {
2495
2496
2497 final double gl0 = ci * ck;
2498
2499
2500 final double scale1 = 1.0 - scale;
2501 final double ereal = gl0 * bn0;
2502 final double efix = scale1 * (gl0 * rr1);
2503 final double e = selfScale * l2 * (ereal - efix);
2504 if (gradientLocal) {
2505 final double gf1 = bn1 * gl0;
2506
2507 double ftm2x = gf1 * xr;
2508 double ftm2y = gf1 * yr;
2509 double ftm2z = gf1 * zr;
2510
2511
2512 if (scale1 != 0.0) {
2513 final double gfr1 = rr3 * gl0;
2514
2515 final double ftm2rx = gfr1 * xr;
2516 final double ftm2ry = gfr1 * yr;
2517 final double ftm2rz = gfr1 * zr;
2518 ftm2x -= scale1 * ftm2rx;
2519 ftm2y -= scale1 * ftm2ry;
2520 ftm2z -= scale1 * ftm2rz;
2521 }
2522 double prefactor = electric * selfScale * l2;
2523 grad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
2524 gxk_local[k] -= prefactor * ftm2x;
2525 gyk_local[k] -= prefactor * ftm2y;
2526 gzk_local[k] -= prefactor * ftm2z;
2527
2528
2529 if (lambdaTerm && soft) {
2530 prefactor = electric * selfScale * dEdLSign * dlPowPerm;
2531 lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
2532 lxk_local[k] -= prefactor * ftm2x;
2533 lyk_local[k] -= prefactor * ftm2y;
2534 lzk_local[k] -= prefactor * ftm2z;
2535 }
2536 }
2537
2538 if (lambdaTermLocal && soft) {
2539 double dRealdL = gl0 * bn1;
2540 double d2RealdL2 = gl0 * bn2;
2541
2542 dUdL += selfScale * (dEdLSign * dlPowPerm * ereal + l2 * dlAlpha * dRealdL);
2543 d2UdL2 += selfScale * (dEdLSign * (d2lPowPerm * ereal
2544 + dlPowPerm * dlAlpha * dRealdL
2545 + dlPowPerm * dlAlpha * dRealdL)
2546 + l2 * d2lAlpha * dRealdL
2547 + l2 * dlAlpha * dlAlpha * d2RealdL2);
2548
2549 double dFixdL = gl0 * rr3;
2550 double d2FixdL2 = gl0 * rr5;
2551 dFixdL *= scale1;
2552 d2FixdL2 *= scale1;
2553 dUdL -= selfScale * (dEdLSign * dlPowPerm * efix + l2 * dlAlpha * dFixdL);
2554 d2UdL2 -= selfScale * (dEdLSign * (d2lPowPerm * efix
2555 + dlPowPerm * dlAlpha * dFixdL
2556 + dlPowPerm * dlAlpha * dFixdL)
2557 + l2 * d2lAlpha * dFixdL
2558 + l2 * dlAlpha * dlAlpha * d2FixdL2);
2559
2560
2561 final double gf1 = bn2 * gl0;
2562
2563 double ftm2x = gf1 * xr;
2564 double ftm2y = gf1 * yr;
2565 double ftm2z = gf1 * zr;
2566
2567
2568 if (scale1 != 0.0) {
2569 final double gfr1 = rr5 * gl0;
2570
2571 final double ftm2rx = gfr1 * xr;
2572 final double ftm2ry = gfr1 * yr;
2573 final double ftm2rz = gfr1 * zr;
2574 ftm2x -= scale1 * ftm2rx;
2575 ftm2y -= scale1 * ftm2ry;
2576 ftm2z -= scale1 * ftm2rz;
2577 }
2578
2579
2580 double prefactor = electric * selfScale * l2 * dlAlpha;
2581 lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
2582 lxk_local[k] -= prefactor * ftm2x;
2583 lyk_local[k] -= prefactor * ftm2y;
2584 lzk_local[k] -= prefactor * ftm2z;
2585 }
2586 return e;
2587 }
2588
2589 private void setMultipoleI(double[] globalMultipolei) {
2590 ci = globalMultipolei[t000];
2591 }
2592
2593 private void setMultipoleK(double[] globalMultipolek) {
2594 ck = globalMultipolek[t000];
2595 }
2596
2597 }
2598
2599 }