text: 
  from __future__ import division
  import mmtbx.solvent.ensemble_ordered_solvent as ensemble_ordered_solvent
  from mmtbx.refinement.ensemble_refinement import ensemble_utils
  from mmtbx.dynamics import ensemble_cd
  import mmtbx.tls.tools as tls_tools
  import mmtbx.command_line
  import mmtbx.utils
  import mmtbx.model
  import mmtbx.maps
  from iotbx.option_parser import iotbx_option_parser
  from iotbx import pdb
  import iotbx.phil
  import iotbx
  from cctbx import geometry_restraints
  from cctbx.array_family import flex
  from cctbx import miller
  from cctbx import adptbx
  from cctbx import xray
  import scitbx.math
  from libtbx.utils import Sorry, user_plus_sys_time, multi_out, show_total_time
  from libtbx import adopt_init_args, slots_getstate_setstate
  from libtbx.str_utils import format_value, make_header
  from libtbx import runtime_utils
  from libtbx import easy_mp
  import libtbx.load_env
  from cStringIO import StringIO
  import cPickle
  import random
  import gzip
  import math
  import time
  import os
  import sys
  
  # these supersede the defaults in included scopes
  customization_params = iotbx.phil.parse("""
  ensemble_refinement.mask.ignore_hydrogens = False
  ensemble_refinement.mask.n_radial_shells = 1
  ensemble_refinement.mask.radial_shell_width = 1.5
  ensemble_refinement.cartesian_dynamics.number_of_steps = 10
  ensemble_refinement.ensemble_ordered_solvent.b_iso_min = 0.0
  ensemble_refinement.ensemble_ordered_solvent.b_iso_max = 100.0
  ensemble_refinement.ensemble_ordered_solvent.find_peaks.map_next_to_model.max_model_peak_dist = 3.0
  ensemble_refinement.ensemble_ordered_solvent.find_peaks.map_next_to_model.use_hydrogens = False
  """)
  
  # the extra fetch() at the end with the customized parameters gives us the
  # actual master phil object
  input_phil = mmtbx.command_line.generate_master_phil_with_inputs(
    phil_string="",
    enable_experimental_phases=True,
    as_phil_string=True)
  master_params = iotbx.phil.parse(input_phil + """
  ensemble_refinement {
    cartesian_dynamics
      .style = menu_item auto_align box
    {
      include scope mmtbx.dynamics.cartesian_dynamics.master_params
      protein_thermostat = True
        .type = bool
        .help = Use protein atoms thermostat
    }
    update_sigmaa_rfree = 0.001
      .type = float
      .help = test function
    ensemble_reduction = True
      .type = bool
      .help = 'Find miminium number of structures to reproduce simulation R-values'
    verbose = -1
      .type = int
    output_file_prefix = None
      .type = str
      .help = 'Prefix for all output files'
  # TODO
  #  write_mmcif_file = False
  #    .type = bool
    gzip_final_model = True
      .type = bool
      .style = hidden
    random_seed = 2679941
      .type = int
      .help = 'Random seed'
    nproc = 1
      .type = int
      .short_caption = Number of processors
      .style = renderer:draw_nproc_widget
    tx = None
      .type = float
      .help = 'Relaxation time (ps)'
      .short_caption = Relaxation time (ps)
    equilibrium_n_tx = 10
      .type = int
      .help = 'Length of equilibration period, n times tx'
      .short_caption = Length of equilibration period
    acquisition_block_n_tx = 2
      .type = int
      .help = 'Length of acquisition block, n times tx'
      .short_caption = Length of acquisition block
    number_of_acquisition_periods = 10
      .type = int
      .help = 'Number of acquisition periods'
    pdb_stored_per_block = 50
      .type = int
      .help = 'Number of model coordinates stored per acquisition block'
      .short_caption = Models stored per acquisition block
    wxray_coupled_tbath = True
      .type = bool
      .help = 'Use temperature control wxray'
      .short_caption = Use temperature control X-ray weight
    wxray_coupled_tbath_offset = 5.0
      .type = float
      .help = 'Temperature offset, increasing offset increases wxray'
      .short_caption = Temperature ofset
    wxray = 1.0
      .type = float
      .help = 'Multiplier for xray weighting; used if wxray_coupled_tbath = Flase'
      .short_caption = X-ray weight
    tls_group_selections = None
      .type = atom_selection
      .multiple = True
      .help = 'TLS groups to use for TLS fitting (TLS details in PDB header not used)'
      .style = use_list
    ptls = 0.80
      .type = floats
      .optional = False
      .help = 'The fraction of atoms to include in TLS fitting'
      .short_caption = Fraction of atoms to include in TLS fitting
      .style = bold
    max_ptls_cycles = 25
      .type = int
      .help = 'Maximum cycles to use in TLS fitting; TLS will stop prior to this if convergence is reached'
      .short_caption = Max. number of cycles of TLS fitting
    isotropic_b_factor_model = False
      .type = bool
      .help = 'Use isotropic B-factor model instead of TLS'
      .short_caption = Use isotropic B-factor model
    pwilson = 0.8
      .type = float
      .help = 'Scale factor for isotropic b-factor model: all atoms = Bwilson * pwilson'
      .short_caption = Scale factor for isotropic b-factor model
    set_occupancies = False
      .type = bool
      .help = 'Set all atoms aoccupancy to 1.0'
      .short_caption = Reset occupancies to 1.0
    target_name = *ml mlhl ls_wunit_k1_fixed ls_wunit_k1
      .type = choice
      .short_caption = Refinement target
      .help = 'Choices for refinement target'
    remove_alt_conf_from_input_pdb = True
      .type = bool
      .help = 'Removes any alternative conformations if present in input PDB model'
    scale_wrt_n_calc_start = False
      .type = bool
      .help = 'Scale <Ncalc> to starting Ncalc'
      .short_caption = Scale <Ncalc> to starting Ncalc
    output_running_kinetic_energy_in_occupancy_column = False
      .type = bool
      .help = 'Output PDB file contains running kinetic energy in place of occupancy'
    ordered_solvent_update = True
      .type = bool
      .help = 'Ordered water molecules automatically updated every nth macro cycle'
    ordered_solvent_update_cycle = 25
      .type = int
      .help = 'Number of macro-cycles per ordered solvent update'
      .short_caption = Solvent update cycles
    harmonic_restraints
      .style = menu_item auto_align box
    {
      selections = None
        .type = atom_selection
        .help = 'Atoms selections to apply harmonic restraints'
      weight = 0.001
        .type = float
        .help = 'Harmonic restraints weight'
      slack = 1.0
        .type = float
        .help = 'Harmonic restraints slack distance'
    }
    electron_density_maps
      .style = menu_item
    {
      apply_default_maps = True
        .type = bool
      include scope mmtbx.maps.map_coeff_params_str
    }
    mask
      .short_caption = Bulk solvent mask
      .style = menu_item auto_align box
    {
      include scope mmtbx.masks.mask_master_params
    }
    ensemble_ordered_solvent
      .style = menu_item auto_align box
    {
      diff_map_cutoff = 2.5
        .type = float
      e_map_cutoff_keep = 0.5
        .type = float
      e_map_cutoff_find = 0.5
        .type = float
      tolerance = 0.9
        .type = float
      ordered_solvent_map_to_model = True
        .type = bool
      include scope mmtbx.solvent.ordered_solvent.output_params_str
      primary_map_type = mFo-DFmodel
        .type=str
      primary_map_cutoff = 3.0
        .type=float
      secondary_map_type = 2mFo-DFmodel
        .type=str
      secondary_map_cutoff_keep = 1.0
        .type=float
      secondary_map_cutoff_find = 1.0
        .type=float
      include scope mmtbx.solvent.ordered_solvent.h_bond_params_str
      include scope mmtbx.solvent.ordered_solvent.adp_occ_params_str
      refine_occupancies = False
        .type = bool
        .help = Refine solvent occupancies.
        .expert_level = 1
      add_hydrogens = False
        .type = bool
        .help = Adds hydrogens to water molecules (except those on special positions)
      refilter = True
        .type = bool
      temperature = 300
        .type = float
        .help = Target temperature for random velocity assignment
      seed = 343534534
        .type = int
        .help = Fixes the random seed for velocity assignment
      preserved_solvent_minimum_distance = 7.0
        .type = float
      find_peaks {
        include scope mmtbx.find_peaks.master_params
      }
    }
    include scope mmtbx.geometry_restraints.external.external_energy_params_str
  }
  refinement.geometry_restraints.edits
    .short_caption = Edit geometry restraints
    .style = menu_item box auto_align
  {
    include scope mmtbx.monomer_library.pdb_interpretation.geometry_restraints_edits_str
  }
  gui
    .help = Phenix GUI parameters, not used in command-line program
  {
    include scope libtbx.phil.interface.tracking_params
    output_dir = None
      .type = path
      .short_caption = Output directory
      .style = output_dir
  }
  extra_restraints_file = None
    .type = path
    .short_caption = Custom geometry restraints
    .help = File containing custom geometry restraints, using the same format \
      as phenix.refine.  On the command line this can be specified directly \
      as a command-line argument, but this parameter is used by the Phenix GUI.
    .style = input_file file_type:phil
  """, process_includes=True).fetch(source=customization_params)
  
  class er_pickle(object):
    def __init__(self,
                 pickle_object,
                 pickle_filename):
      cPickle.dump(pickle_object, gzip.open(pickle_filename, 'wb'))
  
  class ensemble_refinement_data(object):
    def __init__(self, f_calc_running                      = None,
                       f_calc_data_total                   = None,
                       f_calc_data_current                 = None,
                       f_mask_running                      = None,
                       f_mask_current                      = None,
                       f_mask_total                        = None,
                       total_SF_cntr                       = 0,
                       total_SF_cntr_mask                  = 0,
                       fix_scale_factor                    = None,
                       non_solvent_temp                    = None,
                       solvent_temp                        = None,
                       system_temp                         = None,
                       xray_structures                     = [],
                       pdb_hierarchys                      = [],
                       xray_structures_diff_map            = [],
                       seed                                = None,
                       velocities                          = None,
                       ke_protein_running                  = None,
                       ke_pdb                              = [],
                       geo_grad_rms                        = None,
                       xray_grad_rms                       = None,
                       solvent_sel                         = None,
                       all_sel                             = None,
                       er_harmonic_restraints_info         = None,
                       er_harmonic_restraints_weight       = 0.001,
                       er_harmonic_restraints_slack        = 1.0,
                       macro_cycle                         = None,
                       ):
      adopt_init_args(self, locals())
  
  class er_tls_manager(object):
    def __init__(self, tls_selection_strings_no_sol       = None,
                       tls_selection_strings_no_sol_no_hd = None,
                       tls_selections_with_sol            = None,
                       tls_selections_no_sol              = None,
                       tls_selections_no_sol_no_hd        = None,
                       tls_operators                      = None):
      adopt_init_args(self, locals())
  
  class run_ensemble_refinement(object):
    def __init__(self, fmodel,
                       model,
                       log,
                       raw_data,
                       raw_flags,
                       params,
                       ptls,
                       run_number=None):
      adopt_init_args(self, locals())
  #    self.params = params.extract().ensemble_refinement
  
      if self.params.target_name in ['ml', 'mlhl'] :
        self.fix_scale = False
      else:
        self.fix_scale = True
      if not self.params.wxray_coupled_tbath:
        self.params.wxray_coupled_tbath_offset = 0.0
      self.wxray = self.params.wxray
      self.params.ensemble_ordered_solvent.temperature = self.params.cartesian_dynamics.temperature
      self.ensemble_utils = ensemble_utils.manager(ensemble_obj = self)
      self.xray_gradient = None
      self.fc_running_ave = None
      self.macro_cycle = 1
      self.sf_model_ave = None
      self.fmodel_total_block_list = []
      self.reset_velocities = True
      self.cmremove = True
      self.cdp = self.params.cartesian_dynamics
      self.bsp = mmtbx.bulk_solvent.bulk_solvent_and_scaling.master_params.extract()
      if (self.params.target_name == 'mlhl'):
        self.bsp.target = 'ml'
      else :
        self.bsp.target = self.params.target_name
      if self.params.tx == None:
        print >> log, "\nAutomatically set Tx (parameter not defined)"
        print >> log, "Tx          :  2(1/dmin)**2"
        self.params.tx = round(2.0 * ((1.0/self.fmodel.f_obs().d_min())**2),1)
        print >> log, 'Dmin        : ', self.fmodel.f_obs().d_min()
        print >> log, 'Set Tx      : ', self.params.tx
      self.n_mc_per_tx = self.params.tx / (self.cdp.time_step * self.cdp.number_of_steps)
  
      # Set simulation length
      make_header("Simulation length:", out = self.log)
      print >> log, "Number of time steps per macro cycle    : ", self.cdp.number_of_steps
      print >> log, "Tx                                      : ", self.params.tx
      print >> log, "Number macro cycles per Tx period       : ", self.n_mc_per_tx
      self.equilibrium_macro_cycles = int(self.n_mc_per_tx * self.params.equilibrium_n_tx)
      self.acquisition_block_macro_cycles = int(self.n_mc_per_tx * self.params.acquisition_block_n_tx)
      self.total_macro_cycles = int(self.equilibrium_macro_cycles \
                              + (self.acquisition_block_macro_cycles * self.params.number_of_acquisition_periods))
      #
      print >> log, "\nEquilibration"
      print >> log, "    Number Tx periods    : ", self.params.equilibrium_n_tx
      print >> log, "    Number macro cycles  : ", self.equilibrium_macro_cycles
      print >> log, "    Time (ps)            : ", self.equilibrium_macro_cycles \
                                                    * self.cdp.number_of_steps * self.cdp.time_step
      #
      print >> log, "\nAcquisition block"
      print >> log, "    Number Tx periods    : ",  self.params.acquisition_block_n_tx
      print >> log, "    Number macro cycles  : ",  self.acquisition_block_macro_cycles
      print >> log, "    Time (ps)            : ",  self.acquisition_block_macro_cycles \
                                                    * self.cdp.number_of_steps\
                                                    * self.cdp.time_step
      #
      print >> log, "\nSimulation total"
      print >> log, "    Number Tx periods    : ", self.params.equilibrium_n_tx\
                                                  + (self.params.number_of_acquisition_periods\
                                                     * self.params.acquisition_block_n_tx)
      print >> log, "    Number macro cycles  : ", self.total_macro_cycles
      self.total_time = self.total_macro_cycles\
                          * self.cdp.number_of_steps\
                          * self.cdp.time_step
      print >> log, "    Time (ps)            : ", self.total_time
      print >> log, "    Total = Equilibration + nAcquisition"
      # Store block
      self.block_store_cycle_cntr = 0
      self.block_store_cycle = \
          range(self.acquisition_block_macro_cycles + self.equilibrium_macro_cycles,
                self.acquisition_block_macro_cycles + self.total_macro_cycles,
                self.acquisition_block_macro_cycles
                )
      # Store pdb
      self.pdb_store_cycle = max(int(self.acquisition_block_macro_cycles \
                           / self.params.pdb_stored_per_block), 1)
  
      #Setup ensemble_refinement_data_object
      self.er_data = ensemble_refinement_data()
      #Setup fmodels for running average   = refinement target
      #                  total average     = final model
      #                  current model     = model at time point n
      self.fmodel_running = self.fmodel
      self.fmodel_total = None
      self.fmodel_current = None
      self.tls_manager = None
      self.er_data.seed = self.params.random_seed
      self.run_time_stats_dict = {}
  
      #Dummy miller array
      self.copy_ma = self.fmodel_running.f_masks()[0].array(data = self.fmodel_running.f_masks()[0].data()*0).deep_copy()
      #
      self.fmodel_running.xray_structure = self.model.get_xray_structure()
      assert self.fmodel_running.xray_structure is self.model.get_xray_structure()
      self.pdb_hierarchy = self.model.get_hierarchy()
  
      #Atom selections
      self.atom_selections()
  
      self.model.geometry_statistics()
  
      self.setup_bulk_solvent_and_scale()
  
      self.fmodel_running.info(
        free_reflections_per_bin = 100,
        max_number_of_bins       = 999).show_rfactors_targets_in_bins(out = self.log)
  
      if self.params.target_name in ['ml', 'mlhl'] :
        #Must be called before reseting ADPs
        if self.params.scale_wrt_n_calc_start:
          make_header("Calculate Ncalc and restrain to scale kn", out = self.log)
          self.fmodel_running.n_obs_n_calc(update_nobs_ncalc = True)
          n_obs  = self.fmodel_running.n_obs
          n_calc = self.fmodel_running.n_calc
          self.scale_n1_reference = self.scale_helper(target    = n_calc,
                                                      reference = n_obs
                                                      )
          self.scale_n1_target    = self.scale_n1_reference
          self.scale_n1_current   = self.scale_n1_reference
          self.n_calc_reference = self.fmodel_running.n_calc.deep_copy()
          self.n_mc_per_ncalc_update = max(1, int(self.n_mc_per_tx / 10) )
          print >> self.log, "Number macro cycles per tx     : {0:5.0f}".format(self.n_mc_per_tx)
          print >> self.log, "Number macro cycles per update : {0:5.0f}".format(self.n_mc_per_ncalc_update)
          #
          self.fixed_k1_from_start = self.fmodel_running.scale_k1()
          self.target_k1 = self.fmodel_running.scale_k1()
          self.update_normalisation_factors()
        else:
          make_header("Calculate and fix scale of Ncalc", out = self.log)
          self.fmodel_running.n_obs_n_calc(update_nobs_ncalc = True)
          print >> self.log, "Fix Ncalc scale          : True"
          print >> self.log, "Sum current Ncalc        : {0:5.3f}".format(sum(self.fmodel_running.n_calc))
  
      #Set ADP model
      self.tls_manager = er_tls_manager()
      self.setup_tls_selections(tls_group_selection_strings = self.params.tls_group_selections)
      self.fit_tls(input_model = self.model)
      self.assign_solvent_tls_groups()
  
      #Set occupancies to 1.0
      if self.params.set_occupancies:
        make_header("Set occupancies to 1.0", out = self.log)
        self.model.get_xray_structure().set_occupancies(
          value      = 1.0)
        self.model.show_occupancy_statistics(out = self.log)
      #Initiates running average SFs
      self.er_data.f_calc_running = self.fmodel_running.f_calc().data().deep_copy()
      #self.fc_running_ave = self.fmodel_running.f_calc()
      self.fc_running_ave = self.fmodel_running.f_calc().deep_copy()
  
      #Initial sigmaa array, required for ML target function
      #Set eobs and ecalc normalization factors in Fmodel, required for ML
      if self.params.target_name in ['ml', 'mlhl'] :
        self.sigmaa_array = self.fmodel_running.sigmaa().sigmaa().data()
        self.best_r_free = self.fmodel_running.r_free()
        self.fmodel_running.set_sigmaa = self.sigmaa_array
  
      #Harmonic restraints
      if self.params.harmonic_restraints.selections is not None:
        self.add_harmonic_restraints()
  
  ############################## START Simulation ################################
      make_header("Start simulation", out = self.log)
      while self.macro_cycle <= self.total_macro_cycles:
        self.er_data.macro_cycle = self.macro_cycle
        self.time = self.cdp.time_step * self.cdp.number_of_steps * self.macro_cycle
        #XXX Debug
        if False and self.macro_cycle % 10==0:
          print >> self.log, "Sys temp  : ", self.er_data.system_temp
          print >> self.log, "Xray grad : ", self.er_data.xray_grad_rms
          print >> self.log, "Geo grad  : ", self.er_data.geo_grad_rms
          print >> self.log, "Wx        : ", self.wxray
  
        if self.fmodel_running.target_name in ['ml', 'mlhl'] :
          if self.macro_cycle < self.equilibrium_macro_cycles:
            if self.params.scale_wrt_n_calc_start and self.macro_cycle%self.n_mc_per_ncalc_update == 0:
              self.update_normalisation_factors()
            elif self.macro_cycle%int(self.n_mc_per_tx)==0:
              self.update_normalisation_factors()
  
        # Ordered Solvent Update
        if self.params.ordered_solvent_update \
            and (self.macro_cycle == 1\
            or self.macro_cycle%self.params.ordered_solvent_update_cycle == 0):
          self.ordered_solvent_update()
  
        xrs_previous = self.model.get_xray_structure().deep_copy_scatterers()
        assert self.fmodel_running.xray_structure is self.model.get_xray_structure()
  
        if self.cdp.verbose >= 1:
          if self.macro_cycle == 1 or self.macro_cycle%100 == 0:
            cdp_verbose = 1
          else:
            cdp_verbose = -1
        else:
          cdp_verbose = -1
  
        cd_manager = ensemble_cd.cartesian_dynamics(
          structure                   = self.model.get_xray_structure(),
          restraints_manager          = self.model.restraints_manager,
          temperature                 = self.cdp.temperature - self.params.wxray_coupled_tbath_offset,
          protein_thermostat          = self.cdp.protein_thermostat,
          n_steps                     = self.cdp.number_of_steps,
          n_print                     = self.cdp.n_print,
          time_step                   = self.cdp.time_step,
          initial_velocities_zero_fraction = self.cdp.initial_velocities_zero_fraction,
          fmodel                      = self.fmodel_running,
          xray_target_weight          = self.wxray,
          chem_target_weight          = 1.0,
          xray_structure_last_updated = None,
          shift_update                = 0.0,
          xray_gradient               = self.xray_gradient,
          reset_velocities            = self.reset_velocities,
          stop_cm_motion              = self.cmremove,
          update_f_calc               = False,
          er_data                     = self.er_data,
          verbose                     = cdp_verbose,
          log                         = self.log)
  
        self.reset_velocities = False
        self.cmremove = False
  
        #Calc rolling average KE energy
        self.kinetic_energy_running_average()
        #Show KE stats
        if self.params.verbose > 0 and self.macro_cycle % 500 == 0:
          self.ensemble_utils.kinetic_energy_stats()
  
        #Update Fmodel
        self.fmodel_running.update_xray_structure(
          xray_structure      = self.model.get_xray_structure(),
          update_f_calc       = True,
          update_f_mask       = True,
          force_update_f_mask = True)
  
        #Save current Fmask
        self.er_data.f_mask_current = self.fmodel_running.f_masks()[0].data().deep_copy()
  
        #Save current Fcalc
        self.er_data.f_calc_data_current = self.fmodel_running.f_calc().data().deep_copy()
  
        #Total Fmask calculation
        if self.er_data.f_mask_total is None:
          self.er_data.f_mask_total = self.fmodel_running.f_masks()[0].data().deep_copy()
          self.er_data.total_SF_cntr_mask = 1
        else:
          self.er_data.f_mask_total += self.fmodel_running.f_masks()[0].data().deep_copy()
          self.er_data.total_SF_cntr_mask += 1
  
        #Total Fcalc calculation
        if self.er_data.f_calc_data_total is None:
          self.er_data.f_calc_data_total = self.fmodel_running.f_calc().data().deep_copy()
          self.er_data.total_SF_cntr = 1
        else:
          self.er_data.f_calc_data_total += self.fmodel_running.f_calc().data().deep_copy()
          self.er_data.total_SF_cntr += 1
  
        #Running average Fcalc calculation
        if self.params.tx > 0:
          self.a_prime = math.exp(-(self.cdp.time_step * self.cdp.number_of_steps)/self.params.tx)
        else:
          self.a_prime = 0
  
        self.er_data.f_calc_running \
          = (self.a_prime * self.er_data.f_calc_running) + ((1-self.a_prime) * self.fmodel_running.f_calc().data().deep_copy())
        self.fc_running_ave = self.fc_running_ave.array(data = self.er_data.f_calc_running)
  
        #Update running average Fmask
        f_mask = self.fmodel_running.f_masks()[0]
        self.copy_ma, f_mask = self.copy_ma.common_sets(f_mask)
        if self.macro_cycle == 1:
          self.er_data.f_mask_running = f_mask.data().deep_copy()
        else:
          self.er_data.f_mask_running \
            = (self.a_prime * self.er_data.f_mask_running) + ((1-self.a_prime) * f_mask.data())
        self.running_f_mask_update = self.copy_ma.array(data = self.er_data.f_mask_running).deep_copy()
  
        #Update runnning average Fcalc and Fmask
        self.fmodel_running.update(f_calc = self.fc_running_ave,
                                   f_mask = self.running_f_mask_update)
  
        #Update total average Fcalc
        total_f_mask_update \
            = self.copy_ma.array(data = self.er_data.f_mask_total / self.er_data.total_SF_cntr_mask).deep_copy()
  
  
        if self.fmodel_total == None:
          self.fmodel_total = self.fmodel_running.deep_copy()
          self.fmodel_total.update(
            f_calc = self.copy_ma.array(data = self.er_data.f_calc_data_total / self.er_data.total_SF_cntr ),
            f_mask = total_f_mask_update)
  
          if(self.er_data.fix_scale_factor is not None):
            self.fmodel_total.set_scale_switch = self.er_data.fix_scale_factor
        else:
          self.fmodel_total.update(
            f_calc = self.copy_ma.array(data = self.er_data.f_calc_data_total / self.er_data.total_SF_cntr),
            f_mask = total_f_mask_update)
  
        #Update current time-step Fcalc
        current_f_mask_update = self.copy_ma.array(data = self.er_data.f_mask_current)
  
        if self.fmodel_current == None:
          self.fmodel_current = self.fmodel_running.deep_copy()
          self.fmodel_current.update(
            f_calc = self.copy_ma.array(data = self.er_data.f_calc_data_current),
            f_mask = current_f_mask_update)
          if(self.er_data.fix_scale_factor is not None):
            self.fmodel_current.set_scale_switch = self.er_data.fix_scale_factor
        else:
          self.fmodel_current.update(
            f_calc = self.copy_ma.array(data = self.er_data.f_calc_data_current),
            f_mask = current_f_mask_update)
  
        #ML params update
        if self.params.target_name in ['ml', 'mlhl'] :
          if self.macro_cycle < self.equilibrium_macro_cycles:
            if self.fmodel_running.r_free() < (self.best_r_free - self.params.update_sigmaa_rfree):
              self.update_sigmaa()
        # Wxray coupled to temperature bath
        if self.params.wxray_coupled_tbath:
          if self.macro_cycle < 5:
            self.wxray        = 2.5
          elif self.macro_cycle < self.equilibrium_macro_cycles:
            if self.params.tx == 0:
              a_prime_wx = 0
            else:
              wx_tx = min(self.time, self.params.tx)
              a_prime_wx = math.exp(-(self.cdp.time_step * self.cdp.number_of_steps)/wx_tx)
            wxray_t = self.wxray * max(0.01, self.cdp.temperature / self.er_data.non_solvent_temp)
            self.wxray = (a_prime_wx * self.wxray) + ((1-a_prime_wx) * wxray_t)
  
        #Store current structure, current KE
        if self.macro_cycle % self.pdb_store_cycle == 0 \
             and self.macro_cycle >= self.equilibrium_macro_cycles:
          self.er_data.xray_structures.append(self.model.get_xray_structure().deep_copy_scatterers())
          self.er_data.pdb_hierarchys.append(self.model.get_hierarchy().deep_copy())
          if self.er_data.ke_protein_running is None:
            self.er_data.ke_pdb.append(flex.double(self.model.get_xray_structure().sites_cart().size(), 0.0) )
          else:
            ke_expanded = flex.double(self.model.get_sites_cart().size(), 0.0)
            ke_expanded.set_selected(~self.model.solvent_selection(),
                                     self.er_data.ke_protein_running)
            self.er_data.ke_pdb.append(ke_expanded)
  
        #Current structural deviation vs starting structure and previous macro-cycle structure
        if xrs_previous.distances(other = self.model.get_xray_structure()).min_max_mean().mean > 1.0:
          print >> self.log, "\n\nWARNING:"
          print >> self.log, "Macro cycle too long, max atomic deviation w.r.t. previous cycle"
          print >> self.log, "greater than 1.0A"
          print >> self.log, "Reduce params.cartesian_dynamics.number_of_steps"
          print >> self.log, "Max deviation : {0:1.3f}"\
            .format(xrs_previous.distances(other = self.model.get_xray_structure()).min_max_mean().mean)
  
        if self.fmodel_running.r_work() > 0.75:
          raise Sorry("Simulation aborted, running Rfree > 75%")
  
        #Print run time stats
        if self.macro_cycle == 1 or self.macro_cycle%50 == 0:
          print >> self.log, "\n________________________________________________________________________________"
          print >> self.log, "    MC        Time     |  Current  |  Rolling  |   Total   | Temp |  Grad Wxray "
          print >> self.log, "          (ps)     (%) |   Rw   Rf |   Rw   Rf |   Rw   Rf |  (K) |   X/G       "
        print >> self.log, \
            "~{0:5d} {1:7.2f} {2:7.2f} | {3:4.1f} {4:4.1f} | {5:4.1f} {6:4.1f} | {7:4.1f} {8:4.1f} | {9:4.0f} | {10:5.2f} {11:5.2f}"\
            .format(self.macro_cycle,
                    self.time,
                    100 * self.time / self.total_time,
                    100*self.fmodel_current.r_work(),
                    100*self.fmodel_current.r_free(),
                    100*self.fmodel_running.r_work(),
                    100*self.fmodel_running.r_free(),
                    100*self.fmodel_total.r_work(),
                    100*self.fmodel_total.r_free(),
                    self.er_data.non_solvent_temp,
                    self.er_data.xray_grad_rms / self.er_data.geo_grad_rms,
                    self.wxray)
  
        if self.params.verbose > 0:
          if self.macro_cycle == 1\
              or self.macro_cycle%100 == 0\
              or self.macro_cycle == self.total_macro_cycles:
            self.print_fmodels_scale_and_solvent_stats()
  
        if self.params.number_of_acquisition_periods > 1:
          if self.macro_cycle in self.block_store_cycle:
            self.save_multiple_fmodel()
  
        #End of equilibration period, reset total structure factors, atomic cords, kinetic energies
        if self.macro_cycle == self.equilibrium_macro_cycles:
          self.reset_totals()
        #
        assert self.model.get_xray_structure() is cd_manager.structure
        assert self.fmodel_running.xray_structure is cd_manager.structure
        if self.fix_scale == True:
          assert self.fmodel_running.scale_k1() == self.er_data.fix_scale_factor
        self.macro_cycle +=1
  
  ############################## END Simulation ##################################
  
      self.macro_cycle = self.total_macro_cycles
      #Find optimum section of acquisition period
      if self.params.number_of_acquisition_periods > 1:
        self.optimise_multiple_fmodel()
      else:
        self.fmodel_total.set_scale_switch = 0
        self.fmodel_total.update_all_scales(
          log    = self.log,
          remove_outliers=False,
          params = self.bsp)
      #Minimize number of ensemble models
      if self.params.ensemble_reduction:
        self.ensemble_utils.ensemble_reduction()
  
      #Optimise fmodel_total k, b_aniso, k_sol, b_sol
      self.fmodel_total.set_scale_switch = 0
      self.print_fmodels_scale_and_solvent_stats()
      self.fmodel_total.update_all_scales(
        log    = self.log,
        remove_outliers=False,
        params = self.bsp)
      self.print_fmodels_scale_and_solvent_stats()
      print >> self.log, "FINAL Rwork = %6.4f Rfree = %6.4f Rf/Rw = %6.4f"\
          %(self.fmodel_total.r_work(),
            self.fmodel_total.r_free(),
            self.fmodel_total.r_free() / self.fmodel_total.r_work()
            )
      print >> self.log, "Final Twork = %6.4f Tfree = %6.4f Tf/Tw = %6.4f"\
          %(self.fmodel_total.target_w(),
            self.fmodel_total.target_t(),
            self.fmodel_total.target_t() / self.fmodel_total.target_w()
            )
      info = self.fmodel_total.info(free_reflections_per_bin = 100,
                                    max_number_of_bins       = 999
                                    )
      info.show_remark_3(out = self.log)
      info.show_rfactors_targets_in_bins(out = self.log)
  
      self.write_output_files(run_number=run_number)
  
  ############################## END ER ##########################################
  
    def write_output_files(self, run_number=None):
      #PDB output
      prefix = self.params.output_file_prefix
      if (run_number is not None):
        prefix += "_%g" % run_number
      pdb_out = prefix + ".pdb"
      cif_out = prefix + ".cif"
      if (self.params.gzip_final_model):
        pdb_out += ".gz"
        self.write_ensemble_pdb(out = gzip.open(pdb_out, 'wb'))
        # TODO
        if False :#(self.params.write_cif_file):
          self.write_ensemble_mmcif(out=gzip.open(cif_out, 'wb'))
      else :
        self.write_ensemble_pdb(out = open(pdb_out, 'wb'))
        # TODO
        if False :#(self.params.write_cif_file):
          self.write_ensemble_mmcif(out=open(cif_out, 'wb'))
      self.pdb_file = pdb_out
      # Map output
      assert (self.fmodel_total is not None)
      self.mtz_file = write_mtz_file(
        fmodel_total=self.fmodel_total,
        raw_data=self.raw_data,
        raw_flags=self.raw_flags,
        prefix=prefix,
        params=self.params)
  
    def show_overall(self, message = "", fmodel_running = True):
      if fmodel_running:
        message = "Running: " + message
        self.fmodel_running.info().show_rfactors_targets_scales_overall(header = message, out = self.log)
      else:
        message = "Total: " + message
        self.fmodel_total.info().show_rfactors_targets_scales_overall(header = message, out = self.log)
  
    def add_harmonic_restraints(self):
      make_header("Add specific harmonic restraints", out = self.log)
      # ensures all solvent atoms are at the end prior to applying harmonic restraints
      self.ordered_solvent_update()
      hr_selections = mmtbx.utils.get_atom_selections(
          model = self.model,
          selection_strings = self.params.harmonic_restraints.selections)
      pdb_atoms = self.pdb_hierarchy.atoms()
      print >> self.log, "\nAdd atomic harmonic restraints:"
      restraint_info = []
      for i_seq in hr_selections[0]:
        atom_info = pdb_atoms[i_seq].fetch_labels()
        print >> self.log, '    {0} {1} {2} {3} {4}     '.format(
                                     atom_info.name,
                                     atom_info.i_seq+1,
                                     atom_info.resseq,
                                     atom_info.resname,
                                     atom_info.chain_id,
                                     )
        restraint_info.append((i_seq, pdb_atoms[i_seq].xyz))
      self.er_data.er_harmonic_restraints_info = restraint_info
      self.er_data.er_harmonic_restraints_weight = self.params.harmonic_restraints.weight
      self.er_data.er_harmonic_restraints_slack  = self.params.harmonic_restraints.slack
      print >> self.log, "\n|"+"-"*77+"|\n"
  
    def setup_bulk_solvent_and_scale(self):
      make_header("Setup bulk solvent and scale", out = self.log)
      self.show_overall(message = "pre solvent and scale")
      #
      self.fmodel_running.update_all_scales(
        params = self.bsp,
        remove_outliers=False,
        log    = self.log)
      self.fmodel_running.optimize_mask(params = self.bsp)
      #Fixes scale factor for rolling average #ESSENTIAL for LSQ
      if self.fix_scale == True:
        self.er_data.fix_scale_factor = self.fmodel_running.scale_k1()
        self.fmodel_running.set_scale_switch = self.er_data.fix_scale_factor
      self.show_overall(message = "post solvent and scale")
  
    def scale_helper(self, reference, target):
      return flex.sum(reference * target) / flex.sum(flex.pow2(target))
  
    def update_normalisation_factors(self):
      if self.params.scale_wrt_n_calc_start:
        # Adaptive scaling
        # Ncalc_start / Ncalc_current
        make_header("Update Ncalc and restrain to Ncalc ref", out = self.log)
        # Get N_calc current, compare with reference
        n_obs, n_calc =\
          self.fmodel_running.n_obs_n_calc(update_nobs_ncalc = False)
        ref_div_current = self.n_calc_reference / n_calc
  
        n_calc_coeff    = 1.0-math.exp(-self.n_mc_per_ncalc_update/self.n_mc_per_tx)
        n_calc_scaled   = ref_div_current * n_calc_coeff
        n_calc_update   = (self.fmodel_running.n_calc * (1.0-n_calc_coeff) ) + (self.fmodel_running.n_calc * ref_div_current * n_calc_coeff)
  
        # Update with scaled array
        self.fmodel_running.n_calc = n_calc_update
  
      else:
        # Normalise to reference Sum(Ncalc)
        make_header("Update and renormalise Ncalc array", out = self.log)
        eobs_norm_factor, ecalc_norm_factor =\
          self.fmodel_running.n_obs_n_calc(update_nobs_ncalc = False)
        self.scale_n1_current = self.scale_helper(target    = ecalc_norm_factor,
                                                  reference = eobs_norm_factor
                                                  )
        print >> self.log, "Kn current               : {0:5.3f}".format(self.scale_n1_current)
        ecalc_k = sum(self.fmodel_running.n_calc) / sum(ecalc_norm_factor)
        ecalc_k_alt = flex.sum(self.fmodel_running.n_calc * ecalc_norm_factor) / flex.sum(flex.pow2(ecalc_norm_factor) )
        print >> self.log, "Sum current Ncalc        : {0:5.3f}".format(sum(self.fmodel_running.n_calc) )
        print >> self.log, "Sum updated Ncalc        : {0:5.3f}".format(sum(ecalc_norm_factor) )
        print >> self.log, "Rescaling factor         : {0:5.3f}".format(ecalc_k)
        print >> self.log, "Rescaling factor alt     : {0:5.3f}".format(ecalc_k_alt)
        ecalc_norm_factor = ecalc_k * ecalc_norm_factor
        self.fmodel_running.n_calc = ecalc_norm_factor
      print >> self.log, "|"+"-"*77+"|\n"
  
    def update_sigmaa(self):
      make_header("Update sigmaa", out = self.log)
      if self.params.verbose > 0:
        print >> self.log, "Previous best Rfree      : ", self.best_r_free
        print >> self.log, "Current       Rfree      : ", self.fmodel_running.r_free()
        self.print_ml_stats()
        print >> self.log, "  Update sigmaa"
      self.sigmaa_array = self.fmodel_running.sigmaa().sigmaa().data()
      self.fmodel_running.set_sigmaa = self.sigmaa_array
      if self.params.verbose > 0:
        self.print_ml_stats()
      self.best_r_free = self.fmodel_running.r_free()
      print >> self.log, "|"+"-"*77+"|\n"
  
    def setup_tls_selections(self, tls_group_selection_strings):
      make_header("Generating TLS selections from input parameters (not including solvent)", out = self.log)
      model_no_solvent = self.model.deep_copy()
      model_no_solvent = model_no_solvent.remove_solvent()
  
      if len(tls_group_selection_strings) < 1:
        print >> self.log, '\nNo TLS groups supplied - automatic setup'
        # Get chain information
        chains_info = []
        for chain in model_no_solvent.get_hierarchy().chains():
          count_h = 0
          for atom in chain.atoms():
            if atom.element_is_hydrogen(): count_h+=1
          chain_id_non_h = ("'" + chain.id + "'", chain.atoms_size() - count_h)
          chains_info.append(chain_id_non_h)
        # Check all chains > 63 heavy atoms for TLS fitting
        chains_size = flex.int(zip(*chains_info)[1])
        chains_size_ok = flex.bool(chains_size > 63)
        if sum(chains_size) < 63:
          print >> self.log, '\nStructure contains less than 63 atoms (non H/D, non solvent)'
          print >> self.log, '\nUnable to perform TLS fitting, will use isotropic B-factor model'
        elif chains_size_ok.count(False) == 0:
          print >> self.log, '\nTLS selections:'
          print >> self.log, 'Chain, number atoms (non H/D)'
          for chain in chains_info:
            tls_group_selection_strings.append('chain ' + chain[0])
            print >> self.log, chain[0], chain[1]
        else:
          print >> self.log, '\nFollowing chains contain less than 63 atoms (non H/D):'
          tls_group_selection_strings.append('chain ')
          for chain in chains_info:
            tls_group_selection_strings[0] += (chain[0] + ' or chain ')
            if chain[1] < 63:
              print >> self.log, chain[0], chain[1]
          print >> self.log, 'Combining all chains to single TLS group'
          print >> self.log, 'WARNING: this may not be the optimum tls groupings to use'
          print >> self.log, 'TLS selections:'
          tls_group_selection_strings[0] = tls_group_selection_strings[0][0:-10]
          print >> self.log, tls_group_selection_strings[0]
      #
      tls_no_sol_selections =  mmtbx.utils.get_atom_selections(
          model = model_no_solvent,
          selection_strings = tls_group_selection_strings)
      #
      tls_no_hd_selection_strings = []
      for selection_string in tls_group_selection_strings:
        no_hd_string = '(' + selection_string + ') and not (element H or element D)'
        tls_no_hd_selection_strings.append(no_hd_string)
  
      tls_no_sol_no_hd_selections = mmtbx.utils.get_atom_selections(
          model = model_no_solvent,
          selection_strings = tls_no_hd_selection_strings)
  
      #
      assert self.tls_manager is not None
      self.tls_manager.tls_selection_strings_no_sol       = tls_group_selection_strings
      self.tls_manager.tls_selection_strings_no_sol_no_hd = tls_no_hd_selection_strings
      self.tls_manager.tls_selections_no_sol              = tls_no_sol_selections
      self.tls_manager.tls_selections_no_sol_no_hd        = tls_no_sol_no_hd_selections
      self.tls_manager.tls_operators = mmtbx.tls.tools.generate_tlsos(
          selections     = self.tls_manager.tls_selections_no_sol,
          xray_structure = model_no_solvent.get_xray_structure(),
          value          = 0.0)
  
      self.model.tls_groups = mmtbx.tls.tools.tls_groups(
          selection_strings = self.tls_manager.tls_selection_strings_no_sol,
          tlsos             = self.tls_manager.tls_operators)
  
    def fit_tls(self, input_model, verbose = False):
      make_header("Fit TLS from reference model", out = self.log)
      model_copy = input_model.deep_copy()
      model_copy = model_copy.remove_solvent()
      print >> self.log, 'Reference model :'
      model_copy.show_adp_statistics(padded = True, out = self.log)
      start_xrs = model_copy.get_xray_structure().deep_copy_scatterers()
      start_xrs.convert_to_isotropic()
      start_biso = start_xrs.scatterers().extract_u_iso()/adptbx.b_as_u(1)
      model_copy.get_xray_structure().convert_to_anisotropic()
      tls_selection_no_sol_hd            = self.tls_manager.tls_selections_no_sol_no_hd
      tls_selection_no_sol_hd_exclusions = self.tls_manager.tls_selections_no_sol_no_hd
      pre_fitted_mean = 999999.99
      #
      use_isotropic = False
      for group in self.tls_manager.tls_selections_no_sol_no_hd:
        if group.size() < 63:
          self.params.isotropic_b_factor_model = True
        elif self.ptls * group.size() < 63:
          self.ptls = 64.0 / group.size()
          print >> self.log, '\nAutomatically increasing pTLS to : {0:5.3f}'.format(self.ptls)
      if self.params.isotropic_b_factor_model:
        print >> self.log, '\nModel contains less than 63 non-solvent, non-H/D atoms'
        print >> self.log, 'Insufficient to fit TLS model, using isotropic model'
        iso_b  = self.fmodel_running.wilson_b() * self.params.pwilson
        episq = 8.0*(math.pi**2)
        print >> self.log, 'Isotropic translation (B) : {0:5.3f}'.format(iso_b)
        print >> self.log, '  = Wilson b-factor * pwilson'
        iso_u = iso_b / episq
        print >> self.log, 'Isotropic translation (U) : {0:5.3f}'.format(iso_u)
        fit_tlsos = []
        for tls_group in self.tls_manager.tls_operators:
          tls_t_new = (iso_u,
                       iso_u,
                       iso_u,
                       0.0,
                       0.0,
                       0.0)
          tls_l_new = (0.0,0.0,0.0,0.0,0.0,0.0)
          tls_s_new = (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
          fit_tlsos.append(tls_tools.tlso(t      = tls_t_new,
                                        l      = tls_l_new,
                                        s      = tls_s_new,
                                        origin = tls_group.origin))
        self.tls_manager.tls_operators = fit_tlsos
        for tls_group in self.tls_manager.tls_operators:
          mmtbx.tls.tools.show_tls_one_group(tlso = tls_group,
                                             out  = self.log)
  
      else:
        for fit_cycle in xrange(self.params.max_ptls_cycles):
          fit_tlsos = mmtbx.tls.tools.generate_tlsos(
            selections     = tls_selection_no_sol_hd_exclusions,
            xray_structure = model_copy.get_xray_structure(),
            value          = 0.0)
          print >> self.log, '\nFitting cycle : ', fit_cycle+1
          for rt,rl,rs in [[1,0,1],[1,1,1],[0,1,1],
                           [1,0,0],[0,1,0],[0,0,1],[1,1,1],
                           [0,0,1]]*10:
            fit_tlsos = mmtbx.tls.tools.tls_from_uanisos(
              xray_structure               = model_copy.get_xray_structure(),
              selections                   = tls_selection_no_sol_hd_exclusions,
              tlsos_initial                = fit_tlsos,
              number_of_macro_cycles       = 10,
              max_iterations               = 100,
              refine_T                     = rt,
              refine_L                     = rl,
              refine_S                     = rs,
              enforce_positive_definite_TL = True,
              verbose                      = -1,
              out                          = self.log)
          fitted_tls_xrs = model_copy.get_xray_structure().deep_copy_scatterers()
          us_tls = mmtbx.tls.tools.u_cart_from_tls(
                 sites_cart = fitted_tls_xrs.sites_cart(),
                 selections = self.tls_manager.tls_selections_no_sol,
                 tlsos      = fit_tlsos)
          fitted_tls_xrs.set_u_cart(us_tls)
          fitted_tls_xrs.convert_to_isotropic()
          fitted_biso = fitted_tls_xrs.scatterers().extract_u_iso()/adptbx.b_as_u(1)
          mmtbx.tls.tools.show_tls(tlsos = fit_tlsos, out = self.log)
          #For testing
          if verbose:
            pdb_hierarchy = model_copy.pdb_hierarchy
            pdb_atoms = pdb_hierarchy().atoms()
            not_h_selection = pdb_hierarchy().atom_selection_cache().selection('not element H')
            ca_selection = pdb_hierarchy().atom_selection_cache().selection('name ca')
            print >> self.log, '\nCA atoms (Name/res number/res name/chain/atom number/ref biso/fit biso:'
            for i_seq, ca in enumerate(ca_selection):
              if ca:
                atom_info = pdb_atoms[i_seq].fetch_labels()
                print >> self.log, atom_info.name, atom_info.resseq, atom_info.resname, atom_info.chain_id, " | ", i_seq, start_biso[i_seq], fitted_biso[i_seq]
  
          delta_ref_fit = start_biso - fitted_biso
          hd_selection = model_copy.get_hd_selection()
          delta_ref_fit_no_h = delta_ref_fit.select(~hd_selection)
          delta_ref_fit_no_h_basic_stats = scitbx.math.basic_statistics(delta_ref_fit_no_h )
          start_biso_no_hd = start_biso.select(~hd_selection)
          fitted_biso_no_hd = fitted_biso.select(~hd_selection)
  
          if verbose:
            print >> self.log, 'pTLS                                    : ', self.ptls
  
          sorted_delta_ref_fit_no_h = sorted(delta_ref_fit_no_h)
          percentile_cutoff = sorted_delta_ref_fit_no_h[int(len(sorted_delta_ref_fit_no_h) * self.ptls)-1]
          if verbose:
            print >> self.log, 'Cutoff (<)                              : ', percentile_cutoff
  
          print >> self.log, 'Number of atoms (non HD)                : ', delta_ref_fit_no_h.size()
          delta_ref_fit_no_h_include = flex.bool(delta_ref_fit_no_h < percentile_cutoff)
          print >> self.log, 'Number of atoms (non HD) used in fit    : ', delta_ref_fit_no_h_include.count(True)
          print >> self.log, 'Percentage (non HD) used in fit         : {0:5.3f}'.format(delta_ref_fit_no_h_include.count(True) / delta_ref_fit_no_h.size())
  
          # Convergence test
          if fitted_biso_no_hd.min_max_mean().mean == pre_fitted_mean:
            break
          else:
            pre_fitted_mean = fitted_biso_no_hd.min_max_mean().mean
  
          # N.B. map on to full array including hydrogens for i_seqs
          include_array = flex.bool(delta_ref_fit < percentile_cutoff)
          #
          include_i_seq = []
          assert delta_ref_fit.size() == model_copy.get_number_of_atoms()
          assert include_array.size() == model_copy.get_number_of_atoms()
          for i_seq, include_flag in enumerate(include_array):
            if include_flag and not hd_selection[i_seq]:
              include_i_seq.append(i_seq)
          tls_selection_no_sol_hd_exclusions = []
          for group in xrange(len(tls_selection_no_sol_hd)):
            new_group = flex.size_t()
            for x in tls_selection_no_sol_hd[group]:
              if x in include_i_seq:
                new_group.append(x)
            if len(new_group) < 63:
              raise Sorry("Number atoms in TLS too small; increase size of group or reduce cut-off")
            if verbose:
              print >> self.log, 'TLS group ', group+1, ' number atoms ', len(new_group)
            tls_selection_no_sol_hd_exclusions.append(new_group)
  
      print >> self.log, '\nFinal non-solvent b-factor model'
      model_copy.get_xray_structure().convert_to_anisotropic()
  
      us_tls = mmtbx.tls.tools.u_cart_from_tls(
               sites_cart = model_copy.get_sites_cart(),
               selections = self.tls_manager.tls_selections_no_sol_no_hd,
               tlsos      = fit_tlsos)
      model_copy.get_xray_structure().set_u_cart(us_tls)
      model_copy.show_adp_statistics(padded = True, out = self.log)
      del model_copy
  
      #Update TLS params
      self.model.tls_groups.tlsos = fit_tlsos
      self.tls_manager.tls_operators = fit_tlsos
      self.assign_solvent_tls_groups()
  
    def tls_parameters_update(self):
      self.model.get_xray_structure().convert_to_anisotropic()
      #Apply TLS w.r.t. to atomic position
      selections = self.tls_manager.tls_selections_with_sol
      us_tls = mmtbx.tls.tools.u_cart_from_tls(
                 sites_cart = self.model.get_sites_cart(),
                 selections = selections,
                 tlsos      = self.tls_manager.tls_operators)
      for selection in selections:
        self.model.get_xray_structure().set_u_cart(us_tls, selection = selection)
      self.fmodel_running.update_xray_structure(
        xray_structure = self.model.get_xray_structure(),
        update_f_calc  = False,
        update_f_mask  = False)
  
    def assign_solvent_tls_groups(self):
      self.model.get_xray_structure().convert_to_anisotropic(selection =  self.model.solvent_selection())
      self.fmodel_running.update_xray_structure(
        xray_structure  = self.model.get_xray_structure(),
        update_f_calc   = False,
        update_f_mask   = False)
      #
      self.tls_manager.tls_selections_with_sol = []
      for grp in self.tls_manager.tls_selections_no_sol:
        self.tls_manager.tls_selections_with_sol.append(grp.deep_copy())
      #
      if len(self.tls_manager.tls_selections_with_sol) == 1:
        pdb_atoms     = self.pdb_hierarchy.atoms()
        hoh_selection = self.model.solvent_selection()
        for n, atom in enumerate(pdb_atoms):
          if hoh_selection[n]:
            self.tls_manager.tls_selections_with_sol[0].append(n)
      else:
        model             = self.model.deep_copy()
        solvent_selection = model.solvent_selection()
        solvent_xrs       = model.get_xray_structure().select(solvent_selection)
        model             = model.remove_solvent()
        closest_distances = model.get_xray_structure().closest_distances(
                                sites_frac      = solvent_xrs.sites_frac(),
                                use_selection   = ~model.get_xray_structure().hd_selection(),
                                distance_cutoff = 10.0)
        assert len(solvent_xrs.sites_cart()) == len(closest_distances.i_seqs)
        number_non_solvent_atoms = model.get_number_of_atoms()
        for n, i_seq in enumerate(closest_distances.i_seqs):
          for grp in self.tls_manager.tls_selections_with_sol:
            if i_seq in grp:
              grp.append(n+number_non_solvent_atoms)
              break
      #
      self.tls_parameters_update()
  
    def kinetic_energy_running_average(self):
      #Kinetic energy
      atomic_weights = self.model.get_xray_structure().atomic_weights()
      ke = 0.5 * atomic_weights * self.er_data.velocities.dot()
      #Select non-solvent atoms
      ke = ke.select(~self.model.solvent_selection() )
      if self.er_data.ke_protein_running == None:
        self.er_data.ke_protein_running = ke
      else:
        self.er_data.ke_protein_running\
          = (self.a_prime * self.er_data.ke_protein_running) + ( (1-self.a_prime) * ke)
  
    def ordered_solvent_update(self):
      ensemble_ordered_solvent_manager = ensemble_ordered_solvent.manager(
          model             = self.model,
          fmodel            = self.fmodel_running,
          verbose           = self.params.verbose,
          params            = self.params.ensemble_ordered_solvent,
          velocities        = self.er_data.velocities,
          log               = self.log)
      self.model = ensemble_ordered_solvent_manager.model
      self.er_data.velocities = ensemble_ordered_solvent_manager.velocities
      self.fmodel_running.update_xray_structure(
        xray_structure = self.model.get_xray_structure(),
        update_f_calc  = False,
        update_f_mask  = False)
      assert self.fmodel_running.xray_structure is self.model.get_xray_structure()
      self.xray_gradient = None
      #Update atom selections
      self.pdb_hierarchy = self.model.get_hierarchy()
      self.atom_selections()
      #Reset solvent tls groups
      if self.tls_manager is not None:
        self.assign_solvent_tls_groups()
  
    def reset_totals(self):
      make_header("Reseting structure ensemble and total Fmodel", out = self.log)
      self.er_data.xray_structures = []
      self.er_data.xray_structures_diff_map = []
      self.er_data.pdb_hierarchys = []
      self.er_data.ke_pdb = []
      self.er_data.f_calc_data_total = None
      self.er_data.total_SF_cntr = 0
      self.er_data.f_mask_total = None
      self.er_data.total_SF_cntr_mask = 0
  
    #Generates list of atom selections (needed to pass to CD)
    def atom_selections(self):
      self.er_data.all_sel     = flex.bool(self.model.get_number_of_atoms(), True)
      self.er_data.solvent_sel = self.model.solvent_selection()
  
    def save_multiple_fmodel(self):
      make_header("Saving fmodel block", out = self.log)
      #Stores fcalc, fmask, xray structure, pdb hierarchys
      print >> self.log, '{0:<23}: {1:>8} {2:>8} {3:>8} {4:>8}'.format('','MC','Block','Rwork','Rfree')
      print >> self.log, "{0:<23}: {1:8d} {2:8d} {3:8.3f} {4:8.3f}".format(
          'Fmodel block info',
          self.macro_cycle,
          self.block_store_cycle_cntr+1,
          100 * self.fmodel_total.r_work(),
          100 * self.fmodel_total.r_free() )
      fcalc_block  = self.er_data.f_calc_data_total / self.er_data.total_SF_cntr
      fmask_block  = self.er_data.f_mask_total / self.er_data.total_SF_cntr_mask
      xrs_block    = self.er_data.xray_structures
      xrs_dm_block = self.er_data.xray_structures_diff_map
      pdb_h_block  = self.er_data.pdb_hierarchys
      ke_pdb_block = self.er_data.ke_pdb
  
      block_info = (fcalc_block,
                    fmask_block,
                    xrs_block,
                    xrs_dm_block,
                    pdb_h_block,
                    ke_pdb_block)
  
      self.block_store_cycle_cntr+1
      if self.block_store_cycle_cntr+1 == 1:
        self.block_temp_file_list = []
      prefix = self.params.output_file_prefix
      if (self.run_number is not None):
        prefix += "_%g" % self.run_number
      filename = str(self.block_store_cycle_cntr+1)+'_block_'+prefix+'_TEMP.pZ'
      self.block_temp_file_list.append(filename)
      er_pickle(pickle_object = block_info, pickle_filename = filename)
      self.block_store_cycle_cntr += 1
      if self.macro_cycle != self.total_macro_cycles:
        self.reset_totals()
  
    def optimise_multiple_fmodel(self):
      make_header("Block selection by Rwork", out = self.log)
      best_r_work = None
  
      # Load all temp files
      self.fmodel_total_block_list = []
      for filename in self.block_temp_file_list:
        block_info = cPickle.load(gzip.open(filename,'rb'))
        self.fmodel_total_block_list.append(block_info)
        os.remove(filename)
  
      self.fmodel_total.set_scale_switch = 0
      print >> self.log, '  {0:>17} {1:>8} {2:>8}'\
        .format('Block range','Rwork','Rfree','k1')
      for x in xrange(len(self.fmodel_total_block_list)):
        x2 = x+1
        y = len(self.fmodel_total_block_list)
        while y > x:
          fcalc_tot = self.fmodel_total_block_list[x][0].deep_copy()
          fmask_tot = self.fmodel_total_block_list[x][1].deep_copy()
          cntr      = 1
          while x2 < y:
            fcalc_tot += self.fmodel_total_block_list[x2][0].deep_copy()
            fmask_tot += self.fmodel_total_block_list[x2][1].deep_copy()
            x2     += 1
            cntr   += 1
          self.fmodel_total.update(
            f_calc = self.copy_ma.array(data = (fcalc_tot / cntr)),
            f_mask = self.copy_ma.array(data = (fmask_tot / cntr)) )
          self.fmodel_total.update_all_scales(
            params = self.bsp,
            remove_outliers=False,
            log = self.log)
          print >> self.log, "  {0:8d} {1:8d} {2:8.3f} {3:8.3f}"\
            .format(x+1,
                    y,
                    self.fmodel_total.r_work(),
                    self.fmodel_total.r_free(),
                    self.fmodel_total.scale_k1()
                    )
          if best_r_work == None:
            best_r_work = self.fmodel_total.r_work()
            best_r_work_block = [x,y]
            best_r_work_fcalc = (fcalc_tot / cntr)
            best_r_work_fmask = (fmask_tot / cntr)
          elif self.fmodel_total.r_work() < (best_r_work - 0.01):
            best_r_work = self.fmodel_total.r_work()
            best_r_work_block = [x,y]
            best_r_work_fcalc = (fcalc_tot / cntr)
            best_r_work_fmask = (fmask_tot / cntr)
          x2 = x+1
          y -= 1
      self.fmodel_total.update(
        f_calc = self.copy_ma.array(data = best_r_work_fcalc),
        f_mask = self.copy_ma.array(data = best_r_work_fmask) )
      self.fmodel_total.update_all_scales(
            params = self.bsp,
            remove_outliers=False,
            log    = self.log)
  
      print >> self.log, "\nOptimium block :"
      print >> self.log, "  {0:8d} {1:8d} {2:8.3f} {3:8.3f} {4:8.3f} {5:8.3f}"\
        .format(best_r_work_block[0]+1,
                best_r_work_block[1],
                self.fmodel_total.r_work(),
                self.fmodel_total.r_free(),
                self.fmodel_total.scale_k1(),
                self.fmodel_total.fmodel_kbu().k_sols()[0],
                self.fmodel_total.fmodel_kbu().b_sols()[0])
      #Update self.er_data.xray_structures and self.er_data.pdb_hierarchys to correspond to optimum fmodel_total
      self.er_data.xray_structures = []
      self.er_data.xray_structures_diff_map =[]
      self.er_data.pdb_hierarchys  = []
      self.er_data.ke_pdb          = []
      for x in xrange(len(self.fmodel_total_block_list)):
        if x >= best_r_work_block[0] and x < best_r_work_block[1]:
          print  >> self.log, "Block | Number of models in block : ", x+1, " | ", len(self.fmodel_total_block_list[x][2])
          self.er_data.xray_structures.extend(self.fmodel_total_block_list[x][2])
          self.er_data.xray_structures_diff_map.extend(self.fmodel_total_block_list[x][3])
          self.er_data.pdb_hierarchys.extend(self.fmodel_total_block_list[x][4])
          self.er_data.ke_pdb.extend(self.fmodel_total_block_list[x][5])
      assert len(self.er_data.xray_structures) == len(self.er_data.pdb_hierarchys)
      assert len(self.er_data.xray_structures) == len(self.er_data.ke_pdb)
      print >> self.log, "Number of models for PBD          : ", len(self.er_data.xray_structures)
      print >> self.log, "|"+"-"*77+"|\n"
  
    def print_fmodels_scale_and_solvent_stats(self):
      make_header("Fmodel statistics | macrocycle: "+str(self.macro_cycle),
        out = self.log)
      print >> self.log, '{0:<23}: {1:>8} {2:>8} {3:>8} {4:>8}'.format('','MC',
        'k1','Bsol','ksol')
      if self.fmodel_current is not None:
        print >> self.log, "{0:<23}: {1:8d} {2:8.3f} {3:8.3f} {4:8.3f}"\
          .format('Fmodel current',
                  self.macro_cycle,
                  self.fmodel_current.scale_k1(),
                  self.fmodel_current.fmodel_kbu().b_sols()[0],
                  self.fmodel_current.fmodel_kbu().k_sols()[0],
                  )
      if self.fmodel_running is not None:
        print >> self.log, "{0:<23}: {1:8d} {2:8.3f} {3:8.3f} {4:8.3f}"\
          .format('Fmodel running',
                  self.macro_cycle,
                  self.fmodel_running.scale_k1(),
                  self.fmodel_running.fmodel_kbu().b_sols()[0],
                  self.fmodel_running.fmodel_kbu().k_sols()[0] )
      if self.fmodel_total is not None:
        print >> self.log, "{0:<23}: {1:8d} {2:8.3f} {3:8.3f} {4:8.3f}"\
          .format('Fmodel_Total',
                  self.macro_cycle,
                  self.fmodel_total.scale_k1(),
                  self.fmodel_total.fmodel_kbu().b_sols()[0],
                  self.fmodel_total.fmodel_kbu().k_sols()[0] )
      if self.fmodel_current is not None:
        print >> self.log, "Fmodel current bcart   : {0:14.2f} {1:5.2f} {2:5.2f} {3:5.2f} {4:5.2f} {5:5.2f}".format(*self.fmodel_current.fmodel_kbu().b_cart())
      if self.fmodel_running is not None:
        print >> self.log, "Fmodel running bcart   : {0:14.2f} {1:5.2f} {2:5.2f} {3:5.2f} {4:5.2f} {5:5.2f}".format(*self.fmodel_running.fmodel_kbu().b_cart())
      if self.fmodel_total  is not None:
        print >> self.log, "Fmodel total bcart     : {0:14.2f} {1:5.2f} {2:5.2f} {3:5.2f} {4:5.2f} {5:5.2f}".format(*self.fmodel_total.fmodel_kbu().b_cart())
      print >> self.log, "|"+"-"*77+"|\n"
  
    def write_diff_map_ensemble(self, out):
      crystal_symmetry = self.er_data.xray_structures_diff_map[0].crystal_symmetry()
      print >> out, pdb.format_cryst1_record(crystal_symmetry = crystal_symmetry)
      print >> out, pdb.format_scale_records(unit_cell = crystal_symmetry.unit_cell())
      for n,xrs in enumerate(self.er_data.xray_structures_diff_map):
        print >> out, "MODEL %8d"%(n+1)
        print >> out, xrs.as_pdb_file()
        print >> out, "ENDMDL"
      print >> out, "END"
  
    def update_single_hierarchy(self, i_model):
      xrs = self.er_data.xray_structures[i_model]
      scatterers = xrs.scatterers()
      sites_cart = xrs.sites_cart()
      u_isos = xrs.extract_u_iso_or_u_equiv()
      occupancies = scatterers.extract_occupancies()
      u_carts = scatterers.extract_u_cart_plus_u_iso(xrs.unit_cell())
      scat_types = scatterers.extract_scattering_types()
      i_model_pdb_hierarchy = self.er_data.pdb_hierarchys[i_model]
      pdb_atoms = i_model_pdb_hierarchy.atoms()
      i_model_ke = self.er_data.ke_pdb[i_model]
      for j_seq, atom in enumerate(pdb_atoms):
        if j_seq < len(sites_cart):
          atom.xyz = sites_cart[j_seq]
          if self.params.output_running_kinetic_energy_in_occupancy_column:
            #XXX * 0.1 to fit in occ col
            atom.occ = 0.1 * i_model_ke[j_seq]
          else:
            atom.occ = 1.0 / len(self.er_data.xray_structures)
          atom.b = adptbx.u_as_b(u_isos[j_seq])
          e = scat_types[j_seq]
          if (len(e) > 1 and "+-0123456789".find(e[1]) >= 0):
            atom.element = "%2s" % e[:1]
            atom.charge = "%-2s" % e[1:]
          elif (len(e) > 2):
            atom.element = "%2s" % e[:2]
            atom.charge = "%-2s" % e[2:]
          else:
            atom.element = "%2s" % e
            atom.charge = "  "
          if (scatterers[j_seq].flags.use_u_aniso()):
            atom.uij = u_carts[j_seq]
          elif(False):
            atom.uij = self.u_cart
          else:
            atom.uij = (-1,-1,-1,-1,-1,-1)
      return i_model_pdb_hierarchy
  
    def write_ensemble_pdb(self, out):
      crystal_symmetry = self.er_data.xray_structures[0].crystal_symmetry()
      pr = "REMARK   3"
      print >> out, pr
      print >> out,  "REMARK   3 TIME-AVERAGED ENSEMBLE REFINEMENT."
      from phenix import phenix_info # FIXME ???
      ver, tag = phenix_info.version_and_release_tag(f = out)
      if(ver is None):
        prog = "   PROGRAM     : PHENIX (phenix.ensemble_refinement)"
      else:
        if(tag is not None):
          ver = ver+"_"+tag
        prog = "   PROGRAM     : PHENIX (phenix.ensemble_refinement: %s)"%ver
      print >> out,pr+prog
      authors = phenix_info.phenix_developers_last
      l = pr+"   AUTHORS     :"
      j = 0
      i = j
      n = len(l) + 1
      while (j != len(authors)):
        a = len(authors[j]) + 1
        if (n+a > 79):
          print >> out,l, ",".join(authors[i:j]) + ","
          l = pr+"               :"
          i = j
          n = len(l) + 1
        n += a
        j += 1
      if (i != j):
        print >> out,l, ",".join(authors[i:j])
      fmodel_info = self.fmodel_total.info()
      fmodel_info.show_remark_3(out = out)
  #    model_stats = mmtbx.model_statistics.model(model     = self.model,
  #                                               ignore_hd = False)
  #    # set mode_stats.geometry to None as refers to final structure NOT ensemble
  #    model_stats.geometry = None
  #    model_stats.show(out = out, pdb_deposition =True)
      # get mean geometry stats for ensemble
      self.final_geometry_pdb_string = self.ensemble_utils.ensemble_mean_geometry_stats(
          restraints_manager       = self.model.restraints_manager,
          xray_structure           = self.model.get_xray_structure(),
          ensemble_xray_structures = self.er_data.xray_structures,
          ignore_hd                = True,
          verbose                  = False,
          out                      = self.log,
          return_pdb_string        = True)
      print >> out, self.final_geometry_pdb_string
      print >> out, pdb.format_cryst1_record(crystal_symmetry = crystal_symmetry)
      print >> out, pdb.format_scale_records(unit_cell = crystal_symmetry.unit_cell())
      atoms_reset_serial = True
      #
      cntr = 0
      assert len(self.er_data.ke_pdb) == len(self.er_data.xray_structures)
      assert len(self.er_data.pdb_hierarchys) == len(self.er_data.xray_structures)
      for i_model, xrs in enumerate(self.er_data.xray_structures):
        cntr += 1
        print >> out, "MODEL %8d"%cntr
        i_model_pdb_hierarchy = self.update_single_hierarchy(i_model)
        if (atoms_reset_serial):
          atoms_reset_serial_first_value = 1
        else:
          atoms_reset_serial_first_value = None
        out.write(i_model_pdb_hierarchy.as_pdb_string(
          append_end=False,
          atoms_reset_serial_first_value=atoms_reset_serial_first_value))
        #
        print >> out, "ENDMDL"
      print >> out, "END"
  
    def print_ml_stats(self):
      if self.fmodel_running.set_sigmaa is not None:
        self.run_time_stats_dict.update({'Sigma_a':self.fmodel_running.set_sigmaa})
      if self.params.target_name in ['ml', 'mlhl'] :
        self.run_time_stats_dict.update({'Alpha':self.fmodel_running.alpha_beta()[0].data()})
        self.run_time_stats_dict.update({'Beta':self.fmodel_running.alpha_beta()[1].data()})
      if self.fmodel_running.n_obs is not None:
        self.run_time_stats_dict.update({'Eobs(fixed)':self.fmodel_running.n_obs})
      if self.fmodel_running.n_calc is not None:
        self.run_time_stats_dict.update({'Ecalc(fixed)':self.fmodel_running.n_calc})
  
      make_header("ML statistics", out = self.log)
      print >> self.log, '  {0:<23}: {1:>12} {2:>12} {3:>12}'.format('','min','max','mean')
      for key in sorted(self.run_time_stats_dict.keys()):
        info = self.run_time_stats_dict[key].min_max_mean()
        print >> self.log, '  {0:<23}: {1:12.3f} {2:12.3f} {3:12.3f}'.format(
          key,
          info.min,
          info.max,
          info.mean)
      print >> self.log, "|"+"-"*77+"|\n"
  
  ################################################################################
  
  def show_data(fmodel, n_outl, test_flag_value, f_obs_labels, log):
    info = fmodel.info()
    flags_pc = \
     fmodel.r_free_flags().data().count(True)*1./fmodel.r_free_flags().data().size()
    print >> log, "Data statistics"
    try: f_obs_labels = f_obs_labels[:f_obs_labels.index(",")]
    except ValueError: pass
    result = " \n    ".join([
      "data_label                          : %s" % f_obs_labels,
      "high_resolution                     : "+format_value("%-5.2f",info.d_min),
      "low_resolution                      : "+format_value("%-6.2f",info.d_max),
      "completeness_in_range               : " + \
        format_value("%-6.2f",info.completeness_in_range),
      "completeness(d_min-inf)             : " + \
        format_value("%-6.2f",info.completeness_d_min_inf),
      "completeness(6A-inf)                : " + \
        format_value("%-6.2f",info.completeness_6_inf),
      "wilson_b                            : " + \
        format_value("%-6.1f",fmodel.wilson_b()),
      "number_of_reflections               : " + \
        format_value("%-8d",  info.number_of_reflections),
      "test_set_size                       : " + \
        format_value("%-8.4f",flags_pc),
      "test_flag_value                     : " + \
        format_value("%-d",   test_flag_value),
      "number_of_Fobs_outliers             : " + format_value("%-8d",  n_outl),
      "anomalous_flag                      : " + \
        format_value("%-6s",  fmodel.f_obs().anomalous_flag())])
    print >> log, "   ", result
  
  def show_model_vs_data(fmodel, log):
    d_max, d_min = fmodel.f_obs().d_max_min()
    flags_pc = fmodel.r_free_flags().data().count(True)*100./\
      fmodel.r_free_flags().data().size()
    if(flags_pc == 0): r_free = None
    else: r_free = fmodel.r_free()
    k_sol = format_value("%-5.2f",fmodel.fmodel_kbu().k_sols()[0])
    b_sol = format_value("%-7.2f",fmodel.fmodel_kbu().b_sols()[0])
    b_cart = " ".join([("%8.2f"%v).strip() for v in fmodel.fmodel_kbu().b_cart()])
    print >> log, "Model vs data statistics"
    result = " \n    ".join([
      "r_work(re-computed)                 : " + \
        format_value("%-6.4f",fmodel.r_work()),
      "r_free(re-computed)                 : " + format_value("%-6.4f",r_free),
      "scale_k1                            : " + \
        format_value("%-6.4f",fmodel.scale_k1()),
      "bulk_solvent_(k_sol,b_sol)          : %s%s" % (k_sol, b_sol),
      "overall_anisotropic_scale_(b_cart)  : " + format_value("%-s",b_cart)])
    print >> log, "   ", result
  
  def write_mtz_file(fmodel_total, raw_data, raw_flags, prefix, params):
    assert (fmodel_total is not None)
    class labels_decorator:
      def __init__(self, amplitudes_label, phases_label):
        self._amplitudes = amplitudes_label
        self._phases = phases_label
      def amplitudes(self):
        return self._amplitudes
      def phases(self, root_label, anomalous_sign=None):
        assert anomalous_sign is None or not anomalous_sign
        return self._phases
    xray_suffix = "_xray"
    f_obs_label = "F-obs"
    i_obs_label = "I-obs"
    flag_label = "R-free-flags"
    if (raw_data.is_xray_intensity_array()):
      column_root_label = i_obs_label
    else:
      column_root_label = f_obs_label
    mtz_dataset_original = raw_data.as_mtz_dataset(
      column_root_label=column_root_label)
    mtz_dataset_original.add_miller_array(
      miller_array = raw_flags,
      column_root_label=flag_label)
    mtz_dataset_original.set_name("Original-experimental-data")
    new_dataset = mtz_dataset_original.mtz_crystal().add_dataset(
      name = "Experimental-data-used-in-refinement", wavelength=1)
    new_dataset.add_miller_array(
      miller_array = fmodel_total.f_obs(),
      column_root_label="F-obs-filtered"+xray_suffix)
    another_dataset = new_dataset.mtz_crystal().add_dataset(
      name = "Model-structure-factors-(bulk-solvent-and-all-scales-included)",
      wavelength=1)
    another_dataset.add_miller_array(
      miller_array = fmodel_total.f_model_scaled_with_k1_composite_work_free(),
      column_root_label="F-model"+xray_suffix)
    yet_another_dataset = another_dataset.mtz_crystal().add_dataset(
      name = "Fourier-map-coefficients", wavelength=1)
    cmo = mmtbx.maps.compute_map_coefficients(
      fmodel = fmodel_total,
      params = params.electron_density_maps.map_coefficients)
    for ma in cmo.mtz_dataset.mtz_object().as_miller_arrays():
      labels=ma.info().labels
      ld = labels_decorator(amplitudes_label=labels[0], phases_label=labels[1])
      yet_another_dataset.add_miller_array(
        miller_array      = ma,
        column_root_label = labels[0],
        label_decorator   = ld)
    yet_another_dataset.mtz_object().write(
      file_name=prefix+".mtz")
    return prefix + ".mtz"
  
  #-----------------------------------------------------------------------
  def run(args, command_name = "phenix.ensemble_refinement", out=None,
      validate=False, replace_stderr=True):
    if(len(args) == 0): args = ["--help"]
    command_line = (iotbx_option_parser(
      usage="%s reflection_file pdb_file [options]" % command_name,
      description='Example: %s data.mtz model.pdb' % command_name
    ).enable_dry_run().enable_show_defaults()).process(args=args)
    if (out is None):
      out = sys.stdout
    if(command_line.expert_level is not None):
      master_params.show(
        expert_level=command_line.expert_level,
        attributes_level=command_line.attributes_level,
        out=out)
      return
    inputs = mmtbx.command_line.load_model_and_data(
      args=command_line.args,
      master_phil=master_params,
      out=out,
      create_fmodel=False,
      process_pdb_file=False)
    working_phil = inputs.working_phil
    params = working_phil.extract()
    if (params.extra_restraints_file is not None):
      # XXX this is a revolting hack...
      print >> out, "Processing custom geometry restraints in file:"
      print >> out, "  %s" % params.extra_restraints_file
      restraints_phil = iotbx.phil.parse(file_name=params.extra_restraints_file)
      cleanup_phil = iotbx.phil.parse("extra_restraints_file=None")
      working_phil = master_params.fetch(
        sources=[working_phil, restraints_phil, cleanup_phil])
      params = working_phil.extract()
    er_params = params.ensemble_refinement
  
    if er_params.electron_density_maps.apply_default_maps != False\
      or len(er_params.electron_density_maps.map_coefficients) == 0:
      maps_par = libtbx.env.find_in_repositories(
        relative_path=\
          "cctbx_project/mmtbx/refinement/ensemble_refinement/maps.params",
        test=os.path.isfile)
      maps_par_phil = iotbx.phil.parse(file_name=maps_par)
      working_params = mmtbx.refinement.ensemble_refinement.master_params.fetch(
                          sources = [working_phil]+[maps_par_phil])
      er_params = working_params.extract().ensemble_refinement
  
    if er_params.output_file_prefix == None:
      er_params.output_file_prefix = os.path.splitext(
        os.path.basename(inputs.pdb_file_names[0]))[0] + "_ensemble"
    log = multi_out()
    log.register(label="stdout", file_object=out)
    log.register(
      label="log_buffer",
      file_object=StringIO(),
      atexit_send_to=None)
    if (replace_stderr):
      sys.stderr = log
    log_file = open(er_params.output_file_prefix+'.log', "w")
    log.replace_stringio(
        old_label="log_buffer",
        new_label="log",
        new_file_object=log_file)
    timer = user_plus_sys_time()
    mmtbx.utils.print_programs_start_header(log=log, text=command_name)
    make_header("Ensemble refinement parameters", out = log)
    working_phil.show(out = log)
    make_header("Model and data statistics", out = log)
    print >> log, "Data file                               : %s" % \
      format_value("%5s", os.path.basename(params.input.xray_data.file_name))
    print >> log, "Model file                              : %s \n" % \
      (format_value("%5s",os.path.basename(inputs.pdb_file_names[0])))
    print >> log, "\nTLS MUST BE IN ATOM RECORDS OF INPUT PDB\n"
    f_obs = inputs.f_obs
    number_of_reflections = f_obs.indices().size()
  
    r_free_flags = inputs.r_free_flags
    raw_flags = inputs.raw_flags
    raw_data = inputs.raw_data
  
    print >> log, "\nPDB file name : ", inputs.pdb_file_names[0]
  
    # Process PDB file
    cif_objects = inputs.cif_objects
    pdb_file = inputs.pdb_file_names[0]
    pdb_ip = mmtbx.model.manager.get_default_pdb_interpretation_params()
    pdb_ip.pdb_interpretation.clash_guard.nonbonded_distance_threshold = -1.0
    pdb_ip.pdb_interpretation.clash_guard.max_number_of_distances_below_threshold = 100000000
    pdb_ip.pdb_interpretation.clash_guard.max_fraction_of_distances_below_threshold = 1.0
    pdb_ip.pdb_interpretation.proceed_with_excessive_length_bonds=True
  
    # Model
    pdb_inp = iotbx.pdb.input(file_name=pdb_file)
    model = mmtbx.model.manager(
      model_input = pdb_inp,
      restraint_objects = cif_objects,
      pdb_interpretation_params = pdb_ip,
      log = log)
    if model.get_number_of_models() > 1:
      raise Sorry("Multiple models not supported.")
  
    # Remove alternative conformations if present
    n_removed_atoms = 0
    if er_params.remove_alt_conf_from_input_pdb:
      n_removed_atoms = model.remove_alternative_conformations(
          always_keep_one_conformer=True)
  
    if n_removed_atoms > 0:
      pdb_file_removed_alt_confs = pdb_file[0:-4]+'_removed_alt_confs.pdb'
      print >> log, "\nRemoving alternative conformations"
      print >> log, "All occupancies reset to 1.0"
      print >> log, "New PDB : ", pdb_file_removed_alt_confs, "\n"
      pdb_str = model.model_as_pdb()
      f = open(pdb_file_removed_alt_confs, 'w')
      f.write(pdb_str)
      f.close()
      pdb_inp = iotbx.pdb.input(file_name=pdb_file_removed_alt_confs)
      model = mmtbx.model.manager(
        model_input = pdb_inp,
        restraint_objects = cif_objects,
        pdb_interpretation_params = pdb_ip,
        log = log)
  
    # Refinement flags
    class rf:
      def __init__(self, size):
        self.individual_sites     = True
        self.individual_adp       = False
        self.sites_individual     = flex.bool(size, True)
        self.sites_torsion_angles = None
        self.torsion_angles       = None
        self.adp_individual_iso   = None
        self.adp_individual_aniso = None
      def inflate(self, **keywords): pass
  
    refinement_flags = rf(size = model.get_number_of_atoms())
  
    model.set_refinement_flags(refinement_flags)
    model.get_restraints_manager()
  
    # Geometry file
    xray_structure = model.get_xray_structure()
    sites_cart = xray_structure.sites_cart()
    site_labels = xray_structure.scatterers().extract_labels()
    model.restraints_manager.geometry.show_sorted(
      sites_cart=sites_cart,
      site_labels=site_labels,
      f=open(er_params.output_file_prefix+'.geo','w') )
  
    print >> log, "Unit cell                               :", f_obs.unit_cell()
    print >> log, "Space group                             :", \
      f_obs.crystal_symmetry().space_group_info().symbol_and_number()
    print >> log, "Number of symmetry operators            :", \
      f_obs.crystal_symmetry().space_group_info().type().group().order_z()
    print >> log, "Unit cell volume                        : %-15.4f" % \
      f_obs.unit_cell().volume()
    f_obs_labels = f_obs.info().label_string()
  
    if (command_line.options.dry_run):
      return None
  
    fmodel = mmtbx.utils.fmodel_simple(
      f_obs                      = f_obs,
      xray_structures            = [model.get_xray_structure()],
      scattering_table           = "n_gaussian",
      r_free_flags               = r_free_flags,
      target_name                = er_params.target_name,
      bulk_solvent_and_scaling   = False,
      bss_params                 = None,
      mask_params                = None,
      twin_laws                  = None,
      skip_twin_detection        = True,
      twin_switch_tolerance      = 2.0,
      outliers_rejection         = True,
      bulk_solvent_correction    = True,
      anisotropic_scaling        = True,
      log                        = log)
    hl_coeffs = inputs.hl_coeffs
    if (hl_coeffs is not None) and (params.input.use_experimental_phases):
      print >> log, "Using MLHL target with experimental phases"
      er_params.target_name = "mlhl"
      hl_coeffs = hl_coeffs.common_set(other=fmodel.f_obs())
    else :
      hl_coeffs = None
    # XXX is this intentional?
    fmodel = mmtbx.f_model.manager(
      mask_params                  = er_params.mask,
      xray_structure               = model.get_xray_structure(),
      f_obs                        = fmodel.f_obs(),
      r_free_flags                 = fmodel.r_free_flags(),
      target_name                  = er_params.target_name,
      abcd                         = hl_coeffs)
    hd_sel = model.get_hd_selection()
    model.get_xray_structure().set_occupancies(
          value     = 1.0,
          selection = hd_sel)
    model.show_occupancy_statistics(out = log)
  
    fmodel.update_xray_structure(
      xray_structure      = model.get_xray_structure(),
      update_f_calc       = True,
      update_f_mask       = False,
      force_update_f_mask = False)
  
    n_outl = f_obs.data().size() - fmodel.f_obs().data().size()
    show_data(fmodel          = fmodel,
              n_outl          = n_outl,
              test_flag_value = inputs.test_flag_value,
              f_obs_labels    = f_obs_labels,
              log             = log)
    show_model_vs_data(fmodel = fmodel,
                       log    = log)
  
    best_trial = None
    if (len(er_params.ptls) == 1):
      best_trial = run_wrapper(
        fmodel               = fmodel,
        model                = model,
        er_params            = er_params,
        raw_data             = raw_data,
        raw_flags            = raw_flags,
        log                  = log).__call__(
          ptls=er_params.ptls[0],
          buffer_output=False,
          append_ptls=False)
    else :
      driver = run_wrapper(
        fmodel               = fmodel,
        model                = model,
        er_params            = er_params,
        raw_data             = raw_data,
        raw_flags            = raw_flags,
        log                  = log)
      trials = []
      if (er_params.nproc in [1, None]) or (sys.platform == "win32"):
        for ptls in er_params.ptls :
          make_header("Running with pTLS = %g" % ptls, out=log)
          trial_result = driver(ptls, buffer_output=False, write_log=False)
          assert (trial_result is not None)
          trials.append(trial_result)
      else :
        trials = easy_mp.pool_map(
          fixed_func=driver,
          args=er_params.ptls,
          processes=er_params.nproc)
      assert (not None in trials)
      best_trial = min(trials, key=lambda t: t.r_free)
      best_trial.save_final(er_params.output_file_prefix)
  
    show_total_time(out = log)
    return result(
      best_trial=best_trial,
      prefix=er_params.output_file_prefix,
      validate=validate,
      log=log)
  
  class run_wrapper(object):
    def __init__(self, model, fmodel, raw_data, raw_flags, er_params, log):
      adopt_init_args(self, locals())
  
    def __call__(self, ptls, buffer_output=True, write_log=True,
        append_ptls=True):
      out = self.log
      log_out = None
      if (buffer_output):
        out = StringIO()
      run_number = None
      if (append_ptls):
        run_number = ptls
      ensemble_refinement = run_ensemble_refinement(
        fmodel               = self.fmodel.deep_copy(),
        model                = self.model.deep_copy(),
        params               = self.er_params,
        raw_data             = self.raw_data,
        raw_flags            = self.raw_flags,
        run_number           = run_number,
        ptls                 = ptls,
        log                  = out)
      if (buffer_output):
        log_out = out.getvalue()
        if (write_log):
          log_file_name = self.er_params.output_file_prefix + '_ptls-' + \
            str(ptls) + '.log'
          log_file = open(log_file_name, 'w')
          log_file.write(log_out)
      return trial(
        ptls=ptls,
        r_work=ensemble_refinement.fmodel_total.r_work(),
        r_free=ensemble_refinement.fmodel_total.r_free(),
        pdb_file=ensemble_refinement.pdb_file,
        mtz_file=ensemble_refinement.mtz_file,
        log_out=log_out,
        number_of_models=len(ensemble_refinement.er_data.xray_structures))
  
  class trial(slots_getstate_setstate):
    __slots__ = ["r_work", "r_free", "pdb_file", "mtz_file", "number_of_models",
                 "log_out", "ptls"]
    def __init__(self, **kwds):
      kwds = dict(kwds)
      for name in self.__slots__ :
        setattr(self, name, kwds[name])
  
    def save_final(self, prefix):
      pdb_out = prefix + ".pdb"
      if (self.pdb_file.endswith(".gz")):
        pdb_out += ".gz"
      os.rename(self.pdb_file, pdb_out)
      os.rename(self.mtz_file, prefix + ".mtz")
      self.pdb_file = pdb_out
      self.mtz_file = prefix + ".mtz"
  
  ########################################################################
  # Phenix GUI hooks
  class result(slots_getstate_setstate):
    __slots__ = [
      "directory", "r_work", "r_free",
      "number_of_models", "pdb_file", "mtz_file","validation", "chi_angles"
    ]
    def __init__(self,
        best_trial,
        prefix,
        log,
        validate=False):
      for attr in ["r_work", "r_free", "number_of_models", "pdb_file",
                   "mtz_file"] :
        setattr(self, attr, getattr(best_trial, attr))
      self.directory = os.getcwd()
      self.validation = None
      if (validate):
        from mmtbx.command_line import validation_summary
        self.validation = validation_summary.run(
          args=[self.pdb_file],
          out=log)
        assert (type(self.validation).__name__ == 'ensemble')
  
        # calculate chi angles for each model
        # each model is assumed to have the same number of protein residues
        self.chi_angles = list()
        hierarchy = pdb.input(self.pdb_file).construct_hierarchy()
        for model in hierarchy.models():
          self.chi_angles.append(calculate_chi_angles(model))
  
    def get_result_files(self, output_dir=None):
      if (output_dir is None):
        output_dir = self.directory
      return (os.path.join(self.directory, self.pdb_file),
              os.path.join(self.directory, self.mtz_file))
  
    def finish_job(self):
      pdb_file, mtz_file = self.get_result_files()
      return (
        [(pdb_file, "Final ensemble"),
         (mtz_file, "Map coefficients")],
        [("R-work", "%.4f" % self.r_work),
         ("R-free", "%.4f" % self.r_free),
         ("Models", str(self.number_of_models)),]
      )
  
  class launcher(runtime_utils.target_with_save_result):
    def run(self):
      os.mkdir(self.output_dir)
      os.chdir(self.output_dir)
      return run(args=list(self.args),
        out=sys.stdout,
        validate=True)
  
  def validate_params(params):
    mmtbx.command_line.validate_input_params(params)
    if (params.ensemble_refinement.ptls is None):
      raise Sorry("You must specify a fraction of atoms to use for TLS fitting.")
    elif (len(params.input.xray_data.labels[0].split(",")) > 2):
      raise Sorry("Anomalous data are not allowed in this program.")
  
    # check files
    mmtbx.command_line.check_files(
      params.input.pdb.file_name, 'pdb',
      'Please provide a valid structure for the Input model.')
    mmtbx.command_line.check_files(
      params.input.monomers.file_name, 'cif',
      'Please provide valid restraints')
    mmtbx.command_line.check_files(
      params.extra_restraints_file, 'phil',
      'Please provide a valid file containing custom restraints.')
  
    return params
  
  # =============================================================================
  from mmtbx.validation import rotalyze
  from mmtbx.rotamer import sidechain_angles
  
  def calculate_chi_angles(model=None):
    '''
  
    =============================================================================
    Function for calculating all dihedral angles for each protein residue in a
    model.
  
    Parameters:
    -----------
    model - from hierarchy.models()
  
    Return:
    -------
    dictionary - where 'id_str' is a list containing residue ids
                 (atom_group.id_str()) and 'chi_angles' is a list containing
                 a list of dihedral angles for each residue (list of lists)
  
    Notes:
    ------
    The list for an amino acid is of size 0 if it has no defined dihedral angles.
    A dihedral angle in the list for an atom_group can be None if there are
    missing atoms. A protein residue is defined as common_amino_acid.
  
    =============================================================================
  
    '''
  
    id_str = list()
    chi_angles = list()
    xyz = list()
  
    if (model is not None):
      get_class = pdb.common_residue_names_get_class
      angles = sidechain_angles.SidechainAngles(False)
  
      # loop over all atom groups in chain
      for chain in model.chains():
        for rg in chain.residue_groups():
          all_dict = rotalyze.construct_complete_sidechain(rg)
          for ag in rg.atom_groups():
            residue_class = get_class(ag.resname)
            if (residue_class == 'common_amino_acid'):
              atom_dict = all_dict.get(ag.altloc)
              id_str.append(ag.id_str())
              xyz.append(ag.atoms().extract_xyz().mean())
              chi_angles.append(angles.measureChiAngles(
                res=ag,atom_dict=atom_dict))
  
    return { 'id_str': id_str,
             'xyz': xyz,
             'chi_angles': chi_angles }
  