aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/diffraction.c104
-rw-r--r--src/diffraction.h2
-rw-r--r--src/main.c11
-rw-r--r--src/sfac.c74
-rw-r--r--src/sfac.h10
-rw-r--r--src/utils.h51
6 files changed, 203 insertions, 49 deletions
diff --git a/src/diffraction.c b/src/diffraction.c
index 4ed13709..e112e7ce 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -76,55 +76,78 @@ static double lattice_factor(struct threevec q, double ax, double ay, double az,
}
-/* Return structure factor for molecule 'mol' at energy 'en' (J/photon) at
- * scattering vector 'q' */
+/* Look up the structure factor for the nearest Bragg condition */
static double complex molecule_factor(struct molecule *mol, struct threevec q,
- double en)
+ double ax, double ay, double az,
+ double bx, double by, double bz,
+ double cx, double cy, double cz)
{
- int i;
- double complex F = 0.0;
- double s;
+ double hd, kd, ld;
+ signed int h, k, l;
- /* s = sin(theta)/lambda = 1/2d = (1/d)/2.0 */
- s = modulus(q.u, q.v, q.w) / 2.0;
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * cz;
+ h = (signed int)nearbyint(hd);
+ k = (signed int)nearbyint(kd);
+ l = (signed int)nearbyint(ld);
- /* 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];
+ return get_integral(mol->reflections, h, k, l);
+}
- for ( j=0; j<spec->n_atoms; j++ ) {
- double ph;
+static double complex water_factor(struct threevec q, double en)
+{
+ return 0.0;
+}
- ph = q.u*spec->x[j] + q.v*spec->y[j] + q.w*spec->z[j];
- /* Conversion from revolutions to radians is required */
- contrib += cos(2.0*M_PI*ph) + I*sin(2.0*M_PI*ph);
+static void get_reflections_cached(struct molecule *mol, double en)
+{
+ char s[1024];
+ FILE *fh;
+ size_t r;
+
+ snprintf(s, 1023, "reflections-%ieV.cache", (int)J_to_eV(en));
+ fh = fopen(s, "rb");
+ if ( fh == NULL ) {
+ printf("No cache file found (looked for %s)\n", s);
+ goto calc;
+ }
- }
+ mol->reflections = reflist_new();
+ r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh);
+ if ( r < IDIM*IDIM*IDIM ) {
+ printf("Found cache file (%s), but failed to read.\n", s);
+ goto calc;
+ }
- sfac = get_sfac(spec->species, s, en);
- F += sfac * contrib * exp(-2.0 * spec->B[j] * s);
+ printf("Read structure factors (at Bragg positions) from %s\n", s);
+ return;
+calc:
+ printf("Calculating structure factors at Bragg positions...\n");
+ mol->reflections = get_reflections(mol, en);
+ fh = fopen(s, "wb");
+ if ( fh == NULL ) {
+ printf("Failed to write cache file (%s)\n", s);
+ return;
}
- return F;
-}
-
+ r = fwrite(mol->reflections, sizeof(double complex),
+ IDIM*IDIM*IDIM, fh);
+ if ( r < IDIM*IDIM*IDIM ) {
+ printf("Failed to write cache file (%s)\n", s);
+ return;
+ }
+ fclose(fh);
-static double complex water_factor(struct threevec q, double en)
-{
- return 0.0;
+ printf("Successfully saved structure factors at Bragg positions to"
+ " file %s\n", s);
}
-void get_diffraction(struct image *image, UnitCell *cell)
+void get_diffraction(struct image *image)
{
int x, y;
double ax, ay, az;
@@ -139,13 +162,17 @@ void get_diffraction(struct image *image, UnitCell *cell)
if ( image->molecule == NULL ) return;
}
- cell_get_cartesian(cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
+ cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
image->sfacs = malloc(image->width * image->height
* sizeof(double complex));
+ if ( image->molecule->reflections == NULL ) {
+ get_reflections_cached(image->molecule, image->xray_energy);
+ }
+
progress_bar(0, image->width-1);
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
@@ -154,12 +181,13 @@ void get_diffraction(struct image *image, UnitCell *cell)
double complex f_molecule;
double complex f_water;
struct threevec q;
+ double complex val;
q = image->qvecs[x + image->width*y];
f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz);
f_molecule = molecule_factor(image->molecule, q,
- image->xray_energy);
+ ax,ay,az,bx,by,bz,cx,cy,cz);
/* Nasty approximation follows */
if ( y == image->height/2 ) {
@@ -168,8 +196,8 @@ void get_diffraction(struct image *image, UnitCell *cell)
f_water = 0.0;
}
- image->sfacs[x + image->width*y] = (f_lattice * f_molecule)
- + f_water;
+ val = (f_molecule) + f_water;
+ image->sfacs[x + image->width*y] = val;
}
progress_bar(x, image->width-1);
diff --git a/src/diffraction.h b/src/diffraction.h
index c32e62c2..112818bc 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -19,6 +19,6 @@
#include "image.h"
#include "cell.h"
-extern void get_diffraction(struct image *image, UnitCell *cell);
+extern void get_diffraction(struct image *image);
#endif /* DIFFRACTION_H */
diff --git a/src/main.c b/src/main.c
index b7755ce6..f3aeadfc 100644
--- a/src/main.c
+++ b/src/main.c
@@ -41,7 +41,6 @@ static void main_show_help(const char *s)
int main(int argc, char *argv[])
{
int c, done;
- UnitCell *cell;
struct image image;
char filename[1024];
int number = 1;
@@ -59,14 +58,6 @@ int main(int argc, char *argv[])
}
- /* Define unit cell */
- cell = cell_new_from_parameters(28.10e-9,
- 28.10e-9,
- 16.52e-9,
- deg2rad(90.0),
- deg2rad(90.0),
- deg2rad(120.0));
-
/* Define image parameters */
image.width = 1024;
image.height = 1024;
@@ -124,7 +115,7 @@ again:
image.twotheta = NULL;
image.hdr = NULL;
- get_diffraction(&image, cell);
+ get_diffraction(&image);
record_image(&image);
snprintf(filename, 1023, "results/sim-%i.h5", number);
diff --git a/src/sfac.c b/src/sfac.c
index 00d6c69d..c13ed3b8 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -272,6 +272,15 @@ struct molecule *load_molecule()
mol = malloc(sizeof(struct molecule));
if ( mol == NULL ) return NULL;
mol->n_species = 0;
+ mol->reflections = NULL;
+
+ /* FIXME: Read cell from file */
+ mol->cell = cell_new_from_parameters(28.10e-9,
+ 28.10e-9,
+ 16.52e-9,
+ deg2rad(90.0),
+ deg2rad(90.0),
+ deg2rad(120.0));
fh = fopen("molecule.pdb", "r");
if ( fh == NULL ) {
@@ -370,3 +379,68 @@ struct molecule *load_molecule()
return mol;
}
+
+
+double complex *get_reflections(struct molecule *mol, double en)
+{
+ double complex *reflections;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ signed int h, k, l;
+
+ cell_get_reciprocal(mol->cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ reflections = reflist_new();
+
+ for ( h=-INDMAX; h<=INDMAX; h++ ) {
+ for ( k=-INDMAX; k<=INDMAX; k++ ) {
+ for ( l=-INDMAX; l<=INDMAX; l++ ) {
+
+ double complex F;
+ int i;
+ double s;
+
+ /* 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;
+
+ 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 */
+ contrib += cos(-2.0*M_PI*ph)
+ + I*sin(-2.0*M_PI*ph);
+
+ }
+
+ sfac = get_sfac(spec->species, s, en);
+ F += sfac * contrib * exp(-2.0 * spec->B[j] * s);
+
+ }
+
+ integrate_reflection(reflections, h, k, l, F);
+
+ }
+ }
+ progress_bar((h+INDMAX+1), 2*INDMAX);
+ }
+ printf("\n");
+
+ return reflections;
+}
diff --git a/src/sfac.h b/src/sfac.h
index f88729f4..fc4b2bd9 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -18,6 +18,8 @@
#include <complex.h>
+#include "cell.h"
+
struct mol_species
{
@@ -37,6 +39,13 @@ 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;
@@ -45,5 +54,6 @@ struct molecule
extern double complex get_sfac(const char *n, double s, double en);
extern struct molecule *load_molecule(void);
+extern double complex *get_reflections(struct molecule *mol, double en);
#endif /* SFAC_H */
diff --git a/src/utils.h b/src/utils.h
index 7940d77f..ca653973 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -17,6 +17,9 @@
#endif
#include <math.h>
+#include <complex.h>
+#include <string.h>
+#include <stdlib.h>
/* Electron charge in C */
@@ -31,6 +34,10 @@
/* Thomson scattering length (m) */
#define THOMSON_LENGTH (2.81794e-15)
+/* Maxmimum index to go up to */
+#define INDMAX 20
+#define IDIM (INDMAX*2 +1)
+
extern unsigned int biggest(signed int a, signed int b);
extern unsigned int smallest(signed int a, signed int b);
@@ -69,4 +76,48 @@ extern void progress_bar(int val, int total);
/* Joules to eV */
#define J_to_eV(a) ((a)/ELECTRON_CHARGE)
+
+static inline void integrate_reflection(double complex *ref, signed int h,
+ signed int k, signed int l,
+ double complex i)
+{
+ int idx;
+
+ /* Not interested in central beam */
+ if ( (h==0) && (k==0) && (l==0) ) return;
+
+ if ( h < 0 ) h += IDIM;
+ if ( k < 0 ) k += IDIM;
+ if ( l < 0 ) l += IDIM;
+
+ idx = h + (IDIM*k) + (IDIM*IDIM*l);
+ ref[idx] += i;
+}
+
+
+static inline double complex get_integral(double complex *ref, signed int h,
+ signed int k, signed int l)
+{
+ int idx;
+
+ if ( h < 0 ) h += IDIM;
+ if ( k < 0 ) k += IDIM;
+ if ( l < 0 ) l += IDIM;
+
+ idx = h + (IDIM*k) + (IDIM*IDIM*l);
+ return ref[idx];
+}
+
+
+static inline double complex *reflist_new(void)
+{
+ double complex *r;
+ size_t r_size;
+ r_size = IDIM*IDIM*IDIM*sizeof(double complex);
+ r = malloc(r_size);
+ memset(r, 0, r_size);
+ return r;
+}
+
+
#endif /* UTILS_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-29 10:26:28 +0000