From b8ebe9370c6d7df3e2c71665881608bf519276eb Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 15 Mar 2021 16:53:09 +0100 Subject: Move MTZ and XDS export to libcrystfel --- src/gui_export.c | 460 ++----------------------------------------------------- 1 file changed, 15 insertions(+), 445 deletions(-) (limited to 'src/gui_export.c') diff --git a/src/gui_export.c b/src/gui_export.c index c3499827..50bd126c 100644 --- a/src/gui_export.c +++ b/src/gui_export.c @@ -36,12 +36,6 @@ #include #include -#ifdef HAVE_LIBCCP4 -#include -#include -#include -#endif - #include #include #include @@ -65,242 +59,14 @@ struct export_window }; -struct point_group_conversion -{ - char centering; - const char *crystfel; - int friedel; - - int xds_spgnum; - - const char *ccp4; -}; - - -/* Table for converting CrystFEL's knowledge about centering, point group - * and unique axis into something that can be recognised by external programs. - * If xds_spgnum=0, ccp4=NULL, or something is missing form this table entirely, - * it means that there is no way to represent the situation to that program - * without re-indexing the dataset. */ -struct point_group_conversion pg_conversions[] = { - - /* Triclinic */ - {'P', "1", 0, 1, "P 1"}, - {'P', "-1", 1, 1, "P 1"}, - - /* Monoclinic */ - {'P', "2_uaa", 0, 0, "P211"}, - {'P', "m_uaa", 0, 0, "Pm11"}, - {'P', "2/m_uaa", 1, 0, "P211"}, - {'P', "2_uab", 0, 3, "P121"}, - {'P', "m_uab", 0, 0, "P1m1"}, - {'P', "2/m_uab", 1, 3, "P121"}, - {'P', "2_uac", 0, 0, "P112"}, - {'P', "m_uac", 0, 0, "P11m"}, - {'P', "2/m_uac", 1, 0, "P112"}, - {'P', "2", 0, 0, "P121"}, /* unique axis c */ - {'P', "m", 0, 0, "P11m"}, /* unique axis c */ - {'P', "2/m", 1, 0, "P121"}, /* unique axis c */ - - {'A', "2_uab", 0, 0, "A121"}, - {'A', "m_uab", 0, 0, "A1m1"}, - {'A', "2/m_uab", 1, 0, "A121"}, - {'A', "2_uac", 0, 0, "A112"}, - {'A', "m_uac", 0, 0, "A11m"}, - {'A', "2/m_uac", 1, 0, "A112"}, - {'A', "2", 0, 0, "A121"}, /* unique axis c */ - {'A', "m", 0, 0, "A11m"}, /* unique axis c */ - {'A', "2/m", 1, 0, "A121"}, /* unique axis c */ - - {'B', "2_uaa", 0, 0, "B211"}, - {'B', "m_uaa", 0, 0, "Bm11"}, - {'B', "2/m_uaa", 1, 0, "B211"}, - {'B', "2_uac", 0, 0, "B112"}, - {'B', "m_uac", 0, 0, "B11m"}, - {'B', "2/m_uac", 1, 0, "B112"}, - {'B', "2", 0, 0, "B112"}, /* unique axis c */ - {'B', "m", 0, 0, "B11m"}, /* unique axis c */ - {'B', "2/m", 1, 0, "B112"}, /* unique axis c */ - - {'C', "2_uaa", 0, 0, "C211"}, - {'C', "m_uaa", 0, 0, "Cm11"}, - {'C', "2/m_uaa", 1, 0, "C211"}, - {'C', "2_uab", 0, 5, "C121"}, - {'C', "m_uab", 0, 0, "C1m1"}, - {'C', "2/m_uab", 1, 5, "C121"}, - - {'I', "2_uaa", 0, 0, "I211"}, - {'I', "m_uaa", 0, 0, "Im11"}, - {'I', "2/m_uaa", 1, 0, "I211"}, - {'I', "2_uab", 0, 0, "I121"}, - {'I', "m_uab", 0, 0, "I1m1"}, - {'I', "2/m_uab", 1, 0, "I121"}, - {'I', "2_uac", 0, 0, "I112"}, - {'I', "m_uac", 0, 0, "I11m"}, - {'I', "2/m_uac", 1, 0, "I112"}, - {'I', "2", 0, 0, "I121"}, /* unique axis c */ - {'I', "m", 0, 0, "I11m"}, /* unique axis c */ - {'I', "2/m", 1, 0, "I121"}, /* unique axis c */ - - /* Orthorhombic */ - {'P', "222", 0, 16, "P222"}, - {'P', "mmm", 1, 16, "P222"}, - {'P', "mm2", 0, 25, "Pmm2"}, - {'A', "222", 0, 0, "A222"}, - {'A', "mmm", 1, 0, "A222"}, - {'A', "mm2", 0, 38, "Amm2"}, - {'B', "222", 0, 0, "B222"}, - {'B', "mmm", 1, 0, "B222"}, - {'B', "mm2", 0, 0, "Bmm2"}, - {'C', "222", 0, 21, "C222"}, - {'C', "mmm", 1, 21, "C222"}, - {'C', "mm2", 0, 35, "Cmm2"}, - {'F', "222", 0, 22, "F222"}, - {'F', "mmm", 1, 22, "F222"}, - {'F', "mm2", 0, 42, "Fmm2"}, - {'I', "222", 0, 23, "I222"}, - {'I', "mmm", 1, 23, "I222"}, - {'I', "mm2", 0, 45, "Imm2"}, - - /* Tetragonal */ - {'P', "4", 0, 75, "P4"}, /* unique axis c */ - {'P', "4/m", 1, 75, "P4"}, /* unique axis c */ - {'P', "422", 0, 89, "P422"}, /* unique axis c */ - {'P', "4/mmm", 1, 89, "P422"}, /* unique axis c */ - {'P', "4mm", 0, 99, "P4mm"}, /* unique axis c */ - {'P', "-4", 0, 81, "P-4"}, /* unique axis c */ - {'P', "-42m", 0, 111, "P-42m"}, /* unique axis c */ - {'P', "-4m2", 0, 115, "P-4m2"}, /* unique axis c */ - {'P', "4_uac", 0, 75, "P4"}, - {'P', "4/m_uac", 1, 75, "P4"}, - {'P', "422_uac", 0, 89, "P422"}, - {'P', "4/mmm_uac", 1, 89, "P422"}, - {'P', "4mm_uac", 0, 99, "P4mm"}, - {'P', "-4_uac", 0, 81, "P-4"}, - {'P', "-42m_uac", 0, 111, "P-42m"}, - {'P', "-4m2_uac", 0, 115, "P-4m2"}, - {'I', "4", 0, 79, "I4"}, /* unique axis c */ - {'I', "4/m", 1, 79, "I4"}, /* unique axis c */ - {'I', "422", 0, 97, "I422"}, /* unique axis c */ - {'I', "4/mmm", 1, 97, "I422"}, /* unique axis c */ - {'I', "4mm", 0, 107, "I4mm"}, /* unique axis c */ - {'I', "-4", 0, 82, "I-4"}, /* unique axis c */ - {'I', "-42m", 0, 121, "I-42m"}, /* unique axis c */ - {'I', "-4m2", 0, 119, "I-4m2"}, /* unique axis c */ - {'I', "4_uac", 0, 79, "I4"}, - {'I', "4/m_uac", 1, 79, "I4"}, - {'I', "422_uac", 0, 97, "I422"}, - {'I', "4/mmm_uac", 1, 97, "I422"}, - {'I', "4mm_uac", 0, 107, "I4mm"}, - {'I', "-4_uac", 0, 82, "I-4"}, - {'I', "-42m_uac", 0, 121, "I-42m"}, - {'I', "-4m2_uac", 0, 119, "I-4m2"}, - - /* Trigonal (rhombohedral) */ - {'R', "3_R", 0, 0, "R3:R"}, - {'R', "-3_R", 1, 0, "R3:R"}, - {'R', "32_R", 0, 0, "R32:R"}, - {'R', "-3m_R", 1, 0, "R32:R"}, - {'R', "3m_R", 0, 0, "R3m:R"}, - - /* Trigonal (rhombohedral on hexagonal axes) */ - {'H', "3_H", 0, 146, "R3:H"}, - {'H', "-3_H", 1, 146, "R3:H"}, - {'H', "32_H", 0, 155, "R3:H"}, - {'H', "-3m_H", 1, 155, "R3:H"}, - {'H', "3m_H", 0, 0, "R3m:H"}, - - /* Trigonal (hexagonal) */ - {'P', "3_H", 0, 143, "P3"}, - {'P', "-3_H", 1, 143, "P3"}, - {'P', "312_H", 0, 149, "P312"}, - {'P', "-31m_H", 1, 149, "P312"}, - {'P', "321_H", 0, 150, "P321"}, - {'P', "-3m1_H", 1, 150, "P321"}, - {'P', "3m1_H", 0, 156, "P3m1"}, - {'P', "31m_H", 0, 157, "P31m"}, - - /* Hexagonal */ - {'P', "6", 0, 168, "P6"}, - {'P', "6/m", 1, 168, "P6"}, - {'P', "622", 0, 177, "P622"}, - {'P', "6/mmm", 1, 177, "P622"}, - {'P', "6mm", 0, 177, "P6mm"}, - {'P', "-6m2", 0, 187, "P-6m2"}, - {'P', "-62m", 0, 189, "P-62m"}, - - /* Cubic */ - {'P', "23", 0, 195, "P23"}, - {'P', "m-3", 1, 195, "P23"}, - {'P', "432", 0, 207, "P432"}, - {'P', "m-3m", 1, 207, "P432"}, - {'P', "-43m", 0, 215, "P -4 3 m"}, - {'I', "23", 0, 197, "I23"}, - {'I', "m-3", 1, 197, "I23"}, - {'I', "432", 0, 211, "I432"}, - {'I', "m-3m", 1, 211, "I432"}, - {'I', "-43m", 0, 217, "I -4 3 m"}, - {'F', "23", 0, 196, "F23"}, - {'F', "m-3", 1, 196, "F23"}, - {'F', "432", 0, 209, "F432"}, - {'F', "m-3m", 1, 209, "F432"}, - {'F', "-43m", 0, 216, "F -4 3 m"}, - - {'*', NULL, 0, 0, NULL} -}; - - -static int space_group_for_xds(const char *sym_str, char cen) -{ - int i = 0; - do { - if ( (pg_conversions[i].centering == cen) - && (strcmp(sym_str, pg_conversions[i].crystfel) == 0) ) - { - return pg_conversions[i].xds_spgnum; - } - i++; - } while (pg_conversions[i].centering != '*'); - - ERROR("Couldn't derive XDS representation of symmetry.\n"); - return 0; -} - - -#ifdef HAVE_LIBCCP4 -static const char *space_group_for_mtz(const char *sym_str, char cen) -{ - int i = 0; - do { - if ( (pg_conversions[i].centering == cen) - && (strcmp(sym_str, pg_conversions[i].crystfel) == 0) ) - { - return pg_conversions[i].ccp4; - } - i++; - } while (pg_conversions[i].centering != '*'); - - ERROR("Couldn't derive CCP4 representation of symmetry.\n"); - return NULL; -} -#endif - - static int export_to_xds(struct gui_merge_result *result, const char *filename, UnitCell *cell, double min_res, double max_res) { - FILE *fh; RefList *reflist; - RefListIterator *iter; - Reflection *refl; - double a, b, c, al, be,ga; char *sym_str; SymOpList *sym; - int spg; - - fh = fopen(filename, "w"); - if ( fh == NULL ) return 1; + int r; reflist = read_reflections_2(result->hkl, &sym_str); if ( reflist == NULL ) return 1; @@ -309,104 +75,24 @@ static int export_to_xds(struct gui_merge_result *result, sym = get_pointgroup(sym_str); if ( sym == NULL ) return 1; - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - - spg = space_group_for_xds(sym_str, cell_get_centering(cell)); - if ( spg == 0 ) return 1; - - fprintf(fh, "!FORMAT=XDS_ASCII MERGE=TRUE FRIEDEL'S_LAW=%s\n", - is_centrosymmetric(sym) ? "TRUE" : "FALSE"); - fprintf(fh, "!SPACE_GROUP_NUMBER=%i\n", spg); - fprintf(fh, "!UNIT_CELL_CONSTANT= %.2f %.2f %.2f %.2f %.2f %.2f\n", - a*1e10, b*1e10, c*1e10, rad2deg(al), rad2deg(be), rad2deg(ga)); - fprintf(fh, "!NUMBER_OF_ITEMS_IN_EACH_DATA_RECORD=5\n"); - fprintf(fh, "!ITEM_H=1\n"); - fprintf(fh, "!ITEM_K=2\n"); - fprintf(fh, "!ITEM_L=3\n"); - fprintf(fh, "!ITEM_IOBS=4\n"); - fprintf(fh, "!ITEM_SIGMA(IOBS)=5\n"); - fprintf(fh, "!END_OF_HEADER\n"); - - for ( refl = first_refl(reflist, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - double one_over_d; - - get_indices(refl, &h, &k, &l); - - one_over_d = 2.0*resolution(cell, h, k, l); - if ( (one_over_d > min_res) && (one_over_d < max_res) ) { - - fprintf(fh, "%6i %6i %6i %9.2f %9.2f\n", - h, k, l, - get_intensity(refl), - get_esd_intensity(refl)); - - } - } + r = write_to_xds(reflist, sym, cell, min_res, max_res, filename); - fprintf(fh, "!END_OF_DATA\n"); free_symoplist(sym); free(sym_str); reflist_free(reflist); - fclose(fh); - return 0; -} - - -#ifdef HAVE_LIBCCP4 -static CCP4SPG *add_mtz_symmetry_header(MTZ *mtz, const char *spg_name) -{ - CCP4SPG *spg; - float rsymx[192][4][4]; - char ltypex[2]; - int i; - - spg = ccp4spg_load_by_spgname(spg_name); - if ( spg == NULL ) { - ERROR("Couldn't look up CCP4 space group '%s'\n", spg_name); - return NULL; - } - - for ( i=0; insymop; i++ ) { - rotandtrn_to_mat4(rsymx[i], spg->symop[i]); - } - ltypex[0] = spg->symbol_old[0]; - ltypex[1] = '\0'; - - ccp4_lwsymm(mtz, spg->nsymop, spg->nsymop_prim, - rsymx, ltypex, spg->spg_ccp4_num, spg->symbol_old, - spg->point_group); - - return spg; + return r; } -#endif static int export_to_mtz(struct gui_merge_result *result, const char *filename, UnitCell *cell, double min_res, double max_res) { -#ifdef HAVE_LIBCCP4 - MTZ *mtz; - MTZXTAL *cr; - MTZSET *ds; - MTZCOL *columns[7]; - double a, b, c, al, be, ga; - int r; - char tmp[128]; - float cellp[6]; - int refl_i; RefList *reflist; - Reflection *refl; - RefListIterator *iter; - char *sym_str = NULL; - CCP4SPG *spg; - const char *spg_name; + char *sym_str; SymOpList *sym; + int r; reflist = read_reflections_2(result->hkl, &sym_str); if ( reflist == NULL ) return 1; @@ -415,130 +101,14 @@ static int export_to_mtz(struct gui_merge_result *result, sym = get_pointgroup(sym_str); if ( sym == NULL ) return 1; - spg_name = space_group_for_mtz(sym_str, cell_get_centering(cell)); - if ( spg_name == NULL ) { - reflist_free(reflist); - return 1; - } - - mtz = MtzMalloc(0, 0); - - snprintf(tmp, 128, "Data exported via CrystFEL GUI, version %s", - crystfel_version_string()); - ccp4_lwtitl(mtz, tmp, 0); - - mtz->refs_in_memory = 0; - mtz->fileout = MtzOpenForWrite(filename); - - spg = add_mtz_symmetry_header(mtz, spg_name); - if ( spg == NULL ) { - return 1; - } - - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - cellp[0] = a*1e10; - cellp[1] = b*1e10; - cellp[2] = c*1e10; - cellp[3] = rad2deg(al); - cellp[4] = rad2deg(be); - cellp[5] = rad2deg(ga); - - /* FIXME: Proposed labelling: - * title = as above - * project = basename of folder containing crystfel.project - * crystal = name of indexing results run - * dataset = name of merge results run */ - cr = MtzAddXtal(mtz, "Crystal_name", "Project_name", cellp); - ds = MtzAddDataset(mtz, cr, result->name, 0.0); - columns[0] = MtzAddColumn(mtz, ds, "H", "H"); - columns[1] = MtzAddColumn(mtz, ds, "K", "H"); - columns[2] = MtzAddColumn(mtz, ds, "L", "H"); - columns[3] = MtzAddColumn(mtz, ds, "I+", "J"); - columns[4] = MtzAddColumn(mtz, ds, "SIGI+", "Q"); - columns[5] = MtzAddColumn(mtz, ds, "I-", "J"); - columns[6] = MtzAddColumn(mtz, ds, "SIGI-", "Q"); - - refl_i = 1; - for ( refl = first_refl(reflist, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - double one_over_d; - int isym; - - get_indices(refl, &h, &k, &l); - - one_over_d = 2.0*resolution(cell, h, k, l); - if ( (one_over_d > min_res) && (one_over_d < max_res) ) { - - float refldata[7]; - signed int nh, nk, nl; - signed int fh, fk, fl; - Reflection *friedel; - Reflection *refl_plus; - Reflection *refl_minus; - - /* Look for Friedel partner */ - if ( find_equiv_in_list(reflist, -h, -k, -l, - sym, &fh, &fk, &fl) ) - { - friedel = find_refl(reflist, fh, fk, fl); - } else { - friedel = NULL; - } - - /* Move to CCP4's idea of the ASU */ - isym = ccp4spg_put_in_asu(spg, h, k, l, &nh, &nk, &nl); - - /* Ok, do we have an I+ or an I- ? */ - if ( is_odd(isym) ) { - /* I+ */ - refl_plus = refl; - refl_minus = friedel; - } else { - /* I- */ - refl_minus = refl; - refl_plus = friedel; - } - - /* If we are looking at an I-, only write it out now - * if the corresponding I+ if not in 'reflist'. - * If I+ is present, then this I- will get written when - * the Friedel pair is processed. */ - if ( !is_odd(isym) && (refl_plus != NULL) ) continue; - - refldata[0] = nh; - refldata[1] = nk; - refldata[2] = nl; - if ( refl_plus != NULL ) { - refldata[3] = get_intensity(refl_plus); - refldata[4] = get_esd_intensity(refl_plus); - } else { - refldata[3] = NAN; - refldata[4] = NAN; - } - if ( refl_minus != NULL ) { - refldata[5] = get_intensity(refl_minus); - refldata[6] = get_esd_intensity(refl_minus); - } else { - refldata[5] = NAN; - refldata[6] = NAN; - } - - ccp4_lwrefl(mtz, refldata, columns, 7, refl_i++); - - } - } + r = write_to_mtz(reflist, sym, cell, min_res, max_res, filename, + result->name); - r = MtzPut(mtz, " "); - ccp4spg_free(&spg); - MtzFree(mtz); + free_symoplist(sym); + free(sym_str); reflist_free(reflist); - return 1-r; /* Yes, really. MtzPut return values are backwards */ -#else - return 1; -#endif + + return r; } @@ -680,10 +250,10 @@ gint export_sig(GtkWidget *widget, struct crystfelproject *proj) win->format = gtk_combo_box_text_new(); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(win->format), FALSE, FALSE, 4.0); -#ifdef HAVE_LIBCCP4 - gtk_combo_box_text_append(GTK_COMBO_BOX_TEXT(win->format), "mtz", - "MTZ"); -#endif + if ( libcrystfel_can_write_mtz() ) { + gtk_combo_box_text_append(GTK_COMBO_BOX_TEXT(win->format), "mtz", + "MTZ"); + } gtk_combo_box_text_append(GTK_COMBO_BOX_TEXT(win->format), "xds", "XDS ASCII"); gtk_combo_box_set_active(GTK_COMBO_BOX(win->format), 0); -- cgit v1.2.3