from __future__ import absolute_import, division, print_function from cctbx import xray from cctbx.development import random_structure from cctbx.development import debug_utils from cctbx.array_family import flex import sys, random from six.moves import range def shift_u_iso(structure, shift): for sc in structure.scatterers(): if(sc.flags.use_u_iso() and sc.flags.grad_u_iso()): sc.u_iso += (shift * random.random()) def shift_u_aniso(structure, shift): for sc in structure.scatterers(): if(sc.flags.use_u_aniso() and sc.flags.grad_u_aniso()): u_star = list(sc.u_star) for i in range(6): u_star[i] += (shift * random.random()) sc.u_star = u_star def exercise(target_functor, data_type, space_group_info, anomalous_flag, gradient_flags, occupancy_penalty, n_elements=9, d_min=None, shake_sigma=0.1, verbose=0,tan_u_iso=False, param = 0): assert data_type == 'F' or data_type == 'F^2' if (data_type == 'F^2' and not target_functor == xray.unified_least_squares_residual): return structure_ideal = random_structure.xray_structure( space_group_info, elements=(("O","N","C")*(n_elements))[:n_elements],#("Se",)*n_elements, volume_per_atom=200, random_u_iso=True, random_u_cart_scale=.3, random_occupancy=True, use_u_aniso = True) random_structure.random_modify_adp_and_adp_flags( scatterers = structure_ideal.scatterers(), random_u_iso_scale = 0.3, random_u_iso_min = 0.0) xray.set_scatterer_grad_flags(scatterers = structure_ideal.scatterers(), site = gradient_flags.site, u_iso = gradient_flags.u_iso, u_aniso = gradient_flags.u_aniso, occupancy = gradient_flags.occupancy, fp = gradient_flags.fp, fdp = gradient_flags.fdp, tan_u_iso = tan_u_iso, param = param) if(0): print() for sc in structure_ideal.scatterers(): print(sc.flags.use_u_iso(),sc.flags.grad_u_iso(),sc.flags.use_u_aniso(),\ sc.flags.grad_u_aniso(),sc.u_iso, sc.u_star,sc.flags.tan_u_iso(),\ sc.flags.param, sc.occupancy) rnd_f_calc = structure_ideal.structure_factors( anomalous_flag=anomalous_flag, d_min=d_min, algorithm="direct", cos_sin_table=True).f_calc() if data_type == "F": y_obs = abs(rnd_f_calc) elif data_type == "F^2": y_obs = rnd_f_calc.norm() y_obs.set_observation_type_xray_intensity() else: raise "Error: invalid data type: %s" % data_type if (0 or verbose): print("structure_ideal:") structure_ideal.show_summary().show_scatterers() print("n_special_positions:", \ structure_ideal.special_position_indices().size()) print() structure_ideal_cp = structure_ideal.deep_copy_scatterers() structure_shake = structure_ideal if (gradient_flags.site): structure_shake = structure_shake.random_modify_parameters( "site", shake_sigma) if (gradient_flags.occupancy): structure_shake = structure_shake.random_modify_parameters( "occupancy", shake_sigma) if (occupancy_penalty is not None): structure_shake.scatterers()[-1].occupancy = 0 if (gradient_flags.u_aniso): shift_u_aniso(structure_shake, 0.001) if (gradient_flags.u_iso): shift_u_iso(structure_shake, 0.1) assert tuple(structure_ideal.special_position_indices()) \ == tuple(structure_shake.special_position_indices()) if (0 or verbose): print("structure_shake:") structure_shake.show_summary().show_scatterers() print() for i_trial in range(10): try: minimizer = xray.minimization.lbfgs( target_functor=target_functor(y_obs), xray_structure=structure_shake, occupancy_penalty=occupancy_penalty, structure_factor_algorithm="direct") except RuntimeError as e: if (str(e).find("debye_waller_factor_exp: arg_limit exceeded") < 0): raise else: break else: raise RuntimeError("Too many xray.minimization.lbfgs failures.") if (0 or verbose): print("first:", minimizer.first_target_value) print("final:", minimizer.final_target_value) print() assert minimizer.final_target_value < minimizer.first_target_value if (0 or verbose): print("minimized structure_shake:") structure_shake.show_summary().show_scatterers() print() f_final = y_obs.structure_factors_from_scatterers( xray_structure=structure_shake, algorithm="direct", cos_sin_table=True).f_calc() if data_type == 'F': f_final = abs(f_final) elif data_type == 'F^2': f_final = f_final.norm() c = flex.linear_correlation(y_obs.data(), f_final.data()) assert c.is_well_defined() if (0 or verbose): label = gradient_flags.string_of_true() if (anomalous_flag): label += ",anomalous" print("correlation: %10.8f" % c.coefficient(), label) print() c_coefficient = c.coefficient() if(c_coefficient <= 0.999): print(c_coefficient) if data_type == 'F': assert c_coefficient > 0.999 elif data_type == 'F^2': assert c_coefficient > 0.9 def run_call_back(flags, space_group_info): data_type = ('F', 'F^2')[hasattr(flags, 'F_sq')] options = ( ( True,False,False,False), (False, True,False,False), (False,False, True,False), (False,False,False, True), (False, True, True, True)) for target_functor in xray.target_functors.registry().values(): for (fsite, fu_iso, foccupancy, fu_aniso) in options: gradient_flags = xray.structure_factors.gradient_flags( site = fsite, u_iso = fu_iso, u_aniso = fu_aniso, occupancy = foccupancy) for anomalous_flag in (False, True)[:]: #SWITCH u_penalty_types = [None] tan_u_isos = [False] if (gradient_flags.u_iso): tan_u_isos.append(True) occupancy_penalty_types = [None] if (gradient_flags.occupancy): occupancy_penalty_types.append( xray.minimization.occupancy_penalty_exp()) for tan_u_iso in tan_u_isos: if(tan_u_iso): param = 100 else: param = 0 for occupancy_penalty in occupancy_penalty_types: if(0): print(fsite,fu_iso,foccupancy,fu_aniso,anomalous_flag,tan_u_iso) do_exercise = lambda: exercise( target_functor, data_type, space_group_info, anomalous_flag, gradient_flags, occupancy_penalty=occupancy_penalty, verbose=flags.Verbose, tan_u_iso=tan_u_iso, param = param, d_min = 2.5) try: do_exercise() except AssertionError: print("Test did not pass: ruling out a random fluke...") do_exercise() do_exercise() print("Ruled out!") def run(): cmd_args = ['--F', '--F_sq', '--debugging'] extra = [] args = [] for arg in sys.argv[1:]: if arg in cmd_args: extra.append(arg[2:]) else: args.append(arg) extra = tuple(extra) debug = () if 'debugging' in extra: extra = () if 0: extra += ('F_sq',) assert not('F' in extra and 'F_sq' in extra) if 'F' in extra: print('Refinement against F') if 'F_sq' in extra: print('Refinement against F^2') debug_utils.parse_options_loop_space_groups( args, run_call_back, keywords=extra+debug) if (__name__ == "__main__"): run()