1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.xray;
39
40 import edu.rit.pj.ParallelTeam;
41 import ffx.crystal.Crystal;
42 import ffx.crystal.HKL;
43 import ffx.crystal.ReflectionList;
44 import ffx.crystal.Resolution;
45 import ffx.potential.MolecularAssembly;
46 import ffx.potential.Utilities;
47 import ffx.potential.bonded.Atom;
48 import ffx.potential.parameters.ForceField;
49 import ffx.potential.parsers.PDBFilter;
50 import ffx.xray.parallel.GridMethod;
51 import ffx.xray.scatter.NeutronFormFactor;
52 import ffx.xray.solvent.SolventModel;
53 import ffx.xray.parsers.CCP4MapWriter;
54 import ffx.xray.parsers.DiffractionFile;
55 import ffx.xray.parsers.MTZWriter;
56 import ffx.xray.parsers.MTZWriter.MTZType;
57 import ffx.xray.refine.RefinementMode;
58 import ffx.xray.refine.RefinementModel;
59 import ffx.xray.scatter.XRayFormFactor;
60 import org.apache.commons.configuration2.CompositeConfiguration;
61
62 import java.io.BufferedWriter;
63 import java.io.File;
64 import java.io.FileWriter;
65 import java.io.PrintWriter;
66 import java.text.SimpleDateFormat;
67 import java.util.Arrays;
68 import java.util.Date;
69 import java.util.Set;
70 import java.util.logging.Level;
71 import java.util.logging.Logger;
72
73 import static ffx.numerics.math.ScalarMath.b2u;
74 import static ffx.utilities.TinkerUtils.version;
75 import static java.lang.Math.abs;
76 import static java.lang.String.format;
77 import static java.util.Arrays.fill;
78 import static org.apache.commons.io.FilenameUtils.removeExtension;
79
80
81
82
83
84
85
86 public class DiffractionData implements DataContainer {
87
88 private static final Logger logger = Logger.getLogger(DiffractionData.class.getName());
89 private final MolecularAssembly[] assembly;
90 private final String modelName;
91 private final DiffractionFile[] dataFiles;
92 private final int n;
93 private final Crystal[] crystal;
94 private final Resolution[] resolution;
95 private final ReflectionList[] reflectionList;
96 private final DiffractionRefinementData[] refinementData;
97 private final CrystalReciprocalSpace[] crystalReciprocalSpacesFc;
98 private final CrystalReciprocalSpace[] crystalReciprocalSpacesFs;
99 private final SolventModel solventModel;
100 private final RefinementModel refinementModel;
101 private final boolean use_3g;
102 private final double aRadBuff;
103 private final double xrayScaleTol;
104 private final double sigmaATol;
105 private final double bSimWeight;
106 private final double bNonZeroWeight;
107 private final boolean nativeEnvironmentApproximation;
108 private final SigmaAMinimize[] sigmaAMinimize;
109 private final CrystalStats[] crystalStats;
110 private ParallelTeam parallelTeam;
111 private final GridMethod gridMethod;
112 private final boolean[] scaled;
113 private double xWeight;
114
115
116
117 private final boolean gridSearch;
118
119
120
121 private final boolean splineFit;
122
123
124
125
126
127
128
129
130
131
132
133
134 public DiffractionData(MolecularAssembly[] molecularAssemblies, CompositeConfiguration properties,
135 SolventModel solventModel, DiffractionFile... dataFiles) {
136 this.assembly = molecularAssemblies;
137 this.solventModel = solventModel;
138 this.modelName = molecularAssemblies[0].getName();
139 this.dataFiles = dataFiles;
140 this.n = dataFiles.length;
141
142 int rflag = properties.getInt("rfree-flag", -1);
143
144 double fsigfCutoff = properties.getDouble("f-sigf-cutoff", -1.0);
145 gridSearch = properties.getBoolean("solvent-grid-search", false);
146 splineFit = !properties.getBoolean("no-spline-fit", false);
147 use_3g = properties.getBoolean("use-3g", true);
148 aRadBuff = properties.getDouble("scattering-buffer", 0.75);
149 double sampling = properties.getDouble("sampling", 0.6);
150 xrayScaleTol = properties.getDouble("xray-scale-tol", 1e-4);
151 sigmaATol = properties.getDouble("sigmaa-tol", 0.05);
152 xWeight = properties.getDouble("data-weight", 1.0);
153 bSimWeight = properties.getDouble("b-sim-weight", 1.0);
154 bNonZeroWeight = properties.getDouble("b-nonzero-weight", 1.0);
155 boolean residueBFactor = properties.getBoolean("residue-bfactor", false);
156 int nResidueBFactor = properties.getInt("n-residue-bfactor", 1);
157 boolean addAnisou = properties.getBoolean("add-anisou", false);
158 boolean refineMolOcc = properties.getBoolean("refine-mol-occ", false);
159
160 ForceField forceField = molecularAssemblies[0].getForceField();
161 nativeEnvironmentApproximation =
162 forceField.getBoolean("NATIVE_ENVIRONMENT_APPROXIMATION", false);
163
164 crystal = new Crystal[n];
165 resolution = new Resolution[n];
166 reflectionList = new ReflectionList[n];
167 refinementData = new DiffractionRefinementData[n];
168 sigmaAMinimize = new SigmaAMinimize[n];
169 crystalStats = new CrystalStats[n];
170
171 for (int i = 0; i < n; i++) {
172 crystal[i] = Crystal.checkProperties(properties);
173 if (crystal[i] == null) {
174 crystal[i] = molecularAssemblies[i].getCrystal().getUnitCell();
175 }
176 File dataFile = new File(dataFiles[i].getFilename());
177 double res = dataFiles[i].getDiffractionfilter().getResolution(dataFile, crystal[i]);
178 if (dataFiles[i].isNeutron()) {
179 resolution[i] = Resolution.checkProperties(properties, true, res);
180 } else {
181 resolution[i] = Resolution.checkProperties(properties, false, res);
182 }
183 if (resolution[i] == null) {
184 logger.severe(" Resolution could not be determined from property or reflection files.");
185 }
186 reflectionList[i] = new ReflectionList(crystal[i], resolution[i], properties);
187
188 logger.info(resolution[i].toString());
189
190 refinementData[i] = new DiffractionRefinementData(properties, reflectionList[i]);
191 dataFiles[i].getDiffractionfilter().readFile(dataFile, reflectionList[i], refinementData[i], properties);
192 }
193
194 if (logger.isLoggable(Level.INFO)) {
195 StringBuilder sb = new StringBuilder();
196 sb.append("\n X-ray Refinement Settings\n\n");
197 sb.append(" Target Function\n");
198 sb.append(" X-ray refinement weight: ").append(xWeight).append("\n");
199 sb.append(" Use cctbx 3 Gaussians: ").append(use_3g).append("\n");
200 sb.append(" Atomic form factor radius buffer: ").append(aRadBuff).append("\n");
201 sb.append(" Reciprocal space sampling rate: ").append(sampling).append("\n");
202 sb.append(" Resolution dependent spline scale: ").append(splineFit).append("\n");
203 sb.append(" Solvent grid search: ").append(gridSearch).append("\n");
204 sb.append(" X-ray scale fit tolerance: ").append(xrayScaleTol).append("\n");
205 sb.append(" Sigma A fit tolerance: ").append(sigmaATol).append("\n");
206 sb.append(" Native environment approximation: ").append(nativeEnvironmentApproximation)
207 .append("\n");
208 sb.append(" Reflections\n");
209 sb.append(" F/sigF cutoff: ").append(fsigfCutoff).append("\n");
210 sb.append(" R Free flag (-1 auto-determine from the data): ").append(rflag).append("\n");
211 sb.append(" Number of bins: ").append(reflectionList[0].nBins).append("\n");
212 sb.append(" B-Factors\n");
213 sb.append(" Similarity weight: ").append(bSimWeight).append("\n");
214
215
216 sb.append(" Refine by residue: ").append(residueBFactor).append("\n");
217 if (residueBFactor) {
218 sb.append(" Number of residues per B: ").append(nResidueBFactor).append(")\n");
219 }
220 sb.append(" Add ANISOU for refinement: ").append(addAnisou).append("\n");
221 sb.append(" Occupancies\n");
222 sb.append(" Refine on molecules: ").append(refineMolOcc).append("\n");
223
224
225 logger.info(sb.toString());
226 }
227
228
229 refinementModel = new RefinementModel(molecularAssemblies);
230
231
232
233 for (Atom a : refinementModel.getScatteringAtoms()) {
234
235
236
237 XRayFormFactor atomff = new XRayFormFactor(a, use_3g, 2.0);
238
239 if (a.getOccupancy() == 0.0) {
240 a.setFormFactorWidth(1.0);
241 continue;
242 }
243
244 double arad;
245 try {
246 arad = a.getVDWType().radius * 0.5;
247 } catch (NullPointerException ex) {
248 logger.warning(format(" Failure to get van der Waals type for atom %s; ensure the vdW term is enabled!", a));
249 throw ex;
250 }
251 double[] xyz = new double[3];
252 xyz[0] = a.getX() + arad;
253 xyz[1] = a.getY();
254 xyz[2] = a.getZ();
255 double rho = atomff.rho(0.0, 1.0, xyz);
256
257
258 while (rho > 0.001) {
259 arad += 0.1;
260 xyz[0] = a.getX() + arad;
261 rho = atomff.rho(0.0, 1.0, xyz);
262
263 }
264 arad += aRadBuff;
265 a.setFormFactorWidth(arad);
266 }
267
268
269 crystalReciprocalSpacesFc = new CrystalReciprocalSpace[n];
270 crystalReciprocalSpacesFs = new CrystalReciprocalSpace[n];
271
272 parallelTeam = molecularAssemblies[0].getParallelTeam();
273
274 String gridString = properties.getString("grid-method", "SLICE");
275 gridMethod = GridMethod.parse(gridString);
276
277 parallelTeam = new ParallelTeam();
278 for (int i = 0; i < n; i++) {
279
280 crystalReciprocalSpacesFc[i] = new CrystalReciprocalSpace(
281 reflectionList[i], refinementModel.getScatteringAtoms(),
282 parallelTeam, parallelTeam, false,
283 this.dataFiles[i].isNeutron(), SolventModel.NONE, gridMethod);
284 refinementData[i].setCrystalReciprocalSpaceFc(crystalReciprocalSpacesFc[i]);
285 crystalReciprocalSpacesFc[i].setUse3G(use_3g);
286 crystalReciprocalSpacesFc[i].setWeight(this.dataFiles[i].getWeight());
287 crystalReciprocalSpacesFc[i].lambdaTerm = false;
288 crystalReciprocalSpacesFc[i].setNativeEnvironmentApproximation(
289 nativeEnvironmentApproximation);
290
291
292 crystalReciprocalSpacesFs[i] = new CrystalReciprocalSpace(
293 reflectionList[i], refinementModel.getScatteringAtoms(),
294 parallelTeam, parallelTeam, true,
295 this.dataFiles[i].isNeutron(), solventModel, gridMethod);
296 refinementData[i].setCrystalReciprocalSpaceFs(crystalReciprocalSpacesFs[i]);
297 crystalReciprocalSpacesFs[i].setUse3G(use_3g);
298 crystalReciprocalSpacesFs[i].setWeight(this.dataFiles[i].getWeight());
299 crystalReciprocalSpacesFs[i].lambdaTerm = false;
300 crystalReciprocalSpacesFs[i].setNativeEnvironmentApproximation(
301 nativeEnvironmentApproximation);
302 crystalStats[i] = new CrystalStats(reflectionList[i], refinementData[i]);
303 }
304
305 scaled = new boolean[n];
306 fill(scaled, false);
307 }
308
309
310
311
312
313
314
315
316
317 public void AverageFc(MolecularAssembly[] assembly, int index) {
318 RefinementModel tmprefinementmodel = new RefinementModel(assembly);
319
320
321 for (Atom a : tmprefinementmodel.getScatteringAtoms()) {
322 XRayFormFactor atomff = new XRayFormFactor(a, use_3g, 2.0);
323 if (a.getOccupancy() == 0.0) {
324 a.setFormFactorWidth(1.0);
325 continue;
326 }
327
328 double arad = a.getVDWType().radius * 0.5;
329 double[] xyz = new double[3];
330 xyz[0] = a.getX() + arad;
331 xyz[1] = a.getY();
332 xyz[2] = a.getZ();
333 double rho = atomff.rho(0.0, 1.0, xyz);
334 while (rho > 0.001) {
335 arad += 0.1;
336 xyz[0] = a.getX() + arad;
337 rho = atomff.rho(0.0, 1.0, xyz);
338 }
339 arad += aRadBuff;
340 a.setFormFactorWidth(arad);
341 }
342
343
344 for (int i = 0; i < n; i++) {
345 crystalReciprocalSpacesFc[i] =
346 new CrystalReciprocalSpace(
347 reflectionList[i],
348 tmprefinementmodel.getScatteringAtoms(),
349 parallelTeam,
350 parallelTeam,
351 false,
352 dataFiles[i].isNeutron(),
353 SolventModel.NONE,
354 gridMethod);
355 crystalReciprocalSpacesFc[i].setNativeEnvironmentApproximation(
356 nativeEnvironmentApproximation);
357 refinementData[i].setCrystalReciprocalSpaceFc(crystalReciprocalSpacesFc[i]);
358
359 crystalReciprocalSpacesFs[i] =
360 new CrystalReciprocalSpace(
361 reflectionList[i],
362 tmprefinementmodel.getScatteringAtoms(),
363 parallelTeam,
364 parallelTeam,
365 true,
366 dataFiles[i].isNeutron(),
367 solventModel,
368 gridMethod);
369 crystalReciprocalSpacesFs[i].setNativeEnvironmentApproximation(
370 nativeEnvironmentApproximation);
371 refinementData[i].setCrystalReciprocalSpaceFs(crystalReciprocalSpacesFs[i]);
372 }
373
374 int nhkl = refinementData[0].n;
375 double[][] fc = new double[nhkl][2];
376 double[][] fs = new double[nhkl][2];
377
378
379 for (int i = 0; i < n; i++) {
380 crystalReciprocalSpacesFc[i].computeDensity(fc);
381 if (solventModel != SolventModel.NONE) {
382 crystalReciprocalSpacesFs[i].computeDensity(fs);
383 }
384
385
386 for (int j = 0; j < refinementData[i].n; j++) {
387 refinementData[i].fc[j][0] += (fc[j][0] - refinementData[i].fc[j][0]) / index;
388 refinementData[i].fc[j][1] += (fc[j][1] - refinementData[i].fc[j][1]) / index;
389
390 refinementData[i].fs[j][0] += (fs[j][0] - refinementData[i].fs[j][0]) / index;
391 refinementData[i].fs[j][1] += (fs[j][1] - refinementData[i].fs[j][1]) / index;
392 }
393 }
394 }
395
396
397
398
399
400
401
402 public void computeAtomicDensity() {
403 for (int i = 0; i < n; i++) {
404 crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc);
405 if (solventModel != SolventModel.NONE) {
406 crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs);
407 }
408 }
409 }
410
411
412
413
414
415
416 public boolean destroy() {
417 try {
418 boolean underlyingShutdown = true;
419 for (MolecularAssembly assem : assembly) {
420
421 boolean thisShutdown = assem.destroy();
422 underlyingShutdown = underlyingShutdown && thisShutdown;
423 }
424 parallelTeam.shutdown();
425 return underlyingShutdown;
426 } catch (Exception ex) {
427 logger.warning(format(" Exception in shutting down a RealSpaceData: %s", ex));
428 logger.info(Utilities.stackTraceToString(ex));
429 return false;
430 }
431 }
432
433
434
435
436
437
438 public MolecularAssembly[] getAssembly() {
439 return assembly;
440 }
441
442
443
444
445
446
447 public Crystal[] getCrystal() {
448 return crystal;
449 }
450
451
452
453
454
455
456 public CrystalReciprocalSpace[] getCrystalReciprocalSpacesFc() {
457 return crystalReciprocalSpacesFc;
458 }
459
460
461
462
463
464
465 public CrystalReciprocalSpace[] getCrystalReciprocalSpacesFs() {
466 return crystalReciprocalSpacesFs;
467 }
468
469
470
471
472
473
474 public int getN() {
475 return n;
476 }
477
478
479
480
481
482
483 public ParallelTeam getParallelTeam() {
484 return parallelTeam;
485 }
486
487
488
489
490
491
492 public DiffractionRefinementData[] getRefinementData() {
493 return refinementData;
494 }
495
496
497
498
499 @Override
500 public RefinementModel getRefinementModel() {
501 return refinementModel;
502 }
503
504
505
506
507
508
509 public ReflectionList[] getReflectionList() {
510 return reflectionList;
511 }
512
513
514
515
516
517
518 public Resolution[] getResolution() {
519 return resolution;
520 }
521
522
523
524
525
526
527 public boolean[] getScaled() {
528 return scaled;
529 }
530
531
532
533
534 @Override
535 public double getWeight() {
536 return xWeight;
537 }
538
539
540
541
542 @Override
543 public void setWeight(double weight) {
544 this.xWeight = weight;
545 }
546
547
548
549
550 @Override
551 public String printEnergyUpdate() {
552 StringBuilder sb = new StringBuilder();
553 for (int i = 0; i < n; i++) {
554 sb.append(format(
555 " dataset %d (weight: %5.1f): R: %6.2f Rfree: %6.2f chemical energy: %8.2f likelihood: %8.2f free likelihood: %8.2f\n",
556 i + 1,
557 dataFiles[i].getWeight(),
558 crystalStats[i].getR(),
559 crystalStats[i].getRFree(),
560 assembly[0].getPotentialEnergy().getTotalEnergy(),
561 dataFiles[i].getWeight() * sigmaAMinimize[i].calculateLikelihood(),
562 dataFiles[i].getWeight() * refinementData[i].llkF));
563 }
564
565 return sb.toString();
566 }
567
568
569
570
571 @Override
572 public String printOptimizationHeader() {
573 return "R Rfree";
574 }
575
576
577
578
579 @Override
580 public String printOptimizationUpdate() {
581 StringBuilder sb = new StringBuilder();
582 for (int i = 0; i < n; i++) {
583 sb.append(format("%6.2f %6.2f ", crystalStats[i].getR(), crystalStats[i].getRFree()));
584 }
585
586 return sb.toString();
587 }
588
589
590
591
592 public void printStats() {
593 int numberOfScatteringAtoms = 0;
594 int nHeavyScatteringAtoms = 0;
595
596
597 for (Atom a : refinementModel.getScatteringAtoms()) {
598 if (a.getOccupancy() == 0.0) {
599 continue;
600 }
601 numberOfScatteringAtoms++;
602 if (a.getAtomicNumber() == 1) {
603 continue;
604 }
605 nHeavyScatteringAtoms++;
606 }
607
608 for (int i = 0; i < n; i++) {
609 if (!scaled[i]) {
610 scaleBulkFit(i);
611 }
612
613 String sb = format(" Statistics for Data Set %d of %d\n\n"
614 + " Weight: %6.2f\n Neutron data: %4s\n"
615 + " Model: %s\n Data file: %s\n",
616 i + 1, n, dataFiles[i].getWeight(), dataFiles[i].isNeutron(),
617 modelName, dataFiles[i].getFilename());
618 logger.info(sb);
619
620 crystalStats[i].printScaleStats();
621 crystalStats[i].printDPIStats(nHeavyScatteringAtoms, numberOfScatteringAtoms);
622 crystalStats[i].printHKLStats();
623 crystalStats[i].printSNStats();
624 crystalStats[i].printRStats();
625 }
626 }
627
628
629
630
631 public void scaleBulkFit() {
632 for (int i = 0; i < n; i++) {
633 scaleBulkFit(i);
634 }
635 }
636
637
638
639
640
641
642 public void scaleBulkFit(int i) {
643 String sb = format(
644 " Scaling Data Set %d of %d\n\n Weight: %6.2f\n Neutron data: %s\n Model: %s\n Data file: %s",
645 i + 1, n, dataFiles[i].getWeight(), dataFiles[i].isNeutron(),
646 modelName, dataFiles[i].getFilename());
647 logger.info(sb);
648
649
650
651 refinementData[i].modelScaleK = 0.0;
652
653 Arrays.fill(refinementData[i].modelAnisoB, 0.0);
654
655
656 if (dataFiles[i].isNeutron()) {
657
658
659
660 refinementData[i].bulkSolventK = 0.639;
661 } else {
662
663
664 refinementData[i].bulkSolventK = 0.334;
665 }
666 refinementData[i].bulkSolventUeq = b2u(50.0);
667
668
669 refinementData[i].crystalReciprocalSpaceFs.setDefaultSolventAB();
670
671
672 CompositeConfiguration properties = assembly[0].getProperties();
673 if (dataFiles[i].isNeutron()) {
674 if (properties.containsKey("neutron-bulk-solvent")) {
675 String string = properties.getString("neutron-bulk-solvent", "0.639 50.0");
676 String[] split = string.split("\\s+");
677 refinementData[i].bulkSolventK = Double.parseDouble(split[0]);
678 refinementData[i].bulkSolventUeq = b2u(Double.parseDouble(split[1]));
679 refinementData[i].bulkSolventFixed = true;
680 } else {
681 refinementData[i].bulkSolventFixed = false;
682 }
683 } else {
684 if (properties.containsKey("bulk-solvent")) {
685 String string = properties.getString("bulk-solvent", "0.334 50.0");
686 String[] split = string.split("\\s+");
687 refinementData[i].bulkSolventK = Double.parseDouble(split[0]);
688 refinementData[i].bulkSolventUeq = b2u(Double.parseDouble(split[1]));
689 refinementData[i].bulkSolventFixed = true;
690 } else {
691 refinementData[i].bulkSolventFixed = false;
692 }
693 }
694
695
696 crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc);
697 if (solventModel != SolventModel.NONE) {
698 crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs);
699 }
700
701
702 if (solventModel != SolventModel.NONE) {
703 boolean bulkSolventFixed = refinementData[i].bulkSolventFixed;
704
705
706 refinementData[i].bulkSolventFixed = true;
707 ScaleBulkMinimize scaleBulkMinimize = new ScaleBulkMinimize(reflectionList[i],
708 refinementData[i], crystalReciprocalSpacesFs[i], parallelTeam);
709 scaleBulkMinimize.minimize(xrayScaleTol);
710
711
712 if (!bulkSolventFixed) {
713 refinementData[i].bulkSolventFixed = false;
714 scaleBulkMinimize = new ScaleBulkMinimize(reflectionList[i],
715 refinementData[i], crystalReciprocalSpacesFs[i], parallelTeam);
716 }
717
718
719
720 if (gridSearch) {
721
722 scaleBulkMinimize.gridOptimizeBulkSolventModel();
723 }
724
725
726 if (!bulkSolventFixed) {
727 scaleBulkMinimize.gridOptimizeKsBs();
728 }
729
730
731 scaleBulkMinimize.minimize(xrayScaleTol);
732 } else {
733
734 ScaleBulkMinimize scaleBulkMinimize = new ScaleBulkMinimize(
735 reflectionList[i], refinementData[i], crystalReciprocalSpacesFs[i], parallelTeam);
736 scaleBulkMinimize.minimize(xrayScaleTol);
737 }
738
739
740 sigmaAMinimize[i] = new SigmaAMinimize(reflectionList[i], refinementData[i], parallelTeam);
741 sigmaAMinimize[i].minimize(sigmaATol);
742
743 if (splineFit) {
744
745 SplineMinimize splineMinimize = new SplineMinimize(
746 reflectionList[i], refinementData[i], refinementData[i].spline, SplineEnergy.SplineType.FOFC);
747 splineMinimize.minimize(1e-5);
748 }
749
750 scaled[i] = true;
751 }
752
753
754
755
756 public void timings() {
757 logger.info(" Performing 10 Fc calculations for timing...");
758 for (int i = 0; i < n; i++) {
759 for (int j = 0; j < 10; j++) {
760 crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc, true);
761 crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs, true);
762 crystalReciprocalSpacesFc[i].computeAtomicGradients(
763 refinementData[i].dFc, refinementData[i].freeR, refinementData[i].rFreeFlag,
764 RefinementMode.COORDINATES, true);
765 crystalReciprocalSpacesFs[i].computeAtomicGradients(
766 refinementData[i].dFs, refinementData[i].freeR, refinementData[i].rFreeFlag,
767 RefinementMode.COORDINATES, true);
768 }
769 }
770 }
771
772
773
774
775
776
777
778 public void writeData(String filename) {
779 if (n == 1) {
780 writeData(filename, 0);
781 } else {
782 for (int i = 0; i < n; i++) {
783 writeData(removeExtension(filename) + "_" + i + ".mtz", i);
784 }
785 }
786 }
787
788
789
790
791
792
793
794
795 public void writeData(String filename, int i) {
796 MTZWriter mtzwriter;
797 if (scaled[i]) {
798 mtzwriter = new MTZWriter(reflectionList[i], refinementData[i], filename);
799 } else {
800 mtzwriter = new MTZWriter(reflectionList[i], refinementData[i], filename, MTZType.DATAONLY);
801 }
802 mtzwriter.write();
803 }
804
805
806
807
808
809
810 public void writeMaps(String filename, boolean normalize) {
811 if (n == 1) {
812 writeMaps(filename, 0, normalize);
813 } else {
814 for (int i = 0; i < n; i++) {
815 writeMaps(removeExtension(filename) + "_" + i + ".map", i, normalize);
816 }
817 }
818 }
819
820
821
822
823
824
825
826
827 private void writeMaps(String filename, int i, boolean normalize) {
828 if (!scaled[i]) {
829 scaleBulkFit(i);
830 }
831
832
833 crystalReciprocalSpacesFc[i].computeAtomicGradients(refinementData[i].foFc1,
834 refinementData[i].freeR, refinementData[i].rFreeFlag, RefinementMode.COORDINATES);
835 double[] densityGrid = crystalReciprocalSpacesFc[i].getDensityGrid();
836 int extx = (int) crystalReciprocalSpacesFc[i].getXDim();
837 int exty = (int) crystalReciprocalSpacesFc[i].getYDim();
838 int extz = (int) crystalReciprocalSpacesFc[i].getZDim();
839
840 CCP4MapWriter mapwriter = new CCP4MapWriter(extx, exty, extz, crystal[i],
841 removeExtension(filename) + "_fofc.map");
842 mapwriter.write(densityGrid, normalize);
843
844
845 crystalReciprocalSpacesFc[i].computeAtomicGradients(refinementData[i].foFc2,
846 refinementData[i].freeR, refinementData[i].rFreeFlag, RefinementMode.COORDINATES);
847 densityGrid = crystalReciprocalSpacesFc[i].getDensityGrid();
848 extx = (int) crystalReciprocalSpacesFc[i].getXDim();
849 exty = (int) crystalReciprocalSpacesFc[i].getYDim();
850 extz = (int) crystalReciprocalSpacesFc[i].getZDim();
851
852 mapwriter = new CCP4MapWriter(extx, exty, extz, crystal[i],
853 removeExtension(filename) + "_2fofc.map");
854 mapwriter.write(densityGrid, normalize);
855 }
856
857
858
859
860
861
862 public void writeModel(String filename) {
863 StringBuilder remark = new StringBuilder();
864
865 File file = version(new File(filename));
866 CompositeConfiguration properties = assembly[0].getProperties();
867 ForceField forceField = assembly[0].getForceField();
868 PDBFilter pdbFilter = new PDBFilter(file, Arrays.asList(assembly), forceField, properties);
869
870 Date now = new Date();
871 SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss ");
872 remark.append("REMARK FFX output ISO-8601 date: ").append(sdf.format(now)).append("\n");
873 remark.append("REMARK\n");
874 remark.append("REMARK 3\n");
875 remark.append("REMARK 3 REFINEMENT\n");
876 remark.append("REMARK 3 PROGRAM : FORCE FIELD X\n");
877 remark.append("REMARK 3\n");
878 for (int i = 0; i < n; i++) {
879 remark.append("REMARK 3 DATA SET ").append(i + 1).append("\n");
880 if (dataFiles[i].isNeutron()) {
881 remark.append("REMARK 3 DATA SET TYPE : NEUTRON\n");
882 } else {
883 remark.append("REMARK 3 DATA SET TYPE : X-RAY\n");
884 }
885 remark.append("REMARK 3 DATA SET WEIGHT : ").append(dataFiles[i].getWeight()).append("\n");
886 remark.append("REMARK 3\n");
887 remark.append(crystalStats[i].getPDBHeaderString());
888 }
889 for (int i = 0; i < assembly.length; i++) {
890 remark.append("REMARK 3 CHEMICAL SYSTEM ").append(i + 1).append("\n");
891 remark.append(assembly[i].getPotentialEnergy().getPDBHeaderString());
892 }
893 pdbFilter.writeFileWithHeader(file, remark);
894 }
895
896
897
898
899
900
901 public void writeModel(String filename, Set<Atom> excludeAtoms, double pH) {
902 StringBuilder remark = new StringBuilder();
903
904 File file = version(new File(filename));
905 PDBFilter pdbFilter = new PDBFilter(file, Arrays.asList(assembly), null, null);
906 if (pH > 0) {
907 pdbFilter.setRotamerTitration(true);
908 }
909
910 Date now = new Date();
911 SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss ");
912 remark.append("REMARK FFX output ISO-8601 date: ").append(sdf.format(now)).append("\n");
913 remark.append("REMARK\n");
914 remark.append("REMARK 3\n");
915 remark.append("REMARK 3 REFINEMENT\n");
916 remark.append("REMARK 3 PROGRAM : FORCE FIELD X\n");
917 remark.append("REMARK 3\n");
918 for (int i = 0; i < n; i++) {
919 remark.append("REMARK 3 DATA SET ").append(i + 1).append("\n");
920 if (dataFiles[i].isNeutron()) {
921 remark.append("REMARK 3 DATA SET TYPE : NEUTRON\n");
922 } else {
923 remark.append("REMARK 3 DATA SET TYPE : X-RAY\n");
924 }
925 remark.append("REMARK 3 DATA SET WEIGHT : ").append(dataFiles[i].getWeight()).append("\n");
926 remark.append("REMARK 3\n");
927 remark.append(crystalStats[i].getPDBHeaderString());
928 }
929 for (int i = 0; i < assembly.length; i++) {
930 remark.append("REMARK 3 CHEMICAL SYSTEM ").append(i + 1).append("\n");
931 remark.append(assembly[i].getPotentialEnergy().getPDBHeaderString());
932 }
933 remark.append("REMARK 3 TITRATION PH : \n").append(pH).append("\n");
934 String[] remarks = remark.toString().split("\n");
935 pdbFilter.writeFile(file, false, excludeAtoms, true, true, remarks);
936 }
937
938
939
940
941
942
943
944 public void writeSolventMask(String filename) {
945 if (n == 1) {
946 writeSolventMask(filename, 0);
947 } else {
948 for (int i = 0; i < n; i++) {
949 writeSolventMask(removeExtension(filename) + "_" + i + ".map", i);
950 }
951 }
952 }
953
954
955
956
957
958
959 public void writeSolventMaskCNS(String filename) {
960 if (n == 1) {
961 writeSolventMaskCNS(filename, 0);
962 } else {
963 for (int i = 0; i < n; i++) {
964 writeSolventMaskCNS(removeExtension(filename) + "_" + i + ".map", i);
965 }
966 }
967 }
968
969
970
971
972
973
974
975 void updateCoordinates() {
976 for (int i = 0; i < n; i++) {
977 crystalReciprocalSpacesFc[i].updateCoordinates();
978 crystalReciprocalSpacesFs[i].updateCoordinates();
979 }
980 }
981
982
983
984
985
986
987
988
989
990 void computeAtomicGradients(RefinementMode refinementMode) {
991 for (int i = 0; i < n; i++) {
992 crystalReciprocalSpacesFc[i].computeAtomicGradients(refinementData[i].dFc,
993 refinementData[i].freeR, refinementData[i].rFreeFlag, refinementMode);
994 crystalReciprocalSpacesFs[i].computeAtomicGradients(refinementData[i].dFs,
995 refinementData[i].freeR, refinementData[i].rFreeFlag, refinementMode);
996 }
997 }
998
999
1000
1001
1002
1003
1004
1005 double computeLikelihood() {
1006 double e = 0.0;
1007 for (int i = 0; i < n; i++) {
1008 e += dataFiles[i].getWeight() * sigmaAMinimize[i].calculateLikelihood();
1009 }
1010 return e;
1011 }
1012
1013
1014
1015
1016
1017
1018 void setLambdaTerm(boolean lambdaTerm) {
1019 for (int i = 0; i < n; i++) {
1020 crystalReciprocalSpacesFc[i].setLambdaTerm(lambdaTerm);
1021 crystalReciprocalSpacesFs[i].setLambdaTerm(lambdaTerm);
1022 }
1023 }
1024
1025
1026
1027
1028
1029
1030
1031 private void writeSolventMaskCNS(String filename, int i) {
1032 if (solventModel != SolventModel.NONE) {
1033 try {
1034 PrintWriter cnsfile = new PrintWriter(new BufferedWriter(new FileWriter(filename)));
1035 cnsfile.println(" ANOMalous=FALSE");
1036 cnsfile.println(" DECLare NAME=FS DOMAin=RECIprocal TYPE=COMP END");
1037 for (HKL ih : reflectionList[i].hklList) {
1038 int j = ih.getIndex();
1039 cnsfile.printf(" INDE %d %d %d FS= %.4f %.4f\n",
1040 ih.getH(), ih.getK(), ih.getL(),
1041 refinementData[i].fsF(j), Math.toDegrees(refinementData[i].fsPhi(j)));
1042 }
1043 cnsfile.close();
1044 } catch (Exception e) {
1045 String message = "Fatal exception evaluating structure factors.\n";
1046 logger.log(Level.SEVERE, message, e);
1047 System.exit(-1);
1048 }
1049 }
1050 }
1051
1052
1053
1054
1055
1056
1057
1058
1059 private void writeSolventMask(String filename, int i) {
1060 if (solventModel != SolventModel.NONE) {
1061 CCP4MapWriter mapwriter =
1062 new CCP4MapWriter(
1063 (int) crystalReciprocalSpacesFs[i].getXDim(),
1064 (int) crystalReciprocalSpacesFs[i].getYDim(),
1065 (int) crystalReciprocalSpacesFs[i].getZDim(),
1066 crystal[i], filename);
1067 mapwriter.write(crystalReciprocalSpacesFs[i].getSolventGrid());
1068 }
1069 }
1070
1071
1072
1073
1074
1075
1076 double getbSimWeight() {
1077 return bSimWeight;
1078 }
1079
1080
1081
1082
1083
1084
1085 double getbNonZeroWeight() {
1086 return bNonZeroWeight;
1087 }
1088
1089 }