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.xray;
39  
40  import static ffx.crystal.SymOp.applyTransSymRot;
41  import static ffx.numerics.fft.Complex3D.iComplex3D;
42  import static ffx.numerics.math.ScalarMath.mod;
43  import static java.lang.String.format;
44  import static java.lang.System.arraycopy;
45  import static java.util.Arrays.fill;
46  import static org.apache.commons.math3.util.FastMath.abs;
47  import static org.apache.commons.math3.util.FastMath.exp;
48  import static org.apache.commons.math3.util.FastMath.floor;
49  import static org.apache.commons.math3.util.FastMath.max;
50  import static org.apache.commons.math3.util.FastMath.min;
51  import static org.apache.commons.math3.util.FastMath.pow;
52  import static org.apache.commons.math3.util.FastMath.sqrt;
53  
54  import edu.rit.pj.IntegerForLoop;
55  import edu.rit.pj.IntegerSchedule;
56  import edu.rit.pj.ParallelRegion;
57  import edu.rit.pj.ParallelTeam;
58  import edu.rit.pj.reduction.SharedIntegerArray;
59  import ffx.crystal.Crystal;
60  import ffx.crystal.HKL;
61  import ffx.crystal.ReflectionList;
62  import ffx.crystal.Resolution;
63  import ffx.crystal.SymOp;
64  import ffx.numerics.fft.Complex;
65  import ffx.numerics.fft.Complex3DParallel;
66  import ffx.numerics.math.ComplexNumber;
67  import ffx.potential.bonded.Atom;
68  import ffx.potential.nonbonded.RowLoop;
69  import ffx.potential.nonbonded.RowRegion;
70  import ffx.potential.nonbonded.SliceLoop;
71  import ffx.potential.nonbonded.SliceRegion;
72  import ffx.potential.nonbonded.SpatialDensityLoop;
73  import ffx.potential.nonbonded.SpatialDensityRegion;
74  import ffx.xray.RefinementMinimize.RefinementMode;
75  import java.util.ArrayList;
76  import java.util.List;
77  import java.util.logging.Level;
78  import java.util.logging.Logger;
79  
80  /**
81   * Structure factor calculation (including bulk solvent structure factors)
82   *
83   * @author Timothy D. Fenn
84   * @see <a href="http://dx.doi.org/10.1107/S0567739473000458" target="_blank"> L. F. Ten Eyck, Acta
85   *     Cryst. (1973). A29, 183-191.</a>
86   * @see <a href="http://dx.doi.org/10.1107/S0567739477001211" target="_blank"> L. F. Ten Eyck, Acta
87   *     Cryst. (1977). A33, 486-492.</a>
88   * @see <a href="http://dx.doi.org/10.1107/S0365110X55001862" target="_blank"> J. Waser, Acta Cryst.
89   *     (1955). 8, 595.</a>
90   * @see <a href="http://dx.doi.org/10.1107/S0108767388009183" target="_blank"> A. T. Brunger, Acta
91   *     Cryst. (1989). A45, 42-50.</a>
92   * @see <a href="http://dx.doi.org/10.1107/97809553602060000551" target="_blank"> G. Bricogne, Int.
93   *     Tables Cryst. (2006). Vol. B, ch. 1.3, pp. 25-98.</a>
94   * @see <a href="http://dx.doi.org/10.1002/jcc.1032" target="_blank"> J. A. Grant, B. T. Pickup, A.
95   *     Nicholls, J. Comp. Chem. (2001). 22, 608-640</a>
96   * @see <a href="http://dx.doi.org/10.1006/jmbi.1994.1633" target="_blank"> J. S. Jiang, A. T.
97   *     Brunger, JMB (1994) 243, 100-115.</a>
98   * @see <a href="http://dx.doi.org/10.1107/S0907444910031045" target="_blank"> T.D. Fenn, M. J.
99   *     Schnieders, A. T. Brunger, Acta Cryst. (2010). D66, 1024-1031.</a>
100  * @since 1.0
101  */
102 public class CrystalReciprocalSpace {
103 
104   private static final Logger logger = Logger.getLogger(CrystalReciprocalSpace.class.getName());
105   private static final double toSeconds = 1.0e-9;
106   private final boolean neutron;
107   private final double bAdd;
108   private final double fftScale;
109   private final double[] densityGrid;
110   private final double[] solventGrid;
111   private final int nSymm;
112   private final int bulkNSymm;
113   private final int nAtoms;
114   private final int fftX, fftY, fftZ;
115   private final double ifftX, ifftY, ifftZ;
116   private final int halfFFTX;
117   private final int complexFFT3DSpace;
118   private final int threadCount;
119   private final int bufferSize = 3;
120   // Slice Scheduling
121   private final SharedIntegerArray optSliceWeightAtomic;
122   private final SharedIntegerArray optSliceWeightSolvent;
123   private final SharedIntegerArray optSliceWeightBulkSolvent;
124   private final int[] previousOptSliceWeightAtomic;
125   private final int[] previousOptSliceWeightSolvent;
126   private final int[] previousOptSliceWeightBulkSolvent;
127   private final SliceSchedule atomicSliceSchedule;
128   private final SliceSchedule solventSliceSchedule;
129   private final SliceSchedule bulkSolventSliceSchedule;
130   // Row Scheduling
131   private final SharedIntegerArray optRowWeightAtomic;
132   private final SharedIntegerArray optRowWeightSolvent;
133   private final SharedIntegerArray optRowWeightBulkSolvent;
134   private final int[] previousOptRowWeightAtomic;
135   private final int[] previousOptRowWeightSolvent;
136   private final int[] previousOptRowWeightBulkSolvent;
137   private final RowSchedule atomicRowSchedule;
138   private final RowSchedule solventRowSchedule;
139   private final RowSchedule bulkSolventRowSchedule;
140 
141   private final ParallelTeam parallelTeam;
142   private final Atom[] atoms;
143   private final Crystal crystal;
144   private final ReflectionList reflectionList;
145   private final FormFactor[][] atomFormFactors;
146   private final FormFactor[][] solventFormFactors;
147   private final GridMethod gridMethod;
148   private final SharedIntegerArray optAtomicGradientWeight;
149   private final int[] previousOptAtomicGradientWeight;
150   private final GradientSchedule atomicGradientSchedule;
151   /**
152    * Parallelization of putting atomic form factors onto the 3D grid using a 3D spatial
153    * decomposition.
154    */
155   private final SpatialDensityRegion atomicDensityRegion;
156 
157   private final SpatialDensityRegion solventDensityRegion;
158   /**
159    * Parallelization of putting atomic form factors onto the 3D grid using a slice-based
160    * decomposition.
161    */
162   private final SliceRegion atomicSliceRegion;
163 
164   private final AtomicSliceLoop[] atomicSliceLoops;
165   private final SliceRegion solventSliceRegion;
166   private final SolventSliceLoop[] solventSliceLoops;
167   private final BulkSolventSliceLoop[] bulkSolventSliceLoops;
168   /**
169    * Parallelization of putting atomic form factors onto the 3D grid using a row-based
170    * decomposition.
171    */
172   private final RowRegion atomicRowRegion;
173 
174   private final AtomicRowLoop[] atomicRowLoops;
175   private final RowRegion solventRowRegion;
176   private final SolventRowLoop[] solventRowLoops;
177   private final BulkSolventRowLoop[] bulkSolventRowLoops;
178   /** Parallelization over atomic structure factor loops. */
179   private final InitRegion initRegion;
180 
181   private final ExtractRegion extractRegion;
182   private final AtomicScaleRegion atomicScaleRegion;
183   private final AtomicGradientRegion atomicGradientRegion;
184   /** Parallelization over solvent structure factor loops. */
185   private final BabinetRegion babinetRegion;
186 
187   private final SolventGridRegion solventGridRegion;
188   private final SolventScaleRegion solventScaleRegion;
189   private final SolventGradientRegion solventGradientRegion;
190   private final Complex3DParallel complexFFT3D;
191   protected SolventModel solventModel;
192   boolean lambdaTerm = false;
193   double solventA;
194   double solventB;
195   private final boolean solvent;
196   private boolean useThreeGaussians = true;
197   // Not final for purposes of finite differences
198   private double[][][] coordinates;
199   private double weight = 1.0;
200   private final int aRadGrid;
201   /** If the "Native Environment Approximation" is true, the "use" flag is ignored. */
202   private boolean nativeEnvironmentApproximation = false;
203 
204   /**
205    * Crystal Reciprocal Space constructor, assumes this is not a bulk solvent mask and is not a
206    * neutron data set
207    *
208    * @param reflectionList the {@link ffx.crystal.ReflectionList} to fill with structure factors
209    * @param atoms array of {@link ffx.potential.bonded.Atom atoms} for structure factor
210    *     computation
211    * @param fftTeam {@link edu.rit.pj.ParallelTeam} for parallelization
212    * @param parallelTeam {@link edu.rit.pj.ParallelTeam} for parallelization
213    */
214   public CrystalReciprocalSpace(
215       ReflectionList reflectionList,
216       Atom[] atoms,
217       ParallelTeam fftTeam,
218       ParallelTeam parallelTeam) {
219     this(
220         reflectionList,
221         atoms,
222         fftTeam,
223         parallelTeam,
224         false,
225         false,
226         SolventModel.POLYNOMIAL,
227         GridMethod.SLICE);
228   }
229 
230   /**
231    * Crystal Reciprocal Space constructor, assumes this is not a neutron data set and implements a
232    * polynomial bulk solvent mask if needed
233    *
234    * @param reflectionList the {@link ffx.crystal.ReflectionList} to fill with structure factors
235    * @param atoms array of {@link ffx.potential.bonded.Atom atoms} for structure factor
236    *     computation
237    * @param fftTeam {@link edu.rit.pj.ParallelTeam} for parallelization
238    * @param parallelTeam {@link edu.rit.pj.ParallelTeam} for parallelization
239    * @param solventMask true if this is a bulk solvent mask
240    */
241   public CrystalReciprocalSpace(
242       ReflectionList reflectionList,
243       Atom[] atoms,
244       ParallelTeam fftTeam,
245       ParallelTeam parallelTeam,
246       boolean solventMask) {
247     this(
248         reflectionList,
249         atoms,
250         fftTeam,
251         parallelTeam,
252         solventMask,
253         false,
254         SolventModel.POLYNOMIAL,
255         GridMethod.SLICE);
256   }
257 
258   /**
259    * Crystal Reciprocal Space constructor, assumes a polynomial bulk solvent mask if needed
260    *
261    * @param reflectionList the {@link ffx.crystal.ReflectionList} to fill with structure factors
262    * @param atoms array of {@link ffx.potential.bonded.Atom atoms} for structure factor
263    *     computation
264    * @param fftTeam {@link edu.rit.pj.ParallelTeam} for parallelization
265    * @param parallelTeam {@link edu.rit.pj.ParallelTeam} for parallelization
266    * @param solventMask true if this is a bulk solvent mask
267    * @param neutron true if this is a neutron structure
268    */
269   public CrystalReciprocalSpace(
270       ReflectionList reflectionList,
271       Atom[] atoms,
272       ParallelTeam fftTeam,
273       ParallelTeam parallelTeam,
274       boolean solventMask,
275       boolean neutron) {
276     this(
277         reflectionList,
278         atoms,
279         fftTeam,
280         parallelTeam,
281         solventMask,
282         neutron,
283         SolventModel.POLYNOMIAL,
284         GridMethod.SLICE);
285   }
286 
287   /**
288    * Crystal Reciprocal Space constructor, all parameters provided
289    *
290    * @param reflectionlist the {@link ffx.crystal.ReflectionList} to fill with structure factors
291    * @param atoms array of {@link ffx.potential.bonded.Atom atoms} for structure factor
292    *     computation
293    * @param fftTeam {@link edu.rit.pj.ParallelTeam} for parallelization
294    * @param parallelTeam {@link edu.rit.pj.ParallelTeam} for parallelization
295    * @param solventMask true if this is a bulk solvent mask
296    * @param neutron true if this is a neutron structure
297    * @param solventModel bulk solvent model type
298    * @param gridMethod parallel method for putting density onto the FFT grid.
299    * @see CrystalReciprocalSpace.SolventModel
300    */
301   public CrystalReciprocalSpace(
302       ReflectionList reflectionlist,
303       Atom[] atoms,
304       ParallelTeam fftTeam,
305       ParallelTeam parallelTeam,
306       boolean solventMask,
307       boolean neutron,
308       SolventModel solventModel,
309       GridMethod gridMethod) {
310     this.reflectionList = reflectionlist;
311     this.atoms = atoms;
312     this.parallelTeam = parallelTeam;
313     this.solvent = solventMask;
314     this.neutron = neutron;
315     this.solventModel = solventModel;
316     this.gridMethod = gridMethod;
317 
318     nAtoms = atoms.length;
319     crystal = reflectionlist.crystal;
320     nSymm = 1;
321     threadCount = parallelTeam.getThreadCount();
322 
323     // Necssary for the bulk-solvent expansion!
324     bulkNSymm = crystal.spaceGroup.symOps.size();
325     Resolution resolution = reflectionlist.resolution;
326     double gridFactor = resolution.samplingLimit() / 2.0;
327     double gridStep = resolution.resolutionLimit() * gridFactor;
328     // Set default FFT grid size from unit cell dimensions.
329     int nX = (int) Math.floor(crystal.a / gridStep) + 1;
330     int nY = (int) Math.floor(crystal.b / gridStep) + 1;
331     int nZ = (int) Math.floor(crystal.c / gridStep) + 1;
332     if (nX % 2 != 0) {
333       nX += 1;
334     }
335     // Use preferred dimensions.
336     while (!Complex.preferredDimension(nX)) {
337       nX += 2;
338     }
339     while (!Complex.preferredDimension(nY)) {
340       nY += 1;
341     }
342     while (!Complex.preferredDimension(nZ)) {
343       nZ += 1;
344     }
345 
346     fftX = nX;
347     fftY = nY;
348     fftZ = nZ;
349     halfFFTX = nX / 2;
350     ifftX = 1.0 / (double) fftX;
351     ifftY = 1.0 / (double) fftY;
352     ifftZ = 1.0 / (double) fftZ;
353     fftScale = crystal.volume;
354     complexFFT3DSpace = fftX * fftY * fftZ * 2;
355     densityGrid = new double[complexFFT3DSpace];
356 
357     // Slice Method Variables
358     optSliceWeightAtomic = new SharedIntegerArray(fftZ);
359     optSliceWeightSolvent = new SharedIntegerArray(fftZ);
360     optSliceWeightBulkSolvent = new SharedIntegerArray(fftZ);
361     previousOptSliceWeightAtomic = new int[fftZ];
362     previousOptSliceWeightSolvent = new int[fftZ];
363     previousOptSliceWeightBulkSolvent = new int[fftZ];
364 
365     atomicSliceSchedule = new SliceSchedule(threadCount, fftZ);
366     solventSliceSchedule = new SliceSchedule(threadCount, fftZ);
367     bulkSolventSliceSchedule = new SliceSchedule(threadCount, fftZ);
368 
369     // Row Method Variables
370     optRowWeightAtomic = new SharedIntegerArray(fftZ * fftY);
371     optRowWeightSolvent = new SharedIntegerArray(fftZ * fftY);
372     optRowWeightBulkSolvent = new SharedIntegerArray(fftZ * fftY);
373     previousOptRowWeightAtomic = new int[fftZ * fftY];
374     previousOptRowWeightSolvent = new int[fftZ * fftY];
375     previousOptRowWeightBulkSolvent = new int[fftZ * fftY];
376     atomicRowSchedule = new RowSchedule(threadCount, fftZ, fftY);
377     solventRowSchedule = new RowSchedule(threadCount, fftZ, fftY);
378     bulkSolventRowSchedule = new RowSchedule(threadCount, fftZ, fftY);
379 
380     optAtomicGradientWeight = new SharedIntegerArray(nAtoms);
381     previousOptAtomicGradientWeight = new int[nAtoms];
382     atomicGradientSchedule = new GradientSchedule(threadCount, nAtoms);
383 
384     if (solvent) {
385       bAdd = 0.0;
386       atomFormFactors = null;
387       double vdwr;
388       switch (solventModel) {
389         case BINARY:
390           solventA = 1.0;
391           solventB = 1.0;
392           solventGrid = new double[complexFFT3DSpace];
393           solventFormFactors = new SolventBinaryFormFactor[bulkNSymm][nAtoms];
394           for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
395             for (int i = 0; i < nAtoms; i++) {
396               vdwr = atoms[i].getVDWType().radius * 0.5;
397               solventFormFactors[iSymm][i] = new SolventBinaryFormFactor(atoms[i], vdwr + solventA);
398             }
399           }
400           break;
401         case GAUSSIAN:
402           solventA = 11.5;
403           solventB = 0.55;
404           solventGrid = new double[complexFFT3DSpace];
405           solventFormFactors = new SolventGaussFormFactor[bulkNSymm][nAtoms];
406           for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
407             for (int i = 0; i < nAtoms; i++) {
408               vdwr = atoms[i].getVDWType().radius * 0.5;
409               solventFormFactors[iSymm][i] = new SolventGaussFormFactor(atoms[i], vdwr * solventB);
410             }
411           }
412           break;
413         case POLYNOMIAL:
414           solventA = 0.0;
415           solventB = 0.8;
416           solventGrid = new double[complexFFT3DSpace];
417           solventFormFactors = new SolventPolyFormFactor[bulkNSymm][nAtoms];
418           for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
419             for (int i = 0; i < nAtoms; i++) {
420               vdwr = atoms[i].getVDWType().radius * 0.5;
421               solventFormFactors[iSymm][i] =
422                   new SolventPolyFormFactor(atoms[i], vdwr + solventA, solventB);
423             }
424           }
425           break;
426         case NONE:
427         default:
428           solventGrid = null;
429           solventFormFactors = null;
430           break;
431       }
432     } else {
433       bAdd = 2.0;
434       if (neutron) {
435         atomFormFactors = new NeutronFormFactor[bulkNSymm][nAtoms];
436         for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
437           for (int i = 0; i < nAtoms; i++) {
438             atomFormFactors[iSymm][i] = new NeutronFormFactor(atoms[i], bAdd);
439           }
440         }
441       } else {
442         atomFormFactors = new XRayFormFactor[bulkNSymm][nAtoms];
443         for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
444           for (int i = 0; i < nAtoms; i++) {
445             atomFormFactors[iSymm][i] = new XRayFormFactor(atoms[i], useThreeGaussians, bAdd);
446           }
447         }
448       }
449       solventGrid = null;
450       solventFormFactors = null;
451     }
452 
453     // Determine number of grid points to sample density on
454     double aRad = -1.0;
455     for (Atom a : atoms) {
456       double vdwr = a.getVDWType().radius * 0.5;
457       if (!solvent) {
458         aRad = Math.max(aRad, a.getFormFactorWidth());
459       } else {
460         switch (solventModel) {
461           case BINARY:
462             aRad = Math.max(aRad, vdwr + solventA + 0.2);
463             break;
464           case GAUSSIAN:
465             aRad = Math.max(aRad, vdwr * solventB + 2.0);
466             break;
467           case POLYNOMIAL:
468             aRad = Math.max(aRad, vdwr + solventB + 0.2);
469             break;
470           case NONE:
471           default:
472             aRad = 0.0;
473         }
474       }
475     }
476     aRadGrid = (int) Math.floor(aRad * nX / crystal.a) + 1;
477 
478     // local copy of coordinates - note use of bulknsym
479     coordinates = new double[bulkNSymm][3][nAtoms];
480     List<SymOp> symops = crystal.spaceGroup.symOps;
481     double[] xyz = new double[3];
482     for (int i = 0; i < bulkNSymm; i++) {
483       double[] x = coordinates[i][0];
484       double[] y = coordinates[i][1];
485       double[] z = coordinates[i][2];
486       for (int j = 0; j < nAtoms; j++) {
487         Atom aj = atoms[j];
488         crystal.applySymOp(aj.getXYZ(null), xyz, symops.get(i));
489         x[j] = xyz[0];
490         y[j] = xyz[1];
491         z[j] = xyz[2];
492       }
493     }
494 
495     if (logger.isLoggable(Level.INFO)) {
496       StringBuilder sb = new StringBuilder();
497       if (solvent) {
498         sb.append("  Bulk Solvent Grid\n");
499         sb.append(format("  Bulk solvent model type:           %8s\n", solventModel));
500       } else {
501         sb.append("  Atomic Scattering Grid\n");
502       }
503       sb.append(format("  Form factor grid points:             %8d\n", aRadGrid * 2));
504       sb.append(format("  Form factor grid diameter:           %8.3f\n", aRad * 2));
505       sb.append(format("  Grid density:                        %8.3f\n", 1.0 / gridStep));
506       sb.append(format("  Grid dimensions:                (%3d,%3d,%3d)\n", fftX, fftY, fftZ));
507       sb.append(format("  Grid method:                         %8s\n", gridMethod));
508       logger.info(sb.toString());
509     }
510 
511     if (solvent) {
512       int minWork = nSymm;
513       if (solventModel != SolventModel.NONE) {
514         solventGradientRegion =
515             new SolventGradientRegion(RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES);
516         // Null out scattering instances used without bulk solvent.
517         atomicGradientRegion = null;
518         atomicSliceLoops = null;
519         atomicRowLoops = null;
520         switch (gridMethod) {
521           case SPATIAL:
522             atomicDensityRegion =
523                 new SpatialDensityRegion(
524                     fftX,
525                     fftY,
526                     fftZ,
527                     densityGrid,
528                     (aRadGrid + 2) * 2,
529                     nSymm,
530                     minWork,
531                     threadCount,
532                     crystal,
533                     atoms,
534                     coordinates);
535             solventDensityRegion =
536                 new BulkSolventDensityRegion(
537                     fftX,
538                     fftY,
539                     fftZ,
540                     solventGrid,
541                     (aRadGrid + 2) * 2,
542                     bulkNSymm,
543                     minWork,
544                     threadCount,
545                     crystal,
546                     atoms,
547                     coordinates,
548                     4.0,
549                     parallelTeam);
550             if (solventModel == SolventModel.GAUSSIAN) {
551               atomicDensityRegion.setInitValue(0.0);
552             } else {
553               atomicDensityRegion.setInitValue(1.0);
554             }
555             SolventDensityLoop[] solventDensityLoops = new SolventDensityLoop[threadCount];
556             SolventDensityLoop[] bulkSolventDensityLoops = new SolventDensityLoop[threadCount];
557             for (int i = 0; i < threadCount; i++) {
558               solventDensityLoops[i] = new SolventDensityLoop(atomicDensityRegion);
559               bulkSolventDensityLoops[i] = new SolventDensityLoop(solventDensityRegion);
560             }
561             atomicDensityRegion.setDensityLoop(solventDensityLoops);
562             solventDensityRegion.setDensityLoop(bulkSolventDensityLoops);
563 
564             atomicSliceRegion = null;
565             solventSliceRegion = null;
566             bulkSolventSliceLoops = null;
567             solventSliceLoops = null;
568 
569             atomicRowRegion = null;
570             solventRowRegion = null;
571             bulkSolventRowLoops = null;
572             solventRowLoops = null;
573             break;
574           case ROW:
575             atomicRowRegion =
576                 new RowRegion(
577                     fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
578             solventRowRegion =
579                 new BulkSolventRowRegion(
580                     fftX,
581                     fftY,
582                     fftZ,
583                     solventGrid,
584                     bulkNSymm,
585                     threadCount,
586                     crystal,
587                     atoms,
588                     coordinates,
589                     4.0,
590                     parallelTeam);
591             if (solventModel == SolventModel.GAUSSIAN) {
592               atomicRowRegion.setInitValue(0.0);
593             } else {
594               atomicRowRegion.setInitValue(1.0);
595             }
596             solventRowLoops = new SolventRowLoop[threadCount];
597             bulkSolventRowLoops = new BulkSolventRowLoop[threadCount];
598             for (int i = 0; i < threadCount; i++) {
599               solventRowLoops[i] = new SolventRowLoop(atomicRowRegion);
600               bulkSolventRowLoops[i] = new BulkSolventRowLoop(solventRowRegion);
601             }
602             atomicRowRegion.setDensityLoop(solventRowLoops);
603             solventRowRegion.setDensityLoop(bulkSolventRowLoops);
604             atomicDensityRegion = null;
605             solventDensityRegion = null;
606             atomicSliceRegion = null;
607             solventSliceRegion = null;
608             bulkSolventSliceLoops = null;
609             solventSliceLoops = null;
610             break;
611           case SLICE:
612           default:
613             atomicSliceRegion =
614                 new SliceRegion(
615                     fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
616             solventSliceRegion =
617                 new BulkSolventSliceRegion(
618                     fftX,
619                     fftY,
620                     fftZ,
621                     solventGrid,
622                     bulkNSymm,
623                     threadCount,
624                     crystal,
625                     atoms,
626                     coordinates,
627                     4.0,
628                     parallelTeam);
629             if (solventModel == SolventModel.GAUSSIAN) {
630               atomicSliceRegion.setInitValue(0.0);
631             } else {
632               atomicSliceRegion.setInitValue(1.0);
633             }
634             solventSliceLoops = new SolventSliceLoop[threadCount];
635             bulkSolventSliceLoops = new BulkSolventSliceLoop[threadCount];
636             for (int i = 0; i < threadCount; i++) {
637               solventSliceLoops[i] = new SolventSliceLoop(atomicSliceRegion);
638               bulkSolventSliceLoops[i] = new BulkSolventSliceLoop(solventSliceRegion);
639             }
640             atomicSliceRegion.setDensityLoop(solventSliceLoops);
641             solventSliceRegion.setDensityLoop(bulkSolventSliceLoops);
642 
643             atomicDensityRegion = null;
644             solventDensityRegion = null;
645             atomicRowRegion = null;
646             solventRowRegion = null;
647             bulkSolventRowLoops = null;
648             solventRowLoops = null;
649         }
650       } else {
651         atomicDensityRegion = null;
652         solventDensityRegion = null;
653         atomicGradientRegion = null;
654         solventGradientRegion = null;
655         atomicSliceRegion = null;
656         atomicSliceLoops = null;
657         solventSliceRegion = null;
658         bulkSolventSliceLoops = null;
659         solventSliceLoops = null;
660         atomicRowRegion = null;
661         atomicRowLoops = null;
662         solventRowRegion = null;
663         bulkSolventRowLoops = null;
664         solventRowLoops = null;
665       }
666     } else {
667       atomicGradientRegion =
668           new AtomicGradientRegion(RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES);
669 
670       // Null out bulk solvent instances.
671       solventDensityRegion = null;
672       solventSliceRegion = null;
673       solventSliceLoops = null;
674       bulkSolventSliceLoops = null;
675       solventRowRegion = null;
676       solventRowLoops = null;
677       bulkSolventRowLoops = null;
678       solventGradientRegion = null;
679 
680       /*
681        * Create nSymm pieces of work per thread; the empty pieces will be
682        * removed leaving 1 piece of work per thread.
683        */
684       switch (gridMethod) {
685         case SPATIAL:
686           int minWork = nSymm;
687           atomicDensityRegion =
688               new SpatialDensityRegion(
689                   fftX,
690                   fftY,
691                   fftZ,
692                   densityGrid,
693                   (aRadGrid + 2) * 2,
694                   nSymm,
695                   minWork,
696                   threadCount,
697                   crystal,
698                   atoms,
699                   coordinates);
700           atomicDensityRegion.setInitValue(0.0);
701           AtomicDensityLoop[] atomicDensityLoops = new AtomicDensityLoop[threadCount];
702           for (int i = 0; i < threadCount; i++) {
703             atomicDensityLoops[i] = new AtomicDensityLoop(atomicDensityRegion);
704           }
705           atomicDensityRegion.setDensityLoop(atomicDensityLoops);
706           atomicSliceRegion = null;
707           atomicSliceLoops = null;
708           atomicRowRegion = null;
709           atomicRowLoops = null;
710           break;
711 
712         case ROW:
713           atomicRowRegion =
714               new RowRegion(fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
715           atomicRowRegion.setInitValue(0.0);
716           atomicRowLoops = new AtomicRowLoop[threadCount];
717           for (int i = 0; i < threadCount; i++) {
718             atomicRowLoops[i] = new AtomicRowLoop(atomicRowRegion);
719           }
720           atomicRowRegion.setDensityLoop(atomicRowLoops);
721           atomicDensityRegion = null;
722           atomicSliceRegion = null;
723           atomicSliceLoops = null;
724           break;
725 
726         case SLICE:
727         default:
728           atomicSliceRegion =
729               new SliceRegion(
730                   fftX, fftY, fftZ, densityGrid, nSymm, threadCount, atoms, coordinates);
731           atomicSliceRegion.setInitValue(0.0);
732           atomicSliceLoops = new AtomicSliceLoop[threadCount];
733           for (int i = 0; i < threadCount; i++) {
734             atomicSliceLoops[i] = new AtomicSliceLoop(atomicSliceRegion);
735           }
736           atomicSliceRegion.setDensityLoop(atomicSliceLoops);
737           atomicDensityRegion = null;
738           atomicRowRegion = null;
739           atomicRowLoops = null;
740       }
741     }
742 
743     extractRegion = new ExtractRegion(threadCount);
744     babinetRegion = new BabinetRegion(threadCount);
745     solventGridRegion = new SolventGridRegion(threadCount);
746     initRegion = new InitRegion(threadCount);
747     solventScaleRegion = new SolventScaleRegion(threadCount);
748     atomicScaleRegion = new AtomicScaleRegion(threadCount);
749     complexFFT3D = new Complex3DParallel(fftX, fftY, fftZ, fftTeam);
750   }
751 
752   /**
753    * Parse a solvent model string and return a SolventModel enum.
754    *
755    * @param str The solvent model String.
756    * @return The SolventModel to use.
757    */
758   public static SolventModel parseSolventModel(String str) {
759     try {
760       return SolventModel.valueOf(str.toUpperCase().replaceAll("\\s+", ""));
761     } catch (Exception e) {
762       logger.info(format(" Could not parse %s as a Solvent Model; defaulting to Polynomial.", str));
763       return SolventModel.POLYNOMIAL;
764     }
765   }
766 
767   /**
768    * compute inverse FFT to determine atomic gradients
769    *
770    * @param hklData structure factors to apply inverse FFT
771    * @param freer array of free r flags corresponding to hkldata
772    * @param flag Rfree flag value
773    * @param refinementMode {@link RefinementMinimize.RefinementMode refinement mode}
774    * @see RefinementMinimize.RefinementMode
775    * @see DiffractionRefinementData
776    */
777   public void computeAtomicGradients(
778       double[][] hklData, int[] freer, int flag, RefinementMode refinementMode) {
779     computeAtomicGradients(hklData, freer, flag, refinementMode, false);
780   }
781 
782   /**
783    * parallelized computation of bulk solvent structure factors
784    *
785    * @param hklData structure factor list to fill in
786    * @see DiffractionRefinementData
787    */
788   public void computeSolventDensity(double[][] hklData) {
789     computeSolventDensity(hklData, false);
790   }
791 
792   /**
793    * densityNorm
794    *
795    * @param data an array of double.
796    * @param meansd an array of double.
797    * @param norm a boolean.
798    */
799   public void densityNorm(double[] data, double[] meansd, boolean norm) {
800     double mean, sd;
801 
802     mean = sd = 0.0;
803     int n = 0;
804     for (int k = 0; k < fftZ; k++) {
805       for (int j = 0; j < fftY; j++) {
806         for (int i = 0; i < fftX; i++) {
807           int index = iComplex3D(i, j, k, fftX, fftY);
808           n++;
809           mean += (data[index] - mean) / n;
810         }
811       }
812     }
813 
814     n = 0;
815     for (int k = 0; k < fftZ; k++) {
816       for (int j = 0; j < fftY; j++) {
817         for (int i = 0; i < fftX; i++) {
818           int index = iComplex3D(i, j, k, fftX, fftY);
819           sd += pow(data[index] - mean, 2);
820           n++;
821         }
822       }
823     }
824     sd = sqrt(sd / n);
825 
826     if (meansd != null) {
827       meansd[0] = mean;
828       meansd[1] = sd;
829     }
830 
831     if (norm) {
832       double isd = 1.0 / sd;
833       for (int k = 0; k < fftZ; k++) {
834         for (int j = 0; j < fftY; j++) {
835           for (int i = 0; i < fftX; i++) {
836             int index = iComplex3D(i, j, k, fftX, fftY);
837             data[index] = (data[index] - mean) * isd;
838           }
839         }
840       }
841     }
842   }
843 
844   /**
845    * Getter for the field <code>densityGrid</code>.
846    *
847    * @return the densityGrid
848    */
849   public double[] getDensityGrid() {
850     return densityGrid;
851   }
852 
853   /**
854    * return dataset weight
855    *
856    * @return the weight wA
857    */
858   public double getWeight() {
859     return weight;
860   }
861 
862   /**
863    * set the dataset weight
864    *
865    * @param weight desired weight wA
866    */
867   public void setWeight(double weight) {
868     this.weight = weight;
869   }
870 
871   /**
872    * getXDim
873    *
874    * @return a double.
875    */
876   public double getXDim() {
877     return fftX;
878   }
879 
880   /**
881    * getYDim
882    *
883    * @return a double.
884    */
885   public double getYDim() {
886     return fftY;
887   }
888 
889   /**
890    * getZDim
891    *
892    * @return a double.
893    */
894   public double getZDim() {
895     return fftZ;
896   }
897 
898   /**
899    * Set atomic coordinates.
900    *
901    * @param coords new coordinate positions (3 params per atom)
902    */
903   public void setCoordinates(double[] coords) {
904     assert (coords != null);
905     List<SymOp> symops = crystal.spaceGroup.symOps;
906     double[] xyz = new double[3];
907     double[] symXYZ = new double[3];
908     int index = 0;
909     for (int i = 0; i < nAtoms; i++) {
910       if (!atoms[i].isActive()) {
911         // Coordinates for inactive atoms do not change.
912         continue;
913       }
914       xyz[0] = coords[index++];
915       xyz[1] = coords[index++];
916       xyz[2] = coords[index++];
917       for (int j = 0; j < nSymm; j++) {
918         crystal.applySymOp(xyz, symXYZ, symops.get(j));
919         coordinates[j][0][i] = symXYZ[0];
920         coordinates[j][1][i] = symXYZ[1];
921         coordinates[j][2][i] = symXYZ[2];
922       }
923     }
924   }
925 
926   /**
927    * Setter for the field <code>nativeEnvironmentApproximation</code>.
928    *
929    * @param nativeEnvironmentApproximation a boolean.
930    */
931   void setNativeEnvironmentApproximation(boolean nativeEnvironmentApproximation) {
932     this.nativeEnvironmentApproximation = nativeEnvironmentApproximation;
933   }
934 
935   /**
936    * Getter for the field <code>solventGrid</code>.
937    *
938    * @return the solventGrid
939    */
940   double[] getSolventGrid() {
941     return solventGrid;
942   }
943 
944   /**
945    * Setter for the field <code>lambdaTerm</code>.
946    *
947    * @param lambdaTerm a boolean.
948    */
949   void setLambdaTerm(boolean lambdaTerm) {
950     this.lambdaTerm = lambdaTerm;
951   }
952 
953   /**
954    * set the bulk solvent parameters
955    *
956    * @param a added atom width (binary, polynomial only)
957    * @param b falloff of the switch (Gaussian, polynomial only)
958    */
959   void setSolventAB(double a, double b) {
960     double vdwr;
961     this.solventA = a;
962     this.solventB = b;
963     if (!solvent) {
964       return;
965     }
966     switch (solventModel) {
967       case BINARY:
968         for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
969           for (int i = 0; i < nAtoms; i++) {
970             vdwr = atoms[i].getVDWType().radius * 0.5;
971             solventFormFactors[iSymm][i] = new SolventBinaryFormFactor(atoms[i], vdwr + solventA);
972           }
973         }
974         break;
975       case GAUSSIAN:
976         for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
977           for (int i = 0; i < nAtoms; i++) {
978             vdwr = atoms[i].getVDWType().radius * 0.5;
979             solventFormFactors[iSymm][i] = new SolventGaussFormFactor(atoms[i], vdwr * solventB);
980           }
981         }
982         break;
983       case POLYNOMIAL:
984         for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
985           for (int i = 0; i < nAtoms; i++) {
986             vdwr = atoms[i].getVDWType().radius * 0.5;
987             solventFormFactors[iSymm][i] =
988                 new SolventPolyFormFactor(atoms[i], vdwr + solventA, solventB);
989           }
990         }
991         break;
992     }
993   }
994 
995   /**
996    * should the structure factor computation use 3 Gaussians or 6 for atoms?
997    *
998    * @param useThreeGaussians if true, use 3 Gaussian
999    */
1000   void setUse3G(boolean useThreeGaussians) {
1001     this.useThreeGaussians = useThreeGaussians;
1002     if (solvent || neutron) {
1003       return;
1004     }
1005     for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
1006       for (int i = 0; i < nAtoms; i++) {
1007         atomFormFactors[iSymm][i] = new XRayFormFactor(atoms[i], useThreeGaussians, bAdd);
1008       }
1009     }
1010   }
1011 
1012   /**
1013    * offset X coordinates (mostly for finite difference checks)
1014    *
1015    * @param n {@link ffx.potential.bonded.Atom} to apply delta to
1016    * @param delta amount to shift atom by
1017    */
1018   void deltaX(int n, double delta) {
1019     List<SymOp> symops = crystal.spaceGroup.symOps;
1020     double[] xyz = new double[3];
1021     double[] symxyz = new double[3];
1022     atoms[n].getXYZ(xyz);
1023     xyz[0] += delta;
1024     for (int i = 0; i < nSymm; i++) {
1025       crystal.applySymOp(xyz, symxyz, symops.get(i));
1026       coordinates[i][0][n] = symxyz[0];
1027     }
1028   }
1029 
1030   /**
1031    * offset Y coordinates (mostly for finite difference checks)
1032    *
1033    * @param n {@link ffx.potential.bonded.Atom} to apply delta to
1034    * @param delta amount to shift atom by
1035    */
1036   void deltaY(int n, double delta) {
1037     List<SymOp> symops = crystal.spaceGroup.symOps;
1038     double[] xyz = new double[3];
1039     double[] symxyz = new double[3];
1040     atoms[n].getXYZ(xyz);
1041     xyz[1] += delta;
1042     for (int i = 0; i < nSymm; i++) {
1043       crystal.applySymOp(xyz, symxyz, symops.get(i));
1044       coordinates[i][1][n] = symxyz[1];
1045     }
1046   }
1047 
1048   /**
1049    * offset Z coordinates (mostly for finite difference checks)
1050    *
1051    * @param n {@link ffx.potential.bonded.Atom} to apply delta to
1052    * @param delta amount to shift atom by
1053    */
1054   void deltaZ(int n, double delta) {
1055     List<SymOp> symops = crystal.spaceGroup.symOps;
1056     double[] xyz = new double[3];
1057     double[] symxyz = new double[3];
1058     atoms[n].getXYZ(xyz);
1059     xyz[2] += delta;
1060     for (int i = 0; i < nSymm; i++) {
1061       crystal.applySymOp(xyz, symxyz, symops.get(i));
1062       coordinates[i][2][n] = symxyz[2];
1063     }
1064   }
1065 
1066   /**
1067    * parallelized computation of structure factors
1068    *
1069    * @param hklData structure factor list to fill in
1070    * @see DiffractionRefinementData
1071    */
1072   void computeDensity(double[][] hklData) {
1073     computeDensity(hklData, false);
1074   }
1075 
1076   /**
1077    * parallelized computation of structure factors
1078    *
1079    * @param hklData structure factor list to fill in
1080    * @param print if true, print information on timings during the calculation
1081    * @see DiffractionRefinementData
1082    */
1083   void computeDensity(double[][] hklData, boolean print) {
1084     if (solvent) {
1085       if (solventModel != SolventModel.NONE) {
1086         computeSolventDensity(hklData, print);
1087       }
1088     } else {
1089       computeAtomicDensity(hklData, print);
1090     }
1091   }
1092 
1093   /**
1094    * compute inverse FFT to determine atomic gradients
1095    *
1096    * @param hklData structure factors to apply inverse FFT
1097    * @param freer array of free r flags corresponding to hkldata
1098    * @param flag Rfree flag value
1099    * @param refinementMode {@link RefinementMinimize.RefinementMode refinement mode}
1100    * @param print if true, print information on timings during the calculation
1101    * @see RefinementMinimize.RefinementMode
1102    * @see DiffractionRefinementData
1103    */
1104   void computeAtomicGradients(
1105       double[][] hklData, int[] freer, int flag, RefinementMode refinementMode, boolean print) {
1106 
1107     if (solvent && solventModel == SolventModel.NONE) {
1108       return;
1109     }
1110 
1111     // zero out the density
1112     for (int i = 0; i < complexFFT3DSpace; i++) {
1113       densityGrid[i] = 0.0;
1114     }
1115 
1116     int nfree = 0;
1117     StringBuilder sb = new StringBuilder();
1118     long symtime = -System.nanoTime();
1119     int nsym = crystal.spaceGroup.symOps.size();
1120     // int nsym = 1;
1121     List<SymOp> symops = crystal.spaceGroup.symOps;
1122     ComplexNumber c = new ComplexNumber();
1123     ComplexNumber cj = new ComplexNumber();
1124     HKL ij = new HKL();
1125     for (HKL ih : reflectionList.hkllist) {
1126       double[] fc = hklData[ih.getIndex()];
1127       if (Double.isNaN(fc[0])) {
1128         continue;
1129       }
1130 
1131       // cross validation check!!!
1132       if (freer != null) {
1133         if (freer[ih.getIndex()] == flag) {
1134           nfree++;
1135           continue;
1136         }
1137       }
1138 
1139       c.re(fc[0]);
1140       c.im(fc[1]);
1141       // scale
1142       c.timesIP(2.0 / fftScale);
1143 
1144       // apply symmetry
1145       for (int j = 0; j < nsym; j++) {
1146         cj.copy(c);
1147         SymOp symOp = symops.get(j);
1148         applyTransSymRot(ih, ij, symOp);
1149         double shift = symOp.symPhaseShift(ih);
1150 
1151         int h = mod(ij.getH(), fftX);
1152         int k = mod(ij.getK(), fftY);
1153         int l = mod(ij.getL(), fftZ);
1154 
1155         if (h < halfFFTX + 1) {
1156           final int ii = iComplex3D(h, k, l, fftX, fftY);
1157           cj.phaseShiftIP(shift);
1158           densityGrid[ii] += cj.re();
1159           // Added parentheses below on June 6, 2022.
1160           // TODO: Should there be a "minus"?
1161           densityGrid[ii + 1] += (-cj.im());
1162         } else {
1163           h = (fftX - h) % fftX;
1164           k = (fftY - k) % fftY;
1165           l = (fftZ - l) % fftZ;
1166           final int ii = iComplex3D(h, k, l, fftX, fftY);
1167           cj.phaseShiftIP(shift);
1168           densityGrid[ii] += cj.re();
1169           densityGrid[ii + 1] += cj.im();
1170         }
1171       }
1172     }
1173     symtime += System.nanoTime();
1174 
1175     long startTime = System.nanoTime();
1176     complexFFT3D.ifft(densityGrid);
1177     long fftTime = System.nanoTime() - startTime;
1178 
1179     /*
1180     CCP4MapWriter mapout = new CCP4MapWriter(fftX, fftY, fftZ, crystal, "/tmp/foo.map");
1181     mapout.write(densityGrid);
1182     */
1183     startTime = System.nanoTime();
1184     long permanentDensityTime = 0;
1185     try {
1186       if (solvent) {
1187         solventGradientRegion.setRefinementMode(refinementMode);
1188         parallelTeam.execute(solventGradientRegion);
1189       } else {
1190         for (int i = 0; i < nAtoms; i++) {
1191           optAtomicGradientWeight.set(i, 0);
1192         }
1193         atomicGradientSchedule.updateWeights(previousOptAtomicGradientWeight);
1194         atomicGradientRegion.setRefinementMode(refinementMode);
1195         parallelTeam.execute(atomicGradientRegion);
1196         for (int i = 0; i < nAtoms; i++) {
1197           previousOptAtomicGradientWeight[i] = optAtomicGradientWeight.get(i);
1198         }
1199       }
1200       permanentDensityTime = System.nanoTime() - startTime;
1201     } catch (Exception e) {
1202       String message = "Exception computing atomic gradients.";
1203       logger.log(Level.SEVERE, message, e);
1204     }
1205 
1206     if (solvent) {
1207       sb.append(format(" Solvent symmetry: %8.3f\n", symtime * toSeconds));
1208       sb.append(format(" Solvent inverse FFT: %8.3f\n", fftTime * toSeconds));
1209       sb.append(format(" Grid solvent gradients: %8.3f\n", permanentDensityTime * toSeconds));
1210       sb.append(format(" %d reflections ignored (cross validation set)\n", nfree));
1211     } else {
1212       sb.append(format(" Atomic symmetry: %8.3f\n", symtime * toSeconds));
1213       sb.append(format(" Atomic inverse FFT: %8.3f\n", fftTime * toSeconds));
1214       sb.append(format(" Grid atomic gradients: %8.3f\n", permanentDensityTime * toSeconds));
1215       sb.append(format(" %d reflections ignored (cross validation set)\n", nfree));
1216     }
1217 
1218     if (logger.isLoggable(Level.INFO) && print) {
1219       logger.info(sb.toString());
1220     }
1221   }
1222 
1223   /**
1224    * parallelized computation of structure factors (identical to compuateDensity)
1225    *
1226    * @param hklData structure factor list to fill in
1227    * @see DiffractionRefinementData
1228    */
1229   void computeAtomicDensity(double[][] hklData) {
1230     computeAtomicDensity(hklData, false);
1231   }
1232 
1233   /**
1234    * parallelized computation of structure factors (identical to computeDensity)
1235    *
1236    * @param hklData structure factor list to fill in
1237    * @param print if true, print information on timings during the calculation
1238    * @see DiffractionRefinementData
1239    */
1240   private void computeAtomicDensity(double[][] hklData, boolean print) {
1241     // Zero out reflection data.
1242     long initTime = -System.nanoTime();
1243     try {
1244       initRegion.setHKL(hklData);
1245       parallelTeam.execute(initRegion);
1246     } catch (Exception e) {
1247       String message = "Fatal exception zeroing out reflection data.";
1248       logger.log(Level.SEVERE, message, e);
1249     }
1250     initTime += System.nanoTime();
1251 
1252     // Assign atomic electron density to FFT grid.
1253     long atomicGridTime = -System.nanoTime();
1254     try {
1255       switch (gridMethod) {
1256         case SPATIAL:
1257           atomicDensityRegion.assignAtomsToCells();
1258           parallelTeam.execute(atomicDensityRegion);
1259           break;
1260 
1261         case ROW:
1262           for (int i = 0; i < fftZ * fftY; i++) {
1263             optRowWeightAtomic.set(i, 0);
1264           }
1265           atomicRowSchedule.updateWeights(previousOptRowWeightAtomic);
1266           parallelTeam.execute(atomicRowRegion);
1267           for (int i = 0; i < fftZ * fftY; i++) {
1268             previousOptRowWeightAtomic[i] = optRowWeightAtomic.get(i);
1269           }
1270           break;
1271         case SLICE:
1272         default:
1273           for (int i = 0; i < fftZ; i++) {
1274             optSliceWeightAtomic.set(i, 0);
1275           }
1276           atomicSliceSchedule.updateWeights(previousOptSliceWeightAtomic);
1277           parallelTeam.execute(atomicSliceRegion);
1278           for (int i = 0; i < fftZ; i++) {
1279             previousOptSliceWeightAtomic[i] = optSliceWeightAtomic.get(i);
1280           }
1281           break;
1282       }
1283     } catch (Exception e) {
1284       String message = "Fatal exception evaluating atomic electron density.";
1285       logger.log(Level.SEVERE, message, e);
1286     }
1287     atomicGridTime += System.nanoTime();
1288 
1289     // Compute model structure factors via an FFT of the electron density.
1290     long startTime = System.nanoTime();
1291     complexFFT3D.fft(densityGrid);
1292     long fftTime = System.nanoTime() - startTime;
1293 
1294     // Extract and scale structure factors.
1295     long symTime = -System.nanoTime();
1296     try {
1297       extractRegion.setHKL(hklData);
1298       parallelTeam.execute(extractRegion);
1299       double scale = (fftScale) / (fftX * fftY * fftZ);
1300       atomicScaleRegion.setHKL(hklData);
1301       atomicScaleRegion.setScale(scale);
1302       parallelTeam.execute(atomicScaleRegion);
1303     } catch (Exception e) {
1304       String message = "Fatal exception extracting structure factors.";
1305       logger.log(Level.SEVERE, message, e);
1306     }
1307     symTime += System.nanoTime();
1308 
1309     if (logger.isLoggable(Level.INFO) && print) {
1310       StringBuilder sb = new StringBuilder();
1311       sb.append(format("\n Fc Initialization:         %8.4f\n", initTime * toSeconds));
1312       sb.append(format(" Atomic grid density:       %8.4f\n", atomicGridTime * toSeconds));
1313       sb.append(format(" Atomic FFT:                %8.4f\n", fftTime * toSeconds));
1314       sb.append(format(" Atomic symmetry & scaling: %8.4f\n", symTime * toSeconds));
1315       logger.info(sb.toString());
1316 
1317       StringBuilder ASB = new StringBuilder();
1318       switch (gridMethod) {
1319         case ROW:
1320           double atomicRowTotal = 0;
1321           double atomicRowMin = atomicRowLoops[0].getThreadTime();
1322           double atomicRowMax = 0;
1323           double atomicRowWeightTotal = 0;
1324 
1325           for (int i = 0; i < threadCount; i++) {
1326             atomicRowTotal += atomicRowLoops[i].getThreadTime();
1327             atomicRowMax = max(atomicRowLoops[i].getThreadTime(), atomicRowMax);
1328             atomicRowMin = min(atomicRowLoops[i].getThreadTime(), atomicRowMin);
1329             atomicRowWeightTotal += atomicRowLoops[i].getThreadWeight();
1330           }
1331 
1332           // Atomic timing and balance analysis
1333           ASB.append(
1334               format(
1335                   "\n RowLoop (Atomic): %7.5f (sec)                 | Total Weight: %7.0f\n",
1336                   atomicGridTime * toSeconds, atomicRowWeightTotal));
1337           ASB.append(
1338               " Thread     LoopTime    Balance(%)  Normalized   |      Rows       Weight    Balance(%)  Normalized\n");
1339 
1340           // check for no weights issued, then set to 1 so a 0 is printed instead of NaN
1341           if (atomicRowWeightTotal == 0) {
1342             atomicRowWeightTotal = 1;
1343           }
1344 
1345           for (int i = 0; i < threadCount; i++) {
1346             ASB.append(
1347                 format(
1348                     "     %3d     %7.5f     %7.1f     %7.1f     |   %7d     %7d     %7.1f     %7.1f\n",
1349                     i, atomicRowLoops[i].getThreadTime() * toSeconds,
1350                     ((atomicRowLoops[i].getThreadTime()) / (atomicRowTotal)) * 100.00,
1351                     ((atomicRowLoops[i].getThreadTime()) / (atomicRowTotal)) * (100.00
1352                         * threadCount),
1353                     atomicRowLoops[i].getNumberofSlices(),
1354                     atomicRowLoops[i].getThreadWeight(),
1355                     100.00 * (atomicRowLoops[i].getThreadWeight() / atomicRowWeightTotal),
1356                     (100.00 * threadCount) * (atomicRowLoops[i].getThreadWeight()
1357                         / atomicRowWeightTotal)));
1358           }
1359           ASB.append(format("    Min      %7.5f\n", atomicRowMin * toSeconds));
1360           ASB.append(format("    Max      %7.5f\n", atomicRowMax * toSeconds));
1361           ASB.append(format("    Delta    %7.5f\n", (atomicRowMax - atomicRowMin) * toSeconds));
1362           logger.info(ASB.toString());
1363 
1364           break;
1365         case SPATIAL:
1366           break;
1367         case SLICE:
1368           double atomicSliceTotal = 0;
1369           double atomicSliceMin = atomicSliceLoops[0].getThreadTime();
1370           double atomicSliceMax = 0;
1371           double atomicSliceWeightTotal = 0;
1372 
1373           for (int i = 0; i < threadCount; i++) {
1374             atomicSliceTotal += atomicSliceLoops[i].getThreadTime();
1375             atomicSliceMax = max(atomicSliceLoops[i].getThreadTime(), atomicSliceMax);
1376             atomicSliceMin = min(atomicSliceLoops[i].getThreadTime(), atomicSliceMin);
1377             atomicSliceWeightTotal += atomicSliceLoops[i].getThreadWeight();
1378           }
1379 
1380           // Atomic timing and balance analysis
1381           ASB.append(
1382               format(
1383                   "\n SliceLoop (Atomic): %7.5f (sec)               | Total Weight: %7.0f\n",
1384                   atomicGridTime * toSeconds, atomicSliceWeightTotal));
1385           ASB.append(
1386               " Thread     LoopTime    Balance(%)  Normalized   |      Slices    Weight    Balance(%)  Normalized\n");
1387 
1388           // check for no weights issued, then set to 1 so a 0 is printed instead of NaN
1389           if (atomicSliceWeightTotal == 0) {
1390             atomicSliceWeightTotal = 1;
1391           }
1392 
1393           for (int i = 0; i < threadCount; i++) {
1394             ASB.append(
1395                 format(
1396                     "     %3d     %7.5f     %7.1f     %7.1f     |   %7d     %7d     %7.1f     %7.1f\n",
1397                     i,
1398                     atomicSliceLoops[i].getThreadTime() * toSeconds,
1399                     ((atomicSliceLoops[i].getThreadTime()) / (atomicSliceTotal)) * 100.00,
1400                     ((atomicSliceLoops[i].getThreadTime()) / (atomicSliceTotal)) * (100.00
1401                         * threadCount),
1402                     atomicSliceLoops[i].getNumberofSlices(),
1403                     atomicSliceLoops[i].getThreadWeight(),
1404                     100.00 * (atomicSliceLoops[i].getThreadWeight() / atomicSliceWeightTotal),
1405                     (100.00 * threadCount) * (atomicSliceLoops[i].getThreadWeight()
1406                         / atomicSliceWeightTotal)));
1407           }
1408           ASB.append(format("    Min      %7.5f\n", atomicSliceMin * toSeconds));
1409           ASB.append(format("    Max      %7.5f\n", atomicSliceMax * toSeconds));
1410           ASB.append(format("    Delta    %7.5f\n", (atomicSliceMax - atomicSliceMin) * toSeconds));
1411           logger.info(ASB.toString());
1412           break;
1413         default:
1414       }
1415     }
1416   }
1417 
1418   /**
1419    * parallelized computation of bulk solvent structure factors
1420    *
1421    * @param hklData structure factor list to fill in
1422    * @param print if true, print information on timings during the calculation
1423    * @see DiffractionRefinementData
1424    */
1425   private void computeSolventDensity(double[][] hklData, boolean print) {
1426     // Zero out the reflection data.
1427     long initTime = -System.nanoTime();
1428     try {
1429       initRegion.setHKL(hklData);
1430       parallelTeam.execute(initRegion);
1431     } catch (Exception e) {
1432       String message = "Fatal exception initializing structure factors.";
1433       logger.log(Level.SEVERE, message, e);
1434     }
1435     initTime += System.nanoTime();
1436 
1437     long solventGridTime = -System.nanoTime();
1438     try {
1439       switch (gridMethod) {
1440         case SPATIAL:
1441           atomicDensityRegion.assignAtomsToCells();
1442           parallelTeam.execute(atomicDensityRegion);
1443           break;
1444         case ROW:
1445           for (int i = 0; i < fftZ * fftY; i++) {
1446             optRowWeightSolvent.set(i, 0);
1447           }
1448           solventRowSchedule.updateWeights(previousOptRowWeightSolvent);
1449           parallelTeam.execute(atomicRowRegion);
1450           for (int i = 0; i < fftZ * fftY; i++) {
1451             previousOptRowWeightSolvent[i] = optRowWeightSolvent.get(i);
1452           }
1453           break;
1454         case SLICE:
1455         default:
1456           for (int i = 0; i < fftZ; i++) {
1457             optSliceWeightSolvent.set(i, 0);
1458           }
1459           solventSliceSchedule.updateWeights(previousOptSliceWeightSolvent);
1460           parallelTeam.execute(atomicSliceRegion);
1461           for (int i = 0; i < fftZ; i++) {
1462             previousOptSliceWeightSolvent[i] = optSliceWeightSolvent.get(i);
1463           }
1464           break;
1465       }
1466     } catch (Exception e) {
1467       String message = "Fatal exception evaluating solvent electron density.";
1468       logger.log(Level.SEVERE, message, e);
1469     }
1470     solventGridTime += System.nanoTime();
1471 
1472     long expTime = -System.nanoTime();
1473     if (solventModel == SolventModel.BINARY) {
1474       /*
1475        * Need to shrink mask. First, copy densityGrid to solventGrid
1476        * temporarily to store original map.
1477        */
1478       int nmap = densityGrid.length;
1479       arraycopy(densityGrid, 0, solventGrid, 0, nmap);
1480       // Logic to loop within the cutoff box.
1481       double[] xyz = {solventB, solventB, solventB};
1482       double[] uvw = new double[3];
1483       crystal.toFractionalCoordinates(xyz, uvw);
1484       final double frx = fftX * uvw[0];
1485       final int ifrx = (int) frx;
1486       final double fry = fftY * uvw[1];
1487       final int ifry = (int) fry;
1488       final double frz = fftZ * uvw[2];
1489       final int ifrz = (int) frz;
1490       for (int k = 0; k < fftZ; k++) {
1491         for (int j = 0; j < fftY; j++) {
1492           for (int i = 0; i < fftX; i++) {
1493             final int ii = iComplex3D(i, j, k, fftX, fftY);
1494             if (densityGrid[ii] == 1.0) {
1495               continue;
1496             }
1497             for (int iz = k - ifrz; iz <= k + ifrz; iz++) {
1498               int giz = mod(iz, fftZ);
1499               for (int iy = j - ifry; iy <= j + ifry; iy++) {
1500                 int giy = mod(iy, fftY);
1501                 for (int ix = i - ifrx; ix <= i + ifrx; ix++) {
1502                   int gix = mod(ix, fftX);
1503                   final int jj = iComplex3D(gix, giy, giz, fftX, fftY);
1504                   if (densityGrid[jj] == 1.0) {
1505                     solventGrid[ii] = 1.0;
1506                   }
1507                 }
1508               }
1509             }
1510           }
1511         }
1512       }
1513 
1514       // Copy the completed solvent grid back.
1515       arraycopy(solventGrid, 0, densityGrid, 0, nmap);
1516     }
1517 
1518     // Copy the grid over for derivatives.
1519     ArrayCopyRegion arrayCopyRegion = new ArrayCopyRegion();
1520     try {
1521       parallelTeam.execute(arrayCopyRegion);
1522     } catch (Exception e) {
1523       logger.info(e.toString());
1524     }
1525 
1526     // Babinet Principle.
1527     try {
1528       parallelTeam.execute(babinetRegion);
1529     } catch (Exception e) {
1530       String message = "Fatal exception Babinet Principle.";
1531       logger.log(Level.SEVERE, message, e);
1532     }
1533     expTime += System.nanoTime();
1534 
1535     long fftTime = -System.nanoTime();
1536     complexFFT3D.fft(densityGrid);
1537     fftTime += System.nanoTime();
1538 
1539     // Extract and scale structure factors.
1540     long symTime = -System.nanoTime();
1541     try {
1542       extractRegion.setHKL(hklData);
1543       parallelTeam.execute(extractRegion);
1544       final double scale = (fftScale) / (fftX * fftY * fftZ);
1545       solventScaleRegion.setHKL(hklData);
1546       solventScaleRegion.setScale(scale);
1547       parallelTeam.execute(solventScaleRegion);
1548     } catch (Exception e) {
1549       String message = "Fatal exception extracting structure factors.";
1550       logger.log(Level.SEVERE, message, e);
1551     }
1552     symTime += System.nanoTime();
1553 
1554     // Expand derivative mask so it includes nearby local density.
1555     long solventDensityTime = -System.nanoTime();
1556     try {
1557       switch (gridMethod) {
1558         case SPATIAL:
1559           solventDensityRegion.assignAtomsToCells();
1560           parallelTeam.execute(solventDensityRegion);
1561           break;
1562 
1563         case ROW:
1564           for (int i = 0; i < (fftZ * fftY); i++) {
1565             optRowWeightBulkSolvent.set(i, 0);
1566           }
1567           bulkSolventRowSchedule.updateWeights(previousOptRowWeightBulkSolvent);
1568           parallelTeam.execute(solventRowRegion);
1569           for (int i = 0; i < fftZ * fftY; i++) {
1570             previousOptRowWeightBulkSolvent[i] = optRowWeightBulkSolvent.get(i);
1571           }
1572           break;
1573         case SLICE:
1574         default:
1575           for (int i = 0; i < (fftZ); i++) {
1576             optSliceWeightBulkSolvent.set(i, 0);
1577           }
1578 
1579           bulkSolventSliceSchedule.updateWeights(previousOptSliceWeightBulkSolvent);
1580           parallelTeam.execute(solventSliceRegion);
1581           for (int i = 0; i < fftZ; i++) {
1582             previousOptSliceWeightBulkSolvent[i] = optSliceWeightBulkSolvent.get(i);
1583           }
1584           break;
1585       }
1586     } catch (Exception e) {
1587       String message = "Fatal exception evaluating solvent electron density.";
1588       logger.log(Level.SEVERE, message, e);
1589     }
1590     solventDensityTime += System.nanoTime();
1591 
1592     long solventExpTime = 0;
1593     if (solventModel == SolventModel.GAUSSIAN) {
1594       solventExpTime = -System.nanoTime();
1595       try {
1596         parallelTeam.execute(solventGridRegion);
1597       } catch (Exception e) {
1598         String message = "Fatal exception solvent grid region.";
1599         logger.log(Level.SEVERE, message, e);
1600       }
1601       solventExpTime += System.nanoTime();
1602     }
1603 
1604     if (logger.isLoggable(Level.INFO) && print) {
1605       StringBuilder sb = new StringBuilder();
1606       sb.append(format("\n Fc Initialization:                     %8.4f\n", initTime * toSeconds));
1607       sb.append(
1608           format(" Solvent grid density:                  %8.4f\n", solventGridTime * toSeconds));
1609       sb.append(format(" Bulk solvent exponentiation:           %8.4f\n", expTime * toSeconds));
1610       sb.append(format(" Solvent FFT:                           %8.4f\n", fftTime * toSeconds));
1611       sb.append(format(" Solvent symmetry & scaling:            %8.4f\n", symTime * toSeconds));
1612       sb.append(
1613           format(
1614               " Solvent grid expansion:                %8.4f\n", solventDensityTime * toSeconds));
1615       if (solventModel == SolventModel.GAUSSIAN) {
1616         sb.append(format(" Gaussian grid exponentiation: %8.4f\n", solventExpTime * toSeconds));
1617       }
1618       logger.info(sb.toString());
1619 
1620       StringBuilder SSB = new StringBuilder();
1621       StringBuilder BSSB = new StringBuilder();
1622       switch (gridMethod) {
1623         case ROW:
1624           double solventRowTotal = 0;
1625           double solventRowMin = solventRowLoops[0].getThreadTime();
1626           double solventRowMax = 0;
1627           double solventRowWeightTotal = 0;
1628 
1629           double bulkSolventRowTotal = 0;
1630           double bulkSolventRowMin = bulkSolventRowLoops[0].getTimePerThread();
1631           double bulkSolventRowMax = 0;
1632           double bulkSolventRowWeightTotal = 0;
1633 
1634           for (int i = 0; i < threadCount; i++) {
1635             solventRowTotal += solventRowLoops[i].getThreadTime();
1636             solventRowMax = max(solventRowLoops[i].getThreadTime(), solventRowMax);
1637             solventRowMin = min(solventRowLoops[i].getThreadTime(), solventRowMin);
1638             solventRowWeightTotal += solventRowLoops[i].getThreadWeight();
1639 
1640             bulkSolventRowTotal += bulkSolventRowLoops[i].getTimePerThread();
1641             bulkSolventRowMax = max(bulkSolventRowLoops[i].getTimePerThread(), bulkSolventRowMax);
1642             bulkSolventRowMin = min(bulkSolventRowLoops[i].getTimePerThread(), bulkSolventRowMin);
1643             bulkSolventRowWeightTotal += bulkSolventRowLoops[i].getWeightPerThread()[i];
1644           }
1645 
1646           // Solvent timing and balance analysis
1647           SSB.append(
1648               format(
1649                   "\n RowLoop (Solvent): %7.5f (sec)                | Total Weight: %7.0f\n",
1650                   solventGridTime * toSeconds, solventRowWeightTotal));
1651           SSB.append(
1652               " Thread     LoopTime    Balance(%)  Normalized   |      Rows       Weight     Balance(%)  Normalized\n");
1653 
1654           // Check for no weights issued, then set to 1 so a 0 is printed instead of NaN
1655           if (solventRowWeightTotal == 0) {
1656             solventRowWeightTotal = 1;
1657           }
1658 
1659           for (int i = 0; i < threadCount; i++) {
1660             SSB.append(
1661                 format(
1662                     "     %3d     %7.5f     %7.1f     %7.1f     |   %7d     %7d     %7.1f     %7.1f\n",
1663                     i, solventRowLoops[i].getThreadTime() * toSeconds,
1664                     ((solventRowLoops[i].getThreadTime()) / (solventRowTotal)) * 100.00,
1665                     ((solventRowLoops[i].getThreadTime()) / (solventRowTotal))
1666                         * (100.00 * threadCount),
1667                     solventRowLoops[i].getNumberofSlices(),
1668                     solventRowLoops[i].getThreadWeight(),
1669                     100.00 * (solventRowLoops[i].getThreadWeight() / solventRowWeightTotal),
1670                     (100.00 * threadCount)
1671                         * (solventRowLoops[i].getThreadWeight() / solventRowWeightTotal)));
1672           }
1673           SSB.append(format("    Min      %7.5f\n", solventRowMin * toSeconds));
1674           SSB.append(format("    Max      %7.5f\n", solventRowMax * toSeconds));
1675           SSB.append(format("    Delta    %7.5f\n", (solventRowMax - solventRowMin) * toSeconds));
1676           logger.info(SSB.toString());
1677 
1678           // Bulk solvent timing and balance analysis
1679           BSSB.append(
1680               format(
1681                   "\n RowLoop (Bulk Solvent): %7.5f (sec)           | Total Weight: %7.0f\n",
1682                   bulkSolventRowTotal * toSeconds, bulkSolventRowWeightTotal));
1683           BSSB.append(
1684               " Thread     LoopTime    Balance(%)  Normalized   |      Rows       Weight     Balance(%)  Normalized\n");
1685 
1686           // check for no weights issued, then set to 1 so a 0 is printed instead of NaN
1687           if (bulkSolventRowWeightTotal == 0) {
1688             bulkSolventRowWeightTotal = 1;
1689           }
1690 
1691           for (int i = 0; i < threadCount; i++) {
1692             BSSB.append(
1693                 format(
1694                     "     %3d     %7.5f     %7.1f     %7.1f     |   %7d     %7d     %7.1f     %7.1f\n",
1695                     i, bulkSolventRowLoops[i].getTimePerThread() * toSeconds,
1696                     ((bulkSolventRowLoops[i].getTimePerThread()) / (bulkSolventRowTotal)) * 100.00,
1697                     ((bulkSolventRowLoops[i].getTimePerThread()) / (bulkSolventRowTotal))
1698                         * (100.00 * threadCount),
1699                     bulkSolventRowLoops[i].getBoundsPerThread()[i],
1700                     bulkSolventRowLoops[i].getWeightPerThread()[i],
1701                     100.00
1702                         * (bulkSolventRowLoops[i].getWeightPerThread()[i]
1703                         / bulkSolventRowWeightTotal),
1704                     (100.00 * threadCount)
1705                         * (bulkSolventRowLoops[i].getWeightPerThread()[i]
1706                         / bulkSolventRowWeightTotal)));
1707           }
1708           BSSB.append(format("    Min      %7.5f\n", bulkSolventRowMin * toSeconds));
1709           BSSB.append(format("    Max      %7.5f\n", bulkSolventRowMax * toSeconds));
1710           BSSB.append(
1711               format("    Delta    %7.5f\n", (bulkSolventRowMax - bulkSolventRowMin) * toSeconds));
1712           logger.info(BSSB.toString());
1713           break;
1714         case SPATIAL:
1715           break;
1716         case SLICE:
1717           double solventSliceTotal = 0;
1718           double solventSliceMin = solventSliceLoops[0].getThreadTime();
1719           double solventSliceMax = 0;
1720           double solventSliceWeightTotal = 0;
1721 
1722           double bulkSolventSliceTotal = 0;
1723           double bulkSolventSliceMin = bulkSolventSliceLoops[0].getTimePerThread();
1724           double bulkSolventSliceMax = 0;
1725           double bulkSolventSliceWeightTotal = 0;
1726 
1727           for (int i = 0; i < threadCount; i++) {
1728             solventSliceTotal += solventSliceLoops[i].getThreadTime();
1729             solventSliceMax = max(solventSliceLoops[i].getThreadTime(), solventSliceMax);
1730             solventSliceMin = min(solventSliceLoops[i].getThreadTime(), solventSliceMin);
1731             solventSliceWeightTotal += solventSliceLoops[i].getThreadWeight();
1732 
1733             bulkSolventSliceTotal += bulkSolventSliceLoops[i].getTimePerThread();
1734             bulkSolventSliceMax =
1735                 max(bulkSolventSliceLoops[i].getTimePerThread(), bulkSolventSliceMax);
1736             bulkSolventSliceMin =
1737                 min(bulkSolventSliceLoops[i].getTimePerThread(), bulkSolventSliceMin);
1738             bulkSolventSliceWeightTotal += bulkSolventSliceLoops[i].getWeightPerThread()[i];
1739           }
1740 
1741           // Solvent timing and balance analysis
1742           SSB.append(
1743               format(
1744                   "\n SliceLoop (Solvent): %7.5f (sec)               | Total Weight: %7.0f\n",
1745                   solventGridTime * toSeconds, solventSliceWeightTotal));
1746           SSB.append(
1747               " Thread     LoopTime    Balance(%)  Normalized   |      Slices    Weight     Balance(%)  Normalized\n");
1748 
1749           // Check for no weights issued, then set to 1 so a 0 is printed instead of NaN
1750           if (solventSliceWeightTotal == 0) {
1751             solventSliceWeightTotal = 1;
1752           }
1753 
1754           for (int i = 0; i < threadCount; i++) {
1755             SSB.append(
1756                 format(
1757                     "     %3d     %7.5f     %7.1f     %7.1f     |   %7d     %7d     %7.1f     %7.1f\n",
1758                     i, solventSliceLoops[i].getThreadTime() * toSeconds,
1759                     ((solventSliceLoops[i].getThreadTime()) / (solventSliceTotal)) * 100.00,
1760                     ((solventSliceLoops[i].getThreadTime()) / (solventSliceTotal))
1761                         * (100.00 * threadCount),
1762                     solventSliceLoops[i].getNumberofSlices(),
1763                     solventSliceLoops[i].getThreadWeight(),
1764                     100.00 * (solventSliceLoops[i].getThreadWeight() / solventSliceWeightTotal),
1765                     (100.00 * threadCount)
1766                         * (solventSliceLoops[i].getThreadWeight() / solventSliceWeightTotal)));
1767           }
1768           SSB.append(format("    Min      %7.5f\n", solventSliceMin * toSeconds));
1769           SSB.append(format("    Max      %7.5f\n", solventSliceMax * toSeconds));
1770           SSB.append(
1771               format("    Delta    %7.5f\n", (solventSliceMax - solventSliceMin) * toSeconds));
1772           logger.info(SSB.toString());
1773 
1774           // Bulk solvent timing and balance analysis
1775           BSSB.append(
1776               format(
1777                   "\n SliceLoop (Bulk Solvent): %7.5f (sec)          | Total Weight: %7.0f\n",
1778                   bulkSolventSliceTotal * toSeconds, bulkSolventSliceWeightTotal));
1779           BSSB.append(
1780               " Thread     LoopTime    Balance(%)  Normalized   |      Slices     Weight     Balance(%)  Normalized\n");
1781 
1782           // Check for no weights issued, then set to 1 so a 0 is printed instead of NaN
1783           if (bulkSolventSliceWeightTotal == 0) {
1784             bulkSolventSliceWeightTotal = 1;
1785           }
1786 
1787           for (int i = 0; i < threadCount; i++) {
1788             BSSB.append(
1789                 format(
1790                     "     %3d     %7.5f     %7.1f     %7.1f     |   %7d     %7d     %7.1f     %7.1f\n",
1791                     i, bulkSolventSliceLoops[i].getTimePerThread() * toSeconds,
1792                     ((bulkSolventSliceLoops[i].getTimePerThread()) / (bulkSolventSliceTotal))
1793                         * 100.00,
1794                     ((bulkSolventSliceLoops[i].getTimePerThread()) / (bulkSolventSliceTotal))
1795                         * (100.00 * threadCount),
1796                     bulkSolventSliceLoops[i].getBoundsPerThread()[i],
1797                     bulkSolventSliceLoops[i].getWeightPerThread()[i],
1798                     100.00
1799                         * (bulkSolventSliceLoops[i].getWeightPerThread()[i]
1800                         / bulkSolventSliceWeightTotal),
1801                     (100.00 * threadCount)
1802                         * (bulkSolventSliceLoops[i].getWeightPerThread()[i]
1803                         / bulkSolventSliceWeightTotal)));
1804           }
1805           BSSB.append(format("    Min      %7.5f\n", bulkSolventSliceMin * toSeconds));
1806           BSSB.append(format("    Max      %7.5f\n", bulkSolventSliceMax * toSeconds));
1807           BSSB.append(
1808               format(
1809                   "    Delta    %7.5f\n", (bulkSolventSliceMax - bulkSolventSliceMin) * toSeconds));
1810           logger.info(BSSB.toString());
1811           break;
1812         default:
1813       }
1814     }
1815   }
1816 
1817   /** The possible solvent model methods */
1818   public enum SolventModel {
1819 
1820     /** Do not model solvent scattering. */
1821     NONE,
1822     /** The classic binary (0, 1) model. */
1823     BINARY,
1824     /** Smooth the boundary of the classic model using Gaussians. */
1825     GAUSSIAN,
1826     /** Smooth the boundar of the classic model using a cubic polynomial switch (default). */
1827     POLYNOMIAL
1828   }
1829 
1830   public enum GridMethod {
1831     SPATIAL,
1832     SLICE,
1833     ROW
1834   }
1835 
1836   private class AtomicDensityLoop extends SpatialDensityLoop {
1837 
1838     final double[] xyz = new double[3];
1839     final double[] uvw = new double[3];
1840     final double[] xc = new double[3];
1841     final double[] xf = new double[3];
1842     final double[] grid;
1843 
1844     AtomicDensityLoop(SpatialDensityRegion region) {
1845       super(region, region.getNsymm(), region.actualCount);
1846       grid = region.getGrid();
1847     }
1848 
1849     @Override
1850     public void gridDensity(int iSymm, int n) {
1851       if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
1852         return;
1853       }
1854 
1855       if (lambdaTerm && atoms[n].applyLambda()) {
1856         return;
1857       }
1858 
1859       xyz[0] = coordinates[iSymm][0][n];
1860       xyz[1] = coordinates[iSymm][1][n];
1861       xyz[2] = coordinates[iSymm][2][n];
1862       FormFactor atomff = atomFormFactors[iSymm][n];
1863       crystal.toFractionalCoordinates(xyz, uvw);
1864       final int frad =
1865           Math.min(
1866               aRadGrid, (int) Math.floor(atoms[n].getFormFactorWidth() * fftX / crystal.a) + 1);
1867 
1868       // Logic to loop within the cutoff box.
1869       final double frx = fftX * uvw[0];
1870       final int ifrx = (int) frx;
1871       final int ifrxu = ifrx + frad;
1872 
1873       final double fry = fftY * uvw[1];
1874       final int ifry = (int) fry;
1875       final int ifryu = ifry + frad;
1876 
1877       final double frz = fftZ * uvw[2];
1878       final int ifrz = (int) frz;
1879       final int ifrzu = ifrz + frad;
1880 
1881       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
1882         int giz = mod(iz, fftZ);
1883         xf[2] = iz * ifftZ;
1884         for (int iy = ifry - frad; iy <= ifryu; iy++) {
1885           int giy = mod(iy, fftY);
1886           xf[1] = iy * ifftY;
1887           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
1888             int gix = mod(ix, fftX);
1889             xf[0] = ix * ifftX;
1890             crystal.toCartesianCoordinates(xf, xc);
1891             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
1892             grid[ii] = atomff.rho(grid[ii], 1.0, xc);
1893           }
1894         }
1895       }
1896     }
1897   }
1898 
1899   private class SolventDensityLoop extends SpatialDensityLoop {
1900 
1901     final double[] xyz = new double[3];
1902     final double[] uvw = new double[3];
1903     final double[] xc = new double[3];
1904     final double[] xf = new double[3];
1905     final double[] grid;
1906 
1907     SolventDensityLoop(SpatialDensityRegion region) {
1908       super(region, region.getNsymm(), region.actualCount);
1909       grid = region.getGrid();
1910     }
1911 
1912     @Override
1913     public void gridDensity(int iSymm, int n) {
1914       if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
1915         return;
1916       }
1917 
1918       if (lambdaTerm && atoms[n].applyLambda()) {
1919         return;
1920       }
1921 
1922       xyz[0] = coordinates[iSymm][0][n];
1923       xyz[1] = coordinates[iSymm][1][n];
1924       xyz[2] = coordinates[iSymm][2][n];
1925       FormFactor solventff = solventFormFactors[iSymm][n];
1926       crystal.toFractionalCoordinates(xyz, uvw);
1927       double vdwr = atoms[n].getVDWType().radius * 0.5;
1928       int frad = aRadGrid;
1929       switch (solventModel) {
1930         case BINARY:
1931           frad =
1932               Math.min(aRadGrid, (int) Math.floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
1933           break;
1934         case GAUSSIAN:
1935           frad =
1936               Math.min(aRadGrid, (int) Math.floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
1937           break;
1938         case POLYNOMIAL:
1939           frad =
1940               Math.min(aRadGrid, (int) Math.floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
1941           break;
1942         case NONE:
1943         default:
1944           return;
1945       }
1946 
1947       // Logic to loop within the cutoff box.
1948       final double frx = fftX * uvw[0];
1949       final int ifrx = (int) frx;
1950       final int ifrxu = ifrx + frad;
1951 
1952       final double fry = fftY * uvw[1];
1953       final int ifry = (int) fry;
1954       final int ifryu = ifry + frad;
1955 
1956       final double frz = fftZ * uvw[2];
1957       final int ifrz = (int) frz;
1958       final int ifrzu = ifrz + frad;
1959 
1960       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
1961         int giz = mod(iz, fftZ);
1962         xf[2] = iz * ifftZ;
1963         for (int iy = ifry - frad; iy <= ifryu; iy++) {
1964           int giy = mod(iy, fftY);
1965           xf[1] = iy * ifftY;
1966           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
1967             int gix = mod(ix, fftX);
1968             xf[0] = ix * ifftX;
1969             crystal.toCartesianCoordinates(xf, xc);
1970             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
1971             grid[ii] = solventff.rho(grid[ii], 1.0, xc);
1972           }
1973         }
1974       }
1975     }
1976   }
1977 
1978   private class AtomicRowLoop extends RowLoop {
1979 
1980     final double[] xyz = new double[3];
1981     final double[] uvw = new double[3];
1982     final double[] xc = new double[3];
1983     final double[] xf = new double[3];
1984     final double[] grid;
1985     final int[] optLocal;
1986     long timer;
1987     int previousUB, previousLB;
1988     int actualWeight;
1989 
1990     AtomicRowLoop(RowRegion region) {
1991       super(region.getNatoms(), region.getNsymm(), region);
1992       grid = region.getGrid();
1993       optLocal = new int[fftZ * fftY];
1994     }
1995 
1996     public void buildList(int iSymm, int iAtom, int lb, int ub) {
1997       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
1998         return;
1999       }
2000       xyz[0] = coordinates[iSymm][0][iAtom];
2001       xyz[1] = coordinates[iSymm][1][iAtom];
2002       xyz[2] = coordinates[iSymm][2][iAtom];
2003       crystal.toFractionalCoordinates(xyz, uvw);
2004       final int frad =
2005           min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2006 
2007       final double frz = fftZ * uvw[2];
2008       final int ifrz = (int) frz;
2009       final int ifrzu = ifrz + frad;
2010       final int ifrzl = ifrz - frad;
2011 
2012       final double fry = fftY * uvw[1];
2013       final int ifry = (int) fry;
2014       final int ifryu = ifry + frad;
2015       final int ifryl = ifry - frad;
2016 
2017       // Loop over allowed z coordinates for this Loop
2018       // Check if the current atom is close enough
2019       // If so, add to list.
2020       int buff = bufferSize;
2021 
2022       int lbZ = rowRegion.zFromRowIndex(lb);
2023       int ubZ = rowRegion.zFromRowIndex(ub);
2024 
2025       for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2026         int giz = mod(iz, fftZ);
2027         if (lbZ > giz || giz > ubZ) {
2028           continue;
2029         }
2030         int rowLB = rowRegion.rowIndexForYZ(mod(ifryl - buff, fftY), giz);
2031         int rowUB = rowRegion.rowIndexForYZ(mod(ifryu + buff, fftY), giz);
2032         if (lb >= rowLB || rowUB <= ub) {
2033           buildListA.add(iAtom);
2034           buildListS.add(iSymm);
2035           break;
2036         }
2037       }
2038     }
2039 
2040     @Override
2041     public boolean checkList(int[][][] zyAtListBuild, int buff) {
2042       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2043         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2044           if (rowRegion.select[iSymm][iAtom]) {
2045             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2046               continue;
2047             }
2048             xyz[0] = coordinates[iSymm][0][iAtom];
2049             xyz[1] = coordinates[iSymm][1][iAtom];
2050             xyz[2] = coordinates[iSymm][2][iAtom];
2051             crystal.toFractionalCoordinates(xyz, uvw);
2052             final double frz = fftZ * uvw[2];
2053             final int ifrz = (int) frz;
2054             final int previousZ = zyAtListBuild[iSymm][iAtom][0];
2055 
2056             final double fry = fftY * uvw[1];
2057             final int ifry = (int) fry;
2058             final int previousY = zyAtListBuild[iSymm][iAtom][1];
2059             if (abs(ifrz - previousZ) >= buff || abs(ifry - previousY) >= buff) {
2060               return true;
2061             }
2062           }
2063         }
2064       }
2065       return false;
2066     }
2067 
2068     @Override
2069     public void finish() {
2070       for (int i = 0; i < fftZ * fftY; i++) {
2071         optRowWeightAtomic.addAndGet(i, optLocal[i]);
2072       }
2073       timer += System.nanoTime();
2074     }
2075 
2076     @Override
2077     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2078       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2079         return;
2080       }
2081 
2082       if (lambdaTerm && atoms[iAtom].applyLambda()) {
2083         return;
2084       }
2085 
2086       int lbZ = rowRegion.zFromRowIndex(lb);
2087       int ubZ = rowRegion.zFromRowIndex(ub);
2088 
2089       xyz[0] = coordinates[iSymm][0][iAtom];
2090       xyz[1] = coordinates[iSymm][1][iAtom];
2091       xyz[2] = coordinates[iSymm][2][iAtom];
2092       FormFactor atomff = atomFormFactors[iSymm][iAtom];
2093       crystal.toFractionalCoordinates(xyz, uvw);
2094       final int frad =
2095           min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2096 
2097       // Logic to loop within the cutoff box.
2098       final double frx = fftX * uvw[0];
2099       final int ifrx = (int) frx;
2100       final int ifrxu = ifrx + frad;
2101 
2102       final double fry = fftY * uvw[1];
2103       final int ifry = (int) fry;
2104       final int ifryu = ifry + frad;
2105 
2106       final double frz = fftZ * uvw[2];
2107       final int ifrz = (int) frz;
2108       final int ifrzu = ifrz + frad;
2109 
2110       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2111         int giz = mod(iz, fftZ);
2112         if (lbZ > giz || giz > ubZ) {
2113           continue;
2114         }
2115         xf[2] = iz * ifftZ;
2116         for (int iy = ifry - frad; iy <= ifryu; iy++) {
2117           int giy = mod(iy, fftY);
2118           int rowIndex = rowRegion.rowIndexForYZ(giy, giz);
2119           if (lb > rowIndex || rowIndex > ub) {
2120             continue;
2121           }
2122           xf[1] = iy * ifftY;
2123           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2124             int gix = mod(ix, fftX);
2125             xf[0] = ix * ifftX;
2126             crystal.toCartesianCoordinates(xf, xc);
2127             optLocal[rowIndex]++;
2128             actualWeight++;
2129             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
2130             grid[ii] = atomff.rho(grid[ii], 1.0, xc);
2131           }
2132         }
2133       }
2134     }
2135 
2136     @Override
2137     public void run(int lb, int ub) throws Exception {
2138       boolean boundsChange = false;
2139       if (previousLB != lb || previousUB != ub) {
2140         boundsChange = true;
2141       }
2142       previousLB = lb;
2143       previousUB = ub;
2144       if (rebuildList || boundsChange) {
2145         buildListA = new ArrayList<>();
2146         buildListS = new ArrayList<>();
2147         for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2148           for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2149             if (rowRegion.select[iSymm][iAtom]) {
2150               buildList(iSymm, iAtom, lb, ub);
2151             }
2152           }
2153         }
2154       }
2155       for (int i = 0; i < buildListA.size(); i++) {
2156         if (rowRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2157           gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2158         }
2159       }
2160     }
2161 
2162     @Override
2163     public void saveZYValues(int[][][] zyAtListBuild) {
2164       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2165         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2166           zyAtListBuild[iSymm][iAtom][0] = 0;
2167           zyAtListBuild[iSymm][iAtom][1] = 0;
2168           if (rowRegion.select[iSymm][iAtom]) {
2169             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2170               continue;
2171             }
2172             xyz[0] = coordinates[iSymm][0][iAtom];
2173             xyz[1] = coordinates[iSymm][1][iAtom];
2174             xyz[2] = coordinates[iSymm][2][iAtom];
2175             crystal.toFractionalCoordinates(xyz, uvw);
2176             final double frz = fftZ * uvw[2];
2177             final int ifrz = (int) frz;
2178 
2179             final double fry = fftY * uvw[1];
2180             final int ifry = (int) fry;
2181             zyAtListBuild[iSymm][iAtom][0] = ifrz;
2182             zyAtListBuild[iSymm][iAtom][1] = ifry;
2183           }
2184         }
2185       }
2186     }
2187 
2188     @Override
2189     public IntegerSchedule schedule() {
2190       return atomicRowSchedule;
2191     }
2192 
2193     @Override
2194     public void start() {
2195       fill(optLocal, 0);
2196       actualWeight = 0;
2197       timer = -System.nanoTime();
2198     }
2199 
2200     double getThreadTime() {
2201       return timer;
2202     }
2203 
2204     int getNumberofSlices() {
2205       return (previousUB - previousLB + 1);
2206     }
2207 
2208     int getThreadWeight() {
2209       return actualWeight;
2210     }
2211   }
2212 
2213   private class SolventRowLoop extends RowLoop {
2214 
2215     final double[] xyz = new double[3];
2216     final double[] uvw = new double[3];
2217     final double[] xc = new double[3];
2218     final double[] xf = new double[3];
2219     final double[] grid;
2220     final int[] optLocal;
2221     long timer;
2222     int previousUB, previousLB;
2223     int actualWeight;
2224 
2225     public SolventRowLoop(RowRegion region) {
2226       super(region.getNatoms(), region.getNsymm(), region);
2227       grid = region.getGrid();
2228       optLocal = new int[fftZ * fftY];
2229     }
2230 
2231     public void buildList(int iSymm, int iAtom, int lb, int ub) {
2232       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2233         return;
2234       }
2235       xyz[0] = coordinates[iSymm][0][iAtom];
2236       xyz[1] = coordinates[iSymm][1][iAtom];
2237       xyz[2] = coordinates[iSymm][2][iAtom];
2238       crystal.toFractionalCoordinates(xyz, uvw);
2239       final int frad =
2240           min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2241 
2242       final double frz = fftZ * uvw[2];
2243       final int ifrz = (int) frz;
2244       final int ifrzu = ifrz + frad;
2245       final int ifrzl = ifrz - frad;
2246 
2247       final double fry = fftY * uvw[1];
2248       final int ifry = (int) fry;
2249       final int ifryu = ifry + frad;
2250       final int ifryl = ifry - frad;
2251 
2252       // Loop over allowed z coordinates for this Loop
2253       // Check if the current atom is close enough
2254       // If so, add to list.
2255       int buff = bufferSize;
2256 
2257       int lbZ = rowRegion.zFromRowIndex(lb);
2258       int ubZ = rowRegion.zFromRowIndex(ub);
2259 
2260       for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2261         int giz = mod(iz, fftZ);
2262         if (lbZ > giz || giz > ubZ) {
2263           continue;
2264         }
2265         int rowLB = rowRegion.rowIndexForYZ(mod(ifryl - buff, fftY), giz);
2266         int rowUB = rowRegion.rowIndexForYZ(mod(ifryu + buff, fftY), giz);
2267         if (lb >= rowLB || rowUB <= ub) {
2268           buildListA.add(iAtom);
2269           buildListS.add(iSymm);
2270           break;
2271         }
2272       }
2273     }
2274 
2275     @Override
2276     public boolean checkList(int[][][] zyAtListBuild, int buff) {
2277       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2278         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2279           if (rowRegion.select[iSymm][iAtom]) {
2280             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2281               continue;
2282             }
2283             xyz[0] = coordinates[iSymm][0][iAtom];
2284             xyz[1] = coordinates[iSymm][1][iAtom];
2285             xyz[2] = coordinates[iSymm][2][iAtom];
2286             crystal.toFractionalCoordinates(xyz, uvw);
2287             final double frz = fftZ * uvw[2];
2288             final int ifrz = (int) frz;
2289             final int previousZ = zyAtListBuild[iSymm][iAtom][0];
2290 
2291             final double fry = fftY * uvw[1];
2292             final int ifry = (int) fry;
2293             final int previousY = zyAtListBuild[iSymm][iAtom][1];
2294             if (abs(ifrz - previousZ) >= buff || abs(ifry - previousY) >= buff) {
2295               return true;
2296             }
2297           }
2298         }
2299       }
2300       return false;
2301     }
2302 
2303     @Override
2304     public void finish() {
2305       timer += System.nanoTime();
2306       for (int i = 0; i < fftZ * fftY; i++) {
2307         optRowWeightSolvent.addAndGet(i, optLocal[i]);
2308       }
2309     }
2310 
2311     public int getNumberofSlices() {
2312       return (previousUB - previousLB + 1);
2313     }
2314 
2315     public double getThreadTime() {
2316       return timer;
2317     }
2318 
2319     public int getThreadWeight() {
2320       return actualWeight;
2321     }
2322 
2323     @Override
2324     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2325 
2326       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2327         return;
2328       }
2329 
2330       if (lambdaTerm && atoms[iAtom].applyLambda()) {
2331         return;
2332       }
2333 
2334       xyz[0] = coordinates[iSymm][0][iAtom];
2335       xyz[1] = coordinates[iSymm][1][iAtom];
2336       xyz[2] = coordinates[iSymm][2][iAtom];
2337       FormFactor formFactor = solventFormFactors[iSymm][iAtom];
2338       crystal.toFractionalCoordinates(xyz, uvw);
2339       double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
2340       int frad = aRadGrid;
2341       switch (solventModel) {
2342         case BINARY:
2343           frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
2344           break;
2345         case GAUSSIAN:
2346           frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
2347           break;
2348         case POLYNOMIAL:
2349           frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
2350           break;
2351         case NONE:
2352         default:
2353           return;
2354       }
2355 
2356       int ubZ = rowRegion.zFromRowIndex(ub);
2357       int lbZ = rowRegion.zFromRowIndex(lb);
2358 
2359       // Logic to loop within the cutoff box.
2360       final double frx = fftX * uvw[0];
2361       final int ifrx = (int) frx;
2362       final int ifrxu = ifrx + frad;
2363 
2364       final double fry = fftY * uvw[1];
2365       final int ifry = (int) fry;
2366       final int ifryu = ifry + frad;
2367 
2368       final double frz = fftZ * uvw[2];
2369       final int ifrz = (int) frz;
2370       final int ifrzu = ifrz + frad;
2371 
2372       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2373         int giz = mod(iz, fftZ);
2374         if (lbZ > giz || giz > ubZ) {
2375           continue;
2376         }
2377         xf[2] = iz * ifftZ;
2378         for (int iy = ifry - frad; iy <= ifryu; iy++) {
2379           int giy = mod(iy, fftY);
2380           int rowIndex = rowRegion.rowIndexForYZ(giy, giz);
2381           if (lb > rowIndex || rowIndex > ub) {
2382             continue;
2383           }
2384           xf[1] = iy * ifftY;
2385           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2386             int gix = mod(ix, fftX);
2387             xf[0] = ix * ifftX;
2388             crystal.toCartesianCoordinates(xf, xc);
2389             optLocal[rowIndex]++;
2390             actualWeight++;
2391             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
2392             grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
2393           }
2394         }
2395       }
2396     }
2397 
2398     @Override
2399     public void run(int lb, int ub) throws Exception {
2400       boolean boundsChange = false;
2401       if (previousLB != lb || previousUB != ub) {
2402         boundsChange = true;
2403       }
2404       previousLB = lb;
2405       previousUB = ub;
2406       if (rebuildList || boundsChange) {
2407         buildListA = new ArrayList<>();
2408         buildListS = new ArrayList<>();
2409         for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2410           for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2411             if (rowRegion.select[iSymm][iAtom]) {
2412               buildList(iSymm, iAtom, lb, ub);
2413             }
2414           }
2415         }
2416       }
2417       for (int i = 0; i < buildListA.size(); i++) {
2418         if (rowRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2419           gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2420         }
2421       }
2422     }
2423 
2424     @Override
2425     public void saveZYValues(int[][][] zyAtListBuild) {
2426       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2427         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2428           zyAtListBuild[iSymm][iAtom][0] = 0;
2429           zyAtListBuild[iSymm][iAtom][1] = 0;
2430           if (rowRegion.select[iSymm][iAtom]) {
2431             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2432               continue;
2433             }
2434             xyz[0] = coordinates[iSymm][0][iAtom];
2435             xyz[1] = coordinates[iSymm][1][iAtom];
2436             xyz[2] = coordinates[iSymm][2][iAtom];
2437             crystal.toFractionalCoordinates(xyz, uvw);
2438             final double frz = fftZ * uvw[2];
2439             final int ifrz = (int) frz;
2440 
2441             final double fry = fftY * uvw[1];
2442             final int ifry = (int) fry;
2443             zyAtListBuild[iSymm][iAtom][0] = ifrz;
2444             zyAtListBuild[iSymm][iAtom][1] = ifry;
2445           }
2446         }
2447       }
2448     }
2449 
2450     @Override
2451     public IntegerSchedule schedule() {
2452       return solventRowSchedule;
2453     }
2454 
2455     @Override
2456     public void start() {
2457       fill(optLocal, 0);
2458       timer = -System.nanoTime();
2459     }
2460   }
2461 
2462   private class BulkSolventRowLoop extends RowLoop {
2463 
2464     final double[] xyz = new double[3];
2465     final double[] uvw = new double[3];
2466     final double[] xc = new double[3];
2467     final double[] xf = new double[3];
2468     final double[] grid;
2469     final int[] optLocal;
2470     long timer;
2471     int[] threadWeights;
2472     int[] threadBounds;
2473 
2474     public BulkSolventRowLoop(RowRegion region) {
2475       super(region.getNatoms(), region.getNsymm(), region);
2476       grid = region.getGrid();
2477       optLocal = new int[fftZ * fftY];
2478     }
2479 
2480     @Override
2481     public void finish() {
2482       for (int i = 0; i < fftZ * fftY; i++) {
2483         optRowWeightBulkSolvent.addAndGet(i, optLocal[i]);
2484       }
2485       timer += System.nanoTime();
2486       threadBounds = bulkSolventRowSchedule.getLowerBounds().clone();
2487       threadWeights = bulkSolventRowSchedule.getThreadWeights().clone();
2488       for (int i = threadCount - 1; i > 0; i--) {
2489         threadWeights[i] -= threadWeights[i - 1];
2490       }
2491     }
2492 
2493     public int[] getBoundsPerThread() {
2494       return threadBounds;
2495     }
2496 
2497     public double getTimePerThread() {
2498       return timer;
2499     }
2500 
2501     public int[] getWeightPerThread() {
2502       return optLocal;
2503     }
2504 
2505     @Override
2506     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2507       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2508         return;
2509       }
2510 
2511       if (lambdaTerm && atoms[iAtom].applyLambda()) {
2512         return;
2513       }
2514 
2515       xyz[0] = coordinates[iSymm][0][iAtom];
2516       xyz[1] = coordinates[iSymm][1][iAtom];
2517       xyz[2] = coordinates[iSymm][2][iAtom];
2518       FormFactor formFactor = solventFormFactors[iSymm][iAtom];
2519       crystal.toFractionalCoordinates(xyz, uvw);
2520       double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
2521       int frad = aRadGrid;
2522       switch (solventModel) {
2523         case BINARY:
2524           frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
2525           break;
2526         case GAUSSIAN:
2527           frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
2528           break;
2529         case POLYNOMIAL:
2530           frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
2531           break;
2532         case NONE:
2533         default:
2534           return;
2535       }
2536 
2537       int ubZ = rowRegion.zFromRowIndex(ub);
2538       int lbZ = rowRegion.zFromRowIndex(lb);
2539 
2540       // Logic to loop within the cutoff box.
2541       final double frx = fftX * uvw[0];
2542       final int ifrx = (int) frx;
2543       final int ifrxu = ifrx + frad;
2544 
2545       final double fry = fftY * uvw[1];
2546       final int ifry = (int) fry;
2547       final int ifryu = ifry + frad;
2548 
2549       final double frz = fftZ * uvw[2];
2550       final int ifrz = (int) frz;
2551       final int ifrzu = ifrz + frad;
2552 
2553       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2554         int giz = mod(iz, fftZ);
2555         if (lbZ > giz || giz > ubZ) {
2556           continue;
2557         }
2558         xf[2] = iz * ifftZ;
2559         for (int iy = ifry - frad; iy <= ifryu; iy++) {
2560           int giy = mod(iy, fftY);
2561           int rowIndex = rowRegion.rowIndexForYZ(giy, giz);
2562           if (lb > rowIndex || rowIndex > ub) {
2563             continue;
2564           }
2565           xf[1] = iy * ifftY;
2566           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2567             int gix = mod(ix, fftX);
2568             xf[0] = ix * ifftX;
2569             crystal.toCartesianCoordinates(xf, xc);
2570             optLocal[rowIndex]++;
2571             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
2572             grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
2573           }
2574         }
2575       }
2576     }
2577 
2578     @Override
2579     public IntegerSchedule schedule() {
2580       return bulkSolventRowSchedule;
2581     }
2582 
2583     @Override
2584     public void start() {
2585       timer = -System.nanoTime();
2586       fill(optLocal, 0);
2587     }
2588   }
2589 
2590   private class AtomicSliceLoop extends SliceLoop {
2591 
2592     final double[] xyz = new double[3];
2593     final double[] uvw = new double[3];
2594     final double[] xc = new double[3];
2595     final double[] xf = new double[3];
2596     final double[] grid;
2597     final int[] optLocal;
2598     long timer;
2599     int previousUB, previousLB;
2600     int actualWeight;
2601 
2602     public AtomicSliceLoop(SliceRegion region) {
2603       super(region.getNatoms(), region.getNsymm(), region);
2604       grid = region.getGrid();
2605       optLocal = new int[fftZ];
2606     }
2607 
2608     public void buildList(int iSymm, int iAtom, int lb, int ub) {
2609       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2610         return;
2611       }
2612       xyz[0] = coordinates[iSymm][0][iAtom];
2613       xyz[1] = coordinates[iSymm][1][iAtom];
2614       xyz[2] = coordinates[iSymm][2][iAtom];
2615       crystal.toFractionalCoordinates(xyz, uvw);
2616       final int frad =
2617           min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2618       final double frz = fftZ * uvw[2];
2619       final int ifrz = (int) frz;
2620       final int ifrzu = ifrz + frad;
2621       final int ifrzl = ifrz - frad;
2622       // Loop over allowed z coordinates for this Loop
2623       // Check if the current atom is close enough
2624       // If so, add to list.
2625       int buff = bufferSize;
2626       for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2627         int giz = mod(iz, fftZ);
2628         if (giz >= lb && giz <= ub) {
2629           buildListA.add(iAtom);
2630           buildListS.add(iSymm);
2631           break;
2632         }
2633       }
2634     }
2635 
2636     @Override
2637     public boolean checkList(int[][] zAtListBuild, int buff) {
2638       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2639         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2640           if (sliceRegion.select[iSymm][iAtom]) {
2641             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2642               continue;
2643             }
2644             xyz[0] = coordinates[iSymm][0][iAtom];
2645             xyz[1] = coordinates[iSymm][1][iAtom];
2646             xyz[2] = coordinates[iSymm][2][iAtom];
2647             crystal.toFractionalCoordinates(xyz, uvw);
2648             final double frz = fftZ * uvw[2];
2649             final int ifrz = (int) frz;
2650             final int previousZ = zAtListBuild[iSymm][iAtom];
2651             if (abs(ifrz - previousZ) >= buff) {
2652               return true;
2653             }
2654           }
2655         }
2656       }
2657       return false;
2658     }
2659 
2660     @Override
2661     public void finish() {
2662       for (int i = 0; i < fftZ; i++) {
2663         optSliceWeightAtomic.addAndGet(i, optLocal[i]);
2664       }
2665       timer += System.nanoTime();
2666     }
2667 
2668     public int getNumberofSlices() {
2669       return (previousUB - previousLB + 1);
2670     }
2671 
2672     public double getThreadTime() {
2673       return timer;
2674     }
2675 
2676     public int getThreadWeight() {
2677       return actualWeight;
2678     }
2679 
2680     @Override
2681     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2682       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2683         return;
2684       }
2685 
2686       if (lambdaTerm && atoms[iAtom].applyLambda()) {
2687         return;
2688       }
2689 
2690       xyz[0] = coordinates[iSymm][0][iAtom];
2691       xyz[1] = coordinates[iSymm][1][iAtom];
2692       xyz[2] = coordinates[iSymm][2][iAtom];
2693       FormFactor atomff = atomFormFactors[iSymm][iAtom];
2694       crystal.toFractionalCoordinates(xyz, uvw);
2695       final int frad =
2696           min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2697 
2698       // Logic to loop within the cutoff box.
2699       final double frx = fftX * uvw[0];
2700       final int ifrx = (int) frx;
2701       final int ifrxu = ifrx + frad;
2702 
2703       final double fry = fftY * uvw[1];
2704       final int ifry = (int) fry;
2705       final int ifryu = ifry + frad;
2706 
2707       final double frz = fftZ * uvw[2];
2708       final int ifrz = (int) frz;
2709       final int ifrzu = ifrz + frad;
2710 
2711       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2712         int giz = mod(iz, fftZ);
2713         if (lb > giz || giz > ub) {
2714           continue;
2715         }
2716         xf[2] = iz * ifftZ;
2717         for (int iy = ifry - frad; iy <= ifryu; iy++) {
2718           int giy = mod(iy, fftY);
2719           xf[1] = iy * ifftY;
2720           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2721             int gix = mod(ix, fftX);
2722             xf[0] = ix * ifftX;
2723             crystal.toCartesianCoordinates(xf, xc);
2724             optLocal[giz]++;
2725             actualWeight++;
2726             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
2727             grid[ii] = atomff.rho(grid[ii], 1.0, xc);
2728           }
2729         }
2730       }
2731     }
2732 
2733     @Override
2734     public void run(int lb, int ub) throws Exception {
2735       boolean boundsChange = false;
2736       if (previousLB != lb || previousUB != ub) {
2737         boundsChange = true;
2738       }
2739       previousLB = lb;
2740       previousUB = ub;
2741       if (rebuildList || boundsChange) {
2742         buildListA = new ArrayList<>();
2743         buildListS = new ArrayList<>();
2744         for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2745           for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2746             if (sliceRegion.select[iSymm][iAtom]) {
2747               buildList(iSymm, iAtom, lb, ub);
2748             }
2749           }
2750         }
2751       }
2752       for (int i = 0; i < buildListA.size(); i++) {
2753         if (sliceRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2754           gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2755         }
2756       }
2757     }
2758 
2759     @Override
2760     public void saveZValues(int[][] zAtListBuild) {
2761       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2762         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2763           zAtListBuild[iSymm][iAtom] = 0;
2764           if (sliceRegion.select[iSymm][iAtom]) {
2765             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2766               continue;
2767             }
2768             xyz[0] = coordinates[iSymm][0][iAtom];
2769             xyz[1] = coordinates[iSymm][1][iAtom];
2770             xyz[2] = coordinates[iSymm][2][iAtom];
2771             crystal.toFractionalCoordinates(xyz, uvw);
2772             final double frz = fftZ * uvw[2];
2773             final int ifrz = (int) frz;
2774             zAtListBuild[iSymm][iAtom] = ifrz;
2775           }
2776         }
2777       }
2778     }
2779 
2780     @Override
2781     public IntegerSchedule schedule() {
2782       return atomicSliceSchedule;
2783     }
2784 
2785     @Override
2786     public void start() {
2787       fill(optLocal, 0);
2788       actualWeight = 0;
2789       timer = -System.nanoTime();
2790     }
2791   }
2792 
2793   private class SolventSliceLoop extends SliceLoop {
2794 
2795     final double[] xyz = new double[3];
2796     final double[] uvw = new double[3];
2797     final double[] xc = new double[3];
2798     final double[] xf = new double[3];
2799     final double[] grid;
2800     final int[] optLocal;
2801     long timer;
2802     int previousUB, previousLB;
2803     int actualWeight;
2804 
2805     public SolventSliceLoop(SliceRegion region) {
2806       super(region.getNatoms(), region.getNsymm(), region);
2807       grid = region.getGrid();
2808       optLocal = new int[fftZ];
2809     }
2810 
2811     public void buildList(int iSymm, int iAtom, int lb, int ub) {
2812       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2813         return;
2814       }
2815       xyz[0] = coordinates[iSymm][0][iAtom];
2816       xyz[1] = coordinates[iSymm][1][iAtom];
2817       xyz[2] = coordinates[iSymm][2][iAtom];
2818       crystal.toFractionalCoordinates(xyz, uvw);
2819       final int frad =
2820           min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);
2821       final double frz = fftZ * uvw[2];
2822       final int ifrz = (int) frz;
2823       final int ifrzu = ifrz + frad;
2824       final int ifrzl = ifrz - frad;
2825       // Loop over allowed z coordinates for this Loop
2826       // Check if the current atom is close enough
2827       // If so, add to list.
2828       int buff = bufferSize;
2829       for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
2830         int giz = mod(iz, fftZ);
2831         if (giz >= lb && giz <= ub) {
2832           buildListA.add(iAtom);
2833           buildListS.add(iSymm);
2834           break;
2835         }
2836       }
2837     }
2838 
2839     @Override
2840     public boolean checkList(int[][] zAtListBuild, int buff) {
2841       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2842         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2843           if (sliceRegion.select[iSymm][iAtom]) {
2844             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2845               continue;
2846             }
2847             xyz[0] = coordinates[iSymm][0][iAtom];
2848             xyz[1] = coordinates[iSymm][1][iAtom];
2849             xyz[2] = coordinates[iSymm][2][iAtom];
2850             crystal.toFractionalCoordinates(xyz, uvw);
2851             final double frz = fftZ * uvw[2];
2852             final int ifrz = (int) frz;
2853             final int previousZ = zAtListBuild[iSymm][iAtom];
2854             if (abs(ifrz - previousZ) >= buff) {
2855               return true;
2856             }
2857           }
2858         }
2859       }
2860       return false;
2861     }
2862 
2863     @Override
2864     public void finish() {
2865       timer += System.nanoTime();
2866       for (int i = 0; i < fftZ; i++) {
2867         optSliceWeightSolvent.addAndGet(i, optLocal[i]);
2868       }
2869     }
2870 
2871     public int getNumberofSlices() {
2872       return (previousUB - previousLB + 1);
2873     }
2874 
2875     public double getThreadTime() {
2876       return timer;
2877     }
2878 
2879     public int getThreadWeight() {
2880       return actualWeight;
2881     }
2882 
2883     @Override
2884     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
2885       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2886         return;
2887       }
2888 
2889       if (lambdaTerm && atoms[iAtom].applyLambda()) {
2890         return;
2891       }
2892 
2893       xyz[0] = coordinates[iSymm][0][iAtom];
2894       xyz[1] = coordinates[iSymm][1][iAtom];
2895       xyz[2] = coordinates[iSymm][2][iAtom];
2896       FormFactor formFactor = solventFormFactors[iSymm][iAtom];
2897       crystal.toFractionalCoordinates(xyz, uvw);
2898       double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
2899       int frad = aRadGrid;
2900       switch (solventModel) {
2901         case BINARY:
2902           frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
2903           break;
2904         case GAUSSIAN:
2905           frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
2906           break;
2907         case POLYNOMIAL:
2908           frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
2909           break;
2910         case NONE:
2911         default:
2912           return;
2913       }
2914 
2915       // Logic to loop within the cutoff box.
2916       final double frx = fftX * uvw[0];
2917       final int ifrx = (int) frx;
2918       final int ifrxu = ifrx + frad;
2919 
2920       final double fry = fftY * uvw[1];
2921       final int ifry = (int) fry;
2922       final int ifryu = ifry + frad;
2923 
2924       final double frz = fftZ * uvw[2];
2925       final int ifrz = (int) frz;
2926       final int ifrzu = ifrz + frad;
2927 
2928       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
2929         int giz = mod(iz, fftZ);
2930         if (lb > giz || giz > ub) {
2931           continue;
2932         }
2933         xf[2] = iz * ifftZ;
2934         for (int iy = ifry - frad; iy <= ifryu; iy++) {
2935           int giy = mod(iy, fftY);
2936           xf[1] = iy * ifftY;
2937           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
2938             int gix = mod(ix, fftX);
2939             xf[0] = ix * ifftX;
2940             crystal.toCartesianCoordinates(xf, xc);
2941             optLocal[giz]++;
2942             actualWeight++;
2943             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
2944             grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
2945           }
2946         }
2947       }
2948     }
2949 
2950     @Override
2951     public void run(int lb, int ub) throws Exception {
2952       boolean boundsChange = false;
2953       if (previousLB != lb || previousUB != ub) {
2954         boundsChange = true;
2955       }
2956       previousLB = lb;
2957       previousUB = ub;
2958       if (rebuildList || boundsChange) {
2959         buildListA = new ArrayList<>();
2960         buildListS = new ArrayList<>();
2961         for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2962           for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2963             if (sliceRegion.select[iSymm][iAtom]) {
2964               buildList(iSymm, iAtom, lb, ub);
2965             }
2966           }
2967         }
2968       }
2969       for (int i = 0; i < buildListA.size(); i++) {
2970         if (sliceRegion.select[buildListS.get(i)][buildListA.get(i)]) {
2971           gridDensity(buildListS.get(i), buildListA.get(i), lb, ub);
2972         }
2973       }
2974     }
2975 
2976     @Override
2977     public void saveZValues(int[][] zAtListBuild) {
2978       for (int iSymm = 0; iSymm < nSymm; iSymm++) {
2979         for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
2980           zAtListBuild[iSymm][iAtom] = 0;
2981           if (sliceRegion.select[iSymm][iAtom]) {
2982             if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
2983               continue;
2984             }
2985             xyz[0] = coordinates[iSymm][0][iAtom];
2986             xyz[1] = coordinates[iSymm][1][iAtom];
2987             xyz[2] = coordinates[iSymm][2][iAtom];
2988             crystal.toFractionalCoordinates(xyz, uvw);
2989             final double frz = fftZ * uvw[2];
2990             final int ifrz = (int) frz;
2991             zAtListBuild[iSymm][iAtom] = ifrz;
2992           }
2993         }
2994       }
2995     }
2996 
2997     @Override
2998     public IntegerSchedule schedule() {
2999       return solventSliceSchedule;
3000     }
3001 
3002     @Override
3003     public void start() {
3004       fill(optLocal, 0);
3005       actualWeight = 0;
3006       timer = -System.nanoTime();
3007     }
3008   }
3009 
3010   private class BulkSolventSliceLoop extends SliceLoop {
3011 
3012     final double[] xyz = new double[3];
3013     final double[] uvw = new double[3];
3014     final double[] xc = new double[3];
3015     final double[] xf = new double[3];
3016     final double[] grid;
3017     final int[] optLocal;
3018     long timer;
3019     int[] threadBounds;
3020     int[] threadWeights;
3021 
3022     public BulkSolventSliceLoop(SliceRegion region) {
3023       super(region.getNatoms(), region.getNsymm(), region);
3024       grid = region.getGrid();
3025       optLocal = new int[fftZ];
3026     }
3027 
3028     @Override
3029     public void finish() {
3030       timer += System.nanoTime();
3031       for (int i = 0; i < fftZ; i++) {
3032         optSliceWeightBulkSolvent.addAndGet(i, optLocal[i]);
3033       }
3034       threadBounds = bulkSolventSliceSchedule.getLowerBounds().clone();
3035       threadWeights = bulkSolventSliceSchedule.getThreadWeights().clone();
3036       for (int i = threadBounds.length - 1; i > 0; i--) {
3037         threadBounds[i] -= threadBounds[i - 1];
3038       }
3039     }
3040 
3041     public int[] getBoundsPerThread() {
3042       return threadBounds;
3043     }
3044 
3045     public double getTimePerThread() {
3046       return timer;
3047     }
3048 
3049     public int[] getWeightPerThread() {
3050       return threadWeights;
3051     }
3052 
3053     @Override
3054     public void gridDensity(int iSymm, int iAtom, int lb, int ub) {
3055       if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
3056         return;
3057       }
3058 
3059       if (lambdaTerm && atoms[iAtom].applyLambda()) {
3060         return;
3061       }
3062 
3063       xyz[0] = coordinates[iSymm][0][iAtom];
3064       xyz[1] = coordinates[iSymm][1][iAtom];
3065       xyz[2] = coordinates[iSymm][2][iAtom];
3066       FormFactor formFactor = solventFormFactors[iSymm][iAtom];
3067       crystal.toFractionalCoordinates(xyz, uvw);
3068       double vdwr = atoms[iAtom].getVDWType().radius * 0.5;
3069       int frad = aRadGrid;
3070       switch (solventModel) {
3071         case BINARY:
3072           frad = min(aRadGrid, (int) floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
3073           break;
3074         case GAUSSIAN:
3075           frad = min(aRadGrid, (int) floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
3076           break;
3077         case POLYNOMIAL:
3078           frad = min(aRadGrid, (int) floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
3079           break;
3080         case NONE:
3081         default:
3082           return;
3083       }
3084 
3085       // Logic to loop within the cutoff box.
3086       final double frx = fftX * uvw[0];
3087       final int ifrx = (int) frx;
3088       final int ifrxu = ifrx + frad;
3089 
3090       final double fry = fftY * uvw[1];
3091       final int ifry = (int) fry;
3092       final int ifryu = ifry + frad;
3093 
3094       final double frz = fftZ * uvw[2];
3095       final int ifrz = (int) frz;
3096       final int ifrzu = ifrz + frad;
3097 
3098       for (int iz = ifrz - frad; iz <= ifrzu; iz++) {
3099         int giz = mod(iz, fftZ);
3100         if (lb > giz || giz > ub) {
3101           continue;
3102         }
3103         xf[2] = iz * ifftZ;
3104         for (int iy = ifry - frad; iy <= ifryu; iy++) {
3105           int giy = mod(iy, fftY);
3106           xf[1] = iy * ifftY;
3107           for (int ix = ifrx - frad; ix <= ifrxu; ix++) {
3108             int gix = mod(ix, fftX);
3109             xf[0] = ix * ifftX;
3110             crystal.toCartesianCoordinates(xf, xc);
3111             optLocal[giz]++;
3112             final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
3113             grid[ii] = formFactor.rho(grid[ii], 1.0, xc);
3114           }
3115         }
3116       }
3117     }
3118 
3119     @Override
3120     public IntegerSchedule schedule() {
3121       return bulkSolventSliceSchedule;
3122     }
3123 
3124     @Override
3125     public void start() {
3126       fill(optLocal, 0);
3127       timer = -System.nanoTime();
3128     }
3129   }
3130 
3131   private class AtomicGradientRegion extends ParallelRegion {
3132 
3133     private final AtomicGradientLoop[] atomicGradientLoop;
3134     private RefinementMode refinementmode;
3135 
3136     public AtomicGradientRegion(RefinementMode refinementmode) {
3137       this.refinementmode = refinementmode;
3138       atomicGradientLoop = new AtomicGradientLoop[threadCount];
3139       for (int i = 0; i < threadCount; i++) {
3140         atomicGradientLoop[i] = new AtomicGradientLoop();
3141       }
3142     }
3143 
3144     public RefinementMode getRefinementMode() {
3145       return refinementmode;
3146     }
3147 
3148     public void setRefinementMode(RefinementMode refinementmode) {
3149       this.refinementmode = refinementmode;
3150     }
3151 
3152     @Override
3153     public void run() {
3154       try {
3155         execute(0, nAtoms - 1, atomicGradientLoop[getThreadIndex()]);
3156       } catch (Exception e) {
3157         e.printStackTrace();
3158         logger.severe(e.toString());
3159       }
3160     }
3161 
3162     private class AtomicGradientLoop extends IntegerForLoop {
3163 
3164       final double[] xyz = new double[3];
3165       final double[] uvw = new double[3];
3166       final double[] xc = new double[3];
3167       final double[] xf = new double[3];
3168       final int[] optLocal = new int[nAtoms];
3169       long timer;
3170 
3171       // Extra padding to avert cache interference.
3172       long pad0, pad1, pad2, pad3, pad4, pad5, pad6, pad7;
3173       long pad8, pad9, pada, padb, padc, padd, pade, padf;
3174 
3175       @Override
3176       public void finish() {
3177         timer += System.nanoTime();
3178         for (int i = 0; i < nAtoms; i++) {
3179           optAtomicGradientWeight.addAndGet(i, optLocal[i]);
3180         }
3181       }
3182 
3183       @Override
3184       public void run(final int lb, final int ub) {
3185         for (int n = lb; n <= ub; n++) {
3186           if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
3187             continue;
3188           }
3189           if (lambdaTerm && atoms[n].applyLambda()) {
3190             continue;
3191           }
3192 
3193           xyz[0] = coordinates[0][0][n];
3194           xyz[1] = coordinates[0][1][n];
3195           xyz[2] = coordinates[0][2][n];
3196           FormFactor atomff = atomFormFactors[0][n];
3197           atomff.update(xyz, 0.0);
3198           crystal.toFractionalCoordinates(xyz, uvw);
3199           final int dfrad =
3200               Math.min(
3201                   aRadGrid, (int) Math.floor(atoms[n].getFormFactorWidth() * fftX / crystal.a) + 1);
3202 
3203           // Logic to loop within the cutoff box.
3204           final double frx = fftX * uvw[0];
3205           final int ifrx = (int) frx;
3206 
3207           final double fry = fftY * uvw[1];
3208           final int ifry = (int) fry;
3209 
3210           final double frz = fftZ * uvw[2];
3211           final int ifrz = (int) frz;
3212 
3213           for (int ix = ifrx - dfrad; ix <= ifrx + dfrad; ix++) {
3214             int gix = mod(ix, fftX);
3215             xf[0] = ix * ifftX;
3216             for (int iy = ifry - dfrad; iy <= ifry + dfrad; iy++) {
3217               int giy = mod(iy, fftY);
3218               xf[1] = iy * ifftY;
3219               for (int iz = ifrz - dfrad; iz <= ifrz + dfrad; iz++) {
3220                 int giz = mod(iz, fftZ);
3221                 xf[2] = iz * ifftZ;
3222                 crystal.toCartesianCoordinates(xf, xc);
3223                 optLocal[n]++;
3224                 final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
3225                 atomff.rhoGrad(xc, weight * getDensityGrid()[ii], refinementmode);
3226               }
3227             }
3228           }
3229         }
3230       }
3231 
3232       @Override
3233       public IntegerSchedule schedule() {
3234         return atomicGradientSchedule;
3235       }
3236 
3237       @Override
3238       public void start() {
3239         fill(optLocal, 0);
3240         timer = -System.nanoTime();
3241       }
3242     }
3243   }
3244 
3245   private class SolventGradientRegion extends ParallelRegion {
3246 
3247     private final SolventGradientLoop[] solventGradientLoop;
3248     private RefinementMode refinementmode;
3249 
3250     public SolventGradientRegion(RefinementMode refinementmode) {
3251       this.refinementmode = refinementmode;
3252       solventGradientLoop = new SolventGradientLoop[threadCount];
3253       for (int i = 0; i < threadCount; i++) {
3254         solventGradientLoop[i] = new SolventGradientLoop();
3255       }
3256     }
3257 
3258     public RefinementMode getRefinementMode() {
3259       return refinementmode;
3260     }
3261 
3262     public void setRefinementMode(RefinementMode refinementmode) {
3263       this.refinementmode = refinementmode;
3264     }
3265 
3266     @Override
3267     public void run() {
3268       try {
3269         execute(0, nAtoms - 1, solventGradientLoop[getThreadIndex()]);
3270       } catch (Exception e) {
3271         logger.severe(e.toString());
3272       }
3273     }
3274 
3275     private class SolventGradientLoop extends IntegerForLoop {
3276 
3277       double[] xyz = new double[3];
3278       double[] uvw = new double[3];
3279       double[] xc = new double[3];
3280       double[] xf = new double[3];
3281       // Extra padding to avert cache interference.
3282       long pad0, pad1, pad2, pad3, pad4, pad5, pad6, pad7;
3283       long pad8, pad9, pada, padb, padc, padd, pade, padf;
3284 
3285       @Override
3286       public void run(final int lb, final int ub) {
3287         for (int n = lb; n <= ub; n++) {
3288           if (!atoms[n].getUse() && !nativeEnvironmentApproximation) {
3289             continue;
3290           }
3291           if (lambdaTerm && atoms[n].applyLambda()) {
3292             continue;
3293           }
3294           xyz[0] = coordinates[0][0][n];
3295           xyz[1] = coordinates[0][1][n];
3296           xyz[2] = coordinates[0][2][n];
3297           FormFactor solventff = solventFormFactors[0][n];
3298           solventff.update(xyz);
3299           crystal.toFractionalCoordinates(xyz, uvw);
3300           double vdwr = atoms[n].getVDWType().radius * 0.5;
3301           double dfcmult = 1.0;
3302           int dfrad = aRadGrid;
3303           switch (solventModel) {
3304             case BINARY:
3305               dfrad =
3306                   Math.min(
3307                       aRadGrid, (int) Math.floor((vdwr + solventA + 0.2) * fftX / crystal.a) + 1);
3308               break;
3309             case GAUSSIAN:
3310               dfrad =
3311                   Math.min(
3312                       aRadGrid, (int) Math.floor((vdwr * solventB + 2.0) * fftX / crystal.a) + 1);
3313               dfcmult = solventA;
3314               break;
3315             case POLYNOMIAL:
3316               dfrad =
3317                   Math.min(
3318                       aRadGrid, (int) Math.floor((vdwr + solventB + 0.2) * fftX / crystal.a) + 1);
3319               break;
3320             case NONE:
3321               return;
3322           }
3323 
3324           // Logic to loop within the cutoff box.
3325           final double frx = fftX * uvw[0];
3326           final int ifrx = (int) frx;
3327           final int ifrxu = ifrx + dfrad;
3328 
3329           final double fry = fftY * uvw[1];
3330           final int ifry = (int) fry;
3331           final int ifryu = ifry + dfrad;
3332 
3333           final double frz = fftZ * uvw[2];
3334           final int ifrz = (int) frz;
3335           final int ifrzu = ifrz + dfrad;
3336 
3337           for (int ix = ifrx - dfrad; ix <= ifrxu; ix++) {
3338             int gix = mod(ix, fftX);
3339             xf[0] = ix * ifftX;
3340             for (int iy = ifry - dfrad; iy <= ifryu; iy++) {
3341               int giy = mod(iy, fftY);
3342               xf[1] = iy * ifftY;
3343               for (int iz = ifrz - dfrad; iz <= ifrzu; iz++) {
3344                 int giz = mod(iz, fftZ);
3345                 xf[2] = iz * ifftZ;
3346                 crystal.toCartesianCoordinates(xf, xc);
3347                 final int ii = iComplex3D(gix, giy, giz, fftX, fftY);
3348                 solventff.rhoGrad(
3349                     xc,
3350                     weight * getDensityGrid()[ii] * dfcmult * getSolventGrid()[ii],
3351                     refinementmode);
3352               }
3353             }
3354           }
3355         }
3356       }
3357     }
3358   }
3359 
3360   private class BabinetRegion extends ParallelRegion {
3361 
3362     BabinetLoop[] babinetLoops;
3363 
3364     public BabinetRegion(int nThreads) {
3365       babinetLoops = new BabinetLoop[nThreads];
3366     }
3367 
3368     @Override
3369     public void run() throws Exception {
3370       int ti = getThreadIndex();
3371 
3372       if (babinetLoops[ti] == null) {
3373         babinetLoops[ti] = new BabinetLoop();
3374       }
3375 
3376       try {
3377         execute(0, fftZ - 1, babinetLoops[ti]);
3378       } catch (Exception e) {
3379         logger.info(e.toString());
3380       }
3381     }
3382 
3383     private class BabinetLoop extends IntegerForLoop {
3384 
3385       @Override
3386       public void run(int lb, int ub) throws Exception {
3387         switch (solventModel) {
3388           case BINARY:
3389           case POLYNOMIAL:
3390             for (int k = lb; k <= ub; k++) {
3391               for (int j = 0; j < fftY; j++) {
3392                 for (int i = 0; i < fftX; i++) {
3393                   final int ii = iComplex3D(i, j, k, fftX, fftY);
3394                   densityGrid[ii] = 1.0 - getDensityGrid()[ii];
3395                 }
3396               }
3397             }
3398             break;
3399           case GAUSSIAN:
3400             for (int k = lb; k <= ub; k++) {
3401               for (int j = 0; j < fftY; j++) {
3402                 for (int i = 0; i < fftX; i++) {
3403                   final int ii = iComplex3D(i, j, k, fftX, fftY);
3404                   densityGrid[ii] = 1.0 - exp(-solventA * getDensityGrid()[ii]);
3405                 }
3406               }
3407             }
3408             break;
3409           default:
3410         }
3411       }
3412     }
3413   }
3414 
3415   private class SolventGridRegion extends ParallelRegion {
3416 
3417     SolventGridLoop[] solventGridLoops;
3418 
3419     public SolventGridRegion(int nThreads) {
3420       solventGridLoops = new SolventGridLoop[nThreads];
3421     }
3422 
3423     @Override
3424     public void run() throws Exception {
3425       int ti = getThreadIndex();
3426 
3427       if (solventGridLoops[ti] == null) {
3428         solventGridLoops[ti] = new SolventGridLoop();
3429       }
3430 
3431       try {
3432         execute(0, fftZ - 1, solventGridLoops[ti]);
3433       } catch (Exception e) {
3434         logger.info(e.toString());
3435       }
3436     }
3437 
3438     private class SolventGridLoop extends IntegerForLoop {
3439 
3440       @Override
3441       public void run(int lb, int ub) throws Exception {
3442         // case GAUSSIAN:
3443         for (int k = lb; k <= ub; k++) {
3444           for (int j = 0; j < fftY; j++) {
3445             for (int i = 0; i < fftX; i++) {
3446               final int ii = iComplex3D(i, j, k, fftX, fftY);
3447               solventGrid[ii] = exp(-solventA * getSolventGrid()[ii]);
3448             }
3449           }
3450         }
3451       }
3452     }
3453   }
3454 
3455   private class InitRegion extends ParallelRegion {
3456 
3457     InitLoop[] initLoops;
3458     FormFactorUpdateLoop[] formFactorUpdateLoops;
3459     int nHKL = reflectionList.hkllist.size();
3460     double[][] hkldata = null;
3461 
3462     public InitRegion(int nThreads) {
3463       initLoops = new InitLoop[nThreads];
3464       formFactorUpdateLoops = new FormFactorUpdateLoop[nThreads];
3465     }
3466 
3467     @Override
3468     public void run() throws Exception {
3469       int ti = getThreadIndex();
3470 
3471       if (initLoops[ti] == null) {
3472         initLoops[ti] = new InitLoop();
3473         formFactorUpdateLoops[ti] = new FormFactorUpdateLoop();
3474       }
3475 
3476       try {
3477         execute(0, nHKL - 1, initLoops[ti]);
3478         execute(0, nAtoms - 1, formFactorUpdateLoops[ti]);
3479       } catch (Exception e) {
3480         logger.info(e.toString());
3481       }
3482     }
3483 
3484     public void setHKL(double[][] hkldata) {
3485       this.hkldata = hkldata;
3486     }
3487 
3488     private class InitLoop extends IntegerForLoop {
3489 
3490       @Override
3491       public void run(int lb, int ub) throws Exception {
3492         for (int i = lb; i <= ub; i++) {
3493           hkldata[i][0] = 0.0;
3494           hkldata[i][1] = 0.0;
3495         }
3496       }
3497     }
3498 
3499     private class FormFactorUpdateLoop extends IntegerForLoop {
3500 
3501       private final double[] xyz = new double[3];
3502 
3503       @Override
3504       public void run(int lb, int ub) throws Exception {
3505         for (int iSymm = 0; iSymm < bulkNSymm; iSymm++) {
3506           for (int i = lb; i <= ub; i++) {
3507             xyz[0] = coordinates[iSymm][0][i];
3508             xyz[1] = coordinates[iSymm][1][i];
3509             xyz[2] = coordinates[iSymm][2][i];
3510             if (solventFormFactors != null) {
3511               solventFormFactors[iSymm][i].update(xyz);
3512             }
3513             if (atomFormFactors != null) {
3514               atomFormFactors[iSymm][i].update(xyz, bAdd);
3515             }
3516           }
3517         }
3518       }
3519     }
3520   }
3521 
3522   private class AtomicScaleRegion extends ParallelRegion {
3523 
3524     AtomicScaleLoop[] atomicScaleLoops;
3525     int nHKL = reflectionList.hkllist.size();
3526     double[][] hklData = null;
3527     double scale;
3528 
3529     public AtomicScaleRegion(int nThreads) {
3530       atomicScaleLoops = new AtomicScaleLoop[nThreads];
3531     }
3532 
3533     @Override
3534     public void run() throws Exception {
3535       int ti = getThreadIndex();
3536 
3537       if (atomicScaleLoops[ti] == null) {
3538         atomicScaleLoops[ti] = new AtomicScaleLoop();
3539       }
3540 
3541       try {
3542         execute(0, nHKL - 1, atomicScaleLoops[ti]);
3543       } catch (Exception e) {
3544         logger.info(e.toString());
3545       }
3546     }
3547 
3548     public void setHKL(double[][] hkldata) {
3549       this.hklData = hkldata;
3550     }
3551 
3552     public void setScale(double scale) {
3553       this.scale = scale;
3554     }
3555 
3556     private class AtomicScaleLoop extends IntegerForLoop {
3557 
3558       ComplexNumber c;
3559 
3560       public AtomicScaleLoop() {
3561         c = new ComplexNumber();
3562       }
3563 
3564       @Override
3565       public void run(int lb, int ub) throws Exception {
3566         for (int i = lb; i <= ub; i++) {
3567           HKL ih = reflectionList.hkllist.get(i);
3568           double[] fc = hklData[ih.getIndex()];
3569           c.re(fc[0]);
3570           c.im(fc[1]);
3571           // Remove Badd
3572           double s = crystal.invressq(ih);
3573           c.timesIP(scale * exp(0.25 * bAdd * s));
3574           c.conjugateIP();
3575           fc[0] = c.re();
3576           fc[1] = c.im();
3577         }
3578       }
3579     }
3580   }
3581 
3582   private class SolventScaleRegion extends ParallelRegion {
3583 
3584     SolventScaleLoop[] solventScaleLoops;
3585     int nHKL = reflectionList.hkllist.size();
3586     double[][] hkldata = null;
3587     double scale;
3588 
3589     public SolventScaleRegion(int nThreads) {
3590       solventScaleLoops = new SolventScaleLoop[nThreads];
3591     }
3592 
3593     @Override
3594     public void run() throws Exception {
3595       int ti = getThreadIndex();
3596 
3597       if (solventScaleLoops[ti] == null) {
3598         solventScaleLoops[ti] = new SolventScaleLoop();
3599       }
3600 
3601       try {
3602         execute(0, nHKL - 1, solventScaleLoops[ti]);
3603       } catch (Exception e) {
3604         logger.info(e.toString());
3605       }
3606     }
3607 
3608     public void setHKL(double[][] hkldata) {
3609       this.hkldata = hkldata;
3610     }
3611 
3612     public void setScale(double scale) {
3613       this.scale = scale;
3614     }
3615 
3616     private class SolventScaleLoop extends IntegerForLoop {
3617 
3618       ComplexNumber c;
3619 
3620       public SolventScaleLoop() {
3621         c = new ComplexNumber();
3622       }
3623 
3624       @Override
3625       public void run(int lb, int ub) throws Exception {
3626         for (int i = lb; i <= ub; i++) {
3627           HKL ih = reflectionList.hkllist.get(i);
3628           double[] fc = hkldata[ih.getIndex()];
3629           c.re(fc[0]);
3630           c.im(fc[1]);
3631           c.timesIP(scale);
3632           c.conjugateIP();
3633           // negative: babinet
3634           fc[0] = -c.re();
3635           fc[1] = -c.im();
3636         }
3637       }
3638     }
3639   }
3640 
3641   private class ExtractRegion extends ParallelRegion {
3642 
3643     ExtractLoop[] extractLoops;
3644     int nHKL = reflectionList.hkllist.size();
3645     double[][] hkldata = null;
3646 
3647     public ExtractRegion(int nThreads) {
3648       extractLoops = new ExtractLoop[nThreads];
3649     }
3650 
3651     @Override
3652     public void run() throws Exception {
3653       int ti = getThreadIndex();
3654 
3655       if (extractLoops[ti] == null) {
3656         extractLoops[ti] = new ExtractLoop();
3657       }
3658 
3659       try {
3660         execute(0, nHKL - 1, extractLoops[ti]);
3661       } catch (Exception e) {
3662         logger.info(e.toString());
3663       }
3664     }
3665 
3666     public void setHKL(double[][] hkldata) {
3667       this.hkldata = hkldata;
3668     }
3669 
3670     private class ExtractLoop extends IntegerForLoop {
3671 
3672       ComplexNumber c;
3673       int nsym;
3674       List<SymOp> symops;
3675       HKL ij;
3676 
3677       public ExtractLoop() {
3678         c = new ComplexNumber();
3679         nsym = crystal.spaceGroup.symOps.size();
3680         symops = crystal.spaceGroup.symOps;
3681         ij = new HKL();
3682       }
3683 
3684       @Override
3685       public void run(int lb, int ub) throws Exception {
3686         for (int i = lb; i <= ub; i++) {
3687           HKL ih = reflectionList.hkllist.get(i);
3688           double[] fc = hkldata[ih.getIndex()];
3689           // Apply symmetry
3690           for (int j = 0; j < nsym; j++) {
3691             SymOp symOp = symops.get(j);
3692             applyTransSymRot(ih, ij, symOp);
3693             double shift = symOp.symPhaseShift(ih);
3694             int h = mod(ij.getH(), fftX);
3695             int k = mod(ij.getK(), fftY);
3696             int l = mod(ij.getL(), fftZ);
3697             if (h < halfFFTX + 1) {
3698               final int ii = iComplex3D(h, k, l, fftX, fftY);
3699               c.re(getDensityGrid()[ii]);
3700               c.im(getDensityGrid()[ii + 1]);
3701               c.phaseShiftIP(shift);
3702               fc[0] += c.re();
3703               fc[1] += c.im();
3704             } else {
3705               h = (fftX - h) % fftX;
3706               k = (fftY - k) % fftY;
3707               l = (fftZ - l) % fftZ;
3708               final int ii = iComplex3D(h, k, l, fftX, fftY);
3709               c.re(getDensityGrid()[ii]);
3710               c.im(-getDensityGrid()[ii + 1]);
3711               c.phaseShiftIP(shift);
3712               fc[0] += c.re();
3713               fc[1] += c.im();
3714             }
3715           }
3716         }
3717       }
3718     }
3719   }
3720 
3721   private class ArrayCopyRegion extends ParallelRegion {
3722 
3723     ArrayCopyLoop[] arrayCopyLoops;
3724 
3725     public ArrayCopyRegion() {
3726       arrayCopyLoops = new ArrayCopyLoop[threadCount];
3727     }
3728 
3729     @Override
3730     public void run() throws Exception {
3731       int ti = getThreadIndex();
3732 
3733       if (arrayCopyLoops[ti] == null) {
3734         arrayCopyLoops[ti] = new ArrayCopyLoop();
3735       }
3736 
3737       try {
3738         execute(0, complexFFT3DSpace - 1, arrayCopyLoops[ti]);
3739       } catch (Exception e) {
3740         logger.info(e.toString());
3741       }
3742     }
3743 
3744     private class ArrayCopyLoop extends IntegerForLoop {
3745 
3746       @Override
3747       public void run(int lb, int ub) throws Exception {
3748         arraycopy(getDensityGrid(), lb, getSolventGrid(), lb, ub - lb);
3749       }
3750     }
3751   }
3752 }