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