Don't try to render PDBs, part III: tidy up and fix

src/compare_hkl.c

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
 	double *ref1;
 	double *ref2;
 	double *out;
-	struct molecule *mol;
+	UnitCell *cell;
 	char *outfile = NULL;
 	char *afile = NULL;
 	char *bfile = NULL;
@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	mol = load_molecule();
+	cell = load_cell_from_pdb("molecule.pdb");
 	ref1 = read_reflections(afile);
 	if ( ref1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", afile);
@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
 	}
 	}
 
-	write_reflections(outfile, NULL, out, 1, mol->cell);
+	write_reflections(outfile, NULL, out, 1, cell);
 
 	return 0;
 }

src/diffraction.c

@@ -223,7 +223,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 			                              na, nb, nc);
 
 			if ( intensities == NULL ) {
-				I_molecule = 10000.0;
+				I_molecule = 1.0e10;
 			} else {
 				I_molecule = molecule_factor(intensities, q,
 			                            ax,ay,az,bx,by,bz,cx,cy,cz);

src/get_hkl.c

@@ -153,7 +153,7 @@ int main(int argc, char *argv[])
 	}
 
 	mol = load_molecule();
-	get_reflections_cached(mol, eV_to_J(1790.0));
+	ideal_ref = get_reflections(mol, eV_to_J(1790.0));
 
 	counts = new_list_count();
 

src/indexamajig.c

@@ -208,7 +208,7 @@ int main(int argc, char *argv[])
 		{"no-match",           0, &config_nomatch,     1},
 		{"verbose",            0, &config_verbose,     1},
 		{"alternate",          0, &config_alternate,   1},
-		{"intensities",        0, NULL,               'q'},
+		{"intensities",        1, NULL,               'q'},
 		{0, 0, NULL, 0}
 	};
 

src/pattern_sim.c

@@ -175,7 +175,7 @@ int main(int argc, char *argv[])
 		{"no-images",          0, &config_noimages,    1},
 		{"no-water",           0, &config_nowater,     1},
 		{"no-noise",           0, &config_nonoise,     1},
-		{"intensities",        0, NULL,               'i'},
+		{"intensities",        1, NULL,               'i'},
 		{"powder",             0, &config_powder,      1},
 		{0, 0, NULL, 0}
 	};

src/process_hkl.c

@@ -187,7 +187,7 @@ int main(int argc, char *argv[])
 	double *ref, *trueref = NULL;
 	unsigned int *counts;
 	char *rval;
-	struct molecule *mol = NULL;
+	UnitCell *cell;
 	int config_maxonly = 0;
 	int config_every = 1000;
 	int config_rvsq = 0;
@@ -276,8 +276,7 @@ int main(int argc, char *argv[])
 
 	ref = new_list_intensity();
 	counts = new_list_count();
-
-	mol = load_molecule();
+	cell = load_cell_from_pdb("molecule.pdb");
 
 	if ( strcmp(filename, "-") == 0 ) {
 		fh = stdin;
@@ -307,7 +306,7 @@ int main(int argc, char *argv[])
 
 			if (config_every && (n_patterns % config_every == 0)) {
 				process_reflections(ref, trueref, counts,
-				                    n_patterns, mol->cell,
+				                    n_patterns, cell,
 				                    config_rvsq,
 				                    config_zoneaxis);
 			}
@@ -341,20 +340,16 @@ int main(int argc, char *argv[])
 	fclose(fh);
 
 	if ( trueref != NULL ) {
-		process_reflections(ref, trueref, counts, n_patterns, mol->cell,
+		process_reflections(ref, trueref, counts, n_patterns, cell,
 		                    config_rvsq, config_zoneaxis);
 	}
 
 	if ( output != NULL ) {
-		UnitCell *cell = NULL;
-		if ( mol != NULL ) cell = mol->cell;
 		write_reflections(output, counts, ref, 0, cell);
 	}
 
 	if ( config_zoneaxis ) {
 		char name[64];
-		UnitCell *cell = NULL;
-		if ( mol != NULL ) cell = mol->cell;
 		snprintf(name, 63, "ZA-%u.dat", n_patterns);
 		write_reflections(name, counts, ref, 1, cell);
 	}

src/sfac.c

@@ -304,8 +304,8 @@ static void centre_molecule(struct molecule *mol)
 
 	}
 
-	//STATUS("Molecule was shifted by %5.3f, %5.3f, %5.3f nm\n",
-	//       mol->xc*1e9, mol->yc*1e9, mol->zc*1e9);
+	STATUS("The atoms were shifted by %5.3f, %5.3f, %5.3f nm\n",
+	       mol->xc*1e9, mol->yc*1e9, mol->zc*1e9);
 }
 
 
