diff options
Diffstat (limited to 'src/cell.c')
| -rw-r--r-- | src/cell.c | 1164 |
1 files changed, 0 insertions, 1164 deletions
diff --git a/src/cell.c b/src/cell.c deleted file mode 100644 index c31697bc..00000000 --- a/src/cell.c +++ /dev/null @@ -1,1164 +0,0 @@ -/* - * cell.c - * - * Unit Cell Calculations - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <math.h> -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_linalg.h> - -#include "cell.h" -#include "utils.h" -#include "image.h" - - -/** - * SECTION:unitcell - * @short_description: Unit cell - * @title: UnitCell - * @section_id: - * @see_also: - * @include: "cell.h" - * @Image: - * - * This structure represents a unit cell. - */ - - - -/* Weighting factor of lengths relative to angles */ -#define LWEIGHT (10.0e-9) - -typedef enum { - CELL_REP_CRYST, - CELL_REP_CART, - CELL_REP_RECIP -} CellRepresentation; - -struct _unitcell { - - CellRepresentation rep; - - /* Crystallographic representation */ - double a; /* m */ - double b; /* m */ - double c; /* m */ - double alpha; /* Radians */ - double beta; /* Radians */ - double gamma; /* Radians */ - - /* Cartesian representation */ - double ax; double bx; double cx; - double ay; double by; double cy; - double az; double bz; double cz; - - /* Cartesian representation of reciprocal axes */ - double axs; double bxs; double cxs; - double ays; double bys; double cys; - double azs; double bzs; double czs; - - char *pointgroup; - char *spacegroup; -}; - - -/************************** Setters and Constructors **************************/ - - -/** - * cell_new: - * - * Create a new %UnitCell. - * - * Returns: the new unit cell, or NULL on failure. - * - */ -UnitCell *cell_new() -{ - UnitCell *cell; - - cell = malloc(sizeof(UnitCell)); - if ( cell == NULL ) return NULL; - - cell->a = 1.0; - cell->b = 1.0; - cell->c = 1.0; - cell->alpha = M_PI_2; - cell->beta = M_PI_2; - cell->gamma = M_PI_2; - - cell->rep = CELL_REP_CRYST; - - cell->pointgroup = strdup("1"); - cell->spacegroup = strdup("P 1"); - - return cell; -} - - -/** - * cell_free: - * @cell: A %UnitCell to free. - * - * Frees a %UnitCell, and all internal resources concerning that cell. - * - */ -void cell_free(UnitCell *cell) -{ - if ( cell == NULL ) return; - free(cell->pointgroup); - free(cell->spacegroup); - free(cell); -} - - -void cell_set_parameters(UnitCell *cell, double a, double b, double c, - double alpha, double beta, double gamma) -{ - if ( cell == NULL ) return; - - cell->a = a; - cell->b = b; - cell->c = c; - cell->alpha = alpha; - cell->beta = beta; - cell->gamma = gamma; - - cell->rep = CELL_REP_CRYST; -} - - -void cell_set_cartesian(UnitCell *cell, - double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz) -{ - if ( cell == NULL ) return; - - cell->ax = ax; cell->ay = ay; cell->az = az; - cell->bx = bx; cell->by = by; cell->bz = bz; - cell->cx = cx; cell->cy = cy; cell->cz = cz; - - cell->rep = CELL_REP_CART; -} - - -void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az) -{ - if ( cell == NULL ) return; - cell->ax = ax; cell->ay = ay; cell->az = az; - cell->rep = CELL_REP_CART; -} - - -void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz) -{ - if ( cell == NULL ) return; - cell->bx = bx; cell->by = by; cell->bz = bz; - cell->rep = CELL_REP_CART; -} - - -void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz) -{ - if ( cell == NULL ) return; - cell->cx = cx; cell->cy = cy; cell->cz = cz; - cell->rep = CELL_REP_CART; -} - - -UnitCell *cell_new_from_parameters(double a, double b, double c, - double alpha, double beta, double gamma) -{ - UnitCell *cell; - - cell = cell_new(); - if ( cell == NULL ) return NULL; - - cell_set_parameters(cell, a, b, c, alpha, beta, gamma); - - return cell; -} - - -UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, - struct rvec cs) -{ - UnitCell *cell; - - cell = cell_new(); - if ( cell == NULL ) return NULL; - - cell->axs = as.u; cell->ays = as.v; cell->azs = as.w; - cell->bxs = bs.u; cell->bys = bs.v; cell->bzs = bs.w; - cell->cxs = cs.u; cell->cys = cs.v; cell->czs = cs.w; - - cell->rep = CELL_REP_RECIP; - - return cell; -} - - -UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c) -{ - UnitCell *cell; - - cell = cell_new(); - if ( cell == NULL ) return NULL; - - cell->ax = a.u; cell->ay = a.v; cell->az = a.w; - cell->bx = b.u; cell->by = b.v; cell->bz = b.w; - cell->cx = c.u; cell->cy = c.v; cell->cz = c.w; - - cell->rep = CELL_REP_CART; - - return cell; -} - - -UnitCell *cell_new_from_cell(UnitCell *orig) -{ - UnitCell *new; - double ax, ay, az, bx, by, bz, cx, cy, cz; - - new = cell_new(); - - cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz); - cell_set_pointgroup(new, orig->pointgroup); - cell_set_spacegroup(new, orig->spacegroup); - - return new; -} - - -void cell_set_reciprocal(UnitCell *cell, - double asx, double asy, double asz, - double bsx, double bsy, double bsz, - double csx, double csy, double csz) -{ - if ( cell == NULL ) return; - - cell->axs = asx; cell->ays = asy; cell->azs = asz; - cell->bxs = bsx; cell->bys = bsy; cell->bzs = bsz; - cell->cxs = csx; cell->cys = csy; cell->czs = csz; - - cell->rep = CELL_REP_RECIP; -} - - -void cell_set_spacegroup(UnitCell *cell, const char *sym) -{ - free(cell->spacegroup); - cell->spacegroup = strdup(sym); -} - - -void cell_set_pointgroup(UnitCell *cell, const char *sym) -{ - free(cell->pointgroup); - cell->pointgroup = strdup(sym); -} - - -/************************* Getter helper functions ****************************/ - -static int cell_crystallographic_to_cartesian(UnitCell *cell, - double *ax, double *ay, double *az, - double *bx, double *by, double *bz, - double *cx, double *cy, double *cz) -{ - double tmp, V, cosalphastar, cstar; - - /* Firstly: Get a in terms of x, y and z - * +a (cryst) is defined to lie along +x (cart) */ - *ax = cell->a; - *ay = 0.0; - *az = 0.0; - - /* b in terms of x, y and z - * b (cryst) is defined to lie in the xy (cart) plane */ - *bx = cell->b*cos(cell->gamma); - *by = cell->b*sin(cell->gamma); - *bz = 0.0; - - tmp = cos(cell->alpha)*cos(cell->alpha) - + cos(cell->beta)*cos(cell->beta) - + cos(cell->gamma)*cos(cell->gamma) - - 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma); - V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp); - - cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha); - cosalphastar /= sin(cell->beta)*sin(cell->gamma); - - cstar = (cell->a * cell->b * sin(cell->gamma))/V; - - /* c in terms of x, y and z */ - *cx = cell->c*cos(cell->beta); - *cy = -cell->c*sin(cell->beta)*cosalphastar; - *cz = 1.0/cstar; - - return 0; -} - - -/* Why yes, I do enjoy long argument lists...! */ -static int cell_invert(double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - double *asx, double *asy, double *asz, - double *bsx, double *bsy, double *bsz, - double *csx, double *csy, double *csz) -{ - int s; - gsl_matrix *m; - gsl_matrix *inv; - gsl_permutation *perm; - - m = gsl_matrix_alloc(3, 3); - if ( m == NULL ) { - ERROR("Couldn't allocate memory for matrix\n"); - return 1; - } - gsl_matrix_set(m, 0, 0, ax); - gsl_matrix_set(m, 0, 1, bx); - gsl_matrix_set(m, 0, 2, cx); - gsl_matrix_set(m, 1, 0, ay); - gsl_matrix_set(m, 1, 1, by); - gsl_matrix_set(m, 1, 2, cy); - gsl_matrix_set(m, 2, 0, az); - gsl_matrix_set(m, 2, 1, bz); - gsl_matrix_set(m, 2, 2, cz); - - /* Invert */ - perm = gsl_permutation_alloc(m->size1); - if ( perm == NULL ) { - ERROR("Couldn't allocate permutation\n"); - gsl_matrix_free(m); - return 1; - } - inv = gsl_matrix_alloc(m->size1, m->size2); - if ( inv == NULL ) { - ERROR("Couldn't allocate inverse\n"); - gsl_matrix_free(m); - gsl_permutation_free(perm); - return 1; - } - if ( gsl_linalg_LU_decomp(m, perm, &s) ) { - ERROR("Couldn't decompose matrix\n"); - gsl_matrix_free(m); - gsl_permutation_free(perm); - return 1; - } - if ( gsl_linalg_LU_invert(m, perm, inv) ) { - ERROR("Couldn't invert matrix\n"); - gsl_matrix_free(m); - gsl_permutation_free(perm); - return 1; - } - gsl_permutation_free(perm); - gsl_matrix_free(m); - - /* Transpose */ - gsl_matrix_transpose(inv); - - *asx = gsl_matrix_get(inv, 0, 0); - *bsx = gsl_matrix_get(inv, 0, 1); - *csx = gsl_matrix_get(inv, 0, 2); - *asy = gsl_matrix_get(inv, 1, 0); - *bsy = gsl_matrix_get(inv, 1, 1); - *csy = gsl_matrix_get(inv, 1, 2); - *asz = gsl_matrix_get(inv, 2, 0); - *bsz = gsl_matrix_get(inv, 2, 1); - *csz = gsl_matrix_get(inv, 2, 2); - - gsl_matrix_free(inv); - - return 0; -} - - -/********************************** Getters ***********************************/ - -int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, - double *alpha, double *beta, double *gamma) -{ - double ax, ay, az, bx, by, bz, cx, cy, cz; - - if ( cell == NULL ) return 1; - - switch ( cell->rep ) { - - case CELL_REP_CRYST: - /* Direct response */ - *a = cell->a; - *b = cell->b; - *c = cell->c; - *alpha = cell->alpha; - *beta = cell->beta; - *gamma = cell->gamma; - return 0; - - case CELL_REP_CART: - /* Convert cartesian -> crystallographic */ - *a = modulus(cell->ax, cell->ay, cell->az); - *b = modulus(cell->bx, cell->by, cell->bz); - *c = modulus(cell->cx, cell->cy, cell->cz); - - *alpha = angle_between(cell->bx, cell->by, cell->bz, - cell->cx, cell->cy, cell->cz); - *beta = angle_between(cell->ax, cell->ay, cell->az, - cell->cx, cell->cy, cell->cz); - *gamma = angle_between(cell->ax, cell->ay, cell->az, - cell->bx, cell->by, cell->bz); - return 0; - - case CELL_REP_RECIP: - /* Convert reciprocal -> crystallographic. - * Start by converting reciprocal -> cartesian */ - cell_invert(cell->axs, cell->ays, cell->azs, - cell->bxs, cell->bys, cell->bzs, - cell->cxs, cell->cys, cell->czs, - &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - /* Now convert cartesian -> crystallographic */ - *a = modulus(ax, ay, az); - *b = modulus(bx, by, bz); - *c = modulus(cx, cy, cz); - - *alpha = angle_between(bx, by, bz, cx, cy, cz); - *beta = angle_between(ax, ay, az, cx, cy, cz); - *gamma = angle_between(ax, ay, az, bx, by, bz); - return 0; - } - - return 1; -} - - -int cell_get_cartesian(UnitCell *cell, - double *ax, double *ay, double *az, - double *bx, double *by, double *bz, - double *cx, double *cy, double *cz) -{ - if ( cell == NULL ) return 1; - - switch ( cell->rep ) { - - case CELL_REP_CRYST: - /* Convert crystallographic -> cartesian. */ - return cell_crystallographic_to_cartesian(cell, - ax, ay, az, - bx, by, bz, - cx, cy, cz); - - case CELL_REP_CART: - /* Direct response */ - *ax = cell->ax; *ay = cell->ay; *az = cell->az; - *bx = cell->bx; *by = cell->by; *bz = cell->bz; - *cx = cell->cx; *cy = cell->cy; *cz = cell->cz; - return 0; - - case CELL_REP_RECIP: - /* Convert reciprocal -> cartesian */ - return cell_invert(cell->axs, cell->ays, cell->azs, - cell->bxs, cell->bys, cell->bzs, - cell->cxs, cell->cys, cell->czs, - ax, ay, az, bx, by, bz, cx, cy, cz); - - } - - return 1; -} - - -int cell_get_reciprocal(UnitCell *cell, - double *asx, double *asy, double *asz, - double *bsx, double *bsy, double *bsz, - double *csx, double *csy, double *csz) -{ - int r; - double ax, ay, az, bx, by, bz, cx, cy, cz; - if ( cell == NULL ) return 1; - - switch ( cell->rep ) { - - case CELL_REP_CRYST: - /* Convert crystallographic -> reciprocal */ - r = cell_crystallographic_to_cartesian(cell, - &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz); - if ( r ) return r; - return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz, - asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - - case CELL_REP_CART: - /* Convert cartesian -> reciprocal */ - cell_invert(cell->ax, cell->ay, cell->az, - cell->bx, cell->by, cell->bz, - cell->cx, cell->cy, cell->cz, - asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - return 0; - - case CELL_REP_RECIP: - /* Direct response */ - *asx = cell->axs; *asy = cell->ays; *asz = cell->azs; - *bsx = cell->bxs; *bsy = cell->bys; *bsz = cell->bzs; - *csx = cell->cxs; *csy = cell->cys; *csz = cell->czs; - return 0; - - } - - return 1; -} - - -const char *cell_get_pointgroup(UnitCell *cell) -{ - return cell->pointgroup; -} - - -const char *cell_get_spacegroup(UnitCell *cell) -{ - return cell->spacegroup; -} - - - - - -/********************************* Utilities **********************************/ - -static const char *cell_rep(UnitCell *cell) -{ - switch ( cell->rep ) { - case CELL_REP_CRYST: - return "crystallographic, direct space"; - case CELL_REP_CART: - return "cartesian, direct space"; - case CELL_REP_RECIP: - return "cartesian, reciprocal space"; - } - - return "unknown"; -} - - -/** - * cell_rotate: - * @in: A %UnitCell to rotate - * @quat: A %quaternion - * - * Rotate a %UnitCell using a %quaternion. - * - * Returns: a newly allocated rotated copy of @in. - * - */ -UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) -{ - struct rvec a, b, c; - struct rvec an, bn, cn; - UnitCell *out = cell_new_from_cell(in); - - cell_get_cartesian(in, &a.u, &a.v, &a.w, - &b.u, &b.v, &b.w, - &c.u, &c.v, &c.w); - - an = quat_rot(a, quat); - bn = quat_rot(b, quat); - cn = quat_rot(c, quat); - - cell_set_cartesian(out, an.u, an.v, an.w, - bn.u, bn.v, bn.w, - cn.u, cn.v, cn.w); - - return out; -} - - -void cell_print(UnitCell *cell) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double a, b, c, alpha, beta, gamma; - double ax, ay, az, bx, by, bz, cx, cy, cz; - - cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); - - STATUS(" a b c alpha beta gamma\n"); - STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n", - a*1e9, b*1e9, c*1e9, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az); - STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz); - STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz); - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - - STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", - asx, asy, asz, modulus(asx, asy, asz)); - STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", - bsx, bsy, bsz, modulus(bsx, bsy, bsz)); - STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", - csx, csy, csz, modulus(csx, csy, csz)); - - STATUS("Point group: %s\n", cell_get_pointgroup(cell)); - STATUS("Cell representation is %s.\n", cell_rep(cell)); -} - - -#define MAX_CAND (1024) - -static int right_handed(struct rvec a, struct rvec b, struct rvec c) -{ - struct rvec aCb; - double aCb_dot_c; - - /* "a" cross "b" */ - aCb.u = a.v*b.w - a.w*b.v; - aCb.v = - (a.u*b.w - a.w*b.u); - aCb.w = a.u*b.v - a.v*b.u; - - /* "a cross b" dot "c" */ - aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; - - if ( aCb_dot_c > 0.0 ) return 1; - return 0; -} - - -struct cvec { - struct rvec vec; - float na; - float nb; - float nc; - float fom; -}; - - -static int same_vector(struct cvec a, struct cvec b) -{ - if ( a.na != b.na ) return 0; - if ( a.nb != b.nb ) return 0; - if ( a.nc != b.nc ) return 0; - return 1; -} - - -/* Attempt to make 'cell' fit into 'template' somehow */ -UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, - int reduce) -{ - signed int n1l, n2l, n3l; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - int i, j; - double lengths[3]; - double angles[3]; - struct cvec *cand[3]; - UnitCell *new_cell = NULL; - float best_fom = +999999999.9; /* Large number.. */ - int ncand[3] = {0,0,0}; - signed int ilow, ihigh; - float ltl = 5.0; /* percent */ - float angtol = deg2rad(1.5); - - if ( cell_get_reciprocal(template, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell for template.\n"); - return NULL; - } - - lengths[0] = modulus(asx, asy, asz); - lengths[1] = modulus(bsx, bsy, bsz); - lengths[2] = modulus(csx, csy, csz); - - angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); - angles[1] = angle_between(asx, asy, asz, csx, csy, csz); - angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); - - cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); - - if ( cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell.\n"); - return NULL; - } - - if ( reduce ) { - ilow = -2; ihigh = 4; - } else { - ilow = 0; ihigh = 1; - } - - /* Negative values mean 1/n, positive means n, zero means zero */ - for ( n1l=ilow; n1l<=ihigh; n1l++ ) { - for ( n2l=ilow; n2l<=ihigh; n2l++ ) { - for ( n3l=ilow; n3l<=ihigh; n3l++ ) { - - float n1, n2, n3; - signed int b1, b2, b3; - - n1 = (n1l>=0) ? (n1l) : (1.0/n1l); - n2 = (n2l>=0) ? (n2l) : (1.0/n2l); - n3 = (n3l>=0) ? (n3l) : (1.0/n3l); - - if ( !reduce ) { - if ( n1l + n2l + n3l > 1 ) continue; - } - - /* 'bit' values can be +1 or -1 */ - for ( b1=-1; b1<=1; b1+=2 ) { - for ( b2=-1; b2<=1; b2+=2 ) { - for ( b3=-1; b3<=1; b3+=2 ) { - - double tx, ty, tz; - double tlen; - int i; - - n1 *= b1; n2 *= b2; n3 *= b3; - - tx = n1*asx + n2*bsx + n3*csx; - ty = n1*asy + n2*bsy + n3*csy; - tz = n1*asz + n2*bsz + n3*csz; - tlen = modulus(tx, ty, tz); - - /* Test modulus for agreement with moduli of template */ - for ( i=0; i<3; i++ ) { - - if ( !within_tolerance(lengths[i], tlen, ltl) ) - continue; - - cand[i][ncand[i]].vec.u = tx; - cand[i][ncand[i]].vec.v = ty; - cand[i][ncand[i]].vec.w = tz; - cand[i][ncand[i]].na = n1; - cand[i][ncand[i]].nb = n2; - cand[i][ncand[i]].nc = n3; - cand[i][ncand[i]].fom = fabs(lengths[i] - tlen); - if ( ncand[i] == MAX_CAND ) { - ERROR("Too many candidates\n"); - } else { - ncand[i]++; - } - } - - } - } - } - - } - } - } - - if ( verbose ) { - STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); - } - - for ( i=0; i<ncand[0]; i++ ) { - for ( j=0; j<ncand[1]; j++ ) { - - double ang; - int k; - float fom1; - - if ( same_vector(cand[0][i], cand[1][j]) ) continue; - - /* Measure the angle between the ith candidate for axis 0 - * and the jth candidate for axis 1 */ - ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v, - cand[0][i].vec.w, cand[1][j].vec.u, - cand[1][j].vec.v, cand[1][j].vec.w); - - /* Angle between axes 0 and 1 should be angle 2 */ - if ( fabs(ang - angles[2]) > angtol ) continue; - - fom1 = fabs(ang - angles[2]); - - for ( k=0; k<ncand[2]; k++ ) { - - float fom2, fom3; - - if ( same_vector(cand[1][j], cand[2][k]) ) continue; - - /* Measure the angle between the current candidate for - * axis 0 and the kth candidate for axis 2 */ - ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v, - cand[0][i].vec.w, cand[2][k].vec.u, - cand[2][k].vec.v, cand[2][k].vec.w); - - /* ... it should be angle 1 ... */ - if ( fabs(ang - angles[1]) > angtol ) continue; - - fom2 = fom1 + fabs(ang - angles[1]); - - /* Finally, the angle between the current candidate for - * axis 1 and the kth candidate for axis 2 */ - ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, - cand[1][j].vec.w, cand[2][k].vec.u, - cand[2][k].vec.v, cand[2][k].vec.w); - - /* ... it should be angle 0 ... */ - if ( fabs(ang - angles[0]) > angtol ) continue; - - /* Unit cell must be right-handed */ - if ( !right_handed(cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec) ) continue; - - fom3 = fom2 + fabs(ang - angles[0]); - fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom - + cand[2][k].fom); - - if ( fom3 < best_fom ) { - if ( new_cell != NULL ) free(new_cell); - new_cell = cell_new_from_reciprocal_axes( - cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec); - best_fom = fom3; - } - - } - - } - } - - free(cand[0]); - free(cand[1]); - free(cand[2]); - - return new_cell; -} - - - -UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int i; - double lengths[3]; - int used[3]; - struct rvec real_a, real_b, real_c; - struct rvec params[3]; - double alen, blen, clen; - float ltl = 5.0; /* percent */ - int have_real_a; - int have_real_b; - int have_real_c; - - /* Get the lengths to match */ - if ( cell_get_cartesian(template, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell for template.\n"); - return NULL; - } - alen = modulus(ax, ay, az); - blen = modulus(bx, by, bz); - clen = modulus(cx, cy, cz); - - /* Get the lengths from the cell and turn them into anonymous vectors */ - if ( cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell.\n"); - return NULL; - } - lengths[0] = modulus(ax, ay, az); - lengths[1] = modulus(bx, by, bz); - lengths[2] = modulus(cx, cy, cz); - used[0] = 0; used[1] = 0; used[2] = 0; - params[0].u = ax; params[0].v = ay; params[0].w = az; - params[1].u = bx; params[1].v = by; params[1].w = bz; - params[2].u = cx; params[2].v = cy; params[2].w = cz; - - real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; - real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; - real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; - - /* Check each vector against a and b */ - have_real_a = 0; - have_real_b = 0; - for ( i=0; i<3; i++ ) { - if ( within_tolerance(lengths[i], alen, ltl) - && !used[i] && !have_real_a ) - { - used[i] = 1; - memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); - have_real_a = 1; - } - if ( within_tolerance(lengths[i], blen, ltl) - && !used[i] && !have_real_b ) - { - used[i] = 1; - memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); - have_real_b = 1; - } - } - - /* Have we matched both a and b? */ - if ( !(have_real_a && have_real_b) ) return NULL; - - /* "c" is "the other one" */ - have_real_c = 0; - for ( i=0; i<3; i++ ) { - if ( !used[i] ) { - memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); - have_real_c = 1; - } - } - - if ( !have_real_c ) { - ERROR("Huh? Couldn't find the third vector.\n"); - ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); - return NULL; - } - - /* Flip c if not right-handed */ - if ( !right_handed(real_a, real_b, real_c) ) { - real_c.u = -real_c.u; - real_c.v = -real_c.v; - real_c.w = -real_c.w; - } - - return cell_new_from_direct_axes(real_a, real_b, real_c); -} - - -/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ -double resolution(UnitCell *cell, signed int h, signed int k, signed int l) -{ - double a, b, c, alpha, beta, gamma; - - cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); - - const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha) - - cos(beta)*cos(beta) - - cos(gamma)*cos(gamma) - + 2*cos(alpha)*cos(beta)*cos(gamma) ); - - const double S11 = b*b*c*c*sin(alpha)*sin(alpha); - const double S22 = a*a*c*c*sin(beta)*sin(beta); - const double S33 = a*a*b*b*sin(gamma)*sin(gamma); - const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma)); - const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha)); - const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta)); - - const double brackets = S11*h*h + S22*k*k + S33*l*l - + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l; - const double oneoverdsq = brackets / Vsq; - const double oneoverd = sqrt(oneoverdsq); - - return oneoverd / 2; -} - - -static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym) -{ - char *new = NULL; - - if ( strcmp(sym, "P 1") == 0 ) new = "1"; - if ( strcmp(sym, "P 63") == 0 ) new = "6"; - if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222"; - if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222"; - if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422"; - - if ( new != NULL ) { - cell_set_pointgroup(cell, new); - } else { - ERROR("Can't determine point group for '%s'\n", sym); - } -} - - -UnitCell *load_cell_from_pdb(const char *filename) -{ - FILE *fh; - char *rval; - UnitCell *cell = NULL; - - fh = fopen(filename, "r"); - if ( fh == NULL ) { - ERROR("Couldn't open '%s'\n", filename); - return NULL; - } - - do { - - char line[1024]; - - rval = fgets(line, 1023, fh); - - if ( strncmp(line, "CRYST1", 6) == 0 ) { - - float a, b, c, al, be, ga; - char as[10], bs[10], cs[10]; - char als[8], bes[8], gas[8]; - char *sym; - int r; - - memcpy(as, line+6, 9); as[9] = '\0'; - memcpy(bs, line+15, 9); bs[9] = '\0'; - memcpy(cs, line+24, 9); cs[9] = '\0'; - memcpy(als, line+33, 7); als[7] = '\0'; - memcpy(bes, line+40, 7); bes[7] = '\0'; - memcpy(gas, line+47, 7); gas[7] = '\0'; - - r = sscanf(as, "%f", &a); - r += sscanf(bs, "%f", &b); - r += sscanf(cs, "%f", &c); - r += sscanf(als, "%f", &al); - r += sscanf(bes, "%f", &be); - r += sscanf(gas, "%f", &ga); - - if ( r != 6 ) { - STATUS("Couldn't understand CRYST1 line.\n"); - continue; - } - - cell = cell_new_from_parameters(a*1e-10, - b*1e-10, c*1e-10, - deg2rad(al), - deg2rad(be), - deg2rad(ga)); - - if ( strlen(line) > 65 ) { - sym = strndup(line+55, 10); - notrail(sym); - cell_set_pointgroup_from_pdb(cell, sym); - cell_set_spacegroup(cell, sym); - free(sym); - } else { - cell_set_pointgroup_from_pdb(cell, "P 1"); - cell_set_spacegroup(cell, "P 1"); - ERROR("CRYST1 line without space group.\n"); - } - - break; /* Done */ - } - - } while ( rval != NULL ); - - fclose(fh); - - return cell; -} - - -#ifdef GSL_FUDGE -/* Force the linker to bring in CBLAS to make GSL happy */ -void cell_fudge_gslcblas() -{ - STATUS("%p\n", cblas_sgemm); -} -#endif - - -UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot) -{ - UnitCell *out; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double xnew, ynew, znew; - - cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, - &bsz, &csx, &csy, &csz); - - /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */ - xnew = asx*cos(omega) + asy*sin(omega); - ynew = -asx*sin(omega) + asy*cos(omega); - znew = asz; - asx = xnew; asy = ynew; asz = znew; - xnew = bsx*cos(omega) + bsy*sin(omega); - ynew = -bsx*sin(omega) + bsy*cos(omega); - znew = bsz; - bsx = xnew; bsy = ynew; bsz = znew; - xnew = csx*cos(omega) + csy*sin(omega); - ynew = -csx*sin(omega) + csy*cos(omega); - znew = csz; - csx = xnew; csy = ynew; csz = znew; - - /* Rotate by "phi" about +x (not parallel to anything specific) */ - xnew = asx; - ynew = asy*cos(phi) + asz*sin(phi); - znew = -asy*sin(phi) + asz*cos(phi); - asx = xnew; asy = ynew; asz = znew; - xnew = bsx; - ynew = bsy*cos(phi) + bsz*sin(phi); - znew = -bsy*sin(phi) + bsz*cos(phi); - bsx = xnew; bsy = ynew; bsz = znew; - xnew = csx; - ynew = csy*cos(phi) + csz*sin(phi); - znew = -csy*sin(phi) + csz*cos(phi); - csx = xnew; csy = ynew; csz = znew; - - /* Rotate by "rot" about the new +z (in-plane rotation) */ - xnew = asx*cos(rot) + asy*sin(rot); - ynew = -asx*sin(rot) + asy*cos(rot); - znew = asz; - asx = xnew; asy = ynew; asz = znew; - xnew = bsx*cos(rot) + bsy*sin(rot); - ynew = -bsx*sin(rot) + bsy*cos(rot); - znew = bsz; - bsx = xnew; bsy = ynew; bsz = znew; - xnew = csx*cos(rot) + csy*sin(rot); - ynew = -csx*sin(rot) + csy*cos(rot); - znew = csz; - csx = xnew; csy = ynew; csz = znew; - - out = cell_new_from_cell(in); - cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - - return out; -} - - -int cell_is_sensible(UnitCell *cell) -{ - double a, b, c, al, be, ga; - - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - if ( al + be + ga >= 2.0*M_PI ) return 0; - if ( al + be - ga >= 2.0*M_PI ) return 0; - if ( al - be + ga >= 2.0*M_PI ) return 0; - if ( - al + be + ga >= 2.0*M_PI ) return 0; - if ( al + be + ga <= 0.0 ) return 0; - if ( al + be - ga <= 0.0 ) return 0; - if ( al - be + ga <= 0.0 ) return 0; - if ( - al + be + ga <= 0.0 ) return 0; - if ( isnan(al) ) return 0; - if ( isnan(be) ) return 0; - if ( isnan(ga) ) return 0; - return 1; -} |