Skip to content
GitLab
Commit 6a9811f6 authored by Thomas White's avatar Thomas White
Browse files

Don't try to render PDBs, part I

parent a6f27929
Hide whitespace changes
Inline Side-by-side
data/diffraction.cl
View file @ 6a9811f6
......@@ -101,7 +101,7 @@ float lattice_factor(float16 cell, float4 q,
}
float2 get_sfac(global float2 *sfacs, float16 cell, float4 q)
float get_intensity(global float *intensities, float16 cell, float4 q)
{
float hf, kf, lf;
int h, k, l;
......@@ -117,25 +117,25 @@ float2 get_sfac(global float2 *sfacs, float16 cell, float4 q)
/* Return a silly value if indices are out of range */
if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
return (float2)(100000.0, 0.0);
return 1000000.0;
}
h = (h>=0) ? h : h+IDIM;
k = (k>=0) ? k : k+IDIM;
l = (l>=0) ? l : l+IDIM;
if ( (h>=IDIM) || (k>=IDIM) || (l>=IDIM) ) return (float2)(100000.0, 0.0);
if ( (h>=IDIM) || (k>=IDIM) || (l>=IDIM) ) return 1000000.0;
idx = h + (IDIM*k) + (IDIM*IDIM*l);
return sfacs[idx];
return intensities[idx];
}
kernel void diffraction(global float *diff, global float *tt, float klow,
int w, float cx, float cy,
float res, float clen, float16 cell,
global float2 *sfacs, float4 z,
global float *intensities, float4 z,
int xmin, int ymin, int sampling, local float *tmp,
float kstep,
read_only image2d_t func_a,
......@@ -145,8 +145,8 @@ kernel void diffraction(global float *diff, global float *tt, float klow,
float ttv;
const int x = get_global_id(0) + (xmin*sampling);
const int y = get_global_id(1) + (ymin*sampling);
float f_lattice;
float2 f_molecule;
float f_lattice, I_lattice;
float I_molecule;
float4 q;
const int lx = get_local_id(0);
const int ly = get_local_id(1);
......@@ -155,16 +155,15 @@ kernel void diffraction(global float *diff, global float *tt, float klow,
const int ax = x / sampling;
const int ay = y / sampling;
float intensity;
float2 val;
/* Calculate value */
q = get_q(x, y, cx, cy, res, clen, k, &ttv, z, sampling);
f_lattice = lattice_factor(cell, q, func_a, func_b, func_c);
f_molecule = get_sfac(sfacs, cell, q);
I_molecule = get_intensity(intensities, cell, q);
I_lattice = pow(f_lattice, 2.0f);
/* Write the value to local memory */
val = f_molecule * f_lattice;
intensity = pow(val.x, 2.0f) + pow(val.y, 2.0f);
intensity = I_molecule * I_lattice;
tmp[lx+sampling*ly+sampling*sampling*lb] = intensity;
/* Memory fence */
......
src/diffraction-gpu.c
View file @ 6a9811f6
......@@ -38,7 +38,7 @@ struct gpu_context
cl_command_queue cq;
cl_program prog;
cl_kernel kern;
cl_mem sfacs;
cl_mem intensities;
cl_mem tt;
size_t tt_size;
......@@ -183,7 +183,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
ERROR("Couldn't set arg 8: %s\n", clError(err));
return;
}
clSetKernelArg(gctx->kern, 9, sizeof(cl_mem), &gctx->sfacs);
clSetKernelArg(gctx->kern, 9, sizeof(cl_mem), &gctx->intensities);
if ( err != CL_SUCCESS ) {
ERROR("Couldn't set arg 9: %s\n", clError(err));
return;
......@@ -344,7 +344,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
/* Setup the OpenCL stuff, create buffers, load the structure factor table */
struct gpu_context *setup_gpu(int no_sfac, struct image *image,
struct molecule *molecule)
double *intensities)
{
struct gpu_context *gctx;
cl_uint nplat;
......@@ -352,21 +352,11 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
cl_context_properties prop[3];
cl_int err;
cl_device_id dev;
size_t sfac_size;
float *sfac_ptr;
size_t intensities_size;
float *intensities_ptr;
size_t maxwgsize;
int i;
if ( molecule == NULL ) return NULL;
/* Generate structure factors if required */
if ( !no_sfac ) {
if ( molecule->reflections == NULL ) {
get_reflections_cached(molecule,
ph_lambda_to_en(image->lambda));
}
}
STATUS("Setting up GPU..."); fflush(stderr);
err = clGetPlatformIDs(8, platforms, &nplat);
......@@ -411,28 +401,26 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
}
/* Create a single-precision version of the scattering factors */
sfac_size = IDIM*IDIM*IDIM*sizeof(cl_float)*2; /* complex */
sfac_ptr = malloc(sfac_size);
if ( !no_sfac ) {
intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float);
intensities_ptr = malloc(intensities_size);
if ( intensities != NULL ) {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
sfac_ptr[2*i+0] = creal(molecule->reflections[i]);
sfac_ptr[2*i+1] = cimag(molecule->reflections[i]);
intensities_ptr[i] = intensities[i];
}
} else {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
sfac_ptr[2*i+0] = 10000.0;
sfac_ptr[2*i+1] = 0.0;
intensities_ptr[i] = 10000.0;
}
}
gctx->sfacs = clCreateBuffer(gctx->ctx,
gctx->intensities = clCreateBuffer(gctx->ctx,
CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
sfac_size, sfac_ptr, &err);
intensities_size, intensities_ptr, &err);
if ( err != CL_SUCCESS ) {
ERROR("Couldn't allocate sfac memory\n");
ERROR("Couldn't allocate intensities memory\n");
free(gctx);
return NULL;
}
free(sfac_ptr);
free(intensities_ptr);
gctx->tt_size = image->width*image->height*sizeof(cl_float);
gctx->tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, gctx->tt_size,
......@@ -478,7 +466,7 @@ void cleanup_gpu(struct gpu_context *gctx)
clReleaseProgram(gctx->prog);
clReleaseMemObject(gctx->diff);
clReleaseMemObject(gctx->tt);
clReleaseMemObject(gctx->sfacs);
clReleaseMemObject(gctx->intensities);
/* Release LUTs */
for ( i=1; i<=gctx->max_sinc_lut; i++ ) {
......
src/diffraction-gpu.h
View file @ 6a9811f6
......@@ -26,7 +26,7 @@ struct gpu_context;
extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
int na, int nb, int nc);
extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
struct molecule *molecule);
double *intensities);
extern void cleanup_gpu(struct gpu_context *gctx);
#else
......@@ -39,7 +39,7 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
}
static struct gpu_context *setup_gpu(int no_sfac, struct image *image,
struct molecule *molecule)
double *intensities)
{
return NULL;
}
......
src/diffraction.c
View file @ 6a9811f6
......@@ -70,14 +70,14 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
/* Look up the structure factor for the nearest Bragg condition */
static double complex molecule_factor(struct molecule *mol, struct rvec q,
static double molecule_factor(double *intensities, struct rvec q,
double ax, double ay, double az,
double bx, double by, double bz,
double cx, double cy, double cz)
{
double hd, kd, ld;
signed int h, k, l;
double complex r;
double r;
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
......@@ -86,7 +86,7 @@ static double complex molecule_factor(struct molecule *mol, struct rvec q,
k = (signed int)rint(kd);
l = (signed int)rint(ld);
r = lookup_sfac(mol->reflections, h, k, l);
r = lookup_intensity(intensities, h, k, l);
return r;
}
......@@ -170,8 +170,8 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
}
void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
int do_water)
void get_diffraction(struct image *image, int na, int nb, int nc,
double *intensities, int do_water)
{
unsigned int xs, ys;
double ax, ay, az;
......@@ -188,13 +188,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
/* Allocate (and zero) the "diffraction array" */
image->data = calloc(image->width * image->height, sizeof(float));
if ( !no_sfac ) {
if ( image->molecule->reflections == NULL ) {
get_reflections_cached(image->molecule,
ph_lambda_to_en(image->lambda));
}
}
/* Needed later for Lorentz calculation */
image->twotheta = malloc(image->width * image->height * sizeof(double));
......@@ -205,8 +198,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
for ( xs=0; xs<image->width*SAMPLING; xs++ ) {
for ( ys=0; ys<image->height*SAMPLING; ys++ ) {
double f_lattice;
double complex f_molecule;
struct rvec q;
float twotheta;
double sw = 1.0/(SAMPLING*SAMPLING); /* Sample weight */
......@@ -220,8 +211,9 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
float k;
double kw = 1.0/BWSAMPLING;
double complex val;
double intensity;
double f_lattice, I_lattice;
double I_molecule;
/* Calculate k this time round */
k = klow + kstep * bwstep;
......@@ -233,15 +225,17 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
bx, by, bz,
cx, cy, cz,
na, nb, nc);
if ( no_sfac ) {
f_molecule = 10000.0;
if ( intensities == NULL ) {
I_molecule = 10000.0;
} else {
f_molecule = molecule_factor(image->molecule, q,
I_molecule = molecule_factor(intensities, q,
ax,ay,az,bx,by,bz,cx,cy,cz);
}
val = f_molecule * f_lattice;
intensity = sw * kw * pow(cabs(val), 2.0);
I_lattice = pow(f_lattice, 2.0);
intensity = sw * kw * I_lattice * I_molecule;
image->data[x + image->width*y] += intensity;
}
......
src/diffraction.h
View file @ 6a9811f6
......@@ -20,7 +20,8 @@
#include "cell.h"
extern void get_diffraction(struct image *image, int na, int nb, int nc,
int nosfac, int do_water);
double *intensities, int do_water);
extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
unsigned int sampling, float *ttp, float k);
#endif /* DIFFRACTION_H */
src/get_hkl.c
View file @ 6a9811f6
......@@ -110,6 +110,8 @@ int main(int argc, char *argv[])
char *output = NULL;
unsigned int *counts;
signed int h, k, l;
FILE *fh;
char *rval;
/* Long options */
const struct option longopts[] = {
......@@ -154,7 +156,7 @@ int main(int argc, char *argv[])
mol = load_molecule();
get_reflections_cached(mol, eV_to_J(1790.0));
ideal_ref = ideal_intensities(mol->reflections);
counts = new_list_count();
if ( template != NULL ) {
......
src/pattern_sim.c
View file @ 6a9811f6
......@@ -30,6 +30,7 @@
#include "detector.h"
#include "intensities.h"
#include "sfac.h"
#include "reflections.h"
static void show_help(const char *s)
......@@ -50,6 +51,8 @@ static void show_help(const char *s)
" (a new orientation will be used for each image).\n"
" --powder Write a summed pattern of all images simulated by\n"
" this invocation to results/integr.h5.\n"
" -i, --intensities=<file> Specify file containing reflection intensities\n"
" to use.\n"
"\n"
"By default, the simulation aims to be as accurate as possible. For greater\n"
"speed, or for testing, you can choose to disable certain things using the\n"
......@@ -57,7 +60,6 @@ static void show_help(const char *s)
"\n"
" --no-water Do not simulate water background.\n"
" --no-noise Do not calculate Poisson noise.\n"
" --no-sfac Pretend that all structure factors are 1.\n"
);
}
......@@ -78,21 +80,15 @@ static void show_details()
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
"This value is multiplied by the complex structure factor at the nearest\n"
"Bragg position, i.e. the gradient of the shape transform across each\n"
"appearance of the shape transform is not included, for speed of calculation.\n"
"This square modulus of this value is take, and multiplied by the Bragg\n"
"intensitiy at the neares Bragg position. That means that the gradient of\n"
"the underlying molecule transform is not included, for speed of calculation.\n"
"The Bragg intensities are taken from the file specified on the command line\n"
"with the --intensities option.\n"
"\n"
"Complex structure factors are calculated using a combination of the Henke\n"
"and Waasmeier-Kirfel scattering factors. The Henke factors are complex\n"
"and energy dependence, whereas the Waas-Kirf values are real-valued and\n"
"|q|-dependent. The difference between the Waas-Kirf value at the\n"
"appropriate |q| and the same value at |q|=0 is subtracted from the Henke\n"
"value. The Henke values are linearly interpolated from the provided tables\n"
"(note that the interpolation should really be exponential).\n"
"\n"
"The modulus of the structure factor is taken and squared. Intensity from\n"
"water is then added according to the first term of equation 5 from\n"
"Phys Chem Chem Phys 2003 (5) 1981--1991.\n"
"Intensity from water is then added according to the first term of equation\n"
"5 from Phys Chem Chem Phys 2003 (5) 1981--1991. This simulates the\n"
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
......@@ -101,14 +97,12 @@ static void show_details()
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
"where |F(q)|^2 is the value calculated as described above.\n"
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
"algorithm. When the intensity is sufficiently high that Knuth's algorithm\n"
"would result in machine precision problems, a normal distribution with\n"
"standard deviation sqrt(I) is used instead.\n"
"\n"
"The coherent, elastic part of the diffuse scattering from the water jet can\n"
"be simulated.\n"
);
}
......@@ -153,6 +147,8 @@ int main(int argc, char *argv[])
struct image image;
struct gpu_context *gctx = NULL;
double *powder;
char *intfile = NULL;
double *intensities;
int config_simdetails = 0;
int config_nearbragg = 0;
int config_randomquat = 0;
......@@ -178,7 +174,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},
{"no-sfac", 0, &config_nosfac, 1},
{"intensities", 0, NULL, 'i'},
{"powder", 0, &config_powder, 1},
{0, 0, NULL, 0}
};
......@@ -202,6 +198,11 @@ int main(int argc, char *argv[])
break;
}
case 'i' : {
intfile = strdup(optarg);
break;
}
case 0 : {
break;
}
......@@ -224,6 +225,17 @@ int main(int argc, char *argv[])
return 1;
}
if ( intfile == NULL ) {
/* Gentle reminder */
STATUS("You didn't specify the file containing the ");
STATUS("reflection intensities (with --intensities).\n");
STATUS("I'll simulate a flat intensity distribution.\n");
intensities = NULL;
} else {
intensities = read_reflections(intfile);
free(intfile);
}
/* Define image parameters */
image.width = 1024;
image.height = 1024;
......@@ -276,11 +288,11 @@ int main(int argc, char *argv[])
if ( config_gpu ) {
if ( gctx == NULL ) {
gctx = setup_gpu(config_nosfac, &image,
image.molecule);
intensities);
}
get_diffraction_gpu(gctx, &image, na, nb, nc);
} else {
get_diffraction(&image, na, nb, nc, config_nosfac,
get_diffraction(&image, na, nb, nc, intensities,
!config_nowater);
}
if ( image.molecule == NULL ) {
......
src/process_hkl.c
View file @ 6a9811f6
......@@ -37,31 +37,31 @@ static void show_help(const char *s)
printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
" -h, --help Display this help message.\n"
" -i, --input=<filename> Specify input filename (\"-\" for stdin).\n"
" -o, --output=<filename> Specify output filename for merged intensities\n"
" (don't specify for no output).\n"
" -h, --help Display this help message.\n"
" -i, --input=<filename> Specify input filename (\"-\" for stdin).\n"
" -o, --output=<filename> Specify output filename for merged intensities\n"
" (don't specify for no output).\n"
"\n"
" --max-only Take the integrated intensity to be equal to the\n"
" maximum intensity measured for that reflection.\n"
" The default is to use the mean value from all\n"
" measurements.\n"
" -e, --output-every=<n> Analyse figures of merit after every n patterns\n"
" Default: 1000. A value of zero means to do the\n"
" analysis only after reading all the patterns.\n"
" --no-analyse Don't perform any kind of analysis, just merge the\n"
" intensities.\n"
" --sum Sum (rather than average) the intensities for the\n"
" final output list. This is useful for comparing\n"
" results to radially summed powder patterns, but\n"
" will break R-factor analysis.\n"
" -r, --rvsq Output lists of R vs |q| (\"Luzzatti plots\") when\n"
" analysing figures of merit.\n"
" --stop-after=<n> Stop after processing n patterns. Zero means\n"
" keep going until the end of the input, and is the\n"
" default.\n"
" --zone-axis Output an [001] zone axis pattern each time the\n"
" figures of merit are analysed.\n");
" --max-only Take the integrated intensity to be equal to the\n"
" maximum intensity measured for that reflection.\n"
" The default is to use the mean value from all\n"
" measurements.\n"
" --sum Sum (rather than average) the intensities for the\n"
" final output list. This is useful for comparing\n"
" results to radially summed powder patterns, but\n"
" will break R-factor analysis.\n"
" --stop-after=<n> Stop after processing n patterns. Zero means\n"
" keep going until the end of the input, and is\n"
" the default.\n"
" -c, --compare-with=<file> Compare with reflection intensities in this file\n"
"\n"
" -e, --output-every=<n> Analyse figures of merit after every n patterns\n"
" Default: 1000. A value of zero means to do the\n"
" analysis only after reading all the patterns.\n"
" -r, --rvsq Output lists of R vs |q| (\"Luzzatti plots\")\n"
" when analysing figures of merit.\n"
" --zone-axis Output an [001] zone axis pattern each time the\n"
" figures of merit are analysed.\n");
}
......@@ -193,8 +193,8 @@ int main(int argc, char *argv[])
int config_rvsq = 0;
int config_stopafter = 0;
int config_zoneaxis = 0;
int config_noanalyse = 0;
int config_sum = 0;
char *intfile = NULL;
/* Long options */
const struct option longopts[] = {
......@@ -203,10 +203,10 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"max-only", 0, &config_maxonly, 1},
{"output-every", 1, NULL, 'e'},
{"no-analyse", 0, &config_noanalyse, 1},
{"rvsq", 0, NULL, 'r'},
{"stop-after", 1, NULL, 's'},
{"zone-axis", 0, &config_zoneaxis, 1},
{"compare-with", 0, NULL, 'c'},
{"sum", 0, &config_sum, 1},
{0, 0, NULL, 0}
};
......@@ -245,6 +245,11 @@ int main(int argc, char *argv[])
break;
}
case 'c' : {
intfile = strdup(optarg);
break;
}
case 0 : {
break;
}
......@@ -261,14 +266,18 @@ int main(int argc, char *argv[])
return 1;
}
if ( intfile != NULL ) {
STATUS("Comparing against '%s'\n", intfile);
trueref = read_reflections(intfile);
free(intfile);
} else {
trueref = NULL;
}
ref = new_list_intensity();
counts = new_list_count();
mol = load_molecule();
if ( !config_noanalyse ) {
get_reflections_cached(mol, eV_to_J(2.0e3));
trueref = ideal_intensities(mol->reflections);
}
if ( strcmp(filename, "-") == 0 ) {
fh = stdin;
......@@ -331,7 +340,7 @@ int main(int argc, char *argv[])
fclose(fh);
if ( !config_noanalyse ) {
if ( trueref != NULL ) {
process_reflections(ref, trueref, counts, n_patterns, mol->cell,
config_rvsq, config_zoneaxis);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use