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