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;
39
40 import edu.rit.pj.BarrierAction;
41 import edu.rit.pj.IntegerForLoop;
42 import edu.rit.pj.IntegerSchedule;
43 import edu.rit.pj.ParallelRegion;
44 import edu.rit.pj.ParallelTeam;
45 import edu.rit.pj.reduction.SharedBoolean;
46 import edu.rit.pj.reduction.SharedInteger;
47 import edu.rit.util.Range;
48 import ffx.crystal.Crystal;
49 import ffx.potential.MolecularAssembly;
50 import ffx.potential.bonded.Atom;
51 import ffx.potential.utils.PotentialsUtils;
52
53 import java.io.File;
54 import java.util.*;
55 import java.util.logging.Level;
56 import java.util.logging.Logger;
57
58 import static java.lang.Math.max;
59 import static java.lang.String.format;
60 import static java.lang.System.arraycopy;
61 import static java.util.Arrays.*;
62 import static org.apache.commons.math3.util.FastMath.floor;
63 import static org.apache.commons.math3.util.FastMath.min;
64 import static org.apache.commons.math3.util.FastMath.sqrt;
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93 public class NeighborList extends ParallelRegion {
94
95
96 private static final Logger logger = Logger.getLogger(NeighborList.class.getName());
97
98 private final DomainDecomposition domainDecomposition;
99
100 private final MaskingInterface maskingRules;
101
102 private final double cutoff;
103
104
105
106
107 private final double buffer;
108
109 private final double motion2;
110
111 private final double cutoffPlusBuffer;
112
113 private final double cutoffPlusBuffer2;
114
115 private final ParallelTeam parallelTeam;
116
117 private final int threadCount;
118
119
120
121 private boolean forceRebuild = false;
122
123
124
125 private boolean print = false;
126
127
128
129
130 private final SharedBoolean sharedMotion;
131
132
133
134 private final MotionLoop[] motionLoops;
135
136
137
138 private final ListInitBarrierAction listInitBarrierAction;
139
140
141
142 private final AssignAtomsToCellsLoop[] assignAtomsToCellsLoops;
143
144 private final NeighborListLoop[] verletListLoop;
145 private final MxNListLoop[] mxNListLoop;
146
147
148
149 private long motionTime;
150
151
152
153 private long initTime;
154
155
156
157 private long assignAtomsToCellsTime;
158
159
160
161 private long verletListTime;
162
163 private long groupingTime;
164
165 private final SharedInteger sharedCount;
166
167 private final Range[] ranges;
168
169
170
171
172 private Crystal crystal;
173
174 private int nSymm;
175
176 private int nAtoms;
177
178 private double[][] coordinates;
179
180 private double[] previous;
181
182 private int[][][] lists;
183
184 private int[][][] groups;
185 private int nGroups;
186
187 private int[][][] groupLists;
188
189 private int M = -1;
190
191 private int N = -1;
192
193 private int[] listCount;
194
195 private PairwiseSchedule pairwiseSchedule;
196
197 private int asymmetricUnitCount;
198
199 private int symmetryMateCount;
200
201 private int atomsWithIteractions;
202
203 private boolean[] use;
204
205 private Atom[] atoms;
206
207
208
209 private long time;
210
211 private boolean intermolecular = true;
212
213
214
215
216
217 private final boolean includeInactivePairs = true;
218
219 private boolean disableUpdates = false;
220
221
222
223
224
225
226
227
228
229
230
231
232 public NeighborList(MaskingInterface maskingRules, Crystal crystal, Atom[] atoms, double cutoff,
233 double buffer, ParallelTeam parallelTeam) {
234 this.maskingRules = maskingRules;
235 this.crystal = crystal;
236 this.cutoff = cutoff;
237 this.buffer = buffer;
238 this.parallelTeam = new ParallelTeam(parallelTeam.getThreadCount());
239 this.atoms = atoms;
240 nAtoms = atoms.length;
241
242
243 cutoffPlusBuffer = cutoff + buffer;
244 cutoffPlusBuffer2 = cutoffPlusBuffer * cutoffPlusBuffer;
245 motion2 = (buffer / 2.0) * (buffer / 2.0);
246
247
248 threadCount = parallelTeam.getThreadCount();
249 sharedCount = new SharedInteger();
250 ranges = new Range[threadCount];
251
252 sharedMotion = new SharedBoolean();
253 motionLoops = new MotionLoop[threadCount];
254 listInitBarrierAction = new ListInitBarrierAction();
255 verletListLoop = new NeighborListLoop[threadCount];
256 assignAtomsToCellsLoops = new AssignAtomsToCellsLoop[threadCount];
257 mxNListLoop = new MxNListLoop[threadCount];
258 for (int i = 0; i < threadCount; i++) {
259 motionLoops[i] = new MotionLoop();
260 verletListLoop[i] = new NeighborListLoop();
261 assignAtomsToCellsLoops[i] = new AssignAtomsToCellsLoop();
262 mxNListLoop[i] = new MxNListLoop();
263 }
264
265
266 boolean print = logger.isLoggable(Level.FINE);
267 domainDecomposition = new DomainDecomposition(nAtoms, crystal, cutoffPlusBuffer);
268 initNeighborList(print);
269 }
270
271
272
273
274
275
276
277
278
279
280
281
282 public void buildList(final double[][] coordinates, final int[][][] lists, boolean[] use,
283 boolean forceRebuild, boolean print) {
284 if (disableUpdates) {
285 return;
286 }
287 this.coordinates = coordinates;
288 this.lists = lists;
289 this.use = use;
290 this.forceRebuild = forceRebuild;
291 this.print = print;
292
293 try {
294 parallelTeam.execute(this);
295 } catch (Exception e) {
296 String message = "Fatal exception building neighbor list.\n";
297 logger.log(Level.SEVERE, message, e);
298 }
299 }
300
301 public void buildMxNList(int M, int N, final double[][] coordinates, final int[][][] lists, boolean[] use, boolean forceRebuild, boolean print){
302 if (disableUpdates) {
303 return;
304 }
305 this.coordinates = coordinates;
306 this.lists = lists;
307 this.M = M;
308 this.N = N;
309 this.use = use;
310 this.forceRebuild = forceRebuild;
311 this.print = print;
312
313 try {
314 parallelTeam.execute(this);
315 } catch (Exception e) {
316 String message = "Fatal exception building neighbor list.\n";
317 logger.log(Level.SEVERE, message, e);
318 }
319 }
320
321
322
323
324
325
326 public void destroy() throws Exception {
327 parallelTeam.shutdown();
328 }
329
330
331
332
333
334
335
336
337 @Override
338 public void finish() {
339 if (disableUpdates) {
340 return;
341 }
342 if (!forceRebuild && !sharedMotion.get()) {
343 return;
344 }
345
346
347 atomsWithIteractions = 0;
348 asymmetricUnitCount = 0;
349 symmetryMateCount = 0;
350 int[][] list = lists[0];
351 for (int i = 0; i < nAtoms; i++) {
352 asymmetricUnitCount += list[i].length;
353 if (listCount[i] > 0) {
354 atomsWithIteractions++;
355 }
356 }
357 for (int iSymm = 1; iSymm < nSymm; iSymm++) {
358 list = lists[iSymm];
359 for (int i = 0; i < nAtoms; i++) {
360 symmetryMateCount += list[i].length;
361 }
362 }
363
364
365 if (print) {
366 print();
367 }
368
369
370 long scheduleTime = -System.nanoTime();
371 pairwiseSchedule.updateRanges(sharedCount.get(), atomsWithIteractions, listCount);
372 scheduleTime += System.nanoTime();
373
374 if (logger.isLoggable(Level.FINE)) {
375 time = System.nanoTime() - time;
376 StringBuilder sb = new StringBuilder();
377 sb.append(format(" Motion Check: %6.4f sec\n", motionTime * 1e-9));
378 sb.append(format(" List Initialization: %6.4f sec\n", initTime * 1e-9));
379 sb.append(format(" Assign Atoms to Cells: %6.4f sec\n", assignAtomsToCellsTime * 1e-9));
380 sb.append(format(" Create Vertlet Lists: %6.4f sec\n", verletListTime * 1e-9));
381 sb.append(format(" Parallel Schedule: %6.4f sec\n", scheduleTime * 1e-9));
382 if ( M >= 1 || N >= 1) {
383 sb.append(format(" M x N List Total: %6.4f sec\n", groupingTime * 1e-9));
384 }
385 sb.append(format(" Neighbor List Total: %6.4f sec\n", time * 1e-9));
386 logger.fine(sb.toString());
387 }
388 }
389
390
391
392
393
394
395 public double getCutoff() {
396 return cutoff;
397 }
398
399
400
401
402
403
404
405 public void setDisableUpdates(boolean disableUpdate) {
406 this.disableUpdates = disableUpdate;
407 }
408
409
410
411
412
413
414 public int[][][] getNeighborList() {
415 return lists;
416 }
417
418
419
420
421
422
423 public PairwiseSchedule getPairwiseSchedule() {
424 return pairwiseSchedule;
425 }
426
427
428
429
430
431
432
433
434 @Override
435 public void run() {
436 if (disableUpdates) {
437 return;
438 }
439
440 int threadIndex = getThreadIndex();
441 try {
442 if (!forceRebuild) {
443
444 if (threadIndex == 0) {
445 motionTime = -System.nanoTime();
446 }
447 execute(0, nAtoms - 1, motionLoops[threadIndex]);
448 if (threadIndex == 0) {
449 motionTime += System.nanoTime();
450 }
451 }
452 if (forceRebuild || sharedMotion.get()) {
453
454 barrier(listInitBarrierAction);
455
456 if (threadIndex == 0) {
457 assignAtomsToCellsTime = -System.nanoTime();
458 }
459 execute(0, nAtoms - 1, assignAtomsToCellsLoops[threadIndex]);
460 if (threadIndex == 0) {
461 assignAtomsToCellsTime += System.nanoTime();
462 }
463
464 if (threadIndex == 0) {
465 verletListTime = -System.nanoTime();
466 }
467 execute(0, nAtoms - 1, verletListLoop[threadIndex]);
468 if (threadIndex == 0) {
469 verletListTime += System.nanoTime();
470 }
471
472 if ( M >= 1 || N >= 1) {
473 if(threadIndex == 0){
474 groupingTime = -System.nanoTime();
475 }
476
477 if (threadIndex == 0) { groups(); }
478 barrier();
479 execute(0, nGroups - 1, mxNListLoop[threadIndex]);
480 if(threadIndex == 0){
481 groupingTime += System.nanoTime();
482 }
483 }
484 }
485 } catch (Exception e) {
486 String message =
487 "Fatal exception building neighbor list in thread: " + getThreadIndex() + "\n";
488 logger.log(Level.SEVERE, message, e);
489 }
490 }
491
492
493
494
495
496
497
498
499 private void groups() {
500
501 int nA = domainDecomposition.nA;
502 int nB = domainDecomposition.nB;
503 int nC = domainDecomposition.nC;
504
505 groups = new int[nSymm][nAtoms][M];
506 int[] symmGroups = new int[nSymm];
507 for(int i = 0; i < nA; i++){
508 for(int j = 0; j < nB; j++){
509 int[] group = new int[M];
510
511 List<PriorityQueue<AtomIndex>> atomQueue = new ArrayList<>(nSymm);
512 for(int k = 0; k < nSymm; k++){
513 atomQueue.add(new PriorityQueue<>());
514 }
515 for(int k = 0; k < nC; k++){
516
517 Cell cell = domainDecomposition.getCell(i, j, k);
518 assert(cell != null);
519 for(int l = 0; l < cell.list.size(); l++){
520 int iSymm = cell.list.get(l).iSymm;
521
522 PriorityQueue<AtomIndex> queue = atomQueue.get(iSymm);
523 queue.add(cell.list.get(l));
524
525 while(queue.size() >= M){
526 for(int m = 0; m < M; m++){
527 AtomIndex atom = queue.poll();
528 assert(atom != null);
529 group[m] = atom.i;
530 }
531 groups[iSymm][symmGroups[iSymm]++] = group;
532 group = new int[M];
533 }
534 }
535 }
536 for(int k = 0; k < nSymm; k++) {
537 PriorityQueue<AtomIndex> queue = atomQueue.get(k);
538 if (queue.isEmpty()) {
539 continue;
540 }
541
542 for (int l = 0; l < M; l++) {
543 if (!queue.isEmpty()) {
544 AtomIndex atom = queue.poll();
545 assert (atom != null);
546 group[l] = atom.i;
547 } else {
548 group[l] = -1;
549 }
550 }
551 groups[k][symmGroups[k]++] = group;
552 }
553 }
554 }
555 nGroups = symmGroups[0];
556 groupLists = new int[nSymm][nGroups][];
557 }
558
559
560
561
562
563
564 public void setAtoms(Atom[] atoms) {
565 this.atoms = atoms;
566 this.nAtoms = atoms.length;
567 initNeighborList(false);
568 }
569
570
571
572
573
574
575
576 public void setCrystal(Crystal crystal) {
577 this.crystal = crystal;
578 initNeighborList(false);
579 }
580
581
582
583
584
585
586 public void setIntermolecular(boolean intermolecular) {
587 this.intermolecular = intermolecular;
588 }
589
590
591
592
593
594
595
596
597 @Override
598 public void start() {
599 if (disableUpdates) {
600 return;
601 }
602 time = System.nanoTime();
603 motionTime = 0;
604 initTime = 0;
605 assignAtomsToCellsTime = 0;
606 verletListTime = 0;
607 sharedMotion.set(false);
608 }
609
610 private void initNeighborList(boolean print) {
611
612
613 int newNSymm = crystal.spaceGroup.symOps.size();
614 if (nSymm != newNSymm) {
615 nSymm = newNSymm;
616 }
617
618
619 if (previous == null || previous.length < 3 * nAtoms) {
620 previous = new double[3 * nAtoms];
621 listCount = new int[nAtoms];
622 pairwiseSchedule = new PairwiseSchedule(threadCount, nAtoms, ranges);
623 } else {
624 pairwiseSchedule.setAtoms(nAtoms);
625 }
626
627
628 final double sphere = min(min(crystal.interfacialRadiusA, crystal.interfacialRadiusB),
629 crystal.interfacialRadiusC);
630
631
632
633
634
635 if (!crystal.aperiodic()) {
636 assert (sphere >= cutoffPlusBuffer);
637 }
638
639 domainDecomposition.initDomainDecomposition(nAtoms, crystal);
640 if (print) {
641 domainDecomposition.log();
642 }
643 }
644
645 private void print() {
646 StringBuilder sb = new StringBuilder(
647 format(" Buffer: %5.2f (A)\n", buffer));
648 sb.append(format(" Cut-off: %5.2f (A)\n", cutoff));
649 sb.append(format(" Total: %5.2f (A)\n", cutoffPlusBuffer));
650 sb.append(format(" Neighbors in the asymmetric unit:%11d\n", asymmetricUnitCount));
651 if (nSymm > 1) {
652 int num = (asymmetricUnitCount + symmetryMateCount) * nSymm;
653 sb.append(format(" Neighbors in symmetry mates: %12d\n", symmetryMateCount));
654 sb.append(format(" Neighbors in the unit cell: %12d\n", num));
655 }
656 logger.info(sb.toString());
657 }
658
659 private static final int XX = 0;
660 private static final int YY = 1;
661 private static final int ZZ = 2;
662
663
664
665
666
667
668 private class MotionLoop extends IntegerForLoop {
669
670 @Override
671 public void run(int lb, int ub) {
672 double[] current = coordinates[0];
673 for (int i = lb; i <= ub; i++) {
674 int i3 = i * 3;
675 int iX = i3 + XX;
676 int iY = i3 + YY;
677 int iZ = i3 + ZZ;
678 double dx = previous[iX] - current[iX];
679 double dy = previous[iY] - current[iY];
680 double dz = previous[iZ] - current[iZ];
681 double dr2 = crystal.image(dx, dy, dz);
682 if (dr2 > motion2) {
683 if (logger.isLoggable(Level.FINE)) {
684 logger.fine(format(" Motion detected for atom %d (%8.6f A).", i, sqrt(dr2)));
685 }
686 sharedMotion.set(true);
687 }
688 }
689 }
690 }
691
692
693
694
695 private class ListInitBarrierAction extends BarrierAction {
696
697 @Override
698 public void run() {
699 initTime = -System.nanoTime();
700
701 sharedCount.set(0);
702
703 domainDecomposition.clear();
704
705
706 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
707 if (lists[iSymm] == null || lists[iSymm].length < nAtoms) {
708 lists[iSymm] = new int[nAtoms][];
709 }
710 }
711 initTime += System.nanoTime();
712 }
713 }
714
715
716
717
718
719
720 private class AssignAtomsToCellsLoop extends IntegerForLoop {
721
722
723
724
725 private final double[] auCart = new double[3];
726
727
728
729 private final double[] cart = new double[3];
730
731
732
733 private final double[] frac = new double[3];
734
735 @Override
736 public void run(int lb, int ub) throws Exception {
737
738
739 double[] current = coordinates[0];
740 for (int i = lb; i <= ub; i++) {
741 int i3 = i * 3;
742 int iX = i3 + XX;
743 int iY = i3 + YY;
744 int iZ = i3 + ZZ;
745 previous[iX] = current[iX];
746 previous[iY] = current[iY];
747 previous[iZ] = current[iZ];
748 }
749
750 final double[] auXYZ = coordinates[0];
751 final double spCutoff2 = crystal.getSpecialPositionCutoff2();
752 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
753 final double[] xyz = coordinates[iSymm];
754
755 for (int i = lb; i <= ub; i++) {
756 int i3 = i * 3;
757 cart[0] = xyz[i3 + XX];
758 cart[1] = xyz[i3 + YY];
759 cart[2] = xyz[i3 + ZZ];
760 if (iSymm != 0) {
761 auCart[0] = auXYZ[i3 + XX];
762 auCart[1] = auXYZ[i3 + YY];
763 auCart[2] = auXYZ[i3 + ZZ];
764 double dx = auCart[0] - cart[0];
765 double dy = auCart[1] - cart[1];
766 double dz = auCart[2] - cart[2];
767 double dr2 = crystal.image(dx, dy, dz);
768 if (dr2 < spCutoff2) {
769 int molecule = atoms[i].getMoleculeNumber();
770 if (atoms[i].isActive()) {
771 logger.info(format(" Active Atom %s at Special Position in Molecule %d with SymOp %d (%8.6f A).",
772 atoms[i], molecule, iSymm, sqrt(dr2)));
773 } else if (logger.isLoggable(Level.FINE)) {
774 logger.fine(format(" Atom %s at Special Position in Molecule %d with SymOp %d (%8.6f A).",
775 atoms[i], molecule, iSymm, sqrt(dr2)));
776 }
777
778 domainDecomposition.addSpecialPositionExclusion(molecule, iSymm);
779 }
780 }
781
782 crystal.toFractionalCoordinates(cart, frac);
783 domainDecomposition.addAtomToCell(i, iSymm, frac);
784 }
785 }
786 }
787 }
788
789
790
791
792
793
794
795
796 private class NeighborListLoop extends IntegerForLoop {
797
798 private final IntegerSchedule schedule;
799 private int n;
800 private int iSymm;
801 private int atomIndex;
802 private boolean iactive = true;
803 private int count;
804 private double[] xyz;
805 private int[] pairs;
806 private Cell cellForCurrentAtom;
807 private int[] pairCellAtoms;
808 private double[] mask;
809 private boolean[] vdw14;
810
811 NeighborListLoop() {
812 int len = 1000;
813 pairs = new int[len];
814 schedule = IntegerSchedule.dynamic(10);
815 pairCellAtoms = new int[len];
816 }
817
818 @Override
819 public void finish() {
820 sharedCount.addAndGet(count);
821 }
822
823 @Override
824 public void run(final int lb, final int ub) {
825 int nEdge = domainDecomposition.nEdge;
826 int nA = domainDecomposition.nA;
827 int nB = domainDecomposition.nB;
828 int nC = domainDecomposition.nC;
829 for (iSymm = 0; iSymm < nSymm; iSymm++) {
830 int[][] list = lists[iSymm];
831
832 for (atomIndex = lb; atomIndex <= ub; atomIndex++) {
833 n = 0;
834
835 if (iSymm == 0) {
836 listCount[atomIndex] = 0;
837 int molecule = atoms[atomIndex].getMoleculeNumber();
838 List<Integer> symOpList = domainDecomposition.getSpecialPositionSymOps(molecule);
839 atoms[atomIndex].setSpecialPositionSymOps(symOpList);
840 }
841
842 if (use == null || use[atomIndex]) {
843
844 iactive = atoms[atomIndex].isActive();
845
846 cellForCurrentAtom = domainDecomposition.getCellForAtom(atomIndex);
847 int a = cellForCurrentAtom.a;
848 int b = cellForCurrentAtom.b;
849 int c = cellForCurrentAtom.c;
850 int a1 = a + 1;
851 int b1 = b + 1;
852 int c1 = c + 1;
853
854
855
856
857
858
859 int aStart = nA == 1 ? a : a - nEdge;
860 int aStop = nA == 1 ? a : a + nEdge;
861 int bStart = nB == 1 ? b : b - nEdge;
862 int bStop = nB == 1 ? b : b + nEdge;
863 int cStart = nC == 1 ? c : c - nEdge;
864 int cStop = nC == 1 ? c : c + nEdge;
865
866 if (iSymm == 0) {
867
868 atomCellPairs(cellForCurrentAtom);
869
870
871 for (int bi = b1; bi <= bStop; bi++) {
872 atomCellPairs(domainDecomposition.image(a, bi, c));
873 }
874
875 for (int bi = bStart; bi <= bStop; bi++) {
876 for (int ci = c1; ci <= cStop; ci++) {
877 atomCellPairs(domainDecomposition.image(a, bi, ci));
878 }
879 }
880
881 for (int bi = bStart; bi <= bStop; bi++) {
882 for (int ci = cStart; ci <= cStop; ci++) {
883 for (int ai = a1; ai <= aStop; ai++) {
884 atomCellPairs(domainDecomposition.image(ai, bi, ci));
885 }
886 }
887 }
888 } else {
889
890 for (int ai = aStart; ai <= aStop; ai++) {
891 for (int bi = bStart; bi <= bStop; bi++) {
892 for (int ci = cStart; ci <= cStop; ci++) {
893 atomCellPairs(domainDecomposition.image(ai, bi, ci));
894 }
895 }
896 }
897 }
898 }
899
900 list[atomIndex] = new int[n];
901 listCount[atomIndex] += n;
902 count += n;
903 arraycopy(pairs, 0, list[atomIndex], 0, n);
904 }
905 }
906 }
907
908 @Override
909 public IntegerSchedule schedule() {
910 return schedule;
911 }
912
913 @Override
914 public void start() {
915 xyz = coordinates[0];
916 count = 0;
917 if (mask == null || mask.length < nAtoms) {
918 mask = new double[nAtoms];
919 fill(mask, 1.0);
920 vdw14 = new boolean[nAtoms];
921 fill(vdw14, false);
922 }
923 }
924
925 private void atomCellPairs(final Cell cell) {
926 if (pairCellAtoms.length < cell.getCount()) {
927 pairCellAtoms = new int[cell.getCount()];
928 }
929
930
931 int num = cell.getSymOpAtoms(iSymm, pairCellAtoms);
932 if (num == 0) {
933 return;
934 }
935
936
937 if (iSymm == 0 && maskingRules != null) {
938
939 maskingRules.applyMask(atomIndex, vdw14, mask);
940 }
941
942 final int moleculeIndex = atoms[atomIndex].getMoleculeNumber();
943
944 boolean logClosePairs = logger.isLoggable(Level.FINE);
945
946 final int i3 = atomIndex * 3;
947 final double xi = xyz[i3 + XX];
948 final double yi = xyz[i3 + YY];
949 final double zi = xyz[i3 + ZZ];
950 final double[] pair = coordinates[iSymm];
951
952
953 for (int j = 0; j < num; j++) {
954 final int aj = pairCellAtoms[j];
955
956
957 if (iSymm == 0 && cell == cellForCurrentAtom) {
958
959 if (aj <= atomIndex) {
960 continue;
961 }
962 }
963
964 if (use != null && !use[aj]) {
965 continue;
966 }
967
968 if (!includeInactivePairs && !iactive) {
969 if (!atoms[aj].isActive()) {
970 continue;
971 }
972 }
973
974 int moleculeIndexJ = atoms[aj].getMoleculeNumber();
975
976 if (!intermolecular && (moleculeIndex != moleculeIndexJ)) {
977 continue;
978 }
979
980
981 if (iSymm != 0 && (moleculeIndex == moleculeIndexJ)) {
982 if (domainDecomposition.isSpecialPositionExclusion(moleculeIndex, iSymm)) {
983 if (logger.isLoggable(Level.FINEST)) {
984 logger.finest(
985 format(" Excluding Interaction for Atoms %d and %d in Molecule %d for SymOp %d.",
986 atomIndex, aj, moleculeIndex, iSymm));
987 }
988 continue;
989 }
990 }
991
992 if (mask[aj] > 0 && (iSymm == 0 || aj >= atomIndex)) {
993 int aj3 = aj * 3;
994 final double xj = pair[aj3 + XX];
995 final double yj = pair[aj3 + YY];
996 final double zj = pair[aj3 + ZZ];
997 final double xr = xi - xj;
998 final double yr = yi - yj;
999 final double zr = zi - zj;
1000 final double d2 = crystal.image(xr, yr, zr);
1001 if (d2 <= cutoffPlusBuffer2) {
1002 if (logClosePairs && d2 < crystal.getSpecialPositionCutoff2()) {
1003
1004 logger.fine(
1005 format(" Close interaction (%6.3f) between atoms (iSymm = %d):\n %s\n %s\n",
1006 sqrt(d2), iSymm, atoms[atomIndex].toString(), atoms[aj].toString()));
1007 }
1008
1009
1010 try {
1011 pairs[n++] = aj;
1012 } catch (Exception e) {
1013 n = pairs.length;
1014 pairs = copyOf(pairs, n + 100);
1015 pairs[n++] = aj;
1016 }
1017 }
1018 }
1019 }
1020
1021
1022 if (iSymm == 0 && maskingRules != null) {
1023 maskingRules.removeMask(atomIndex, vdw14, mask);
1024 }
1025 }
1026
1027 }
1028
1029 private class MxNListLoop extends IntegerForLoop {
1030
1031 @Override
1032 public void run(final int lb, final int ub) throws Exception {
1033
1034 for(int i = 0; i < nSymm; i++){
1035 for(int j = lb; j < ub; j++){
1036
1037 int[] pairs = new int[1000*N];
1038 int n = 0;
1039 boolean[] exists = new boolean[nAtoms];
1040 for(int k = 0; k < M; k++){
1041 int atomID = groups[i][j][k];
1042 if (atomID == -1) {
1043 continue;
1044 }
1045 int[] neighborList = lists[i][atomID];
1046 for(int l = 0; l < neighborList.length; l++){
1047 boolean newInteraction = !exists[neighborList[l]];
1048 if(newInteraction){
1049 try{
1050 Arrays.fill(pairs, n, n+N, neighborList[l]);
1051 n+=N;
1052 } catch (Exception e){
1053 pairs = copyOf(pairs, n + 100*N);
1054 Arrays.fill(pairs, n, n+N, neighborList[l]);
1055 n+=N;
1056 }
1057 exists[neighborList[l]] = true;
1058 }
1059 }
1060 }
1061 groupLists[i][j] = copyOf(pairs, n);
1062 }
1063 }
1064 }
1065 }
1066
1067
1068
1069
1070 private static class DomainDecomposition {
1071
1072
1073
1074
1075 private Crystal crystal;
1076
1077
1078
1079
1080 private final double targetInterfacialRadius;
1081
1082
1083
1084
1085
1086
1087
1088
1089 private final int nEdge;
1090
1091
1092
1093
1094 private final int nSearch;
1095
1096 private int nA;
1097
1098 private int nB;
1099
1100
1101
1102 private int nC;
1103
1104
1105
1106 private int nAtoms;
1107
1108 private int[] cellA;
1109
1110 private int[] cellB;
1111
1112 private int[] cellC;
1113
1114
1115
1116 private Cell[][][] cells;
1117
1118
1119
1120 private Map<Integer, List<Integer>> specialPositionExclusionMap;
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130 public DomainDecomposition(int nAtoms, Crystal crystal, double cutoffPlusBuffer) {
1131 this.nEdge = 2;
1132 this.nSearch = 2 * nEdge + 1;
1133 targetInterfacialRadius = cutoffPlusBuffer / (double) nEdge;
1134 initDomainDecomposition(nAtoms, crystal);
1135 }
1136
1137
1138
1139
1140 public void initDomainDecomposition(int nAtoms, Crystal crystal) {
1141 this.nAtoms = nAtoms;
1142 this.crystal = crystal;
1143
1144
1145 if (cellA == null || cellA.length < nAtoms) {
1146 cellA = new int[nAtoms];
1147 cellB = new int[nAtoms];
1148 cellC = new int[nAtoms];
1149 }
1150
1151
1152 int maxA = (int) floor(2 * crystal.interfacialRadiusA / targetInterfacialRadius);
1153 int maxB = (int) floor(2 * crystal.interfacialRadiusB / targetInterfacialRadius);
1154 int maxC = (int) floor(2 * crystal.interfacialRadiusC / targetInterfacialRadius);
1155
1156
1157
1158
1159
1160 if (!crystal.aperiodic()) {
1161 assert (maxA >= nSearch - 1);
1162 assert (maxB >= nSearch - 1);
1163 assert (maxC >= nSearch - 1);
1164 }
1165
1166
1167
1168
1169
1170 nA = maxA < nSearch ? 1 : maxA;
1171 nB = maxB < nSearch ? 1 : maxB;
1172 nC = maxC < nSearch ? 1 : maxC;
1173
1174
1175 if (cells == null || cells.length != nA || cells[0].length != nB || cells[0][0].length != nC) {
1176 cells = new Cell[nA][nB][nC];
1177 for (int i = 0; i < nA; i++) {
1178 for (int j = 0; j < nB; j++) {
1179 for (int k = 0; k < nC; k++) {
1180 cells[i][j][k] = new Cell(i, j, k);
1181 }
1182 }
1183 }
1184 } else {
1185 clear();
1186 }
1187
1188 }
1189
1190
1191
1192
1193 public void clear() {
1194
1195 for (int i = 0; i < nA; i++) {
1196 for (int j = 0; j < nB; j++) {
1197 for (int k = 0; k < nC; k++) {
1198 cells[i][j][k].clear();
1199 }
1200 }
1201 }
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211 }
1212
1213
1214
1215
1216
1217
1218
1219 public void addSpecialPositionExclusion(int molecule, int symOpIndex) {
1220
1221 if (specialPositionExclusionMap == null) {
1222 specialPositionExclusionMap = new HashMap<>();
1223 }
1224
1225 List<Integer> list = specialPositionExclusionMap.get(molecule);
1226 if (list == null) {
1227 list = new ArrayList<>();
1228 specialPositionExclusionMap.put(molecule, list);
1229 }
1230
1231 if (!list.contains(symOpIndex)) {
1232 list.add(symOpIndex);
1233 }
1234 }
1235
1236
1237
1238
1239
1240
1241
1242
1243 public boolean isSpecialPositionExclusion(int molecule, int symOpIndex) {
1244 if (specialPositionExclusionMap == null) {
1245 return false;
1246 }
1247 List<Integer> list = specialPositionExclusionMap.get(molecule);
1248 if (list == null) {
1249 return false;
1250 }
1251 return list.contains(symOpIndex);
1252 }
1253
1254
1255
1256
1257
1258
1259
1260 public List<Integer> getSpecialPositionSymOps(int molecule) {
1261 if (specialPositionExclusionMap == null) {
1262 return null;
1263 }
1264 if (specialPositionExclusionMap.containsKey(molecule)) {
1265 List<Integer> list = specialPositionExclusionMap.get(molecule);
1266 if (list.isEmpty()) {
1267 return null;
1268 } else {
1269 return list;
1270 }
1271 }
1272 return null;
1273 }
1274
1275
1276
1277
1278 public void log() {
1279 int nCells = nA * nB * nC;
1280 int nSymm = crystal.spaceGroup.getNumberOfSymOps();
1281 StringBuilder sb = new StringBuilder(" Neighbor List Builder\n");
1282 sb.append(format(" Total Cells: %8d\n", nCells));
1283 if (nCells > 1) {
1284 sb.append(format(" Interfacial radius %8.3f A\n", targetInterfacialRadius));
1285 int searchA = nA == 1 ? 1 : nSearch;
1286 int searchB = nB == 1 ? 1 : nSearch;
1287 int searchC = nC == 1 ? 1 : nSearch;
1288 sb.append(format(" Domain Decomposition: %8d %4d %4d\n", nA, nB, nC));
1289 sb.append(format(" Neighbor Search: %8d x%3d x%3d\n", searchA, searchB, searchC));
1290 sb.append(format(" Mean Atoms per Cell: %8d", nAtoms * nSymm / nCells));
1291 }
1292 logger.info(sb.toString());
1293 }
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305 public Cell image(int i, int j, int k) {
1306 if (i >= nA) {
1307 i -= nA;
1308 } else if (i < 0) {
1309 i += nA;
1310 }
1311 if (j >= nB) {
1312 j -= nB;
1313 } else if (j < 0) {
1314 j += nB;
1315 }
1316 if (k >= nC) {
1317 k -= nC;
1318 } else if (k < 0) {
1319 k += nC;
1320 }
1321 return cells[i][j][k];
1322 }
1323
1324
1325
1326
1327
1328
1329
1330
1331 public void addAtomToCell(int i, int iSymm, double[] frac) {
1332 double xu = frac[0];
1333 double yu = frac[1];
1334 double zu = frac[2];
1335
1336 while (xu < 0.0) {
1337 xu += 1.0;
1338 }
1339 while (xu >= 1.0) {
1340 xu -= 1.0;
1341 }
1342 while (yu < 0.0) {
1343 yu += 1.0;
1344 }
1345 while (yu >= 1.0) {
1346 yu -= 1.0;
1347 }
1348 while (zu < 0.0) {
1349 zu += 1.0;
1350 }
1351 while (zu >= 1.0) {
1352 zu -= 1.0;
1353 }
1354
1355 final int a = (int) floor(xu * nA);
1356 final int b = (int) floor(yu * nB);
1357 final int c = (int) floor(zu * nC);
1358
1359 if (iSymm == 0) {
1360
1361 cellA[i] = a;
1362 cellB[i] = b;
1363 cellC[i] = c;
1364 }
1365 cells[a][b][c].add(i, iSymm, zu);
1366 }
1367
1368
1369
1370
1371
1372
1373
1374 public Cell getCellForAtom(int i) {
1375 return cells[cellA[i]][cellB[i]][cellC[i]];
1376 }
1377
1378 public Cell getCell(int a, int b, int c) {
1379 if (a < 0 || a >= nA || b < 0 || b >= nB || c < 0 || c >= nC) {
1380 return null;
1381 }
1382 return cells[a][b][c];
1383 }
1384 }
1385
1386
1387
1388
1389 public static class AtomIndex implements Comparable<AtomIndex> {
1390
1391 public final int iSymm;
1392 public final int i;
1393 public final double z;
1394
1395 public AtomIndex(int iSymm, int i, double z) {
1396 this.iSymm = iSymm;
1397 this.i = i;
1398 this.z = z;
1399 }
1400
1401 @Override
1402 public int compareTo(AtomIndex o) {
1403 return Double.compare(this.z, o.z);
1404 }
1405 }
1406
1407
1408
1409
1410 public static class Cell {
1411
1412
1413
1414
1415 final int a;
1416
1417
1418
1419 final int b;
1420
1421
1422
1423 final int c;
1424
1425
1426
1427
1428 final List<AtomIndex> list;
1429
1430
1431
1432 final Set<Integer> groupList;
1433
1434 public Cell(int a, int b, int c) {
1435 this.a = a;
1436 this.b = b;
1437 this.c = c;
1438 list = Collections.synchronizedList(new ArrayList<>());
1439 groupList = Collections.synchronizedSet(new HashSet<>());
1440 }
1441
1442
1443
1444
1445
1446
1447
1448 public void add(int atomIndex, int symOpIndex, double z) {
1449 list.add(new AtomIndex(symOpIndex, atomIndex, z));
1450 }
1451
1452 public void groupAdd(int groupIndex){
1453 groupList.add(groupIndex);
1454 }
1455
1456 public AtomIndex get(int index) {
1457 if (index >= list.size()) {
1458 return null;
1459 }
1460 return list.get(index);
1461 }
1462
1463 public int getCount() {
1464 return list.size();
1465 }
1466
1467 public int getGroupCount(){
1468 return groupList.size();
1469 }
1470
1471
1472
1473
1474
1475
1476
1477
1478 public int getSymOpAtoms(int symOpIndex, int[] index) {
1479 int count = 0;
1480 for (AtomIndex atomIndex : list) {
1481 if (atomIndex.iSymm == symOpIndex) {
1482 if (count >= index.length) {
1483 throw new IndexOutOfBoundsException(
1484 "Index out of bounds: count " + count + " " + index.length + " " + list.size());
1485 }
1486 index[count++] = atomIndex.i;
1487 }
1488 }
1489 return count;
1490 }
1491
1492
1493
1494
1495 public void clear() {
1496 list.clear();
1497 }
1498
1499 }
1500
1501
1502
1503
1504
1505 public static void main(String[] args){
1506 PotentialsUtils potentialsUtils = new PotentialsUtils();
1507 File xyzFile = new File("./examples/waterbox.xyz");
1508 if(!xyzFile.exists()){
1509 System.out.println("File does not exist");
1510 }
1511 MolecularAssembly molecularAssembly = potentialsUtils.open(xyzFile);
1512 molecularAssembly.getPotentialEnergy().energy(null, true);
1513 }
1514 }