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