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.xray;
39  
40  import static ffx.utilities.TinkerUtils.version;
41  import static ffx.xray.CrystalReciprocalSpace.SolventModel.POLYNOMIAL;
42  import static java.lang.String.format;
43  import static java.util.Arrays.fill;
44  import static org.apache.commons.io.FilenameUtils.removeExtension;
45  
46  import edu.rit.pj.ParallelTeam;
47  import ffx.xray.parsers.CCP4MapWriter;
48  import ffx.crystal.Crystal;
49  import ffx.crystal.HKL;
50  import ffx.crystal.ReflectionList;
51  import ffx.crystal.Resolution;
52  import ffx.potential.MolecularAssembly;
53  import ffx.potential.Utilities;
54  import ffx.potential.bonded.Atom;
55  import ffx.potential.bonded.Molecule;
56  import ffx.potential.bonded.Residue;
57  import ffx.potential.parameters.ForceField;
58  import ffx.potential.parsers.PDBFilter;
59  import ffx.xray.CrystalReciprocalSpace.SolventModel;
60  import ffx.xray.RefinementMinimize.RefinementMode;
61  import ffx.xray.parsers.DiffractionFile;
62  import ffx.xray.parsers.MTZWriter;
63  import ffx.xray.parsers.MTZWriter.MTZType;
64  import java.io.BufferedWriter;
65  import java.io.File;
66  import java.io.FileWriter;
67  import java.io.PrintWriter;
68  import java.text.SimpleDateFormat;
69  import java.util.Arrays;
70  import java.util.Date;
71  import java.util.List;
72  import java.util.Set;
73  import java.util.logging.Level;
74  import java.util.logging.Logger;
75  import org.apache.commons.configuration2.CompositeConfiguration;
76  
77  /**
78   * DiffractionData class.
79   *
80   * @author Timothy D. Fenn
81   * @since 1.0
82   */
83  public class DiffractionData implements DataContainer {
84  
85    private static final Logger logger = Logger.getLogger(DiffractionData.class.getName());
86    private final MolecularAssembly[] assembly;
87    private final String modelName;
88    private final DiffractionFile[] dataFiles;
89    private final int n;
90    private final Crystal[] crystal;
91    private final Resolution[] resolution;
92    private final ReflectionList[] reflectionList;
93    private final DiffractionRefinementData[] refinementData;
94    private final CrystalReciprocalSpace[] crystalReciprocalSpacesFc;
95    private final CrystalReciprocalSpace[] crystalReciprocalSpacesFs;
96    private final SolventModel solventModel;
97    private final RefinementModel refinementModel;
98    // Settings
99    private final double fsigfCutoff;
100   private final boolean use_3g;
101   private final double aRadBuff;
102   private final double xrayScaleTol;
103   private final double sigmaATol;
104   private final double bSimWeight;
105   private final double bNonZeroWeight;
106 
107 
108   private final double bMass;
109   private final boolean residueBFactor;
110   private final int nResidueBFactor;
111   private final boolean addAnisou;
112   private final boolean refineMolOcc;
113   private final double occMass;
114   private final boolean nativeEnvironmentApproximation;
115   private final ScaleBulkMinimize[] scaleBulkMinimize;
116   private final SigmaAMinimize[] sigmaAMinimize;
117   private final SplineMinimize[] splineMinimize;
118   private final CrystalStats[] crystalStats;
119   private ParallelTeam parallelTeam;
120   private CrystalReciprocalSpace.GridMethod gridMethod;
121   private final boolean[] scaled;
122   private double xWeight;
123   /** If true, perform a grid search for bulk solvent parameters. */
124   private final boolean gridSearch;
125   /** If true, fit a scaling spline between Fo and Fc. */
126   private final boolean splineFit;
127 
128   /**
129    * construct a diffraction data assembly, assumes an X-ray data set with a weight of 1.0 using the
130    * same name as the molecular assembly
131    *
132    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object, used as
133    *     the atomic model for comparison against the data
134    * @param properties system properties file
135    */
136   public DiffractionData(MolecularAssembly assembly, CompositeConfiguration properties) {
137     this(new MolecularAssembly[] {assembly}, properties, POLYNOMIAL, new DiffractionFile(assembly));
138   }
139 
140   /**
141    * construct a diffraction data assembly
142    *
143    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object, used as
144    *     the atomic model for comparison against the data
145    * @param properties system properties file
146    * @param datafile one or more {@link ffx.xray.parsers.DiffractionFile} to be refined against
147    */
148   public DiffractionData(
149       MolecularAssembly assembly, CompositeConfiguration properties, DiffractionFile... datafile) {
150     this(new MolecularAssembly[] {assembly}, properties, POLYNOMIAL, datafile);
151   }
152 
153   /**
154    * construct a diffraction data assembly, assumes an X-ray data set with a weight of 1.0 using the
155    * same name as the molecular assembly
156    *
157    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object, used as
158    *     the atomic model for comparison against the data
159    * @param properties system properties file
160    * @param solventmodel the type of solvent model desired - see {@link
161    *     CrystalReciprocalSpace.SolventModel bulk solvent model} selections
162    */
163   public DiffractionData(
164       MolecularAssembly assembly, CompositeConfiguration properties, SolventModel solventmodel) {
165     this(
166         new MolecularAssembly[] {assembly},
167         properties,
168         solventmodel,
169         new DiffractionFile(assembly));
170   }
171 
172   /**
173    * construct a diffraction data assembly
174    *
175    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object, used as
176    *     the atomic model for comparison against the data
177    * @param properties system properties file
178    * @param solventmodel the type of solvent model desired - see {@link
179    *     CrystalReciprocalSpace.SolventModel bulk solvent model} selections
180    * @param datafile one or more {@link ffx.xray.parsers.DiffractionFile} to be refined against
181    */
182   public DiffractionData(
183       MolecularAssembly assembly,
184       CompositeConfiguration properties,
185       SolventModel solventmodel,
186       DiffractionFile... datafile) {
187     this(new MolecularAssembly[] {assembly}, properties, solventmodel, datafile);
188   }
189 
190   /**
191    * construct a diffraction data assembly, assumes an X-ray data set with a weight of 1.0 using the
192    * same name as the molecular assembly
193    *
194    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object array
195    *     (typically containing alternate conformer assemblies), used as the atomic model for
196    *     comparison against the data
197    * @param properties system properties file
198    */
199   public DiffractionData(MolecularAssembly[] assembly, CompositeConfiguration properties) {
200     this(assembly, properties, POLYNOMIAL, new DiffractionFile(assembly[0]));
201   }
202 
203   /**
204    * construct a diffraction data assembly
205    *
206    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object array
207    *     (typically containing alternate conformer assemblies), used as the atomic model for
208    *     comparison against the data
209    * @param properties system properties file
210    * @param datafile one or more {@link ffx.xray.parsers.DiffractionFile} to be refined against
211    */
212   public DiffractionData(
213       MolecularAssembly[] assembly,
214       CompositeConfiguration properties,
215       DiffractionFile... datafile) {
216     this(assembly, properties, POLYNOMIAL, datafile);
217   }
218 
219   /**
220    * construct a diffraction data assembly
221    *
222    * @param assembly {@link ffx.potential.MolecularAssembly molecular assembly} object array
223    *     (typically containing alternate conformer assemblies), used as the atomic model for
224    *     comparison against the data
225    * @param properties system properties file
226    * @param solventmodel the type of solvent model desired - see {@link
227    *     CrystalReciprocalSpace.SolventModel bulk solvent model} selections
228    * @param datafile one or more {@link ffx.xray.parsers.DiffractionFile} to be refined against
229    */
230   public DiffractionData(
231       MolecularAssembly[] assembly,
232       CompositeConfiguration properties,
233       SolventModel solventmodel,
234       DiffractionFile... datafile) {
235 
236     this.assembly = assembly;
237     this.solventModel = solventmodel;
238     this.modelName = assembly[0].getFile().getName();
239     this.dataFiles = datafile;
240     this.n = datafile.length;
241 
242     int rflag = properties.getInt("rfree-flag", -1);
243     fsigfCutoff = properties.getDouble("f-sigf-cutoff", -1.0);
244     gridSearch = properties.getBoolean("solvent-grid-search", false);
245     splineFit = !properties.getBoolean("no-spline-fit", false);
246     use_3g = properties.getBoolean("use-3g", true);
247     aRadBuff = properties.getDouble("scattering-buffer", 0.75);
248     double sampling = properties.getDouble("sampling", 0.6);
249     xrayScaleTol = properties.getDouble("xray-scale-tol", 1e-4);
250     sigmaATol = properties.getDouble("sigmaa-tol", 0.05);
251     xWeight = properties.getDouble("data-weight", 1.0);
252     bSimWeight = properties.getDouble("b-sim-weight", 1.0);
253     bNonZeroWeight = properties.getDouble("b-nonzero-weight", 1.0);
254     bMass = properties.getDouble("bfactor-mass", 5.0);
255     residueBFactor = properties.getBoolean("residue-bfactor", false);
256     nResidueBFactor = properties.getInt("n-residue-bfactor", 1);
257     addAnisou = properties.getBoolean("add-anisou", false);
258     refineMolOcc = properties.getBoolean("refine-mol-occ", false);
259     occMass = properties.getDouble("occmass", 10.0);
260 
261     ForceField forceField = assembly[0].getForceField();
262     nativeEnvironmentApproximation =
263         forceField.getBoolean("NATIVE_ENVIRONMENT_APPROXIMATION", false);
264 
265     crystal = new Crystal[n];
266     resolution = new Resolution[n];
267     reflectionList = new ReflectionList[n];
268     refinementData = new DiffractionRefinementData[n];
269     scaleBulkMinimize = new ScaleBulkMinimize[n];
270     sigmaAMinimize = new SigmaAMinimize[n];
271     splineMinimize = new SplineMinimize[n];
272     crystalStats = new CrystalStats[n];
273 
274     // Read in Fo/sigFo/FreeR
275     File tmp;
276     Crystal crystalinit = Crystal.checkProperties(properties);
277     Resolution resolutioninit = Resolution.checkProperties(properties);
278     if (crystalinit == null || resolutioninit == null) {
279       for (int i = 0; i < n; i++) {
280         tmp = new File(datafile[i].getFilename());
281         reflectionList[i] = datafile[i].getDiffractionfilter().getReflectionList(tmp, properties);
282         if (reflectionList[i] == null) {
283           logger.info(" Crystal information from the PDB or property file will be used.");
284           crystalinit = assembly[i].getCrystal().getUnitCell();
285           double res = datafile[i].getDiffractionfilter().getResolution(tmp, crystalinit);
286           if (res < 0.0) {
287             logger.severe("MTZ/CIF/CNS file does not contain full crystal information!");
288           } else {
289             resolutioninit = new Resolution(res, sampling);
290             reflectionList[i] = new ReflectionList(crystalinit, resolutioninit, properties);
291           }
292         }
293       }
294     } else {
295       for (int i = 0; i < n; i++) {
296         reflectionList[i] = new ReflectionList(crystalinit, resolutioninit, properties);
297       }
298     }
299 
300     for (int i = 0; i < n; i++) {
301       crystal[i] = reflectionList[i].crystal;
302       resolution[i] = reflectionList[i].resolution;
303       refinementData[i] = new DiffractionRefinementData(properties, reflectionList[i]);
304       tmp = new File(datafile[i].getFilename());
305       datafile[i].getDiffractionfilter()
306           .readFile(tmp, reflectionList[i], refinementData[i], properties);
307     }
308 
309     if (!crystal[0].getUnitCell().equals(assembly[0].getCrystal().getUnitCell())) {
310       logger.info("\n The PDB and reflection file crystal information do not match.");
311       logger.info(" PDB File:" + assembly[0].getCrystal().getUnitCell().toString());
312       logger.info(" Reflection File:" + crystal[0].getUnitCell().toString());
313       logger.severe(" Please check the concordance of the PDB CRYST1 record with the diffraction file.");
314     }
315 
316     if (logger.isLoggable(Level.INFO)) {
317       StringBuilder sb = new StringBuilder();
318       sb.append("\n X-ray Refinement Settings\n\n");
319       sb.append("  Target Function\n");
320       sb.append("   X-ray refinement weight: ").append(xWeight).append("\n");
321       sb.append("   Use cctbx 3 Gaussians: ").append(use_3g).append("\n");
322       sb.append("   Atomic form factor radius buffer: ").append(aRadBuff).append("\n");
323       sb.append("   Reciprocal space sampling rate: ").append(sampling).append("\n");
324       sb.append("   Resolution dependent spline scale: ").append(splineFit).append("\n");
325       sb.append("   Solvent grid search: ").append(gridSearch).append("\n");
326       sb.append("   X-ray scale fit tolerance: ").append(xrayScaleTol).append("\n");
327       sb.append("   Sigma A fit tolerance: ").append(sigmaATol).append("\n");
328       sb.append("   Native environment approximation: ").append(nativeEnvironmentApproximation)
329           .append("\n");
330       sb.append("  Reflections\n");
331       sb.append("   F/sigF cutoff: ").append(fsigfCutoff).append("\n");
332       sb.append("   R Free flag (-1 auto-determine from the data): ").append(rflag).append("\n");
333       sb.append("   Number of bins: ").append(reflectionList[0].nBins).append("\n");
334       sb.append("  B-Factors\n");
335       sb.append("   Similarity weight: ").append(bSimWeight).append("\n");
336       // sb.append("  Non-zero weight (bnonzeroweight): ").append(bNonZeroWeight).append("\n");
337       // sb.append("  Lagrangian mass (bmass): ").append(bMass).append("\n");
338       sb.append("   Refine by residue: ").append(residueBFactor).append("\n");
339       if (residueBFactor) {
340         sb.append("   Number of residues per B: ").append(nResidueBFactor).append(")\n");
341       }
342       sb.append("   Add ANISOU for refinement: ").append(addAnisou).append("\n");
343       sb.append("  Occupancies\n");
344       sb.append("   Refine on molecules: ").append(refineMolOcc).append("\n");
345       // sb.append("  Lagrangian mass (occmass): ").append(occMass).append("\n");
346 
347       logger.info(sb.toString());
348     }
349 
350     // now set up the refinement model
351     refinementModel = new RefinementModel(assembly, refineMolOcc);
352 
353     // initialize atomic form factors
354     for (Atom a : refinementModel.getTotalAtomArray()) {
355       a.setFormFactorIndex(-1);
356       XRayFormFactor atomff = new XRayFormFactor(a, use_3g, 2.0);
357       a.setFormFactorIndex(atomff.ffIndex);
358 
359       if (a.getOccupancy() == 0.0) {
360         a.setFormFactorWidth(1.0);
361         continue;
362       }
363 
364       double arad;
365       try {
366         arad = a.getVDWType().radius * 0.5;
367       } catch (NullPointerException ex) {
368         logger.warning(format(" Failure to get van der Waals type for atom %s; ensure the vdW term is enabled!", a));
369         throw ex;
370       }
371       double[] xyz = new double[3];
372       xyz[0] = a.getX() + arad;
373       xyz[1] = a.getY();
374       xyz[2] = a.getZ();
375       double rho = atomff.rho(0.0, 1.0, xyz);
376       while (rho > 0.001) {
377         arad += 0.1;
378         xyz[0] = a.getX() + arad;
379         rho = atomff.rho(0.0, 1.0, xyz);
380       }
381       arad += aRadBuff;
382       a.setFormFactorWidth(arad);
383     }
384 
385     // set up FFT and run it
386     crystalReciprocalSpacesFc = new CrystalReciprocalSpace[n];
387     crystalReciprocalSpacesFs = new CrystalReciprocalSpace[n];
388 
389     parallelTeam = assembly[0].getParallelTeam();
390 
391     String gridString = properties.getString("grid-method", "SLICE").toUpperCase();
392     try {
393       gridMethod = CrystalReciprocalSpace.GridMethod.valueOf(gridString);
394     } catch (Exception e) {
395       gridMethod = CrystalReciprocalSpace.GridMethod.SLICE;
396     }
397 
398     parallelTeam = new ParallelTeam();
399     for (int i = 0; i < n; i++) {
400       // Atomic Scattering
401       crystalReciprocalSpacesFc[i] =
402           new CrystalReciprocalSpace(
403               reflectionList[i],
404               refinementModel.getTotalAtomArray(),
405               parallelTeam,
406               parallelTeam,
407               false,
408               dataFiles[i].isNeutron(),
409               SolventModel.NONE,
410               gridMethod);
411       refinementData[i].setCrystalReciprocalSpaceFc(crystalReciprocalSpacesFc[i]);
412       crystalReciprocalSpacesFc[i].setUse3G(use_3g);
413       crystalReciprocalSpacesFc[i].setWeight(dataFiles[i].getWeight());
414       crystalReciprocalSpacesFc[i].lambdaTerm = false;
415       crystalReciprocalSpacesFc[i].setNativeEnvironmentApproximation(
416           nativeEnvironmentApproximation);
417 
418       // Bulk Solvent Scattering
419       crystalReciprocalSpacesFs[i] =
420           new CrystalReciprocalSpace(
421               reflectionList[i],
422               refinementModel.getTotalAtomArray(),
423               parallelTeam,
424               parallelTeam,
425               true,
426               dataFiles[i].isNeutron(),
427               solventmodel,
428               gridMethod);
429       refinementData[i].setCrystalReciprocalSpaceFs(crystalReciprocalSpacesFs[i]);
430       crystalReciprocalSpacesFs[i].setUse3G(use_3g);
431       crystalReciprocalSpacesFs[i].setWeight(dataFiles[i].getWeight());
432       crystalReciprocalSpacesFs[i].lambdaTerm = false;
433       crystalReciprocalSpacesFs[i].setNativeEnvironmentApproximation(
434           nativeEnvironmentApproximation);
435       crystalStats[i] = new CrystalStats(reflectionList[i], refinementData[i]);
436     }
437 
438     scaled = new boolean[n];
439     fill(scaled, false);
440   }
441 
442   /**
443    * read in a different assembly to average in structure factors
444    *
445    * @param assembly the {@link ffx.potential.MolecularAssembly molecular assembly} object array
446    *     (typically containing alternate conformer assemblies), used as the atomic model to average
447    *     in with previous assembly
448    * @param index the current data index (for cumulative average purposes)
449    */
450   public void AverageFc(MolecularAssembly[] assembly, int index) {
451     RefinementModel tmprefinementmodel = new RefinementModel(assembly, refineMolOcc);
452 
453     // initialize atomic form factors
454     for (Atom a : tmprefinementmodel.getTotalAtomArray()) {
455       a.setFormFactorIndex(-1);
456       XRayFormFactor atomff = new XRayFormFactor(a, use_3g, 2.0);
457       a.setFormFactorIndex(atomff.ffIndex);
458 
459       if (a.getOccupancy() == 0.0) {
460         a.setFormFactorWidth(1.0);
461         continue;
462       }
463 
464       double arad = a.getVDWType().radius * 0.5;
465       double[] xyz = new double[3];
466       xyz[0] = a.getX() + arad;
467       xyz[1] = a.getY();
468       xyz[2] = a.getZ();
469       double rho = atomff.rho(0.0, 1.0, xyz);
470       while (rho > 0.001) {
471         arad += 0.1;
472         xyz[0] = a.getX() + arad;
473         rho = atomff.rho(0.0, 1.0, xyz);
474       }
475       arad += aRadBuff;
476       a.setFormFactorWidth(arad);
477     }
478 
479     // set up FFT and run it
480     for (int i = 0; i < n; i++) {
481       crystalReciprocalSpacesFc[i] =
482           new CrystalReciprocalSpace(
483               reflectionList[i],
484               tmprefinementmodel.getTotalAtomArray(),
485               parallelTeam,
486               parallelTeam,
487               false,
488               dataFiles[i].isNeutron(),
489               SolventModel.NONE,
490               gridMethod);
491       crystalReciprocalSpacesFc[i].setNativeEnvironmentApproximation(
492           nativeEnvironmentApproximation);
493       refinementData[i].setCrystalReciprocalSpaceFc(crystalReciprocalSpacesFc[i]);
494 
495       crystalReciprocalSpacesFs[i] =
496           new CrystalReciprocalSpace(
497               reflectionList[i],
498               tmprefinementmodel.getTotalAtomArray(),
499               parallelTeam,
500               parallelTeam,
501               true,
502               dataFiles[i].isNeutron(),
503               solventModel,
504               gridMethod);
505       crystalReciprocalSpacesFs[i].setNativeEnvironmentApproximation(
506           nativeEnvironmentApproximation);
507       refinementData[i].setCrystalReciprocalSpaceFs(crystalReciprocalSpacesFs[i]);
508     }
509 
510     int nhkl = refinementData[0].n;
511     double[][] fc = new double[nhkl][2];
512     double[][] fs = new double[nhkl][2];
513 
514     // Run FFTs.
515     for (int i = 0; i < n; i++) {
516       crystalReciprocalSpacesFc[i].computeDensity(fc);
517       if (solventModel != SolventModel.NONE) {
518         crystalReciprocalSpacesFs[i].computeDensity(fs);
519       }
520 
521       // Average in with current data.
522       for (int j = 0; j < refinementData[i].n; j++) {
523         refinementData[i].fc[j][0] += (fc[j][0] - refinementData[i].fc[j][0]) / index;
524         refinementData[i].fc[j][1] += (fc[j][1] - refinementData[i].fc[j][1]) / index;
525 
526         refinementData[i].fs[j][0] += (fs[j][0] - refinementData[i].fs[j][0]) / index;
527         refinementData[i].fs[j][1] += (fs[j][1] - refinementData[i].fs[j][1]) / index;
528       }
529     }
530   }
531 
532   /**
533    * Parallelized call to compute atomic density on a grid, followed by FFT to compute structure
534    * factors.
535    *
536    * @see CrystalReciprocalSpace#computeDensity(double[][], boolean)
537    */
538   public void computeAtomicDensity() {
539     for (int i = 0; i < n; i++) {
540       crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc);
541       if (solventModel != SolventModel.NONE) {
542         crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs);
543       }
544     }
545   }
546 
547   /**
548    * Similar to Potential.destroy(), frees up resources associated with this RealSpaceData.
549    *
550    * @return If assets successfully freed.
551    */
552   public boolean destroy() {
553     try {
554       boolean underlyingShutdown = true;
555       for (MolecularAssembly assem : assembly) {
556         // Continue trying to shut assemblies down even if one fails to shut down.
557         boolean thisShutdown = assem.destroy();
558         underlyingShutdown = underlyingShutdown && thisShutdown;
559       }
560       parallelTeam.shutdown();
561       return underlyingShutdown;
562     } catch (Exception ex) {
563       logger.warning(format(" Exception in shutting down a RealSpaceData: %s", ex));
564       logger.info(Utilities.stackTraceToString(ex));
565       return false;
566     }
567   }
568 
569   /**
570    * getActiveAtomArray.
571    *
572    * @return an array of {@link ffx.potential.bonded.Atom} objects.
573    */
574   public Atom[] getActiveAtomArray() {
575     return getAtomArray();
576   }
577 
578   /** {@inheritDoc} */
579   @Override
580   public List<List<Molecule>> getAltMolecules() {
581     return refinementModel.getAltMolecules();
582   }
583 
584   /** {@inheritDoc} */
585   @Override
586   public List<List<Residue>> getAltResidues() {
587     return refinementModel.getAltResidues();
588   }
589 
590   /**
591    * Getter for the field <code>assembly</code>.
592    *
593    * @return the assembly
594    */
595   public MolecularAssembly[] getAssembly() {
596     return assembly;
597   }
598 
599   /**
600    * {@inheritDoc}
601    *
602    * <p>return the atom array for the model associated with this data
603    */
604   @Override
605   public Atom[] getAtomArray() {
606     return refinementModel.getTotalAtomArray();
607   }
608 
609   /**
610    * Getter for the field <code>crystal</code>.
611    *
612    * @return the crystal
613    */
614   public Crystal[] getCrystal() {
615     return crystal;
616   }
617 
618   /**
619    * Getter for the field <code>crs_fc</code>.
620    *
621    * @return the crs_fc
622    */
623   public CrystalReciprocalSpace[] getCrystalReciprocalSpacesFc() {
624     return crystalReciprocalSpacesFc;
625   }
626 
627   /**
628    * Getter for the field <code>crs_fs</code>.
629    *
630    * @return the crs_fs
631    */
632   public CrystalReciprocalSpace[] getCrystalReciprocalSpacesFs() {
633     return crystalReciprocalSpacesFs;
634   }
635 
636   /**
637    * Getter for the field <code>crystalStats</code>.
638    *
639    * @return the crystalStats
640    */
641   public CrystalStats[] getCrystalStats() {
642     return crystalStats;
643   }
644 
645   /**
646    * Getter for the field <code>dataFiles</code>.
647    *
648    * @return the dataFiles
649    */
650   public DiffractionFile[] getDataFiles() {
651     return dataFiles;
652   }
653 
654   /**
655    * Getter for the field <code>fsigfCutoff</code>.
656    *
657    * @return the fsigfCutoff
658    */
659   public double getFsigfCutoff() {
660     return fsigfCutoff;
661   }
662 
663   /**
664    * Getter for the field <code>modelName</code>.
665    *
666    * @return the modelName
667    */
668   public String getModelName() {
669     return modelName;
670   }
671 
672   /** {@inheritDoc} */
673   @Override
674   public MolecularAssembly[] getMolecularAssemblies() {
675     return assembly;
676   }
677 
678   /**
679    * Getter for the field <code>n</code>.
680    *
681    * @return the n
682    */
683   public int getN() {
684     return n;
685   }
686 
687   /**
688    * Getter for the field <code>parallelTeam</code>.
689    *
690    * @return the parallelTeam
691    */
692   public ParallelTeam getParallelTeam() {
693     return parallelTeam;
694   }
695 
696   /**
697    * Return R value for OST x-ray minimization.
698    *
699    * @return a double.
700    */
701   public double getRCrystalStat() {
702     return crystalStats[0].getR();
703   }
704 
705   /**
706    * Getter for the field <code>refinementData</code>.
707    *
708    * @return the refinementData
709    */
710   public DiffractionRefinementData[] getRefinementData() {
711     return refinementData;
712   }
713 
714   /** {@inheritDoc} */
715   @Override
716   public RefinementModel getRefinementModel() {
717     return refinementModel;
718   }
719 
720   /**
721    * Getter for the field <code>reflectionList</code>.
722    *
723    * @return the reflectionList
724    */
725   public ReflectionList[] getReflectionList() {
726     return reflectionList;
727   }
728 
729   /**
730    * Getter for the field <code>resolution</code>.
731    *
732    * @return the resolution
733    */
734   public Resolution[] getResolution() {
735     return resolution;
736   }
737 
738   /**
739    * Getter for the field <code>scaleBulkMinimize</code>.
740    *
741    * @return the scaleBulkMinimize
742    */
743   public ScaleBulkMinimize[] getScaleBulkMinimize() {
744     return scaleBulkMinimize;
745   }
746 
747   /**
748    * Getter for the field <code>scaled</code>.
749    *
750    * @return the scaled
751    */
752   public boolean[] getScaled() {
753     return scaled;
754   }
755 
756   /**
757    * Getter for the field <code>sigmaAMinimize</code>.
758    *
759    * @return the sigmaAMinimize
760    */
761   public SigmaAMinimize[] getSigmaAMinimize() {
762     return sigmaAMinimize;
763   }
764 
765   /**
766    * Getter for the field <code>sigmaATol</code>.
767    *
768    * @return the sigmaATol
769    */
770   public double getSigmaATol() {
771     return sigmaATol;
772   }
773 
774   /**
775    * Getter for the field <code>solventModel</code>.
776    *
777    * @return the solventModel
778    */
779   public SolventModel getSolventModel() {
780     return solventModel;
781   }
782 
783   /**
784    * Getter for the field <code>splineMinimize</code>.
785    *
786    * @return the splineMinimize
787    */
788   public SplineMinimize[] getSplineMinimize() {
789     return splineMinimize;
790   }
791 
792   /** {@inheritDoc} */
793   @Override
794   public double getWeight() {
795     return xWeight;
796   }
797 
798   /** {@inheritDoc} */
799   @Override
800   public void setWeight(double weight) {
801     this.xWeight = weight;
802   }
803 
804   /**
805    * Getter for the field <code>xrayScaleTol</code>.
806    *
807    * @return the xrayScaleTol
808    */
809   public double getXrayScaleTol() {
810     return xrayScaleTol;
811   }
812 
813   /**
814    * Getter for the field <code>aRadBuff</code>.
815    *
816    * @return the aRadBuff
817    */
818   public double getaRadBuff() {
819     return aRadBuff;
820   }
821 
822   /**
823    * Getter for the field <code>xWeight</code>.
824    *
825    * @return the xWeight
826    */
827   public double getxWeight() {
828     return xWeight;
829   }
830 
831   /**
832    * isRefineMolOcc.
833    *
834    * @return the refineMolOcc
835    */
836   public boolean isRefineMolOcc() {
837     return refineMolOcc;
838   }
839 
840   /**
841    * isUse_3g.
842    *
843    * @return the use_3g
844    */
845   public boolean isUse_3g() {
846     return use_3g;
847   }
848 
849   /** {@inheritDoc} */
850   @Override
851   public String printEnergyUpdate() {
852     StringBuilder sb = new StringBuilder();
853     for (int i = 0; i < n; i++) {
854       sb.append(
855           format(
856               "     dataset %d (weight: %5.1f): R: %6.2f Rfree: %6.2f chemical energy: %8.2f likelihood: %8.2f free likelihood: %8.2f\n",
857               i + 1,
858               dataFiles[i].getWeight(),
859               crystalStats[i].getR(),
860               crystalStats[i].getRFree(),
861               assembly[0].getPotentialEnergy().getTotalEnergy(),
862               dataFiles[i].getWeight() * sigmaAMinimize[i].calculateLikelihood(),
863               dataFiles[i].getWeight() * refinementData[i].llkF));
864     }
865 
866     return sb.toString();
867   }
868 
869   /** {@inheritDoc} */
870   @Override
871   public String printOptimizationHeader() {
872     return "R  Rfree";
873   }
874 
875   /** {@inheritDoc} */
876   @Override
877   public String printOptimizationUpdate() {
878     StringBuilder sb = new StringBuilder();
879     for (int i = 0; i < n; i++) {
880       sb.append(format("%6.2f %6.2f ", crystalStats[i].getR(), crystalStats[i].getRFree()));
881     }
882 
883     return sb.toString();
884   }
885 
886   /** Print scale and R statistics for all datasets associated with the model. */
887   public void printScaleAndR() {
888     for (int i = 0; i < n; i++) {
889       if (!scaled[i]) {
890         scaleBulkFit(i);
891       }
892       crystalStats[i].printScaleStats();
893       crystalStats[i].printRStats();
894     }
895   }
896 
897   /** Print all statistics for all datasets associated with the model. */
898   public void printStats() {
899     int nat = 0;
900     int nnonh = 0;
901 
902     // Note - we are including inactive and/or un-used atoms in the following loop.
903     for (Atom a : refinementModel.getTotalAtomList()) {
904       if (a.getOccupancy() == 0.0) {
905         continue;
906       }
907       nat++;
908       if (a.getAtomicNumber() == 1) {
909         continue;
910       }
911       nnonh++;
912     }
913 
914     for (int i = 0; i < n; i++) {
915       if (!scaled[i]) {
916         scaleBulkFit(i);
917       }
918 
919       String sb =
920           format(
921               " Statistics for Data Set %d of %d\n\n"
922                   + "  Weight:     %6.2f\n  Neutron data: %4s\n"
923                   + "  Model:        %s\n  Data file:    %s\n",
924               i + 1,
925               n,
926               dataFiles[i].getWeight(),
927               dataFiles[i].isNeutron(),
928               modelName,
929               dataFiles[i].getFilename());
930       logger.info(sb);
931 
932       crystalStats[i].printScaleStats();
933       crystalStats[i].printDPIStats(nnonh, nat);
934       crystalStats[i].printHKLStats();
935       crystalStats[i].printSNStats();
936       crystalStats[i].printRStats();
937     }
938   }
939 
940   /** Scale model and fit bulk solvent to all data. */
941   public void scaleBulkFit() {
942     for (int i = 0; i < n; i++) {
943       scaleBulkFit(i);
944     }
945   }
946 
947   /**
948    * Scale model and fit bulk solvent to dataset i of n.
949    *
950    * @param i a int.
951    */
952   public void scaleBulkFit(int i) {
953     String sb =
954         format(
955             " Scaling Data Set %d of %d\n\n  Weight: %6.2f\n  Neutron data: %s\n  Model: %s\n  Data file: %s",
956             i + 1,
957             n,
958             dataFiles[i].getWeight(),
959             dataFiles[i].isNeutron(),
960             modelName,
961             dataFiles[i].getFilename());
962     logger.info(sb);
963 
964     // reset some values
965     refinementData[i].bulkSolventK = 0.33;
966     refinementData[i].bulkSolventUeq = 50.0 / (8.0 * Math.PI * Math.PI);
967     refinementData[i].modelScaleK = 0.0;
968     for (int j = 0; j < 6; j++) {
969       refinementData[i].modelAnisoB[j] = 0.0;
970     }
971 
972     // run FFTs
973     crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc);
974     if (solventModel != SolventModel.NONE) {
975       crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs);
976     }
977 
978     // initialize minimizers
979     scaleBulkMinimize[i] =
980         new ScaleBulkMinimize(
981             reflectionList[i], refinementData[i], crystalReciprocalSpacesFs[i], parallelTeam);
982     splineMinimize[i] =
983         new SplineMinimize(
984             reflectionList[i], refinementData[i], refinementData[i].spline, SplineEnergy.Type.FOFC);
985 
986     // minimize
987     if (solventModel != SolventModel.NONE && gridSearch) {
988       scaleBulkMinimize[i].minimize(6, xrayScaleTol);
989       scaleBulkMinimize[i].gridOptimize();
990     }
991     scaleBulkMinimize[i].minimize(6, xrayScaleTol);
992 
993     // sigmaA / LLK calculation
994     sigmaAMinimize[i] = new SigmaAMinimize(reflectionList[i], refinementData[i], parallelTeam);
995     sigmaAMinimize[i].minimize(7, sigmaATol);
996 
997     if (splineFit) {
998       splineMinimize[i].minimize(7, 1e-5);
999     }
1000 
1001     scaled[i] = true;
1002   }
1003 
1004   /**
1005    * set the bulk solvent parameters for a given bulk solvent model
1006    *
1007    * @param a typically the width of the atom
1008    * @param b typically the rate with which the atom transitions to bulk
1009    * @see CrystalReciprocalSpace#setSolventAB(double, double)
1010    */
1011   public void setSolventAB(double a, double b) {
1012     for (int i = 0; i < n; i++) {
1013       if (solventModel != SolventModel.NONE) {
1014         refinementData[i].solventA = a;
1015         refinementData[i].solventB = b;
1016         crystalReciprocalSpacesFs[i].setSolventAB(a, b);
1017       }
1018     }
1019   }
1020 
1021   /** Perform 10 Fc calculations for the purposes of timings. */
1022   public void timings() {
1023     logger.info(" Performing 10 Fc calculations for timing...");
1024     for (int i = 0; i < n; i++) {
1025       for (int j = 0; j < 10; j++) {
1026         crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc, true);
1027         crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs, true);
1028         crystalReciprocalSpacesFc[i].computeAtomicGradients(
1029             refinementData[i].dFc,
1030             refinementData[i].freeR,
1031             refinementData[i].rFreeFlag,
1032             RefinementMode.COORDINATES,
1033             true);
1034         crystalReciprocalSpacesFs[i].computeAtomicGradients(
1035             refinementData[i].dFs,
1036             refinementData[i].freeR,
1037             refinementData[i].rFreeFlag,
1038             RefinementMode.COORDINATES,
1039             true);
1040       }
1041     }
1042   }
1043 
1044   /**
1045    * write current datasets to MTZ files
1046    *
1047    * @param filename output filename, or filename root for multiple datasets
1048    * @see MTZWriter#write()
1049    */
1050   public void writeData(String filename) {
1051     if (n == 1) {
1052       writeData(filename, 0);
1053     } else {
1054       for (int i = 0; i < n; i++) {
1055         writeData(removeExtension(filename) + "_" + i + ".mtz", i);
1056       }
1057     }
1058   }
1059 
1060   /**
1061    * write dataset i to MTZ file
1062    *
1063    * @param filename output filename
1064    * @param i dataset to write out
1065    * @see MTZWriter#write()
1066    */
1067   public void writeData(String filename, int i) {
1068     MTZWriter mtzwriter;
1069     if (scaled[i]) {
1070       mtzwriter = new MTZWriter(reflectionList[i], refinementData[i], filename);
1071     } else {
1072       mtzwriter = new MTZWriter(reflectionList[i], refinementData[i], filename, MTZType.DATAONLY);
1073     }
1074     mtzwriter.write();
1075   }
1076 
1077   /**
1078    * write 2Fo-Fc and Fo-Fc maps for all datasets
1079    *
1080    * @param filename output root filename for Fo-Fc and 2Fo-Fc maps
1081    */
1082   public void writeMaps(String filename) {
1083     if (n == 1) {
1084       writeMaps(filename, 0);
1085     } else {
1086       for (int i = 0; i < n; i++) {
1087         writeMaps(removeExtension(filename) + "_" + i + ".map", i);
1088       }
1089     }
1090   }
1091 
1092   /**
1093    * write 2Fo-Fc and Fo-Fc maps for a datasets
1094    *
1095    * @param filename output root filename for Fo-Fc and 2Fo-Fc maps
1096    * @param i a int.
1097    */
1098   private void writeMaps(String filename, int i) {
1099     if (!scaled[i]) {
1100       scaleBulkFit(i);
1101     }
1102 
1103     // Fo-Fc
1104     crystalReciprocalSpacesFc[i].computeAtomicGradients(
1105         refinementData[i].foFc1,
1106         refinementData[i].freeR,
1107         refinementData[i].rFreeFlag,
1108         RefinementMode.COORDINATES);
1109     double[] densityGrid = crystalReciprocalSpacesFc[i].getDensityGrid();
1110     int extx = (int) crystalReciprocalSpacesFc[i].getXDim();
1111     int exty = (int) crystalReciprocalSpacesFc[i].getYDim();
1112     int extz = (int) crystalReciprocalSpacesFc[i].getZDim();
1113 
1114     CCP4MapWriter mapwriter = new CCP4MapWriter(extx, exty, extz, crystal[i],
1115         removeExtension(filename) + "_fofc.map");
1116     mapwriter.write(densityGrid);
1117 
1118     // 2Fo-Fc
1119     crystalReciprocalSpacesFc[i].computeAtomicGradients(
1120         refinementData[i].foFc2,
1121         refinementData[i].freeR,
1122         refinementData[i].rFreeFlag,
1123         RefinementMode.COORDINATES);
1124     densityGrid = crystalReciprocalSpacesFc[i].getDensityGrid();
1125     extx = (int) crystalReciprocalSpacesFc[i].getXDim();
1126     exty = (int) crystalReciprocalSpacesFc[i].getYDim();
1127     extz = (int) crystalReciprocalSpacesFc[i].getZDim();
1128 
1129     mapwriter = new CCP4MapWriter(extx, exty, extz, crystal[i],
1130         removeExtension(filename) + "_2fofc.map");
1131     mapwriter.write(densityGrid);
1132   }
1133 
1134   /**
1135    * Write current model to PDB file.
1136    *
1137    * @param filename output PDB filename
1138    */
1139   public void writeModel(String filename) {
1140     StringBuilder remark = new StringBuilder();
1141 
1142     File file = version(new File(filename));
1143     CompositeConfiguration properties = assembly[0].getProperties();
1144     ForceField forceField = assembly[0].getForceField();
1145     PDBFilter pdbFilter = new PDBFilter(file, Arrays.asList(assembly), forceField, properties);
1146 
1147     Date now = new Date();
1148     SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss ");
1149     remark.append("REMARK FFX output ISO-8601 date: ").append(sdf.format(now)).append("\n");
1150     remark.append("REMARK\n");
1151     remark.append("REMARK   3\n");
1152     remark.append("REMARK   3 REFINEMENT\n");
1153     remark.append("REMARK   3   PROGRAM     : FORCE FIELD X\n");
1154     remark.append("REMARK   3\n");
1155     for (int i = 0; i < n; i++) {
1156       remark.append("REMARK   3  DATA SET ").append(i + 1).append("\n");
1157       if (dataFiles[i].isNeutron()) {
1158         remark.append("REMARK   3   DATA SET TYPE   : NEUTRON\n");
1159       } else {
1160         remark.append("REMARK   3   DATA SET TYPE   : X-RAY\n");
1161       }
1162       remark.append("REMARK   3   DATA SET WEIGHT : ").append(dataFiles[i].getWeight()).append("\n");
1163       remark.append("REMARK   3\n");
1164       remark.append(crystalStats[i].getPDBHeaderString());
1165     }
1166     for (int i = 0; i < assembly.length; i++) {
1167       remark.append("REMARK   3  CHEMICAL SYSTEM ").append(i + 1).append("\n");
1168       remark.append(assembly[i].getPotentialEnergy().getPDBHeaderString());
1169     }
1170     pdbFilter.writeFileWithHeader(file, remark);
1171   }
1172 
1173   /**
1174    * Write current model to PDB file.
1175    *
1176    * @param filename output PDB filename
1177    */
1178   public void writeModel(String filename, Set<Atom> excludeAtoms, double pH) {
1179     StringBuilder remark = new StringBuilder();
1180 
1181     File file = version(new File(filename));
1182     PDBFilter pdbFilter = new PDBFilter(file, Arrays.asList(assembly),null, null);
1183     if(pH > 0){
1184       pdbFilter.setRotamerTitration(true);
1185     }
1186 
1187     Date now = new Date();
1188     SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss ");
1189     remark.append("REMARK FFX output ISO-8601 date: ").append(sdf.format(now)).append("\n");
1190     remark.append("REMARK\n");
1191     remark.append("REMARK   3\n");
1192     remark.append("REMARK   3 REFINEMENT\n");
1193     remark.append("REMARK   3   PROGRAM     : FORCE FIELD X\n");
1194     remark.append("REMARK   3\n");
1195     for (int i = 0; i < n; i++) {
1196       remark.append("REMARK   3  DATA SET ").append(i + 1).append("\n");
1197       if (dataFiles[i].isNeutron()) {
1198         remark.append("REMARK   3   DATA SET TYPE   : NEUTRON\n");
1199       } else {
1200         remark.append("REMARK   3   DATA SET TYPE   : X-RAY\n");
1201       }
1202       remark.append("REMARK   3   DATA SET WEIGHT : ").append(dataFiles[i].getWeight()).append("\n");
1203       remark.append("REMARK   3\n");
1204       remark.append(crystalStats[i].getPDBHeaderString());
1205     }
1206     for (int i = 0; i < assembly.length; i++) {
1207       remark.append("REMARK   3  CHEMICAL SYSTEM ").append(i + 1).append("\n");
1208       remark.append(assembly[i].getPotentialEnergy().getPDBHeaderString());
1209     }
1210     remark.append("REMARK   3   TITRATION PH   : \n").append(pH).append("\n");
1211     String[] remarks = remark.toString().split("\n");
1212     pdbFilter.writeFile(file, false, excludeAtoms,true, true, remarks);
1213   }
1214 
1215   /**
1216    * Write bulk solvent mask for all datasets to a CCP4 map file.
1217    *
1218    * @param filename output filename, or output root filename for multiple datasets
1219    * @see CCP4MapWriter#write(double[], boolean)
1220    */
1221   public void writeSolventMask(String filename) {
1222     if (n == 1) {
1223       writeSolventMask(filename, 0);
1224     } else {
1225       for (int i = 0; i < n; i++) {
1226         writeSolventMask(removeExtension(filename) + "_" + i + ".map", i);
1227       }
1228     }
1229   }
1230 
1231   /**
1232    * Write bulk solvent mask for all datasets to a CNS map file.
1233    *
1234    * @param filename output filename, or output root filename for multiple datasets
1235    */
1236   public void writeSolventMaskCNS(String filename) {
1237     if (n == 1) {
1238       writeSolventMaskCNS(filename, 0);
1239     } else {
1240       for (int i = 0; i < n; i++) {
1241         writeSolventMaskCNS(removeExtension(filename) + "_" + i + ".map", i);
1242       }
1243     }
1244   }
1245 
1246   /**
1247    * Move the atomic coordinates for the FFT calculation - this is independent of the {@link
1248    * ffx.potential.bonded.Atom} object coordinates as it uses a linearized array
1249    *
1250    * @param x array of coordinates to move atoms to
1251    * @see CrystalReciprocalSpace#setCoordinates(double[])
1252    */
1253   void setFFTCoordinates(double[] x) {
1254     for (int i = 0; i < n; i++) {
1255       crystalReciprocalSpacesFc[i].setCoordinates(x);
1256       crystalReciprocalSpacesFs[i].setCoordinates(x);
1257     }
1258   }
1259 
1260   /**
1261    * Determine the total atomic gradients due to the diffraction data - performs an inverse FFT
1262    *
1263    * @param refinementMode the {@link RefinementMinimize.RefinementMode refinement mode}
1264    *     requested
1265    * @see CrystalReciprocalSpace#computeAtomicGradients(double[][], int[], int,
1266    *     ffx.xray.RefinementMinimize.RefinementMode, boolean) computeAtomicGradients
1267    */
1268   void computeAtomicGradients(RefinementMode refinementMode) {
1269     for (int i = 0; i < n; i++) {
1270       crystalReciprocalSpacesFc[i].computeAtomicGradients(
1271           refinementData[i].dFc,
1272           refinementData[i].freeR,
1273           refinementData[i].rFreeFlag,
1274           refinementMode);
1275       crystalReciprocalSpacesFs[i].computeAtomicGradients(
1276           refinementData[i].dFs,
1277           refinementData[i].freeR,
1278           refinementData[i].rFreeFlag,
1279           refinementMode);
1280     }
1281   }
1282 
1283   /**
1284    * compute total crystallographic likelihood target
1285    *
1286    * @return the total -log likelihood
1287    * @see SigmaAMinimize#calculateLikelihood()
1288    */
1289   double computeLikelihood() {
1290     double e = 0.0;
1291     for (int i = 0; i < n; i++) {
1292       e += dataFiles[i].getWeight() * sigmaAMinimize[i].calculateLikelihood();
1293     }
1294     return e;
1295   }
1296 
1297   /**
1298    * setLambdaTerm.
1299    *
1300    * @param lambdaTerm a boolean.
1301    */
1302   void setLambdaTerm(boolean lambdaTerm) {
1303     for (int i = 0; i < n; i++) {
1304       crystalReciprocalSpacesFc[i].setLambdaTerm(lambdaTerm);
1305       crystalReciprocalSpacesFs[i].setLambdaTerm(lambdaTerm);
1306     }
1307   }
1308 
1309   /**
1310    * Write bulk solvent mask for dataset i to a CNS map file.
1311    *
1312    * @param filename output filename
1313    * @param i dataset to write out
1314    */
1315   private void writeSolventMaskCNS(String filename, int i) {
1316     if (solventModel != SolventModel.NONE) {
1317       try {
1318         PrintWriter cnsfile = new PrintWriter(new BufferedWriter(new FileWriter(filename)));
1319         cnsfile.println(" ANOMalous=FALSE");
1320         cnsfile.println(" DECLare NAME=FS DOMAin=RECIprocal TYPE=COMP END");
1321         for (HKL ih : reflectionList[i].hklList) {
1322           int j = ih.getIndex();
1323           cnsfile.printf(
1324               " INDE %d %d %d FS= %.4f %.4f\n",
1325               ih.getH(),
1326               ih.getK(),
1327               ih.getL(),
1328               refinementData[i].fsF(j),
1329               Math.toDegrees(refinementData[i].fsPhi(j)));
1330         }
1331         cnsfile.close();
1332       } catch (Exception e) {
1333         String message = "Fatal exception evaluating structure factors.\n";
1334         logger.log(Level.SEVERE, message, e);
1335         System.exit(-1);
1336       }
1337     }
1338   }
1339 
1340   /**
1341    * Write bulk solvent mask for dataset i to a CCP4 map file.
1342    *
1343    * @param filename output filename
1344    * @param i dataset to write out
1345    * @see CCP4MapWriter#write(double[], boolean)
1346    */
1347   private void writeSolventMask(String filename, int i) {
1348     if (solventModel != SolventModel.NONE) {
1349       CCP4MapWriter mapwriter =
1350           new CCP4MapWriter(
1351               (int) crystalReciprocalSpacesFs[i].getXDim(),
1352               (int) crystalReciprocalSpacesFs[i].getYDim(),
1353               (int) crystalReciprocalSpacesFs[i].getZDim(),
1354               crystal[i],
1355               filename);
1356       mapwriter.write(crystalReciprocalSpacesFs[i].getSolventGrid());
1357     }
1358   }
1359 
1360   /**
1361    * Getter for the field <code>bSimWeight</code>.
1362    *
1363    * @return the bSimWeight
1364    */
1365   double getbSimWeight() {
1366     return bSimWeight;
1367   }
1368 
1369   /**
1370    * Getter for the field <code>bNonZeroWeight</code>.
1371    *
1372    * @return the bNonZeroWeight
1373    */
1374   double getbNonZeroWeight() {
1375     return bNonZeroWeight;
1376   }
1377 
1378   /**
1379    * Getter for the field <code>bMass</code>.
1380    *
1381    * @return the bMass
1382    */
1383   double getbMass() {
1384     return bMass;
1385   }
1386 
1387   /**
1388    * isResidueBFactor.
1389    *
1390    * @return the residueBFactor
1391    */
1392   boolean isResidueBFactor() {
1393     return residueBFactor;
1394   }
1395 
1396   /**
1397    * Getter for the field <code>nResidueBFactor</code>.
1398    *
1399    * @return the nResidueBFactor
1400    */
1401   int getnResidueBFactor() {
1402     return nResidueBFactor;
1403   }
1404 
1405   /**
1406    * isAddAnisou.
1407    *
1408    * @return the addAnisou
1409    */
1410   boolean isAddAnisou() {
1411     return addAnisou;
1412   }
1413 
1414   /**
1415    * Getter for the field <code>occMass</code>.
1416    *
1417    * @return the occMass
1418    */
1419   double getOccMass() {
1420     return occMass;
1421   }
1422 }