Remove all remaining PDB rendering stuff

13 files changed, 30 insertions, 691 deletions

diff --git a/src/check_hkl.c b/src/check_hkl.c
index b36150b2..fb6d02d5 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -22,7 +22,6 @@
 #include <getopt.h>
 
 #include "utils.h"
-#include "sfac.h"
 #include "reflections.h"
 #include "statistics.h"
 #include "symmetry.h"
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index dca4a298..abe8bc6a 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -22,7 +22,6 @@
 #include <getopt.h>
 
 #include "utils.h"
-#include "sfac.h"
 #include "reflections.h"
 #include "statistics.h"
 #include "symmetry.h"
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 74eca439..c0814e79 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -29,7 +29,6 @@
 #include "utils.h"
 #include "cell.h"
 #include "diffraction.h"
-#include "sfac.h"
 #include "cl-utils.h"
 #include "beam-parameters.h"
 
diff --git a/src/diffraction.c b/src/diffraction.c
index e3d7435c..cd372a28 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -22,7 +22,6 @@
 #include "utils.h"
 #include "cell.h"
 #include "diffraction.h"
-#include "sfac.h"
 #include "beam-parameters.h"
 #include "symmetry.h"
 
diff --git a/src/dirax.c b/src/dirax.c
index 35a2b1ce..291dc85b 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -36,7 +36,6 @@
 #include "image.h"
 #include "dirax.h"
 #include "utils.h"
-#include "sfac.h"
 #include "peaks.h"
 
 
diff --git a/src/get_hkl.c b/src/get_hkl.c
index cddb54b8..f8a43ddb 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -1,9 +1,9 @@
 /*
  * get_hkl.c
  *
- * Small program to write out a list of h,k,l,I values given a structure
+ * Small program to manipulate reflection lists
  *
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
  *
  * Part of CrystFEL - crystallography with a FEL
  *
@@ -22,7 +22,6 @@
 #include <getopt.h>
 
 #include "utils.h"
-#include "sfac.h"
 #include "reflections.h"
 #include "symmetry.h"
 #include "beam-parameters.h"
@@ -32,27 +31,34 @@ static void show_help(const char *s)
 {
 	printf("Syntax: %s [options]\n\n", s);
 	printf(
-"Create reflections lists.\n"
+"Manipulate reflection lists.\n"
 "\n"
 "  -h, --help                 Display this help message.\n"
 "\n"
-"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
+"  -i, --input=<file>         Read reflections from <file>.\n"
+"  -y, --symmetry=<sym>       The symmetry of the input reflection list.\n"
+"\n"
+"You can add noise to the reflections with either of:\n"
 "      --poisson              Simulate Poisson samples.\n"
 "      --noise                Add 10%% random noise.\n"
-"  -y, --symmetry=<sym>       The symmetry of the input file (-i).\n"
+"\n"
+"To calculate Poisson samples accurately, you must also give:\n"
+"  -b, --beam=<file>          Get beam parameters from file.\n"
+"\n"
+"You can artificially 'twin' the reflections, or expand them out:\n"
 "  -w, --twin=<sym>           Generate twinned data according to the given\n"
 "                              point group.\n"
 "  -e, --expand=<sym>         Expand reflections to this point group.\n"
-"  -o, --output=<filename>    Output filename (default: stdout).\n"
-"  -i, --intensities=<file>   Read intensities from file instead of\n"
-"                              calculating them from scratch.  You might use\n"
-"                              this if you need to apply noise or twinning.\n"
-"  -p, --pdb=<file>           PDB file from which to get the structure.\n"
-"      --no-phases            Do not try to use phases in the input file.\n"
+"\n"
+"You can restrict which reflections are written out:\n"
+"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
+"\n"
+"You might sometimes need to do this:\n"
 "      --multiplicity         Multiply intensities by the number of\n"
 "                              equivalent reflections.\n"
-"  -b, --beam=<file>          Get beam parameters from file (used for sigmas).\n"
-"      --max-res=<d>          Calculate structure factors out to d=<d> nm.\n"
+"\n"
+"Don't forget to specify the output filename:\n"
+"  -o, --output=<filename>    Output filename (default: stdout).\n"
 );
 }
 
@@ -253,7 +259,6 @@ int main(int argc, char *argv[])
 	double *ideal_ref;
 	double *phases;
 	double *esds;
-	struct molecule *mol;
 	char *template = NULL;
 	int config_noise = 0;
 	int config_poisson = 0;
@@ -269,10 +274,7 @@ int main(int argc, char *argv[])
 	ReflItemList *write_items;
 	UnitCell *cell = NULL;
 	char *beamfile = NULL;
-	char *rval;
-	struct beam_params *beam;  /* Beam parameters for SF calculation */
-	int have_max_res = 0;
-	double max_res = 0.0;
+	struct beam_params *beam;  /* Beam parameters for Poisson calculation */
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -285,16 +287,13 @@ int main(int argc, char *argv[])
 		{"twin",               1, NULL,               'w'},
 		{"expand",             1, NULL,               'e'},
 		{"intensities",        1, NULL,               'i'},
