Move MTZ and XDS export to libcrystfel

6 files changed, 504 insertions, 1 deletions

diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index baada975..f1bd3c14 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -9,6 +9,7 @@ find_package(FLEX REQUIRED)
 find_package(BISON REQUIRED)
 find_package(Doxygen)
 pkg_search_module(FFTW fftw3)
+pkg_search_module(LIBCCP4 libccp4c)
 
 configure_file(src/libcrystfel-version.c.cmake.in libcrystfel-version.c)
 
@@ -17,6 +18,7 @@ set(HAVE_FFTW ${FFTW_FOUND})
 set(HAVE_XGANDALF ${XGANDALF_FOUND})
 set(HAVE_PINKINDEXER ${PINKINDEXER_FOUND})
 set(HAVE_FDIP ${FDIP_FOUND})
+set(HAVE_LIBCCP4 ${LIBCCP4_FOUND})
 
 # Recent enough version of zlib?
 set(CMAKE_REQUIRED_LIBRARIES "-lz")
@@ -161,6 +163,11 @@ if (CURSES_FOUND)
   target_link_libraries(${PROJECT_NAME} PRIVATE ${CURSES_LIBRARIES})
 endif (CURSES_FOUND)
 
+if (LIBCCP4_FOUND)
+	target_include_directories(${PROJECT_NAME} PRIVATE ${LIBCCP4_INCLUDES})
+	target_link_libraries(${PROJECT_NAME} PRIVATE ${LIBCCP4_LIBRARIES})
+endif (LIBCCP4_FOUND)
+
 target_compile_options(${PROJECT_NAME} PRIVATE -Wall)
 set_target_properties(${PROJECT_NAME} PROPERTIES PUBLIC_HEADER "${LIBCRYSTFEL_HEADERS}")
 
diff --git a/libcrystfel/config.h.cmake.in b/libcrystfel/config.h.cmake.in
index 7f0dcd94..63211779 100644
--- a/libcrystfel/config.h.cmake.in
+++ b/libcrystfel/config.h.cmake.in
@@ -8,6 +8,7 @@
 #cmakedefine HAVE_CURSES
 #cmakedefine HAVE_GZBUFFER
 #cmakedefine HAVE_GDKPIXBUF
+#cmakedefine HAVE_LIBCCP4
 
 #cmakedefine HAVE_FORKPTY_PTY_H
 #cmakedefine HAVE_FORKPTY_UTIL_H
diff --git a/libcrystfel/config.h.in b/libcrystfel/config.h.in
index 1fdd2e13..c00f559e 100644
--- a/libcrystfel/config.h.in
+++ b/libcrystfel/config.h.in
@@ -6,6 +6,7 @@
 #mesondefine HAVE_FDIP
 #mesondefine HAVE_CURSES
 #mesondefine HAVE_GZBUFFER
+#mesondefine HAVE_LIBCCP4
 
 #mesondefine HAVE_FORKPTY_PTY_H
 #mesondefine HAVE_FORKPTY_UTIL_H
diff --git a/libcrystfel/meson.build b/libcrystfel/meson.build
index 8b8eb77d..38b8002d 100644
--- a/libcrystfel/meson.build
+++ b/libcrystfel/meson.build
@@ -32,6 +32,11 @@ if fdipdep.found()
   conf_data.set10('HAVE_FDIP', 1)
 endif
 
+ccp4dep = dependency('libccp4c', required: false)
+if ccp4dep.found()
+  conf_data.set10('HAVE_LIBCCP4', 1)
+endif
+
 libcrystfel_versionc = vcs_tag(input: 'src/libcrystfel-version.c.in',
                                output: 'libcrystfel-version.c')
 
@@ -122,7 +127,8 @@ libcrystfel = library('crystfel', [libcrystfel_sources, libcrystfel_versionc],
                       include_directories : [libcrystfel_includes, libcrystfel_conf_inc],
                       dependencies : [mdep, utildep, fftwdep, gsldep, zlibdep,
                                       hdf5dep, pthreaddep, ncursesdep,
-                                      xgandalfdep, pinkindexerdep, fdipdep],
+                                      xgandalfdep, pinkindexerdep, fdipdep,
+                                      ccp4dep],
                       install : true)
 
 libcrystfeldep = declare_dependency(include_directories : libcrystfel_includes,
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 42c3f06d..90f7b91b 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -35,6 +35,12 @@
 #include <stdio.h>
 #include <assert.h>
 
+#ifdef HAVE_LIBCCP4
+#include <ccp4/cmtzlib.h>
+#include <ccp4/csymlib.h>
+#include <ccp4/ccp4_parser.h>
+#endif
+
 #include "reflist.h"
 #include "cell.h"
 #include "cell-utils.h"
@@ -707,3 +713,475 @@ void reflist_add_command_and_version(RefList *list, int argc, char *argv[])
 	reflist_add_notes(list, tmp);
 	free(tmp);
 }
+
+
+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
+
+
+int write_to_xds(RefList *reflist,
+                 SymOpList *sym,
+                 UnitCell *cell,
+                 double min_res,
+                 double max_res,
+                 const char *filename)
+{
+	FILE *fh;
+	RefListIterator *iter;
+	Reflection *refl;
+	double a, b, c, al, be,ga;
+	int spg;
+
+	cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+
+	spg = space_group_for_xds(symmetry_name(sym),
+	                          cell_get_centering(cell));
+	if ( spg == 0 ) return 1;
+
+	fh = fopen(filename, "w");
+	if ( fh == NULL ) 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));
+
+		}
+	}
+
+	fprintf(fh, "!END_OF_DATA\n");
+
+	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; i<spg->nsymop; 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;
+}
+#endif
+
+
+int write_to_mtz(RefList *reflist,
+                 SymOpList *sym,
+                 UnitCell *cell,
+                 double min_res,
+                 double max_res,
+                 const char *filename,
+                 const char *dataset_name)
+{
+#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;
+	Reflection *refl;
+	RefListIterator *iter;
+	CCP4SPG *spg;
+	const char *spg_name;
+
+	spg_name = space_group_for_mtz(symmetry_name(sym),
+	                               cell_get_centering(cell));
+	if ( spg_name == NULL ) {
+		reflist_free(reflist);
+		return 1;
+	}
+
+	mtz = MtzMalloc(0, 0);
+
+	snprintf(tmp, 128, "Data exported via CrystFEL version %s",
+	         libcrystfel_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, dataset_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 = MtzPut(mtz, " ");
+	ccp4spg_free(&spg);
+	MtzFree(mtz);
+	return 1-r; /* Yes, really.  MtzPut return values are backwards */
+#else
+	return 1;
+#endif
+}
+
+
+int libcrystfel_can_write_mtz()
+{
+#ifdef HAVE_LIBCCP4
+	return 1;
+#else
+	return 0;
+#endif
+}
diff --git a/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h
index f5177cee..bd60c642 100644
--- a/libcrystfel/src/reflist-utils.h
+++ b/libcrystfel/src/reflist-utils.h
@@ -79,6 +79,16 @@ extern void free_contribs(RefList *list);
 extern void reflist_add_command_and_version(RefList *list,
                                             int argcv, char *argv[]);
 
+extern int write_to_mtz(RefList *reflist, SymOpList *sym, UnitCell *cell,
+                        double min_res, double max_res,
+                        const char *filename, const char *dataset_name);
+
+extern int write_to_xds(RefList *reflist, SymOpList *sym, UnitCell *cell,
+                        double min_res, double max_res, const char *filename);
+
+
+extern int libcrystfel_can_write_mtz(void);
+
 #ifdef __cplusplus
 }
 #endif