from __future__ import division from scitbx.array_family import flex import sys from mmtbx import bulk_solvent from cctbx import adptbx import boost.python ext = boost.python.import_ext("mmtbx_f_model_ext") from cctbx import sgtbx from mmtbx.bulk_solvent import bulk_solvent_and_scaling import mmtbx.f_model import math from libtbx import group_args import scitbx.math from cctbx import miller import mmtbx.arrays import scitbx.math import cctbx.xray from libtbx import group_args from libtbx.test_utils import approx_equal def moving_average(x): x_ = x_ = [x[0]] + list(x) + [x[len(x)-1]] for cycle in xrange(5): result = x_[:] selection = flex.bool(len(result), False) for i, s in enumerate(selection): if(i!=0 and i!=len(result)-1): if((result[i-1]result[i] and result[i+1]>result[i])): selection[i]=True for i in xrange(len(result)): if(i!=0 and i!=len(result)-1 and selection[i]): result[i] = (x_[i-1]+x_[i]+x_[i+1])/3. x_ = result[:] return result[1:len(result)-1] class run(object): def __init__(self, f_obs, f_calc, # can be a sum: f_calc=f_hydrogens+f_calc+f_part f_mask, # only one shell is supported r_free_flags, ss, bin_selections=None, scale_method="combo", number_of_cycles=20, # termination occures much earlier auto_convergence_tolerance = 1.e-4, log=None, auto=True, auto_convergence=True, bulk_solvent = True, try_poly = True, try_expanal = True, try_expmin = False, verbose=False): assert f_obs.indices().all_eq(r_free_flags.indices()) self.log = log self.scale_method = scale_method self.verbose = verbose self.r_free_flags = r_free_flags self.ss = ss self.bulk_solvent = bulk_solvent self.try_poly = try_poly self.try_expanal = try_expanal self.try_expmin = try_expmin self.d_spacings = f_obs.d_spacings() self.r_low = None self.poly_approx_cutoff = None self.f_obs = f_obs self.r_free_flags = r_free_flags self.auto_convergence = auto_convergence self.auto_convergence_tolerance = auto_convergence_tolerance self.scale_matrices = None self.auto = auto self.bin_selections = bin_selections self.k_exp_overall, self.b_exp_overall = None,None if(self.bin_selections is None): self.bin_selections = self.f_obs.log_binning() # If R-free flags are bad - discard them and use all reflections instead. ifg = self.is_flags_good() if(ifg): self.selection_work = miller.array( miller_set = self.f_obs, data = ~self.r_free_flags.data()) else: self.selection_work = miller.array( miller_set = self.f_obs, data = flex.bool(self.f_obs.data().size(), True)) # self.ss = ss def init_result(): return group_args( k_mask_bin_orig = None, k_mask_bin_smooth = None, k_mask = None, k_isotropic = None, k_mask_fit_params = None) self.bss_result = init_result() if(log is None): log = sys.stdout if(verbose): print >> log, "-"*80 print >> log, \ "Overall, iso- and anisotropic scaling and bulk-solvent modeling:" point_group = sgtbx.space_group_info( symbol=f_obs.space_group().type().lookup_symbol() ).group().build_derived_point_group() self.adp_constraints = sgtbx.tensor_rank_2_constraints( space_group=point_group, reciprocal_space=True) self.core = mmtbx.arrays.init(f_calc = f_calc, f_masks = f_mask) if(abs(self.core.f_mask()).data().all_eq(0)): self.bulk_solvent=False self.cores_and_selections = [] self.low_resolution_selection = self._low_resolution_selection() self.high_resolution_selection = self._high_resolution_selection() if(verbose): print >> log, " Using %d resolution bins"%len(self.bin_selections) self.ss_bin_values=[] sel_positive = self.f_obs.data()>0 self.selection_work = self.selection_work.customized_copy( data = self.selection_work.data() & sel_positive) for i_sel, sel in enumerate(self.bin_selections): core_selected = self.core.select(selection=sel) sel_use = self.selection_work.data().select(sel) sel_work = sel & self.selection_work.data() self.cores_and_selections.append([sel, core_selected, sel_use, sel_work]) ss = self.ss.select(sel) self.ss_bin_values.append([ flex.min(ss), flex.max(ss), flex.mean(ss)]) for cycle in xrange(number_of_cycles): r_start = self.r_factor(use_scale=True) r_start0 = r_start use_scale_r=False if(cycle==0): use_scale_r=True if(verbose): print >> log, " cycle %d:"%cycle print >> log, " r(start): %6.4f"%(r_start) # bulk-solvent and overall isotropic scale if(self.bulk_solvent): if(cycle==0): r_start = self.set_k_isotropic_exp(r_start = r_start, use_scale_r=use_scale_r, verbose=verbose) r_start = self.k_mask_grid_search(r_start=r_start) if(verbose): print >> self.log, " r(bulk_solvent_grid_search): %6.4f"%r_start r_start = self.set_k_isotropic_exp(r_start = r_start, use_scale_r=use_scale_r, verbose=verbose) else: r_start = self.bulk_solvent_scaling(r_start = r_start) # anisotropic scale if([try_poly, try_expanal, try_expmin].count(True)): if(verbose): print >> log, " anisotropic scaling:" self.anisotropic_scaling(r_start = r_start) if(self.auto_convergence and self.is_converged(r_start=r_start0, tolerance=self.auto_convergence_tolerance)): break self.apply_overall_scale() if(verbose): print >> log, " r(final): %6.4f"%(self.r_factor()) self.show() # self.r_low = self._r_low() self.r_high = self._r_high() if(verbose): d = self.d_spacings.data().select(self.low_resolution_selection) d1 = ("%7.4f"%flex.min(d)).strip() d2 = ("%7.4f"%flex.max(d)).strip() n = d.size() print >> self.log, "r(low-resolution: %s-%s A; %d reflections): %6.4f"%( d2,d1,n,self.r_low) print >> log, "-"*80 self.r_final = self.r_factor() def is_flags_good(self): """ This function detects inadequate R-free flags. """ result = True fd = self.r_free_flags.data() for i_sel, sel in enumerate(self.bin_selections): fd_ = fd.select(sel) if(fd_.count(True)==0): result = False break return result def set_k_isotropic_exp(self, r_start,use_scale_r, verbose, b_lower_limit = -100): arrays = mmtbx.arrays.init( f_calc = self.core.f_calc, f_masks = self.core.f_mask(), k_isotropic = self.core.k_isotropic, k_anisotropic = self.core.k_anisotropic, k_masks = self.core.k_mask()) sel = self.selection_work.data() rf = scitbx.math.gaussian_fit_1d_analytical( x = flex.sqrt(self.ss).select(sel), y = self.f_obs.data().select(sel), z = abs(arrays.f_model).data().select(sel)) if(rf.b < b_lower_limit): return r_start k = rf.a * flex.exp(-self.ss * rf.b) r = self.try_scale(k_isotropic_exp = k, use_scale_r=use_scale_r) if(r> self.log, " r(set_k_isotropic_exp): %6.4f"%r self.core = self.core.update(k_isotropic_exp = k) r_start = self.r_factor() else: if(verbose): print >> self.log, " r(set_k_isotropic_exp): %6.4f (rejected)"%r return r_start def try_scale(self, k_isotropic_exp=None, k_isotropic=None, k_mask=None, k_anisotropic=None, selection=None, use_scale_r=False): if(k_isotropic_exp is None): k_isotropic_exp = self.core.k_isotropic_exp if(k_isotropic is None): k_isotropic = self.core.k_isotropic if(k_mask is None): k_mask = self.core.k_mask() if(k_anisotropic is None): k_anisotropic = self.core.k_anisotropic c = mmtbx.arrays.init( f_calc = self.core.f_calc, f_masks = self.core.f_mask(), k_isotropic_exp = k_isotropic_exp, k_isotropic = k_isotropic, k_anisotropic = k_anisotropic, k_masks = k_mask) sel = self.selection_work.data() if(not use_scale_r): if(selection is not None): sel = selection & sel return bulk_solvent.r_factor(self.f_obs.data(), c.f_model.data(), sel) # XXX else: return cctbx.xray.r_factor( self.f_obs.data().select(sel), c.f_model.data().select(sel)).value() def r_factor(self, use_scale=False): if(use_scale): return bulk_solvent.r_factor(self.f_obs.data(), self.core.f_model.data(), self.selection_work.data()) else: return bulk_solvent.r_factor(self.f_obs.data(), self.core.f_model.data(), self.selection_work.data(), 1) def _r_high(self): return bulk_solvent.r_factor(self.f_obs.data(), self.core.f_model.data(), self.high_resolution_selection, 1) def _r_low(self): return bulk_solvent.r_factor(self.f_obs.data(), self.core.f_model.data(), self.low_resolution_selection, 1) def _high_resolution_selection(self): return self.bin_selections[len(self.bin_selections)-1] & \ self.selection_work.data() def _low_resolution_selection(self): d = self.d_spacings.data() result = None if((d>=8).count(True)>=500): result = d>=8 else: sel = self.bin_selections[0].deep_copy() if(sel.count(True)>=500): result = sel else: for i in range(1, len(self.bin_selections)): #XXX sel |= self.bin_selections[i].deep_copy() if(sel.count(True)>=300): break result = sel return result & self.selection_work.data() def populate_bin_to_individual_k_mask_linear_interpolation(self, k_mask_bin): assert len(k_mask_bin) == len(self.cores_and_selections) def linear_interpolation(x1,x2,y1,y2): k=0 if(x1!=x2): k=(y2-y1)/(x2-x1) b = y1-k*x1 return k,b result1 = flex.double(self.f_obs.size(), -1) result2 = flex.double(self.f_obs.size(), -1) result = flex.double(self.f_obs.size(), -1) for i, cas in enumerate(self.cores_and_selections): selection, zzz, zzz, zzz = cas x1,x2 = self.ss_bin_values[i][0], self.ss_bin_values[i][1] y1 = k_mask_bin[i] if(i==len(k_mask_bin)-1): y2 = k_mask_bin[i-1] else: y2 = k_mask_bin[i+1] k,b = linear_interpolation(x1=x1,x2=x2,y1=y1,y2=y2) bulk_solvent.set_to_linear_interpolated(self.ss,k,b,selection,result1) result2.set_selected(selection, y1) r1 = self.try_scale(k_mask = result1, selection=selection) # XXX inefficient r2 = self.try_scale(k_mask = result2, selection=selection) # XXX inefficient if(r12): x=flex.double() y=flex.double() for i_sel, cas in enumerate(self.cores_and_selections): selection, core, selection_use, sel_work = cas sel = sel_work ss_ = self.ss_bin_values[i_sel][2] k_isotropic_ = flex.mean(self.core.k_isotropic.select(sel)) x.append(ss_) y.append(k_isotropic_) import scitbx.math r = scitbx.math.gaussian_fit_1d_analytical(x = flex.sqrt(x), y = y) r_start = self.r_factor() k_isotropic = r.a*flex.exp(-self.ss*r.b) r = self.try_scale(k_isotropic = k_isotropic) if(rself.r_final and abs(r_start-self.r_final)> self.log, " Statistics in resolution bins:" #assert k_mask.size() == len(self.bin_selections) fmt=" %7.5f %6.2f -%6.2f %5.1f %5d %-6s %-6s %-6s %6.3f %6.3f %8.2f %6.4f" f_model = self.core.f_model.data() print >> self.log, " s^2 Resolution Compl Nrefl k_mask k_iso k_ani R" print >> self.log, " (A) (%) orig smooth average" k_mask_bin_orig_ = str(None) k_mask_bin_smooth_ = str(None) k_mask_bin_approx_ = str(None) for i_sel, cas in enumerate(self.cores_and_selections): selection, core, selection_use, sel_work = cas sel = sel_work ss_ = self.ss_bin_values[i_sel][2] if(b is not None and self.bss_result.k_mask_bin_orig is not None): k_mask_bin_orig_ = "%6.4f"%self.bss_result.k_mask_bin_orig[i_sel] if(b is not None and self.bss_result.k_mask_bin_smooth is not None): k_mask_bin_smooth_ = "%6.4f"%self.bss_result.k_mask_bin_smooth[i_sel] k_mask_bin_averaged_ = "%6.4f"%flex.mean(self.core.k_mask().select(sel)) d_ = self.d_spacings.data().select(sel) d_min_ = flex.min(d_) d_max_ = flex.max(d_) n_ref_ = d_.size() f_obs_ = self.f_obs.select(sel) f_obs_mean_ = flex.mean(f_obs_.data()) k_isotropic_ = flex.mean(self.core.k_isotropic.select(sel)) k_anisotropic_ = flex.mean(self.core.k_anisotropic.select(sel)) cmpl_ = f_obs_.completeness(d_max=d_max_)*100. r_ = bulk_solvent.r_factor(f_obs_.data(),f_model.select(sel),1) print >> self.log, fmt%(ss_, d_max_, d_min_, cmpl_, n_ref_, k_mask_bin_orig_, k_mask_bin_smooth_,k_mask_bin_averaged_, k_isotropic_, k_anisotropic_, f_obs_mean_, r_) def _k_isotropic_as_scale_k1(self, r_start, k_mask=None): k_isotropic = flex.double(self.ss.size(), -1) if(k_mask is None): k_mask = self.core.k_mask() core_data = mmtbx.arrays.init( f_calc = self.core.f_calc, f_masks = self.core.f_mask(), k_isotropic_exp = self.core.k_isotropic_exp, k_anisotropic = self.core.k_anisotropic, k_masks = k_mask).data for i_cas, cas in enumerate(self.cores_and_selections): selection, core, selection_use, sel_work = cas scale_k1 = bulk_solvent.scale(self.f_obs.data(), core_data.f_model, sel_work) k_isotropic = k_isotropic.set_selected(selection, scale_k1) assert k_isotropic.count(-1.) == 0 # r = self.try_scale(k_isotropic = k_isotropic, k_mask=k_mask) if(r<=r_start): self.core = self.core.update(k_isotropic = k_isotropic, k_masks=k_mask) r_start = r return k_isotropic def estimate_scale_k1(self, cutoff=4, width=1, min_reflections=500): cutoff = min(cutoff, self.f_obs.d_min()+width) sel_high = self.d_spacings.data()min_reflections): core = self.core.select(selection = sel_high) f_obs = self.f_obs.select(sel_high) fm = core.k_isotropic_exp * core.k_anisotropic * (core.f_calc.data() + core.k_mask() * core.f_mask().data()) scale_k1 = bulk_solvent.scale(f_obs.data(), fm) return scale_k1 def bulk_solvent_scaling(self, r_start): k_mask = flex.double(self.f_obs.size(), -1) k_mask_bin = flex.double() k_mask_trial_range = flex.double([i/1000. for i in range(0,1000,10)]) scale_k1 = self.estimate_scale_k1() def get_k_mask_trial_range(x, shift=0.05): if(x > 1): x = 1 result = flex.double() inc = max(0,x-shift) while inc<=x+shift+1.e-3: result.append(inc) inc+=0.01 return result for i_cas, cas in enumerate(self.cores_and_selections): selection, core, selection_use, sel_work = cas scale = self.core.k_anisotropic.select(selection)*scale_k1*\ self.core.k_isotropic_exp.select(selection) f_obs = self.f_obs.select(selection).data() f_calc = core.f_calc.data()*scale f_mask = core.f_mask().data()*scale if(self.scale_method == "k_iso_k_mask_anal"): obj = bulk_solvent.overall_and_bulk_solvent_scale_coefficients_analytical( f_obs = f_obs, f_calc = f_calc, f_mask = f_mask, selection = selection_use) k_mask_bin.append(obj.x_best) k_mask.set_selected(selection, obj.x_best) elif(self.scale_method == "k_mask_anal"): obj = bulk_solvent.bulk_solvent_scale_coefficients_analytical( f_obs = f_obs, f_calc = f_calc, f_mask = f_mask, selection = selection_use) k_mask_bin.append(obj.x_best) k_mask.set_selected(selection, obj.x_best) elif(self.scale_method == "k_mask_r_grid_search"): k_mask_bin_, k_isotropic_bin_ = \ bulk_solvent.k_mask_and_k_overall_grid_search( f_obs, f_calc, f_mask, k_mask_trial_range, selection_use) k_mask_bin.append(k_mask_bin_) k_mask.set_selected(selection, k_mask_bin_) elif(self.scale_method == "combo"): r = flex.double() k = flex.double() x0 = self.bss_result.k_mask_bin_orig[i_cas] k_mask.set_selected(selection, x0) r0 = self.try_scale(k_mask = k_mask, selection=selection) r.append(r0) k.append(x0) # obj1 = bulk_solvent.overall_and_bulk_solvent_scale_coefficients_analytical( f_obs = f_obs, f_calc = f_calc, f_mask = f_mask, selection = selection_use) k_mask.set_selected(selection, obj1.x_best) r.append(self.try_scale(k_mask = k_mask, selection=selection)) k.append(obj1.x_best) # obj2 = bulk_solvent.bulk_solvent_scale_coefficients_analytical( f_obs = f_obs, f_calc = f_calc, f_mask = f_mask, selection = selection_use) k_mask.set_selected(selection, obj2.x_best) r.append(self.try_scale(k_mask = k_mask, selection=selection)) k.append(obj2.x_best) s = flex.sort_permutation(r) x = k.select(s)[0] # fine-sample k_mask around minimum of LS to fall into minimum of R k_mask_bin_, k_isotropic_bin_ = \ bulk_solvent.k_mask_and_k_overall_grid_search( f_obs, f_calc, f_mask, get_k_mask_trial_range(x = x), selection_use) k_mask_bin.append(k_mask_bin_) k_mask.set_selected(selection, k_mask_bin_) else: assert 0 # k_mask_bin_smooth = self.smooth(k_mask_bin) k_mask = self.populate_bin_to_individual_k_mask_linear_interpolation( k_mask_bin = k_mask_bin_smooth) k_isotropic = self._k_isotropic_as_scale_k1(r_start=r_start,k_mask = k_mask) # r = self.try_scale(k_mask = k_mask, k_isotropic = k_isotropic) if(r<=r_start or (r>r_start and abs(r-r_start)*100<0.5)): # may be 0.5? self.core = self.core.update(k_isotropic = k_isotropic, k_masks = k_mask) self.bss_result.k_mask_bin_orig = k_mask_bin self.bss_result.k_mask_bin_smooth = k_mask_bin_smooth self.bss_result.k_mask = k_mask self.bss_result.k_isotropic = k_isotropic r_start = r if(self.verbose): print >> self.log, " %s: %6.4f"%("bulk-solvent", r_start) else: if(self.verbose): print >> self.log, " %s: %6.4f (%s)"%("bulk-solvent", r, "result rejected") return self.r_factor() def smooth(self, x): result = moving_average(x = x) result_ = flex.double(len(result), 0) for i, r in enumerate(result): d = 1/math.sqrt(self.ss_bin_values[i][1])/2 if(r==0 and d<3): break result_[i]=r return result_ def format_scale_matrix(self, m=None, log=None): sm = m if(sm is None): sm = self.scale_matrices out = log if(sm is None): print >> log, " k_anisotropic=1" return if(len(sm)<=6): print >> out, " b_cart(11,22,33,12,13,23):",\ ",".join([str("%8.4f"%i).strip() for i in sm]) else: v0=[] v1=[] for i, a in enumerate(sm): if(i in [0,2,4,6,8,10]): v1.append(a) else: v0.append(a) print >> out, " V0:",\ ",".join([str("%8.4f"%i).strip() for i in v0]) print >> out, " V1:",\ ",".join([str("%8.4f"%i).strip() for i in v1]) def anisotropic_scaling(self, r_start): r_expanal, r_poly, r_expmin = None,None,None k_anisotropic_expanal, k_anisotropic_poly, \ k_anisotropic_expmin = None, None, None scale_matrix_expanal, scale_matrix_poly, scale_matrix_expmin= None,None,None sel = self.selection_work.data() f_model = self.core.f_model_no_aniso_scale.data().select(sel) f_obs = self.f_obs.data().select(sel) mi = self.f_obs.indices().select(sel) uc = self.f_obs.unit_cell() mi_all = self.f_obs.indices() # try exp_anal if(self.try_expanal): obj = bulk_solvent.aniso_u_scaler( f_model = f_model, f_obs = f_obs, miller_indices = mi, adp_constraint_matrix = self.adp_constraints.gradient_sum_matrix()) u_star = self.adp_constraints.all_params(tuple(obj.u_star_independent)) scale_matrix_expanal = adptbx.u_as_b(adptbx.u_star_as_u_cart(uc, u_star)) k_anisotropic_expanal = ext.k_anisotropic(mi_all, u_star) r_expanal = self.try_scale(k_anisotropic = k_anisotropic_expanal) if(self.verbose): print >> self.log, " r_expanal: %6.4f"%r_expanal # try poly if(self.try_poly): obj = bulk_solvent.aniso_u_scaler( f_model = f_model, f_obs = f_obs, miller_indices = mi, unit_cell = uc) scale_matrix_poly = obj.a k_anisotropic_poly = ext.k_anisotropic(mi_all, obj.a, uc) r_poly = self.try_scale(k_anisotropic = k_anisotropic_poly) if(self.verbose): print >> self.log, " r_poly : %6.4f"%r_poly # pre-analyze force_to_use_expmin=False if(k_anisotropic_poly is not None and self.auto and r_poly0): force_to_use_expmin = True self.try_expmin = True # try expmin if(self.try_expmin): zero = self.f_obs.select(sel).customized_copy(data = flex.complex_double(f_obs.size(), 0)) fm = mmtbx.f_model.manager_kbu( f_obs = self.f_obs.select(sel), f_calc = self.core.f_model_no_aniso_scale.select(sel), f_masks = [zero], f_part1 = zero, f_part2 = zero, ss = self.ss) obj = bulk_solvent_and_scaling.u_star_minimizer( fmodel_core_data = fm, f_obs = self.f_obs.select(sel), u_initial = [0,0,0,0,0,0], refine_u = True, min_iterations = 500, max_iterations = 500, symmetry_constraints_on_b_cart = True, u_min_max = 500., u_min_min =-500.) u_star = obj.u_min scale_matrix_expmin = adptbx.u_as_b(adptbx.u_star_as_u_cart(uc, u_star)) k_anisotropic_expmin = ext.k_anisotropic(mi_all, u_star) r_expmin = self.try_scale(k_anisotropic = k_anisotropic_expmin) if(self.verbose): print >> self.log, " r_expmin : %6.4f"%r_expmin if(force_to_use_expmin): self.core = self.core.update(k_anisotropic = k_anisotropic_expmin) if(self.verbose): self.format_scale_matrix(m=scale_matrix_expmin) return # select best r = [(r_expanal, k_anisotropic_expanal, scale_matrix_expanal), (r_poly, k_anisotropic_poly, scale_matrix_poly), (r_expmin, k_anisotropic_expmin, scale_matrix_expmin)] r_best = r_start k_anisotropic_best = None scale_matrix_best = None for result in r: r_factor, k_anisotropic, scale_matrix = result if(r_factor is not None and r_factor < r_best): r_best = r_factor k_anisotropic_best = k_anisotropic.deep_copy() scale_matrix_best = scale_matrix[:] if(scale_matrix_best is None): if(self.verbose): print >> self.log, " result rejected due to r-factor increase" else: self.scale_matrices = scale_matrix_best self.core = self.core.update(k_anisotropic = k_anisotropic_best) r_aniso = self.r_factor() if(self.verbose): self.format_scale_matrix() print >> self.log, " r_final : %6.4f"%r_aniso def overall_isotropic_kb_estimate(self): k_total = self.core.k_isotropic * self.core.k_anisotropic * \ self.core.k_isotropic_exp r = scitbx.math.gaussian_fit_1d_analytical(x=flex.sqrt(self.ss), y=k_total) return r.a, r.b def apply_back_trace_of_overall_exp_scale_matrix(self, xray_structure=None): result = None k,b=self.overall_isotropic_kb_estimate() k_total = self.core.k_isotropic * self.core.k_anisotropic * \ self.core.k_isotropic_exp k,b,r = mmtbx.bulk_solvent.fit_k_exp_b_to_k_total(k_total, self.ss, k, b) if(r<0.5): self.k_exp_overall,self.b_exp_overall = k,b if(xray_structure is not None): if([self.k_exp_overall,self.b_exp_overall].count(None)==0 and k != 0): bs1 = xray_structure.extract_u_iso_or_u_equiv()*adptbx.u_as_b(1.) def split(b_trace, xray_structure): b_min = xray_structure.min_u_cart_eigenvalue()*adptbx.u_as_b(1.) b_res = min(0, b_min + b_trace+1.e-6) b_adj = b_trace-b_res xray_structure.shift_us(b_shift = b_adj) return b_adj, b_res b_adj,b_res=split(b_trace=self.b_exp_overall,xray_structure=xray_structure) k_new = self.k_exp_overall*flex.exp(-self.ss*b_adj) result = group_args( xray_structure= xray_structure, k_isotropic = self.core.k_isotropic/k_new*self.core.k_isotropic_exp, k_mask = [m*flex.exp(-self.ss*b_adj) for m in self.core.k_masks], k_anisotropic = self.core.k_anisotropic) bs2 = xray_structure.extract_u_iso_or_u_equiv()*adptbx.u_as_b(1.) diff = bs2-bs1 assert approx_equal(flex.min(diff), flex.max(diff)) assert approx_equal(flex.max(diff), b_adj) return result