From 1a25445aea5e028c6a6f2cefc8946d9fc4c04a0f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 23 Feb 2018 16:16:03 +0100 Subject: partialator: Add -w and --operator options --- src/partialator.c | 61 +++++++++++++++++++++++++++-- src/post-refinement.c | 105 +++++++++++++++++++++++++------------------------- src/post-refinement.h | 4 +- 3 files changed, 112 insertions(+), 58 deletions(-) (limited to 'src') 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 Apparent point group for resolving ambiguities.\n" +" --operator= 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; idxno_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); -- cgit v1.2.3