From 70cb465c18a8b14796607fa32b2d945af500560b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 24 Oct 2019 15:30:27 +0200 Subject: get_hkl: Clean up resolution cutoffs and add --lowres --- doc/man/get_hkl.1 | 3 ++ src/get_hkl.c | 101 ++++++++++++++++++++++++++++-------------------------- 2 files changed, 56 insertions(+), 48 deletions(-) diff --git a/doc/man/get_hkl.1 b/doc/man/get_hkl.1 index 90637a8a..8d76c154 100644 --- a/doc/man/get_hkl.1 +++ b/doc/man/get_hkl.1 @@ -79,9 +79,12 @@ The intensities of the reflections will be multiplied by their symmetric multipl .PD 0 .IP \fB--cutoff-angstroms=\fR\fIn\fR .IP \fB--cutoff-angstroms=\fR\fIn1,n2,n3\fR +.IP \fB--highres=\fR\fIn\fR +.IP \fB--lowres=\fR\fIn\fR .PD In the first form, reflections with d (=lamba/2*sin(theta)) < \fIn\fR will be removed. In the second form, anisotropic truncation will be performed with separate resolution limits \fIn1\fR, \fIn2\fR and \fIn3\fR along a*, b* and c* respectively. You must also specify \fB-p\fR or \fB--pdb\fR. +The option \fB--lowres\fR will remove reflections with d > \fIn\fR. \fB--highres\fR is a synonym for \fB--cutoff-angstroms\fR. .SH REINDEXING REFLECTIONS .PD 0 diff --git a/src/get_hkl.c b/src/get_hkl.c index 2a17aab5..057a5ce7 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -84,6 +84,8 @@ static void show_help(const char *s) "You can restrict which reflections are written out:\n" " -t, --template= Only include reflections mentioned in file.\n" " --cutoff-angstroms= Only include reflections with d < n Angstroms.\n" +" --highres= Synonym for --cutoff-angstroms.\n" +" --lowres= Only include reflections with d > n Angstroms.\n" "\n" "You might sometimes need to do this:\n" " --multiplicity Multiply intensities by the number of\n" @@ -404,6 +406,37 @@ static RefList *trim_centrics(RefList *in, const SymOpList *sym) } +static RefList *apply_resolution_cutoff(RefList *input, double lowres, + double highres, UnitCell *cell) +{ + RefList *n; + Reflection *refl; + RefListIterator *iter; + + n = reflist_new(); + copy_notes(n, input); + + for ( refl = first_refl(input, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + double res; + get_indices(refl, &h, &k, &l); + res = 2.0 * resolution(cell, h, k, l); + if ( (res < 1e10/highres) && (res > 1e10/lowres) ) { + Reflection *a; + a = add_refl(n, h, k, l); + copy_data(a, refl); + } + } + + cell_free(cell); + reflist_free(input); + return n; +} + + int main(int argc, char *argv[]) { int c; @@ -425,14 +458,15 @@ int main(int argc, char *argv[]) RefList *input; double adu_per_photon = 0.0; int have_adu_per_photon = 0; - int have_cutoff_iso = 0; int have_cutoff_aniso = 0; char *cutoff_str = NULL; - double cutiso = 0.0; float cutn1, cutn2, cutn3; char *cellfile = NULL; char *reindex_str = NULL; SymOpList *reindex = NULL; + float lowres = 0.0; + double highres = INFINITY; + UnitCell *cell = NULL; /* Long options */ const struct option longopts[] = { @@ -452,7 +486,9 @@ int main(int argc, char *argv[]) {"no-need-all-parts", 0, &config_nap, 0}, {"adu-per-photon", 1, NULL, 2}, {"cutoff-angstroms", 1, NULL, 3}, + {"highres", 1, NULL, 3}, {"reindex", 1, NULL, 4}, + {"lowres", 1, NULL, 6}, {0, 0, NULL, 0} }; @@ -512,6 +548,13 @@ int main(int argc, char *argv[]) reindex_str = strdup(optarg); break; + case 6 : + if (sscanf(optarg, "%f", &lowres) != 1) { + ERROR("Invalid value for --lowres\n"); + return 1; + } + break; + case 0 : break; @@ -538,19 +581,20 @@ int main(int argc, char *argv[]) /* Convert Angstroms -> m */ cutn1 /= 1e10; cutn2 /= 1e10; cutn3 /= 1e10; + /* Set isotropic cutoff to the largest */ + highres = biggest(cutn1, biggest(cutn2, cutn3)); + } else { char *rval; errno = 0; - cutiso = strtod(cutoff_str, &rval); + highres = strtod(cutoff_str, &rval); if ( *rval != '\0' ) { ERROR("Invalid value for --cutoff-angstroms.\n"); return 1; } - have_cutoff_iso = 1; - } free(cutoff_str); @@ -729,16 +773,11 @@ int main(int argc, char *argv[]) } - if ( have_cutoff_iso ) { - - RefList *n; - Reflection *refl; - RefListIterator *iter; - UnitCell *cell; + if ( (highres < INFINITY) || have_cutoff_aniso || (lowres > 0.0) ) { if ( cellfile == NULL ) { ERROR("You must provide a unit cell when using " - "--cutoff-angstroms.\n"); + "--cutoff-angstroms, --highres or --lowres.\n"); return 1; } @@ -749,54 +788,20 @@ int main(int argc, char *argv[]) } free(cellfile); - n = reflist_new(); - copy_notes(n, input); - - for ( refl = first_refl(input, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - double res; - get_indices(refl, &h, &k, &l); - res = 2.0 * resolution(cell, h, k, l); - if ( res < 1e10 / cutiso ) { - Reflection *a; - a = add_refl(n, h, k, l); - copy_data(a, refl); - } - } - - cell_free(cell); - reflist_free(input); - input = n; - } + input = apply_resolution_cutoff(input, lowres, highres, cell); + if ( have_cutoff_aniso ) { RefList *n; Reflection *refl; RefListIterator *iter; - UnitCell *cell; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double as, bs, cs; - if ( cellfile == NULL ) { - ERROR("You must provide a unit cell when using " - "--cutoff-angstroms.\n"); - return 1; - } - - cell = load_cell_from_file(cellfile); - if ( cell == NULL ) { - ERROR("Failed to load cell from '%s'\n", cellfile); - return 1; - } - free(cellfile); - cell_get_reciprocal(cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); -- cgit v1.2.3