partialator: Add --max-rel-B

4 files changed, 26 insertions, 6 deletions

diff --git a/doc/man/partialator.1 b/doc/man/partialator.1
index 11e535c7..640376d9 100644
--- a/doc/man/partialator.1
+++ b/doc/man/partialator.1
@@ -122,6 +122,11 @@ Disable cross-validation (for testing only).
 .PD
 Read a list of filenames, event IDs and dataset IDs from \fIfilename\fR.  See the section on \fBCUSTOM DATASET SPLITTING\fR below.
 
+.PD 0
+.IP \fB--max-rel-B=\fIn\fR
+.PD
+Reject crystals if the absolute values of their relative Debye-Waller ("B") factors are more than \fIn\fR square Angstroms.  The default is \fB--max-rel-B=100\fR.
+
 .SH PARTIALITY MODELS
 
 The available partiality models are:
diff --git a/src/partialator.c b/src/partialator.c
index d6606ef4..079d4d46 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -310,7 +310,8 @@ static void show_help(const char *s)
 "      --push-res=<n>         Merge higher than apparent resolution cutoff.\n"
 "  -j <n>                     Run <n> analyses in parallel.\n"
 "      --no-free              Disable cross-validation (testing only).\n"
-"      --custom-split         List of files for custom dataset splitting.\n");
+"      --custom-split         List of files for custom dataset splitting.\n"
+"      --max-rel-B            Maximum allowable relative |B| factor.\n");
 }
 
 
@@ -620,6 +621,7 @@ int main(int argc, char *argv[])
 	int no_free = 0;
 	char *csplit_fn = NULL;
 	struct custom_split *csplit = NULL;
+	double max_B = 1e-18;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -640,6 +642,8 @@ int main(int argc, char *argv[])
 		{"push-res",           1, NULL,                5},
 		{"res-push",           1, NULL,                5}, /* compat */
 		{"custom-split",       1, NULL,                6},
+		{"max-rel-B"   ,       1, NULL,                7},
+		{"max-rel-b"   ,       1, NULL,                7}, /* compat */
 
 		{"no-scale",           0, &no_scale,           1},
 		{"no-pr",              0, &no_pr,              1},
@@ -754,6 +758,16 @@ int main(int argc, char *argv[])
 			csplit_fn = strdup(optarg);
 			break;
 
+			case 7 :
+			errno = 0;
+			max_B = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --max-rel-B.\n");
+				return 1;
+			}
+			max_B = max_B * 1e-20;
+			break;
+
 			case 0 :
 			break;
 
@@ -969,7 +983,7 @@ int main(int argc, char *argv[])
 	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
 	                         min_measurements, push_res);
 
-	check_rejection(crystals, n_crystals, full);
+	check_rejection(crystals, n_crystals, full, max_B);
 
 	write_pgraph(full, crystals, n_crystals, 0);
 
@@ -991,7 +1005,7 @@ int main(int argc, char *argv[])
 		STATUS("Overall free residual: initial = %e, final = %e\n",
 		       init_free_dev, final_free_dev);
 
-		check_rejection(crystals, n_crystals, full);
+		check_rejection(crystals, n_crystals, full, max_B);
 		normalise_scales(crystals, n_crystals);
 
 		/* Re-estimate all the full intensities */
diff --git a/src/rejection.c b/src/rejection.c
index cea43701..ef185034 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -193,7 +193,7 @@ static void show_duds(Crystal **crystals, int n_crystals)
 }
 
 
-void check_rejection(Crystal **crystals, int n, RefList *full)
+void check_rejection(Crystal **crystals, int n, RefList *full, double max_B)
 {
 	int i;
 	int n_acc = 0;
@@ -204,7 +204,7 @@ void check_rejection(Crystal **crystals, int n, RefList *full)
 	for ( i=0; i<n; i++ ) {
 
 		/* Reject if B factor modulus is very large */
-		if ( fabs(crystal_get_Bfac(crystals[i])) > 1e-18 ) {
+		if ( fabs(crystal_get_Bfac(crystals[i])) > max_B ) {
 			crystal_set_user_flag(crystals[i], PRFLAG_BIGB);
 		}
 
diff --git a/src/rejection.h b/src/rejection.h
index ec529941..8979313b 100644
--- a/src/rejection.h
+++ b/src/rejection.h
@@ -38,6 +38,7 @@
 #include "crystal.h"
 
 extern void early_rejection(Crystal **crystals, int n);
-extern void check_rejection(Crystal **crystals, int n, RefList *full);
+extern void check_rejection(Crystal **crystals, int n, RefList *full,
+                            double max_B);
 
 #endif	/* REJECTION_H */