diff options
-rw-r--r-- | doc/man/partialator.1 | 18 | ||||
-rw-r--r-- | src/partialator.c | 61 | ||||
-rw-r--r-- | src/post-refinement.c | 105 | ||||
-rw-r--r-- | src/post-refinement.h | 4 |
4 files changed, 128 insertions, 60 deletions
diff --git a/doc/man/partialator.1 b/doc/man/partialator.1 index 980af00f..1a598e89 100644 --- a/doc/man/partialator.1 +++ b/doc/man/partialator.1 @@ -1,7 +1,7 @@ .\" .\" partialator man page .\" -.\" Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, +.\" Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, .\" a research centre of the Helmholtz Association. .\" .\" Part of CrystFEL - crystallography with a FEL @@ -137,6 +137,20 @@ Write out the per-crystal parameters and reflection lists after every cycle of r .PD Do not write the extensive log files needed for plotting contour maps and spectrum graphs. This makes the process a lot faster, but you probably do want these logs to check that post-refinement is working reasonably. +.PD 0 +.IP "\fB-w\fR \fIpg\fR" +.PD +Set the apparent symmetry of the crystals. The ambiguity operator will be determined by comparing this to the actual symmetry. +.IP +If you prefer (or the scenario demands it), you can specify the ambiguity operator directly using \fB--operator\fR. + +.PD 0 +.IP \fB--operator=\fR\fIop\fR +.PD +Specify the indexing ambiguity operator. Example: \fB--operator=k,h,-l\fR. +.IP +If you prefer, you can specify the ambiguity operator by specifying the apparent symmetry using \fB-w\fR. + .SH PARTIALITY MODELS The available partiality models are: @@ -221,7 +235,7 @@ This page was written by Thomas White. Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>. .SH COPYRIGHT AND DISCLAIMER -Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. +Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. .P partialator, and this manual, are part of CrystFEL. .P diff --git a/src/partialator.c b/src/partialator.c index 5c751732..baab84c4 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -327,7 +327,9 @@ static void show_help(const char *s) " --no-free Disable cross-validation (testing only).\n" " --custom-split List of files for custom dataset splitting.\n" " --max-rel-B Maximum allowable relative |B| factor.\n" -" --no-logs Do not write extensive log files.\n"); +" --no-logs Do not write extensive log files.\n" +" -w <pg> Apparent point group for resolving ambiguities.\n" +" --operator=<op> Indexing ambiguity operator for resolving.\n"); } @@ -820,6 +822,8 @@ int main(int argc, char *argv[]) char *outfile = NULL; char *sym_str = NULL; SymOpList *sym; + SymOpList *amb; + SymOpList *w_sym; int nthreads = 1; int i; int n_iter = 10; @@ -852,6 +856,8 @@ int main(int argc, char *argv[]) char *rfile = NULL; RefList *reference = NULL; int no_logs = 0; + char *w_sym_str = NULL; + char *operator = NULL; /* Long options */ const struct option longopts[] = { @@ -875,6 +881,7 @@ int main(int argc, char *argv[]) {"max-rel-B", 1, NULL, 7}, {"max-rel-b", 1, NULL, 7}, /* compat */ {"reference", 1, NULL, 8}, /* ssshhh! */ + {"operator", 1, NULL, 9}, {"no-scale", 0, &no_scale, 1}, {"no-pr", 0, &no_pr, 1}, @@ -896,7 +903,7 @@ int main(int argc, char *argv[]) } /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:", + while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:", longopts, NULL)) != -1) { @@ -955,6 +962,10 @@ int main(int argc, char *argv[]) pmodel_str = strdup(optarg); break; + case 'w' : + w_sym_str = strdup(optarg); + break; + case 2 : errno = 0; min_measurements = strtod(optarg, &rval); @@ -1005,6 +1016,10 @@ int main(int argc, char *argv[]) rfile = strdup(optarg); break; + case 9 : + operator = strdup(optarg); + break; + case 0 : break; @@ -1044,6 +1059,44 @@ int main(int argc, char *argv[]) free(sym_str); if ( sym == NULL ) return 1; + if ( (w_sym_str != NULL) && (operator != NULL) ) { + ERROR("Specify the apparent symmetry (-w) or the operator, " + "not both.\n"); + return 1; + } + + if ( w_sym_str == NULL ) { + w_sym = NULL; + amb = NULL; + } else { + pointgroup_warning(w_sym_str); + w_sym = get_pointgroup(w_sym_str); + free(w_sym_str); + if ( w_sym == NULL ) return 1; + amb = get_ambiguities(w_sym, sym); + if ( amb == NULL ) { + ERROR("Couldn't find ambiguity operator.\n"); + ERROR("Check that your values for -y and -w are " + "correct.\n"); + return 1; + } + + } + + if ( operator ) { + amb = parse_symmetry_operations(operator); + if ( amb == NULL ) return 1; + set_symmetry_name(amb, "Ambiguity"); + } + + if ( amb != NULL ) { + STATUS("Indexing ambiguity resolution enabled. " + "The ambiguity operation(s) are:\n"); + describe_symmetry(amb); + /* In contrast to ambigator, partialator can deal with multiple + * ambiguities at once */ + } + if ( pmodel_str != NULL ) { if ( strcmp(pmodel_str, "unity") == 0 ) { pmodel = PMODEL_UNITY; @@ -1264,7 +1317,7 @@ int main(int argc, char *argv[]) if ( !no_pr ) { refine_all(crystals, n_crystals, full, nthreads, pmodel, - 0, i+1, no_logs); + 0, i+1, no_logs, sym, amb); } else if ( !no_scale ) { scale_all_to_reference(crystals, n_crystals, full, nthreads); } @@ -1357,7 +1410,7 @@ int main(int argc, char *argv[]) crystal_free(crystals[i]); } reflist_free(full); - free(sym); + free_symoplist(sym); free(outfile); free(crystals); free(infile); diff --git a/src/post-refinement.c b/src/post-refinement.c index 1c33e211..7e5420d3 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -389,7 +389,7 @@ static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial, static RefList *reindex_reflections(RefList *input, SymOpList *amb, - SymOpList *sym) + SymOpList *sym, int idx) { RefList *n; Reflection *refl; @@ -405,13 +405,13 @@ static RefList *reindex_reflections(RefList *input, SymOpList *amb, Reflection *rn; get_indices(refl, &h, &k, &l); - get_equiv(amb, NULL, 0, h, k, l, &h, &k, &l); + get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l); get_asymm(sym, h, k, l, &h, &k, &l); rn = add_refl(n, h, k, l); copy_data(rn, refl); get_symmetric_indices(rn, &h, &k, &l); - get_equiv(amb, NULL, 0, h, k, l, &h, &k, &l); + get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l); set_symmetric_indices(rn, h, k, l); } @@ -419,7 +419,7 @@ static RefList *reindex_reflections(RefList *input, SymOpList *amb, } -static void reindex_cell(UnitCell *cell, SymOpList *amb) +static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx) { signed int h, k, l; struct rvec na, nb, nc; @@ -429,17 +429,17 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb) &bs.u, &bs.v, &bs.w, &cs.u, &cs.v, &cs.w); - get_equiv(amb, NULL, 0, 1, 0, 0, &h, &k, &l); + get_equiv(amb, NULL, idx, 1, 0, 0, &h, &k, &l); na.u = as.u*h + bs.u*k + cs.u*l; na.v = as.v*h + bs.v*k + cs.v*l; na.w = as.w*h + bs.w*k + cs.w*l; - get_equiv(amb, NULL, 0, 0, 1, 0, &h, &k, &l); + get_equiv(amb, NULL, idx, 0, 1, 0, &h, &k, &l); nb.u = as.u*h + bs.u*k + cs.u*l; nb.v = as.v*h + bs.v*k + cs.v*l; nb.w = as.w*h + bs.w*k + cs.w*l; - get_equiv(amb, NULL, 0, 0, 0, 1, &h, &k, &l); + get_equiv(amb, NULL, idx, 0, 0, 1, &h, &k, &l); nc.u = as.u*h + bs.u*k + cs.u*l; nc.v = as.v*h + bs.v*k + cs.v*l; nc.w = as.w*h + bs.w*k + cs.w*l; @@ -450,45 +450,50 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb) } -void try_reindex(Crystal *crin, const RefList *full) +void try_reindex(Crystal *crin, const RefList *full, + SymOpList *sym, SymOpList *amb) { RefList *list; Crystal *cr; UnitCell *cell; - SymOpList *amb; - SymOpList *sym; double residual_original, residual_flipped; + int idx, n; - amb = parse_symmetry_operations("k,h,-l"); - sym = get_pointgroup("m-3"); + if ( sym == NULL || amb == NULL ) return; residual_original = residual(crin, full, 0, NULL, NULL, 0); cr = crystal_copy(crin); - cell = cell_new_from_cell(crystal_get_cell(crin)); - if ( cell == NULL ) return; - reindex_cell(cell, amb); - crystal_set_cell(cr, cell); + n = num_equivs(amb, NULL); - list = reindex_reflections(crystal_get_reflections(crin), amb, sym); - crystal_set_reflections(cr, list); + for ( idx=0; idx<n; idx++ ) { - update_predictions(cr); - calculate_partialities(cr, PMODEL_XSPHERE); + cell = cell_new_from_cell(crystal_get_cell(crin)); + if ( cell == NULL ) return; + reindex_cell(cell, amb, idx); + crystal_set_cell(cr, cell); - residual_flipped = residual(cr, full, 0, NULL, NULL, 1); + list = reindex_reflections(crystal_get_reflections(crin), + amb, sym, idx); + crystal_set_reflections(cr, list); - if ( residual_flipped < residual_original ) { - crystal_set_cell(crin, cell); - crystal_set_reflections(crin, list); - } else { - cell_free(crystal_get_cell(cr)); - reflist_free(crystal_get_reflections(cr)); - crystal_free(cr); + update_predictions(cr); + calculate_partialities(cr, PMODEL_XSPHERE); + + residual_flipped = residual(cr, full, 0, NULL, NULL, 1); + + if ( residual_flipped < residual_original ) { + crystal_set_cell(crin, cell); + crystal_set_reflections(crin, list); + residual_original = residual_flipped; + } else { + cell_free(crystal_get_cell(cr)); + reflist_free(crystal_get_reflections(cr)); + } } - free_symoplist(amb); + crystal_free(cr); } @@ -686,7 +691,8 @@ void write_gridscan(Crystal *cr, const RefList *full, static void do_pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel, int verbose, int serial, - int cycle, int write_logs) + int cycle, int write_logs, + SymOpList *sym, SymOpList *amb) { gsl_multimin_fminimizer *min; struct rf_priv priv; @@ -697,7 +703,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full, double residual_start, residual_free_start; FILE *fh = NULL; - try_reindex(cr, full); + try_reindex(cr, full, sym, amb); residual_start = residual(cr, full, 0, NULL, NULL, 0); residual_free_start = residual(cr, full, 1, NULL, NULL, 0); @@ -858,24 +864,6 @@ static void do_pr_refine(Crystal *cr, const RefList *full, } -static struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose, int serial, - int cycle, int write_logs) -{ - struct prdata prdata; - - prdata.refined = 0; - - do_pr_refine(cr, full, pmodel, verbose, serial, cycle, write_logs); - - if ( crystal_get_user_flag(cr) == 0 ) { - prdata.refined = 1; - } - - return prdata; -} - - struct refine_args { RefList *full; @@ -886,6 +874,8 @@ struct refine_args int verbose; int cycle; int no_logs; + SymOpList *sym; + SymOpList *amb; }; @@ -906,9 +896,15 @@ static void refine_image(void *task, int id) int write_logs = 0; write_logs = !pargs->no_logs && (pargs->serial % 20 == 0); - pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel, - pargs->verbose, pargs->serial, pargs->cycle, - write_logs); + pargs->prdata.refined = 0; + + do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose, + pargs->serial, pargs->cycle, write_logs, + pargs->sym, pargs->amb); + + if ( crystal_get_user_flag(cr) == 0 ) { + pargs->prdata.refined = 1; + } } @@ -942,7 +938,8 @@ static void done_image(void *vqargs, void *task) void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, - int verbose, int cycle, int no_logs) + int verbose, int cycle, int no_logs, + SymOpList *sym, SymOpList *amb) { struct refine_args task_defaults; struct queue_args qargs; @@ -954,6 +951,8 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.verbose = verbose; task_defaults.cycle = cycle; task_defaults.no_logs = no_logs; + task_defaults.sym = sym; + task_defaults.amb = amb; qargs.task_defaults = task_defaults; qargs.n_started = 0; diff --git a/src/post-refinement.h b/src/post-refinement.h index 3c4ac7ef..a27127ff 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -41,6 +41,7 @@ #include "utils.h" #include "crystal.h" #include "geometry.h" +#include "symmetry.h" enum prflag @@ -59,7 +60,8 @@ extern const char *str_prflag(enum prflag flag); extern void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, - int verbose, int cycle, int no_logs); + int verbose, int cycle, int no_logs, + SymOpList *sym, SymOpList *amb); extern void write_gridscan(Crystal *cr, const RefList *full, int cycle, int serial); |