Commit 74cfb2b2 authored by Thomas White's avatar Thomas White
Browse files

process_hkl: Take input ("apparent") symmetry on command line

parent 30012bfe
......@@ -66,6 +66,7 @@ static void show_help(const char *s)
" --scale Scale each pattern for best fit with the current\n"
" model.\n"
" -y, --symmetry=<sym> Merge according to point group <sym>.\n"
" -a, --input-symmetry=<a> Specify the apparent (input) symmetry.\n"
" --reference=<file> Compare against intensities from <file> when\n"
" scaling or resolving ambiguities.\n"
" The symmetry of the reference list must be the\n"
......@@ -552,11 +553,11 @@ int main(int argc, char *argv[])
int config_rmerge = 0;
unsigned int n_total_patterns;
char *sym = NULL;
char *in_sym = NULL;
char *pdb = NULL;
ReflItemList *twins;
ReflItemList *observed;
int i;
const char *holo = NULL;
char *histo = NULL;
signed int hist_h, hist_k, hist_l;
double *hist_vals = NULL;
......@@ -580,10 +581,11 @@ int main(int argc, char *argv[])
{"sum", 0, &config_sum, 1},
{"scale", 0, &config_scale, 1},
{"symmetry", 1, NULL, 'y'},
{"input-symmetry", 1, NULL, 'a'},
{"pdb", 1, NULL, 'p'},
{"histogram", 1, NULL, 'g'},
{"rmerge", 0, &config_rmerge, 1},
{"outstream", 1, NULL, 'a'},
{"outstream", 1, NULL, 2},
{"reference", 1, NULL, 'r'},
{"beam", 1, NULL, 'b'},
{0, 0, NULL, 0}
......@@ -630,10 +632,14 @@ int main(int argc, char *argv[])
reference = strdup(optarg);
break;
case 'a' :
case 2 :
outstream = strdup(optarg);
break;
case 'a' :
in_sym = strdup(optarg);
break;
case 'b' :
beam = get_beam_parameters(optarg);
if ( beam == NULL ) {
......@@ -670,33 +676,35 @@ int main(int argc, char *argv[])
cell = load_cell_from_pdb(pdb);
free(pdb);
if ( sym == NULL ) sym = strdup("1");
/* Show useful symmetry information */
if ( sym != NULL ) {
holo = get_holohedral(sym);
int np = num_general_equivs(holo) / num_general_equivs(sym);
if ( (sym != NULL) && (in_sym != NULL) ) {
int np = num_general_equivs(in_sym) / num_general_equivs(sym);
if ( np > 1 ) {
STATUS("Resolving point group %s into %s "
"(%i possibilities)\n",
holo, sym, np);
in_sym, sym, np);
/* Get the list of twin/Bijvoet possibilities */
twins = get_twin_possibilities(holo, sym);
STATUS("Twin/inversion operation indices from %s are:",
holo);
twins = get_twin_possibilities(in_sym, sym);
STATUS("Twin operation indices from %s are:", in_sym);
for ( i=0; i<num_items(twins); i++ ) {
STATUS(" %i", get_item(twins, i)->op);
}
STATUS("\n");
} else {
STATUS("No twin/inversion resolution necessary.\n");
STATUS("No resolution required to get from %s to %s\n",
in_sym, sym);
twins = NULL;
}
} else {
STATUS("Not performing any twin/inversion resolution.\n");
STATUS("No twin resolution requested (use -a otherwise).\n");
twins = NULL;
sym = strdup("1");
holo = strdup("1");
in_sym = strdup(sym);
}
if ( histo != NULL ) {
......@@ -758,7 +766,7 @@ int main(int argc, char *argv[])
merge_all(fh, &model, &observed, &counts,
config_maxonly, config_scale, config_sum,
config_startafter, config_stopafter,
twins, holo, sym, n_total_patterns,
twins, in_sym, sym, n_total_patterns,
reference_items, reference_i,
hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL,
outfh);
......@@ -801,7 +809,8 @@ int main(int argc, char *argv[])
rewind(fh);
merge_all(fh, &model, &observed, &counts,
config_maxonly, config_scale, 0,
config_startafter, config_stopafter, twins, holo, sym,
config_startafter, config_stopafter,
twins, in_sym, sym,
n_total_patterns, reference_items, reference_i,
NULL, 0, 0, 0, NULL, devs, tots, model, NULL);
......
......@@ -328,31 +328,6 @@ void get_asymm(signed int h, signed int k, signed int l,
}
const char *get_holohedral(const char *sym)
{
/* Triclinic */
if ( strcmp(sym, "1") == 0 ) return "-1";
if ( strcmp(sym, "-1") == 0 ) return "-1";
if ( strcmp(sym, "222") == 0 ) return "mmm";
if ( strcmp(sym, "mmm") == 0 ) return "mmm";
/* Tetragonal */
if ( strcmp(sym, "422") == 0 ) return "4/mmm";
if ( strcmp(sym, "4/mmm") == 0 ) return "4/mmm";
/* Hexagonal */
if ( strcmp(sym, "6") == 0 ) return "6/mmm";
if ( strcmp(sym, "6/m") == 0 ) return "6/mmm";
if ( strcmp(sym, "6/mmm") == 0 ) return "6/mmm";
/* TODO: Add more groups here */
ERROR("Couldn't find holohedral point group for '%s'\n", sym);
abort();
}
/* This is kind of like a "numerical" left coset decomposition.
* Given a reflection index and a point group, it returns the list of twinning
* possibilities.
......
......@@ -35,8 +35,6 @@ extern void get_general_equiv(signed int h, signed int k, signed int l,
signed int *he, signed int *ke, signed int *le,
const char *sym, int idx);
extern const char *get_holohedral(const char *sym);
extern ReflItemList *get_twins(ReflItemList *items,
const char *holo, const char *mero);
......
......@@ -113,7 +113,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
double nominal_photon_energy)
{
struct _indexingprivate_template *priv;
const char *holo;
const char *sym;
double omega_max, phi_max;
int n_templates;
const double omega_step = deg2rad(0.5);
......@@ -137,24 +137,23 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
priv = calloc(1, sizeof(struct _indexingprivate_template));
priv->base.indm = INDEXING_TEMPLATE;
/* We can only distinguish orientations within the holohedral cell */
holo = get_holohedral(cell_get_pointgroup(cell));
sym = cell_get_pointgroup(cell);
/* These define the orientation in space */
if ( is_polyhedral(holo) ) {
ERROR("WARNING: Holohedral point group is polyhedral.\n");
if ( is_polyhedral(sym) ) {
ERROR("WARNING: Point group is polyhedral.\n");
ERROR("This means I can't properly determine the orientation");
ERROR(" ranges for template matching. Expect trouble.\n");
}
omega_max = 2.0*M_PI / rotational_order(holo);
if ( has_bisecting_mirror_or_diad(holo) ) omega_max /= 2.0;
omega_max = 2.0*M_PI / rotational_order(sym);
if ( has_bisecting_mirror_or_diad(sym) ) omega_max /= 2.0;
phi_max = M_PI;
if ( has_perpendicular_mirror(holo) ) phi_max /= 2.0;
if ( has_perpendicular_mirror(sym) ) phi_max /= 2.0;
/* One more axis would define the rotation in the plane of the image */
STATUS("Orientation ranges in %s: %.0f-%.0f, %.0f-%.0f deg.\n",
holo, 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max));
sym, 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max));
n_templates = (omega_max * phi_max)/(omega_step * phi_step);
STATUS("%i templates to be calculated.\n", n_templates);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment