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