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 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
79
80
81
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
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
124 private final boolean gridSearch;
125
126 private final boolean splineFit;
127
128
129
130
131
132
133
134
135
136 public DiffractionData(MolecularAssembly assembly, CompositeConfiguration properties) {
137 this(new MolecularAssembly[] {assembly}, properties, POLYNOMIAL, new DiffractionFile(assembly));
138 }
139
140
141
142
143
144
145
146
147
148 public DiffractionData(
149 MolecularAssembly assembly, CompositeConfiguration properties, DiffractionFile... datafile) {
150 this(new MolecularAssembly[] {assembly}, properties, POLYNOMIAL, datafile);
151 }
152
153
154
155
156
157
158
159
160
161
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
174
175
176
177
178
179
180
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
192
193
194
195
196
197
198
199 public DiffractionData(MolecularAssembly[] assembly, CompositeConfiguration properties) {
200 this(assembly, properties, POLYNOMIAL, new DiffractionFile(assembly[0]));
201 }
202
203
204
205
206
207
208
209
210
211
212 public DiffractionData(
213 MolecularAssembly[] assembly,
214 CompositeConfiguration properties,
215 DiffractionFile... datafile) {
216 this(assembly, properties, POLYNOMIAL, datafile);
217 }
218
219
220
221
222
223
224
225
226
227
228
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
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
337
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
346
347 logger.info(sb.toString());
348 }
349
350
351 refinementModel = new RefinementModel(assembly, refineMolOcc);
352
353
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
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
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
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
444
445
446
447
448
449
450 public void AverageFc(MolecularAssembly[] assembly, int index) {
451 RefinementModel tmprefinementmodel = new RefinementModel(assembly, refineMolOcc);
452
453
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
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
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
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
534
535
536
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
549
550
551
552 public boolean destroy() {
553 try {
554 boolean underlyingShutdown = true;
555 for (MolecularAssembly assem : assembly) {
556
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
571
572
573
574 public Atom[] getActiveAtomArray() {
575 return getAtomArray();
576 }
577
578
579 @Override
580 public List<List<Molecule>> getAltMolecules() {
581 return refinementModel.getAltMolecules();
582 }
583
584
585 @Override
586 public List<List<Residue>> getAltResidues() {
587 return refinementModel.getAltResidues();
588 }
589
590
591
592
593
594
595 public MolecularAssembly[] getAssembly() {
596 return assembly;
597 }
598
599
600
601
602
603
604 @Override
605 public Atom[] getAtomArray() {
606 return refinementModel.getTotalAtomArray();
607 }
608
609
610
611
612
613
614 public Crystal[] getCrystal() {
615 return crystal;
616 }
617
618
619
620
621
622
623 public CrystalReciprocalSpace[] getCrystalReciprocalSpacesFc() {
624 return crystalReciprocalSpacesFc;
625 }
626
627
628
629
630
631
632 public CrystalReciprocalSpace[] getCrystalReciprocalSpacesFs() {
633 return crystalReciprocalSpacesFs;
634 }
635
636
637
638
639
640
641 public CrystalStats[] getCrystalStats() {
642 return crystalStats;
643 }
644
645
646
647
648
649
650 public DiffractionFile[] getDataFiles() {
651 return dataFiles;
652 }
653
654
655
656
657
658
659 public double getFsigfCutoff() {
660 return fsigfCutoff;
661 }
662
663
664
665
666
667
668 public String getModelName() {
669 return modelName;
670 }
671
672
673 @Override
674 public MolecularAssembly[] getMolecularAssemblies() {
675 return assembly;
676 }
677
678
679
680
681
682
683 public int getN() {
684 return n;
685 }
686
687
688
689
690
691
692 public ParallelTeam getParallelTeam() {
693 return parallelTeam;
694 }
695
696
697
698
699
700
701 public double getRCrystalStat() {
702 return crystalStats[0].getR();
703 }
704
705
706
707
708
709
710 public DiffractionRefinementData[] getRefinementData() {
711 return refinementData;
712 }
713
714
715 @Override
716 public RefinementModel getRefinementModel() {
717 return refinementModel;
718 }
719
720
721
722
723
724
725 public ReflectionList[] getReflectionList() {
726 return reflectionList;
727 }
728
729
730
731
732
733
734 public Resolution[] getResolution() {
735 return resolution;
736 }
737
738
739
740
741
742
743 public ScaleBulkMinimize[] getScaleBulkMinimize() {
744 return scaleBulkMinimize;
745 }
746
747
748
749
750
751
752 public boolean[] getScaled() {
753 return scaled;
754 }
755
756
757
758
759
760
761 public SigmaAMinimize[] getSigmaAMinimize() {
762 return sigmaAMinimize;
763 }
764
765
766
767
768
769
770 public double getSigmaATol() {
771 return sigmaATol;
772 }
773
774
775
776
777
778
779 public SolventModel getSolventModel() {
780 return solventModel;
781 }
782
783
784
785
786
787
788 public SplineMinimize[] getSplineMinimize() {
789 return splineMinimize;
790 }
791
792
793 @Override
794 public double getWeight() {
795 return xWeight;
796 }
797
798
799 @Override
800 public void setWeight(double weight) {
801 this.xWeight = weight;
802 }
803
804
805
806
807
808
809 public double getXrayScaleTol() {
810 return xrayScaleTol;
811 }
812
813
814
815
816
817
818 public double getaRadBuff() {
819 return aRadBuff;
820 }
821
822
823
824
825
826
827 public double getxWeight() {
828 return xWeight;
829 }
830
831
832
833
834
835
836 public boolean isRefineMolOcc() {
837 return refineMolOcc;
838 }
839
840
841
842
843
844
845 public boolean isUse_3g() {
846 return use_3g;
847 }
848
849
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
870 @Override
871 public String printOptimizationHeader() {
872 return "R Rfree";
873 }
874
875
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
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
898 public void printStats() {
899 int nat = 0;
900 int nnonh = 0;
901
902
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
941 public void scaleBulkFit() {
942 for (int i = 0; i < n; i++) {
943 scaleBulkFit(i);
944 }
945 }
946
947
948
949
950
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
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
973 crystalReciprocalSpacesFc[i].computeDensity(refinementData[i].fc);
974 if (solventModel != SolventModel.NONE) {
975 crystalReciprocalSpacesFs[i].computeDensity(refinementData[i].fs);
976 }
977
978
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
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
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
1006
1007
1008
1009
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
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
1046
1047
1048
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
1062
1063
1064
1065
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
1079
1080
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
1094
1095
1096
1097
1098 private void writeMaps(String filename, int i) {
1099 if (!scaled[i]) {
1100 scaleBulkFit(i);
1101 }
1102
1103
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
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
1136
1137
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
1175
1176
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
1217
1218
1219
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
1233
1234
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
1248
1249
1250
1251
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
1262
1263
1264
1265
1266
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
1285
1286
1287
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
1299
1300
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
1311
1312
1313
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
1342
1343
1344
1345
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
1362
1363
1364
1365 double getbSimWeight() {
1366 return bSimWeight;
1367 }
1368
1369
1370
1371
1372
1373
1374 double getbNonZeroWeight() {
1375 return bNonZeroWeight;
1376 }
1377
1378
1379
1380
1381
1382
1383 double getbMass() {
1384 return bMass;
1385 }
1386
1387
1388
1389
1390
1391
1392 boolean isResidueBFactor() {
1393 return residueBFactor;
1394 }
1395
1396
1397
1398
1399
1400
1401 int getnResidueBFactor() {
1402 return nResidueBFactor;
1403 }
1404
1405
1406
1407
1408
1409
1410 boolean isAddAnisou() {
1411 return addAnisou;
1412 }
1413
1414
1415
1416
1417
1418
1419 double getOccMass() {
1420 return occMass;
1421 }
1422 }