View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * The BulkSolventList class builds a list of atoms in symmetry mates that are within a cutoff
58   * distance of an atom in the asymmetric unit. This is done in parallel via a spatial decomposition.
59   *
60   * <ol>
61   *   <li>The unit cell is partitioned into <code>nA * nB * nC</code> smaller axis-aligned cells,
62   *       where nA, nB and nC are chosen as large as possible subject to the criteria that the length
63   *       of each side of a sub-volume (rCellA, rCellB, rCellC) multiplied by (nEdgeA, nEdgeB,
64   *       nEdgeC), respectively, must be greater than the cutoff distance <code>Rcut</code> plus a
65   *       buffer distance <code>delta</code>: <br>
66   *       <code>rCellA * nEdgeA .GE. (Rcut + delta)</code> <br>
67   *       <code>rCellB * nEdgeB .GE. (Rcut + delta)</code> <br>
68   *       <code>rCellC * nEdgeC .GE. (Rcut + delta)</code> All neighbors of an atom are in a block of
69   *       (2*nEdgeA+1)(2*nEdgeB+1)(2*nEdgeC+1) neighborCells.
70   * </ol>
71   *
72   * @author Michael J. Schnieders
73   * @since 1.0
74   */
75  public class BulkSolventList extends ParallelRegion {
76  
77    private static final Logger logger = Logger.getLogger(BulkSolventList.class.getName());
78    /** The crystal object defines the unit cell dimensions and spacegroup. */
79    private final Crystal crystal;
80    /** The number of asymmetric units in the unit cell. */
81    private final int nSymm;
82    /** The number of atoms in the asymmetric unit. */
83    private final int nAtoms;
84    /** The selected atoms. [nsymm][natom] */
85    private final SharedBooleanArray[] sharedSelect;
86  
87    private final int nEdgeA, nEdgeB, nEdgeC;
88    /** The number of cells in one plane (nDivisions^2). */
89    private final int nAB;
90    /** The number of cells (nDivisions^3). */
91    private final int nCells;
92    /** A temporary array that holds the index of the cell each atom is assigned to. */
93    private final int[][] cellIndex;
94    /** The cell indices of each atom along a A-axis. */
95    private final int[] cellA;
96    /** The cell indices of each atom along a B-axis. */
97    private final int[] cellB;
98    /** The cell indices of each atom along a C-axis. */
99    private final int[] cellC;
100   /** The list of atoms in each cell. [nsymm][natom] = atom index */
101   private final int[][] cellList;
102   /**
103    * The offset of each atom from the start of the cell. The first atom atom in the cell has 0
104    * offset. [nsymm][natom] = offset of the atom
105    */
106   private final int[][] cellOffset;
107   /** The number of atoms in each cell. [nsymm][ncell] */
108   private final int[][] cellCount;
109   /** The index of the first atom in each cell. [nsymm][ncell] */
110   private final int[][] cellStart;
111   /** The cutoff beyound which the pairwise energy is zero. */
112   private final double cutoff;
113   /** For distance comparisons without taking a sqrt. */
114   private final double cutoff2;
115   /** The array of fractional "a", "b", and "c" coordinates. */
116   private final double[][] frac;
117   /** The ParallelTeam coordinates use of threads and their schedules. */
118   private final ParallelTeam parallelTeam;
119   /** A Verlet list loop for each thread. */
120   private final SelectionListLoop[] selectionListLoop;
121 
122   private final double toSeconds = 1.0e-9;
123   /** Reduced coordinates for each symmetry copy. [nsymm][3][natom] */
124   private double[][][] coordinates;
125 
126   private boolean[][] selected;
127   /** Number of selected atoms. */
128   private int nSelected;
129   /** The number of divisions along the A-axis. */
130   private int nA;
131   /** The number of divisions along the B-axis. */
132   private int nB;
133   /** The number of divisions along the C-Axis. */
134   private int nC;
135 
136   private long time;
137   private long cellTime, selectTime, totalTime;
138 
139   /**
140    * Constructor for the NeighborList class.
141    *
142    * @param crystal Definition of the unit cell and space group.
143    * @param atoms The atoms to generate Verlet lists for.
144    * @param cutoff The cutoff distance.
145    * @param parallelTeam Specifies the parallel environment.
146    * @since 1.0
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     // nEdgeA, nEdgeB and nEdgeC must be >= 1.
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     // Parallel constructs.
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    * This method can be called as necessary to build/rebuild the neighbor lists.
208    *
209    * @param coordinates The coordinates of each atom [nSymm][nAtoms*3].
210    * @param selected The list of selected atoms [nSymm][nAtoms].
211    * @param log a boolean.
212    * @since 1.0
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    * {@inheritDoc}
234    *
235    * <p>This is method should not be called; it is invoked by Parallel Java.
236    *
237    * <p>since 0.1
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    * {@inheritDoc}
249    *
250    * <p>This is method should not be called; it is invoked by Parallel Java.
251    *
252    * @since 1.0
253    */
254   @Override
255   public void run() {}
256 
257   /**
258    * {@inheritDoc}
259    *
260    * <p>This is method should not be called; it is invoked by Parallel Java.
261    *
262    * @since 1.0
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    * Assign asymmetric and symmetry mate atoms to cells. This is very fast; there is little to be
286    * gained from parallelizing it at this point.
287    *
288    * @since 1.0
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       // Zero out the cell counts.
298       for (int i = 0; i < nCells; i++) {
299         cellCounts[i] = 0;
300       }
301       // Convert to fractional coordinates.
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       // Assign each atom to a cell using fractional coordinates.
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         // Move the atom into the range 0.0 <= x < 1.0
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         // The cell indices of this atom.
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         // The cell index of this atom.
340         final int index = a + b * nA + c * nAB;
341         cellIndexs[i] = index;
342         // The offset of this atom from the beginning of the cell.
343         cellOffsets[i] = cellCounts[index]++;
344       }
345       // Define the starting indices.
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       // Move atom locations into a list ordered by cell.
352       for (int i = 0; i < nAtoms; i++) {
353         final int index = cellIndexs[i];
354         cellLists[cellStarts[index]++] = i;
355       }
356       // Define the starting indices again.
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    * Execute the parallel atom selection builder.
367    *
368    * @since 1.0
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    * The SelectionListLoop class encapsulates thread local variables and methods for selecting atoms
401    * based on a spatial decomposition of the unit cell.
402    *
403    * @author Michael J. Schnieders
404    * @since 1.0
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         // Loop over all atoms.
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            * If the number of divisions is 1 in any direction then set
438            * the loop limits to the current cell value.
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           // Check atoms in symmetry mate cells.
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      * If the index is >= to nX, it is mapped back into the periodic unit cell by subtracting nX. If
477      * the index is < 0, it is mapped into the periodic unit cell by adding nX. The Neighbor list
478      * algorithm never requires multiple additions or subtractions of nX.
479      *
480      * @param i The index along the a-axis.
481      * @param j The index along the b-axis.
482      * @param k The index along the c-axis.
483      * @return The pointer into the 1D cell array.
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       // Loop over atoms in the pair cell.
513       for (int j = start; j < pairStop; j++) {
514         final int aj = pairList[j];
515         // If this atom is already selected, continue.
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 }