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 ffx.crystal.Crystal;
41 import ffx.crystal.HKL;
42 import ffx.crystal.ReflectionList;
43 import ffx.crystal.ReflectionSpline;
44 import ffx.numerics.math.ComplexNumber;
45 import ffx.xray.solvent.SolventModel;
46
47 import java.util.logging.Logger;
48
49 import static ffx.numerics.math.ScalarMath.u2b;
50 import static java.lang.Double.isNaN;
51 import static java.lang.String.format;
52 import static org.apache.commons.math3.util.FastMath.abs;
53 import static org.apache.commons.math3.util.FastMath.exp;
54 import static org.apache.commons.math3.util.FastMath.pow;
55 import static org.apache.commons.math3.util.FastMath.sqrt;
56
57
58
59
60
61
62
63 public class CrystalStats {
64
65 private static final Logger logger = Logger.getLogger(CrystalStats.class.getName());
66 private final ReflectionList reflectionList;
67 private final DiffractionRefinementData refinementData;
68 private final Crystal crystal;
69 private final int nBins;
70 private final double[][] fSigF;
71 private final int[] freeR;
72 private final double[][] fcTot;
73 private final double[][] fomPhi;
74 private final ReflectionSpline spline;
75
76 private double blowDPI = -1.0;
77 private boolean print = true;
78
79 private int nObsHKL, highnObsHKL, nObsRFree, highnObsRFree;
80 private double resHigh, resLow, highResHigh, highResLow;
81 private double completeness, highCompleteness;
82 private double rall, r, rfree, highR, highRfree;
83 private double blowDPIH, cruickDPI, cruickDPIH;
84
85
86
87
88
89
90
91 CrystalStats(ReflectionList reflectionList, DiffractionRefinementData refinementData) {
92 this.reflectionList = reflectionList;
93 this.refinementData = refinementData;
94 this.crystal = reflectionList.crystal;
95 this.nBins = refinementData.nBins;
96 this.fSigF = refinementData.fSigF;
97 this.freeR = refinementData.freeR;
98 this.fcTot = refinementData.fcTot;
99 this.fomPhi = refinementData.fomPhi;
100 this.spline = new ReflectionSpline(reflectionList, refinementData.spline.length);
101 }
102
103
104
105
106
107
108 public double getR() {
109 double numer;
110 double denom;
111 double sum = 0.0;
112 double sumfo = 0.0;
113 double sumall = 0.0;
114 double sumfoall = 0.0;
115 for (HKL ih : reflectionList.hklList) {
116 int i = ih.getIndex();
117
118
119 if (isNaN(fcTot[i][0]) || isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
120 continue;
121 }
122
123
124 double ss = crystal.invressq(ih);
125 double fh = spline.f(ss, refinementData.spline);
126
127 ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
128 numer = abs(abs(fSigF[i][0]) - fh * abs(c.abs()));
129 denom = abs(fSigF[i][0]);
130 sumall += numer;
131 sumfoall += denom;
132 if (!refinementData.isFreeR(i)) {
133 sum += numer;
134 sumfo += denom;
135 }
136 }
137
138 rall = (sumall / sumfoall) * 100.0;
139 r = (sum / sumfo) * 100.0;
140 return r;
141 }
142
143
144
145
146
147
148 public double getSigmaA() {
149 double sum = 0.0;
150 int nhkl = 0;
151 ReflectionSpline sigmaaspline =
152 new ReflectionSpline(reflectionList, refinementData.sigmaA.length);
153
154 for (HKL ih : reflectionList.hklList) {
155 int i = ih.getIndex();
156
157
158 if (isNaN(fcTot[i][0]) || isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
159 continue;
160 }
161
162
163 double ss = crystal.invressq(ih);
164 spline.f(ss, refinementData.spline);
165 double sa = sigmaaspline.f(ss, refinementData.sigmaA);
166
167 nhkl++;
168 sum += (sa - sum) / nhkl;
169 }
170
171 return sum;
172 }
173
174
175
176
177
178
179 String getPDBHeaderString() {
180 print = false;
181 printHKLStats();
182 printRStats();
183 print = true;
184
185 StringBuilder sb = new StringBuilder();
186
187 sb.append("REMARK 3 DATA USED IN REFINEMENT\n");
188 sb.append(format("REMARK 3 RESOLUTION RANGE HIGH (ANGSTROMS) : %6.2f\n", resHigh));
189 sb.append(format("REMARK 3 RESOLUTION RANGE LOW (ANGSTROMS) : %6.2f\n", resLow));
190 if (refinementData.fSigFCutoff > 0.0) {
191 sb.append(
192 format(
193 "REMARK 3 DATA CUTOFF (SIGMA(F)) : %6.2f\n",
194 refinementData.fSigFCutoff));
195 } else {
196 sb.append("REMARK 3 DATA CUTOFF (SIGMA(F)) : NONE\n");
197 }
198 sb.append(format("REMARK 3 COMPLETENESS FOR RANGE (%%) : %6.2f\n", completeness));
199 sb.append("REMARK 3 NUMBER OF REFLECTIONS : ").append(nObsHKL).append("\n");
200 sb.append("REMARK 3\n");
201 sb.append("REMARK 3 FIT TO DATA USED IN REFINEMENT\n");
202 sb.append("REMARK 3 CROSS-VALIDATION METHOD : THROUGHOUT\n");
203 sb.append("REMARK 3 FREE R VALUE TEST SET SELECTION : RANDOM\n");
204 sb.append(format("REMARK 3 R VALUE (OBSERVED) : %8.6f\n", rall / 100.0));
205 sb.append(format("REMARK 3 R VALUE (WORKING SET) : %8.6f\n", r / 100.0));
206 sb.append(format("REMARK 3 FREE R VALUE : %8.6f\n", rfree / 100.0));
207 sb.append(
208 format(
209 "REMARK 3 FREE R VALUE TEST SET SIZE (%%) : %6.1f\n",
210 (((double) nObsRFree) / nObsHKL) * 100.0));
211 sb.append("REMARK 3 FREE R VALUE TEST SET COUNT : ").append(nObsRFree).append("\n");
212 sb.append("REMARK 3\n");
213 sb.append("REMARK 3 FIT IN THE HIGHEST RESOLUTION BIN\n");
214 sb.append("REMARK 3 TOTAL NUMBER OF BINS USED : ").append(nBins).append("\n");
215 sb.append(format("REMARK 3 BIN RESOLUTION RANGE HIGH : %6.2f\n", highResHigh));
216 sb.append(format("REMARK 3 BIN RESOLUTION RANGE LOW : %6.2f\n", highResLow));
217 sb.append("REMARK 3 REFLECTION IN BIN (WORKING SET) : ").append(highnObsHKL).append("\n");
218 sb.append(format("REMARK 3 BIN COMPLETENESS (WORKING+TEST) (%%) : %6.2f\n", highCompleteness));
219 sb.append(format("REMARK 3 BIN R VALUE (WORKING SET) : %8.6f\n", highR / 100.0));
220 sb.append("REMARK 3 BIN FREE R VALUE SET COUNT : ").append(highnObsRFree)
221 .append("\n");
222 sb.append(format("REMARK 3 BIN FREE R VALUE : %8.6f\n", highRfree / 100.0));
223 sb.append("REMARK 3\n");
224
225 sb.append("REMARK 3 OVERALL SCALE FACTORS\n");
226 sb.append(format("REMARK 3 SCALE: %4.2f\n", exp(0.25 * refinementData.modelScaleK)));
227 sb.append("REMARK 3 ANISOTROPIC SCALE TENSOR:\n");
228 sb.append(
229 format(
230 "REMARK 3 %g %g %g\n",
231 refinementData.modelAnisoB[0],
232 refinementData.modelAnisoB[3],
233 refinementData.modelAnisoB[4]));
234 sb.append(
235 format(
236 "REMARK 3 %g %g %g\n",
237 refinementData.modelAnisoB[3],
238 refinementData.modelAnisoB[1],
239 refinementData.modelAnisoB[5]));
240 sb.append(
241 format(
242 "REMARK 3 %g %g %g\n",
243 refinementData.modelAnisoB[4],
244 refinementData.modelAnisoB[5],
245 refinementData.modelAnisoB[2]));
246 sb.append("REMARK 3\n");
247
248 CrystalReciprocalSpace crystalReciprocalSpaceFs = refinementData.crystalReciprocalSpaceFs;
249 SolventModel solventModel = crystalReciprocalSpaceFs.getSolventModel();
250 if (solventModel != SolventModel.NONE) {
251 sb.append("REMARK 3 BULK SOLVENT MODELLING\n");
252 double solventA = crystalReciprocalSpaceFs.getSolventA();
253 double solventB = crystalReciprocalSpaceFs.getSolventB();
254 switch (solventModel) {
255 case BINARY:
256 sb.append("REMARK 3 METHOD USED: BINARY MASK\n");
257 sb.append(format("REMARK 3 PROBE RADIUS : %g\n", solventA));
258 sb.append(format("REMARK 3 SHRINK RADIUS : %g\n", solventB));
259 break;
260 case POLYNOMIAL:
261 sb.append("REMARK 3 METHOD USED: POLYNOMIAL SWITCH\n");
262 sb.append(format("REMARK 3 ATOMIC RADIUS BUFFER : %g\n", solventA));
263 sb.append(format("REMARK 3 SWITCH RADIUS : %g\n", solventB));
264 break;
265 case GAUSSIAN:
266 sb.append("REMARK 3 METHOD USED: GAUSSIAN\n");
267 sb.append(format("REMARK 3 ATOMIC RADIUS BUFFER : %g\n", solventA));
268 sb.append(format("REMARK 3 STD DEV SCALE : %g\n", solventB));
269 break;
270 }
271 sb.append(format("REMARK 3 K_SOL: %g\n", refinementData.bulkSolventK));
272 sb.append(
273 format(
274 "REMARK 3 B_SOL: %g\n",
275 refinementData.bulkSolventUeq * 8.0 * Math.PI * Math.PI));
276 sb.append("REMARK 3\n");
277 }
278
279 if (blowDPI > 0.0) {
280 sb.append("REMARK 3 ERROR ESTIMATES\n");
281 sb.append("REMARK 3 ACTA CRYST (1999) D55, 583-601\n");
282 sb.append("REMARK 3 ACTA CRYST (2002) D58, 792-797\n");
283 sb.append(format("REMARK 3 BLOW DPI ALL ATOMS (EQN 7) : %7.4f\n", blowDPIH));
284 sb.append(format("REMARK 3 BLOW DPI NONH ATOMS (EQN 7) : %7.4f\n", blowDPI));
285 sb.append(format("REMARK 3 CRUICKSHANK DPI ALL ATOMS (EQN 27) : %7.4f\n", cruickDPIH));
286 sb.append(format("REMARK 3 CRUICKSHANK DPI NONH ATOMS (EQN 27) : %7.4f\n", cruickDPI));
287 sb.append("REMARK 3\n");
288 }
289
290 sb.append("REMARK 3 DATA TARGET\n");
291 sb.append("REMARK 3 METHOD USED: MAXIMUM LIKELIHOOD\n");
292 sb.append(format("REMARK 3 -LOG LIKELIHOOD : %g\n", refinementData.llkR));
293 sb.append(format("REMARK 3 -LOG LIKELIHOOD (FREE SET) : %g\n", refinementData.llkF));
294 sb.append("REMARK 3\n");
295
296 return sb.toString();
297 }
298
299
300
301
302
303
304 double getRFree() {
305 double sum = 0.0;
306 double sumfo = 0.0;
307 for (HKL ih : reflectionList.hklList) {
308 int i = ih.getIndex();
309
310
311 if (isNaN(fcTot[i][0]) || isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
312 continue;
313 }
314
315
316 double ss = crystal.invressq(ih);
317 double fh = spline.f(ss, refinementData.spline);
318
319 ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
320 if (refinementData.isFreeR(i)) {
321 sum += abs(abs(fSigF[i][0]) - fh * abs(c.abs()));
322 sumfo += abs(fSigF[i][0]);
323 }
324 }
325
326 rfree = (sum / sumfo) * 100.0;
327 return rfree;
328 }
329
330
331
332
333
334
335 double getSigmaW() {
336 double sum = 0.0;
337 int nhkl = 0;
338 ReflectionSpline sigmaaspline =
339 new ReflectionSpline(reflectionList, refinementData.sigmaA.length);
340
341 for (HKL ih : reflectionList.hklList) {
342 int i = ih.getIndex();
343
344
345 if (isNaN(fcTot[i][0]) || isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
346 continue;
347 }
348
349
350 double ss = crystal.invressq(ih);
351 spline.f(ss, refinementData.spline);
352 double wa = sigmaaspline.f(ss, refinementData.sigmaW);
353
354 nhkl++;
355 sum += (wa - sum) / nhkl;
356 }
357
358 return sum;
359 }
360
361
362
363
364
365
366
367
368
369
370
371 void printDPIStats(int nAtoms, int nNonHydrogenAtoms) {
372 int nhkli = 0;
373 int nhklo = refinementData.n;
374 double rfreefrac = getRFree() * 0.01;
375 double res = reflectionList.resolution.resolutionLimit();
376 for (HKL ih : reflectionList.hklList) {
377 int i = ih.getIndex();
378
379
380 if (isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
381 continue;
382 }
383 nhkli++;
384 }
385
386 double va = pow(crystal.volume / crystal.spaceGroup.getNumberOfSymOps(), 0.3333);
387 blowDPIH = 1.28 * sqrt(nAtoms) * va * pow(nhkli, -0.8333) * rfreefrac;
388 blowDPI = 1.28 * sqrt(nNonHydrogenAtoms) * va * pow(nhkli, -0.8333) * rfreefrac;
389 double natni = sqrt((double) nAtoms / nhkli);
390 double noni = pow((double) nhkli / nhklo, -0.3333);
391 cruickDPIH = natni * noni * rfreefrac * res;
392 natni = sqrt((double) nNonHydrogenAtoms / nhkli);
393 cruickDPI = natni * noni * rfreefrac * res;
394
395 StringBuilder sb = new StringBuilder("\n");
396 sb.append(format(" Blow DPI for all / non-H atoms: %7.4f / %7.4f\n", blowDPIH, blowDPI));
397 sb.append(
398 format(" Cruickshank DPI for all / non-H atoms: %7.4f / %7.4f", cruickDPIH, cruickDPI));
399
400 if (print) {
401 logger.info(sb.toString());
402 }
403 }
404
405
406
407
408 void printHKLStats() {
409 double[][] res = new double[nBins][2];
410 int[][] nhkl = new int[nBins][3];
411
412 for (int i = 0; i < nBins; i++) {
413 res[i][0] = Double.NEGATIVE_INFINITY;
414 res[i][1] = Double.POSITIVE_INFINITY;
415 }
416
417 for (HKL ih : reflectionList.hklList) {
418 int i = ih.getIndex();
419 int b = ih.getBin();
420
421
422 if (isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
423 nhkl[b][2]++;
424 continue;
425 }
426
427
428 double rh = crystal.res(ih);
429 if (rh > res[b][0]) {
430 res[b][0] = rh;
431 }
432 if (rh < res[b][1]) {
433 res[b][1] = rh;
434 }
435
436
437 if (freeR[i] == refinementData.rFreeFlag) {
438 nhkl[b][1]++;
439 } else {
440 nhkl[b][0]++;
441 }
442 }
443
444 StringBuilder sb =
445 new StringBuilder(
446 format(
447 "\n %15s | %8s|%9s| %7s | %7s | %s\n",
448 "Res. Range", " HKL (R)", " HKL (cv)", " Bin", " Miss", "Complete (%)"));
449 for (int i = 0; i < nBins; i++) {
450 sb.append(format(" %7.3f %7.3f | ", res[i][0], res[i][1]));
451 sb.append(
452 format(
453 "%7d | %7d | %7d | %7d | ",
454 nhkl[i][0], nhkl[i][1], nhkl[i][0] + nhkl[i][1], nhkl[i][2]));
455 sb.append(
456 format(
457 "%6.2f\n",
458 (((double) nhkl[i][0] + nhkl[i][1]) / (nhkl[i][0] + nhkl[i][1] + nhkl[i][2]))
459 * 100.0));
460 }
461 sb.append(format(" %7.3f %7.3f | ", res[0][0], res[nBins - 1][1]));
462 int sum1 = 0;
463 int sum2 = 0;
464 int sum3 = 0;
465 for (int i = 0; i < nBins; i++) {
466 sum1 += nhkl[i][0];
467 sum2 += nhkl[i][1];
468 sum3 += nhkl[i][2];
469 }
470 sb.append(format("%7d | %7d | %7d | %7d | ", sum1, sum2, sum1 + sum2, sum3));
471 sb.append(format("%6.2f\n", (((double) sum1 + sum2) / (sum1 + sum2 + sum3)) * 100.0));
472 sb.append(format(" Number of reflections if complete: %10d", refinementData.n));
473
474 nObsHKL = sum1 + sum2;
475 highnObsHKL = nhkl[nBins - 1][0] + nhkl[nBins - 1][1];
476 nObsRFree = sum2;
477 highnObsRFree = nhkl[nBins - 1][1];
478 completeness = (((double) sum1 + sum2) / (sum1 + sum2 + sum3)) * 100.0;
479 highCompleteness =
480 (((double) nhkl[nBins - 1][0] + nhkl[nBins - 1][1])
481 / (nhkl[nBins - 1][0] + nhkl[nBins - 1][1] + nhkl[nBins - 1][2]))
482 * 100.0;
483
484 if (print) {
485 logger.info(sb.toString());
486 }
487 }
488
489
490
491
492 void printRStats() {
493 double[][] res = new double[nBins][2];
494 double[] nhkl = new double[nBins + 1];
495 double[][] rb = new double[nBins + 1][2];
496 double[][] sumfo = new double[nBins + 1][2];
497 double[][] s = new double[nBins + 1][4];
498 double numer;
499 double denom;
500 double sumall = 0.0;
501 double sumfoall = 0.0;
502 ReflectionSpline sigmaASpline =
503 new ReflectionSpline(reflectionList, refinementData.sigmaA.length);
504
505 for (int i = 0; i < nBins; i++) {
506 res[i][0] = Double.NEGATIVE_INFINITY;
507 res[i][1] = Double.POSITIVE_INFINITY;
508 }
509
510 for (HKL ih : reflectionList.hklList) {
511 int i = ih.getIndex();
512 int b = ih.getBin();
513
514
515 if (isNaN(fcTot[i][0]) || isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
516 continue;
517 }
518
519
520 double ss = crystal.invressq(ih);
521 double fh = spline.f(ss, refinementData.spline);
522 double sa = sigmaASpline.f(ss, refinementData.sigmaA);
523 double wa = sigmaASpline.f(ss, refinementData.sigmaW);
524 double eoscale = sigmaASpline.f(ss, refinementData.esqFo);
525
526
527 double rs = crystal.res(ih);
528 if (rs > res[b][0]) {
529 res[b][0] = rs;
530 }
531 if (rs < res[b][1]) {
532 res[b][1] = rs;
533 }
534
535 ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
536 numer = abs(abs(fSigF[i][0]) - fh * abs(c.abs()));
537 denom = abs(fSigF[i][0]);
538 if (refinementData.isFreeR(i)) {
539 rb[b][1] += numer;
540 sumfo[b][1] += denom;
541 rb[nBins][1] += numer;
542 sumfo[nBins][1] += denom;
543 sumall += numer;
544 sumfoall += denom;
545 } else {
546 rb[b][0] += numer;
547 sumfo[b][0] += denom;
548 rb[nBins][0] += numer;
549 sumfo[nBins][0] += denom;
550 sumall += numer;
551 sumfoall += denom;
552 }
553
554 nhkl[b]++;
555 nhkl[nBins]++;
556 s[b][0] += (sa - s[b][0]) / nhkl[b];
557 s[b][1] += (wa - s[b][1]) / nhkl[b];
558 s[b][2] += ((wa / sqrt(eoscale)) - s[b][2]) / nhkl[b];
559 s[b][3] += (fomPhi[i][0] - s[b][3]) / nhkl[b];
560 s[nBins][0] += (sa - s[nBins][0]) / nhkl[nBins];
561 s[nBins][1] += (wa - s[nBins][1]) / nhkl[nBins];
562 s[nBins][2] += ((wa / Math.sqrt(eoscale)) - s[nBins][2]) / nhkl[nBins];
563 s[nBins][3] += (fomPhi[i][0] - s[nBins][3]) / nhkl[nBins];
564 }
565
566 StringBuilder sb =
567 new StringBuilder(format("\n %15s | %7s | %7s | %7s | %7s | %7s | %7s\n",
568 "Res. Range", " R", "Rfree", "s", "w(E)", "w(F)", "FOM"));
569 for (int i = 0; i < nBins; i++) {
570 sb.append(format(" %7.3f %7.3f | ", res[i][0], res[i][1]));
571 sb.append(format("%7.2f | %7.2f | %7.4f | %7.4f | %7.2f | %7.4f\n",
572 (rb[i][0] / sumfo[i][0]) * 100.0,
573 (rb[i][1] / sumfo[i][1]) * 100.0,
574 s[i][0], s[i][1], s[i][2], s[i][3]));
575 }
576
577 sb.append(format(" %7.3f %7.3f | ", res[0][0], res[nBins - 1][1]));
578 sb.append(format("%7.2f | %7.2f | %7.4f | %7.4f | %7.2f | %7.4f\n",
579 (rb[nBins][0] / sumfo[nBins][0]) * 100.0,
580 (rb[nBins][1] / sumfo[nBins][1]) * 100.0,
581 s[nBins][0], s[nBins][1], s[nBins][2], s[nBins][3]));
582
583 resLow = res[0][0];
584 resHigh = res[nBins - 1][1];
585 highResLow = res[nBins - 1][0];
586 highResHigh = res[nBins - 1][1];
587 r = (rb[nBins][0] / sumfo[nBins][0]) * 100.0;
588 rfree = (rb[nBins][1] / sumfo[nBins][1]) * 100.0;
589 rall = (sumall / sumfoall) * 100.0;
590 highR = (rb[nBins - 1][0] / sumfo[nBins - 1][0]) * 100.0;
591 highRfree = (rb[nBins - 1][1] / sumfo[nBins - 1][1]) * 100.0;
592
593 if (print) {
594 logger.info(sb.toString());
595 }
596 }
597
598
599
600
601 void printScaleStats() {
602 int[] nhkl = new int[nBins];
603 double[] scale = new double[nBins];
604
605 for (HKL ih : reflectionList.hklList) {
606 int i = ih.getIndex();
607 int b = ih.getBin();
608
609
610 if (isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
611 continue;
612 }
613
614
615 double ss = crystal.invressq(ih);
616 double fh = spline.f(ss, refinementData.spline);
617
618 nhkl[b]++;
619 scale[b] += (fh - scale[b]) / nhkl[b];
620 }
621
622 StringBuilder sb = new StringBuilder(
623 format(" K_Overall (Fc to Fo Scale): %4.2f\n", exp(0.25 * refinementData.modelScaleK)));
624 sb.append(" Fc to Fo Spline Scale: ");
625 for (int i = 0; i < nBins; i++) {
626 sb.append(format("%4.2f ", scale[i]));
627 }
628 sb.append("\n B_Overall Anisotropic Tensor:\n");
629 sb.append(format(" %10.4f %10.4f %10.4f\n", refinementData.modelAnisoB[0],
630 refinementData.modelAnisoB[3], refinementData.modelAnisoB[4]));
631 sb.append(format(" %10.4f %10.4f %10.4f\n", refinementData.modelAnisoB[3],
632 refinementData.modelAnisoB[1], refinementData.modelAnisoB[5]));
633 sb.append(format(" %10.4f %10.4f %10.4f\n", refinementData.modelAnisoB[4],
634 refinementData.modelAnisoB[5], refinementData.modelAnisoB[2]));
635
636 CrystalReciprocalSpace crystalReciprocalSpaceFs = refinementData.crystalReciprocalSpaceFs;
637 SolventModel solventModel = crystalReciprocalSpaceFs.getSolventModel();
638 if (solventModel != SolventModel.NONE) {
639 double solventA = crystalReciprocalSpaceFs.getSolventA();
640 double solventB = crystalReciprocalSpaceFs.getSolventB();
641 switch (solventModel) {
642 case BINARY -> {
643 sb.append(" Bulk Solvent Model: Binary Mask\n");
644 sb.append(format(" Probe Radius: %8.4f\n", solventA));
645 sb.append(format(" Shrink Radius: %8.4f\n", solventB));
646 }
647 case POLYNOMIAL -> {
648 sb.append(" Bulk Solvent Model: Polynomial Switch\n");
649 sb.append(format(" Probe Radius: %8.4f\n", solventA));
650 sb.append(format(" Window Width: %8.4f\n", solventB));
651 }
652 case GAUSSIAN -> {
653 sb.append(" Bulk Solvent Model: Gaussian\n");
654 sb.append(format(" Exponential Scale A: %8.4f\n", solventA));
655 sb.append(format(" Radius Scale Factor: %8.4f\n", solventB));
656 }
657 case NONE -> sb.append(" Bulk Solvent Model: None\n");
658 }
659 sb.append(format(" Scattering Density (k): %8.4f\n", refinementData.bulkSolventK));
660 sb.append(format(" B-Factor: %8.4f\n", u2b(refinementData.bulkSolventUeq)));
661 }
662 sb.append(format("\n -log Likelihood: %14.3f (free set: %14.3f)", refinementData.llkR, refinementData.llkF));
663
664 if (print) {
665 logger.info(sb.toString());
666 }
667 }
668
669
670
671
672 void printSNStats() {
673 double[][] res = new double[nBins][2];
674 double[] nhkl = new double[nBins + 1];
675 double[][] sn = new double[nBins + 1][3];
676
677 for (int i = 0; i < nBins; i++) {
678 res[i][0] = Double.NEGATIVE_INFINITY;
679 res[i][1] = Double.POSITIVE_INFINITY;
680 }
681
682 for (HKL ih : reflectionList.hklList) {
683 int i = ih.getIndex();
684 int b = ih.getBin();
685
686
687 if (isNaN(fSigF[i][0]) || fSigF[i][1] <= 0.0) {
688 continue;
689 }
690
691
692 double rs = crystal.res(ih);
693 if (rs > res[b][0]) {
694 res[b][0] = rs;
695 }
696 if (rs < res[b][1]) {
697 res[b][1] = rs;
698 }
699
700
701 nhkl[b]++;
702 nhkl[nBins]++;
703 sn[b][0] += (fSigF[i][0] - sn[b][0]) / nhkl[b];
704 sn[b][1] += (fSigF[i][1] - sn[b][1]) / nhkl[b];
705 sn[b][2] += ((fSigF[i][0] / fSigF[i][1]) - sn[b][2]) / nhkl[b];
706
707 sn[nBins][0] += (fSigF[i][0] - sn[nBins][0]) / nhkl[b];
708 sn[nBins][1] += (fSigF[i][1] - sn[nBins][1]) / nhkl[b];
709 sn[nBins][2] += ((fSigF[i][0] / fSigF[i][1]) - sn[nBins][2]) / nhkl[nBins];
710 }
711
712 StringBuilder sb =
713 new StringBuilder(
714 format("\n %15s | %7s | %7s | %7s \n", "Res. Range", "Signal", "Sigma", "S/N"));
715 for (int i = 0; i < nBins; i++) {
716 sb.append(format(" %7.3f %7.3f | ", res[i][0], res[i][1]));
717 sb.append(format("%7.2f | %7.2f | %7.2f\n", sn[i][0], sn[i][1], sn[i][2]));
718 }
719
720 sb.append(format(" %7.3f %7.3f | ", res[0][0], res[nBins - 1][1]));
721 sb.append(format("%7.2f | %7.2f | %7.2f", sn[nBins][0], sn[nBins][1], sn[nBins][2]));
722
723 if (print) {
724 logger.info(sb.toString());
725 }
726 }
727 }