From 193c900d095203fb74671abc1e0b52f0b1d6d62a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 16 Mar 2011 14:17:40 +0100 Subject: Remove all remaining PDB rendering stuff --- src/check_hkl.c | 1 - src/compare_hkl.c | 1 - src/diffraction-gpu.c | 1 - src/diffraction.c | 1 - src/dirax.c | 1 - src/get_hkl.c | 105 +++------- src/index.c | 1 - src/indexamajig.c | 1 - src/mosflm.c | 1 - src/pattern_sim.c | 1 - src/process_hkl.c | 1 - src/sfac.c | 542 -------------------------------------------------- src/sfac.h | 64 ------ 13 files changed, 30 insertions(+), 691 deletions(-) delete mode 100644 src/sfac.c delete mode 100644 src/sfac.h (limited to 'src') 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 #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 #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 + * (c) 2006-2011 Thomas White * * Part of CrystFEL - crystallography with a FEL * @@ -22,7 +22,6 @@ #include #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= Only include reflections mentioned in file.\n" +" -i, --input= Read reflections from .\n" +" -y, --symmetry= 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= The symmetry of the input file (-i).\n" +"\n" +"To calculate Poisson samples accurately, you must also give:\n" +" -b, --beam= Get beam parameters from file.\n" +"\n" +"You can artificially 'twin' the reflections, or expand them out:\n" " -w, --twin= Generate twinned data according to the given\n" " point group.\n" " -e, --expand= Expand reflections to this point group.\n" -" -o, --output= Output filename (default: stdout).\n" -" -i, --intensities= 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= 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= 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= Get beam parameters from file (used for sigmas).\n" -" --max-res= Calculate structure factors out to d= nm.\n" +"\n" +"Don't forget to specify the output filename:\n" +" -o, --output= 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 - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#include -#include -#include -#include -#include -#include - -#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; in_species; i++ ) { - - struct mol_species *spec; - int j; - - spec = mol->species[i]; - - for ( j=0; jn_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; in_species; i++ ) { - - struct mol_species *spec; - int j; - - spec = mol->species[i]; - - for ( j=0; jn_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; jn_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; in_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; in_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; in_species; i++ ) { - - double complex sfac; - double complex contrib = 0.0; - struct mol_species *spec; - int j; - - spec = mol->species[i]; - - for ( j=0; jn_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 - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#ifndef SFAC_H -#define SFAC_H - -#include - -#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 */ -- cgit v1.2.3