from __future__ import absolute_import, division, print_function from scitbx.lstbx import normal_eqns_solving from cctbx import miller from cctbx import euclidean_model_matching as emma from cctbx.array_family import flex from cctbx.xray import structure from libtbx.test_utils import approx_equal, show_diff import libtbx.utils from smtbx import masks from smtbx.refinement import constraints import smtbx.utils from six.moves import cStringIO as StringIO from six.moves import range def exercise_masks(): mt = flex.mersenne_twister(seed=0) xs_ref = structure.from_shelx( file=StringIO(YAKRUY_ins)) mi = xs_ref.crystal_symmetry().build_miller_set( d_min=0.5, anomalous_flag=False) fo = mi.structure_factors_from_scatterers( xs_ref, algorithm="direct").f_calc().as_amplitude_array() k = 0.05 + 10 * mt.random_double() fo = fo.customized_copy(data=fo.data()*k) fo2 = fo.f_as_f_sq() acetonitrile_sel = xs_ref.label_selection( 'N4', 'C20', 'C21', 'H211', 'H212', 'H213') xs_no_sol = xs_ref.deep_copy_scatterers().select( acetonitrile_sel, negate=True) # check what happens when no voids are found mask = masks.mask(xs_ref, fo2) mask.compute(solvent_radius=1.2, shrink_truncation_radius=1.2, resolution_factor=1/2, atom_radii_table={'C':1.70, 'B':1.63, 'N':1.55, 'O':1.52}) assert mask.structure_factors() is None assert mask.n_voids() == 0 assert mask.n_solvent_grid_points() == 0 assert mask.f_mask() is None assert mask.f_model() is None assert mask.modified_intensities() is None assert mask.f_000 is None s = StringIO() mask.show_summary(log=s) assert not show_diff(s.getvalue(), """\ use_set_completion: False solvent_radius: 1.20 shrink_truncation_radius: 1.20 van der Waals radii: B C H N O 1.63 1.70 1.20 1.55 1.52 Total solvent accessible volume / cell = 0.0 Ang^3 [0.0%] gridding: (30,45,54) """) # and now with some voids fo2_complete = fo2.sort() fo2_missing_1 = fo2.select_indices(flex.miller_index([(0,0,1), ]), negate=True) mt = flex.mersenne_twister(seed=0) fo2_incomplete = fo2.select(mt.random_bool(fo2.size(), 0.95)) for fo2, use_space_group_symmetry in zip( (fo2_complete, fo2_complete, fo2_missing_1, fo2_incomplete), (True, False, True, True)): if fo2 is fo2_complete: use_set_completion=False else: use_set_completion=True mask = masks.mask(xs_no_sol, fo2, use_set_completion=use_set_completion) mask.compute(solvent_radius=1.2, shrink_truncation_radius=1.2, resolution_factor=1/3, #atom_radii_table={'C':1.70, 'B':1.63, 'N':1.55, 'O':1.52}, use_space_group_symmetry=use_space_group_symmetry) n_voids = flex.max(mask.mask.data) - 1 f_mask = mask.structure_factors() f_model = mask.f_model() modified_fo = mask.modified_intensities().as_amplitude_array() f_obs_minus_f_model = fo.common_set(f_model).f_obs_minus_f_calc(f_obs_factor=1/k, f_calc=f_model) diff_map = miller.fft_map(mask.crystal_gridding, f_obs_minus_f_model) diff_map.apply_volume_scaling() stats = diff_map.statistics() assert n_voids == 2 assert approx_equal(n_voids, mask.n_voids()) assert mask.n_solvent_grid_points() == 42148 if fo2 is fo2_complete: # check the difference map has no large peaks/holes assert max(stats.max(), abs(stats.min())) < 0.11 # expected electron count: 44 assert approx_equal(mask.f_000_s, 44, eps=1) assert modified_fo.r1_factor(mask.f_calc.common_set(modified_fo), k) < 0.006 assert fo.common_set(fo2).r1_factor(f_model.common_set(fo2), k) < 0.006 s = StringIO() mask.show_summary(log=s) assert not show_diff(s.getvalue(), """\ use_set_completion: True solvent_radius: 1.20 shrink_truncation_radius: 1.20 van der Waals radii: C H N O 1.77 1.20 1.50 1.45 Total solvent accessible volume / cell = 146.5 Ang^3 [16.3%] Total electron count / cell = 43.2 gridding: (45,72,80) Void #Grid points Vol/A^3 Vol/% Centre of mass (frac) Eigenvectors (frac) 1 21074 73.3 8.1 ( 0.267, 0.461, 0.672) 1 ( 0.982, 0.126, 0.142) 2 (-0.166, 0.206, 0.964) 3 (-0.092, 0.970,-0.223) 2 21074 73.3 8.1 (-0.267, 0.539, 0.328) 1 ( 0.982, 0.126, 0.142) 2 (-0.166, 0.206, 0.964) 3 (-0.092, 0.970,-0.223) Void Vol/Ang^3 #Electrons 1 73.3 21.6 2 73.3 21.6 """) cif_block = mask.as_cif_block() fo2 = fo.f_as_f_sq() # this bit is necessary until we have constraints, as # otherwise the hydrogens just disappear into the ether. xs = xs_no_sol.deep_copy_scatterers() h_selection = xs.element_selection('H') orig_flags = xs.scatterer_flags() flags = orig_flags.deep_copy() for flag, is_h in zip(flags, h_selection): if is_h: flag.set_grads(False) xs.set_scatterer_flags(flags) # first refine with no mask xs = exercise_least_squares(xs, fo2, mask=None) xs.set_scatterer_flags(orig_flags) for i in range(1): # compute improved mask/f_mask mask = masks.mask(xs, fo2) mask.compute(solvent_radius=1.2, shrink_truncation_radius=1.2, atom_radii_table={'C':1.70, 'B':1.63, 'N':1.55, 'O':1.52}, resolution_factor=1/3) mask.structure_factors() xs = exercise_least_squares(xs, fo2, mask) # again exclude hydrogens from tests because of lack of constraints emma_ref = xs_no_sol.select(h_selection, negate=True).as_emma_model() match = emma.model_matches(emma_ref, xs.select( h_selection, negate=True).as_emma_model()).refined_matches[0] assert approx_equal(match.rms, 0, eps=1e-3) def exercise_least_squares(xray_structure, fo_sq, mask=None): from smtbx.refinement import least_squares fo_sq = fo_sq.customized_copy(sigmas=flex.double(fo_sq.size(),1.)) xs = xray_structure.deep_copy_scatterers() if mask is not None: f_mask = mask.f_mask() else: f_mask = None connectivity_table = smtbx.utils.connectivity_table(xs) reparametrisation = constraints.reparametrisation( structure=xs, constraints=[], connectivity_table=connectivity_table) obs = fo_sq.as_xray_observations() ls = least_squares.crystallographic_ls( obs, reparametrisation, f_mask=f_mask, weighting_scheme="default") cycles = normal_eqns_solving.naive_iterations(ls, n_max_iterations=3) return xs # durham structure 99srv215, CSD code YAKRUY YAKRUY_ins = """ CELL 0.71073 7.086 10.791 12.850 104.16 105.87 95.86 ZERR 2.00 0.001 0.003 0.002 0.02 0.02 0.01 LATT 1 SFAC C H N O UNIT 42 40 8 4 O1 4 0.263832 0.309851 0.916785 11.00000 0.03486 0.02047 = 0.02294 0.00910 0.01221 0.00848 O2 4 0.259857 0.121826 0.750368 11.00000 0.03493 0.02585 = 0.01911 0.00907 0.01212 0.00849 N1 3 0.244867 -0.219643 1.065706 11.00000 0.02111 0.01983 = 0.02052 0.00529 0.00767 0.00528 N2 3 0.238512 -0.124207 1.256893 11.00000 0.02403 0.01605 = 0.01797 0.00354 0.00780 0.00428 N3 3 0.236863 -0.003326 1.323510 11.00000 0.02340 0.01743 = 0.02104 0.00455 0.00803 0.00515 N4 3 0.313652 0.621206 0.720884 11.00000 0.05212 0.03488 = 0.03702 0.00851 0.01917 0.00007 C1 1 0.258790 0.192265 0.938230 11.00000 0.01801 0.02150 = 0.02273 0.00950 0.00637 0.00518 C2 1 0.258684 0.087261 0.845526 11.00000 0.01829 0.02676 = 0.01813 0.00919 0.00712 0.00512 C3 1 0.258107 -0.035539 0.857494 11.00000 0.01763 0.02320 = 0.01795 0.00465 0.00660 0.00526 H3 2 0.265419 -0.106192 0.797272 11.00000 0.02235 C4 1 0.252521 -0.058937 0.960660 11.00000 0.01442 0.02188 = 0.01933 0.00569 0.00576 0.00362 C5 1 0.249654 0.044792 1.051770 11.00000 0.01428 0.02074 = 0.01927 0.00546 0.00489 0.00333 C6 1 0.254907 0.171444 1.038876 11.00000 0.01823 0.01960 = 0.01938 0.00441 0.00607 0.00410 H6 2 0.260661 0.244361 1.101313 11.00000 0.02028 C7 1 0.244340 0.011338 1.151454 11.00000 0.01528 0.01760 = 0.02039 0.00560 0.00559 0.00358 C8 1 0.242078 -0.118214 1.151436 11.00000 0.01632 0.02066 = 0.01923 0.00609 0.00594 0.00402 C9 1 0.250274 -0.186707 0.973843 11.00000 0.01810 0.02164 = 0.01942 0.00361 0.00684 0.00426 H9 2 0.252942 -0.255996 0.908935 11.00000 0.01373 C10 1 0.238915 0.077502 1.260637 11.00000 0.01646 0.02046 = 0.02005 0.00660 0.00460 0.00293 C11 1 0.241179 -0.230214 1.303905 11.00000 0.01923 0.01988 = 0.02379 0.01016 0.00832 0.00634 C12 1 0.301101 -0.204248 1.421081 11.00000 0.02960 0.02212 = 0.02187 0.00694 0.00858 0.00529 H12 2 0.347647 -0.113158 1.469110 11.00000 0.03067 C13 1 0.297066 -0.306624 1.468541 11.00000 0.03346 0.03026 = 0.02258 0.01182 0.00920 0.00828 H13 2 0.338439 -0.286559 1.549032 11.00000 0.03493 C14 1 0.238170 -0.433815 1.401003 11.00000 0.03483 0.02587 = 0.03168 0.01598 0.01317 0.01089 H14 2 0.239203 -0.504609 1.435015 11.00000 0.03255 C15 1 0.182288 -0.458453 1.284537 11.00000 0.03573 0.01878 = 0.02946 0.00736 0.01196 0.00646 H15 2 0.141861 -0.547777 1.233813 11.00000 0.03479 C16 1 0.182000 -0.357827 1.235343 11.00000 0.02803 0.02211 = 0.02063 0.00598 0.00799 0.00558 H16 2 0.140443 -0.375280 1.154110 11.00000 0.02660 C17 1 0.228172 0.216669 1.308393 11.00000 0.02843 0.01824 = 0.02238 0.00458 0.00869 0.00462 H171 2 0.343696 0.277860 1.305773 11.00000 0.03122 H172 2 0.233564 0.235853 1.390463 11.00000 0.03214 H173 2 0.103775 0.237732 1.266193 11.00000 0.03633 C18 1 0.262699 0.418567 1.006264 11.00000 0.03735 0.01968 = 0.02610 0.00669 0.01069 0.00556 H181 2 0.384186 0.436044 1.074408 11.00000 0.02881 H182 2 0.136744 0.404648 1.027622 11.00000 0.02787 H183 2 0.269188 0.493767 0.976527 11.00000 0.02735 C19 1 0.246327 0.019380 0.652743 11.00000 0.03131 0.03014 = 0.01841 0.00608 0.01024 0.00667 H191 2 0.369862 -0.025013 0.666318 11.00000 0.03203 H192 2 0.122807 -0.042860 0.628185 11.00000 0.03308 H193 2 0.249416 0.060472 0.594291 11.00000 0.03051 C20 1 0.295027 0.510549 0.696308 11.00000 0.03015 0.03545 = 0.02324 0.01128 0.01118 0.00361 C21 1 0.271551 0.369674 0.665707 11.00000 0.03896 0.03186 = 0.04785 0.01623 0.01900 0.00949 H211 2 0.372427 0.342541 0.630638 11.00000 0.06494 H212 2 0.140789 0.331654 0.610485 11.00000 0.06386 H213 2 0.292381 0.339113 0.730897 11.00000 0.08983 HKLF 4 END """ def run(): libtbx.utils.show_times_at_exit() exercise_masks() if __name__ == '__main__': run()