diff options
author | Thomas White <taw@physics.org> | 2024-01-05 16:39:43 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2024-01-05 16:39:43 +0100 |
commit | cc2717022fca0fd6b24b1dba3aaa5173672a3223 (patch) | |
tree | d328fe9f0474e9dacbffffe7cdbe04df578b919f | |
parent | 73992281c975a1a631d2efe7b9b04fccb5459751 (diff) |
get_hkl: Read MTZ files
There are still some rough edges, e.g. it only works with a simple
I/SIGI column (not I+/I-), and can't yet interpret the symmetry
information in the file. However, it's still better than the old
mtz2hkl script.
Closes: https://gitlab.desy.de/thomas.white/crystfel/-/issues/7
-rw-r--r-- | libcrystfel/src/reflist-utils.c | 144 | ||||
-rw-r--r-- | libcrystfel/src/reflist-utils.h | 1 | ||||
-rwxr-xr-x | scripts/mtz2hkl | 44 |
3 files changed, 137 insertions, 52 deletions
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c index f2229541..66a4996c 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -159,6 +159,85 @@ int find_equiv_in_list(RefList *list, signed int h, signed int k, } +RefList *read_mtz(const char *filename, char **psym_name, UnitCell **pcell) +{ +#ifdef HAVE_LIBCCP4 + int nspg; + MTZ *mtz; + int done; + int i; + const MTZCOL *columns[5]; + MTZXTAL *xtal; + RefList *refls; + CCP4SPG *spg; + UnitCell *cell; + + mtz = MtzGet(filename, 0); + if ( mtz == NULL ) return NULL; + + columns[0] = MtzColLookup(mtz, "H"); + columns[1] = MtzColLookup(mtz, "K"); + columns[2] = MtzColLookup(mtz, "L"); + columns[3] = MtzColLookup(mtz, "I"); + columns[4] = MtzColLookup(mtz, "SIGI"); + + if ( columns[3] == NULL ) { + /* FIXME: Try I+/I- */ + STATUS("Couldn't find intensity column in MTZ file\n"); + return NULL; + } + + xtal = MtzIxtal(mtz, 0); + cell = cell_new_from_parameters(xtal->cell[0]*1e-10, + xtal->cell[1]*1e-10, + xtal->cell[2]*1e-10, + deg2rad(xtal->cell[3]), + deg2rad(xtal->cell[4]), + deg2rad(xtal->cell[5])); + + nspg = MtzSpacegroupNumber(mtz); + spg = ccp4spg_load_by_ccp4_num(nspg); + /* FIXME: Convert to CrystFEL (oriented) point group name */ + ccp4spg_free(&spg); + + i = 1; + refls = reflist_new(); + do { + + float res; + float vals[5]; + int flags[5]; + Reflection *refl; + signed int h, k, l; + + done = ccp4_lrreff(mtz, &res, vals, flags, columns, 5, i++); + if ( done ) continue; + + h = vals[0]; + k = vals[1]; + l = vals[2]; + + refl = add_refl(refls, h, k, l); + set_intensity(refl, vals[3]); + set_esd_intensity(refl, vals[4]); + set_redundancy(refl, 1); + + } while ( !done ); + + MtzFree(mtz); + + if ( pcell != NULL ) { + *pcell = cell; + } else { + cell_free(cell); + } + return refls; +#else + return NULL; +#endif +} + + /* * Write the actual reflections to the file, no headers etc. * Reflections which have a redundancy of zero will not be written. @@ -422,36 +501,71 @@ static RefList *read_reflections_from_file(FILE *fh, char **sym) /** - * read_reflections_2: + * read_reflections_3: * \param filename: Filename to read from * \param sym: Pointer to a "char *" at which to store the symmetry + * \param cell: Pointer to a "UnitCell *" at which to store the unit cell * * This function reads a reflection list from a file, including the - * symmetry from the header (e.g. "Symmetry: 4/mmm"). + * point group name from the header (e.g. "4/mmm"). + * + * The file can be a CrystFEL reflection data file, or an MTZ file provided + * that CrystFEL has been compiled with libCCP4 available. The unit cell will + * only be returned when reading from an MTZ file - CrystFEL reflection data + * files don't contain this information. * * Returns: A %RefList read from the file, or NULL on error */ -RefList *read_reflections_2(const char *filename, char **sym) +RefList *read_reflections_3(const char *filename, char **sym, UnitCell **cell) { - FILE *fh; - RefList *out; + const char *ext; if ( filename == NULL ) { - fh = stdout; - } else { - fh = fopen(filename, "r"); + return read_reflections_from_file(stdin, sym); } - if ( fh == NULL ) { - ERROR("Couldn't open input file '%s'.\n", filename); - return NULL; - } + ext = filename_extension(filename, NULL); + if ( strcmp(ext, ".mtz") == 0 ) { + return read_mtz(filename, sym, cell); + } else { - out = read_reflections_from_file(fh, sym); + RefList *out; + FILE *fh = fopen(filename, "r"); - fclose(fh); + if ( fh == NULL ) { + ERROR("Couldn't open input file '%s'.\n", filename); + return NULL; + } - return out; + out = read_reflections_from_file(fh, sym); + if ( cell != NULL ) *cell = NULL; + + fclose(fh); + + return out; + + } +} + + +/** + * read_reflections_2: + * \param filename: Filename to read from + * \param sym: Pointer to a "char *" at which to store the symmetry + * + * This function reads a reflection list from a file, including the + * point group name from the header (e.g. "4/mmm"). + * + * The file can be a CrystFEL reflection data file, or an MTZ file provided + * that CrystFEL has been compiled with libCCP4 available. MTZ files contain + * the unit cell parameters - if you want this information, use + * %read_reflections_2 instead. + * + * Returns: A %RefList read from the file, or NULL on error + */ +RefList *read_reflections_2(const char *filename, char **sym) +{ + return read_reflections_3(filename, sym, NULL); } diff --git a/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h index e853698f..e91a48d5 100644 --- a/libcrystfel/src/reflist-utils.h +++ b/libcrystfel/src/reflist-utils.h @@ -56,6 +56,7 @@ extern int write_reflist_2(const char *filename, RefList *list, SymOpList *sym); extern RefList *read_reflections(const char *filename); extern RefList *read_reflections_2(const char *filename, char **sym); +extern RefList *read_reflections_3(const char *filename, char **sym, UnitCell **cell); extern int check_list_symmetry(RefList *list, const SymOpList *sym); extern int find_equiv_in_list(RefList *list, signed int h, signed int k, diff --git a/scripts/mtz2hkl b/scripts/mtz2hkl index e1c98251..5a7409c9 100755 --- a/scripts/mtz2hkl +++ b/scripts/mtz2hkl @@ -1,39 +1,9 @@ #!/bin/sh -mtz2various hklin $1 hklout $2.temp <<EOF -LABIN H=H K=K L=L I=IMEAN SIGI=SIGIMEAN -OUTPUT USER '(3I4,2F15.1)' -EOF - -perl < $2.temp > $2 << WIBBLE -use strict; - -my \$line; -open(FILE, "$2.temp"); - -printf("CrystFEL reflection list version 2.0\n"); -printf("Symmetry: 1\n"); -printf(" h k l I phase sigma(I) nmeas\n"); - -while ( \$line = <FILE> ) { - - if ( \$line =~ /^\s*([\d\-]+)\s+([\d\-]+)\s+([\d\-]+)\s+([\d\-\.]+)\s+([\d\-\.]+)/ ) { - - my \$h = \$1; - my \$k = \$2; - my \$l = \$3; - my \$intensity = \$4; - my \$sigi = \$5; - - printf("%4i %4i %4i %10.2f %s %10.2f %7i\n", - \$h, \$k, \$l, \$intensity, " -", \$sigi, 1); - - } else { - printf(STDERR "Couldn't understand line '%s'\n", \$line); - } - -} -close(FILE); -printf("End of reflections\n"); -WIBBLE -exit +echo ------------------------------------------------------------- +echo +echo In this CrystFEL version, MTZ import is now done via get_hkl: +echo get_hkl -i file.mtz -o out.hkl +echo +echo ------------------------------------------------------------- +exit 1 |