@@ -462,7 +462,7 @@ struct molecule *load_molecule()
 
 	centre_molecule(mol);
 
-	STATUS("There are %i species\n", mol->n_species); fflush(stderr);
+	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);
@@ -490,31 +490,9 @@ void free_molecule(struct molecule *mol)
 }
 
 
-struct molecule *copy_molecule(struct molecule *orig)
+double *get_reflections(struct molecule *mol, double en)
 {
-	struct molecule *new;
-	int i;
-
-	new = malloc(sizeof(*new));
-	*new = *orig;
-
-	/* Now sort out pointers */
-	for ( i=0; i<orig->n_species; i++ ) {
-		new->species[i] = malloc(sizeof(struct mol_species));
-		memcpy(new->species[i], orig->species[i],
-		       sizeof(struct mol_species));
-	}
-
-	new->cell = cell_new_from_cell(orig->cell);
-	new->reflections = NULL;
-
-	return new;
-}
-
-
-double complex *get_reflections(struct molecule *mol, double en)
-{
-	double complex *reflections;
+	double *reflections;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
@@ -525,7 +503,7 @@ double complex *get_reflections(struct molecule *mol, double en)
 	                               &bsx, &bsy, &bsz,
 	                               &csx, &csy, &csz);
 
-	reflections = new_list_sfac();
+	reflections = new_list_intensity();
 
 	for ( h=-INDMAX; h<=INDMAX; h++ ) {
 	for ( k=-INDMAX; k<=INDMAX; k++ ) {
@@ -572,61 +550,13 @@ double complex *get_reflections(struct molecule *mol, double en)
 
 		}
 
-		set_sfac(reflections, h, k, l, F);
+		set_intensity(reflections, h, k, l, pow(cabs(F), 2.0));
 
 	}
 	progress_bar((k+INDMAX)+IDIM*(h+INDMAX), IDIM*IDIM-1,
-	             "Calculating structure factors");
+	             "Calculating reflection intensities");
 	}
 	}
 
 	return reflections;
 }
-
-
-void get_reflections_cached(struct molecule *mol, double en)
-{
-	char s[1024];
-	FILE *fh;
-	size_t r;
-
-	/* Add 0.5 to improve rounding */
-	snprintf(s, 1023, "reflections-%ieV.cache", (int)(J_to_eV(en)+0.5));
-	fh = fopen(s, "rb");
-	if ( fh == NULL ) {
-		STATUS("No cache file found (looked for %s)\n", s);
-		goto calc;
-	}
-
-	mol->reflections = new_list_sfac();
-	r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh);
-	if ( r < IDIM*IDIM*IDIM ) {
-		STATUS("Found cache file (%s), but failed to read"
-		       " (got %lli items, wanted %i).\n",
-		       s, (long long)r, IDIM*IDIM*IDIM);
-		goto calc;
-	}
-
-	STATUS("Read structure factors (at Bragg positions) from %s\n", s);
-	return;
-
-calc:
-	STATUS("Calculating structure factors at Bragg positions...\n");
-	mol->reflections = get_reflections(mol, en);
-	fh = fopen(s, "wb");
-	if ( fh == NULL ) {
-		STATUS("Failed to write cache file (%s)\n", s);
-		return;
-	}
-
-	r = fwrite(mol->reflections, sizeof(double complex),
-	           IDIM*IDIM*IDIM, fh);
-	if ( r < IDIM*IDIM*IDIM ) {
-		STATUS("Failed to write cache file (%s)\n", s);
-		return;
-	}
-	fclose(fh);
-
-	STATUS("Successfully saved structure factors at Bragg positions to"
-	       " file %s\n", s);
-}

src/sfac.h

@@ -55,12 +55,12 @@ struct molecule
 };
 
 
+/* This is so that the water background calculation can use it */
 extern double complex get_sfac(const char *n, double s, double en);
+
 extern struct molecule *load_molecule(void);
 extern void free_molecule(struct molecule *mol);
-extern struct molecule *copy_molecule(struct molecule *orig);
-extern double complex *get_reflections(struct molecule *mol, double en);
-extern void get_reflections_cached(struct molecule *mol, double en);
+extern double *get_reflections(struct molecule *mol, double en);
 
 
 #endif	/* SFAC_H */