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