-		{"pdb",                1, NULL,               'p'},
-		{"no-phases",          0, &config_nophase,     1},
 		{"multiplicity",       0, &config_multi,       1},
 		{"beam",               1, NULL,               'b'},
-		{"max-res",            1, NULL,                2},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "ht:o:i:p:w:y:e:b:",
+	while ((c = getopt_long(argc, argv, "ht:o:i:w:y:e:b:",
 	                        longopts, NULL)) != -1) {
 
 		switch (c) {
@@ -314,10 +313,6 @@ int main(int argc, char *argv[])
 			input = strdup(optarg);
 			break;
 
-		case 'p' :
-			filename = strdup(optarg);
-			break;
-
 		case 'y' :
 			mero = strdup(optarg);
 			break;
@@ -334,16 +329,6 @@ int main(int argc, char *argv[])
 			beamfile = strdup(optarg);
 			break;
 
-		case 2 :
-			max_res = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid maximum resolution.\n");
-				return 1;
-			}
-			max_res = 1.0 / (max_res * 1.0e-9);
-			have_max_res = 1;
-			break;
-
 		case 0 :
 			break;
 
@@ -363,7 +348,6 @@ int main(int argc, char *argv[])
 		exit(1);
 	}
 
-	mol = load_molecule(filename);
 	cell = load_cell_from_pdb(filename);
 	if ( !config_nophase ) {
 		phases = new_list_phase();
@@ -371,44 +355,15 @@ int main(int argc, char *argv[])
 		phases = NULL;
 	}
 	esds = new_list_sigma();
