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-2021.
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 static ffx.numerics.math.ScalarMath.mod;
41  import static ffx.numerics.fft.Complex3D.iComplex3D;
42  import static ffx.numerics.spline.UniformBSpline.bSpline;
43  import static ffx.numerics.spline.UniformBSpline.bSplineDerivatives;
44  import static ffx.potential.parameters.MultipoleType.t000;
45  import static ffx.potential.parameters.MultipoleType.t001;
46  import static ffx.potential.parameters.MultipoleType.t002;
47  import static ffx.potential.parameters.MultipoleType.t003;
48  import static ffx.potential.parameters.MultipoleType.t010;
49  import static ffx.potential.parameters.MultipoleType.t011;
50  import static ffx.potential.parameters.MultipoleType.t012;
51  import static ffx.potential.parameters.MultipoleType.t020;
52  import static ffx.potential.parameters.MultipoleType.t021;
53  import static ffx.potential.parameters.MultipoleType.t030;
54  import static ffx.potential.parameters.MultipoleType.t100;
55  import static ffx.potential.parameters.MultipoleType.t101;
56  import static ffx.potential.parameters.MultipoleType.t102;
57  import static ffx.potential.parameters.MultipoleType.t110;
58  import static ffx.potential.parameters.MultipoleType.t111;
59  import static ffx.potential.parameters.MultipoleType.t120;
60  import static ffx.potential.parameters.MultipoleType.t200;
61  import static ffx.potential.parameters.MultipoleType.t201;
62  import static ffx.potential.parameters.MultipoleType.t210;
63  import static ffx.potential.parameters.MultipoleType.t300;
64  import static java.lang.String.format;
65  import static java.lang.System.arraycopy;
66  import static org.apache.commons.math3.util.FastMath.PI;
67  import static org.apache.commons.math3.util.FastMath.cos;
68  import static org.apache.commons.math3.util.FastMath.exp;
69  import static org.apache.commons.math3.util.FastMath.floor;
70  import static org.apache.commons.math3.util.FastMath.max;
71  import static org.apache.commons.math3.util.FastMath.min;
72  import static org.apache.commons.math3.util.FastMath.pow;
73  import static org.apache.commons.math3.util.FastMath.round;
74  import static org.apache.commons.math3.util.FastMath.sin;
75  import static org.apache.commons.math3.util.FastMath.sqrt;
76  
77  import edu.rit.pj.IntegerForLoop;
78  import edu.rit.pj.IntegerSchedule;
79  import edu.rit.pj.ParallelRegion;
80  import edu.rit.pj.ParallelTeam;
81  import ffx.crystal.Crystal;
82  import ffx.numerics.fft.Complex;
83  import ffx.numerics.fft.Complex3DParallel;
84  import ffx.potential.bonded.Atom;
85  import ffx.potential.parameters.ForceField;
86  import java.nio.DoubleBuffer;
87  import java.util.logging.Level;
88  import java.util.logging.Logger;
89  import org.apache.commons.configuration2.CompositeConfiguration;
90  
91  /**
92   * The Reciprocal Space class computes the reciprocal space contribution to {@link ParticleMeshEwald}
93   * for the AMOEBA force field.
94   *
95   * <ol>
96   *   <li>Assignment of polarizable multipole charge density to the 3D grid, via b-Splines, is
97   *       parallelized using a spatial decomposition.
98   *   <li>The convolution depends on methods of the {@link ffx.numerics.fft.Real3DParallel} and
99   *       {@link ffx.numerics.fft.Complex3DParallel} classes.
100  *   <li>Finally, the electric potential and its gradients are collected, in parallel, off the grid
101  *       using b-Splines.
102  * </ol>
103  *
104  * @author Michael J. Schnieders
105  * @since 1.0
106  */
107 public class ReciprocalSpace {
108 
109   private static final Logger logger = Logger.getLogger(ReciprocalSpace.class.getName());
110   /** First lookup index to pack a 2D tensor into a 1D array. */
111   private static final int[] qi1 = {0, 1, 2, 0, 0, 1};
112   /** Second lookup index to pack a 2D tensor into a 1D array. */
113   private static final int[] qi2 = {0, 1, 2, 1, 2, 2};
114 
115   private static final double toSeconds = 1.0e-9;
116   private final ParticleMeshEwald particleMeshEwald;
117   private final ForceField forceField;
118   /** Number of unit cell symmetry operators. */
119   private final int nSymm;
120   /** Number of threads. */
121   private final int threadCount;
122   /** The b-Spline order to use for discretization to/from the reciprocal grid. */
123   private final int bSplineOrder;
124   /**
125    * Three derivatives of the potential are needed for AMOEBA. Specifically, the field gradient is
126    * used to compute the energy of the quadrupole moment and the 2nd field gradient (3rd derivative
127    * of the reciprocal potential) the energy gradient.
128    */
129   private final int derivOrder = 3;
130   /** Ewald convergence parameter. */
131   private final double aEwald;
132   /** Parallel team instance. */
133   private final ParallelTeam parallelTeam;
134   private final ParallelTeam fftTeam;
135   private final BSplineRegion bSplineRegion;
136   private final SpatialPermanentLoop[] spatialPermanentLoops;
137   private final SpatialInducedLoop[] spatialInducedLoops;
138   private final SlicePermanentLoop[] slicePermanentLoops;
139   private final SliceInducedLoop[] sliceInducedLoops;
140   private final RowPermanentLoop[] rowPermanentLoops;
141   private final RowInducedLoop[] rowInducedLoops;
142   /**
143    * Number of atoms for a given symmetry operator that a given thread is responsible for applying to
144    * the FFT grid. gridAtomCount[nSymm][nThread]
145    */
146   private final int[][] gridAtomCount;
147   /**
148    * Atom indices for a given symmetry operator that a given thread is responsible for applying to
149    * the FFT grid. gridAtomCount[nSymm][nThread][nAtoms]
150    */
151   private final int[][][] gridAtomList;
152   private final PermanentPhiRegion permanentPhiRegion;
153   private final InducedPhiRegion polarizationPhiRegion;
154   private final IntegerSchedule recipSchedule;
155   /** Convolution variables. */
156   private final double[][] transformFieldMatrix = new double[10][10];
157   private final double[][] transformMultipoleMatrix = new double[10][10];
158   private final double[][] a = new double[3][3];
159   /**
160    * The unit cell for the simulation. A ReplicatesCrystal will give equivalent results up to the
161    * limits of double precision, but is more expensive.
162    */
163   private Crystal crystal;
164   /** Array of atoms in the asymmetric unit. */
165   private Atom[] atoms;
166   /** Number of atoms in the asymmetric unit. */
167   private int nAtoms;
168   /** The X-dimension of the FFT grid. */
169   private int fftX;
170   /** The Y-dimension of the FFT grid. */
171   private int fftY;
172   /** The Z-dimension of the FFT grid. */
173   private int fftZ;
174   /** Number of doubles needed to compute a complex to complex 3D FFT (fftX * fftY * fftZ * 2). */
175   private int fftSpace;
176   /** Reciprocal space grid. [fftSpace] */
177   private double[] splineGrid;
178   /** Wraps the splineGrid. */
179   private DoubleBuffer splineBuffer;
180   /** Reference to atomic coordinate array. */
181   private double[][][] coordinates;
182   // Spline moments and perform the convolution.
183   private SpatialDensityRegion spatialDensityRegion;
184   private SliceRegion sliceRegion;
185   private RowRegion rowRegion;
186   private Complex3DParallel complex3DFFT;
187   private GridMethod gridMethod;
188   // Timing variables.
189   private long bSplineTotal, splinePermanentTotal, splineInducedTotal;
190   private long permanentPhiTotal, inducedPhiTotal, convTotal;
191   /** Timing variables. */
192   private final long[] bSplineTime;
193   private final long[] splinePermanentTime;
194   private final long[] permanentPhiTime;
195   private final long[] splineInducedTime;
196   private final int[] splineCount;
197   private final long[] inducedPhiTime;
198 
199 
200   /**
201    * Reciprocal Space PME contribution.
202    *
203    * @param particleMeshEwald a {@link ParticleMeshEwald} object.
204    * @param crystal a {@link ffx.crystal.Crystal} object.
205    * @param forceField a {@link ffx.potential.parameters.ForceField} object.
206    * @param atoms an array of {@link ffx.potential.bonded.Atom} objects.
207    * @param aewald the Ewald parameter.
208    * @param fftTeam a {@link edu.rit.pj.ParallelTeam} object.
209    * @param parallelTeam a {@link edu.rit.pj.ParallelTeam} object.
210    */
211   public ReciprocalSpace(
212       ParticleMeshEwald particleMeshEwald,
213       Crystal crystal,
214       ForceField forceField,
215       Atom[] atoms,
216       double aewald,
217       ParallelTeam fftTeam,
218       ParallelTeam parallelTeam) {
219 
220     this.particleMeshEwald = particleMeshEwald;
221     this.crystal = crystal.getUnitCell();
222     this.forceField = forceField;
223     this.atoms = atoms;
224     this.nAtoms = atoms.length;
225     this.aEwald = aewald;
226     this.fftTeam = fftTeam;
227     this.parallelTeam = parallelTeam;
228 
229     coordinates = particleMeshEwald.getCoordinates();
230     threadCount = parallelTeam.getThreadCount();
231     nSymm = this.crystal.spaceGroup.getNumberOfSymOps();
232 
233     // Construct the parallel convolution object.
234     boolean available;
235     String recipStrategy = null;
236     try {
237       recipStrategy = forceField.getString("RECIP_SCHEDULE");
238       IntegerSchedule.parse(recipStrategy);
239       available = true;
240     } catch (Exception e) {
241       available = false;
242     }
243     if (available) {
244       recipSchedule = IntegerSchedule.parse(recipStrategy);
245       logger.log(Level.INFO, "  Convolution schedule {0}", recipStrategy);
246     } else {
247       recipSchedule = IntegerSchedule.fixed();
248     }
249 
250     CompositeConfiguration properties = forceField.getProperties();
251     String gridString = properties.getString("grid-method", "SPATIAL").toUpperCase();
252     try {
253       gridMethod = GridMethod.valueOf(gridString);
254     } catch (Exception e) {
255       gridMethod = GridMethod.SPATIAL;
256     }
257 
258     bSplineOrder = forceField.getInteger("PME_ORDER", 5);
259 
260     // Initialize convolution objects that may be re-allocated during NPT simulations.
261     double density = initConvolution();
262 
263     // Log PME reciprocal space parameters.
264     if (logger.isLoggable(Level.INFO)) {
265       StringBuilder sb = new StringBuilder();
266       sb.append(format("    B-Spline Order:                    %8d\n", bSplineOrder));
267       sb.append(format("    Mesh Density:                      %8.3f\n", density));
268       sb.append(format("    Mesh Dimensions:              (%3d,%3d,%3d)\n", fftX, fftY, fftZ));
269       sb.append(format("    Grid Method:                       %8s\n", gridMethod.toString()));
270       logger.info(sb.toString());
271     }
272 
273     // Construct the parallel BSplineRegion, DensityLoops and Phi objects.
274     bSplineRegion = new BSplineRegion();
275     switch (gridMethod) {
276       case SPATIAL:
277         spatialPermanentLoops = new SpatialPermanentLoop[threadCount];
278         spatialInducedLoops = new SpatialInducedLoop[threadCount];
279         for (int i = 0; i < threadCount; i++) {
280           spatialPermanentLoops[i] = new SpatialPermanentLoop(spatialDensityRegion, bSplineRegion);
281           spatialInducedLoops[i] = new SpatialInducedLoop(spatialDensityRegion, bSplineRegion);
282         }
283         slicePermanentLoops = null;
284         sliceInducedLoops = null;
285         rowPermanentLoops = null;
286         rowInducedLoops = null;
287         gridAtomCount = null;
288         gridAtomList = null;
289         break;
290 
291       case ROW:
292         rowPermanentLoops = new RowPermanentLoop[threadCount];
293         rowInducedLoops = new RowInducedLoop[threadCount];
294         for (int i = 0; i < threadCount; i++) {
295           rowPermanentLoops[i] = new RowPermanentLoop(rowRegion, bSplineRegion);
296           rowInducedLoops[i] = new RowInducedLoop(rowRegion, bSplineRegion);
297         }
298         gridAtomCount = new int[nSymm][threadCount];
299         gridAtomList = new int[nSymm][threadCount][nAtoms];
300         spatialPermanentLoops = null;
301         spatialInducedLoops = null;
302         slicePermanentLoops = null;
303         sliceInducedLoops = null;
304         break;
305 
306       case SLICE:
307       default:
308         slicePermanentLoops = new SlicePermanentLoop[threadCount];
309         sliceInducedLoops = new SliceInducedLoop[threadCount];
310         for (int i = 0; i < threadCount; i++) {
311           slicePermanentLoops[i] = new SlicePermanentLoop(sliceRegion, bSplineRegion);
312           sliceInducedLoops[i] = new SliceInducedLoop(sliceRegion, bSplineRegion);
313         }
314         gridAtomCount = new int[nSymm][threadCount];
315         gridAtomList = new int[nSymm][threadCount][nAtoms];
316         spatialPermanentLoops = null;
317         spatialInducedLoops = null;
318         rowPermanentLoops = null;
319         rowInducedLoops = null;
320     }
321     permanentPhiRegion = new PermanentPhiRegion(bSplineRegion);
322     polarizationPhiRegion = new InducedPhiRegion(bSplineRegion);
323 
324     // Initialize timing variables.
325     bSplineTime = new long[threadCount];
326     splinePermanentTime = new long[threadCount];
327     splineInducedTime = new long[threadCount];
328     splineCount = new int[threadCount];
329     permanentPhiTime = new long[threadCount];
330     inducedPhiTime = new long[threadCount];
331   }
332 
333   /**
334    * Computes the modulus of the discrete Fourier Transform of "bsarray" and stores it in "bsmod".
335    *
336    * @param bsmod B-Spline modulus.
337    * @param bsarray B-Spline array.
338    * @param nfft FFT dimension.
339    * @param order B-Spline order.
340    */
341   private static void discreteFTMod(double[] bsmod, double[] bsarray, int nfft, int order) {
342     // Get the modulus of the discrete Fourier fft.
343     double factor = 2.0 * PI / nfft;
344     for (int i = 0; i < nfft; i++) {
345       double sum1 = 0.0;
346       double sum2 = 0.0;
347       for (int j = 0; j < nfft; j++) {
348         double arg = factor * (i * j);
349         sum1 = sum1 + bsarray[j] * cos(arg);
350         sum2 = sum2 + bsarray[j] * sin(arg);
351       }
352       bsmod[i] = sum1 * sum1 + sum2 * sum2;
353     }
354 
355     // Fix for exponential Euler spline interpolation failure.
356     double eps = 1.0e-7;
357     if (bsmod[0] < eps) {
358       bsmod[0] = 0.5 * bsmod[1];
359     }
360     for (int i = 1; i < nfft - 1; i++) {
361       if (bsmod[i] < eps) {
362         bsmod[i] = 0.5 * (bsmod[i - 1] + bsmod[i + 1]);
363       }
364     }
365     if (bsmod[nfft - 1] < eps) {
366       bsmod[nfft - 1] = 0.5 * bsmod[nfft - 2];
367     }
368 
369     // Compute and apply the optimal zeta coefficient.
370     int jcut = 50;
371     int order2 = 2 * order;
372     for (int i = 0; i < nfft; i++) {
373       int k = i;
374       double zeta;
375       if (i > nfft / 2) {
376         k = k - nfft;
377       }
378       if (k == 0) {
379         zeta = 1.0;
380       } else {
381         double sum1 = 1.0;
382         double sum2 = 1.0;
383         factor = PI * k / nfft;
384         for (int j = 0; j < jcut; j++) {
385           double arg = factor / (factor + PI * (j + 1));
386           sum1 = sum1 + pow(arg, order);
387           sum2 = sum2 + pow(arg, order2);
388         }
389         for (int j = 0; j < jcut; j++) {
390           double arg = factor / (factor - PI * (j + 1));
391           sum1 = sum1 + pow(arg, order);
392           sum2 = sum2 + pow(arg, order2);
393         }
394         zeta = sum2 / sum1;
395       }
396       bsmod[i] = bsmod[i] * zeta * zeta;
397     }
398   }
399 
400   /**
401    * cartToFracInducedDipoles
402    *
403    * @param inducedDipole Input Cartesian induced dipole.
404    * @param inducedDipoleCR Input Cartesian induced dipole CR.
405    * @param fracInducedDipole Output fractional induced dipole.
406    * @param fracInducedDipoleCR Output fractional induced dipole CR.
407    */
408   public void cartToFracInducedDipole(
409       double[] inducedDipole,
410       double[] inducedDipoleCR,
411       double[] fracInducedDipole,
412       double[] fracInducedDipoleCR) {
413     double[] in = inducedDipole;
414     fracInducedDipole[0] = a[0][0] * in[0] + a[0][1] * in[1] + a[0][2] * in[2];
415     fracInducedDipole[1] = a[1][0] * in[0] + a[1][1] * in[1] + a[1][2] * in[2];
416     fracInducedDipole[2] = a[2][0] * in[0] + a[2][1] * in[1] + a[2][2] * in[2];
417     in = inducedDipoleCR;
418     fracInducedDipoleCR[0] = a[0][0] * in[0] + a[0][1] * in[1] + a[0][2] * in[2];
419     fracInducedDipoleCR[1] = a[1][0] * in[0] + a[1][1] * in[1] + a[1][2] * in[2];
420     fracInducedDipoleCR[2] = a[2][0] * in[0] + a[2][1] * in[1] + a[2][2] * in[2];
421   }
422 
423   /** computeBSplines */
424   public void computeBSplines() {
425     try {
426       bSplineTotal -= System.nanoTime();
427       parallelTeam.execute(bSplineRegion);
428       bSplineTotal += System.nanoTime();
429     } catch (Exception e) {
430       String message = " Fatal exception evaluating b-Splines.";
431       logger.log(Level.SEVERE, message, e);
432     }
433   }
434 
435   /**
436    * computeInducedPhi
437    *
438    * @param cartInducedDipolePhi an array of double.
439    * @param cartInducedDipoleCRPhi an array of double.
440    * @param fracInducedDipolePhi an array of double.
441    * @param fracInducedDipoleCRPhi an array of double.
442    */
443   public void computeInducedPhi(
444       double[][] cartInducedDipolePhi, double[][] cartInducedDipoleCRPhi,
445       double[][] fracInducedDipolePhi, double[][] fracInducedDipoleCRPhi) {
446     inducedPhiTotal -= System.nanoTime();
447     try {
448       polarizationPhiRegion.setInducedDipolePhi(
449           cartInducedDipolePhi, cartInducedDipoleCRPhi,
450           fracInducedDipolePhi, fracInducedDipoleCRPhi);
451       parallelTeam.execute(polarizationPhiRegion);
452     } catch (Exception e) {
453       String message = "Fatal exception evaluating induced reciprocal space potential.";
454       logger.log(Level.SEVERE, message, e);
455     }
456     inducedPhiTotal += System.nanoTime();
457   }
458 
459   /**
460    * Compute the potential Phi and its derivatives for all atoms.
461    *
462    * @param cartPermanentPhi an array of double.
463    * @param fracPermanentPhi an array of double.
464    */
465   public void computePermanentPhi(double[][] cartPermanentPhi, double[][] fracPermanentPhi) {
466     permanentPhiTotal -= System.nanoTime();
467     try {
468       permanentPhiRegion.setPermanentPhi(cartPermanentPhi, fracPermanentPhi);
469       parallelTeam.execute(permanentPhiRegion);
470     } catch (Exception e) {
471       String message = " Fatal exception evaluating permanent reciprocal space potential.";
472       logger.log(Level.SEVERE, message, e);
473     }
474     permanentPhiTotal += System.nanoTime();
475   }
476 
477   /**
478    * getXDim
479    *
480    * @return a double.
481    */
482   public int getXDim() {
483     return fftX;
484   }
485 
486   /**
487    * getYDim
488    *
489    * @return a double.
490    */
491   public int getYDim() {
492     return fftY;
493   }
494 
495   /**
496    * getZDim
497    *
498    * @return a double.
499    */
500   public int getZDim() {
501     return fftZ;
502   }
503 
504   /**
505    * globalToFracDipole.
506    *
507    * @param globalDipole an array of {@link double} objects.
508    * @return an array of {@link double} objects.
509    */
510   public double[] globalToFracDipole(double[] globalDipole) {
511     final double[] fracDipole = new double[3];
512     // Dipole
513     for (int j = 0; j < 3; j++) {
514       fracDipole[j] = 0.0;
515       for (int k = 0; k < 3; k++) {
516         fracDipole[j] = fracDipole[j] + transformMultipoleMatrix[j + 1][k + 1] * globalDipole[k];
517       }
518     }
519     return fracDipole;
520   }
521 
522   /**
523    * globalToFracMultipole.
524    *
525    * @param globalMultipole an array of {@link double} objects.
526    * @return an array of {@link double} objects.
527    */
528   public double[] globalToFracMultipole(double[] globalMultipole) {
529     double[] fracMultipole = new double[10];
530     // Charge
531     fracMultipole[0] = globalMultipole[0];
532     // Dipole
533     for (int j = 1; j < 4; j++) {
534       fracMultipole[j] = 0.0;
535       for (int k = 1; k < 4; k++) {
536         fracMultipole[j] = fracMultipole[j] + transformMultipoleMatrix[j][k] * globalMultipole[k];
537       }
538     }
539     // Quadrupole
540     for (int j = 4; j < 10; j++) {
541       fracMultipole[j] = 0.0;
542       for (int k = 4; k < 7; k++) {
543         fracMultipole[j] = fracMultipole[j] + transformMultipoleMatrix[j][k] * globalMultipole[k];
544       }
545       for (int k = 7; k < 10; k++) {
546         fracMultipole[j] =
547             fracMultipole[j] + transformMultipoleMatrix[j][k] * 2.0 * globalMultipole[k];
548       }
549       // Apply the oneThird factor for quadrupole components.
550       fracMultipole[j] = fracMultipole[j] / 3.0;
551     }
552     return fracMultipole;
553   }
554 
555   /**
556    * Note that the Java function "signum" and the FORTRAN version have different definitions for an
557    * argument of zero.
558    *
559    * <p>JAVA: Math.signum(0.0) == 0.0
560    *
561    * <p>FORTRAN: signum(0.0) .eq. 1.0
562    */
563   public void inducedDipoleConvolution() {
564     convTotal -= System.nanoTime();
565     try {
566       complex3DFFT.convolution(splineGrid);
567     } catch (Exception e) {
568       String message = "Fatal exception evaluating induced convolution.";
569       logger.log(Level.SEVERE, message, e);
570     }
571     convTotal += System.nanoTime();
572   }
573 
574   /** initTimings. */
575   public void initTimings() {
576     // Reset total timings.
577     bSplineTotal = 0;
578     splinePermanentTotal = 0;
579     splineInducedTotal = 0;
580     permanentPhiTotal = 0;
581     inducedPhiTotal = 0;
582     convTotal = 0;
583 
584     // Reset timing arrays.
585     for (int i = 0; i < threadCount; i++) {
586       bSplineTime[i] = 0;
587       splinePermanentTime[i] = 0;
588       splineInducedTime[i] = 0;
589       permanentPhiTime[i] = 0;
590       inducedPhiTime[i] = 0;
591       splineCount[i] = 0;
592     }
593 
594     // Reset 3D convolution timing.
595     if (complex3DFFT != null) {
596       complex3DFFT.initTiming();
597     }
598   }
599 
600   /** permanentMultipoleConvolution */
601   public void permanentMultipoleConvolution() {
602     convTotal -= System.nanoTime();
603     try {
604       complex3DFFT.convolution(splineGrid);
605     } catch (Exception e) {
606       String message = " Fatal exception evaluating permanent convolution.";
607       logger.log(Level.SEVERE, message, e);
608     }
609     convTotal += System.nanoTime();
610   }
611 
612   /** printTimings. */
613   public void printTimings() {
614     if (logger.isLoggable(Level.FINE)) {
615       if (complex3DFFT != null) {
616         double total =
617             (bSplineTotal
618                 + convTotal
619                 + splinePermanentTotal
620                 + permanentPhiTotal
621                 + splineInducedTotal
622                 + inducedPhiTotal)
623                 * toSeconds;
624 
625         logger.fine(String.format("\n Reciprocal Space: %7.4f (sec)", total));
626         long[] convTime = complex3DFFT.getTimings();
627         logger.fine("                           Direct Field    SCF Field");
628         logger.fine(" Thread  B-Spline  3DConv  Spline  Phi     Spline  Phi      Count");
629 
630         long minBSpline = Long.MAX_VALUE;
631         long maxBSpline = 0;
632         long minConv = Long.MAX_VALUE;
633         long maxConv = 0;
634         long minBSPerm = Long.MAX_VALUE;
635         long maxBSPerm = 0;
636         long minPhiPerm = Long.MAX_VALUE;
637         long maxPhiPerm = 0;
638         long minBSInduced = Long.MAX_VALUE;
639         long maxBSInduced = 0;
640         long minPhiInduced = Long.MAX_VALUE;
641         long maxPhiInduced = 0;
642         int minCount = Integer.MAX_VALUE;
643         int maxCount = 0;
644 
645         for (int i = 0; i < threadCount; i++) {
646           logger.fine(
647               String.format(
648                   "    %3d   %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f  %6d",
649                   i,
650                   bSplineTime[i] * toSeconds,
651                   convTime[i] * toSeconds,
652                   splinePermanentTime[i] * toSeconds,
653                   permanentPhiTime[i] * toSeconds,
654                   splineInducedTime[i] * toSeconds,
655                   inducedPhiTime[i] * toSeconds,
656                   splineCount[i]));
657           minBSpline = min(bSplineTime[i], minBSpline);
658           maxBSpline = max(bSplineTime[i], maxBSpline);
659           minConv = min(convTime[i], minConv);
660           maxConv = max(convTime[i], maxConv);
661           minBSPerm = min(splinePermanentTime[i], minBSPerm);
662           maxBSPerm = max(splinePermanentTime[i], maxBSPerm);
663           minPhiPerm = min(permanentPhiTime[i], minPhiPerm);
664           maxPhiPerm = max(permanentPhiTime[i], maxPhiPerm);
665           minBSInduced = min(splineInducedTime[i], minBSInduced);
666           maxBSInduced = max(splineInducedTime[i], maxBSInduced);
667           minPhiInduced = min(inducedPhiTime[i], minPhiInduced);
668           maxPhiInduced = max(inducedPhiTime[i], maxPhiInduced);
669           minCount = min(splineCount[i], minCount);
670           maxCount = max(splineCount[i], maxCount);
671         }
672         logger.fine(
673             String.format(
674                 " Min      %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f  %6d",
675                 minBSpline * toSeconds,
676                 minConv * toSeconds,
677                 minBSPerm * toSeconds,
678                 minPhiPerm * toSeconds,
679                 minBSInduced * toSeconds,
680                 minPhiInduced * toSeconds,
681                 minCount));
682         logger.fine(
683             String.format(
684                 " Max      %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f  %6d",
685                 maxBSpline * toSeconds,
686                 maxConv * toSeconds,
687                 maxBSPerm * toSeconds,
688                 maxPhiPerm * toSeconds,
689                 maxBSInduced * toSeconds,
690                 maxPhiInduced * toSeconds,
691                 maxCount));
692         logger.fine(
693             String.format(
694                 " Delta    %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f  %6d",
695                 (maxBSpline - minBSpline) * toSeconds,
696                 (maxConv - minConv) * toSeconds,
697                 (maxBSPerm - minBSPerm) * toSeconds,
698                 (maxPhiPerm - minPhiPerm) * toSeconds,
699                 (maxBSInduced - minBSInduced) * toSeconds,
700                 (maxPhiInduced - minPhiInduced) * toSeconds,
701                 (maxCount - minCount)));
702         logger.fine(
703             String.format(
704                 " Actual   %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f  %6d",
705                 bSplineTotal * toSeconds,
706                 convTotal * toSeconds,
707                 splinePermanentTotal * toSeconds,
708                 permanentPhiTotal * toSeconds,
709                 splineInducedTotal * toSeconds,
710                 inducedPhiTotal * toSeconds,
711                 nAtoms * nSymm));
712       }
713     }
714   }
715 
716   /**
717    * Setter for the field <code>atoms</code>.
718    *
719    * @param atoms an array of {@link ffx.potential.bonded.Atom} objects.
720    */
721   public void setAtoms(Atom[] atoms) {
722     this.atoms = atoms;
723     nAtoms = atoms.length;
724     coordinates = particleMeshEwald.getCoordinates();
725     switch (gridMethod) {
726       case SPATIAL:
727         spatialDensityRegion.setAtoms(atoms);
728         spatialDensityRegion.coordinates = coordinates;
729         break;
730       case ROW:
731         rowRegion.setAtoms(atoms);
732         rowRegion.coordinates = coordinates;
733         break;
734       case SLICE:
735         sliceRegion.setAtoms(atoms);
736         sliceRegion.coordinates = coordinates;
737     }
738   }
739 
740   /**
741    * Setter for the field <code>crystal</code>.
742    *
743    * @param crystal a {@link ffx.crystal.Crystal} object.
744    */
745   public void setCrystal(Crystal crystal) {
746     // Check if the number of symmetry operators has changed.
747     this.crystal = crystal.getUnitCell();
748     if (nSymm != this.crystal.spaceGroup.getNumberOfSymOps()) {
749       logger.info(this.crystal.toString());
750       logger.info(crystal.toString());
751       logger.severe(
752           " The reciprocal space class does not currently allow changes in the number of symmetry operators.");
753     }
754     this.coordinates = particleMeshEwald.getCoordinates();
755     initConvolution();
756   }
757 
758   /**
759    * Place the induced dipoles onto the FFT grid for the atoms in use.
760    *
761    * @param inducedDipole Induced dipoles.
762    * @param inducedDipoleCR Chain rule term for induced dipole gradient.
763    * @param use The atoms in use.
764    */
765   public void splineInducedDipoles(
766       double[][][] inducedDipole, double[][][] inducedDipoleCR, boolean[] use) {
767     splineInducedTotal -= System.nanoTime();
768 
769     switch (gridMethod) {
770       case SPATIAL:
771         spatialDensityRegion.setDensityLoop(spatialInducedLoops);
772         for (int i = 0; i < threadCount; i++) {
773           spatialInducedLoops[i].setInducedDipoles(inducedDipole, inducedDipoleCR);
774           spatialInducedLoops[i].setUse(use);
775           spatialInducedLoops[i].setRegion(spatialDensityRegion);
776         }
777         try {
778           parallelTeam.execute(spatialDensityRegion);
779         } catch (Exception e) {
780           String message = " Fatal exception evaluating induced density.\n";
781           logger.log(Level.SEVERE, message, e);
782         }
783         break;
784       case ROW:
785         rowRegion.setDensityLoop(rowInducedLoops);
786         for (int i = 0; i < threadCount; i++) {
787           rowInducedLoops[i].setInducedDipoles(inducedDipole, inducedDipoleCR);
788           rowInducedLoops[i].setUse(use);
789         }
790         try {
791           parallelTeam.execute(rowRegion);
792         } catch (Exception e) {
793           String message = " Fatal exception evaluating induced density.";
794           logger.log(Level.SEVERE, message, e);
795         }
796         break;
797       default:
798         sliceRegion.setDensityLoop(sliceInducedLoops);
799         for (int i = 0; i < threadCount; i++) {
800           sliceInducedLoops[i].setInducedDipoles(inducedDipole, inducedDipoleCR);
801           sliceInducedLoops[i].setUse(use);
802         }
803         try {
804           parallelTeam.execute(sliceRegion);
805         } catch (Exception e) {
806           String message = " Fatal exception evaluating induced density.";
807           logger.log(Level.SEVERE, message, e);
808         }
809         break;
810     }
811     splineInducedTotal += System.nanoTime();
812   }
813 
814   /**
815    * Use b-Splines to place the permanent multipoles onto the FFT grid for the atoms in use.
816    *
817    * @param globalMultipoles an array of double.
818    * @param fracMultipoles an array of double.
819    * @param use an array of boolean.
820    */
821   public void splinePermanentMultipoles(double[][][] globalMultipoles, double[][][] fracMultipoles, boolean[] use) {
822     splinePermanentTotal -= System.nanoTime();
823 
824     switch (gridMethod) {
825       case SPATIAL:
826         spatialDensityRegion.setCrystal(crystal.getUnitCell(), fftX, fftY, fftZ);
827         spatialDensityRegion.assignAtomsToCells();
828         spatialDensityRegion.setDensityLoop(spatialPermanentLoops);
829         for (int i = 0; i < threadCount; i++) {
830           spatialPermanentLoops[i].setPermanent(globalMultipoles, fracMultipoles);
831           spatialPermanentLoops[i].setUse(use);
832           spatialPermanentLoops[i].setRegion(spatialDensityRegion);
833         }
834         try {
835           parallelTeam.execute(spatialDensityRegion);
836         } catch (Exception e) {
837           String message = " Fatal exception evaluating permanent multipole density.";
838           logger.log(Level.SEVERE, message, e);
839         }
840         break;
841       case ROW:
842         rowRegion.setCrystal(crystal.getUnitCell(), fftX, fftY, fftZ);
843         rowRegion.setDensityLoop(rowPermanentLoops);
844         for (int i = 0; i < threadCount; i++) {
845           rowPermanentLoops[i].setPermanent(globalMultipoles, fracMultipoles);
846           rowPermanentLoops[i].setUse(use);
847         }
848         try {
849           parallelTeam.execute(rowRegion);
850         } catch (Exception e) {
851           String message = " Fatal exception evaluating permanent multipole density.";
852           logger.log(Level.SEVERE, message, e);
853         }
854         break;
855       case SLICE:
856       default:
857         sliceRegion.setCrystal(crystal.getUnitCell(), fftX, fftY, fftZ);
858         sliceRegion.setDensityLoop(slicePermanentLoops);
859         for (int i = 0; i < threadCount; i++) {
860           slicePermanentLoops[i].setPermanent(globalMultipoles, fracMultipoles);
861           slicePermanentLoops[i].setUse(use);
862         }
863         try {
864           parallelTeam.execute(sliceRegion);
865         } catch (Exception e) {
866           String message = " Fatal exception evaluating permanent multipole density.";
867           logger.log(Level.SEVERE, message, e);
868         }
869         break;
870     }
871 
872     splinePermanentTotal += System.nanoTime();
873   }
874 
875   private double initConvolution() {
876 
877     // Store the current reciprocal space grid dimensions.
878     int fftXCurrent = fftX;
879     int fftYCurrent = fftY;
880     int fftZCurrent = fftZ;
881 
882     double density = forceField.getDouble("PME_MESH_DENSITY", 1.2);
883     int nX = forceField.getInteger("PME_GRID_X", -1);
884     if (nX < 2) {
885       nX = (int) floor(crystal.a * density) + 1;
886       if (nX % 2 != 0) {
887         nX += 1;
888       }
889       while (!Complex.preferredDimension(nX)) {
890         nX += 2;
891       }
892     }
893     int nY = forceField.getInteger("PME_GRID_Y", -1);
894     if (nY < 2) {
895       nY = (int) floor(crystal.b * density) + 1;
896       if (nY % 2 != 0) {
897         nY += 1;
898       }
899       while (!Complex.preferredDimension(nY)) {
900         nY += 2;
901       }
902     }
903     int nZ = forceField.getInteger("PME_GRID_Z", -1);
904     if (nZ < 2) {
905       nZ = (int) floor(crystal.c * density) + 1;
906       if (nZ % 2 != 0) {
907         nZ += 1;
908       }
909       while (!Complex.preferredDimension(nZ)) {
910         nZ += 2;
911       }
912     }
913 
914     int minGrid = forceField.getInteger("PME_GRID_MIN", 16);
915     nX = max(nX, minGrid);
916     nY = max(nY, minGrid);
917     nZ = max(nZ, minGrid);
918 
919     fftX = nX;
920     fftY = nY;
921     fftZ = nZ;
922 
923     // Populate the matrix that fractionalizes multipoles.
924     transformMultipoleMatrix();
925 
926     // Populate the matrix that convert fractional potential components into orthogonal Cartesian
927     // coordinates.
928     transformFieldMatrix();
929 
930     // Compute the Cartesian to fractional matrix.
931     for (int i = 0; i < 3; i++) {
932       a[0][i] = fftX * crystal.A[i][0];
933       a[1][i] = fftY * crystal.A[i][1];
934       a[2][i] = fftZ * crystal.A[i][2];
935     }
936 
937     fftSpace = fftX * fftY * fftZ * 2;
938     boolean dimChanged = fftX != fftXCurrent || fftY != fftYCurrent || fftZ != fftZCurrent;
939 
940     if (complex3DFFT == null || dimChanged) {
941       complex3DFFT = new Complex3DParallel(fftX, fftY, fftZ, fftTeam, recipSchedule);
942       if (splineGrid == null || splineGrid.length < fftSpace) {
943         splineGrid = new double[fftSpace];
944       }
945       splineBuffer = DoubleBuffer.wrap(splineGrid);
946     }
947     complex3DFFT.setRecip(generalizedInfluenceFunction());
948 
949     switch (gridMethod) {
950       case SPATIAL:
951         if (spatialDensityRegion == null || dimChanged) {
952           spatialDensityRegion =
953               new SpatialDensityRegion(
954                   fftX,
955                   fftY,
956                   fftZ,
957                   splineGrid,
958                   bSplineOrder,
959                   nSymm,
960                   10,
961                   threadCount,
962                   crystal,
963                   atoms,
964                   coordinates);
965         } else {
966           spatialDensityRegion.setCrystal(crystal, fftX, fftY, fftZ);
967           spatialDensityRegion.coordinates = coordinates;
968         }
969         break;
970       case ROW:
971         if (rowRegion == null || dimChanged) {
972           rowRegion = new RowRegion(fftX, fftY, fftZ, splineGrid, nSymm, threadCount, atoms,
973               coordinates);
974         } else {
975           rowRegion.setCrystal(crystal, fftX, fftY, fftZ);
976           rowRegion.coordinates = coordinates;
977         }
978         break;
979       case SLICE:
980       default:
981         if (sliceRegion == null || dimChanged) {
982           sliceRegion = new SliceRegion(fftX, fftY, fftZ, splineGrid, nSymm, threadCount, atoms,
983               coordinates);
984         } else {
985           sliceRegion.setCrystal(crystal, fftX, fftY, fftZ);
986           sliceRegion.coordinates = coordinates;
987         }
988     }
989 
990     return density;
991   }
992 
993   private int RowIndexZ(int i) {
994     return i / fftY;
995   }
996 
997   private int RowIndexY(int i) {
998     return i % fftY;
999   }
1000 
1001   private double[] generalizedInfluenceFunction() {
1002 
1003     double[] influenceFunction = new double[fftSpace / 2];
1004 
1005     double[] bsModX = new double[fftX];
1006     double[] bsModY = new double[fftY];
1007     double[] bsModZ = new double[fftZ];
1008     int maxfft = max(max(max(fftX, fftY), fftZ), bSplineOrder + 1);
1009     double[] bsArray = new double[maxfft];
1010     double[] c = new double[bSplineOrder];
1011 
1012     bSpline(0.0, bSplineOrder, c);
1013     arraycopy(c, 0, bsArray, 1, bSplineOrder);
1014 
1015     discreteFTMod(bsModX, bsArray, fftX, bSplineOrder);
1016     discreteFTMod(bsModY, bsArray, fftY, bSplineOrder);
1017     discreteFTMod(bsModZ, bsArray, fftZ, bSplineOrder);
1018 
1019     double r00 = crystal.A[0][0];
1020     double r01 = crystal.A[0][1];
1021     double r02 = crystal.A[0][2];
1022     double r10 = crystal.A[1][0];
1023     double r11 = crystal.A[1][1];
1024     double r12 = crystal.A[1][2];
1025     double r20 = crystal.A[2][0];
1026     double r21 = crystal.A[2][1];
1027     double r22 = crystal.A[2][2];
1028     int ntot = fftX * fftY * fftZ;
1029     double piTerm = (PI / aEwald) * (PI / aEwald);
1030     double volTerm = PI * crystal.volume;
1031     int nfXY = fftX * fftY;
1032     int nX_2 = (fftX + 1) / 2;
1033     int nY_2 = (fftY + 1) / 2;
1034     int nZ_2 = (fftZ + 1) / 2;
1035 
1036     for (int i = 0; i < ntot - 1; i++) {
1037       int kZ = (i + 1) / nfXY;
1038       int j = i - kZ * nfXY + 1;
1039       int kY = j / fftX;
1040       int kX = j - kY * fftX;
1041       int h = kX;
1042       int k = kY;
1043       int l = kZ;
1044       if (kX >= nX_2) {
1045         h -= fftX;
1046       }
1047       if (kY >= nY_2) {
1048         k -= fftY;
1049       }
1050       if (kZ >= nZ_2) {
1051         l -= fftZ;
1052       }
1053       double sX = r00 * h + r01 * k + r02 * l;
1054       double sY = r10 * h + r11 * k + r12 * l;
1055       double sZ = r20 * h + r21 * k + r22 * l;
1056       double sSquared = sX * sX + sY * sY + sZ * sZ;
1057       double term = -piTerm * sSquared;
1058       double expterm = 0.0;
1059       if (term > -50.0) {
1060         double denom = sSquared * volTerm * bsModX[kX] * bsModY[kY] * bsModZ[kZ];
1061         expterm = exp(term) / denom;
1062         if (crystal.aperiodic()) {
1063           expterm *= (1.0 - cos(PI * crystal.a * sqrt(sSquared)));
1064         }
1065       }
1066       int ii = iComplex3D(kX, kY, kZ, fftX, fftY) / 2;
1067       influenceFunction[ii] = expterm;
1068     }
1069 
1070     // Account for the zeroth grid point for a periodic system.
1071     influenceFunction[0] = 0.0;
1072     if (crystal.aperiodic()) {
1073       influenceFunction[0] = 0.5 * PI / crystal.a;
1074     }
1075 
1076     return influenceFunction;
1077   }
1078 
1079   private void transformMultipoleMatrix() {
1080     double[][] a = new double[3][3];
1081     for (int i = 0; i < 3; i++) {
1082       a[0][i] = fftX * crystal.A[i][0];
1083       a[1][i] = fftY * crystal.A[i][1];
1084       a[2][i] = fftZ * crystal.A[i][2];
1085     }
1086 
1087     for (int i = 0; i < 10; i++) {
1088       for (int j = 0; j < 10; j++) {
1089         transformMultipoleMatrix[i][j] = 0.0;
1090       }
1091     }
1092 
1093     // Charge
1094     transformMultipoleMatrix[0][0] = 1.0;
1095     // Dipole
1096     for (int i = 1; i < 4; i++) {
1097       arraycopy(a[i - 1], 0, transformMultipoleMatrix[i], 1, 3);
1098     }
1099     // Quadrupole
1100     for (int i1 = 0; i1 < 3; i1++) {
1101       int k = qi1[i1];
1102       for (int i2 = 0; i2 < 6; i2++) {
1103         int i = qi1[i2];
1104         int j = qi2[i2];
1105         transformMultipoleMatrix[i1 + 4][i2 + 4] = a[k][i] * a[k][j];
1106       }
1107     }
1108     for (int i1 = 3; i1 < 6; i1++) {
1109       int k = qi1[i1];
1110       int l = qi2[i1];
1111       for (int i2 = 0; i2 < 6; i2++) {
1112         int i = qi1[i2];
1113         int j = qi2[i2];
1114         transformMultipoleMatrix[i1 + 4][i2 + 4] = a[k][i] * a[l][j] + a[k][j] * a[l][i];
1115       }
1116     }
1117   }
1118 
1119   private void transformFieldMatrix() {
1120     double[][] a = new double[3][3];
1121 
1122     for (int i = 0; i < 3; i++) {
1123       a[i][0] = fftX * crystal.A[i][0];
1124       a[i][1] = fftY * crystal.A[i][1];
1125       a[i][2] = fftZ * crystal.A[i][2];
1126     }
1127     for (int i = 0; i < 10; i++) {
1128       for (int j = 0; j < 10; j++) {
1129         transformFieldMatrix[i][j] = 0.0;
1130       }
1131     }
1132     transformFieldMatrix[0][0] = 1.0;
1133     for (int i = 1; i < 4; i++) {
1134       arraycopy(a[i - 1], 0, transformFieldMatrix[i], 1, 3);
1135     }
1136     for (int i1 = 0; i1 < 3; i1++) {
1137       int k = qi1[i1];
1138       for (int i2 = 0; i2 < 3; i2++) {
1139         int i = qi1[i2];
1140         transformFieldMatrix[i1 + 4][i2 + 4] = a[k][i] * a[k][i];
1141       }
1142       for (int i2 = 3; i2 < 6; i2++) {
1143         int i = qi1[i2];
1144         int j = qi2[i2];
1145         transformFieldMatrix[i1 + 4][i2 + 4] = 2.0 * a[k][i] * a[k][j];
1146       }
1147     }
1148     for (int i1 = 3; i1 < 6; i1++) {
1149       int k = qi1[i1];
1150       int n = qi2[i1];
1151       for (int i2 = 0; i2 < 3; i2++) {
1152         int i = qi1[i2];
1153         transformFieldMatrix[i1 + 4][i2 + 4] = a[k][i] * a[n][i];
1154       }
1155       for (int i2 = 3; i2 < 6; i2++) {
1156         int i = qi1[i2];
1157         int j = qi2[i2];
1158         transformFieldMatrix[i1 + 4][i2 + 4] = a[k][i] * a[n][j] + a[n][i] * a[k][j];
1159       }
1160     }
1161   }
1162 
1163   public enum GridMethod {
1164     SPATIAL,
1165     SLICE,
1166     ROW
1167   }
1168 
1169   /**
1170    * The class computes b-Splines that are used to spline multipoles and induced dipoles onto the PME
1171    * grid. Following convolution, the b-Splines are then used to obtain the reciprocal space
1172    * potential and fields from the PME grid. This class automatically updates itself to be consistent
1173    * with the current Crystal boundary conditions.
1174    *
1175    * <p>An external ParallelRegion can be used as follows: <code>
1176    * start() { bSplineRegion.start(); } run() { execute(0, nAtoms - 1,
1177    * bSplineRegion.bSplineLoop[threadID]); }
1178    * </code>
1179    */
1180   public class BSplineRegion extends ParallelRegion {
1181 
1182     final BSplineLoop[] bSplineLoop;
1183     double[][][][] splineX;
1184     double[][][][] splineY;
1185     double[][][][] splineZ;
1186     int[][][] initGrid;
1187     private double r00;
1188     private double r01;
1189     private double r02;
1190     private double r10;
1191     private double r11;
1192     private double r12;
1193     private double r20;
1194     private double r21;
1195     private double r22;
1196 
1197     BSplineRegion() {
1198       bSplineLoop = new BSplineLoop[threadCount];
1199       for (int i = 0; i < threadCount; i++) {
1200         bSplineLoop[i] = new BSplineLoop();
1201       }
1202     }
1203 
1204     @Override
1205     public void run() {
1206       /*
1207        Currently this condition would indicate a programming bug, since
1208        the space group is not allowed to change and the ReciprocalSpace
1209        class should be operating on a unit cell with a fixed number of
1210        space group symmetry operators (and not a replicated unit cell).
1211       */
1212       if (splineX.length < nSymm) {
1213         logger.severe(
1214             " Programming Error: the number of reciprocal space symmetry operators changed.");
1215       }
1216       try {
1217         int threadID = getThreadIndex();
1218         execute(0, nAtoms - 1, bSplineLoop[threadID]);
1219       } catch (Exception e) {
1220         e.printStackTrace();
1221         logger.severe(e.toString());
1222       }
1223     }
1224 
1225     @Override
1226     public void start() {
1227       r00 = crystal.A[0][0];
1228       r01 = crystal.A[0][1];
1229       r02 = crystal.A[0][2];
1230       r10 = crystal.A[1][0];
1231       r11 = crystal.A[1][1];
1232       r12 = crystal.A[1][2];
1233       r20 = crystal.A[2][0];
1234       r21 = crystal.A[2][1];
1235       r22 = crystal.A[2][2];
1236       if (splineX == null || splineX[0].length < nAtoms) {
1237         initGrid = new int[nSymm][nAtoms][];
1238         splineX = new double[nSymm][nAtoms][][];
1239         splineY = new double[nSymm][nAtoms][][];
1240         splineZ = new double[nSymm][nAtoms][][];
1241       }
1242     }
1243 
1244     public class BSplineLoop extends IntegerForLoop {
1245 
1246       private final double[][] bSplineWork;
1247       private final IntegerSchedule schedule = IntegerSchedule.fixed();
1248       // Extra padding to avert cache interference.
1249       private long pad0, pad1, pad2, pad3, pad4, pad5, pad6, pad7;
1250       private long pad8, pad9, pada, padb, padc, padd, pade, padf;
1251 
1252       BSplineLoop() {
1253         bSplineWork = new double[bSplineOrder][bSplineOrder];
1254       }
1255 
1256       @Override
1257       public void finish() {
1258         bSplineTime[getThreadIndex()] += System.nanoTime();
1259       }
1260 
1261       @Override
1262       public void run(int lb, int ub) {
1263         for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
1264           final double[] x = coordinates[iSymOp][0];
1265           final double[] y = coordinates[iSymOp][1];
1266           final double[] z = coordinates[iSymOp][2];
1267           final int[][] initgridi = initGrid[iSymOp];
1268           final double[][][] splineXi = splineX[iSymOp];
1269           final double[][][] splineYi = splineY[iSymOp];
1270           final double[][][] splineZi = splineZ[iSymOp];
1271           for (int i = lb; i <= ub; i++) {
1272             if (initgridi[i] == null) {
1273               initgridi[i] = new int[3];
1274               splineXi[i] = new double[bSplineOrder][derivOrder + 1];
1275               splineYi[i] = new double[bSplineOrder][derivOrder + 1];
1276               splineZi[i] = new double[bSplineOrder][derivOrder + 1];
1277             }
1278             final double xi = x[i];
1279             final double yi = y[i];
1280             final double zi = z[i];
1281             final int[] grd = initgridi[i];
1282             // X-dimension
1283             final double wx = xi * r00 + yi * r10 + zi * r20;
1284             final double ux = wx - round(wx) + 0.5;
1285             final double frx = fftX * ux;
1286             final int ifrx = (int) frx;
1287             final double bx = frx - ifrx;
1288             grd[0] = ifrx - bSplineOrder;
1289             bSplineDerivatives(bx, bSplineOrder, derivOrder, splineXi[i], bSplineWork);
1290             // Y-dimension
1291             final double wy = xi * r01 + yi * r11 + zi * r21;
1292             final double uy = wy - round(wy) + 0.5;
1293             final double fry = fftY * uy;
1294             final int ifry = (int) fry;
1295             final double by = fry - ifry;
1296             grd[1] = ifry - bSplineOrder;
1297             bSplineDerivatives(by, bSplineOrder, derivOrder, splineYi[i], bSplineWork);
1298             // Z-dimension
1299             final double wz = xi * r02 + yi * r12 + zi * r22;
1300             final double uz = wz - round(wz) + 0.5;
1301             final double frz = fftZ * uz;
1302             final int ifrz = (int) frz;
1303             final double bz = frz - ifrz;
1304             grd[2] = ifrz - bSplineOrder;
1305             bSplineDerivatives(bz, bSplineOrder, derivOrder, splineZi[i], bSplineWork);
1306           }
1307         }
1308       }
1309 
1310       @Override
1311       public IntegerSchedule schedule() {
1312         return schedule;
1313       }
1314 
1315       @Override
1316       public void start() {
1317         bSplineTime[getThreadIndex()] -= System.nanoTime();
1318       }
1319     }
1320   }
1321 
1322   private class SpatialPermanentLoop extends SpatialDensityLoop {
1323 
1324     private final BSplineRegion bSplines;
1325     private double[][][] globalMultipoles = null;
1326     private double[][][] fracMultipoles = null;
1327     private boolean[] use = null;
1328     private int threadIndex;
1329 
1330     SpatialPermanentLoop(SpatialDensityRegion region, BSplineRegion splines) {
1331       super(region, region.nSymm, region.actualCount);
1332       this.bSplines = splines;
1333     }
1334 
1335     @Override
1336     public void finish() {
1337       splinePermanentTime[threadIndex] += System.nanoTime();
1338     }
1339 
1340     @Override
1341     public void gridDensity(int iSymm, int n) {
1342       splineCount[threadIndex]++;
1343 
1344       // Convert Cartesian multipoles in the global frame to fractional multipoles.
1345       final double[] gm = globalMultipoles[iSymm][n];
1346       final double[] fm = fracMultipoles[iSymm][n];
1347 
1348       // Charge.
1349       fm[0] = gm[0];
1350 
1351       // Dipole.
1352       for (int j = 1; j < 4; j++) {
1353         fm[j] = 0.0;
1354         for (int k = 1; k < 4; k++) {
1355           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
1356         }
1357       }
1358 
1359       // Quadrupole.
1360       for (int j = 4; j < 10; j++) {
1361         fm[j] = 0.0;
1362         for (int k = 4; k < 7; k++) {
1363           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
1364         }
1365         for (int k = 7; k < 10; k++) {
1366           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * 2.0 * gm[k];
1367         }
1368         /*
1369          Fractional quadrupole components are pre-multiplied by a
1370          factor of 1/3 that arises in their potential.
1371         */
1372         fm[j] = fm[j] / 3.0;
1373       }
1374 
1375       // Some atoms are not used during Lambda dynamics.
1376       if (use != null && !use[n]) {
1377         return;
1378       }
1379 
1380       final double[][] splx = bSplines.splineX[iSymm][n];
1381       final double[][] sply = bSplines.splineY[iSymm][n];
1382       final double[][] splz = bSplines.splineZ[iSymm][n];
1383       final int igrd0 = bSplines.initGrid[iSymm][n][0];
1384       final int jgrd0 = bSplines.initGrid[iSymm][n][1];
1385       int k0 = bSplines.initGrid[iSymm][n][2];
1386       final double c = fm[t000];
1387       final double dx = fm[t100];
1388       final double dy = fm[t010];
1389       final double dz = fm[t001];
1390       final double qxx = fm[t200];
1391       final double qyy = fm[t020];
1392       final double qzz = fm[t002];
1393       final double qxy = fm[t110];
1394       final double qxz = fm[t101];
1395       final double qyz = fm[t011];
1396       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1397         final double[] splzi = splz[ith3];
1398         final double v0 = splzi[0];
1399         final double v1 = splzi[1];
1400         final double v2 = splzi[2];
1401         final double c0 = c * v0;
1402         final double dx0 = dx * v0;
1403         final double dy0 = dy * v0;
1404         final double dz1 = dz * v1;
1405         final double qxx0 = qxx * v0;
1406         final double qyy0 = qyy * v0;
1407         final double qzz2 = qzz * v2;
1408         final double qxy0 = qxy * v0;
1409         final double qxz1 = qxz * v1;
1410         final double qyz1 = qyz * v1;
1411         final int k = mod(++k0, fftZ);
1412         int j0 = jgrd0;
1413         for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
1414           final double[] splyi = sply[ith2];
1415           final double u0 = splyi[0];
1416           final double u1 = splyi[1];
1417           final double u2 = splyi[2];
1418           final double term0 = (c0 + dz1 + qzz2) * u0 + (dy0 + qyz1) * u1 + qyy0 * u2;
1419           final double term1 = (dx0 + qxz1) * u0 + qxy0 * u1;
1420           final double term2 = qxx0 * u0;
1421           final int j = mod(++j0, fftY);
1422           int i0 = igrd0;
1423           for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
1424             final int i = mod(++i0, fftX);
1425             final int ii = iComplex3D(i, j, k, fftX, fftY);
1426             final double[] splxi = splx[ith1];
1427             final double add = splxi[0] * term0 + splxi[1] * term1 + splxi[2] * term2;
1428             final double current = splineBuffer.get(ii);
1429             splineBuffer.put(ii, current + add);
1430           }
1431         }
1432       }
1433     }
1434 
1435     @Override
1436     public void start() {
1437       threadIndex = getThreadIndex();
1438       splinePermanentTime[threadIndex] -= System.nanoTime();
1439     }
1440 
1441     void setPermanent(double[][][] globalMultipoles, double[][][] fracMultipoles) {
1442       this.globalMultipoles = globalMultipoles;
1443       this.fracMultipoles = fracMultipoles;
1444     }
1445 
1446     private void setUse(boolean[] use) {
1447       this.use = use;
1448     }
1449   }
1450 
1451   private class SpatialInducedLoop extends SpatialDensityLoop {
1452 
1453     private final double[][] m = new double[3][3];
1454     private final BSplineRegion bSplineRegion;
1455     private double[][][] inducedDipole = null;
1456     private double[][][] inducedDipoleCR = null;
1457     private boolean[] use = null;
1458     private double a00, a01, a02;
1459     private double a10, a11, a12;
1460     private double a20, a21, a22;
1461 
1462     SpatialInducedLoop(SpatialDensityRegion region, BSplineRegion splines) {
1463       super(region, region.nSymm, region.actualCount);
1464       this.bSplineRegion = splines;
1465     }
1466 
1467     @Override
1468     public void finish() {
1469       splineInducedTime[getThreadIndex()] += System.nanoTime();
1470     }
1471 
1472     @Override
1473     public void gridDensity(int iSymm, int n) {
1474 
1475       // Convert Cartesian induced dipole to fractional induced dipole.
1476       double[] ind = inducedDipole[iSymm][n];
1477       double dx = ind[0];
1478       double dy = ind[1];
1479       double dz = ind[2];
1480       final double ux = a00 * dx + a01 * dy + a02 * dz;
1481       final double uy = a10 * dx + a11 * dy + a12 * dz;
1482       final double uz = a20 * dx + a21 * dy + a22 * dz;
1483       // final double[] find = fracInducedDipole[iSymm][n];
1484       // find[0] = ux;
1485       // find[1] = uy;
1486       // find[2] = uz;
1487 
1488       // Convert Cartesian induced dipole CR term to fractional induced dipole CR term.
1489       double[] indCR = inducedDipoleCR[iSymm][n];
1490       dx = indCR[0];
1491       dy = indCR[1];
1492       dz = indCR[2];
1493       final double px = a00 * dx + a01 * dy + a02 * dz;
1494       final double py = a10 * dx + a11 * dy + a12 * dz;
1495       final double pz = a20 * dx + a21 * dy + a22 * dz;
1496       // final double[] findCR = fracInducedDipoleCR[iSymm][n];
1497       // findCR[0] = px;
1498       // findCR[1] = py;
1499       // findCR[2] = pz;
1500       if (use != null && !use[n]) {
1501         return;
1502       }
1503       final double[][] splx = bSplineRegion.splineX[iSymm][n];
1504       final double[][] sply = bSplineRegion.splineY[iSymm][n];
1505       final double[][] splz = bSplineRegion.splineZ[iSymm][n];
1506       final int igrd0 = bSplineRegion.initGrid[iSymm][n][0];
1507       final int jgrd0 = bSplineRegion.initGrid[iSymm][n][1];
1508       int k0 = bSplineRegion.initGrid[iSymm][n][2];
1509       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1510         final double[] splzi = splz[ith3];
1511         final double v0 = splzi[0];
1512         final double v1 = splzi[1];
1513         final double dx0 = ux * v0;
1514         final double dy0 = uy * v0;
1515         final double dz1 = uz * v1;
1516         final double px0 = px * v0;
1517         final double py0 = py * v0;
1518         final double pz1 = pz * v1;
1519         final int k = mod(++k0, fftZ);
1520         int j0 = jgrd0;
1521         for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
1522           final double[] splyi = sply[ith2];
1523           final double u0 = splyi[0];
1524           final double u1 = splyi[1];
1525           final double term0 = dz1 * u0 + dy0 * u1;
1526           final double term1 = dx0 * u0;
1527           final double termp0 = pz1 * u0 + py0 * u1;
1528           final double termp1 = px0 * u0;
1529           final int j = mod(++j0, fftY);
1530           int i0 = igrd0;
1531           for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
1532             final int i = mod(++i0, fftX);
1533             final int ii = iComplex3D(i, j, k, fftX, fftY);
1534             final double[] splxi = splx[ith1];
1535             final double add = splxi[0] * term0 + splxi[1] * term1;
1536             final double addi = splxi[0] * termp0 + splxi[1] * termp1;
1537             final double current = splineBuffer.get(ii);
1538             final double currenti = splineBuffer.get(ii + 1);
1539             splineBuffer.put(ii, current + add);
1540             splineBuffer.put(ii + 1, currenti + addi);
1541             // splineGrid[ii] += add;
1542             // splineGrid[ii + 1] += addi;
1543           }
1544         }
1545       }
1546     }
1547 
1548     public void setUse(boolean[] use) {
1549       this.use = use;
1550     }
1551 
1552     @Override
1553     public void start() {
1554       splineInducedTime[getThreadIndex()] -= System.nanoTime();
1555       for (int i = 0; i < 3; i++) {
1556         m[0][i] = fftX * crystal.A[i][0];
1557         m[1][i] = fftY * crystal.A[i][1];
1558         m[2][i] = fftZ * crystal.A[i][2];
1559       }
1560       a00 = m[0][0];
1561       a01 = m[0][1];
1562       a02 = m[0][2];
1563       a10 = m[1][0];
1564       a11 = m[1][1];
1565       a12 = m[1][2];
1566       a20 = m[2][0];
1567       a21 = m[2][1];
1568       a22 = m[2][2];
1569     }
1570 
1571     void setInducedDipoles(double[][][] inducedDipole, double[][][] inducedDipoleCR) {
1572       this.inducedDipole = inducedDipole;
1573       this.inducedDipoleCR = inducedDipoleCR;
1574     }
1575   }
1576 
1577   private class RowPermanentLoop extends RowLoop {
1578 
1579     private final double[] fracMPole = new double[10];
1580     private final BSplineRegion bSplines;
1581     private double[][][] globalMultipoles = null;
1582     private double[][][] fracMultipoles = null;
1583     private boolean[] use = null;
1584     private int threadIndex;
1585 
1586     RowPermanentLoop(RowRegion region, BSplineRegion splines) {
1587       super(region.nAtoms, region.nSymm, region);
1588       this.bSplines = splines;
1589     }
1590 
1591     @Override
1592     public void finish() {
1593       splinePermanentTime[threadIndex] += System.nanoTime();
1594     }
1595 
1596     @Override
1597     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
1598       boolean atomContributes = false;
1599       int k0 = bSplines.initGrid[iSymm][iAtom][2];
1600       int lbZ = RowIndexZ(lb);
1601       // int lbY = RowIndexY(lb);
1602       int ubZ = RowIndexZ(ub);
1603       // int ubY = RowIndexY(ub);
1604       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1605         final int k = mod(++k0, fftZ);
1606         if (lbZ <= k && k <= ubZ) {
1607           atomContributes = true;
1608           break;
1609         }
1610       }
1611 
1612       if (!atomContributes) {
1613         return;
1614       }
1615 
1616       splineCount[threadIndex]++;
1617 
1618       // Add atom n to the list for the current symmOp and thread.
1619       int index = gridAtomCount[iSymm][threadIndex];
1620       gridAtomList[iSymm][threadIndex][index] = iAtom;
1621       gridAtomCount[iSymm][threadIndex]++;
1622 
1623       // Convert Cartesian multipoles in the global frame to fractional multipoles.
1624       final double[] gm = globalMultipoles[iSymm][iAtom];
1625       final double[] fm = fracMPole;
1626 
1627       // Charge.
1628       fm[0] = gm[0];
1629 
1630       // Dipole.
1631       for (int j = 1; j < 4; j++) {
1632         fm[j] = 0.0;
1633         for (int k = 1; k < 4; k++) {
1634           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
1635         }
1636       }
1637 
1638       // Quadrupole.
1639       for (int j = 4; j < 10; j++) {
1640         fm[j] = 0.0;
1641         for (int k = 4; k < 7; k++) {
1642           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
1643         }
1644         for (int k = 7; k < 10; k++) {
1645           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * 2.0 * gm[k];
1646         }
1647         // Fractional quadrupole components are pre-multiplied by a
1648         // factor of 1/3 that arises in their potential.
1649         fm[j] = fm[j] / 3.0;
1650       }
1651 
1652       arraycopy(fm, 0, fracMultipoles[iSymm][iAtom], 0, 10);
1653 
1654       // Some atoms are not used during Lambda dynamics.
1655       if (use != null && !use[iAtom]) {
1656         return;
1657       }
1658 
1659       final double[][] splx = bSplines.splineX[iSymm][iAtom];
1660       final double[][] sply = bSplines.splineY[iSymm][iAtom];
1661       final double[][] splz = bSplines.splineZ[iSymm][iAtom];
1662       final int igrd0 = bSplines.initGrid[iSymm][iAtom][0];
1663       final int jgrd0 = bSplines.initGrid[iSymm][iAtom][1];
1664       k0 = bSplines.initGrid[iSymm][iAtom][2];
1665       final double c = fm[t000];
1666       final double dx = fm[t100];
1667       final double dy = fm[t010];
1668       final double dz = fm[t001];
1669       final double qxx = fm[t200];
1670       final double qyy = fm[t020];
1671       final double qzz = fm[t002];
1672       final double qxy = fm[t110];
1673       final double qxz = fm[t101];
1674       final double qyz = fm[t011];
1675       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1676         final double[] splzi = splz[ith3];
1677         final double v0 = splzi[0];
1678         final double v1 = splzi[1];
1679         final double v2 = splzi[2];
1680         final double c0 = c * v0;
1681         final double dx0 = dx * v0;
1682         final double dy0 = dy * v0;
1683         final double dz1 = dz * v1;
1684         final double qxx0 = qxx * v0;
1685         final double qyy0 = qyy * v0;
1686         final double qzz2 = qzz * v2;
1687         final double qxy0 = qxy * v0;
1688         final double qxz1 = qxz * v1;
1689         final double qyz1 = qyz * v1;
1690         final int k = mod(++k0, fftZ);
1691         if (k < lbZ || k > ubZ) {
1692           continue;
1693         }
1694         int j0 = jgrd0;
1695         for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
1696           final double[] splyi = sply[ith2];
1697           final double u0 = splyi[0];
1698           final double u1 = splyi[1];
1699           final double u2 = splyi[2];
1700           // Pieces of a multipole
1701           final double term0 = (c0 + dz1 + qzz2) * u0 + (dy0 + qyz1) * u1 + qyy0 * u2;
1702           final double term1 = (dx0 + qxz1) * u0 + qxy0 * u1;
1703           final double term2 = qxx0 * u0;
1704           final int j = mod(++j0, fftY);
1705           int rowIndex = rowRegion.rowIndexForYZ(j, k);
1706           if (lb > rowIndex || rowIndex > ub) {
1707             continue;
1708           }
1709           int i0 = igrd0;
1710           for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
1711             final int i = mod(++i0, fftX);
1712             final int ii = iComplex3D(i, j, k, fftX, fftY);
1713             final double[] splxi = splx[ith1];
1714             final double add = splxi[0] * term0 + splxi[1] * term1 + splxi[2] * term2;
1715             final double current = splineBuffer.get(ii);
1716             splineBuffer.put(ii, current + add);
1717           }
1718         }
1719       }
1720     }
1721 
1722     @Override
1723     public void start() {
1724       threadIndex = getThreadIndex();
1725       splinePermanentTime[threadIndex] -= System.nanoTime();
1726       for (int i = 0; i < nSymm; i++) {
1727         gridAtomCount[i][threadIndex] = 0;
1728       }
1729     }
1730 
1731     void setPermanent(double[][][] globalMultipoles, double[][][] fracMultipoles) {
1732       this.globalMultipoles = globalMultipoles;
1733       this.fracMultipoles = fracMultipoles;
1734     }
1735 
1736     private void setUse(boolean[] use) {
1737       this.use = use;
1738     }
1739   }
1740 
1741   private class RowInducedLoop extends RowLoop {
1742 
1743     private final double[][] m = new double[3][3];
1744     private final BSplineRegion bSplineRegion;
1745     private double[][][] inducedDipole = null;
1746     private double[][][] inducedDipoleCR = null;
1747     private boolean[] use = null;
1748     private double a00, a01, a02;
1749     private double a10, a11, a12;
1750     private double a20, a21, a22;
1751 
1752     RowInducedLoop(RowRegion region, BSplineRegion splines) {
1753       super(region.nAtoms, region.nSymm, region);
1754       this.bSplineRegion = splines;
1755     }
1756 
1757     @Override
1758     public void finish() {
1759       int threadIndex = getThreadIndex();
1760       splineInducedTime[threadIndex] += System.nanoTime();
1761     }
1762 
1763     @Override
1764     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
1765 
1766       // Convert Cartesian induced dipole to fractional induced dipole.
1767       int lbZ = RowIndexZ(lb);
1768       int ubZ = RowIndexZ(ub);
1769       double[] ind = inducedDipole[iSymm][iAtom];
1770       double dx = ind[0];
1771       double dy = ind[1];
1772       double dz = ind[2];
1773       final double ux = a00 * dx + a01 * dy + a02 * dz;
1774       final double uy = a10 * dx + a11 * dy + a12 * dz;
1775       final double uz = a20 * dx + a21 * dy + a22 * dz;
1776       // final double[] find = fracInducedDipole[iSymm][iAtom];
1777       // find[0] = ux;
1778       // find[1] = uy;
1779       // find[2] = uz;
1780 
1781       // Convert Cartesian induced dipole CR term to fractional induced dipole CR term.
1782       double[] indCR = inducedDipoleCR[iSymm][iAtom];
1783       dx = indCR[0];
1784       dy = indCR[1];
1785       dz = indCR[2];
1786       final double px = a00 * dx + a01 * dy + a02 * dz;
1787       final double py = a10 * dx + a11 * dy + a12 * dz;
1788       final double pz = a20 * dx + a21 * dy + a22 * dz;
1789       // final double[] findCR = fracInducedDipoleCR[iSymm][iAtom];
1790       // findCR[0] = px;
1791       // findCR[1] = py;
1792       // findCR[2] = pz;
1793       if (use != null && !use[iAtom]) {
1794         return;
1795       }
1796       final double[][] splx = bSplineRegion.splineX[iSymm][iAtom];
1797       final double[][] sply = bSplineRegion.splineY[iSymm][iAtom];
1798       final double[][] splz = bSplineRegion.splineZ[iSymm][iAtom];
1799       final int igrd0 = bSplineRegion.initGrid[iSymm][iAtom][0];
1800       final int jgrd0 = bSplineRegion.initGrid[iSymm][iAtom][1];
1801       int k0 = bSplineRegion.initGrid[iSymm][iAtom][2];
1802       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1803         final double[] splzi = splz[ith3];
1804         final double v0 = splzi[0];
1805         final double v1 = splzi[1];
1806         final double dx0 = ux * v0;
1807         final double dy0 = uy * v0;
1808         final double dz1 = uz * v1;
1809         final double px0 = px * v0;
1810         final double py0 = py * v0;
1811         final double pz1 = pz * v1;
1812         final int k = mod(++k0, fftZ);
1813         if (k < lbZ || k > ubZ) {
1814           continue;
1815         }
1816         int j0 = jgrd0;
1817         for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
1818           final double[] splyi = sply[ith2];
1819           final double u0 = splyi[0];
1820           final double u1 = splyi[1];
1821           final double term0 = dz1 * u0 + dy0 * u1;
1822           final double term1 = dx0 * u0;
1823           final double termp0 = pz1 * u0 + py0 * u1;
1824           final double termp1 = px0 * u0;
1825           final int j = mod(++j0, fftY);
1826           int rowIndex = rowRegion.rowIndexForYZ(j, k);
1827           if (lb > rowIndex || rowIndex > ub) {
1828             continue;
1829           }
1830           int i0 = igrd0;
1831           for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
1832             final int i = mod(++i0, fftX);
1833             final int ii = iComplex3D(i, j, k, fftX, fftY);
1834             final double[] splxi = splx[ith1];
1835             final double add = splxi[0] * term0 + splxi[1] * term1;
1836             final double addi = splxi[0] * termp0 + splxi[1] * termp1;
1837             final double current = splineBuffer.get(ii);
1838             final double currenti = splineBuffer.get(ii + 1);
1839             splineBuffer.put(ii, current + add);
1840             splineBuffer.put(ii + 1, currenti + addi);
1841             // splineGrid[ii] += add;
1842             // splineGrid[ii + 1] += addi;
1843           }
1844         }
1845       }
1846     }
1847 
1848     @Override
1849     public void run(int lb, int ub) throws Exception {
1850       int ti = getThreadIndex();
1851       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
1852         int[] list = gridAtomList[iSymm][ti];
1853         int n = gridAtomCount[iSymm][ti];
1854         for (int i = 0; i < n; i++) {
1855           int iAtom = list[i];
1856           gridDensity(iSymm, iAtom, lb, ub);
1857         }
1858       }
1859     }
1860 
1861     public void setUse(boolean[] use) {
1862       this.use = use;
1863     }
1864 
1865     @Override
1866     public void start() {
1867       int threadIndex = getThreadIndex();
1868       splineInducedTime[threadIndex] -= System.nanoTime();
1869       for (int i = 0; i < 3; i++) {
1870         m[0][i] = fftX * crystal.A[i][0];
1871         m[1][i] = fftY * crystal.A[i][1];
1872         m[2][i] = fftZ * crystal.A[i][2];
1873       }
1874       a00 = m[0][0];
1875       a01 = m[0][1];
1876       a02 = m[0][2];
1877       a10 = m[1][0];
1878       a11 = m[1][1];
1879       a12 = m[1][2];
1880       a20 = m[2][0];
1881       a21 = m[2][1];
1882       a22 = m[2][2];
1883     }
1884 
1885     void setInducedDipoles(double[][][] inducedDipole, double[][][] inducedDipoleCR) {
1886       this.inducedDipole = inducedDipole;
1887       this.inducedDipoleCR = inducedDipoleCR;
1888     }
1889   }
1890 
1891   private class SlicePermanentLoop extends SliceLoop {
1892 
1893     private final double[] fracMPole = new double[10];
1894     private final BSplineRegion bSplines;
1895     private double[][][] globalMultipoles = null;
1896     private double[][][] fracMultipoles = null;
1897     private boolean[] use = null;
1898     private int threadIndex;
1899 
1900     SlicePermanentLoop(SliceRegion region, BSplineRegion splines) {
1901       super(region.nAtoms, region.nSymm, region);
1902       this.bSplines = splines;
1903     }
1904 
1905     @Override
1906     public void finish() {
1907       splinePermanentTime[threadIndex] += System.nanoTime();
1908     }
1909 
1910     @Override
1911     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
1912       boolean atomContributes = false;
1913       int k0 = bSplines.initGrid[iSymm][iAtom][2];
1914       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1915         final int k = mod(++k0, fftZ);
1916         if (lb <= k && k <= ub) {
1917           atomContributes = true;
1918           break;
1919         }
1920       }
1921       if (!atomContributes) {
1922         return;
1923       }
1924 
1925       splineCount[threadIndex]++;
1926 
1927       // Add atom n to the list for the current symmOp and thread.
1928       int index = gridAtomCount[iSymm][threadIndex];
1929       gridAtomList[iSymm][threadIndex][index] = iAtom;
1930       gridAtomCount[iSymm][threadIndex]++;
1931 
1932       // Convert Cartesian multipoles in the global frame to fractional multipoles.
1933       final double[] gm = globalMultipoles[iSymm][iAtom];
1934       final double[] fm = fracMPole;
1935 
1936       // Charge
1937       fm[0] = gm[0];
1938       // Dipole
1939       for (int j = 1; j < 4; j++) {
1940         fm[j] = 0.0;
1941         for (int k = 1; k < 4; k++) {
1942           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
1943         }
1944       }
1945       // Quadrupole
1946       for (int j = 4; j < 10; j++) {
1947         fm[j] = 0.0;
1948         for (int k = 4; k < 7; k++) {
1949           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
1950         }
1951         for (int k = 7; k < 10; k++) {
1952           fm[j] = fm[j] + transformMultipoleMatrix[j][k] * 2.0 * gm[k];
1953         }
1954         /*
1955         Fractional quadrupole components are pre-multiplied by a
1956         factor of 1/3 that arises in their potential.
1957         */
1958         fm[j] = fm[j] / 3.0;
1959       }
1960 
1961       arraycopy(fm, 0, fracMultipoles[iSymm][iAtom], 0, 10);
1962 
1963       // Some atoms are not used during Lambda dynamics.
1964       if (use != null && !use[iAtom]) {
1965         return;
1966       }
1967 
1968       final double[][] splx = bSplines.splineX[iSymm][iAtom];
1969       final double[][] sply = bSplines.splineY[iSymm][iAtom];
1970       final double[][] splz = bSplines.splineZ[iSymm][iAtom];
1971       final int igrd0 = bSplines.initGrid[iSymm][iAtom][0];
1972       final int jgrd0 = bSplines.initGrid[iSymm][iAtom][1];
1973       k0 = bSplines.initGrid[iSymm][iAtom][2];
1974       final double c = fm[t000];
1975       final double dx = fm[t100];
1976       final double dy = fm[t010];
1977       final double dz = fm[t001];
1978       final double qxx = fm[t200];
1979       final double qyy = fm[t020];
1980       final double qzz = fm[t002];
1981       final double qxy = fm[t110];
1982       final double qxz = fm[t101];
1983       final double qyz = fm[t011];
1984       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
1985         final double[] splzi = splz[ith3];
1986         final double v0 = splzi[0];
1987         final double v1 = splzi[1];
1988         final double v2 = splzi[2];
1989         final double c0 = c * v0;
1990         final double dx0 = dx * v0;
1991         final double dy0 = dy * v0;
1992         final double dz1 = dz * v1;
1993         final double qxx0 = qxx * v0;
1994         final double qyy0 = qyy * v0;
1995         final double qzz2 = qzz * v2;
1996         final double qxy0 = qxy * v0;
1997         final double qxz1 = qxz * v1;
1998         final double qyz1 = qyz * v1;
1999         final int k = mod(++k0, fftZ);
2000         if (k < lb || k > ub) {
2001           continue;
2002         }
2003         int j0 = jgrd0;
2004         for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
2005           final double[] splyi = sply[ith2];
2006           final double u0 = splyi[0];
2007           final double u1 = splyi[1];
2008           final double u2 = splyi[2];
2009           // Pieces of a multipole
2010           final double term0 = (c0 + dz1 + qzz2) * u0 + (dy0 + qyz1) * u1 + qyy0 * u2;
2011           final double term1 = (dx0 + qxz1) * u0 + qxy0 * u1;
2012           final double term2 = qxx0 * u0;
2013           final int j = mod(++j0, fftY);
2014           int i0 = igrd0;
2015           for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
2016             final int i = mod(++i0, fftX);
2017             final int ii = iComplex3D(i, j, k, fftX, fftY);
2018             final double[] splxi = splx[ith1];
2019             final double add = splxi[0] * term0 + splxi[1] * term1 + splxi[2] * term2;
2020             final double current = splineBuffer.get(ii);
2021             splineBuffer.put(ii, current + add);
2022             /*
2023             if (n == 0) {
2024             logger.info(String.format(" %d %16.8f", ii, current + add));
2025             } */
2026             // splineGrid[ii] += add;
2027           }
2028         }
2029       }
2030     }
2031 
2032     @Override
2033     public void start() {
2034       threadIndex = getThreadIndex();
2035       splinePermanentTime[threadIndex] -= System.nanoTime();
2036       for (int i = 0; i < nSymm; i++) {
2037         gridAtomCount[i][threadIndex] = 0;
2038       }
2039     }
2040 
2041     void setPermanent(double[][][] globalMultipoles, double[][][] fracMultipoles) {
2042       this.globalMultipoles = globalMultipoles;
2043       this.fracMultipoles = fracMultipoles;
2044     }
2045 
2046     private void setUse(boolean[] use) {
2047       this.use = use;
2048     }
2049   }
2050 
2051   private class SliceInducedLoop extends SliceLoop {
2052 
2053     private final double[][] m = new double[3][3];
2054     private final BSplineRegion bSplineRegion;
2055     private double[][][] inducedDipole = null;
2056     private double[][][] inducedDipoleCR = null;
2057     private boolean[] use = null;
2058     private double a00, a01, a02;
2059     private double a10, a11, a12;
2060     private double a20, a21, a22;
2061 
2062     SliceInducedLoop(SliceRegion region, BSplineRegion splines) {
2063       super(region.nAtoms, region.nSymm, region);
2064       this.bSplineRegion = splines;
2065     }
2066 
2067     @Override
2068     public void finish() {
2069       int threadIndex = getThreadIndex();
2070       splineInducedTime[threadIndex] += System.nanoTime();
2071     }
2072 
2073     @Override
2074     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2075       // Convert Cartesian induced dipole to fractional induced dipole.
2076       double[] ind = inducedDipole[iSymm][iAtom];
2077       double dx = ind[0];
2078       double dy = ind[1];
2079       double dz = ind[2];
2080       final double ux = a00 * dx + a01 * dy + a02 * dz;
2081       final double uy = a10 * dx + a11 * dy + a12 * dz;
2082       final double uz = a20 * dx + a21 * dy + a22 * dz;
2083       // final double[] find = fracInducedDipole[iSymm][iAtom];
2084       // find[0] = ux;
2085       // find[1] = uy;
2086       // find[2] = uz;
2087 
2088       // Convert Cartesian induced dipole CR term to fractional induced dipole CR term.
2089       double[] indCR = inducedDipoleCR[iSymm][iAtom];
2090       dx = indCR[0];
2091       dy = indCR[1];
2092       dz = indCR[2];
2093       final double px = a00 * dx + a01 * dy + a02 * dz;
2094       final double py = a10 * dx + a11 * dy + a12 * dz;
2095       final double pz = a20 * dx + a21 * dy + a22 * dz;
2096       // final double[] findCR = fracInducedDipoleCR[iSymm][iAtom];
2097       // findCR[0] = px;
2098       // findCR[1] = py;
2099       // findCR[2] = pz;
2100 
2101       if (use != null && !use[iAtom]) {
2102         return;
2103       }
2104       final double[][] splx = bSplineRegion.splineX[iSymm][iAtom];
2105       final double[][] sply = bSplineRegion.splineY[iSymm][iAtom];
2106       final double[][] splz = bSplineRegion.splineZ[iSymm][iAtom];
2107       final int igrd0 = bSplineRegion.initGrid[iSymm][iAtom][0];
2108       final int jgrd0 = bSplineRegion.initGrid[iSymm][iAtom][1];
2109       int k0 = bSplineRegion.initGrid[iSymm][iAtom][2];
2110       for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
2111         final double[] splzi = splz[ith3];
2112         final double v0 = splzi[0];
2113         final double v1 = splzi[1];
2114         final double dx0 = ux * v0;
2115         final double dy0 = uy * v0;
2116         final double dz1 = uz * v1;
2117         final double px0 = px * v0;
2118         final double py0 = py * v0;
2119         final double pz1 = pz * v1;
2120         final int k = mod(++k0, fftZ);
2121         if (k < lb || k > ub) {
2122           continue;
2123         }
2124         int j0 = jgrd0;
2125         for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
2126           final double[] splyi = sply[ith2];
2127           final double u0 = splyi[0];
2128           final double u1 = splyi[1];
2129           final double term0 = dz1 * u0 + dy0 * u1;
2130           final double term1 = dx0 * u0;
2131           final double termp0 = pz1 * u0 + py0 * u1;
2132           final double termp1 = px0 * u0;
2133           final int j = mod(++j0, fftY);
2134           int i0 = igrd0;
2135           for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
2136             final int i = mod(++i0, fftX);
2137             final int ii = iComplex3D(i, j, k, fftX, fftY);
2138             final double[] splxi = splx[ith1];
2139             final double add = splxi[0] * term0 + splxi[1] * term1;
2140             final double addi = splxi[0] * termp0 + splxi[1] * termp1;
2141             final double current = splineBuffer.get(ii);
2142             final double currenti = splineBuffer.get(ii + 1);
2143             splineBuffer.put(ii, current + add);
2144             splineBuffer.put(ii + 1, currenti + addi);
2145             // splineGrid[ii] += add;
2146             // splineGrid[ii + 1] += addi;
2147           }
2148         }
2149       }
2150     }
2151 
2152     @Override
2153     public void run(int lb, int ub) throws Exception {
2154       int ti = getThreadIndex();
2155       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2156         int[] list = gridAtomList[iSymm][ti];
2157         int n = gridAtomCount[iSymm][ti];
2158         for (int i = 0; i < n; i++) {
2159           int iAtom = list[i];
2160           gridDensity(iSymm, iAtom, lb, ub);
2161         }
2162       }
2163     }
2164 
2165     public void setUse(boolean[] use) {
2166       this.use = use;
2167     }
2168 
2169     @Override
2170     public void start() {
2171       int threadIndex = getThreadIndex();
2172       splineInducedTime[threadIndex] -= System.nanoTime();
2173       for (int i = 0; i < 3; i++) {
2174         m[0][i] = fftX * crystal.A[i][0];
2175         m[1][i] = fftY * crystal.A[i][1];
2176         m[2][i] = fftZ * crystal.A[i][2];
2177       }
2178       a00 = m[0][0];
2179       a01 = m[0][1];
2180       a02 = m[0][2];
2181       a10 = m[1][0];
2182       a11 = m[1][1];
2183       a12 = m[1][2];
2184       a20 = m[2][0];
2185       a21 = m[2][1];
2186       a22 = m[2][2];
2187     }
2188 
2189     void setInducedDipoles(double[][][] inducedDipole, double[][][] inducedDipoleCR) {
2190       this.inducedDipole = inducedDipole;
2191       this.inducedDipoleCR = inducedDipoleCR;
2192     }
2193   }
2194 
2195   /**
2196    * This class is used to obtain the reciprocal space potential and fields due to permanent
2197    * multipoles from the PME grid.
2198    *
2199    * <p>An external ParallelRegion can be used as follows: <code>
2200    * start() { permanentPhiRegion.setCartPermanentPhi(...); }
2201    * <p>
2202    * run() { execute(0, nAtoms - 1, permanentPhiLoops[threadID]); }
2203    * </code>
2204    */
2205   private class PermanentPhiRegion extends ParallelRegion {
2206 
2207     final PermanentPhiLoop[] permanentPhiLoop;
2208 
2209     private final BSplineRegion bSplineRegion;
2210     private double[][] cartPermPhi;
2211     private double[][] fracPermPhi;
2212 
2213     PermanentPhiRegion(BSplineRegion bSplineRegion) {
2214       this.bSplineRegion = bSplineRegion;
2215       permanentPhiLoop = new PermanentPhiLoop[threadCount];
2216       for (int i = 0; i < threadCount; i++) {
2217         permanentPhiLoop[i] = new PermanentPhiLoop();
2218       }
2219     }
2220 
2221     @Override
2222     public void run() {
2223       try {
2224         int threadID = getThreadIndex();
2225         execute(0, nAtoms - 1, permanentPhiLoop[threadID]);
2226       } catch (Exception e) {
2227         logger.severe(e.toString());
2228       }
2229     }
2230 
2231     void setPermanentPhi(double[][] cartPermanentPhi, double[][] fracPermanentPhi) {
2232       this.cartPermPhi = cartPermanentPhi;
2233       this.fracPermPhi = fracPermanentPhi;
2234     }
2235 
2236     public class PermanentPhiLoop extends IntegerForLoop {
2237 
2238       @Override
2239       public void finish() {
2240         int threadIndex = getThreadIndex();
2241         permanentPhiTime[threadIndex] += System.nanoTime();
2242       }
2243 
2244       @Override
2245       public void run(final int lb, final int ub) {
2246         for (int n = lb; n <= ub; n++) {
2247           final double[][] splx = bSplineRegion.splineX[0][n];
2248           final double[][] sply = bSplineRegion.splineY[0][n];
2249           final double[][] splz = bSplineRegion.splineZ[0][n];
2250           final int[] igrd = bSplineRegion.initGrid[0][n];
2251           final int igrd0 = igrd[0];
2252           final int jgrd0 = igrd[1];
2253           int k0 = igrd[2];
2254           double tuv000 = 0.0;
2255           double tuv100 = 0.0;
2256           double tuv010 = 0.0;
2257           double tuv001 = 0.0;
2258           double tuv200 = 0.0;
2259           double tuv020 = 0.0;
2260           double tuv002 = 0.0;
2261           double tuv110 = 0.0;
2262           double tuv101 = 0.0;
2263           double tuv011 = 0.0;
2264           double tuv300 = 0.0;
2265           double tuv030 = 0.0;
2266           double tuv003 = 0.0;
2267           double tuv210 = 0.0;
2268           double tuv201 = 0.0;
2269           double tuv120 = 0.0;
2270           double tuv021 = 0.0;
2271           double tuv102 = 0.0;
2272           double tuv012 = 0.0;
2273           double tuv111 = 0.0;
2274           for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
2275             final int k = mod(++k0, fftZ);
2276             int j0 = jgrd0;
2277             double tu00 = 0.0;
2278             double tu10 = 0.0;
2279             double tu01 = 0.0;
2280             double tu20 = 0.0;
2281             double tu11 = 0.0;
2282             double tu02 = 0.0;
2283             double tu30 = 0.0;
2284             double tu21 = 0.0;
2285             double tu12 = 0.0;
2286             double tu03 = 0.0;
2287             for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
2288               final int j = mod(++j0, fftY);
2289               int i0 = igrd0;
2290               double t0 = 0.0;
2291               double t1 = 0.0;
2292               double t2 = 0.0;
2293               double t3 = 0.0;
2294               for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
2295                 final int i = mod(++i0, fftX);
2296                 final int ii = iComplex3D(i, j, k, fftX, fftY);
2297                 final double tq = splineBuffer.get(ii);
2298                 final double[] splxi = splx[ith1];
2299                 t0 += tq * splxi[0];
2300                 t1 += tq * splxi[1];
2301                 t2 += tq * splxi[2];
2302                 t3 += tq * splxi[3];
2303               }
2304               final double[] splyi = sply[ith2];
2305               final double u0 = splyi[0];
2306               final double u1 = splyi[1];
2307               final double u2 = splyi[2];
2308               final double u3 = splyi[3];
2309               tu00 += t0 * u0;
2310               tu10 += t1 * u0;
2311               tu01 += t0 * u1;
2312               tu20 += t2 * u0;
2313               tu11 += t1 * u1;
2314               tu02 += t0 * u2;
2315               tu30 += t3 * u0;
2316               tu21 += t2 * u1;
2317               tu12 += t1 * u2;
2318               tu03 += t0 * u3;
2319             }
2320             final double[] splzi = splz[ith3];
2321             final double v0 = splzi[0];
2322             final double v1 = splzi[1];
2323             final double v2 = splzi[2];
2324             final double v3 = splzi[3];
2325             tuv000 += tu00 * v0;
2326             tuv100 += tu10 * v0;
2327             tuv010 += tu01 * v0;
2328             tuv001 += tu00 * v1;
2329             tuv200 += tu20 * v0;
2330             tuv020 += tu02 * v0;
2331             tuv002 += tu00 * v2;
2332             tuv110 += tu11 * v0;
2333             tuv101 += tu10 * v1;
2334             tuv011 += tu01 * v1;
2335             tuv300 += tu30 * v0;
2336             tuv030 += tu03 * v0;
2337             tuv003 += tu00 * v3;
2338             tuv210 += tu21 * v0;
2339             tuv201 += tu20 * v1;
2340             tuv120 += tu12 * v0;
2341             tuv021 += tu02 * v1;
2342             tuv102 += tu10 * v2;
2343             tuv012 += tu01 * v2;
2344             tuv111 += tu11 * v1;
2345           }
2346           double[] out = fracPermPhi[n];
2347           out[t000] = tuv000;
2348           out[t100] = tuv100;
2349           out[t010] = tuv010;
2350           out[t001] = tuv001;
2351           out[t200] = tuv200;
2352           out[t020] = tuv020;
2353           out[t002] = tuv002;
2354           out[t110] = tuv110;
2355           out[t101] = tuv101;
2356           out[t011] = tuv011;
2357           out[t300] = tuv300;
2358           out[t030] = tuv030;
2359           out[t003] = tuv003;
2360           out[t210] = tuv210;
2361           out[t201] = tuv201;
2362           out[t120] = tuv120;
2363           out[t021] = tuv021;
2364           out[t102] = tuv102;
2365           out[t012] = tuv012;
2366           out[t111] = tuv111;
2367           double[] in = out;
2368           out = cartPermPhi[n];
2369           out[0] = transformFieldMatrix[0][0] * in[0];
2370           for (int j = 1; j < 4; j++) {
2371             out[j] = 0.0;
2372             for (int k = 1; k < 4; k++) {
2373               out[j] += transformFieldMatrix[j][k] * in[k];
2374             }
2375           }
2376           for (int j = 4; j < 10; j++) {
2377             out[j] = 0.0;
2378             for (int k = 4; k < 10; k++) {
2379               out[j] += transformFieldMatrix[j][k] * in[k];
2380             }
2381           }
2382         }
2383       }
2384 
2385       @Override
2386       public IntegerSchedule schedule() {
2387         return IntegerSchedule.fixed();
2388       }
2389 
2390       @Override
2391       public void start() {
2392         int threadIndex = getThreadIndex();
2393         permanentPhiTime[threadIndex] -= System.nanoTime();
2394       }
2395     }
2396   }
2397 
2398   /**
2399    * This class is used to obtain the reciprocal space potential and fields due to induced dipoles
2400    * from the PME grid.
2401    *
2402    * <p>An external ParallelRegion can be used as follows: <code>
2403    * start() { inducedPhiRegion.setCartInducedDipolePhi(...); }
2404    * <p>
2405    * run() { execute(0, nAtoms - 1, inducedPhiLoops[threadID]); }
2406    * </code>
2407    */
2408   private class InducedPhiRegion extends ParallelRegion {
2409 
2410     final InducedPhiLoop[] inducedPhiLoops;
2411 
2412     private final BSplineRegion bSplineRegion;
2413     private double[][] cartInducedDipolePhi;
2414     private double[][] cartInducedDipolePhiCR;
2415     private double[][] fracInducedDipolePhi;
2416     private double[][] fracInducedDipolePhiCR;
2417 
2418     InducedPhiRegion(BSplineRegion bSplineRegion) {
2419       this.bSplineRegion = bSplineRegion;
2420       inducedPhiLoops = new InducedPhiLoop[threadCount];
2421       for (int i = 0; i < threadCount; i++) {
2422         inducedPhiLoops[i] = new InducedPhiLoop();
2423       }
2424     }
2425 
2426     @Override
2427     public void run() {
2428       try {
2429         int threadID = getThreadIndex();
2430         execute(0, nAtoms - 1, inducedPhiLoops[threadID]);
2431       } catch (Exception e) {
2432         logger.severe(e.toString());
2433       }
2434     }
2435 
2436     void setInducedDipolePhi(
2437         double[][] cartInducedDipolePhi,
2438         double[][] cartInducedDipolePhiCR,
2439         double[][] fracInducedDipolePhi,
2440         double[][] fracInducedDipolePhiCR) {
2441       this.cartInducedDipolePhi = cartInducedDipolePhi;
2442       this.cartInducedDipolePhiCR = cartInducedDipolePhiCR;
2443       this.fracInducedDipolePhi = fracInducedDipolePhi;
2444       this.fracInducedDipolePhiCR = fracInducedDipolePhiCR;
2445     }
2446 
2447     public class InducedPhiLoop extends IntegerForLoop {
2448 
2449       @Override
2450       public void finish() {
2451         int threadIndex = getThreadIndex();
2452         inducedPhiTime[threadIndex] += System.nanoTime();
2453       }
2454 
2455       @Override
2456       public void run(int lb, int ub) {
2457         for (int n = lb; n <= ub; n++) {
2458           final double[][] splx = bSplineRegion.splineX[0][n];
2459           final double[][] sply = bSplineRegion.splineY[0][n];
2460           final double[][] splz = bSplineRegion.splineZ[0][n];
2461           final int[] igrd = bSplineRegion.initGrid[0][n];
2462           final int igrd0 = igrd[0];
2463           final int jgrd0 = igrd[1];
2464           int k0 = igrd[2];
2465           double tuv000 = 0.0;
2466           double tuv100 = 0.0;
2467           double tuv010 = 0.0;
2468           double tuv001 = 0.0;
2469           double tuv200 = 0.0;
2470           double tuv020 = 0.0;
2471           double tuv002 = 0.0;
2472           double tuv110 = 0.0;
2473           double tuv101 = 0.0;
2474           double tuv011 = 0.0;
2475           double tuv300 = 0.0;
2476           double tuv030 = 0.0;
2477           double tuv003 = 0.0;
2478           double tuv210 = 0.0;
2479           double tuv201 = 0.0;
2480           double tuv120 = 0.0;
2481           double tuv021 = 0.0;
2482           double tuv102 = 0.0;
2483           double tuv012 = 0.0;
2484           double tuv111 = 0.0;
2485           double tuv000p = 0.0;
2486           double tuv100p = 0.0;
2487           double tuv010p = 0.0;
2488           double tuv001p = 0.0;
2489           double tuv200p = 0.0;
2490           double tuv020p = 0.0;
2491           double tuv002p = 0.0;
2492           double tuv110p = 0.0;
2493           double tuv101p = 0.0;
2494           double tuv011p = 0.0;
2495           double tuv300p = 0.0;
2496           double tuv030p = 0.0;
2497           double tuv003p = 0.0;
2498           double tuv210p = 0.0;
2499           double tuv201p = 0.0;
2500           double tuv120p = 0.0;
2501           double tuv021p = 0.0;
2502           double tuv102p = 0.0;
2503           double tuv012p = 0.0;
2504           double tuv111p = 0.0;
2505           for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
2506             final int k = mod(++k0, fftZ);
2507             int j0 = jgrd0;
2508             double tu00 = 0.0;
2509             double tu10 = 0.0;
2510             double tu01 = 0.0;
2511             double tu20 = 0.0;
2512             double tu11 = 0.0;
2513             double tu02 = 0.0;
2514             double tu30 = 0.0;
2515             double tu21 = 0.0;
2516             double tu12 = 0.0;
2517             double tu03 = 0.0;
2518             double tu00p = 0.0;
2519             double tu10p = 0.0;
2520             double tu01p = 0.0;
2521             double tu20p = 0.0;
2522             double tu11p = 0.0;
2523             double tu02p = 0.0;
2524             double tu30p = 0.0;
2525             double tu21p = 0.0;
2526             double tu12p = 0.0;
2527             double tu03p = 0.0;
2528             for (int ith2 = 0; ith2 < bSplineOrder; ith2++) {
2529               final int j = mod(++j0, fftY);
2530               int i0 = igrd0;
2531               double t0 = 0.0;
2532               double t1 = 0.0;
2533               double t2 = 0.0;
2534               double t3 = 0.0;
2535               double t0p = 0.0;
2536               double t1p = 0.0;
2537               double t2p = 0.0;
2538               double t3p = 0.0;
2539               for (int ith1 = 0; ith1 < bSplineOrder; ith1++) {
2540                 final int i = mod(++i0, fftX);
2541                 final int ii = iComplex3D(i, j, k, fftX, fftY);
2542                 // final double tq = splineGrid[ii];
2543                 // final double tp = splineGrid[ii + 1];
2544                 final double tq = splineBuffer.get(ii);
2545                 final double tp = splineBuffer.get(ii + 1);
2546                 final double[] splxi = splx[ith1];
2547                 t0 += tq * splxi[0];
2548                 t1 += tq * splxi[1];
2549                 t2 += tq * splxi[2];
2550                 t3 += tq * splxi[3];
2551                 t0p += tp * splxi[0];
2552                 t1p += tp * splxi[1];
2553                 t2p += tp * splxi[2];
2554                 t3p += tp * splxi[3];
2555               }
2556               final double[] splyi = sply[ith2];
2557               final double u0 = splyi[0];
2558               final double u1 = splyi[1];
2559               final double u2 = splyi[2];
2560               final double u3 = splyi[3];
2561               tu00 += t0 * u0;
2562               tu10 += t1 * u0;
2563               tu01 += t0 * u1;
2564               tu20 += t2 * u0;
2565               tu11 += t1 * u1;
2566               tu02 += t0 * u2;
2567               tu30 += t3 * u0;
2568               tu21 += t2 * u1;
2569               tu12 += t1 * u2;
2570               tu03 += t0 * u3;
2571               tu00p += t0p * u0;
2572               tu10p += t1p * u0;
2573               tu01p += t0p * u1;
2574               tu20p += t2p * u0;
2575               tu11p += t1p * u1;
2576               tu02p += t0p * u2;
2577               tu30p += t3p * u0;
2578               tu21p += t2p * u1;
2579               tu12p += t1p * u2;
2580               tu03p += t0p * u3;
2581             }
2582             final double[] splzi = splz[ith3];
2583             final double v0 = splzi[0];
2584             final double v1 = splzi[1];
2585             final double v2 = splzi[2];
2586             final double v3 = splzi[3];
2587             tuv000 += tu00 * v0;
2588             tuv100 += tu10 * v0;
2589             tuv010 += tu01 * v0;
2590             tuv001 += tu00 * v1;
2591             tuv200 += tu20 * v0;
2592             tuv020 += tu02 * v0;
2593             tuv002 += tu00 * v2;
2594             tuv110 += tu11 * v0;
2595             tuv101 += tu10 * v1;
2596             tuv011 += tu01 * v1;
2597             tuv300 += tu30 * v0;
2598             tuv030 += tu03 * v0;
2599             tuv003 += tu00 * v3;
2600             tuv210 += tu21 * v0;
2601             tuv201 += tu20 * v1;
2602             tuv120 += tu12 * v0;
2603             tuv021 += tu02 * v1;
2604             tuv102 += tu10 * v2;
2605             tuv012 += tu01 * v2;
2606             tuv111 += tu11 * v1;
2607             tuv000p += tu00p * v0;
2608             tuv100p += tu10p * v0;
2609             tuv010p += tu01p * v0;
2610             tuv001p += tu00p * v1;
2611             tuv200p += tu20p * v0;
2612             tuv020p += tu02p * v0;
2613             tuv002p += tu00p * v2;
2614             tuv110p += tu11p * v0;
2615             tuv101p += tu10p * v1;
2616             tuv011p += tu01p * v1;
2617             tuv300p += tu30p * v0;
2618             tuv030p += tu03p * v0;
2619             tuv003p += tu00p * v3;
2620             tuv210p += tu21p * v0;
2621             tuv201p += tu20p * v1;
2622             tuv120p += tu12p * v0;
2623             tuv021p += tu02p * v1;
2624             tuv102p += tu10p * v2;
2625             tuv012p += tu01p * v2;
2626             tuv111p += tu11p * v1;
2627           }
2628           double[] out = fracInducedDipolePhi[n];
2629           out[t000] = tuv000;
2630           out[t100] = tuv100;
2631           out[t010] = tuv010;
2632           out[t001] = tuv001;
2633           out[t200] = tuv200;
2634           out[t020] = tuv020;
2635           out[t002] = tuv002;
2636           out[t110] = tuv110;
2637           out[t101] = tuv101;
2638           out[t011] = tuv011;
2639           out[t300] = tuv300;
2640           out[t030] = tuv030;
2641           out[t003] = tuv003;
2642           out[t210] = tuv210;
2643           out[t201] = tuv201;
2644           out[t120] = tuv120;
2645           out[t021] = tuv021;
2646           out[t102] = tuv102;
2647           out[t012] = tuv012;
2648           out[t111] = tuv111;
2649           double[] in = out;
2650           out = cartInducedDipolePhi[n];
2651           out[0] = transformFieldMatrix[0][0] * in[0];
2652           for (int j = 1; j < 4; j++) {
2653             out[j] = 0.0;
2654             for (int k = 1; k < 4; k++) {
2655               out[j] += transformFieldMatrix[j][k] * in[k];
2656             }
2657           }
2658           for (int j = 4; j < 10; j++) {
2659             out[j] = 0.0;
2660             for (int k = 4; k < 10; k++) {
2661               out[j] += transformFieldMatrix[j][k] * in[k];
2662             }
2663           }
2664 
2665           out = fracInducedDipolePhiCR[n];
2666           out[t000] = tuv000p;
2667           out[t100] = tuv100p;
2668           out[t010] = tuv010p;
2669           out[t001] = tuv001p;
2670           out[t200] = tuv200p;
2671           out[t020] = tuv020p;
2672           out[t002] = tuv002p;
2673           out[t110] = tuv110p;
2674           out[t101] = tuv101p;
2675           out[t011] = tuv011p;
2676           out[t300] = tuv300p;
2677           out[t030] = tuv030p;
2678           out[t003] = tuv003p;
2679           out[t210] = tuv210p;
2680           out[t201] = tuv201p;
2681           out[t120] = tuv120p;
2682           out[t021] = tuv021p;
2683           out[t102] = tuv102p;
2684           out[t012] = tuv012p;
2685           out[t111] = tuv111p;
2686           in = out;
2687           out = cartInducedDipolePhiCR[n];
2688           out[0] = transformFieldMatrix[0][0] * in[0];
2689           for (int j = 1; j < 4; j++) {
2690             out[j] = 0.0;
2691             for (int k = 1; k < 4; k++) {
2692               out[j] += transformFieldMatrix[j][k] * in[k];
2693             }
2694           }
2695           for (int j = 4; j < 10; j++) {
2696             out[j] = 0.0;
2697             for (int k = 4; k < 10; k++) {
2698               out[j] += transformFieldMatrix[j][k] * in[k];
2699             }
2700           }
2701         }
2702       }
2703 
2704       @Override
2705       public IntegerSchedule schedule() {
2706         return IntegerSchedule.fixed();
2707       }
2708 
2709       @Override
2710       public void start() {
2711         int threadIndex = getThreadIndex();
2712         inducedPhiTime[threadIndex] -= System.nanoTime();
2713       }
2714     }
2715   }
2716 }