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-2024.
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.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   * The NeighborList class builds Verlet lists in parallel via a spatial decomposition. <br>
69   *
70   * <ol>
71   *   <li>The unit cell is partitioned into <code>nA * nB * nC</code> smaller axis-aligned cells,
72   *       where {nA, nB, nC} are chosen as large as possible subject to the criteria that the length
73   *       of each side of a sub-volume (rCellA, rCellB, rCellC) multiplied by nEdge
74   *       must be greater than the cutoff distance <code>Rcut</code> plus a
75   *       buffer distance <code>delta</code>: <br>
76   *       <code>rCellA * nEdge .GE. (Rcut + delta)</code> <br>
77   *       <code>rCellB * nEdge .GE. (Rcut + delta)</code> <br>
78   *       <code>rCellC * nEdge .GE. (Rcut + delta)</code> <br>
79   *       All neighbors of an atom are in a block of (2*nEdge+1)^3 neighborCells.
80   *   <li>Interactions between an atom and neighbors in the asymmetric unit require only half the
81   *       neighboring cells to be searched to avoid double counting. However, enumeration of
82   *       interactions between an atom in the asymmetric unit and its neighbors in a symmetry mate
83   *       require all cells to be searched.
84   *   <li>Verlet lists from the search are stored, which reduces the number of neighbors whose
85   *       distances must be calculated by a factor of approximately: <br>
86   *       <code>(4/3*Pi*Rcut^3)/(neighborCells*Vcell)</code>
87   *       About 1/3 as many interactions are contained in the Verlet lists compared to
88   *       the total amount in the neighboring cells.
89   * </ol>
90   *
91   * @author Michael J. Schnieders
92   * @since 1.0
93   */
94  public class NeighborList extends ParallelRegion {
95  
96    /** The logger. */
97    private static final Logger logger = Logger.getLogger(NeighborList.class.getName());
98  
99    private final DomainDecomposition domainDecomposition;
100   /** The masking rules to apply when building the neighbor list. */
101   private final MaskingInterface maskingRules;
102   /** The cutoff beyond which the pairwise energy is zero. */
103   private final double cutoff;
104   /**
105    * A buffer, which is added to the cutoff distance, such that the Verlet lists do not need to be
106    * calculated for all coordinate changes.
107    */
108   private final double buffer;
109   /** The maximum squared displacement allowed before list rebuild. */
110   private final double motion2;
111   /** The sum of the cutoff + buffer. */
112   private final double cutoffPlusBuffer;
113   /** Total^2 for distance comparisons without taking a sqrt. */
114   private final double cutoffPlusBuffer2;
115   /** The ParallelTeam coordinates use of threads and their schedules. */
116   private final ParallelTeam parallelTeam;
117   /** Number of threads used by the parallelTeam. */
118   private final int threadCount;
119   /**
120    * If this flag is set, then the Verlet lists will be rebuilt even in the absence of motion.
121    */
122   private boolean forceRebuild = false;
123   /**
124    * If this flag is set then the Verlet lists will be printed to the logger.
125    */
126   private boolean print = false;
127   /**
128    * Shared Boolean used to detect list rebuild. If any thread detects motion in the MotionLoop, the
129    * value is set to true.
130    */
131   private final SharedBoolean sharedMotion;
132   /**
133    * Detect motion in the asymmetric unit.
134    */
135   private final MotionLoop[] motionLoops;
136   /**
137    * Perform some initialization tasks prior to rebuilding the list.
138    */
139   private final ListInitBarrierAction listInitBarrierAction;
140   /**
141    * Assign atoms to cells in parallel. This loop is executed once per list rebuild.
142    */
143   private final AssignAtomsToCellsLoop[] assignAtomsToCellsLoops;
144   /** A Verlet list loop for each thread. */
145   private final NeighborListLoop[] verletListLoop;
146   /**
147    * Time to check for motion.
148    */
149   private long motionTime;
150   /**
151    * Time for some initializations.
152    */
153   private long initTime;
154   /**
155    * Time to assign atoms to cells.
156    */
157   private long assignAtomsToCellsTime;
158   /**
159    * Time to build the Verlet lists.
160    */
161   private long verletListTime;
162   /** Total number of interactions. */
163   private final SharedInteger sharedCount;
164   /** Optimal pairwise ranges. */
165   private final Range[] ranges;
166   /**
167    * The crystal object defines the unit cell dimensions and space group. If replicates of the unit
168    * cell are being used, this should be the overall replicates crystal and not the small unit cell.
169    */
170   private Crystal crystal;
171   /** The number of asymmetric units in the unit cell. */
172   private int nSymm;
173   /** The number of atoms in the asymmetric unit. */
174   private int nAtoms;
175   /** Reduced coordinates for each symmetry copy. [nSymm][3*nAtoms] */
176   private double[][] coordinates;
177   /** The Cartesian coordinates of the asymmetric unit when the list was last rebuilt. */
178   private double[] previous;
179   /** The Verlet lists. [nSymm][nAtoms][nNeighbors] */
180   private int[][][] lists;
181   /** Number of interactions per atom. */
182   private int[] listCount;
183   /** Pairwise ranges for load balancing. */
184   private PairwiseSchedule pairwiseSchedule;
185   /** Number of interactions between atoms in the asymmetric unit. */
186   private int asymmetricUnitCount;
187   /** Number of interactions between atoms in the asymmetric unit and atoms in symmetry mates. */
188   private int symmetryMateCount;
189   /** Number of atoms in the asymmetric unit with interactions lists greater than 0. */
190   private int atomsWithIteractions;
191   /** Atoms being used list. */
192   private boolean[] use;
193   /** List of atoms. */
194   private Atom[] atoms;
195   /**
196    * Time building the neighbor list.
197    */
198   private long time;
199   /** Include intermolecular interactions. */
200   private boolean intermolecular = true;
201   /**
202    * If true, interactions between two inactive atoms are included in the Neighborlist. Set to true
203    * to match OpenMM behavior (i.e. interactions between two atoms with zero mass are included). Set
204    * to false for more efficient "pure Java" optimizations.
205    */
206   private final boolean includeInactivePairs = true;
207   /** Disable updates to the NeighborList; use with caution. */
208   private boolean disableUpdates = false;
209 
210   /**
211    * Constructor for the NeighborList class.
212    *
213    * @param maskingRules This parameter may be null.
214    * @param crystal Definition of the unit cell and space group.
215    * @param atoms The atoms to generate Verlet lists for.
216    * @param cutoff The cutoff distance.
217    * @param buffer The buffer distance.
218    * @param parallelTeam Specifies the parallel environment.
219    * @since 1.0
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     // Configure the neighbor cutoff and list rebuilding criteria.
232     cutoffPlusBuffer = cutoff + buffer;
233     cutoffPlusBuffer2 = cutoffPlusBuffer * cutoffPlusBuffer;
234     motion2 = (buffer / 2.0) * (buffer / 2.0);
235 
236     // Initialize parallel constructs.
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     // Initialize the neighbor list builder subcells.
253     boolean print = logger.isLoggable(Level.FINE);
254     domainDecomposition = new DomainDecomposition(nAtoms, crystal, cutoffPlusBuffer);
255     initNeighborList(print);
256   }
257 
258   /**
259    * This method can be called as necessary to build/rebuild the neighbor lists.
260    *
261    * @param coordinates The coordinates of each atom [nSymm][nAtoms*3].
262    * @param lists The neighbor lists [nSymm][nAtoms][nPairs].
263    * @param forceRebuild If true, the list is rebuilt even if no atom has moved half the buffer
264    *     size.
265    * @param use an array of boolean.
266    * @param print a boolean.
267    * @since 1.0
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    * destroy.
290    *
291    * @throws java.lang.Exception if any.
292    */
293   public void destroy() throws Exception {
294     parallelTeam.shutdown();
295   }
296 
297   /**
298    * {@inheritDoc}
299    *
300    * <p>This is method should not be called; it is invoked by Parallel Java.
301    *
302    * @since 1.0
303    */
304   @Override
305   public void finish() {
306     if (disableUpdates) {
307       return;
308     }
309     if (!forceRebuild && !sharedMotion.get()) {
310       return;
311     }
312 
313     // Collect interactions.
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     // Print the interactions counts.
332     if (print) {
333       print();
334     }
335 
336     // Update the pair-wise schedule.
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    * Returns the cutoff distance used internally by NeighborList.
356    *
357    * @return Cutoff distance in Angstroms.
358    */
359   public double getCutoff() {
360     return cutoff;
361   }
362 
363   /**
364    * If disableUpdates true, disable updating the neighbor list upon motion. Use with caution; best
365    * recommendation is to only use if all atoms have a coordinate restraint.
366    *
367    * @param disableUpdate Disable updating the neighbor list
368    */
369   public void setDisableUpdates(boolean disableUpdate) {
370     this.disableUpdates = disableUpdate;
371   }
372 
373   /**
374    * Return the Verlet list.
375    *
376    * @return The Verlet list of size [nSymm][nAtoms][nNeighbors].
377    */
378   public int[][][] getNeighborList() {
379     return lists;
380   }
381 
382   /**
383    * Getter for the field <code>pairwiseSchedule</code>.
384    *
385    * @return a {@link ffx.potential.nonbonded.PairwiseSchedule} object.
386    */
387   public PairwiseSchedule getPairwiseSchedule() {
388     return pairwiseSchedule;
389   }
390 
391   /**
392    * {@inheritDoc}
393    *
394    * <p>This is method should not be called; it is invoked by Parallel Java.
395    *
396    * @since 1.0
397    */
398   @Override
399   public void run() {
400     if (disableUpdates) {
401       return;
402     }
403 
404     int threadIndex = getThreadIndex();
405     try {
406       if (!forceRebuild) {
407         // Check for motion.
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         // Complete some initializations.
418         barrier(listInitBarrierAction);
419         // Assign atoms to cells.
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         // Build the Verlet list.
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    * The NeighborList will be re-configured, if necessary, for the supplied atom list.
445    *
446    * @param atoms A new list of atoms to operate on.
447    */
448   public void setAtoms(Atom[] atoms) {
449     this.atoms = atoms;
450     this.nAtoms = atoms.length;
451     initNeighborList(false);
452   }
453 
454   /**
455    * The NeighborList will be re-configured, if necessary, for the supplied Crystal. Changes to both
456    * unit cell parameters and number of symmetry operators are also acceptable.
457    *
458    * @param crystal A crystal defining boundary conditions and symmetry.
459    */
460   public void setCrystal(Crystal crystal) {
461     this.crystal = crystal;
462     initNeighborList(false);
463   }
464 
465   /**
466    * Setter for the field <code>intermolecular</code>.
467    *
468    * @param intermolecular If true, the NeighborList will include intermolecular interactions.
469    */
470   public void setIntermolecular(boolean intermolecular) {
471     this.intermolecular = intermolecular;
472   }
473 
474   /**
475    * {@inheritDoc}
476    *
477    * <p>This is method should not be called; it is invoked by Parallel Java.
478    *
479    * @since 1.0
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     // Set the number of symmetry operators.
497     int newNSymm = crystal.spaceGroup.symOps.size();
498     if (nSymm != newNSymm) {
499       nSymm = newNSymm;
500     }
501 
502     // Allocate memory for fractional coordinates and subcell pointers for each atom.
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     // Find the largest sphere that is enclosed by the unit cell.
512     final double sphere = min(min(crystal.interfacialRadiusA, crystal.interfacialRadiusB),
513         crystal.interfacialRadiusC);
514 
515     /*
516      Assert that the boundary conditions defined by the crystal allow use
517      of the minimum image condition.
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    * Detect if any atom has moved 1/2 the buffer size.
549    *
550    * @since 1.0
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    * Perform some initialization tasks prior to list rebuild.
578    */
579   private class ListInitBarrierAction extends BarrierAction {
580 
581     @Override
582     public void run() {
583       initTime = -System.nanoTime();
584       // Clear the list count.
585       sharedCount.set(0);
586 
587       domainDecomposition.clear();
588 
589       // Allocate memory for neighbor lists.
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    * Assign atoms to cells.
601    *
602    * @since 1.0
603    */
604   private class AssignAtomsToCellsLoop extends IntegerForLoop {
605 
606     /**
607      * The cartesian coordinates of the atom in the asymmetric unit.
608      */
609     private final double[] auCart = new double[3];
610     /**
611      * The cartesian coordinates of the atom after symmetry application.
612      */
613     private final double[] cart = new double[3];
614     /**
615      * The fractional coordinates of the atom after symmetry application.
616      */
617     private final double[] frac = new double[3];
618 
619     @Override
620     public void run(int lb, int ub) throws Exception {
621 
622       // Save the current coordinates.
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         // Assign each atom to a cell using fractional coordinates.
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               // Exclude molecule-molecule interactions between the asymmetric unit and the iSymm symmetry mate.
664               domainDecomposition.addSpecialPositionExclusion(molecule, iSymm);
665             }
666           }
667           // Convert to fractional coordinates.
668           crystal.toFractionalCoordinates(cart, frac);
669           domainDecomposition.addAtomToCell(i, iSymm, frac);
670         }
671       }
672     }
673   }
674 
675   /**
676    * The VerletListLoop class encapsulates thread local variables and methods for building Verlet
677    * lists based on a spatial decomposition of the unit cell.
678    *
679    * @author Michael J. Schnieders
680    * @since 1.0
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         // Loop over all atoms.
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              If the number of divisions is 1 in any direction then
742              set the loop limits to the current cell value.
743              Otherwise search from a - nEdge to a + nEdge.
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               // Interactions within the "self-volume".
754               atomCellPairs(cellForCurrentAtom);
755               // Search half of the neighboring volumes to avoid double counting.
756               // (a, b+1..b+nE, c)
757               for (int bi = b1; bi <= bStop; bi++) {
758                 atomCellPairs(domainDecomposition.image(a, bi, c));
759               }
760               // (a, b-nE..b+nE, c+1..c+nE)
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               // (a+1..a+nE, b-nE..b+nE, c-nE..c+nE)
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               // Interactions with all adjacent symmetry mate cells.
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       // Load the cell's atom indices into an array for the specified SymOp.
817       int num = cell.getSymOpAtoms(iSymm, pairCellAtoms);
818       if (num == 0) {
819         return;
820       }
821 
822       // Check if this pair search is over atoms in the asymmetric unit.
823       if (iSymm == 0 && maskingRules != null) {
824         // Interactions between atoms in the asymmetric unit may be masked.
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       // Loop over atoms in the "pair" cell.
839       for (int j = 0; j < num; j++) {
840         final int aj = pairCellAtoms[j];
841 
842         // If the self-volume is being searched for pairs, avoid double counting.
843         if (iSymm == 0 && cell == cellForCurrentAtom) {
844           // Only consider atoms with a higher index than the current atom.
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         // Check if this pair search is only intra-molecular.
862         if (!intermolecular && (moleculeIndex != moleculeIndexJ)) {
863           continue;
864         }
865 
866         // Check for a special position exclusion between two atoms in the same molecule.
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               // Warn about close interactions.
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             // Add the pair to the list, reallocating the array size if necessary.
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       // Interactions between atoms in the asymmetric unit may be masked.
908       if (iSymm == 0 && maskingRules != null) {
909         maskingRules.removeMask(atomIndex, vdw14, mask);
910       }
911     }
912 
913   }
914 
915   /**
916    * Break down the simulation domain into a 3D grid of cells.
917    */
918   private static class DomainDecomposition {
919 
920     /**
921      * The crystal that defines the simulation domain.
922      */
923     private Crystal crystal;
924     /**
925      * The targetInterfacialRadius is the smallest interfacial radius for a subcell, given a search
926      * of nEdge subcells along an axis, that assures all neighbors will be found.
927      */
928     private final double targetInterfacialRadius;
929     /**
930      * The number of subcells that must be searched along an axis to find all neighbors within
931      * the cutoff + buffer distance.
932      *
933      * <p>If the nEdge == 1 for {X=A,B,C} then all neighbors will be found in 3x3x3 = 27 cells.
934      * If each nEdge == 2, then all neighbors will be found in 5x5x5 = 125 cells (in this case the
935      * cells are smaller).
936      */
937     private final int nEdge;
938     /**
939      * The number of subcells that must be searched along an axis to find all neighbors within the
940      * cutoff + buffer distance. nSearch = 2 * nEdge + 1
941      */
942     private final int nSearch;
943     /** The number of divisions along the A-axis. */
944     private int nA;
945     /** The number of divisions along the B-axis. */
946     private int nB;
947     /**
948      * The number of divisions along the C-Axis.
949      */
950     private int nC;
951     /**
952      * The number of atoms.
953      */
954     private int nAtoms;
955     /** The cell indices of each atom along the A-axis. */
956     private int[] cellA;
957     /** The cell indices of each atom along the B-axis. */
958     private int[] cellB;
959     /** The cell indices of each atom along the C-axis. */
960     private int[] cellC;
961     /**
962      * The fractional sub-volumes of the unit cell.
963      */
964     private Cell[][][] cells;
965     /**
966      * Special position exclusions. Map<Molecule, List<SymOpIndex>>
967      */
968     private Map<Integer, List<Integer>> specialPositionExclusionMap;
969 
970 
971     /**
972      * DomainDecomposition Constructor.
973      *
974      * @param nAtoms Number of atoms.
975      * @param crystal The crystal.
976      * @param cutoffPlusBuffer The cutoff plus buffer distance.
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      * Initialize the domain decomposition.
987      */
988     public void initDomainDecomposition(int nAtoms, Crystal crystal) {
989       this.nAtoms = nAtoms;
990       this.crystal = crystal;
991 
992       // Allocate memory for fractional coordinates and subcell pointers for each atom.
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       // Determine the number of subcells along each axis.
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         To satisfy the cutoff + buffer distance, the number of subcells (maxA, maxB, maxC)
1006         is guaranteed to be at least nSearch - 1 for periodic crystals.
1007         */
1008       if (!crystal.aperiodic()) {
1009         assert (maxA >= nSearch - 1);
1010         assert (maxB >= nSearch - 1);
1011         assert (maxC >= nSearch - 1);
1012       }
1013 
1014       /*
1015         If the number of subcells less than nSearch, then turn off the overhead of the
1016         subcell search by setting the number of cells along that axis to 1.
1017        */
1018       nA = maxA < nSearch ? 1 : maxA;
1019       nB = maxB < nSearch ? 1 : maxB;
1020       nC = maxC < nSearch ? 1 : maxC;
1021 
1022       // Allocate memory for the subcells.
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      * Clear the sub-cells of the atom list.
1040      */
1041     public void clear() {
1042       // Clear cell contents.
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       // Clear special position exclusions.
1052       // if (specialPositionExclusionMap != null) {
1053       // Clear the lists, but don't remove the keys.
1054       //   for (List<Integer> list : specialPositionExclusionMap.values()) {
1055       //     list.clear();
1056       //   }
1057       //   specialPositionExclusionMap.clear();
1058       // }
1059     }
1060 
1061     /**
1062      * Add a special position exclusion.
1063      *
1064      * @param molecule The molecule.
1065      * @param symOpIndex The symmetry operation index.
1066      */
1067     public void addSpecialPositionExclusion(int molecule, int symOpIndex) {
1068       // Initialize the map.
1069       if (specialPositionExclusionMap == null) {
1070         specialPositionExclusionMap = new HashMap<>();
1071       }
1072       // Initialize the list for the molecule.
1073       List<Integer> list = specialPositionExclusionMap.get(molecule);
1074       if (list == null) {
1075         list = new ArrayList<>();
1076         specialPositionExclusionMap.put(molecule, list);
1077       }
1078       // Add the symOpIndex to the list.
1079       if (!list.contains(symOpIndex)) {
1080         list.add(symOpIndex);
1081       }
1082     }
1083 
1084     /**
1085      * Check if a special position exclusion exists.
1086      *
1087      * @param molecule The molecule.
1088      * @param symOpIndex The symmetry operation index.
1089      * @return True if the special position exclusion exists.
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      * Get the special position SymOps for a molecule.
1104      *
1105      * @param molecule The molecule.
1106      * @return The special position SymOps.
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      * Log information about the domain decomposition.
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      * If the index is >= to nX, it is mapped back into the periodic unit cell by subtracting nX. If
1145      * the index is less than 0, it is mapped into the periodic unit cell by adding nX. The Neighbor
1146      * list algorithm never requires multiple additions or subtractions of nX.
1147      *
1148      * @param i The index along the a-axis.
1149      * @param j The index along the b-axis.
1150      * @param k The index along the c-axis.
1151      * @return The requested Cell.
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      * Add an atom to a sub-cell.
1174      *
1175      * @param i The index of the atom.
1176      * @param iSymm The index of the symmetry operator.
1177      * @param frac The fractional coordinates of the atom.
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       // Move the atom into the range 0.0 <= x < 1.0
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       // The cell indices of this atom.
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         // Set the sub-cell indices for an asymmetric unit atom.
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      * Get the sub-cell for an asymmetric unit atom.
1218      *
1219      * @param i The index of the atom.
1220      * @return The sub-cell for the atom.
1221      */
1222     public Cell getCellForAtom(int i) {
1223       return cells[cellA[i]][cellB[i]][cellC[i]];
1224     }
1225   }
1226 
1227   /**
1228    * Hold the atom index and its symmetry operator.
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    * Hold the atoms in each cell.
1243    */
1244   public static class Cell {
1245 
1246     /**
1247      * The cell index along the A-axis.
1248      */
1249     final int a;
1250     /**
1251      * The cell index along the B-axis.
1252      */
1253     final int b;
1254     /**
1255      * The cell index along the C-axis.
1256      */
1257     final int c;
1258 
1259     /**
1260      * The list of atoms in the cell, together with their symmetry operator.
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      * Add an atom to the cell.
1273      *
1274      * @param atomIndex The atom index.
1275      * @param symOpIndex The symmetry operator index.
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      * Return the number of atoms in the cell for a given symmetry operator.
1294      *
1295      * @param symOpIndex The symmetry operator index.
1296      * @param index The list of indexes for the given symmetry operator.
1297      * @return The number of atoms in the cell for the symmetry operator.
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      * Clear the list of atoms in the cell.
1315      */
1316     public void clear() {
1317       list.clear();
1318     }
1319 
1320   }
1321 }