-	if ( input == NULL ) {
-
-		if ( beamfile == NULL ) {
-			ERROR("To calculate structure factors, you must"
-			      " provide a beam parameters file (use -b)\n");
-			return 1;
-		}
-
-		beam = get_beam_parameters(beamfile);
-		if ( beam == NULL ) {
-			ERROR("Failed to read beam parameters from '%s'\n", beamfile);
-			return 1;
-		}
-		free(beamfile);
-
-		if ( !have_max_res ) {
-			STATUS("You didn't specify the maximum resolution to"
-			       " calculate structure factors.  I'll go to"
-			       " d = 0.5 nm.\n");
-			max_res = 1.0/0.5e-9;
-		}
-
-		input_items = new_items();
-		ideal_ref = get_reflections(mol, eV_to_J(beam->photon_energy),
-		                            max_res, phases, input_items);
-
-	} else {
-
-		ideal_ref = new_list_intensity();
-		input_items = read_reflections(input, ideal_ref, phases,
-		                               NULL, esds);
-		free(input);
-		if ( check_symmetry(input_items, mero) ) {
-			ERROR("The input reflection list does not appear to"
-			      " have symmetry %s\n", mero);
-			return 1;
-		}
 
+	ideal_ref = new_list_intensity();
+	input_items = read_reflections(input, ideal_ref, phases,
+	                               NULL, esds);
+	free(input);
+	if ( check_symmetry(input_items, mero) ) {
+		ERROR("The input reflection list does not appear to"
+		      " have symmetry %s\n", mero);
+		return 1;
 	}
 
 	if ( config_poisson ) poisson_reflections(ideal_ref, input_items);
diff --git a/src/index.c b/src/index.c
index d24b5852..317c9200 100644
--- a/src/index.c
+++ b/src/index.c
@@ -26,7 +26,6 @@
 #include "peaks.h"
 #include "dirax.h"
 #include "mosflm.h"
-#include "sfac.h"
 #include "detector.h"
 #include "index.h"
 #include "index-priv.h"
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 80dd7e50..f53e4d85 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -30,7 +30,6 @@
 #include "index.h"
 #include "peaks.h"
 #include "detector.h"
-#include "sfac.h"
 #include "filters.h"
 #include "reflections.h"
 #include "thread-pool.h"
diff --git a/src/mosflm.c b/src/mosflm.c
index fe6cfca7..593b04eb 100644
--- a/src/mosflm.c
+++ b/src/mosflm.c
@@ -61,7 +61,6 @@
 #include "image.h"
 #include "mosflm.h"
 #include "utils.h"
-#include "sfac.h"
 #include "peaks.h"
 
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 18b009eb..6d3c07ab 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -29,7 +29,6 @@
 #include "hdf5-file.h"
 #include "detector.h"
 #include "peaks.h"
-#include "sfac.h"
 #include "reflections.h"
 #include "beam-parameters.h"
 #include "symmetry.h"
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 90428426..6ff40159 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -23,7 +23,6 @@
 
 #include "utils.h"
 #include "statistics.h"
-#include "sfac.h"
 #include "reflist-utils.h"
 #include "symmetry.h"
 #include "stream.h"
diff --git a/src/sfac.c b/src/sfac.c
deleted file mode 100644
index 8b69f2e3..00000000
--- a/src/sfac.c
+++ /dev/null
@@ -1,542 +0,0 @@
-/*
- * sfac.c
- *
- * Scattering factors
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <complex.h>
-#include <string.h>
-#include <ctype.h>
-
-#include "utils.h"
-#include "sfac.h"
-
-
-#define N_MEMO 1024
-
-
-/* Look up f1 and f2 for this atom at this energy (in J/photon) */
-static double complex get_f1f2(const char *n, double en)
-{
-	FILE *fh;
-	char filename[64];
-	char line[1024];
-	char *rval;
-	double last_E, last_f1, last_f2;
-	static char *memo_n[N_MEMO];
-	static int memo_eV[N_MEMO];
-	static double complex memo_res[N_MEMO];
-	static int n_memo = 0;
-	int eV;
-	int i;
-
-	eV = (int)rint(J_to_eV(en));
-
-	for ( i=0; i<n_memo; i++ ) {
-		if ( (memo_eV[i] == eV) && (strcmp(memo_n[i], n) == 0) ) {
-			return memo_res[i];
-		}
-	}
-
-	snprintf(filename, 63, DATADIR"/crystfel/%s.nff", n);
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-		ERROR("Couldn't open file '%s'\n", filename);
-		return 0.0;
-	}
-
-	en = J_to_eV(en);
-
-	/* Discard first line */
-	fgets(line, 1023, fh);
-
-	last_E = 0.0;
-	last_f1 = 0.0;
-	last_f2 = 0.0;
-	do {
-
-		int r;
-		double E, f1, f2;
-		float E_f, f1_f, f2_f;
-
-		rval = fgets(line, 1023, fh);
-
-		r = sscanf(line, "%f %f %f", &E_f, &f1_f, &f2_f);
-		if ( r != 3 ) {
-			STATUS("I couldn't understand a line in the f1f2 "
-			       "tables\n");
-			continue;
-		}
-		/* Promote to double precision */
-		E = E_f;  f1 = f1_f;  f2 = f2_f;
-
-		/* Find the first energy greater than the required value */
-		if ( E < en ) {
-			/* Store old values ready for interpolation*/
-			last_E = E;
-			last_f1 = f1;
-			last_f2 = f2;
-		} else {
-
-			/* Perform (linear) interpolation */
-			double f;
-			double actual_f1, actual_f2;
-			double complex res;
-
-			f = (en - last_E) / (E - last_E);
-
-			actual_f1 = last_f1 + f * (f1 - last_f1);
-			actual_f2 = last_f2 + f * (f2 - last_f2);
-
-			fclose(fh);
-
-			res = actual_f1 + I*actual_f2;
-
-			memo_n[n_memo] = strdup(n);
-			memo_eV[n_memo] = eV;
-			memo_res[n_memo++] = res;
-			n_memo = n_memo % N_MEMO;
-
-			return res;
-
-		}
-
-	} while ( rval != NULL );
-
-	fclose(fh);
-
-	ERROR("Couldn't find scattering factors for '%s' at %f eV!\n", n, en);
-
-	return 0.0;
-}
-
-
-struct waas_kirf_factors {
-	char *n;
-	float a1;
-	float a2;
-	float a3;
-	float a4;
-	float a5;
-	float b1;
-	float b2;
-	float b3;
-	float b4;
-	float b5;
-	float c;
-};
-
-/* s = sin(theta)/lambda in metres^-1*/
-static double get_waas_kirf(const char *n, double s)
-{
-	FILE *fh;
-	char *rval;
-	double f;
-	float a1, a2, a3, a4, a5, c, b1, b2, b3, b4, b5;
-	double s2;
-	int i;
-	static struct waas_kirf_factors waaskirfcache[N_MEMO];
-	int found;
-	static int n_waaskirfcache = 0;
-
-	found = 0;
-	for ( i=0; i<n_waaskirfcache; i++ ) {
-		if ( strcmp(waaskirfcache[i].n, n) == 0 ) {
-
-			found = 1;
-
-			a1 = waaskirfcache[n_waaskirfcache].a1;
-			a2 = waaskirfcache[n_waaskirfcache].a2;
-			a3 = waaskirfcache[n_waaskirfcache].a3;
-			a4 = waaskirfcache[n_waaskirfcache].a4;
-			a5 = waaskirfcache[n_waaskirfcache].a5;
-			b1 = waaskirfcache[n_waaskirfcache].b1;
-			b2 = waaskirfcache[n_waaskirfcache].b2;
-			b3 = waaskirfcache[n_waaskirfcache].b3;
-			b4 = waaskirfcache[n_waaskirfcache].b4;
-			b5 = waaskirfcache[n_waaskirfcache].b5;
-			c = waaskirfcache[n_waaskirfcache].c;
-
-			break;
-		}
-	}
-
-	if ( !found ) {
-
-		fh = fopen(DATADIR"/crystfel/f0_WaasKirf.dat", "r");
-		if ( fh == NULL ) {
-			ERROR("Couldn't open f0_WaasKirf.dat\n");
-			return 0.0;
-		}
-
-		do {
-
-			int r;
-			char line[1024];
-			char sp[1024];
-			int Z;
-
-			rval = fgets(line, 1023, fh);
-
-			if ( (line[0] != '#') || (line[1] != 'S') ) continue;
-
-			r = sscanf(line, "#S  %i  %s", &Z, sp);
-			if ( (r != 2) || (strcmp(sp, n) != 0) ) continue;
-
-			/* Skip two lines */
-			fgets(line, 1023, fh);
-			fgets(line, 1023, fh);
-
-			/* Read scattering coefficients */
-			rval = fgets(line, 1023, fh);
-			r = sscanf(line,
-			             "  %f  %f  %f  %f  %f  %f  %f %f %f %f %f",
-			            &a1, &a2, &a3, &a4, &a5, &c,
-			            &b1, &b2, &b3, &b4, &b5);
-			if ( r != 11 ) {
-				ERROR("Couldn't read scattering "
-				      "factors (WaasKirf)\n");
-				return 0.0;
-			}
-
-			break;
-
-		} while ( rval != NULL );
-
-		fclose(fh);
-
-		waaskirfcache[n_waaskirfcache].a1 = a1;
-		waaskirfcache[n_waaskirfcache].a2 = a2;
-		waaskirfcache[n_waaskirfcache].a3 = a3;
-		waaskirfcache[n_waaskirfcache].a4 = a4;
-		waaskirfcache[n_waaskirfcache].a5 = a5;
-		waaskirfcache[n_waaskirfcache].c = c;
-		waaskirfcache[n_waaskirfcache].b1 = b1;
-		waaskirfcache[n_waaskirfcache].b2 = b2;
-		waaskirfcache[n_waaskirfcache].b3 = b3;
-		waaskirfcache[n_waaskirfcache].b4 = b4;
-		waaskirfcache[n_waaskirfcache].b5 = b5;
-		waaskirfcache[n_waaskirfcache++].n = strdup(n);
-		n_waaskirfcache = n_waaskirfcache % N_MEMO;  /* Unlikely */
-
-	}
-
-	s2 = pow(s/1e10, 2.0);  /* s2 is s squared in Angstroms^-2 */
-	f = c + a1*exp(-b1*s2) + a2*exp(-b2*s2) + a3*exp(-b3*s2)
-	      + a4*exp(-b4*s2) + a5*exp(-b5*s2);
-
-	return f;
-}
-
-
-/* Get complex scattering factors for element 'n' at energy 'en' (J/photon),
- * at resolution 's' = sin(theta)/lambda (in m^-1) */
-static double complex get_sfac(const char *n, double s, double en)
-{
-	double complex f1f2;
-	double fq;
-
-	/* Anomalous part (point-like K-shell electrons only) */
-	f1f2 = get_f1f2(n, en);
-
-	/* Atomic form factor part (not complex-valued) */
-	fq = get_waas_kirf(n, s) - get_waas_kirf(n, 0.0);
-
-	return fq + f1f2;
-}
-
-
-static void centre_molecule(struct molecule *mol)
-{
-	int i;
-	double tx = 0.0;
-	double ty = 0.0;
-	double tz = 0.0;
-	double total = 0.0;
-
-	/* Atoms are grouped by species for faster calculation */
-	for ( i=0; i<mol->n_species; i++ ) {
-
-		struct mol_species *spec;
-		int j;
-
-		spec = mol->species[i];
-
-		for ( j=0; j<spec->n_atoms; j++ ) {
-
-			double sfac = get_waas_kirf(spec->species, 0.0);
-
-			tx += spec->x[j] * sfac;
-			ty += spec->y[j] * sfac;
-			tz += spec->z[j] * sfac;
-			total += sfac;
-
-		}
-
-	}
-
-	mol->xc = tx / total;
-	mol->yc = ty / total;
-	mol->zc = tz / total;
-
-	for ( i=0; i<mol->n_species; i++ ) {
-
-		struct mol_species *spec;
-		int j;
-
-		spec = mol->species[i];
-
-		for ( j=0; j<spec->n_atoms; j++ ) {
-
-			spec->x[j] -= mol->xc;
-			spec->y[j] -= mol->yc;
-			spec->z[j] -= mol->zc;
-
-		}
-
-	}
-
-	STATUS("The atoms were shifted by %5.3f, %5.3f, %5.3f nm\n",
-	       mol->xc*1e9, mol->yc*1e9, mol->zc*1e9);
-}
-
-
-/* Load PDB file into a memory format suitable for efficient(ish) structure
- * factor calculation */
-struct molecule *load_molecule(const char *filename)
-{
-	struct molecule *mol;
-	FILE *fh;
-	char line[1024];
-	char *rval;
-	int i;
-
-	mol = malloc(sizeof(struct molecule));
-	if ( mol == NULL ) return NULL;
-	mol->n_species = 0;
-	mol->reflections = NULL;
-	mol->cell = NULL;
-
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-                ERROR("Couldn't open PDB file\n");
-                return NULL;
-        }
-
-	do {
-
-		int j, r;
-		int done = 0;
-		char xs[8];
-		char ys[8];
-		char zs[8];
-		char occs[6];
-		char Bs[6];
-		char el[3];
-		float xf, yf, zf, occf, Bf;
-		double x, y, z, occ, B;
-
-		rval = fgets(line, 1023, fh);
-
-		/* Only interested in atoms */
-		if ( (strncmp(line, "HETATM", 6) != 0)
-		  && (strncmp(line, "ATOM  ", 6) != 0) ) continue;
-
-		chomp(line);
-		if ( strlen(line) != 80 ) {
-			STATUS("Line does not have the correct length (%i):\n",
-			       (int)strlen(line));
-			STATUS("'%s'\n", line);
-			continue;
-		}
-
-		/* Separate out fixed-width fields */
-		memcpy(xs, line+30, 5);    xs[7] = '\0';
-		memcpy(ys, line+38, 5);    ys[7] = '\0';
-		memcpy(zs, line+46, 5);    zs[7] = '\0';
-		memcpy(occs, line+54, 5);  occs[5] = '\0';
-		memcpy(Bs, line+60, 5);    Bs[5] = '\0';
-		memcpy(el, line+76, 2);    el[2] = '\0';
-
-		/* Convert fields */
-		r = sscanf(xs, "%f", &xf);
-		r += sscanf(ys, "%f", &yf);
-		r += sscanf(zs, "%f", &zf);
-		r += sscanf(occs, "%f", &occf);
-		r += sscanf(Bs, "%f", &Bf);
-		if ( el[0] == ' ' ) {
-			el[0] = el[1];
-			el[1] = '\0';
-		} else {
-			el[1] = tolower(el[1]);
-		}
-
-		/* Promote to double precision */
-		x = xf;  y = yf;  z = zf;  occ = occf;  B = Bf;
-
-		for ( j=0; j<mol->n_species; j++ ) {
-
-			struct mol_species *spec;
-			int n;
-
-			spec = mol->species[j];
-
-			if ( strcmp(spec->species, el) != 0 ) continue;
-
-			n = mol->species[j]->n_atoms;
-
-			spec->x[n] = x*1.0e-10;  /* Convert to m */
-			spec->y[n] = y*1.0e-10;
-			spec->z[n] = z*1.0e-10;
-			spec->occ[n] = occ;
-			spec->B[n] = B*1.0e-20;  /* Convert to m^2 */
-			mol->species[j]->n_atoms++;
-			if ( mol->species[j]->n_atoms == MAX_ATOMS ) {
-				ERROR("Too many atoms of type '%s'!\n", el);
-				exit(1);
-			}
-
-			done = 1;
-
-		}
-
-		if ( !done ) {
-
-			/* Need to create record for this species */
-			struct mol_species *spec;
-
-			spec = malloc(sizeof(struct mol_species));
-
-			memcpy(spec->species, el, 1+strlen(el));
-			spec->x[0] = x*1.0e-10;  /* Convert to m */
-			spec->y[0] = y*1.0e-10;
-			spec->z[0] = z*1.0e-10;
-			spec->occ[0] = occ;
-			spec->B[0] = B*1.0e-20;  /* Convert to m^2 */
-			spec->n_atoms = 1;
-
-			mol->species[mol->n_species] = spec;
-			mol->n_species++;
-
-		}
-
-	} while ( rval != NULL );
-
-	fclose(fh);
-
-	centre_molecule(mol);
-
-	STATUS("There are %i species:\n", mol->n_species); fflush(stderr);
-	for ( i=0; i<mol->n_species; i++ ) {
-		STATUS("%3s : %6i\n", mol->species[i]->species,
-		       mol->species[i]->n_atoms);
-	}
-
-	mol->cell = load_cell_from_pdb(filename);
-	if ( mol->cell == NULL ) {
-		ERROR("No unit cell found in PDB file\n");
-		return NULL;
-	}
-
-	return mol;
-}
-
-
-void free_molecule(struct molecule *mol)
-{
-	int i;
-
-	for ( i=0; i<mol->n_species; i++ ) {
-		free(mol->species[i]);
-	}
-
-	free(mol->cell);
-	free(mol);
-}
-
-
-double *get_reflections(struct molecule *mol, double en, double res,
-                        double *phases, ReflItemList *items)
-{
-	double *reflections;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	signed int h, k, l;
-	const int do_thermal = 1;
-
-	cell_get_reciprocal(mol->cell, &asx, &asy, &asz,
-	                               &bsx, &bsy, &bsz,
-	                               &csx, &csy, &csz);
-
-	reflections = new_list_intensity();
-
-	for ( h=-INDMAX; h<=INDMAX; h++ ) {
-	for ( k=-INDMAX; k<=INDMAX; k++ ) {
-	for ( l=-INDMAX; l<=INDMAX; l++ ) {
-
-		double complex F = 0.0;
-		int i;
-		double s, oneoverd;
-
-		/* We need sin(theta)/lambda = 1/2d */
-		s = resolution(mol->cell, h, k, l);
-		oneoverd = 2.0 * s;
-		if ( oneoverd > res ) continue;
-
-		/* Atoms are grouped by species for faster calculation */
-		for ( i=0; i<mol->n_species; i++ ) {
-
-			double complex sfac;
-			double complex contrib = 0.0;
-			struct mol_species *spec;
-			int j;
-
-			spec = mol->species[i];
-
-			for ( j=0; j<spec->n_atoms; j++ ) {
-
-				double ph, u, v, w;
-				double complex cpart;
-
-				u = h*asx + k*bsx + l*csx;
-				v = h*asy + k*bsy + l*csy;
-				w = h*asz + k*bsz + l*csz;
-
-				ph = u*spec->x[j] + v*spec->y[j] + w*spec->z[j];
-
-				/* Conversion from revolutions to radians */
-				cpart = cexp(-2.0*M_PI*I*ph);
-				if ( do_thermal ) {
-					cpart *= exp(-2.0 * spec->B[j] * s * s);
-				}
-				contrib += cpart;
-
-			}
-
-			sfac = get_sfac(spec->species, s, en);
-			F += sfac * contrib;
-
-		}
-
-		set_intensity(reflections, h, k, l, pow(cabs(F), 2.0));
-		if ( phases != NULL ) set_phase(phases, h, k, l, carg(F));
-		if ( items != NULL ) add_item(items, h, k, l);
-
-	}
-	progress_bar((k+INDMAX)+IDIM*(h+INDMAX), IDIM*IDIM-1,
-	             "Calculating reflection intensities");
-	}
-	}
-
-	return reflections;
-}
diff --git a/src/sfac.h b/src/sfac.h
deleted file mode 100644
index e05b3ddf..00000000
--- a/src/sfac.h
+++ /dev/null
@@ -1,64 +0,0 @@
-/*
- * sfac.h
- *
- * Scattering factors
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef SFAC_H
-#define SFAC_H
-
-#include <complex.h>
-
-#include "cell.h"
-#include "utils.h"
-
-
-#define MAX_ATOMS (128*1024)
-
-
-struct mol_species
-{
-	char species[4];    /* Species name */
-	int n_atoms;        /* Number of atoms of this species */
-
-	double x[MAX_ATOMS];
-	double y[MAX_ATOMS];
-	double z[MAX_ATOMS];
-	double occ[MAX_ATOMS];
-	double B[MAX_ATOMS];
-};
-
-
-struct molecule
-{
-	int n_species;
-	struct mol_species *species[32];
-
-	/* Unit cell */
-	UnitCell *cell;
-
-	/* Reflection intensities at Bragg positions */
-	double complex *reflections;
-
-	/* Offset to molecule's centre of scattering power */
-	double xc;
-	double yc;
-	double zc;
-};
-
-extern struct molecule *load_molecule(const char *filename);
-extern void free_molecule(struct molecule *mol);
-extern double *get_reflections(struct molecule *mol, double en, double res,
-                               double *phases, ReflItemList *items);
-
-
-#endif	/* SFAC_H */