From 0f07feaec6e2c96769ded00c7c31ed61bf7e019d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 10 May 2015 23:32:59 +0200 Subject: Refine B factor as well --- src/post-refinement.c | 41 +++++++++++++++++++++++++++++++---------- 1 file changed, 31 insertions(+), 10 deletions(-) (limited to 'src/post-refinement.c') diff --git a/src/post-refinement.c b/src/post-refinement.c index 3008b7af..57d7c5d9 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -174,8 +174,8 @@ static double volume_fraction(double rlow, double rhigh, double pr, } -/* Return the gradient of pG wrt parameter 'k' given the current status - * of the crystal. */ +/* Return the gradient of pGexp(-Bs^2) wrt parameter 'k' given the current + * status of the crystal. */ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { double glow, ghigh; @@ -184,11 +184,20 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) double R = crystal_get_profile_radius(cr); double osf = crystal_get_osf(cr); double gr; + signed int hi, ki, li; + double s, B; + get_indices(refl, &hi, &ki, &li); + s = resolution(crystal_get_cell(cr), hi, ki, li); + B = crystal_get_Bfac(cr); get_partial(refl, &rlow, &rhigh, &p); if ( k == GPARAM_OSF ) { - return p; + return p*exp(-B*s*s); + } + + if ( k == GPARAM_BFAC ) { + return -s*s*osf*p*exp(-B*s*s); } if ( k == GPARAM_R ) { @@ -283,6 +292,12 @@ static void apply_shift(Crystal *cr, int k, double shift) crystal_set_profile_radius(cr, t); break; + case GPARAM_BFAC : + t = crystal_get_Bfac(cr); + t += shift; + crystal_set_Bfac(cr, t); + break; + case GPARAM_OSF : if ( isnan(shift) ) { ERROR("Refusing nan shift of OSF\n"); @@ -335,7 +350,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, const int verbose = 0; int num_params = 0; enum gparam rv[32]; - double G; + double G, B; *n_filtered = 0; @@ -356,7 +371,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, /* If we are scaling, refine scale factors (duh) */ if ( !no_scale ) { rv[num_params++] = GPARAM_OSF; - //rv[num_params++] = GPARAM_BFAC; FIXME: Next + rv[num_params++] = GPARAM_BFAC; } reflections = crystal_get_reflections(cr); @@ -365,6 +380,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, v = gsl_vector_calloc(num_params); G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); /* Construct the equations, one per reflection in this image */ for ( refl = first_refl(reflections, &iter); @@ -376,7 +392,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, double w; double I_partial; int k; - double p, L; + double p, L, s; Reflection *match; double gradients[num_params]; @@ -393,6 +409,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, p = get_partiality(refl); L = get_lorentz(refl); + s = resolution(crystal_get_cell(cr), ha, ka, la); /* Actual measurement of this reflection from this pattern? */ I_partial = get_intensity(refl); @@ -429,7 +446,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, } - delta_I = I_partial - (G * p * I_full / L); + delta_I = I_partial - G*exp(-B*s*s)*p*I_full / L; v_c = w * delta_I * gradients[k]; v_curr = gsl_vector_get(v, k); gsl_vector_set(v, k, v_curr + v_c); @@ -477,6 +494,10 @@ static double pr_iterate(Crystal *cr, const RefList *full, static double guide_dev(Crystal *cr, const RefList *full) { double dev = 0.0; + double G, B; + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); /* For each reflection */ Reflection *refl; @@ -486,7 +507,7 @@ static double guide_dev(Crystal *cr, const RefList *full) refl != NULL; refl = next_refl(refl, iter) ) { - double G, p, L; + double p, L, s; signed int h, k, l; Reflection *full_version; double I_full, I_partial; @@ -503,16 +524,16 @@ static double guide_dev(Crystal *cr, const RefList *full) * scale_intensities() might not yet have been called, so the * full version may not have been calculated yet. */ - G = crystal_get_osf(cr); p = get_partiality(refl); L = get_lorentz(refl); I_partial = get_intensity(refl); I_full = get_intensity(full_version); + s = resolution(crystal_get_cell(cr), h, k, l); //STATUS("%3i %3i %3i %5.2f %5.2f %5.2f %5.2f %5.2f\n", // h, k, l, G, p, I_partial, I_full, // I_partial - p*G*I_full); - dev += pow(I_partial - G*p*I_full/L, 2.0); + dev += pow(I_partial - G*exp(-B*s*s)*p*I_full/L, 2.0); } -- cgit v1.2.3