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.ParallelSection;
44 import edu.rit.pj.ParallelTeam;
45 import edu.rit.pj.reduction.SharedInteger;
46 import edu.rit.util.Range;
47 import ffx.crystal.Crystal;
48 import ffx.crystal.SymOp;
49 import ffx.numerics.atomic.AtomicDoubleArray3D;
50 import ffx.potential.bonded.Atom;
51 import ffx.potential.nonbonded.MaskingInterface;
52 import ffx.potential.nonbonded.ReciprocalSpace;
53 import ffx.potential.parameters.ForceField;
54
55 import java.util.List;
56 import java.util.logging.Level;
57 import java.util.logging.Logger;
58
59 import static ffx.numerics.special.Erf.erfc;
60 import static ffx.potential.parameters.MultipoleType.t000;
61 import static ffx.potential.parameters.MultipoleType.t001;
62 import static ffx.potential.parameters.MultipoleType.t002;
63 import static ffx.potential.parameters.MultipoleType.t010;
64 import static ffx.potential.parameters.MultipoleType.t011;
65 import static ffx.potential.parameters.MultipoleType.t020;
66 import static ffx.potential.parameters.MultipoleType.t100;
67 import static ffx.potential.parameters.MultipoleType.t101;
68 import static ffx.potential.parameters.MultipoleType.t110;
69 import static ffx.potential.parameters.MultipoleType.t200;
70 import static java.util.Arrays.copyOf;
71 import static java.util.Arrays.fill;
72 import static org.apache.commons.math3.util.FastMath.exp;
73 import static org.apache.commons.math3.util.FastMath.min;
74 import static org.apache.commons.math3.util.FastMath.sqrt;
75
76
77
78
79
80
81
82
83
84
85
86
87 public class PermanentFieldRegion extends ParallelRegion implements MaskingInterface {
88
89 private static final Logger logger = Logger.getLogger(PermanentFieldRegion.class.getName());
90
91
92
93 private static final double oneThird = 1.0 / 3.0;
94
95
96
97 private final boolean intermolecularSoftcore;
98
99
100
101 private final boolean intramolecularSoftcore;
102
103
104
105 public double[][][] inducedDipole;
106 public double[][][] inducedDipoleCR;
107
108
109
110 protected int[][] ip11;
111
112
113
114 private Atom[] atoms;
115
116
117
118 private Crystal crystal;
119
120
121
122 private double[][][] coordinates;
123
124
125
126 private double[][][] globalMultipole;
127
128
129
130 private int[][][] neighborLists;
131
132
133
134
135 private int[][][] preconditionerLists;
136
137
138
139 private int[][] preconditionerCounts;
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158 private boolean[] use;
159
160
161
162 private int[] molecule;
163 private double[] ipdamp;
164 private double[] thole;
165
166
167
168 private int[][] mask12;
169 private int[][] mask13;
170 private int[][] mask14;
171
172
173
174 private LambdaMode lambdaMode = LambdaMode.OFF;
175
176
177
178 private ReciprocalSpace reciprocalSpace;
179 private boolean reciprocalSpaceTerm;
180 private double off2;
181 private double preconditionerCutoff;
182 private double an0, an1, an2;
183 private double aewald;
184
185
186
187 private int[][][] realSpaceLists;
188
189
190
191 private int[][] realSpaceCounts;
192 private Range[] realSpaceRanges;
193 private IntegerSchedule permanentSchedule;
194
195
196
197 private AtomicDoubleArray3D field;
198
199
200
201 private AtomicDoubleArray3D fieldCR;
202 private ScaleParameters scaleParameters;
203 private final PermanentRealSpaceFieldSection permanentRealSpaceFieldSection;
204 private final PermanentReciprocalSection permanentReciprocalSection;
205
206
207
208 private PMETimings pmeTimings;
209
210 public PermanentFieldRegion(ParallelTeam pt, ForceField forceField, boolean lambdaTerm) {
211 permanentRealSpaceFieldSection = new PermanentRealSpaceFieldSection(pt);
212 permanentReciprocalSection = new PermanentReciprocalSection();
213
214
215 if (lambdaTerm) {
216 intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
217 intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
218 } else {
219 intermolecularSoftcore = false;
220 intramolecularSoftcore = false;
221 }
222 }
223
224
225
226
227
228
229
230
231 @Override
232 public void applyMask(int i, boolean[] is14, double[]... masks) {
233 if (ip11[i] != null) {
234 double[] energyMask = masks[0];
235 var m12 = mask12[i];
236 for (int value : m12) {
237 energyMask[value] = scaleParameters.p12scale;
238 }
239 var m13 = mask13[i];
240 for (int value : m13) {
241 energyMask[value] = scaleParameters.p13scale;
242 }
243 var m14 = mask14[i];
244 for (int value : m14) {
245 energyMask[value] = scaleParameters.p14scale;
246 for (int k : ip11[i]) {
247 if (k == value) {
248 energyMask[value] = scaleParameters.intra14Scale * scaleParameters.p14scale;
249 break;
250 }
251 }
252 }
253
254 double[] inductionMask = masks[1];
255 for (int index : ip11[i]) {
256 inductionMask[index] = scaleParameters.d11scale;
257 }
258 }
259 }
260
261 public void initTimings() {
262 permanentRealSpaceFieldSection.initTimings();
263 }
264
265 public long getRealSpacePermTime() {
266 return permanentRealSpaceFieldSection.time;
267 }
268
269 public long getInitTime(int threadId) {
270 return permanentRealSpaceFieldSection.permanentRealSpaceFieldRegion.initializationLoop[threadId].time;
271 }
272
273 public long getPermTime(int threadId) {
274 return permanentRealSpaceFieldSection.permanentRealSpaceFieldRegion.permanentRealSpaceFieldLoop[threadId].time;
275 }
276
277 public void init(
278 Atom[] atoms,
279 Crystal crystal,
280 double[][][] coordinates,
281 double[][][] globalMultipole,
282 double[][][] inducedDipole,
283 double[][][] inducedDipoleCR,
284 int[][][] neighborLists,
285 ScaleParameters scaleParameters,
286 boolean[] use,
287 int[] molecule,
288 double[] ipdamp,
289 double[] thole,
290 int[][] ip11,
291 int[][] mask12,
292 int[][] mask13,
293 int[][] mask14,
294 LambdaMode lambdaMode,
295 boolean reciprocalSpaceTerm,
296 ReciprocalSpace reciprocalSpace,
297 EwaldParameters ewaldParameters,
298 PCGSolver pcgSolver,
299 IntegerSchedule permanentSchedule,
300 RealSpaceNeighborParameters realSpaceNeighborParameters,
301 AtomicDoubleArray3D field,
302 AtomicDoubleArray3D fieldCR) {
303 this.atoms = atoms;
304 this.crystal = crystal;
305 this.coordinates = coordinates;
306 this.globalMultipole = globalMultipole;
307 this.inducedDipole = inducedDipole;
308 this.inducedDipoleCR = inducedDipoleCR;
309 this.neighborLists = neighborLists;
310 this.scaleParameters = scaleParameters;
311 this.use = use;
312 this.molecule = molecule;
313 this.ipdamp = ipdamp;
314 this.thole = thole;
315 this.ip11 = ip11;
316 this.mask12 = mask12;
317 this.mask13 = mask13;
318 this.mask14 = mask14;
319 this.lambdaMode = lambdaMode;
320 this.reciprocalSpaceTerm = reciprocalSpaceTerm;
321 this.reciprocalSpace = reciprocalSpace;
322 if (pcgSolver != null) {
323 this.preconditionerCutoff = pcgSolver.getPreconditionerCutoff();
324 this.preconditionerLists = pcgSolver.getPreconditionerLists();
325 this.preconditionerCounts = pcgSolver.getPreconditionerCounts();
326 }
327 this.aewald = ewaldParameters.aewald;
328 this.an0 = ewaldParameters.an0;
329 this.an1 = ewaldParameters.an1;
330 this.an2 = ewaldParameters.an2;
331 this.off2 = ewaldParameters.off2;
332 this.permanentSchedule = permanentSchedule;
333 this.realSpaceLists = realSpaceNeighborParameters.realSpaceLists;
334 this.realSpaceCounts = realSpaceNeighborParameters.realSpaceCounts;
335 this.realSpaceRanges = realSpaceNeighborParameters.realSpaceRanges;
336 this.field = field;
337 this.fieldCR = fieldCR;
338 }
339
340
341
342
343
344
345
346
347 @Override
348 public void removeMask(int i, boolean[] is14, double[]... masks) {
349 if (ip11[i] != null) {
350 double[] energyMask = masks[0];
351 var m12 = mask12[i];
352 for (int value : m12) {
353 energyMask[value] = 1.0;
354 }
355 var m13 = mask13[i];
356 for (int value : m13) {
357 energyMask[value] = 1.0;
358 }
359 var m14 = mask14[i];
360 for (int value : m14) {
361 energyMask[value] = 1.0;
362 for (int k : ip11[i]) {
363 if (k == value) {
364 energyMask[value] = 1.0;
365 break;
366 }
367 }
368 }
369
370 double[] inductionMask = masks[1];
371 for (int index : ip11[i]) {
372 inductionMask[index] = 1.0;
373 }
374 }
375 }
376
377 @Override
378 public void run() {
379 try {
380 execute(permanentRealSpaceFieldSection, permanentReciprocalSection);
381 } catch (RuntimeException e) {
382 String message = "Runtime exception computing the permanent multipole field.\n";
383 logger.log(Level.WARNING, message, e);
384 throw e;
385 } catch (Exception e) {
386 String message = "Fatal exception computing the permanent multipole field.\n";
387 logger.log(Level.SEVERE, message, e);
388 }
389 }
390
391
392
393
394 private class PermanentRealSpaceFieldSection extends ParallelSection {
395
396 private final PermanentRealSpaceFieldRegion permanentRealSpaceFieldRegion;
397 private final ParallelTeam parallelTeam;
398 protected long time;
399
400 PermanentRealSpaceFieldSection(ParallelTeam pt) {
401 this.parallelTeam = pt;
402 int nt = pt.getThreadCount();
403 permanentRealSpaceFieldRegion = new PermanentRealSpaceFieldRegion(nt);
404 }
405
406 public void initTimings() {
407 time = 0;
408 permanentRealSpaceFieldRegion.initTimings();
409 }
410
411 @Override
412 public void run() {
413 try {
414 time -= System.nanoTime();
415 parallelTeam.execute(permanentRealSpaceFieldRegion);
416 time += System.nanoTime();
417 } catch (RuntimeException e) {
418 String message = "Fatal exception computing the real space field.\n";
419 logger.log(Level.WARNING, message, e);
420 } catch (Exception e) {
421 String message = "Fatal exception computing the real space field.\n";
422 logger.log(Level.SEVERE, message, e);
423 }
424 }
425 }
426
427
428
429
430
431
432 private class PermanentReciprocalSection extends ParallelSection {
433
434 @Override
435 public void run() {
436 if (reciprocalSpaceTerm && aewald > 0.0) {
437 reciprocalSpace.performConvolution();
438 }
439 }
440 }
441
442 private class PermanentRealSpaceFieldRegion extends ParallelRegion {
443
444 private final InitializationLoop[] initializationLoop;
445 private final PermanentRealSpaceFieldLoop[] permanentRealSpaceFieldLoop;
446 private final SharedInteger sharedCount;
447 private final int threadCount;
448
449 PermanentRealSpaceFieldRegion(int nt) {
450 threadCount = nt;
451 initializationLoop = new InitializationLoop[threadCount];
452 permanentRealSpaceFieldLoop = new PermanentRealSpaceFieldLoop[threadCount];
453 sharedCount = new SharedInteger();
454 for (int i = 0; i < threadCount; i++) {
455 initializationLoop[i] = new InitializationLoop();
456 permanentRealSpaceFieldLoop[i] = new PermanentRealSpaceFieldLoop();
457 }
458 }
459
460 public void initTimings() {
461 for (int i = 0; i < threadCount; i++) {
462 permanentRealSpaceFieldLoop[i].time = 0;
463 initializationLoop[i].time = 0;
464 }
465 }
466
467 @Override
468 public void finish() {
469 if (realSpaceRanges == null) {
470 logger.severe(" RealSpaceRange array is null");
471 }
472
473 int nAtoms = atoms.length;
474
475
476 int id = 0;
477 int goal = sharedCount.get() / threadCount;
478 int num = 0;
479 int start = 0;
480
481 for (int i = 0; i < nAtoms; i++) {
482 List<SymOp> symOps = crystal.spaceGroup.symOps;
483 int nSymm = symOps.size();
484 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
485 num += realSpaceCounts[iSymm][i];
486 }
487 if (num >= goal) {
488
489 if (id == threadCount - 1) {
490 realSpaceRanges[id] = new Range(start, nAtoms - 1);
491 break;
492 }
493
494 realSpaceRanges[id] = new Range(start, i);
495
496
497 num = 0;
498
499
500 id++;
501
502
503 start = i + 1;
504
505
506 if (start == nAtoms) {
507 for (int j = id; j < threadCount; j++) {
508 realSpaceRanges[j] = null;
509 }
510 break;
511 }
512 } else if (i == nAtoms - 1) {
513
514
515 realSpaceRanges[id] = new Range(start, nAtoms - 1);
516 for (int j = id + 1; j < threadCount; j++) {
517 realSpaceRanges[j] = null;
518 }
519 }
520 }
521 }
522
523 @Override
524 public void run() {
525 int threadIndex = getThreadIndex();
526 try {
527 int nAtoms = atoms.length;
528 execute(0, nAtoms - 1, initializationLoop[threadIndex]);
529 execute(0, nAtoms - 1, permanentRealSpaceFieldLoop[threadIndex]);
530 } catch (RuntimeException e) {
531 String message = "Runtime exception computing the real space field.\n";
532 logger.log(Level.SEVERE, message, e);
533 } catch (Exception e) {
534 String message =
535 "Fatal exception computing the real space field in thread " + getThreadIndex() + "\n";
536 logger.log(Level.SEVERE, message, e);
537 }
538 }
539
540 @Override
541 public void start() {
542 sharedCount.set(0);
543 }
544
545 private class InitializationLoop extends IntegerForLoop {
546
547 protected long time;
548
549 @Override
550 public void start() {
551 time -= System.nanoTime();
552 }
553
554 @Override
555 public void finish() {
556 time += System.nanoTime();
557 }
558
559 @Override
560 public void run(int lb, int ub) {
561
562 List<SymOp> symOps = crystal.spaceGroup.symOps;
563 int nSymm = symOps.size();
564 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
565 double[][] ind0 = inducedDipole[0];
566 double[][] indCR0 = inducedDipoleCR[0];
567 for (int i = lb; i <= ub; i++) {
568 double[] ind = ind0[i];
569 double[] indCR = indCR0[i];
570 ind[0] = 0.0;
571 ind[1] = 0.0;
572 ind[2] = 0.0;
573 indCR[0] = 0.0;
574 indCR[1] = 0.0;
575 indCR[2] = 0.0;
576 }
577 }
578 }
579
580 @Override
581 public IntegerSchedule schedule() {
582 return IntegerSchedule.fixed();
583 }
584
585 }
586
587 private class PermanentRealSpaceFieldLoop extends IntegerForLoop {
588
589 protected long time;
590 private final double[] dx_local;
591 private final double[][] transOp;
592 private int threadID;
593 private double[] inductionMaskLocal;
594 private double[] energyMaskLocal;
595 private int count;
596
597 PermanentRealSpaceFieldLoop() {
598 super();
599 dx_local = new double[3];
600 transOp = new double[3][3];
601 }
602
603 @Override
604 public void finish() {
605 sharedCount.addAndGet(count);
606 time += System.nanoTime();
607 }
608
609 @Override
610 public void run(int lb, int ub) {
611 int[][] lists = neighborLists[0];
612 int[][] ewalds = realSpaceLists[0];
613 int[] counts = realSpaceCounts[0];
614 int[][] preLists = preconditionerLists[0];
615 int[] preCounts = preconditionerCounts[0];
616 final double[] x = coordinates[0][0];
617 final double[] y = coordinates[0][1];
618 final double[] z = coordinates[0][2];
619 final double[][] mpole = globalMultipole[0];
620
621 for (int i = lb; i <= ub; i++) {
622 if (!use[i]) {
623 continue;
624 }
625
626
627 double fix = 0.0;
628 double fiy = 0.0;
629 double fiz = 0.0;
630 double fixCR = 0.0;
631 double fiyCR = 0.0;
632 double fizCR = 0.0;
633
634 final int moleculei = molecule[i];
635 final double pdi = ipdamp[i];
636 final double pti = thole[i];
637 final double xi = x[i];
638 final double yi = y[i];
639 final double zi = z[i];
640 final double[] globalMultipolei = mpole[i];
641 final double ci = globalMultipolei[0];
642 final double dix = globalMultipolei[t100];
643 final double diy = globalMultipolei[t010];
644 final double diz = globalMultipolei[t001];
645 final double qixx = globalMultipolei[t200] * oneThird;
646 final double qiyy = globalMultipolei[t020] * oneThird;
647 final double qizz = globalMultipolei[t002] * oneThird;
648 final double qixy = globalMultipolei[t110] * oneThird;
649 final double qixz = globalMultipolei[t101] * oneThird;
650 final double qiyz = globalMultipolei[t011] * oneThird;
651
652
653 applyMask(i, null, energyMaskLocal, inductionMaskLocal);
654
655
656 final int[] list = lists[i];
657 counts[i] = 0;
658 preCounts[i] = 0;
659 int[] ewald = ewalds[i];
660 int[] preList = preLists[i];
661 for (int k : list) {
662 if (!use[k]) {
663 continue;
664 }
665 if (lambdaMode == LambdaMode.VAPOR) {
666 boolean sameMolecule = (moleculei == molecule[k]);
667 if ((intermolecularSoftcore && !sameMolecule)
668 || (intramolecularSoftcore && sameMolecule)) {
669 continue;
670 }
671 }
672 final double xk = x[k];
673 final double yk = y[k];
674 final double zk = z[k];
675 dx_local[0] = xk - xi;
676 dx_local[1] = yk - yi;
677 dx_local[2] = zk - zi;
678 final double r2 = crystal.image(dx_local);
679 if (r2 <= off2) {
680 count++;
681
682 if (ewald.length <= counts[i]) {
683 int len = ewald.length;
684 ewalds[i] = copyOf(ewald, len + 10);
685 ewald = ewalds[i];
686 }
687 ewald[counts[i]++] = k;
688 final double xr = dx_local[0];
689 final double yr = dx_local[1];
690 final double zr = dx_local[2];
691 final double pdk = ipdamp[k];
692 final double ptk = thole[k];
693 final double[] globalMultipolek = mpole[k];
694 final double ck = globalMultipolek[t000];
695 final double dkx = globalMultipolek[t100];
696 final double dky = globalMultipolek[t010];
697 final double dkz = globalMultipolek[t001];
698 final double qkxx = globalMultipolek[t200] * oneThird;
699 final double qkyy = globalMultipolek[t020] * oneThird;
700 final double qkzz = globalMultipolek[t002] * oneThird;
701 final double qkxy = globalMultipolek[t110] * oneThird;
702 final double qkxz = globalMultipolek[t101] * oneThird;
703 final double qkyz = globalMultipolek[t011] * oneThird;
704 double r = sqrt(r2);
705
706
707 if (r < preconditionerCutoff) {
708 if (preList.length <= preCounts[i]) {
709 int len = preList.length;
710 preLists[i] = copyOf(preList, len + 10);
711 preList = preLists[i];
712 }
713 preList[preCounts[i]++] = k;
714 }
715
716
717 final double ralpha = aewald * r;
718 final double exp2a = exp(-ralpha * ralpha);
719 final double rr1 = 1.0 / r;
720 final double rr2 = rr1 * rr1;
721 final double bn0 = erfc(ralpha) * rr1;
722 final double bn1 = (bn0 + an0 * exp2a) * rr2;
723 final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
724 final double bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
725
726
727 double scale3 = 1.0;
728 double scale5 = 1.0;
729 double scale7 = 1.0;
730 double damp = pdi * pdk;
731 final double pgamma = min(pti, ptk);
732 final double rdamp = r * damp;
733 damp = -pgamma * rdamp * rdamp * rdamp;
734 if (damp > -50.0) {
735 double expdamp = exp(damp);
736 scale3 = 1.0 - expdamp;
737 scale5 = 1.0 - expdamp * (1.0 - damp);
738 scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
739 }
740 final double scale = inductionMaskLocal[k];
741 final double scalep = energyMaskLocal[k];
742
743 final double dsc3 = scale3 * scale;
744 final double dsc5 = scale5 * scale;
745 final double dsc7 = scale7 * scale;
746
747 final double psc3 = scale3 * scalep;
748 final double psc5 = scale5 * scalep;
749 final double psc7 = scale7 * scalep;
750 final double rr3 = rr1 * rr2;
751 final double rr5 = 3.0 * rr3 * rr2;
752 final double rr7 = 5.0 * rr5 * rr2;
753
754 final double drr3 = (1.0 - dsc3) * rr3;
755 final double drr5 = (1.0 - dsc5) * rr5;
756 final double drr7 = (1.0 - dsc7) * rr7;
757
758 final double prr3 = (1.0 - psc3) * rr3;
759 final double prr5 = (1.0 - psc5) * rr5;
760 final double prr7 = (1.0 - psc7) * rr7;
761 final double dir = dix * xr + diy * yr + diz * zr;
762 final double qix = 2.0 * (qixx * xr + qixy * yr + qixz * zr);
763 final double qiy = 2.0 * (qixy * xr + qiyy * yr + qiyz * zr);
764 final double qiz = 2.0 * (qixz * xr + qiyz * yr + qizz * zr);
765 final double qir = (qix * xr + qiy * yr + qiz * zr) * 0.5;
766
767 final double bn123i = bn1 * ci + bn2 * dir + bn3 * qir;
768 final double fkmx = xr * bn123i - bn1 * dix - bn2 * qix;
769 final double fkmy = yr * bn123i - bn1 * diy - bn2 * qiy;
770 final double fkmz = zr * bn123i - bn1 * diz - bn2 * qiz;
771
772 final double ddr357i = drr3 * ci + drr5 * dir + drr7 * qir;
773 final double fkdx = xr * ddr357i - drr3 * dix - drr5 * qix;
774 final double fkdy = yr * ddr357i - drr3 * diy - drr5 * qiy;
775 final double fkdz = zr * ddr357i - drr3 * diz - drr5 * qiz;
776 field.add(threadID, k, fkmx - fkdx, fkmy - fkdy, fkmz - fkdz);
777
778 final double prr357i = prr3 * ci + prr5 * dir + prr7 * qir;
779 final double fkpx = xr * prr357i - prr3 * dix - prr5 * qix;
780 final double fkpy = yr * prr357i - prr3 * diy - prr5 * qiy;
781 final double fkpz = zr * prr357i - prr3 * diz - prr5 * qiz;
782 fieldCR.add(threadID, k, fkmx - fkpx, fkmy - fkpy, fkmz - fkpz);
783 final double dkr = dkx * xr + dky * yr + dkz * zr;
784 final double qkx = 2.0 * (qkxx * xr + qkxy * yr + qkxz * zr);
785 final double qky = 2.0 * (qkxy * xr + qkyy * yr + qkyz * zr);
786 final double qkz = 2.0 * (qkxz * xr + qkyz * yr + qkzz * zr);
787 final double qkr = (qkx * xr + qky * yr + qkz * zr) * 0.5;
788 final double bn123k = bn1 * ck - bn2 * dkr + bn3 * qkr;
789
790 final double fimx = -xr * bn123k - bn1 * dkx + bn2 * qkx;
791 final double fimy = -yr * bn123k - bn1 * dky + bn2 * qky;
792 final double fimz = -zr * bn123k - bn1 * dkz + bn2 * qkz;
793
794 final double drr357k = drr3 * ck - drr5 * dkr + drr7 * qkr;
795 final double fidx = -xr * drr357k - drr3 * dkx + drr5 * qkx;
796 final double fidy = -yr * drr357k - drr3 * dky + drr5 * qky;
797 final double fidz = -zr * drr357k - drr3 * dkz + drr5 * qkz;
798 fix += fimx - fidx;
799 fiy += fimy - fidy;
800 fiz += fimz - fidz;
801
802 final double prr357k = prr3 * ck - prr5 * dkr + prr7 * qkr;
803 final double fipx = -xr * prr357k - prr3 * dkx + prr5 * qkx;
804 final double fipy = -yr * prr357k - prr3 * dky + prr5 * qky;
805 final double fipz = -zr * prr357k - prr3 * dkz + prr5 * qkz;
806 fixCR += fimx - fipx;
807 fiyCR += fimy - fipy;
808 fizCR += fimz - fipz;
809 }
810 }
811
812 field.add(threadID, i, fix, fiy, fiz);
813 fieldCR.add(threadID, i, fixCR, fiyCR, fizCR);
814
815 removeMask(i, null, energyMaskLocal, inductionMaskLocal);
816 }
817
818
819 List<SymOp> symOps = crystal.spaceGroup.symOps;
820 int nSymm = symOps.size();
821 for (int iSymm = 1; iSymm < nSymm; iSymm++) {
822 SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
823 crystal.getTransformationOperator(symOp, transOp);
824 lists = neighborLists[iSymm];
825 ewalds = realSpaceLists[iSymm];
826 counts = realSpaceCounts[iSymm];
827 preLists = preconditionerLists[iSymm];
828 preCounts = preconditionerCounts[iSymm];
829 double[] xs = coordinates[iSymm][0];
830 double[] ys = coordinates[iSymm][1];
831 double[] zs = coordinates[iSymm][2];
832 double[][] mpoles = globalMultipole[iSymm];
833
834
835 for (int i = lb; i <= ub; i++) {
836 if (!use[i]) {
837 continue;
838 }
839
840 double fix = 0.0;
841 double fiy = 0.0;
842 double fiz = 0.0;
843 final double pdi = ipdamp[i];
844 final double pti = thole[i];
845 final double[] multipolei = mpole[i];
846 final double ci = multipolei[t000];
847 final double dix = multipolei[t100];
848 final double diy = multipolei[t010];
849 final double diz = multipolei[t001];
850 final double qixx = multipolei[t200] * oneThird;
851 final double qiyy = multipolei[t020] * oneThird;
852 final double qizz = multipolei[t002] * oneThird;
853 final double qixy = multipolei[t110] * oneThird;
854 final double qixz = multipolei[t101] * oneThird;
855 final double qiyz = multipolei[t011] * oneThird;
856 final double xi = x[i];
857 final double yi = y[i];
858 final double zi = z[i];
859
860
861 final int[] list = lists[i];
862 counts[i] = 0;
863 preCounts[i] = 0;
864 int[] ewald = ewalds[i];
865 int[] preList = preLists[i];
866 for (int k : list) {
867 if (!use[k]) {
868 continue;
869 }
870 final double xk = xs[k];
871 final double yk = ys[k];
872 final double zk = zs[k];
873 dx_local[0] = xk - xi;
874 dx_local[1] = yk - yi;
875 dx_local[2] = zk - zi;
876 final double r2 = crystal.image(dx_local);
877 if (r2 <= off2) {
878 count++;
879
880 if (ewald.length <= counts[i]) {
881 int len = ewald.length;
882 ewalds[i] = copyOf(ewald, len + 10);
883 ewald = ewalds[i];
884 }
885 ewald[counts[i]++] = k;
886 double selfScale = 1.0;
887 if (i == k) {
888 selfScale = 0.5;
889 }
890 final double xr = dx_local[0];
891 final double yr = dx_local[1];
892 final double zr = dx_local[2];
893 final double pdk = ipdamp[k];
894 final double ptk = thole[k];
895 final double[] multipolek = mpoles[k];
896 final double ck = multipolek[t000];
897 final double dkx = multipolek[t100];
898 final double dky = multipolek[t010];
899 final double dkz = multipolek[t001];
900 final double qkxx = multipolek[t200] * oneThird;
901 final double qkyy = multipolek[t020] * oneThird;
902 final double qkzz = multipolek[t002] * oneThird;
903 final double qkxy = multipolek[t110] * oneThird;
904 final double qkxz = multipolek[t101] * oneThird;
905 final double qkyz = multipolek[t011] * oneThird;
906 final double r = sqrt(r2);
907 if (r < preconditionerCutoff) {
908 if (preList.length <= preCounts[i]) {
909 int len = preList.length;
910 preLists[i] = copyOf(preList, len + 10);
911 preList = preLists[i];
912 }
913 preList[preCounts[i]++] = k;
914 }
915
916
917 final double ralpha = aewald * r;
918 final double exp2a = exp(-ralpha * ralpha);
919 final double rr1 = 1.0 / r;
920 final double rr2 = rr1 * rr1;
921 final double bn0 = erfc(ralpha) * rr1;
922 final double bn1 = (bn0 + an0 * exp2a) * rr2;
923 final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
924 final double bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
925
926
927 double scale3 = 1.0;
928 double scale5 = 1.0;
929 double scale7 = 1.0;
930 double damp = pdi * pdk;
931 final double pgamma = min(pti, ptk);
932 final double rdamp = r * damp;
933 damp = -pgamma * rdamp * rdamp * rdamp;
934 if (damp > -50.0) {
935 double expdamp = exp(damp);
936 scale3 = 1.0 - expdamp;
937 scale5 = 1.0 - expdamp * (1.0 - damp);
938 scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
939 }
940
941 final double dsc3 = scale3;
942 final double dsc5 = scale5;
943 final double dsc7 = scale7;
944 final double rr3 = rr1 * rr2;
945 final double rr5 = 3.0 * rr3 * rr2;
946 final double rr7 = 5.0 * rr5 * rr2;
947 final double drr3 = (1.0 - dsc3) * rr3;
948 final double drr5 = (1.0 - dsc5) * rr5;
949 final double drr7 = (1.0 - dsc7) * rr7;
950
951 final double dkr = dkx * xr + dky * yr + dkz * zr;
952 final double qkx = 2.0 * (qkxx * xr + qkxy * yr + qkxz * zr);
953 final double qky = 2.0 * (qkxy * xr + qkyy * yr + qkyz * zr);
954 final double qkz = 2.0 * (qkxz * xr + qkyz * yr + qkzz * zr);
955 final double qkr = (qkx * xr + qky * yr + qkz * zr) * 0.5;
956 final double bn123k = bn1 * ck - bn2 * dkr + bn3 * qkr;
957 final double drr357k = drr3 * ck - drr5 * dkr + drr7 * qkr;
958 final double fimx = -xr * bn123k - bn1 * dkx + bn2 * qkx;
959 final double fimy = -yr * bn123k - bn1 * dky + bn2 * qky;
960 final double fimz = -zr * bn123k - bn1 * dkz + bn2 * qkz;
961 final double fidx = -xr * drr357k - drr3 * dkx + drr5 * qkx;
962 final double fidy = -yr * drr357k - drr3 * dky + drr5 * qky;
963 final double fidz = -zr * drr357k - drr3 * dkz + drr5 * qkz;
964
965 final double dir = dix * xr + diy * yr + diz * zr;
966 final double qix = 2.0 * (qixx * xr + qixy * yr + qixz * zr);
967 final double qiy = 2.0 * (qixy * xr + qiyy * yr + qiyz * zr);
968 final double qiz = 2.0 * (qixz * xr + qiyz * yr + qizz * zr);
969 final double qir = (qix * xr + qiy * yr + qiz * zr) * 0.5;
970 final double bn123i = bn1 * ci + bn2 * dir + bn3 * qir;
971 final double ddr357i = drr3 * ci + drr5 * dir + drr7 * qir;
972 final double fkmx = xr * bn123i - bn1 * dix - bn2 * qix;
973 final double fkmy = yr * bn123i - bn1 * diy - bn2 * qiy;
974 final double fkmz = zr * bn123i - bn1 * diz - bn2 * qiz;
975 final double fkdx = xr * ddr357i - drr3 * dix - drr5 * qix;
976 final double fkdy = yr * ddr357i - drr3 * diy - drr5 * qiy;
977 final double fkdz = zr * ddr357i - drr3 * diz - drr5 * qiz;
978 fix += selfScale * (fimx - fidx);
979 fiy += selfScale * (fimy - fidy);
980 fiz += selfScale * (fimz - fidz);
981 final double xc = selfScale * (fkmx - fkdx);
982 final double yc = selfScale * (fkmy - fkdy);
983 final double zc = selfScale * (fkmz - fkdz);
984 final double fkx = xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0];
985 final double fky = xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1];
986 final double fkz = xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2];
987 field.add(threadID, k, fkx, fky, fkz);
988 fieldCR.add(threadID, k, fkx, fky, fkz);
989 }
990 }
991 field.add(threadID, i, fix, fiy, fiz);
992 fieldCR.add(threadID, i, fix, fiy, fiz);
993 }
994 }
995 }
996
997 @Override
998 public IntegerSchedule schedule() {
999 return permanentSchedule;
1000 }
1001
1002 @Override
1003 public void start() {
1004 time = -System.nanoTime();
1005 threadID = getThreadIndex();
1006 count = 0;
1007 int nAtoms = atoms.length;
1008 if (inductionMaskLocal == null || inductionMaskLocal.length < nAtoms) {
1009 inductionMaskLocal = new double[nAtoms];
1010 energyMaskLocal = new double[nAtoms];
1011 fill(inductionMaskLocal, 1.0);
1012 fill(energyMaskLocal, 1.0);
1013 }
1014 }
1015 }
1016 }
1017 }