From 189da15810deabd739d7c11c6e95fea55739fe60 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 1 Aug 2020 15:13:49 +0200 Subject: Initial import from archive --- src/model.c | 523 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 523 insertions(+) create mode 100644 src/model.c (limited to 'src/model.c') diff --git a/src/model.c b/src/model.c new file mode 100644 index 0000000..08c1bb2 --- /dev/null +++ b/src/model.c @@ -0,0 +1,523 @@ +/* + * model.c + * + * Atomic Model Structures + * + * (c) 2006-2007 Thomas White + * + * synth2d - Two-Dimensional Crystallographic Fourier Synthesis + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "data.h" +#include "elements.h" +#include "main.h" +#include "statistics.h" +#include "symmetry.h" +#include "model.h" +#include "displaywindow.h" +#include "model-editor.h" +#include "multislice.h" +#include "dpsynth.h" +#include "normalise.h" + +/* This is the SPOT for the "current atomic model", but don't assume it's the only model around... */ +AtomicModel *model_current = NULL; + +/* --------------------------------------------------- Loading and saving of models --------------------------------------------------- */ + +/* NB This doesn't exactly replicate Perl's chomp()... */ +static void chomp(char *line) { + + size_t i; + + for ( i=0; iatoms[model->n_atoms].x = x; + model->atoms[model->n_atoms].y = y; + model->atoms[model->n_atoms].z = z; + model->atoms[model->n_atoms].B = B; + model->atoms[model->n_atoms].occ = occ; + model->atoms[model->n_atoms].ref = elements_lookup(tmp); + model->atoms[model->n_atoms].active = (active == 'A')?1:0; /* Anything other than "A" means "inactive" */ + model->atoms[model->n_atoms].refine = (refine == 'R')?1:0; /* Anything other than "R" means "don't refine" */ + model->n_atoms++; + } else if ( strncmp(line, "symmetry ", 9) == 0 ) { + model->sym = symmetry_encode(line+9); + } else if ( sscanf(line, "thickness %f nm", &thickness) == 1 ) { + model->thickness = thickness; + } else { + fprintf(stderr, "Unrecognised line in model file (line %i)\n", i); + } + + i++; + + } + fclose(fh); + + return model; + +} + +void model_load_as_current(const char *filename) { + + AtomicModel *new; + double new_thickness; + Symmetry new_sym; + + new = model_load(filename); + new_thickness = new->thickness; + new_sym = new->sym; + + model_move(model_current, new); /* Preserves lock_count, editor, refinementwindow, thickness and sym */ + model_free(new); + + /* Restore thickness and symmetry */ + model_current->thickness = new_thickness; + model_current->sym = new_sym; + + model_notify_update(model_current); + +} + +void model_save(const char *filename, AtomicModel *model) { + + FILE *fh; + size_t i; + + fh = fopen(filename, "w"); + fprintf(fh, "symmetry %s\n", symmetry_decode(model->sym)); + fprintf(fh, "thickness %.2f nm\n", model->thickness); + for ( i=0; in_atoms; i++ ) { + + fprintf(fh, "%f %f %f %s %f %f %c %c\n", model->atoms[i].x, model->atoms[i].y, model->atoms[i].z, + elements[model->atoms[i].ref].element_name, model->atoms[i].B, + model->atoms[i].occ, + model->atoms[i].active?'A':'a', model->atoms[i].refine?'R':'r'); + + } + + fclose(fh); + +} + +/* --------------------------------------------------- Basic model handling functions ------------------------------------------------- */ + +AtomicModel *model_new() { + + AtomicModel *model; + + model = malloc(sizeof(AtomicModel)); + + model->n_atoms = 0; + model->sym = PLANEGROUP_P1; + model->thickness = 0.0; + model->point_atoms = 0; + + model->editor = NULL; + model->refine_window = NULL; + model->lock_count = 0; + + return model; + +} + +AtomicModel *model_copy(const AtomicModel *a) { + AtomicModel *model; + model = malloc(sizeof(AtomicModel)); + memcpy(model, a, sizeof(AtomicModel)); + return model; +} + +void model_move(AtomicModel *to, AtomicModel *from) { + + ModelEditor *old_editor; + RefinementWindow *old_refine; + int old_lock_count; + Symmetry old_sym; + double old_thickness; + + old_editor = to->editor; + old_refine = to->refine_window; + old_lock_count = to->lock_count; + old_sym = to->sym; + old_thickness = to->thickness; + + memcpy(to, from, sizeof(AtomicModel)); + + to->editor = old_editor; + to->refine_window = old_refine; + to->lock_count = old_lock_count; + to->sym = old_sym; + to->thickness = old_thickness; + +} + +void model_free(AtomicModel *model) { + free(model); +} + +/* --------------------------------------------------- Structure-Factor Calculations -------------------------------------------------- */ + +/* Generate a template for writing structure factors when there's nothing better to use */ +extern void model_generate_template(ReflectionList *reflections, signed int layer) { + + signed int h; + signed int k; + for ( h=-40; h<=40; h++ ) { + for ( k=-40; k<=40; k++ ) { + reflist_addref(reflections, h, k, layer); + } + } + +} + +/* Return f(j, hkl) */ +static double model_sfac(AtomicModel *model, unsigned int j, signed int h, signed int k, signed int l) { + + unsigned int ref; + double x, y, z; + double sfac; + double s; + double a = data_a() * 10; + double b = data_b() * 10; + double c = data_c() * 10; /* Change to Angstroms */ + double gamma = data_gamma(); + double B; + + if ( model->point_atoms ) return 1; + + x = model->atoms[j].x; y = model->atoms[j].y; z = model->atoms[j].z; + ref = model->atoms[j].ref; B = model->atoms[j].B; + + s = resolution(h, k, l, a, b, c, gamma); + + /* Calculate f(j,hkl) */ + sfac = elements[ref].sfac_c; + sfac += elements[ref].sfac_a1 * exp(-elements[ref].sfac_b1 * s * s); + sfac += elements[ref].sfac_a2 * exp(-elements[ref].sfac_b2 * s * s); + sfac += elements[ref].sfac_a3 * exp(-elements[ref].sfac_b3 * s * s); + sfac += elements[ref].sfac_a4 * exp(-elements[ref].sfac_b4 * s * s); + + /* Thermal factor */ + sfac = sfac * exp(- B * s * s); + + return sfac; + +} + +/* Real part of a particular atom's (and it's equivalents') contribution to F */ +static double model_f_contribution_re(AtomicModel *model, unsigned int j, signed int h, signed int k, signed int l) { + + double sfac, cont; + AtomicModel *equivalents; + size_t i; + + sfac = model_sfac(model, j, h, k, l); + equivalents = symmetry_generate_equivalent_atoms(model, j, MODEL_TARGET_CALCULATION); + + cont = 0; + for ( i=0; in_atoms; i++ ) { + double dot, x, y, z; + x = equivalents->atoms[i].x; y = equivalents->atoms[i].y; z = equivalents->atoms[i].z; + dot = h * (1-x) + k * (1-y) + l * z; + cont += model->atoms[j].occ * sfac * cos(-2 * M_PI * dot); + } + + model_free(equivalents); + + return cont; + +} + +/* Imaginary part of a particular atom's (and it's equivalents') contribution to F */ +static double model_f_contribution_im(AtomicModel *model, unsigned int j, signed int h, signed int k, signed int l) { + + double sfac, cont; + AtomicModel *equivalents; + size_t i; + + sfac = model_sfac(model, j, h, k, l); + equivalents = symmetry_generate_equivalent_atoms(model, j, MODEL_TARGET_CALCULATION); + + cont = 0; + for ( i=0; in_atoms; i++ ) { + double dot, x, y, z; + x = equivalents->atoms[i].x; y = equivalents->atoms[i].y; z = equivalents->atoms[i].z; + dot = h * (1-x) + k * (1-y) + l * z; + cont += model->atoms[j].occ * sfac * sin(-2 * M_PI * dot); + } + + model_free(equivalents); + + return cont; + +} + +double model_mod_f(AtomicModel *model, signed int h, signed int k, signed int l) { + + double sf_re, sf_im; + unsigned int j; + + sf_re = 0; sf_im = 0; + /* Work out the total contribution */ + for ( j=0; jn_atoms; j++ ) { + if ( model->atoms[j].active ) { + sf_re += model_f_contribution_re(model, j, h, k, l); + sf_im += model_f_contribution_im(model, j, h, k, l); + } + } + + return sqrt((sf_re*sf_re) + (sf_im*sf_im)); + +} + +/* Calculate kinematical structure factor amplitudes. + * "layer" is meaningless and is ignored if "given_template" is not NULL */ +ReflectionList *model_calculate_f(ReflectionList *given_template, AtomicModel *given_model, signed int layer) { + + ReflectionList *model_reflections; + unsigned int i; + AtomicModel *model; + ReflectionList *template; + + model_reflections = reflist_new(); + + if ( !given_template ) { + template = reflist_new(); + model_generate_template(template, layer); + } else { + template = given_template; + } + + if ( !given_model ) { + model = model_get_current(); + } else { + model = given_model; + } + + if ( model->thickness > 0.0001 ) { + + /* Dynamical diffraction amplitudes */ + model_reflections = multislice_calculate_f_dyn(model, template); + + } else { + + /* Kinematical structure factors */ + for ( i=1; in_reflections; i++ ) { /* 'hkl' loop */ + + signed int h = template->refs[i].h; + signed int k = template->refs[i].k; + signed int l = template->refs[i].l; + double sf_re = 0; + double sf_im = 0; + unsigned int j; + + /* Work out the total contribution */ + for ( j=0; jn_atoms; j++ ) { + if ( model->atoms[j].active ) { + sf_re += model_f_contribution_re(model, j, h, k, l); + sf_im += model_f_contribution_im(model, j, h, k, l); + } + } + + //printf("%3i %3i %3i am=%f ph=%f re=%f im=%f\n", h, k, l, sqrt((sf_re*sf_re) + (sf_im*sf_im)), atan2(sf_im, sf_re), sf_re, sf_im); + reflist_addref_am_ph(model_reflections, h, k, l, sqrt((sf_re*sf_re) + (sf_im*sf_im)), atan2(sf_im, sf_re)); + reflist_set_components(model_reflections, h, k, l, sf_re, sf_im); + + } + + } + + if ( given_template == NULL ) free(template); + + return model_reflections; + +} + +#if 0 +/* Apply a smoothing filter to the amplitudes */ +static void model_smooth(ReflectionList *reflections) { + + unsigned int i; + float a = data_a(); + float b = data_b(); + float c = data_c(); + double sigma = 5.19; + + for ( i=1; in_reflections; i++ ) { + + double s; + signed int h, k, l; + + h = reflections->refs[i].h; + k = reflections->refs[i].k; + l = reflections->refs[i].l; + + s = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)) + ((l*l)/(c*c))); + + reflections->refs[i].amplitude *= exp(-(s*s)/(2*sigma*sigma)); + + } + + +} +#endif + +/* For the difference synthesis */ +void model_calculate_difference_coefficients(ReflectionList *reflections) { + + unsigned int i; + double scale; + ReflectionList *model_reflections; + + if ( reflections->n_reflections == 0 ) return; /* No reflections */ + model_reflections = model_calculate_f(reflections, NULL, 69); + if ( !model_reflections ) return; + + scale = stat_scale(reflections, model_reflections); + //printf("Scale for difference coefficients = %f\n", scale); + + for ( i=0; in_reflections; i++ ) { /* 'hkl' loop */ + + if ( model_reflections->refs[i].amplitude > 0.0000001 ) { + reflections->refs[i].amplitude = reflections->refs[i].amplitude - (scale * model_reflections->refs[i].amplitude); + reflections->refs[i].phase_known = model_reflections->refs[i].phase_known; + reflections->refs[i].phase_known_set = 1; + } else { + reflections->refs[i].phase_known_set = 0; + } + + } + + //model_smooth(reflections); + + free(model_reflections); + +} + +/* For the Fourier refinement synthesis */ +void model_calculate_refinement_coefficients(ReflectionList *reflections) { + + unsigned int i; + ReflectionList *model_reflections; + + if ( reflections->n_reflections == 0 ) return; /* No reflections */ + model_reflections = model_calculate_f(reflections, NULL, 69); + if ( !model_reflections ) return; + + for ( i=0; in_reflections; i++ ) { /* 'hkl' loop */ + + reflections->refs[i].amplitude = reflections->refs[i].amplitude; + reflections->refs[i].phase_known = model_reflections->refs[i].phase_known; + + } + + //model_smooth(reflections); + + free(model_reflections); + +} + +/* ---------------------------------------------------------------- Misc -------------------------------------------------------------- */ + +void model_open_editor() { + model_current->editor = model_editor_open(model_current); +} + +AtomicModel *model_get_current() { + return model_current; +} + +/* Notify that a model has been changed. + * Special version for when a ModelEditor initiates the update + * (don't try and update the editor in this case!) */ +void model_notify_update_editor(AtomicModel *model) { + + if ( model == model_current ) { + displaywindow_switchview(); + } + +} + +/* Notify that a model has been changed */ +void model_notify_update(AtomicModel *model) { + + model_notify_update_editor(model); + + if ( model->editor ) { + model_editor_get_model(model->editor); + } + +} + +/* Initial setting of the current model. Many functions assume that + * a model structure exists. */ +void model_default() { + model_current = model_new(); +} + +void model_lock(AtomicModel *model) { + model->lock_count++; + if ( model->editor ) { + model_editor_lock(model->editor); + } +} + +void model_unlock(AtomicModel *model) { + model->lock_count--; + if ( model->lock_count == 0 ) { + if ( model->editor ) { + model_editor_unlock(model->editor); + } + } +} + +int model_current_is_blank() { + if ( model_get_current()->n_atoms) { + return 0; + } + return 1; +} -- cgit v1.2.3