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.xray;
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.SharedBooleanArray;
45 import ffx.crystal.Crystal;
46 import ffx.potential.bonded.Atom;
47
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50
51 import static java.lang.String.format;
52 import static java.util.Arrays.fill;
53 import static org.apache.commons.math3.util.FastMath.floor;
54 import static org.apache.commons.math3.util.FastMath.min;
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75 public class BulkSolventList extends ParallelRegion {
76
77 private static final Logger logger = Logger.getLogger(BulkSolventList.class.getName());
78
79 private final Crystal crystal;
80
81 private final int nSymm;
82
83 private final int nAtoms;
84
85 private final SharedBooleanArray[] sharedSelect;
86
87 private final int nEdgeA, nEdgeB, nEdgeC;
88
89 private final int nAB;
90
91 private final int nCells;
92
93 private final int[][] cellIndex;
94
95 private final int[] cellA;
96
97 private final int[] cellB;
98
99 private final int[] cellC;
100
101 private final int[][] cellList;
102
103
104
105
106 private final int[][] cellOffset;
107
108 private final int[][] cellCount;
109
110 private final int[][] cellStart;
111
112 private final double cutoff;
113
114 private final double cutoff2;
115
116 private final double[][] frac;
117
118 private final ParallelTeam parallelTeam;
119
120 private final SelectionListLoop[] selectionListLoop;
121
122 private final double toSeconds = 1.0e-9;
123
124 private double[][][] coordinates;
125
126 private boolean[][] selected;
127
128 private int nSelected;
129
130 private int nA;
131
132 private int nB;
133
134 private int nC;
135
136 private long time;
137 private long cellTime, selectTime, totalTime;
138
139
140
141
142
143
144
145
146
147
148 public BulkSolventList(Crystal crystal, Atom[] atoms, double cutoff, ParallelTeam parallelTeam) {
149 this.crystal = crystal;
150 this.cutoff = cutoff;
151 this.parallelTeam = new ParallelTeam(parallelTeam.getThreadCount());
152 nAtoms = atoms.length;
153 nSymm = crystal.spaceGroup.symOps.size();
154
155 cutoff2 = cutoff * cutoff;
156 final double side = min(min(crystal.a, crystal.b), crystal.c);
157
158 assert (side > 2.0 * cutoff);
159
160
161 nEdgeA = 2;
162 nEdgeB = nEdgeA;
163 nEdgeC = nEdgeA;
164 double minLengthA = cutoff / (double) nEdgeA;
165 double minLengthB = cutoff / (double) nEdgeB;
166 double minLengthC = cutoff / (double) nEdgeC;
167
168 nA = (int) floor(crystal.a / minLengthA);
169 nB = (int) floor(crystal.b / minLengthB);
170 nC = (int) floor(crystal.c / minLengthC);
171
172 if (nA < nEdgeA * 2 + 1) {
173 nA = 1;
174 }
175 if (nB < nEdgeB * 2 + 1) {
176 nB = 1;
177 }
178 if (nC < nEdgeC * 2 + 1) {
179 nC = 1;
180 }
181
182 nAB = nA * nB;
183 nCells = nAB * nC;
184
185 cellList = new int[nSymm][nAtoms];
186 cellIndex = new int[nSymm][nAtoms];
187 cellOffset = new int[nSymm][nAtoms];
188 cellStart = new int[nSymm][nCells];
189 cellCount = new int[nSymm][nCells];
190 cellA = new int[nAtoms];
191 cellB = new int[nAtoms];
192 cellC = new int[nAtoms];
193 frac = new double[3][nAtoms];
194
195 int threadCount = parallelTeam.getThreadCount();
196 selectionListLoop = new SelectionListLoop[threadCount];
197 for (int i = 0; i < threadCount; i++) {
198 selectionListLoop[i] = new SelectionListLoop();
199 }
200 sharedSelect = new SharedBooleanArray[nSymm];
201 for (int i = 0; i < nSymm; i++) {
202 sharedSelect[i] = new SharedBooleanArray(nAtoms);
203 }
204 }
205
206
207
208
209
210
211
212
213
214 public void buildList(final double[][][] coordinates, final boolean[][] selected, boolean log) {
215 this.coordinates = coordinates;
216 this.selected = selected;
217
218 cellTime = -System.nanoTime();
219 assignAtomsToCells();
220 cellTime += System.nanoTime();
221
222 selectTime = -System.nanoTime();
223 selectAtoms();
224 selectTime += System.nanoTime();
225 totalTime = cellTime + selectTime;
226
227 if (log) {
228 log();
229 }
230 }
231
232
233
234
235
236
237
238
239 @Override
240 public void finish() {
241 if (logger.isLoggable(Level.FINE)) {
242 logger.fine(
243 format("Parallel Neighbor List: %10.3f seconds", (System.nanoTime() - time) * 1e-9));
244 }
245 }
246
247
248
249
250
251
252
253
254 @Override
255 public void run() {}
256
257
258
259
260
261
262
263
264 @Override
265 public void start() {
266 time = System.nanoTime();
267 }
268
269 private void log() {
270 StringBuilder sb = new StringBuilder(format(" The cutoff is %5.2f angstroms.\n", cutoff));
271 sb.append(format(" Assignment to cells: %8.3f\n", cellTime * toSeconds));
272 sb.append(format(" Selection: %8.3f\n", selectTime * toSeconds));
273 sb.append(format(" Total: %8.3f (sec)\n", totalTime * toSeconds));
274 if (nSymm > 1) {
275 double speedup = (double) (nAtoms * nSymm) / (double) (nAtoms + nSelected);
276 sb.append(format(" Atoms in the asymmetric unit: %12d\n", nAtoms));
277 sb.append(format(" Selected symmetry mate atoms: %12d\n", nSelected));
278 sb.append(format(" Babinet speed up factor: %12.3f\n", speedup));
279 }
280 sb.append("\n");
281 logger.info(sb.toString());
282 }
283
284
285
286
287
288
289
290 private void assignAtomsToCells() {
291 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
292 final int[] cellIndexs = cellIndex[iSymm];
293 final int[] cellCounts = cellCount[iSymm];
294 final int[] cellStarts = cellStart[iSymm];
295 final int[] cellLists = cellList[iSymm];
296 final int[] cellOffsets = cellOffset[iSymm];
297
298 for (int i = 0; i < nCells; i++) {
299 cellCounts[i] = 0;
300 }
301
302 final double[][] xyz = coordinates[iSymm];
303
304 crystal.toFractionalCoordinates(nAtoms, xyz[0], xyz[1], xyz[2], frac[0], frac[1], frac[2]);
305
306
307 for (int i = 0; i < nAtoms; i++) {
308 double xu = frac[0][i];
309 double yu = frac[1][i];
310 double zu = frac[2][i];
311
312 while (xu >= 1.0) {
313 xu -= 1.0;
314 }
315 while (xu < 0.0) {
316 xu += 1.0;
317 }
318 while (yu >= 1.0) {
319 yu -= 1.0;
320 }
321 while (yu < 0.0) {
322 yu += 1.0;
323 }
324 while (zu >= 1.0) {
325 zu -= 1.0;
326 }
327 while (zu < 0.0) {
328 zu += 1.0;
329 }
330
331 final int a = (int) floor(xu * nA);
332 final int b = (int) floor(yu * nB);
333 final int c = (int) floor(zu * nC);
334 if (iSymm == 0) {
335 cellA[i] = a;
336 cellB[i] = b;
337 cellC[i] = c;
338 }
339
340 final int index = a + b * nA + c * nAB;
341 cellIndexs[i] = index;
342
343 cellOffsets[i] = cellCounts[index]++;
344 }
345
346 cellStarts[0] = 0;
347 for (int i = 1; i < nCells; i++) {
348 final int i1 = i - 1;
349 cellStarts[i] = cellStarts[i1] + cellCounts[i1];
350 }
351
352 for (int i = 0; i < nAtoms; i++) {
353 final int index = cellIndexs[i];
354 cellLists[cellStarts[index]++] = i;
355 }
356
357 cellStarts[0] = 0;
358 for (int i = 1; i < nCells; i++) {
359 final int i1 = i - 1;
360 cellStarts[i] = cellStarts[i1] + cellCounts[i1];
361 }
362 }
363 }
364
365
366
367
368
369
370 private void selectAtoms() {
371 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
372 if (selected[iSymm] == null) {
373 selected[iSymm] = new boolean[nAtoms];
374 }
375 for (int i = 0; i < nAtoms; i++) {
376 sharedSelect[iSymm].set(i, false);
377 }
378 }
379 try {
380 parallelTeam.execute(this);
381 } catch (Exception e) {
382 String message = "Fatal exception building neighbor list.\n";
383 logger.log(Level.SEVERE, message, e);
384 }
385 nSelected = 0;
386 fill(selected[0], false);
387 for (int iSymm = 1; iSymm < nSymm; iSymm++) {
388 boolean[] select = selected[iSymm];
389 SharedBooleanArray shared = sharedSelect[iSymm];
390 for (int i = 0; i < nAtoms; i++) {
391 select[i] = shared.get(i);
392 if (select[i]) {
393 nSelected++;
394 }
395 }
396 }
397 }
398
399
400
401
402
403
404
405
406 private class SelectionListLoop extends IntegerForLoop {
407
408 private final IntegerSchedule schedule;
409 private int iSymm;
410 private int atomIndex;
411 private double[][] xyz;
412 private SharedBooleanArray select;
413
414 public SelectionListLoop() {
415 super();
416 schedule = IntegerSchedule.dynamic(10);
417 }
418
419 @Override
420 public void run(final int lb, final int ub) {
421 for (iSymm = 1; iSymm < nSymm; iSymm++) {
422 select = sharedSelect[iSymm];
423
424 for (atomIndex = lb; atomIndex <= ub; atomIndex++) {
425 final int a = cellA[atomIndex];
426 final int b = cellB[atomIndex];
427 final int c = cellC[atomIndex];
428
429 int aStart = a - nEdgeA;
430 int aStop = a + nEdgeA;
431 int bStart = b - nEdgeB;
432 int bStop = b + nEdgeB;
433 int cStart = c - nEdgeC;
434 int cStop = c + nEdgeC;
435
436
437
438
439
440 if (nA == 1) {
441 aStart = a;
442 aStop = a;
443 }
444 if (nB == 1) {
445 bStart = b;
446 bStop = b;
447 }
448 if (nC == 1) {
449 cStart = c;
450 cStop = c;
451 }
452
453
454 for (int ai = aStart; ai <= aStop; ai++) {
455 for (int bi = bStart; bi <= bStop; bi++) {
456 for (int ci = cStart; ci <= cStop; ci++) {
457 selectAsymmetricAtoms(image(ai, bi, ci));
458 }
459 }
460 }
461 }
462 }
463 }
464
465 @Override
466 public IntegerSchedule schedule() {
467 return schedule;
468 }
469
470 @Override
471 public void start() {
472 xyz = coordinates[0];
473 }
474
475
476
477
478
479
480
481
482
483
484
485 private int image(int i, int j, int k) {
486 if (i >= nA) {
487 i -= nA;
488 } else if (i < 0) {
489 i += nA;
490 }
491 if (j >= nB) {
492 j -= nB;
493 } else if (j < 0) {
494 j += nB;
495 }
496 if (k >= nC) {
497 k -= nC;
498 } else if (k < 0) {
499 k += nC;
500 }
501 return i + j * nA + k * nAB;
502 }
503
504 private void selectAsymmetricAtoms(final int pairCellIndex) {
505 final double xi = xyz[0][atomIndex];
506 final double yi = xyz[1][atomIndex];
507 final double zi = xyz[2][atomIndex];
508 final int[] pairList = cellList[iSymm];
509 int start = cellStart[iSymm][pairCellIndex];
510 final int pairStop = start + cellCount[iSymm][pairCellIndex];
511 final double[][] pair = coordinates[iSymm];
512
513 for (int j = start; j < pairStop; j++) {
514 final int aj = pairList[j];
515
516 if (select.get(aj)) {
517 continue;
518 }
519 final double xj = pair[0][aj];
520 final double yj = pair[1][aj];
521 final double zj = pair[2][aj];
522 final double xr = xi - xj;
523 final double yr = yi - yj;
524 final double zr = zi - zj;
525 final double d2 = crystal.image(xr, yr, zr);
526 if (d2 <= cutoff2) {
527 select.set(aj, true);
528 }
529 }
530 }
531 }
532 }