from __future__ import absolute_import, division, print_function from cctbx import miller from cctbx import maptbx import mmtbx.masks class target_map(object): """ Data structure for real-space refinement. Holds map and its reciprocal-space equivalent, and its masked version that can be updated as model changes. """ def __init__(self, map_data, xray_structure, d_min, atom_radius, map_manager=None): self.d_min = d_min self.map_data = map_data self.map_manager = map_manager # XXX switching to map_manager from map_data self.atom_radius = atom_radius self.f_map_diff = None # XXX rudimentary left-over, remove later self.crystal_gridding = maptbx.crystal_gridding( #XXX Likewise, remove later unit_cell = xray_structure.unit_cell(), pre_determined_n_real = map_data.all(), space_group_info = xray_structure.space_group_info()) self.miller_array = None while d_min %.6g" % ( # rotamer_id, refined.rs_f_start, refined.rs_f_final) # return refined # def refine_restrained(): # refined = residue_refine_restrained( # pdb_hierarchy=pdb_hierarchy, # residue=residue, # density_map=density_map, # geometry_restraints_manager=geometry_restraints_manager, # real_space_target_weight=real_space_target_weight, # real_space_gradients_delta=real_space_gradients_delta, # lbfgs_termination_params=lbfgs_termination_params) # print residue.id_str(), "restr. refined(%s): %.6g -> %.6g" % ( # rotamer_id, refined.rs_f_start, refined.rs_f_final) # return refined # def refine(): # residue.atoms().set_xyz(new_xyz=refine_constrained().sites_cart_residue) # return refine_restrained() # for model in pdb_hierarchy.models(): # for chain in model.chains(): # for residue in chain.only_conformer().residues(): # if (get_class(residue.resname) != "common_amino_acid"): # n_other_residues += 1 # else: # rotamer_iterator = mmtbx.rotamer.iterator( # mon_lib_srv=mon_lib_srv, # residue=residue, # atom_selection_bool=atom_selection_bool) # if (rotamer_iterator is None): # n_amino_acids_ignored += 1 # else: # best = None # rotamer_id = "as_given" # if not force_rotamer: # best = group_args(rotamer_id=rotamer_id, refined=refine()) # n_amino_acids_scored += 1 # for rotamer,rotamer_sites_cart in rotamer_iterator: # residue.atoms().set_xyz(new_xyz=rotamer_sites_cart) # rotamer_id=rotamer.id # trial = group_args(rotamer_id=rotamer.id, refined=refine()) # if best is None: # best = trial # if (trial.refined.rs_f_final > best.refined.rs_f_final): # best = trial # print residue.id_str(), "best rotamer:", best.rotamer_id # residue.atoms().set_xyz(new_xyz=best.refined.sites_cart_residue) # print # # print "number of amino acid residues scored:", n_amino_acids_scored # print "number of amino acid residues ignored:", n_amino_acids_ignored # print "number of other residues:", n_other_residues # print # sys.stdout.flush() # #def get_best_rotamer(processed_pdb_file, # fft_map, # d_max, # d_min): # master_phil = get_master_phil() # work_params = master_phil.extract() # get_class = iotbx.pdb.common_residue_names_get_class # mon_lib_srv = mmtbx.monomer_library.server.server() # grm = processed_pdb_file.geometry_restraints_manager( # params_edits=work_params.geometry_restraints.edits, # params_remove=work_params.geometry_restraints.remove) # density_map = fft_map.real_map() # atom_selection_bool=get_atom_selection_bool( # scope_extract=work_params.rotamer_score_and_choose_best, # attr="atom_selection", # processed_pdb_file=processed_pdb_file, # work_params=work_params) # n_other_residues = 0 # n_amino_acids_ignored = 0 # n_amino_acids_scored = 0 # pdb_hierarchy = processed_pdb_file.all_chain_proxies.pdb_hierarchy # for model in pdb_hierarchy.models(): # for chain in model.chains(): # for residue in chain.only_conformer().residues(): # if (get_class(residue.resname) != "common_amino_acid"): # n_other_residues += 1 # else: # rotamer_iterator = mmtbx.rotamer.iterator( # mon_lib_srv=mon_lib_srv, # residue=residue, # atom_selection_bool=atom_selection_bool) # if (rotamer_iterator is None): # n_amino_acids_ignored += 1 # else: # best = None # best_sites_cart = None # n_amino_acids_scored += 1 # for rotamer,rotamer_sites_cart in rotamer_iterator: # residue.atoms().set_xyz(new_xyz=rotamer_sites_cart) # trial = group_args(rotamer_id=rotamer.id, # rs_f=compute_functional_lite( # pdb_hierarchy=pdb_hierarchy, # residue=residue, # density_map=density_map, # geometry_restraints_manager=grm, # real_space_target_weight=work_params.real_space_target_weight)) # #print residue.id_str(),trial.rotamer_id, trial.rs_f # if best is None: # best = trial # best_sites_cart = rotamer_sites_cart # if (trial.rs_f > best.rs_f): # best = trial # best_sites_cart = rotamer_sites_cart # print residue.id_str(), "best rotamer:", best.rotamer_id # residue.atoms().set_xyz(new_xyz=best_sites_cart) # #print # #real_space_params_phil_str = """\ #real_space_target_weight = 10 # .type = float #real_space_gradients_delta_resolution_factor = 1/3 # .type = float #""" # #coordinate_refinement_export_phil_str = """\ #finishing_geometry_minimization #{ # cycles_max = 100 # .type = int # first_weight_scale = 0.1 # .type = float # cycle_weight_multiplier = 1.0 # .type = float # superpose_cycle_end_with_cycle_start = False # .type = bool # dihedral_restraints = False # .type = bool # output_file = Auto # .type = str #} #real_space_target_weights { # first_sample = 10 # .type = float # sampling_step = 30 # .type = float # number_of_samples = 10 # .type = int # bond_rmsd_target = 0.03 # .type = float # worst_acceptable_bond_rmsd { # pool_size = 10 # .type = int # max_pool_average = 0.1 # .type = float # } #} #lbfgs_max_iterations = 500 # .type = int #""" # #coordinate_refinement_params_phil_str = """\ #coordinate_refinement { # run = False # .type = bool # atom_selection = Auto # .type = str # compute_final_correlation = True # .type = bool # home_restraints # .multiple = True # { # selection = None # .type = str # sigma = 0.05 # .type = float # slack = 0 # .type = float # } # %(coordinate_refinement_export_phil_str)s #} #""" % vars() # #def get_master_phil(): # return iotbx.phil.parse( # input_string="""\ #atom_selection = None # .type = str #strict_processing = True # .type = bool # #map { # coeff_labels { # f = 2FOFCWT,PH2FOFCWT # .type = str # phases = None # .type = str # weights = None # .type = str # } # std_dev_weight_coeff_labels { # f = None # .type = str # phases = None # .type = str # weights = None # .type = str # } # low_resolution = None # .type = float # high_resolution = None # .type = float # grid_resolution_factor = 1/3 # .type = float #} # #output_file = Auto # .type = str #final_geo_file = Auto # .type = str # #%s #local_standard_deviations_radius = None # .type = float #weight_map_normalization = *fourier sigma # .type = choice #weight_map_scale_factor = None # .type = float # #%s # #rotamer_score_and_choose_best { # run = False # .type = bool # atom_selection = Auto # .type = str # lbfgs_max_iterations = 50 # .type = int # output_file = Auto # .type = str #} # #include scope mmtbx.monomer_library.pdb_interpretation.grand_master_phil_str #""" % (real_space_params_phil_str, coordinate_refinement_params_phil_str), # process_includes=True) # #def extract_map_coeffs(params, optional, miller_arrays): # def find(labels): # for miller_array in miller_arrays: # if (",".join(miller_array.info().labels) == labels): # return miller_array # matching_array = None # for miller_array in miller_arrays: # if (",".join(miller_array.info().labels).lower() == labels.lower()): # if (matching_array is not None): # return None # matching_array = miller_array # return matching_array # def raise_sorry(msg_intro, name): # msg = [ # msg_intro, # " %s = %s" % params.__phil_path_and_value__(object_name=name), # " List of available labels:"] # for miller_array in miller_arrays: # msg.append(" %s" % ",".join(miller_array.info().labels)) # raise Sorry("\n".join(msg)) # if (params.f is None): # if (not optional): # raise_sorry(msg_intro="Missing assignment:", name="f") # return None # f = find(labels=params.f) # if (f is None): # raise_sorry(msg_intro="Cannot find map coefficients:", name="f") # if (not f.is_complex_array() or params.phases is not None): # if (params.phases is None): # raise_sorry(msg_intro="Missing assignment:", name="phases") # phases = find(labels=params.phases) # if (phases is None): # raise_sorry( # msg_intro="Cannot find map coefficient phases:", name="phases") # if (phases.is_hendrickson_lattman_array() # and phases.data().all_eq((0,0,0,0))): # raise Sorry( # "All Hendrickson-Lattman coefficient zero:\n" # " %s = %s" % params.__phil_path_and_value__(object_name="phases")) # cf, cp = f.common_sets(other=phases) # if (cf.indices().size() != f.indices().size()): # raise Sorry( # "Number of missing map coefficient phases: %d" % ( # f.indices().size() - cf.indices().size())) # f = cf.phase_transfer(phase_source=cp, deg=True) # if (params.weights is not None): # weights = find(labels=params.weights) # if (weights is None): # raise_sorry( # msg_intro="Cannot find map coefficient weights:", # name="weights") # cf, cw = f.common_sets(other=weights) # if (cf.indices().size() != f.indices().size()): # raise Sorry( # "Number of missing map coefficient weights: %d" % ( # f.indices().size() - cf.indices().size())) # f = cf.customized_copy(data=cw.data()*cf.data()) # return f # #class home_restraints(object): # # __slots__ = ["iselection", "weight", "slack"] # # def __init__(O, iselection, weight, slack): # O.iselection = iselection # O.weight = weight # O.slack = slack # #def process_home_restraints_params(work_params, processed_pdb_file): # result = [] # for params in work_params: # sigma = params.sigma # if (sigma is None or sigma <= 0): # continue # bsel = processed_pdb_file.all_chain_proxies.phil_atom_selection( # cache=None, # scope_extract=params, # attr="selection", # allow_none=True, # allow_auto=False) # if (bsel is not None): # slack = params.slack # if (slack is None or slack <= 0): # slack = 0 # result.append(home_restraints( # iselection=bsel.iselection(), weight=1/sigma**2, slack=slack)) # return result # #class geometry_restraints_manager_plus(object): # # __slots__ = [ # "manager", # "home_sites_cart", # "home_restraints_list", # "home_dihedral_proxies", # "crystal_symmetry", # "pair_proxies", # "angle_proxies", # "dihedral_proxies", # "site_symmetry_table", # "plain_pair_sym_table", # "plain_pairs_radius"] # # def __init__(O, # manager, # home_sites_cart=None, # home_restraints_list=None, # home_dihedral_proxies=None): # O.manager = manager # #O.sites_cart = sites_cart # O.home_sites_cart = home_sites_cart # O.home_restraints_list = home_restraints_list # O.home_dihedral_proxies = home_dihedral_proxies # O.crystal_symmetry = manager.crystal_symmetry # O.pair_proxies = manager.pair_proxies # O.angle_proxies = manager.angle_proxies # O.dihedral_proxies = manager.dihedral_proxies # O.site_symmetry_table = manager.site_symmetry_table # O.plain_pair_sym_table = manager.plain_pair_sym_table # O.plain_pairs_radius = manager.plain_pairs_radius # # def energies_add(O, energies_obj): # if (O.manager.site_symmetry_table is not None): # site_symmetry_table_indices = O.manager.site_symmetry_table.indices() # else: # site_symmetry_table_indices = None # from cctbx.geometry_restraints import bond # n_restraints = 0 # r_sum = 0 # hr_summation = cctbx.geometry_restraints \ # .home_restraints_summation_skip_special_positions # for hr in O.home_restraints_list: # r_sum += hr_summation( # sites_cart=energies_obj.sites_cart, # gradients=energies_obj.gradients, # site_symmetry_table_indices=site_symmetry_table_indices, # home_sites_cart=O.home_sites_cart, # iselection=hr.iselection, # weight=hr.weight, # slack=hr.slack) # n_restraints += hr.iselection.size() # energies_obj.n_home_restraints = n_restraints # energies_obj.home_restraints_residual_sum = r_sum # energies_obj.number_of_restraints += n_restraints # energies_obj.residual_sum += r_sum # r_sum=cctbx.geometry_restraints.dihedral_residual_sum( # sites_cart=energies_obj.sites_cart, # proxies=O.home_dihedral_proxies, # gradient_array=energies_obj.gradients) # if O.home_dihedral_proxies is not None: # n_restraints = len(O.home_dihedral_proxies) # else: # n_restraints = 0 # energies_obj.n_dihedral_restraints = n_restraints # energies_obj.dihedral_restraints_residual_sum = r_sum # energies_obj.residual_sum += r_sum # energies_obj.number_of_restraints += n_restraints # # def energies_show(O, energies_obj, f, prefix): # print >> f, prefix+" home_restraints_residual_sum (n=%d): %.6g" % ( # energies_obj.n_home_restraints, # energies_obj.home_restraints_residual_sum) # # def dihedral_energies_show(O, energies_obj, f, prefix): # print >> f, prefix+" dihedral_restraints_residual_sum (n=%d): %.6g" % ( # energies_obj.n_dihedral_restraints, # energies_obj.dihedral_restraints_residual_sum) # # def energies_sites(O, # sites_cart, # flags=None, # custom_nonbonded_function=None, # compute_gradients=False, # gradients=None, # disable_asu_cache=False, # normalization=False): # return O.manager.energies_sites( # sites_cart=sites_cart, # flags=flags, # custom_nonbonded_function=custom_nonbonded_function, # compute_gradients=compute_gradients, # gradients=gradients, # disable_asu_cache=disable_asu_cache, # normalization=normalization, # extension_objects=[O]) # # def show_sorted(O, # flags=None, # sites_cart=None, # site_labels=None, # f=None): # return O.manager.show_sorted( # flags=flags, # sites_cart=sites_cart, # site_labels=site_labels, # f=f) # # def select(O, selection=None, iselection=None): # return O.manager.select(selection=selection, iselection=iselection) # #def fit_rotamers_simple(processed_pdb_file, # fft_map, # d_max, # d_min): # master_phil = get_master_phil() # work_params = master_phil.extract() # mon_lib_srv = mmtbx.monomer_library.server.server() # #ener_lib = mmtbx.monomer_library.server.ener_lib() # grm = processed_pdb_file.geometry_restraints_manager( # params_edits=work_params.geometry_restraints.edits, # params_remove=work_params.geometry_restraints.remove) # density_map = fft_map.real_map() # real_space_gradients_delta = \ # d_min * work_params.real_space_gradients_delta_resolution_factor # rotamer_score_and_choose_best( # mon_lib_srv=mon_lib_srv, # density_map=density_map, # pdb_hierarchy=processed_pdb_file.all_chain_proxies.pdb_hierarchy, # geometry_restraints_manager=grm, # atom_selection_bool=get_atom_selection_bool( # scope_extract=work_params.rotamer_score_and_choose_best, # attr="atom_selection", # processed_pdb_file=processed_pdb_file, # work_params=work_params), # real_space_target_weight=work_params.real_space_target_weight, # real_space_gradients_delta=real_space_gradients_delta, # lbfgs_termination_params=scitbx.lbfgs.termination_parameters( # max_iterations=work_params # .rotamer_score_and_choose_best.lbfgs_max_iterations), # force_rotamer=True) # #def get_atom_selection_bool(scope_extract, # attr, # processed_pdb_file, # work_params): # common_atom_selection_bool_cache=[] # result = processed_pdb_file.all_chain_proxies \ # .phil_atom_selection( # cache=None, # scope_extract=scope_extract, # attr="atom_selection", # allow_none=True, # allow_auto=True) # if (result is None or result is not Auto): # return result # if (len(common_atom_selection_bool_cache) == 0): # common_atom_selection_bool_cache.append( # processed_pdb_file.all_chain_proxies # .phil_atom_selection( # cache=None, # scope_extract=work_params, # attr="atom_selection", # allow_none=True)) # return common_atom_selection_bool_cache[0] # #def compose_output_file_name(input_pdb_file_name, new_suffix): # result = op.basename(input_pdb_file_name) # if ( result.endswith(".pdb") # or result.endswith(".ent")): # result = result[:-4] # result += new_suffix # return result # #def write_pdb( # file_name, # input_pdb_file_name, # processed_pdb_file, # grm, # new_suffix): # if (file_name is not None): # if (file_name is Auto): # file_name = compose_output_file_name( # input_pdb_file_name=input_pdb_file_name, # new_suffix=new_suffix) # pdb_hierarchy = processed_pdb_file.all_chain_proxies.pdb_hierarchy # print "Writing file: %s" % show_string(file_name) # sys.stdout.flush() # pdb_hierarchy.write_pdb_file( # file_name=file_name, # crystal_symmetry=grm.crystal_symmetry) # print # #def run_coordinate_refinement_driver( # processed_pdb_file, # geometry_restraints_manager, # fft_map, # real_space_gradients_delta, # work_params, # std_dev_weight_fft_map=None, # write_pdb_callback=None, # log=None): # if (log is None): log = null_out() # atom_selection_bool = get_atom_selection_bool( # scope_extract=work_params.coordinate_refinement, # attr="atom_selection", # processed_pdb_file=processed_pdb_file, # work_params=work_params) # home_restraints_list = process_home_restraints_params( # work_params=work_params.coordinate_refinement.home_restraints, # processed_pdb_file=processed_pdb_file) # if (not work_params.coordinate_refinement.compute_final_correlation): # work_scatterers = None # else: # work_scatterers = processed_pdb_file.xray_structure( # show_summary=False).scatterers() # if (atom_selection_bool is None): # work_scatterers = work_scatterers.deep_copy() # else: # work_scatterers = work_scatterers.select(atom_selection_bool) # if (std_dev_weight_fft_map is None): # weight_map = None # else: # assert std_dev_weight_fft_map.n_real() == fft_map.n_real() # weight_map = std_dev_weight_fft_map.real_map() # return run_coordinate_refinement( # pdb_atoms=processed_pdb_file.all_chain_proxies.pdb_atoms, # geometry_restraints_manager=geometry_restraints_manager, # selection_variable=atom_selection_bool, # density_map=fft_map.real_map(), # real_space_gradients_delta=real_space_gradients_delta, # work_params=work_params, # home_restraints_list=home_restraints_list, # work_scatterers=work_scatterers, # unit_cell=fft_map.unit_cell(), # d_min=fft_map.d_min(), # weight_map=weight_map, # write_pdb_callback=write_pdb_callback, # log=log) # #def run_coordinate_refinement( # geometry_restraints_manager, # selection_variable, # density_map, # real_space_gradients_delta, # work_params, # pdb_atoms=None, # sites_cart=None, # home_restraints_list=[], # work_scatterers=None, # unit_cell=None, # d_min=None, # weight_map=None, # write_pdb_callback=None, # log=None): # assert [pdb_atoms, sites_cart].count(None)==1 # if (log is None): log = null_out() # if (work_scatterers is not None): # assert unit_cell is not None # assert d_min is not None # best_info = None # if (pdb_atoms is not None): # sites_cart_start = pdb_atoms.extract_xyz() # site_labels = [atom.id_str() for atom in pdb_atoms] # else: # sites_cart_start = sites_cart # site_labels = None # grmp = geometry_restraints_manager_plus( # manager=geometry_restraints_manager, # home_sites_cart=sites_cart_start, # home_restraints_list=home_restraints_list) # print >> log, "Before coordinate refinement:" # grmp.energies_sites(sites_cart=sites_cart_start).show(f=log) # print >> log # log.flush() # rstw_params = work_params.coordinate_refinement.real_space_target_weights # if (rstw_params.number_of_samples is None): # rstw_list = [work_params.real_space_target_weight] # else: # rstw_list = [rstw_params.first_sample + i * rstw_params.sampling_step # for i in range(rstw_params.number_of_samples)] # lbfgs_termination_params = scitbx.lbfgs.termination_parameters( # max_iterations=work_params.coordinate_refinement # .lbfgs_max_iterations) # lbfgs_exception_handling_params = \ # scitbx.lbfgs.exception_handling_parameters( # ignore_line_search_failed_step_at_lower_bound=True) # bond_rmsd_list = [] # for rstw in rstw_list: # refined = maptbx.real_space_refinement_simple.lbfgs( # sites_cart=sites_cart_start, # density_map=density_map, # weight_map=weight_map, # selection_variable=selection_variable, # geometry_restraints_manager=grmp, # real_space_target_weight=rstw, # real_space_gradients_delta=real_space_gradients_delta, # local_standard_deviations_radius= # work_params.local_standard_deviations_radius, # weight_map_scale_factor=work_params.weight_map_scale_factor, # lbfgs_termination_params=lbfgs_termination_params, # lbfgs_exception_handling_params=lbfgs_exception_handling_params) # print >> log, "After coordinate refinement" \ # " with real-space target weight %.1f:" % rstw # grmp.energies_sites(sites_cart=refined.sites_cart).show(f=log) # bond_proxies = grmp.pair_proxies().bond_proxies # bond_proxies.show_sorted( # by_value="residual", # sites_cart=refined.sites_cart, # site_labels=site_labels, # f=log, # prefix=" ", # max_items=3) # print >> log # print >> log, " number_of_function_evaluations:", \ # refined.number_of_function_evaluations # print >> log, " real+geo target start: %.6g" % refined.f_start # print >> log, " real+geo target final: %.6g" % refined.f_final # deltas = bond_proxies.deltas(sites_cart=refined.sites_cart) # deltas_abs = flex.abs(deltas) # deltas_abs_sorted = deltas_abs.select( # flex.sort_permutation(deltas_abs), reverse=True) # wabr_params = rstw_params.worst_acceptable_bond_rmsd # acceptable = ( # flex.mean_default(deltas_abs_sorted[:wabr_params.pool_size],0) # < wabr_params.max_pool_average) # if(deltas.size()>0): # bond_rmsd = flex.mean_sq(deltas)**0.5 # else: bond_rmsd=0 # bond_rmsd_list.append(bond_rmsd) # print >> log, " Bond RMSD: %.3f" % bond_rmsd # # # if (work_scatterers is None): # region_cc = None # else: # region_cc = maptbx.region_density_correlation( # large_unit_cell=unit_cell, # large_d_min=d_min, # large_density_map=density_map, # sites_cart=refined.sites_cart_variable, # site_radii=flex.double(refined.sites_cart_variable.size(), 1), # work_scatterers=work_scatterers) # # # rmsd_diff = abs(rstw_params.bond_rmsd_target - bond_rmsd) # if ( best_info is None # or best_info.rmsd_diff > rmsd_diff # and (acceptable or not best_info.acceptable)): # best_info = group_args( # rstw=rstw, # refined=refined, # acceptable=acceptable, # rmsd_diff=rmsd_diff, # region_cc=region_cc, # sites_cart=refined.sites_cart) # print >> log # if (best_info is not None): # print >> log, "Table of real-space target weights vs. bond RMSD:" # print >> log, " weight RMSD" # for w,d in zip(rstw_list, bond_rmsd_list): # print >> log, " %6.1f %5.3f" % (w,d) # print >> log, "Best real-space target weight: %.1f" % best_info.rstw # print >> log, \ # "Associated refined final value: %.6g" % best_info.refined.f_final # if (best_info.region_cc is not None): # print >> log, "Associated region correlation: %.4f" % best_info.region_cc # print >> log # if (pdb_atoms is not None): # pdb_atoms.set_xyz(new_xyz=refined.sites_cart) # if (write_pdb_callback is not None): # write_pdb_callback(situation="after_best_weight_determination") # # # fgm_params = work_params.coordinate_refinement \ # .finishing_geometry_minimization # if (fgm_params.cycles_max is not None and fgm_params.cycles_max > 0): # print >> log, "As previously obtained with target weight %.1f:" \ # % best_info.rstw # grmp.energies_sites(sites_cart=refined.sites_cart).show( # f=log, prefix=" ") # print >> log, "Finishing refinement to idealize geometry:" # print >> log, " number of function" # print >> log, " weight evaluations cycle RMSD" # number_of_fgm_cycles = 0 # rstw = best_info.rstw * fgm_params.first_weight_scale # sites_cart_start = best_info.refined.sites_cart.deep_copy() # for i_cycle in range(fgm_params.cycles_max): # fgm_refined = maptbx.real_space_refinement_simple.lbfgs( # sites_cart=sites_cart_start, # density_map=density_map, # selection_variable=selection_variable, # geometry_restraints_manager=grmp, # energies_sites_flags=cctbx.geometry_restraints.flags.flags( # default=True, dihedral=fgm_params.dihedral_restraints), # real_space_target_weight=rstw, # real_space_gradients_delta=real_space_gradients_delta, # local_standard_deviations_radius= # work_params.local_standard_deviations_radius, # weight_map_scale_factor=work_params.weight_map_scale_factor, # lbfgs_termination_params=lbfgs_termination_params, # lbfgs_exception_handling_params=lbfgs_exception_handling_params) # cycle_rmsd = sites_cart_start.rms_difference(fgm_refined.sites_cart) # print >> log, " %6.1f %10d %6.3f" % ( # rstw, fgm_refined.number_of_function_evaluations, cycle_rmsd) # number_of_fgm_cycles += 1 # rstw *= fgm_params.cycle_weight_multiplier # if (cycle_rmsd < 1.e-4): # break # if (fgm_params.superpose_cycle_end_with_cycle_start): # fit = scitbx.math.superpose.least_squares_fit( # reference_sites=best_info.refined.sites_cart, # other_sites=fgm_refined.sites_cart) # fgm_refined.sites_cart = fit.other_sites_best_fit() # sites_cart_start = fgm_refined.sites_cart.deep_copy() # print >> log, "After %d refinements to idealize geometry:" % ( # number_of_fgm_cycles) # grmp.energies_sites(sites_cart=fgm_refined.sites_cart).show( # f=log, prefix=" ") # if (work_scatterers is None): # fgm_region_cc = None # else: # fgm_region_cc = maptbx.region_density_correlation( # large_unit_cell=unit_cell, # large_d_min=d_min, # large_density_map=density_map, # sites_cart=fgm_refined.sites_cart_variable, # site_radii=flex.double(fgm_refined.sites_cart_variable.size(), 1), # work_scatterers=work_scatterers) # print >> log, " Associated region correlation: %.4f" % fgm_region_cc # print >> log # best_info.fgm_refined = fgm_refined # best_info.fgm_region_cc = fgm_region_cc # if (pdb_atoms is not None): # pdb_atoms.set_xyz(new_xyz=fgm_refined.sites_cart) # best_info.sites_cart = fgm_refined.sites_cart # if (write_pdb_callback is not None): # write_pdb_callback(situation="after_finishing_geo_min") # else: # best_info.fgm_refined = None # best_info.fgm_region_cc = None # return best_info # #def run(args,log=None): # if log is None: # log = sys.stdout # # show_times = libtbx.utils.show_times(time_start="now") # master_phil = get_master_phil() # import iotbx.utils # input_objects = iotbx.utils.process_command_line_inputs( # args=args, # master_phil=master_phil, # input_types=("mtz", "pdb", "cif")) # work_phil = master_phil.fetch(sources=input_objects["phil"]) # work_phil.show() # print >> log # print >> log, "#phil __OFF__" # print >> log # work_params = work_phil.extract() # # # assert len(input_objects["mtz"]) == 1 # miller_arrays = input_objects["mtz"][0].file_content.as_miller_arrays() # map_coeffs = extract_map_coeffs( # miller_arrays=miller_arrays, # params=work_params.map.coeff_labels, # optional=False) # std_dev_weight_coeffs = extract_map_coeffs( # miller_arrays=miller_arrays, # params=work_params.map.std_dev_weight_coeff_labels, # optional=True) # del miller_arrays # # # mon_lib_srv = mmtbx.monomer_library.server.server() # ener_lib = mmtbx.monomer_library.server.ener_lib() # for file_obj in input_objects["cif"]: # print >> log, "Processing CIF file: %s" % show_string(file_obj.file_name) # for srv in [mon_lib_srv, ener_lib]: # srv.process_cif_object( # cif_object=file_obj.file_content, # file_name=file_obj.file_name) # if (len(input_objects["cif"]) != 0): # print >> log # # # assert len(input_objects["pdb"]) == 1 # TODO not implemented # file_obj = input_objects["pdb"][0] # input_pdb_file_name = file_obj.file_name # print >> log, "Crystal symmetry (from map file):" # map_coeffs.crystal_symmetry().show_summary(prefix=" ") # pdb_crystal_symmetry = file_obj.file_content.crystal_symmetry() # if ( pdb_crystal_symmetry is not None # and not pdb_crystal_symmetry.is_similar_symmetry( # map_coeffs.crystal_symmetry())): # print >>log, " NOTE: Crystal symmetry from PDB file is different:" # pdb_crystal_symmetry.show_summary(prefix=" Ignored: ") # print >> log # processed_pdb_file = mmtbx.monomer_library.pdb_interpretation.process( # mon_lib_srv=mon_lib_srv, # ener_lib=ener_lib, # params=work_params.pdb_interpretation, # file_name=file_obj.file_name, # pdb_inp=file_obj.file_content, # strict_conflict_handling=work_params.strict_processing, # substitute_non_crystallographic_unit_cell_if_necessary=True, # crystal_symmetry=map_coeffs.crystal_symmetry(), # force_symmetry=True, # log=sys.stdout) # if (work_params.strict_processing): # msg = processed_pdb_file.all_chain_proxies.fatal_problems_message() # if (msg is not None): # raise Sorry(msg) # # # grm = processed_pdb_file.geometry_restraints_manager( # params_edits=work_params.geometry_restraints.edits, # params_remove=work_params.geometry_restraints.remove) # print >> log # sys.stdout.flush() # # # def show_completeness(annotation): # print >>log, "Completeness of %s map coefficients:" % annotation # map_coeffs.setup_binner(auto_binning=True) # if (map_coeffs.binner().n_bins_used() > 12): # map_coeffs.setup_binner(n_bins=12) # map_coeffs.completeness(use_binning=True).show(prefix=" ") # print >> log # sys.stdout.flush() # show_completeness("input") # map_coeffs_input = map_coeffs # # # low_res = work_params.map.low_resolution # high_res = work_params.map.high_resolution # d_max, d_min = map_coeffs.d_max_min() # d_max_apply, d_min_apply = None, None # if (low_res is not None and low_res < d_max): # d_max_apply = low_res # print >> log , "Applying low resolution cutoff to map coefficients:" \ # " d_max=%.6g" % d_max_apply # if (high_res is not None and high_res > d_min): # d_min_apply = high_res # print >> log, "Applying high resolution cutoff to map coefficients:" \ # " d_min=%.6g" % d_min_apply # if (d_max_apply is not None or d_min_apply is not None): # map_coeffs = map_coeffs.resolution_filter( # d_max=d_max_apply, d_min=d_min_apply) # if (std_dev_weight_coeffs is not None): # std_dev_weight_coeffs = std_dev_weight_coeffs.resolution_filter( # d_max=d_max_apply, d_min=d_min_apply) # if (d_min_apply is not None): # d_min = d_min_apply # print >> log # sys.stdout.flush() # # # if (map_coeffs is not map_coeffs_input): # show_completeness("final") # # # fft_map = map_coeffs.fft_map( # d_min=d_min, # resolution_factor=work_params.map.grid_resolution_factor) # fft_map.apply_sigma_scaling() # std_dev_weight_fft_map = None # if (std_dev_weight_coeffs is not None): # std_dev_weight_fft_map = miller.fft_map( # crystal_gridding=fft_map, # fourier_coefficients=std_dev_weight_coeffs) # if (work_params.weight_map_normalization == "fourier"): # std_dev_weight_fft_map.apply_fourier_scaling() # elif (work_params.weight_map_normalization == "sigma"): # std_dev_weight_fft_map.apply_sigma_scaling() # else: # assert work_params.weight_map_normalization is None # # # real_space_gradients_delta = \ # d_min * work_params.real_space_gradients_delta_resolution_factor # print >> log, "real_space_gradients_delta: %.6g" % real_space_gradients_delta # print >> log # sys.stdout.flush() # # # if (work_params.rotamer_score_and_choose_best.run): # rotamer_score_and_choose_best( # mon_lib_srv=mon_lib_srv, # density_map=fft_map.real_map(), # pdb_hierarchy=processed_pdb_file.all_chain_proxies.pdb_hierarchy, # geometry_restraints_manager=grm, # atom_selection_bool=get_atom_selection_bool( # scope_extract=work_params.rotamer_score_and_choose_best, # attr="atom_selection", # processed_pdb_file=processed_pdb_file, # work_params=work_params), # real_space_target_weight=work_params.real_space_target_weight, # real_space_gradients_delta=real_space_gradients_delta, # lbfgs_termination_params=scitbx.lbfgs.termination_parameters( # max_iterations=work_params # .rotamer_score_and_choose_best.lbfgs_max_iterations)) # write_pdb( # file_name=work_params.rotamer_score_and_choose_best.output_file, # input_pdb_file_name=input_pdb_file_name, # processed_pdb_file=processed_pdb_file, # grm=grm, # new_suffix="_lockit_rotamer_score_and_choose_best.pdb") # # # best_info = None # if (work_params.coordinate_refinement.run): # def write_pdb_callback(situation): # if (situation == "after_best_weight_determination"): # write_pdb( # file_name=work_params.output_file, # input_pdb_file_name=input_pdb_file_name, # processed_pdb_file=processed_pdb_file, # grm=grm, # new_suffix="_lockit_best_weight.pdb") # elif (situation == "after_finishing_geo_min"): # write_pdb( # file_name=work_params.coordinate_refinement \ # .finishing_geometry_minimization.output_file, # input_pdb_file_name=input_pdb_file_name, # processed_pdb_file=processed_pdb_file, # grm=grm, # new_suffix="_lockit_finishing_geo_min.pdb") # else: # raise AssertionError # best_info = run_coordinate_refinement_driver( # processed_pdb_file=processed_pdb_file, # geometry_restraints_manager=grm, # fft_map=fft_map, # std_dev_weight_fft_map=std_dev_weight_fft_map, # real_space_gradients_delta=real_space_gradients_delta, # work_params=work_params, # write_pdb_callback=write_pdb_callback, # log=sys.stdout) # else: # write_pdb( # file_name=work_params.output_file, # input_pdb_file_name=input_pdb_file_name, # processed_pdb_file=processed_pdb_file, # grm=grm, # new_suffix="_lockit_best_weight.pdb") # # # file_name = work_params.final_geo_file # if (file_name is not None): # if (file_name is Auto): # file_name = compose_output_file_name( # input_pdb_file_name=input_pdb_file_name, # new_suffix="_lockit_final.geo") # pdb_atoms = processed_pdb_file.all_chain_proxies.pdb_atoms # sites_cart = pdb_atoms.extract_xyz() # site_labels = [atom.id_str() for atom in pdb_atoms] # print >> log, "Writing file: %s" % show_string(file_name) # sys.stdout.flush() # f = open(file_name, "w") # grm.show_sorted(sites_cart=sites_cart, site_labels=site_labels, f=f) # del f # print >> log # # # show_times() # sys.stdout.flush() # return best_info # #if (__name__ == "__main__"): # run(args=sys.argv[1:],log=None) #