diff options
-rw-r--r-- | .gitignore | 2 | ||||
-rw-r--r-- | Makefile | 16 | ||||
-rw-r--r-- | crystal.c | 110 | ||||
-rw-r--r-- | crystal.h | 33 | ||||
-rw-r--r-- | main.c | 334 | ||||
-rw-r--r-- | triclinator.run | 32 | ||||
-rw-r--r-- | util.c | 39 | ||||
-rw-r--r-- | util.h | 32 |
8 files changed, 598 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6003c91 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.o +triclinator diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..3878741 --- /dev/null +++ b/Makefile @@ -0,0 +1,16 @@ +triclinator: main.o crystal.o util.o + gcc main.o crystal.o util.o -o triclinator -lgsl -lgslcblas + +main.o: main.c + gcc -g -c -I/usr/include main.c -o main.o + +crystal.o: crystal.c + gcc -g -c -I/usr/include crystal.c -o crystal.o + +util.o: util.c + gcc -g -c -I/usr/include util.c -o util.o + +clean: + rm -f triclinator main.o crystal.o util.o + +.PHONY: clean diff --git a/crystal.c b/crystal.c new file mode 100644 index 0000000..065236a --- /dev/null +++ b/crystal.c @@ -0,0 +1,110 @@ +/* + * crystal.c + * + * Crystallographic stuff + * + * (c) 2009 Thomas White <taw27@cam.ac.uk> + * + * Triclinator - solve nasty triclinic unit cells + * + */ + + +#include <math.h> + +#include "crystal.h" +#include "util.h" + + +/* Some attempt at preserving sanity */ +#define S11 s11(cell) +#define S22 s22(cell) +#define S33 s33(cell) +#define S12 s12(cell) +#define S23 s23(cell) +#define S13 s13(cell) + + +static double s11(Cell cell) +{ + return cell.b*cell.b*cell.c*cell.c*sin(cell.al)*sin(cell.al); +} + +static double s22(Cell cell) +{ + return cell.a*cell.a*cell.c*cell.c*sin(cell.be)*sin(cell.be); +} + +static double s33(Cell cell) +{ + return cell.a*cell.a*cell.b*cell.b*sin(cell.ga)*sin(cell.ga); +} + +static double s12(Cell cell) +{ + return cell.a*cell.b*cell.c*cell.c*(cos(cell.al)*cos(cell.be) - cos(cell.ga)); +} + +static double s23(Cell cell) +{ + return cell.a*cell.a*cell.b*cell.c*(cos(cell.be)*cos(cell.ga) - cos(cell.al)); +} + +static double s13(Cell cell) +{ + return cell.a*cell.b*cell.b*cell.c*(cos(cell.ga)*cos(cell.al) - cos(cell.be)); +} + + +static double volume(Cell cell) +{ + return cell.a*cell.b*cell.c*sqrt(1 - cos(cell.al)*cos(cell.al) + - cos(cell.be)*cos(cell.be) + - cos(cell.ga)*cos(cell.ga) + + 2*cos(cell.al)*cos(cell.be)*cos(cell.ga) ); +} + + +static double dspacing(Cell cell, double h, double k, double l) +{ + double sum_S_terms, one_over_Vsq, one_over_LHS; + + sum_S_terms = S11*h*h + S22*k*k + S33*l*l + + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l; + + one_over_Vsq = 1 / pow(volume(cell), 2); + + one_over_LHS = 1/( one_over_Vsq*sum_S_terms ); + + return sqrt(one_over_LHS); +} + + +static double plane_angle(Cell cell, double h1, double k1, double l1, + double h2, double k2, double l2) +{ + double sum_S_terms, d1, d2, dd_over_Vsq; + + sum_S_terms = S11*h1*h2 + S22*k1*k2 * S33*l1*l2 + + (S23*(k1*l2 + k2*l1)) + + (S13*(l1*h2 + l2*h1)) + + (S12*(h1*k2 + h2*k1)); + + d1 = dspacing(cell, h1, k1, l1); + d2 = dspacing(cell, h2, k2, l2); + dd_over_Vsq = d1*d2 / pow(volume(cell), 2); + + return acos(dd_over_Vsq * sum_S_terms); +} + + +/* Return what the measurement 'val' would have been if the cell were 'cell' */ +double crystal_calc(MVal val, Cell cell) +{ + if ( is_dspacing(val) ) { + return dspacing(cell, val.h1, val.k1, val.l1); + } else { + return plane_angle(cell, val.h1, val.k1, val.l1, + val.h2, val.k2, val.l2); + } +} diff --git a/crystal.h b/crystal.h new file mode 100644 index 0000000..5c38da8 --- /dev/null +++ b/crystal.h @@ -0,0 +1,33 @@ +/* + * crystal.h + * + * Crystallographic stuff + * + * (c) 2009 Thomas White <taw27@cam.ac.uk> + * + * Triclinator - solve nasty triclinic unit cells + * + */ + +#ifndef CRYSTAL_H +#define CRYSTAL_H + +#define RAD2DEG(a) ((a)*180/M_PI) +#define DEG2RAD(a) ((a)*M_PI/180) + +#include "util.h" + +typedef struct +{ + double a; + double b; + double c; + double al; + double be; + double ga; +} Cell; + +/* Return what the measurement 'val' would have been if the cell were 'cell' */ +extern double crystal_calc(MVal val, Cell cell); + +#endif /* CRYSTAL_H */ @@ -0,0 +1,334 @@ +/* + * main.c + * + * The Top Level Source File + * + * (c) 2009 Thomas White <taw27@cam.ac.uk> + * + * Triclinator - solve nasty triclinic unit cells + * + */ + + +#include <stdlib.h> +#include <stdio.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_randist.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_multifit_nlin.h> +#include <math.h> + +#include "crystal.h" +#include "util.h" + +#define DSTEP 0.0001 +#define ASTEP 0.00001 + +struct data { + size_t n; /* Number of measurements */ + MVal *vals; /* The measurements themselves */ +}; + + +int expb_f(const gsl_vector *x, void *data, gsl_vector *f) +{ + size_t n = ((struct data *)data)->n; + MVal *vals = ((struct data *) data)->vals; + Cell cell; + size_t i; + + cell.a = gsl_vector_get(x, 0); + cell.b = gsl_vector_get(x, 1); + cell.c = gsl_vector_get(x, 2); + cell.al = gsl_vector_get(x, 3); + cell.be = gsl_vector_get(x, 4); + cell.ga = gsl_vector_get(x, 5); + + for ( i=0; i<n; i++ ) { + double calc = crystal_calc(vals[i], cell); + gsl_vector_set(f, i, (calc - vals[i].meas)/vals[i].esd); + } + + return GSL_SUCCESS; +} + + +int expb_df(const gsl_vector *x, void *data, gsl_matrix *J) +{ + size_t n = ((struct data *)data)->n; + MVal *vals = ((struct data *) data)->vals; + Cell cell; + size_t i; + + cell.a = gsl_vector_get(x, 0); + cell.b = gsl_vector_get(x, 1); + cell.c = gsl_vector_get(x, 2); + cell.al = gsl_vector_get(x, 3); + cell.be = gsl_vector_get(x, 4); + cell.ga = gsl_vector_get(x, 5); + + for ( i=0; i<n; i++ ) { + + double calc, calc2, s; + Cell cell2; + MVal val; + + val = vals[i]; + + calc = crystal_calc(vals[i], cell); + s = vals[i].esd; + + cell2 = cell; cell2.a += DSTEP; calc2 = crystal_calc(val, cell2); + gsl_matrix_set(J, i, 0, (calc2 - calc)/(DSTEP*s)); /* df/da */ + + cell2 = cell; cell2.b += DSTEP; calc2 = crystal_calc(val, cell2); + gsl_matrix_set(J, i, 1, (calc2 - calc)/(DSTEP*s)); /* df/db */ + + cell2 = cell; cell2.c += DSTEP; calc2 = crystal_calc(val, cell2); + gsl_matrix_set(J, i, 2, (calc2 - calc)/(DSTEP*s)); /* df/dc */ + + cell2 = cell; cell2.al += ASTEP; calc2 = crystal_calc(val, cell2); + gsl_matrix_set(J, i, 3, (calc2 - calc)/(ASTEP*s)); /* df/dal */ + + cell2 = cell; cell2.be += ASTEP; calc2 = crystal_calc(val, cell2); + gsl_matrix_set(J, i, 4, (calc2 - calc)/(ASTEP*s)); /* df/dbe */ + + cell2 = cell; cell2.ga += ASTEP; calc2 = crystal_calc(val, cell2); + gsl_matrix_set(J, i, 5, (calc2 - calc)/(ASTEP*s)); /* df/dga */ + + } + + return GSL_SUCCESS; +} + + +int expb_fdf(const gsl_vector *x, void *data, gsl_vector *f, gsl_matrix *J) +{ + expb_f(x, data, f); + expb_df(x, data, J); + return GSL_SUCCESS; +} + + +void print_state (size_t iter, gsl_multifit_fdfsolver *s) +{ + printf ("%3u %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f", iter, + gsl_vector_get(s->x, 0), + gsl_vector_get(s->x, 1), + gsl_vector_get(s->x, 2), + gsl_vector_get(s->x, 3), + gsl_vector_get(s->x, 4), + gsl_vector_get(s->x, 5)); + +} + + +int main(void) +{ + const gsl_multifit_fdfsolver_type *T; + gsl_multifit_fdfsolver *s; + int status; + unsigned int i, iter, done; + size_t n; + gsl_matrix *covar; + MVal *vals; + struct data d; + gsl_multifit_function_fdf f; + double x_init[6]; + gsl_vector_view x = gsl_vector_view_array(x_init, 6); + Cell cell; + + /* Get the data to be fitted */ + done = 0; n = 0; + vals = NULL; + while ( !done ) { + + signed int h, k, l; + char buf[64]; + int dspacing = 0; + +// printf("Enter the indicies of the first set of planes, " +// " in the format 'h k l'.\n"); +// if ( n != 0 ) printf("If you have no further planes, " +// "just hit enter.\n"); + printf(" First plane indicies: "); + if ( fgets(buf, 63, stdin) != buf ) { + fprintf(stderr, "Error reading from input\n"); + } + if ( sscanf(buf, "%i %i %i", &h, &k, &l) != 3 ) { + if ( buf[0] == '\n' ) { + done = 1; + continue; + } else { + printf("Invalid input, try again.\n"); + } + } + + /* Allocate space for the new reading, + now that we know it exists */ + vals = realloc(vals, (n+1)*sizeof(MVal)); + vals[n].h1 = h; + vals[n].k1 = k; + vals[n].l1 = l; + +// printf("Enter the indicies of the second set of planes in the " +// "same format.\n"); +// if ( n != 0 ) printf("If you are trying to give a d-spacing " +// "only, just hit enter here.\n"); + printf(" Second plane indicies: "); + if ( fgets(buf, 63, stdin) != buf ) { + fprintf(stderr, "Error reading from input\n"); + } + if ( sscanf(buf, "%i %i %i", &h, &k, &l) != 3 ) { + if ( buf[0] == '\n' ) { + dspacing = 1; + } else { + printf("Invalid input, try again.\n"); + } + } + + if ( dspacing ) { + vals[n].h2 = 0; + vals[n].k2 = 0; + vals[n].l2 = 0; + vals[n].meas = read_value(" D-spacing: "); + } else { + vals[n].h2 = h; + vals[n].k2 = k; + vals[n].l2 = l; + vals[n].meas = DEG2RAD(read_value(" Angle between planes: ")); + } + + vals[n].esd = read_value("Estimated standard deviation: "); + + printf("--------------------------------------\n"); + n++; + + } + + x_init[0] = read_value(" Initial guess for a: "); + x_init[1] = read_value(" Initial guess for b: "); + x_init[2] = read_value(" Initial guess for c: "); + x_init[3] = read_value("Initial guess for alpha: "); + x_init[4] = read_value(" Initial guess for beta: "); + x_init[5] = read_value("Initial guess for gamma: "); + + printf("-----------------------------------------\n\n"); + + printf("Inputted values:\n"); + printf("--------------------------------------------------------------\n"); + printf(" h1 k1 l1 h2 k2 l2 Spacing Angle ESD\n"); + printf("--------------------------------------------------------------\n"); + for ( i=0; i<n; i++ ) { + printf("%3i %3i %3i ", vals[i].h1, vals[i].k1, vals[i].l1); + if ( is_dspacing(vals[i]) ) { + printf(" - - - %8.5f - ", vals[i].meas); + } else { + printf("%3i %3i %3i - %8.5f", vals[i].h2, + vals[i].k2, vals[i].l2, RAD2DEG(vals[i].meas)); + } + printf(" +/- %8.5f", vals[i].esd); + if ( !is_dspacing(vals[i]) ) { + printf(" deg\n"); + } else { + printf(" nm\n"); + } + } + printf("--------------------------------------------------------------\n"); + + d.vals = vals; + d.n = n; + + f.f = &expb_f; + f.df = &expb_df; + f.fdf = &expb_fdf; + f.n = n; + f.p = 6; + f.params = &d; + + T = gsl_multifit_fdfsolver_lmsder; + s = gsl_multifit_fdfsolver_alloc(T, n, 6); + gsl_multifit_fdfsolver_set(s, &f, &x.vector); + + iter = 0; + printf("\n"); + printf("iter a b c alpha beta gamma\n"); + printf("---------------------------------------------------------\n"); + print_state(iter, s); + + do { + + iter++; + status = gsl_multifit_fdfsolver_iterate(s); + + printf(" status = %s\n", gsl_strerror(status)); + print_state(iter, s); + if ( status ) break; + + status = gsl_multifit_test_delta(s->dx, s->x, 1e-6, 1e-6); + + } while ( (status == GSL_CONTINUE) && (iter < 500) ); + + covar = gsl_matrix_alloc(6, 6); + gsl_multifit_covar(s->J, 0.0, covar); + + #define FIT(i) gsl_vector_get(s->x, i) + #define ERR(i) sqrt(gsl_matrix_get(covar,i,i)) + + printf("\n---------------------------------------------------------\n"); + + { + double chi = gsl_blas_dnrm2(s->f); + double dof = n - 6; + double c = GSL_MAX_DBL(1, chi / sqrt(dof)); + + printf("chisq/dof = %g / %g = %g\n", pow(chi, 2.0), dof, + pow(chi, 2.0) / dof); + + printf ("a = %8.5f +/- %8.5f\n", FIT(0), c*ERR(0)); + printf ("b = %8.5f +/- %8.5f\n", FIT(1), c*ERR(1)); + printf ("c = %8.5f +/- %8.5f\n", FIT(2), c*ERR(2)); + printf ("alpha = %8.5f +/- %8.5f\n", FIT(3), c*ERR(3)); + printf (" beta = %8.5f +/- %8.5f\n", FIT(4), c*ERR(4)); + printf ("gamma = %8.5f +/- %8.5f\n", FIT(5), c*ERR(5)); + } + + printf ("status = %s\n", gsl_strerror (status)); + + gsl_multifit_fdfsolver_free(s); + gsl_matrix_free(covar); + + cell.a = FIT(0); + cell.b = FIT(1); + cell.c = FIT(2); + cell.al = FIT(3); + cell.be = FIT(4); + cell.ga = FIT(5); + + printf("\nValues given by fitted cell:\n"); + printf("--------------------------------------------------------------\n"); + printf(" h1 k1 l1 h2 k2 l2 Spacing Angle ESD\n"); + printf("--------------------------------------------------------------\n"); + for ( i=0; i<n; i++ ) { + printf("%3i %3i %3i ", vals[i].h1, vals[i].k1, vals[i].l1); + if ( is_dspacing(vals[i]) ) { + printf(" - - - %8.5f - ", + crystal_calc(vals[i], cell)); + } else { + printf("%3i %3i %3i - %8.5f", vals[i].h2, + vals[i].k2, vals[i].l2, + RAD2DEG(crystal_calc(vals[i], cell))); + } + printf(" +/- %8.5f", vals[i].esd); + if ( !is_dspacing(vals[i]) ) { + printf(" deg\n"); + } else { + printf(" nm\n"); + } + } + printf("--------------------------------------------------------------\n"); + + + return 0; +} diff --git a/triclinator.run b/triclinator.run new file mode 100644 index 0000000..2a6d503 --- /dev/null +++ b/triclinator.run @@ -0,0 +1,32 @@ +1 0 0 + +0.799 +0.04 +0 1 0 + +0.548 +0.03 +0 0 1 + +0.475 +0.02 +0 1 1 + +0.376 +0.02 +1 0 0 +0 1 0 +90.5 +0.3 +1 0 0 +0 0 1 +84.2 +0.3 + +0.821 +0.493 +0.585 +90.0 +94.7 +90.0 + @@ -0,0 +1,39 @@ +/* + * util.c + * + * Utility stuff + * + * (c) 2009 Thomas White <taw27@cam.ac.uk> + * + * Triclinator - solve nasty triclinic unit cells + * + */ + +#include <stdio.h> + +#include "util.h" + +double read_value(const char *text) +{ + while ( 1 ) { + + float d; + char buf[64]; + + printf("%s", text); + if ( fgets(buf, 63, stdin) != buf ) { + fprintf(stderr, "Error reading from input\n"); + } + if ( sscanf(buf, "%f", &d) != 1 ) { + printf("Invalid input, try again.\n"); + } else { + return d; + } + } +} + +int is_dspacing(MVal val) +{ + if ( (val.h2 == 0) && (val.k2 == 0) && (val.l2 == 0) ) return 1; + return 0; +} @@ -0,0 +1,32 @@ +/* + * util.h + * + * Utility stuff + * + * (c) 2009 Thomas White <taw27@cam.ac.uk> + * + * Triclinator - solve nasty triclinic unit cells + * + */ + +#ifndef UTIL_H +#define UTIL_H + +typedef struct +{ + signed int h1; + signed int k1; + signed int l1; + + signed int h2; + signed int k2; + signed int l2; /* Or zero */ + + double meas; + double esd; +} MVal; + +extern double read_value(const char *text); +extern int is_dspacing(MVal val); + +#endif /* UTIL_H */ |