View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.realspace;
39  
40  import static ffx.numerics.math.ScalarMath.mod;
41  import static java.lang.String.format;
42  import static java.util.Arrays.fill;
43  import static org.apache.commons.math3.util.FastMath.floor;
44  
45  import edu.rit.pj.IntegerForLoop;
46  import edu.rit.pj.ParallelRegion;
47  import edu.rit.pj.ParallelTeam;
48  import edu.rit.pj.reduction.SharedDouble;
49  import ffx.crystal.Crystal;
50  import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
51  import ffx.numerics.atomic.AtomicDoubleArray3D;
52  import ffx.numerics.spline.TriCubicSpline;
53  import ffx.potential.MolecularAssembly;
54  import ffx.potential.Utilities;
55  import ffx.potential.bonded.Atom;
56  import ffx.potential.bonded.Molecule;
57  import ffx.potential.bonded.Residue;
58  import ffx.realspace.parsers.RealSpaceFile;
59  import ffx.xray.DataContainer;
60  import ffx.xray.DiffractionData;
61  import ffx.xray.RefinementMinimize.RefinementMode;
62  import ffx.xray.RefinementModel;
63  import java.util.List;
64  import java.util.logging.Level;
65  import java.util.logging.Logger;
66  import org.apache.commons.configuration2.CompositeConfiguration;
67  
68  /**
69   * RealSpaceData class.
70   *
71   * @author Timothy D. Fenn
72   * @since 1.0
73   */
74  public class RealSpaceData implements DataContainer {
75  
76    private static final Logger logger = Logger.getLogger(RealSpaceData.class.getName());
77    /** Parallel Team instance. */
78    private final ParallelTeam parallelTeam;
79    /** Calculate the real space target in parallel. */
80    private final RealSpaceRegion realSpaceRegion;
81    /** The real space data files. */
82    private final RealSpaceFile[] realSpaceFile;
83    /** The number of real space data sources to consider. */
84    private final int nRealSpaceData;
85    /** The refinement model. */
86    private final RefinementModel refinementModel;
87    /** A flag to indicate if the lambda term is active. */
88    private final boolean lambdaTerm;
89    /** The collection of MolecularAssembly instances that model the real space data. */
90    private final MolecularAssembly[] molecularAssemblies;
91    /** The collection of crystals that describe the PBC and symmetry of each data source. */
92    private Crystal[] crystal;
93    /** The real space refinement data for each data source. */
94    private RealSpaceRefinementData[] refinementData;
95    /** The total real space energy. */
96    private double realSpaceEnergy = 0.0;
97    /** The partial derivative of the real space energy with respect to lambda. */
98    private double realSpacedUdL = 0.0;
99    /** The real space gradient. */
100   private double[] realSpaceGradient;
101   /** The partial derivative of the real space gradient with respect to lambda. */
102   private double[] realSpacedUdXdL;
103   /** The weight of the data relative to the weight of the force field. */
104   private double xweight;
105   /** The current value of the state variable lambda. */
106   private double lambda = 1.0;
107 
108   /**
109    * Construct a real space data molecularAssemblies, assumes a real space map with a weight of 1.0
110    * using the same name as the molecular molecularAssemblies.
111    *
112    * @param molecularAssembly {@link ffx.potential.MolecularAssembly molecular molecularAssemblies}
113    *     object, used as the atomic model for comparison against the data
114    * @param properties system properties file
115    * @param parallelTeam a {@link edu.rit.pj.ParallelTeam} object.
116    * @param diffractionData {@link ffx.xray.DiffractionData diffraction data}
117    */
118   public RealSpaceData(
119       MolecularAssembly molecularAssembly,
120       CompositeConfiguration properties,
121       ParallelTeam parallelTeam,
122       DiffractionData diffractionData) {
123     this(new MolecularAssembly[] {molecularAssembly}, properties, parallelTeam, diffractionData);
124   }
125 
126   /**
127    * Construct a real space data molecularAssemblies, assumes a real space map with a weight of 1.0
128    * using the same name as the molecularAssemblies.
129    *
130    * @param molecularAssemblies {@link ffx.potential.MolecularAssembly molecular
131    *     molecularAssemblies} object, used as the atomic model for comparison against the data
132    * @param properties system properties file
133    * @param parallelTeam a {@link edu.rit.pj.ParallelTeam} object.
134    * @param diffractionData {@link ffx.xray.DiffractionData diffraction data}
135    */
136   public RealSpaceData(
137       MolecularAssembly[] molecularAssemblies,
138       CompositeConfiguration properties,
139       ParallelTeam parallelTeam,
140       DiffractionData diffractionData) {
141     this.molecularAssemblies = molecularAssemblies;
142     this.parallelTeam = parallelTeam;
143     this.realSpaceFile = null;
144     this.nRealSpaceData = 1;
145     crystal = new Crystal[nRealSpaceData];
146     crystal[0] = diffractionData.getCrystal()[0];
147     refinementData = new RealSpaceRefinementData[nRealSpaceData];
148     refinementData[0] = new RealSpaceRefinementData();
149     refinementData[0].setPeriodic(true);
150 
151     xweight = properties.getDouble("xweight", 1.0);
152     lambdaTerm = properties.getBoolean("lambdaterm", false);
153 
154     if (logger.isLoggable(Level.INFO)) {
155       StringBuilder sb = new StringBuilder(" Real Space Refinement Settings\n");
156       sb.append(format("  Refinement weight (xweight): %5.3f", xweight));
157       logger.info(sb.toString());
158     }
159 
160     if (!diffractionData.getScaled()[0]) {
161       diffractionData.scaleBulkFit();
162       diffractionData.printStats();
163     }
164 
165     // get Fo-Fc difference density
166     diffractionData.getCrystalReciprocalSpacesFc()[0].computeAtomicGradients(
167         diffractionData.getRefinementData()[0].foFc1,
168         diffractionData.getRefinementData()[0].freeR,
169         diffractionData.getRefinementData()[0].rFreeFlag,
170         RefinementMode.COORDINATES);
171 
172     refinementData[0].setOrigin(0, 0, 0);
173     int extx = (int) diffractionData.getCrystalReciprocalSpacesFc()[0].getXDim();
174     int exty = (int) diffractionData.getCrystalReciprocalSpacesFc()[0].getYDim();
175     int extz = (int) diffractionData.getCrystalReciprocalSpacesFc()[0].getZDim();
176     refinementData[0].setExtent(extx, exty, extz);
177     refinementData[0].setNI(extx, exty, extz);
178     refinementData[0].setData(new double[extx * exty * extz]);
179     for (int k = 0; k < extz; k++) {
180       for (int j = 0; j < exty; j++) {
181         for (int i = 0; i < extx; i++) {
182           int index1 = (i + extx * (j + exty * k));
183           int index2 = 2 * index1;
184           refinementData[0].getData()[index1] =
185               diffractionData.getCrystalReciprocalSpacesFc()[0].getDensityGrid()[index2];
186         }
187       }
188     }
189 
190     // Initialize the refinement model.
191     refinementModel = new RefinementModel(molecularAssemblies);
192 
193     // Initialize the RealSpaceRegion.
194     int nAtoms = refinementModel.getTotalAtomArray().length;
195     realSpaceRegion =
196         new RealSpaceRegion(parallelTeam.getThreadCount(), nAtoms, refinementData.length);
197   }
198 
199   /**
200    * Construct a real space data molecularAssemblies, assumes a real space map with a weight of 1.0
201    * using the same name as the molecular molecularAssemblies.
202    *
203    * @param assembly {@link ffx.potential.MolecularAssembly molecular molecularAssemblies} object,
204    *     used as the atomic model for comparison against the data
205    * @param properties system properties file
206    * @param parallelTeam a {@link edu.rit.pj.ParallelTeam} object.
207    */
208   public RealSpaceData(
209       MolecularAssembly assembly, CompositeConfiguration properties, ParallelTeam parallelTeam) {
210     this(new MolecularAssembly[] {assembly}, properties, parallelTeam, new RealSpaceFile(assembly));
211   }
212 
213   /**
214    * Construct a real space data molecularAssemblies.
215    *
216    * @param assembly {@link ffx.potential.MolecularAssembly molecular molecularAssemblies} object,
217    *     used as the atomic model for comparison against the data
218    * @param properties system properties file
219    * @param parallelTeam a {@link edu.rit.pj.ParallelTeam} object.
220    * @param datafile one or more {@link ffx.realspace.parsers.RealSpaceFile} to be refined against
221    */
222   public RealSpaceData(
223       MolecularAssembly assembly,
224       CompositeConfiguration properties,
225       ParallelTeam parallelTeam,
226       RealSpaceFile... datafile) {
227     this(new MolecularAssembly[] {assembly}, properties, parallelTeam, datafile);
228   }
229 
230   /**
231    * Construct a real space data molecularAssemblies.
232    *
233    * @param molecularAssemblies {@link ffx.potential.MolecularAssembly molecular
234    *     molecularAssemblies} object array (typically containing alternate conformer assemblies),
235    *     used as the atomic model for comparison against the data
236    * @param properties system properties file
237    * @param parallelTeam a {@link edu.rit.pj.ParallelTeam} object.
238    * @param dataFile one or more {@link ffx.realspace.parsers.RealSpaceFile} to be refined against
239    */
240   public RealSpaceData(
241       MolecularAssembly[] molecularAssemblies,
242       CompositeConfiguration properties,
243       ParallelTeam parallelTeam,
244       RealSpaceFile... dataFile) {
245 
246     this.molecularAssemblies = molecularAssemblies;
247     this.parallelTeam = parallelTeam;
248     this.realSpaceFile = dataFile;
249     this.nRealSpaceData = dataFile.length;
250     crystal = new Crystal[nRealSpaceData];
251     refinementData = new RealSpaceRefinementData[nRealSpaceData];
252 
253     xweight = properties.getDouble("xweight", 1.0);
254     lambdaTerm = properties.getBoolean("lambdaterm", false);
255 
256     for (int i = 0; i < nRealSpaceData; i++) {
257       crystal[i] =
258           dataFile[i].getRealSpaceFileFilter().getCrystal(dataFile[i].getFilename(), properties);
259       if (crystal[i] == null) {
260         logger.severe(" CCP4 map file does not contain full crystal information!");
261       }
262     }
263 
264     for (int i = 0; i < nRealSpaceData; i++) {
265       refinementData[i] = new RealSpaceRefinementData();
266       dataFile[i]
267           .getRealSpaceFileFilter()
268           .readFile(dataFile[i].getFilename(), refinementData[i], properties);
269 
270       if (refinementData[i].getOrigin()[0] == 0
271           && refinementData[i].getOrigin()[1] == 0
272           && refinementData[i].getOrigin()[2] == 0
273           && refinementData[i].getExtent()[0] == refinementData[i].getNi()[0]
274           && refinementData[i].getExtent()[1] == refinementData[i].getNi()[1]
275           && refinementData[i].getExtent()[2] == refinementData[i].getNi()[2]) {
276         refinementData[i].setPeriodic(true);
277       }
278     }
279 
280     if (logger.isLoggable(Level.INFO)) {
281       StringBuilder sb = new StringBuilder(" Real Space Refinement Settings\n");
282       sb.append(format("  Refinement weight (xweight): %5.3f", xweight));
283       logger.info(sb.toString());
284     }
285 
286     // now set up the refinement model
287     refinementModel = new RefinementModel(molecularAssemblies);
288 
289     // Initialize the RealSpaceRegion.
290     int nAtoms = refinementModel.getTotalAtomArray().length;
291     realSpaceRegion =
292         new RealSpaceRegion(parallelTeam.getThreadCount(), nAtoms, refinementData.length);
293   }
294 
295   /**
296    * Similar to Potential.destroy(), frees up resources associated with this RealSpaceData.
297    *
298    * @return If assets successfully freed.
299    */
300   public boolean destroy() {
301     try {
302       boolean underlyingShutdown = true;
303       for (MolecularAssembly assembly : molecularAssemblies) {
304         // Continue trying to shut assemblies down even if one fails to shut down.
305         boolean thisShutdown = assembly.destroy();
306         underlyingShutdown = underlyingShutdown && thisShutdown;
307       }
308       parallelTeam.shutdown();
309       return underlyingShutdown;
310     } catch (Exception ex) {
311       logger.warning(format(" Exception in shutting down a RealSpaceData: %s", ex));
312       logger.info(Utilities.stackTraceToString(ex));
313       return false;
314     }
315   }
316 
317   /** {@inheritDoc} */
318   @Override
319   public Atom[] getActiveAtomArray() {
320     return refinementModel.getActiveAtomArray();
321   }
322 
323   /** {@inheritDoc} */
324   @Override
325   public List<List<Molecule>> getAltMolecules() {
326     return refinementModel.getAltMolecules();
327   }
328 
329   /** {@inheritDoc} */
330   @Override
331   public List<List<Residue>> getAltResidues() {
332     return refinementModel.getAltResidues();
333   }
334 
335   /** {@inheritDoc} */
336   @Override
337   public Atom[] getAtomArray() {
338     return refinementModel.getTotalAtomArray();
339   }
340 
341   /**
342    * Getter for the field <code>crystal</code>.
343    *
344    * @return the crystal
345    */
346   public Crystal[] getCrystal() {
347     return crystal;
348   }
349 
350   /**
351    * Setter for the field <code>crystal</code>.
352    *
353    * @param crystal the crystal to set
354    */
355   public void setCrystal(Crystal[] crystal) {
356     this.crystal = crystal;
357   }
358 
359   /**
360    * Getter for the field <code>lambda</code>.
361    *
362    * @return the lambda
363    */
364   public double getLambda() {
365     return lambda;
366   }
367 
368   /**
369    * Set the current value of the state variable.
370    *
371    * @param lambda a double.
372    */
373   public void setLambda(double lambda) {
374     this.lambda = lambda;
375   }
376 
377   /** {@inheritDoc} */
378   @Override
379   public MolecularAssembly[] getMolecularAssemblies() {
380     return molecularAssemblies;
381   }
382 
383   /**
384    * Getter for the field <code>realSpaceGradient</code>.
385    *
386    * @param gradient an array of {@link double} objects.
387    * @return an array of {@link double} objects.
388    */
389   public double[] getRealSpaceGradient(double[] gradient) {
390     int nAtoms = refinementModel.getTotalAtomArray().length;
391     int nActiveAtoms = 0;
392     for (int i = 0; i < nAtoms; i++) {
393       if (refinementModel.getTotalAtomArray()[i].isActive()) {
394         nActiveAtoms++;
395       }
396     }
397 
398     if (gradient == null || gradient.length < nActiveAtoms * 3) {
399       gradient = new double[nActiveAtoms * 3];
400     }
401     for (int i = 0; i < nActiveAtoms * 3; i++) {
402       gradient[i] += realSpaceGradient[i];
403     }
404     return gradient;
405   }
406 
407   /**
408    * Getter for the field <code>refinementData</code>.
409    *
410    * @return the refinementData
411    */
412   public RealSpaceRefinementData[] getRefinementData() {
413     return refinementData;
414   }
415 
416   /**
417    * Setter for the field <code>refinementData</code>.
418    *
419    * @param refinementData the refinementData to set
420    */
421   public void setRefinementData(RealSpaceRefinementData[] refinementData) {
422     this.refinementData = refinementData;
423   }
424 
425   /** {@inheritDoc} */
426   @Override
427   public RefinementModel getRefinementModel() {
428     return refinementModel;
429   }
430 
431   /** {@inheritDoc} */
432   @Override
433   public double getWeight() {
434     return xweight;
435   }
436 
437   /** {@inheritDoc} */
438   @Override
439   public void setWeight(double weight) {
440     this.xweight = weight;
441   }
442 
443   /** {@inheritDoc} */
444   @Override
445   public String printEnergyUpdate() {
446     StringBuilder sb = new StringBuilder();
447     for (int i = 0; i < nRealSpaceData; i++) {
448       sb.append(
449           format(
450               "     dataset %d (weight: %5.1f): chemical energy: %8.2f density score: %8.2f\n",
451               i + 1,
452               realSpaceFile[i].getWeight(),
453               molecularAssemblies[0].getPotentialEnergy().getTotalEnergy(),
454               realSpaceFile[i].getWeight() * refinementData[i].getDensityScore()));
455     }
456     return sb.toString();
457   }
458 
459   /** {@inheritDoc} */
460   @Override
461   public String printOptimizationHeader() {
462     return "Density score";
463   }
464 
465   /** {@inheritDoc} */
466   @Override
467   public String printOptimizationUpdate() {
468     StringBuilder sb = new StringBuilder();
469     for (int i = 0; i < nRealSpaceData; i++) {
470       sb.append(format("%6.2f ", refinementData[i].getDensityScore()));
471     }
472     return sb.toString();
473   }
474 
475   /**
476    * compute real space target value, also fills in atomic derivatives
477    *
478    * @return target value sum over all data sets.
479    */
480   double computeRealSpaceTarget() {
481 
482     long time = -System.nanoTime();
483     // Zero out the realSpaceTarget energy.
484     realSpaceEnergy = 0.0;
485     // Zero out the realSpacedUdL energy.
486     realSpacedUdL = 0.0;
487     // Initialize gradient to zero; allocate space if necessary.
488     int nActive = 0;
489     int nAtoms = refinementModel.getTotalAtomArray().length;
490     for (int i = 0; i < nAtoms; i++) {
491       if (refinementModel.getTotalAtomArray()[i].isActive()) {
492         nActive++;
493       }
494     }
495 
496     int nGrad = nActive * 3;
497     if (realSpaceGradient == null || realSpaceGradient.length < nGrad) {
498       realSpaceGradient = new double[nGrad];
499     } else {
500       fill(realSpaceGradient, 0.0);
501     }
502 
503     // Initialize dUdXdL to zero; allocate space if necessary.
504     if (realSpacedUdXdL == null || realSpacedUdXdL.length < nGrad) {
505       realSpacedUdXdL = new double[nGrad];
506     } else {
507       fill(realSpacedUdXdL, 0.0);
508     }
509 
510     try {
511       parallelTeam.execute(realSpaceRegion);
512     } catch (Exception e) {
513       String message = " Exception computing real space energy";
514       logger.log(Level.SEVERE, message, e);
515     }
516 
517     time += System.nanoTime();
518     if (logger.isLoggable(Level.FINE)) {
519       logger.fine(format(" Real space energy time: %16.8f (sec).", time * 1.0e-9));
520     }
521 
522     return realSpaceEnergy;
523   }
524 
525   /**
526    * getdEdL.
527    *
528    * @return a double.
529    */
530   double getdEdL() {
531     return realSpacedUdL;
532   }
533 
534   /**
535    * getdEdXdL.
536    *
537    * @param gradient an array of {@link double} objects.
538    * @return an array of {@link double} objects.
539    */
540   double[] getdEdXdL(double[] gradient) {
541     int nAtoms = refinementModel.getTotalAtomArray().length;
542     int nActiveAtoms = 0;
543     for (int i = 0; i < nAtoms; i++) {
544       if (refinementModel.getTotalAtomArray()[i].isActive()) {
545         nActiveAtoms++;
546       }
547     }
548 
549     if (gradient == null || gradient.length < nActiveAtoms * 3) {
550       gradient = new double[nActiveAtoms * 3];
551     }
552 
553     for (int i = 0; i < nActiveAtoms * 3; i++) {
554       gradient[i] += realSpacedUdXdL[i];
555     }
556     return gradient;
557   }
558 
559   /**
560    * Getter for the field <code>realSpaceFile</code>.
561    *
562    * @return the realSpaceFile
563    */
564   private RealSpaceFile[] getRealSpaceFile() {
565     return realSpaceFile;
566   }
567 
568   /**
569    * Getter for the field <code>nRealSpaceData</code>.
570    *
571    * @return the nRealSpaceData
572    */
573   private int getnRealSpaceData() {
574     return nRealSpaceData;
575   }
576 
577   private class RealSpaceRegion extends ParallelRegion {
578 
579     private final AtomicDoubleArray3D gradient;
580     private final AtomicDoubleArray3D lambdaGrad;
581     private final InitializationLoop[] initializationLoops;
582     private final RealSpaceLoop[] realSpaceLoops;
583     private final SharedDouble[] sharedTarget;
584     private final SharedDouble shareddUdL;
585     private final int nAtoms;
586     private final int nData;
587 
588     RealSpaceRegion(int nThreads, int nAtoms, int nData) {
589       this.nAtoms = nAtoms;
590       this.nData = nData;
591       initializationLoops = new InitializationLoop[nThreads];
592       realSpaceLoops = new RealSpaceLoop[nThreads];
593       sharedTarget = new SharedDouble[nData];
594       for (int i = 0; i < nData; i++) {
595         sharedTarget[i] = new SharedDouble();
596       }
597       shareddUdL = new SharedDouble();
598       gradient = new AtomicDoubleArray3D(AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
599       lambdaGrad = new AtomicDoubleArray3D(AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
600     }
601 
602     @Override
603     public void finish() {
604       // Load final values.
605       realSpaceEnergy = 0;
606       for (int i = 0; i < nData; i++) {
607         refinementData[i].setDensityScore(sharedTarget[i].get());
608         realSpaceEnergy += getRefinementData()[i].getDensityScore();
609       }
610       realSpacedUdL = shareddUdL.get();
611       int index = 0;
612       for (int i = 0; i < nAtoms; i++) {
613         Atom atom = refinementModel.getTotalAtomArray()[i];
614         if (atom.isActive()) {
615           int ii = index * 3;
616           double gx = gradient.getX(i);
617           double gy = gradient.getY(i);
618           double gz = gradient.getZ(i);
619           realSpaceGradient[ii] = gx;
620           realSpaceGradient[ii + 1] = gy;
621           realSpaceGradient[ii + 2] = gz;
622           atom.setXYZGradient(gx, gy, gz);
623           gx = lambdaGrad.getX(i);
624           gy = lambdaGrad.getY(i);
625           gz = lambdaGrad.getZ(i);
626           realSpacedUdXdL[ii] = gx;
627           realSpacedUdXdL[ii + 1] = gy;
628           realSpacedUdXdL[ii + 2] = gz;
629           atom.setLambdaXYZGradient(gx, gy, gz);
630           index++;
631         }
632       }
633     }
634 
635     @Override
636     public void run() throws Exception {
637       int threadID = getThreadIndex();
638 
639       if (initializationLoops[threadID] == null) {
640         initializationLoops[threadID] = new InitializationLoop();
641       }
642       execute(0, nAtoms - 1, initializationLoops[threadID]);
643 
644       if (realSpaceLoops[threadID] == null) {
645         realSpaceLoops[threadID] = new RealSpaceLoop();
646       }
647       execute(0, nAtoms - 1, realSpaceLoops[threadID]);
648     }
649 
650     @Override
651     public void start() {
652       for (int i = 0; i < nData; i++) {
653         sharedTarget[i].set(0.0);
654       }
655       shareddUdL.set(0.0);
656       gradient.alloc(nAtoms);
657       lambdaGrad.alloc(nAtoms);
658     }
659 
660     /** Initialize gradient and lambda gradient arrays. */
661     private class InitializationLoop extends IntegerForLoop {
662 
663       @Override
664       public void run(int lb, int ub) {
665         int threadID = getThreadIndex();
666         gradient.reset(threadID, lb, ub);
667         lambdaGrad.reset(threadID, lb, ub);
668         for (int i = lb; i <= ub; i++) {
669           Atom a = refinementModel.getTotalAtomArray()[i];
670           a.setXYZGradient(0.0, 0.0, 0.0);
671           a.setLambdaXYZGradient(0.0, 0.0, 0.0);
672         }
673       }
674     }
675 
676     private class RealSpaceLoop extends IntegerForLoop {
677 
678       double[] target = new double[nData];
679       double localdUdL;
680 
681       @Override
682       public void finish() {
683         for (int i = 0; i < nData; i++) {
684           sharedTarget[i].addAndGet(target[i]);
685         }
686         shareddUdL.addAndGet(localdUdL);
687       }
688 
689       @Override
690       public void run(int first, int last) throws Exception {
691 
692         int threadID = getThreadIndex();
693         double[] xyz = new double[3];
694         double[] uvw = new double[3];
695         double[] grad = new double[3];
696         double[][][] scalar = new double[4][4][4];
697         TriCubicSpline spline = new TriCubicSpline();
698 
699         for (int i = 0; i < getnRealSpaceData(); i++) {
700 
701           // Define the extent of this real space data sources.
702           int extX = getRefinementData()[i].getExtent()[0];
703           int extY = getRefinementData()[i].getExtent()[1];
704           int extZ = getRefinementData()[i].getExtent()[2];
705           int nX = getRefinementData()[i].getNi()[0];
706           int nY = getRefinementData()[i].getNi()[1];
707           int nZ = getRefinementData()[i].getNi()[2];
708           int originX = getRefinementData()[i].getOrigin()[0];
709           int originY = getRefinementData()[i].getOrigin()[1];
710           int originZ = getRefinementData()[i].getOrigin()[2];
711 
712           for (int ia = first; ia <= last; ia++) {
713             Atom a = refinementModel.getTotalAtomArray()[ia];
714             // Only include atoms in the target function that have
715             // their use flag set to true and are Active.
716             if (!a.getUse()) {
717               continue;
718             }
719 
720             double lambdai = 1.0;
721             double dUdL = 0.0;
722             if (lambdaTerm && a.applyLambda()) {
723               lambdai = getLambda();
724               dUdL = 1.0;
725             }
726             a.getXYZ(xyz);
727             getCrystal()[i].toFractionalCoordinates(xyz, uvw);
728 
729             // Logic to find atom in 3d scalar field box.
730             final double frx = nX * uvw[0];
731             final int ifrx = ((int) floor(frx)) - originX;
732             final double dfrx = frx - floor(frx);
733 
734             final double fry = nY * uvw[1];
735             final int ifry = ((int) floor(fry)) - originY;
736             final double dfry = fry - floor(fry);
737 
738             final double frz = nZ * uvw[2];
739             final int ifrz = ((int) floor(frz)) - originZ;
740             final double dfrz = frz - floor(frz);
741 
742             if (!refinementData[i].isPeriodic()) {
743               if (ifrx - 1 < 0
744                   || ifrx + 2 > extX
745                   || ifry - 1 < 0
746                   || ifry + 2 > extY
747                   || ifrz - 1 < 0
748                   || ifrz + 2 > extZ) {
749                 String message = format(" Atom %s is outside the density will be ignored.", a);
750                 logger.warning(message);
751                 continue;
752               }
753             }
754 
755             // Fill in scalar 4x4 array for interpolation.
756             if (getRefinementData()[i].isPeriodic()) {
757               for (int ui = ifrx - 1; ui < ifrx + 3; ui++) {
758                 int uii = ui - (ifrx - 1);
759                 int pui = mod(ui, extX);
760                 for (int vi = ifry - 1; vi < ifry + 3; vi++) {
761                   int vii = vi - (ifry - 1);
762                   int pvi = mod(vi, extY);
763                   for (int wi = ifrz - 1; wi < ifrz + 3; wi++) {
764                     int wii = wi - (ifrz - 1);
765                     int pwi = mod(wi, extZ);
766                     scalar[uii][vii][wii] = getRefinementData()[i].getDataIndex(pui, pvi, pwi);
767                   }
768                 }
769               }
770             } else {
771               for (int ui = ifrx - 1; ui < ifrx + 3; ui++) {
772                 int uii = ui - (ifrx - 1);
773                 for (int vi = ifry - 1; vi < ifry + 3; vi++) {
774                   int vii = vi - (ifry - 1);
775                   for (int wi = ifrz - 1; wi < ifrz + 3; wi++) {
776                     int wii = wi - (ifrz - 1);
777                     scalar[uii][vii][wii] = getRefinementData()[i].getDataIndex(ui, vi, wi);
778                   }
779                 }
780               }
781             }
782 
783             // Scale and interpolate.
784             double scale;
785             double scaledUdL;
786             double atomicWeight = a.getAtomType().atomicWeight;
787             if (getRealSpaceFile() != null) {
788               atomicWeight *= getRealSpaceFile()[i].getWeight();
789             }
790             scale = -1.0 * lambdai * atomicWeight;
791             scaledUdL = -1.0 * dUdL * atomicWeight;
792 
793             double val = spline.spline(dfrx, dfry, dfrz, scalar, grad);
794             target[i] += scale * val;
795             localdUdL += scaledUdL * val;
796 
797             if (a.isActive()) {
798               grad[0] = grad[0] * nX;
799               grad[1] = grad[1] * nY;
800               grad[2] = grad[2] * nZ;
801               // transpose of toFractional
802               xyz[0] =
803                   grad[0] * getCrystal()[i].A00
804                       + grad[1] * getCrystal()[i].A01
805                       + grad[2] * getCrystal()[i].A02;
806               xyz[1] =
807                   grad[0] * getCrystal()[i].A10
808                       + grad[1] * getCrystal()[i].A11
809                       + grad[2] * getCrystal()[i].A12;
810               xyz[2] =
811                   grad[0] * getCrystal()[i].A20
812                       + grad[1] * getCrystal()[i].A21
813                       + grad[2] * getCrystal()[i].A22;
814               gradient.add(threadID, ia, scale * xyz[0], scale * xyz[1], scale * xyz[2]);
815               lambdaGrad.add(
816                   threadID, ia, scaledUdL * xyz[0], scaledUdL * xyz[1], scaledUdL * xyz[2]);
817             }
818           }
819         }
820         gradient.reduce(first, last);
821         lambdaGrad.reduce(first, last);
822       }
823 
824       @Override
825       public void start() {
826         for (int i = 0; i < nData; i++) {
827           target[i] = 0;
828         }
829         localdUdL = 0;
830       }
831     }
832   }
833 }