#!/usr/bin/env python # XDSCorrect.py # Copyright (C) 2006 CCLRC, Graeme Winter # # This code is distributed under the BSD license, a copy of which is # included in the root directory of this package. # # A wrapper to handle the JOB=CORRECT module in XDS. # # 04/JAN/07 FIXME - need to know if we want anomalous pairs separating # in this module... import os import sys import shutil import math from cctbx.uctbx import unit_cell if not os.environ.has_key('XIA2CORE_ROOT'): raise RuntimeError, 'XIA2CORE_ROOT not defined' if not os.environ.has_key('XIA2_ROOT'): raise RuntimeError, 'XIA2_ROOT not defined' if not os.path.join(os.environ['XIA2CORE_ROOT'], 'Python') in sys.path: sys.path.append(os.path.join(os.environ['XIA2CORE_ROOT'], 'Python')) if not os.environ['XIA2_ROOT'] in sys.path: sys.path.append(os.environ['XIA2_ROOT']) from Driver.DriverFactory import DriverFactory # interfaces that this inherits from ... from Schema.Interfaces.FrameProcessor import FrameProcessor # generic helper stuff from XDS import imageset_to_xds, xds_check_version_supported, xds_check_error # specific helper stuff from XDSCorrectHelpers import _parse_correct_lp from Experts.ResolutionExperts import xds_integrate_hkl_to_list, \ bin_o_tron, digest # global flags from Handlers.Flags import Flags from Handlers.Streams import Debug from libtbx.phil import parse master_params = parse(""" refine = *DISTANCE *BEAM *AXIS *ORIENTATION *CELL .type = choice(multi = True) .help = 'what to refine in the CORRECT step' """) def XDSCorrect(DriverType = None, params=None): DriverInstance = DriverFactory.Driver(DriverType) class XDSCorrectWrapper(DriverInstance.__class__, FrameProcessor): '''A wrapper for wrapping XDS in correct mode.''' def __init__(self, params=None): # set up the object ancestors... DriverInstance.__class__.__init__(self) FrameProcessor.__init__(self) # phil parameters if not params: params = master_params.extract() self._params = params # now set myself up... self._parallel = Flags.get_parallel() self.set_cpu_threads(self._parallel) if self._parallel <= 1: self.set_executable('xds') else: self.set_executable('xds_par') # generic bits self._data_range = (0, 0) self._spot_range = [] self._background_range = (0, 0) self._resolution_range = (0, 0) self._resolution_high = 0.0 self._resolution_low = 40.0 # specific information self._cell = None self._spacegroup_number = None self._anomalous = False self._polarization = 0.0 self._reindex_matrix = None self._reindex_used = None self._input_data_files = { } self._output_data_files = { } self._input_data_files_list = [] self._output_data_files_list = ['GXPARM.XDS'] self._ice = 0 self._excluded_regions = [] # the following input files are also required: # # INTEGRATE.HKL # REMOVE.HKL # # and XDS_ASCII.HKL is produced. # in self._integrate_hkl = None self._remove_hkl = None # out self._xds_ascii_hkl = None self._results = None self._remove = [] return # getter and setter for input / output data def set_anomalous(self, anomalous): self._anomalous = anomalous return def get_anomalous(self): return self._anomalous def set_spacegroup_number(self, spacegroup_number): self._spacegroup_number = spacegroup_number return def set_cell(self, cell): self._cell = cell return def set_ice(self, ice): self._ice = ice return def set_excluded_regions(self, excluded_regions): self._excluded_regions = excluded_regions return def set_polarization(self, polarization): if polarization > 1.0 or polarization < 0.0: raise RuntimeError, 'bad value for polarization: %.2f' % \ polarization self._polarization = polarization def set_reindex_matrix(self, reindex_matrix): if not len(reindex_matrix) == 12: raise RuntimeError, 'reindex matrix must be 12 numbers' self._reindex_matrix = reindex_matrix return def get_reindex_used(self): return self._reindex_used def set_resolution_high(self, resolution_high): self._resolution_high = resolution_high def set_resolution_low(self, resolution_low): self._resolution_low = resolution_low def set_input_data_file(self, name, data): self._input_data_files[name] = data return def get_output_data_file(self, name): return self._output_data_files[name] def set_integrate_hkl(self, integrate_hkl): self._integrate_hkl = integrate_hkl return def set_remove_hkl(self, remove_hkl): self._remove_hkl = remove_hkl return def get_remove(self): return self._remove def get_xds_ascii_hkl(self): return self._xds_ascii_hkl # this needs setting up from setup_from_image in FrameProcessor def set_data_range(self, start, end): self._data_range = (start, end) def add_spot_range(self, start, end): self._spot_range.append((start, end)) def set_background_range(self, start, end): self._background_range = (start, end) def get_result(self, name): if not self._results: raise RuntimeError, 'no results' if not self._results.has_key(name): raise RuntimeError, 'result name "%s" unknown' % name return self._results[name] def run(self): '''Run correct.''' # this is ok... # if not self._cell: # raise RuntimeError, 'cell not set' # if not self._spacegroup_number: # raise RuntimeError, 'spacegroup not set' #image_header = self.get_header() ## crank through the header dictionary and replace incorrect ## information with updated values through the indexer ## interface if available... ## need to add distance, wavelength - that should be enough... #if self.get_distance(): #image_header['distance'] = self.get_distance() #if self.get_wavelength(): #image_header['wavelength'] = self.get_wavelength() #if self.get_two_theta(): #image_header['two_theta'] = self.get_two_theta() header = imageset_to_xds(self.get_imageset()) xds_inp = open(os.path.join(self.get_working_directory(), 'XDS.INP'), 'w') # what are we doing? xds_inp.write('JOB=CORRECT\n') xds_inp.write('MAXIMUM_NUMBER_OF_PROCESSORS=%d\n' % \ self._parallel) # check to see if we are excluding ice rings if self._ice != 0: Debug.write('Excluding ice rings') for record in open(os.path.join( os.environ['XIA2_ROOT'], 'Data', 'ice-rings.dat')).readlines(): resol = tuple(map(float, record.split()[:2])) xds_inp.write('EXCLUDE_RESOLUTION_RANGE= %.2f %.2f\n' % \ resol) # exclude requested resolution ranges if len(self._excluded_regions) != 0: Debug.write('Excluding regions: %s' % `self._excluded_regions`) for upper, lower in self._excluded_regions: xds_inp.write('EXCLUDE_RESOLUTION_RANGE= %.2f %.2f\n' % \ (upper, lower)) # postrefine everything to give better values to the # next INTEGRATE run xds_inp.write( 'REFINE(CORRECT)=%s\n' %self._params.refine) if self._polarization > 0.0: xds_inp.write('FRACTION_OF_POLARIZATION=%.2f\n' % \ self._polarization) for record in header: xds_inp.write('%s\n' % record) name_template = os.path.join(self.get_directory(), self.get_template().replace('#', '?')) record = 'NAME_TEMPLATE_OF_DATA_FRAMES=%s\n' % \ name_template xds_inp.write(record) xds_inp.write('DATA_RANGE=%d %d\n' % self._data_range) # xds_inp.write('MINIMUM_ZETA=0.1\n') # include the resolution range, perhaps if self._resolution_high or self._resolution_low: xds_inp.write('INCLUDE_RESOLUTION_RANGE=%.2f %.2f\n' % \ (self._resolution_low, self._resolution_high)) if self._anomalous: xds_inp.write('FRIEDEL\'S_LAW=FALSE\n') xds_inp.write('STRICT_ABSORPTION_CORRECTION=TRUE\n') else: xds_inp.write('FRIEDEL\'S_LAW=TRUE\n') if self._spacegroup_number: if not self._cell: raise RuntimeError, \ 'cannot set spacegroup without unit cell' xds_inp.write('SPACE_GROUP_NUMBER=%d\n' % \ self._spacegroup_number) if self._cell: xds_inp.write('UNIT_CELL_CONSTANTS=') xds_inp.write('%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n' % \ tuple(self._cell)) if self._reindex_matrix: xds_inp.write('REIDX=%d %d %d %d %d %d %d %d %d %d %d %d' % \ tuple(map(int, self._reindex_matrix))) xds_inp.close() # copy the input file... shutil.copyfile(os.path.join(self.get_working_directory(), 'XDS.INP'), os.path.join(self.get_working_directory(), '%d_CORRECT.INP' % self.get_xpid())) # write the input data files... for file_name in self._input_data_files_list: src = self._input_data_files[file_name] dst = os.path.join( self.get_working_directory(), file_name) if src != dst: shutil.copyfile(src, dst) self.start() self.close_wait() xds_check_version_supported(self.get_all_output()) xds_check_error(self.get_all_output()) # look for errors # like this perhaps # !!! ERROR !!! ILLEGAL SPACE GROUP NUMBER OR UNIT CELL # copy the LP file shutil.copyfile(os.path.join(self.get_working_directory(), 'CORRECT.LP'), os.path.join(self.get_working_directory(), '%d_CORRECT.LP' % self.get_xpid())) # gather the output files for file in self._output_data_files_list: self._output_data_files[file] = os.path.join( self.get_working_directory(), file) self._xds_ascii_hkl = os.path.join( self.get_working_directory(), 'XDS_ASCII.HKL') # do some parsing of the correct output... self._results = _parse_correct_lp(os.path.join( self.get_working_directory(), 'CORRECT.LP')) # check that the unit cell is comparable to what went in i.e. # the volume is the same to within a factor of 10 (which is # extremely generous and should only spot gross errors) original = unit_cell(self._cell) refined = unit_cell(self._results['cell']) if original.volume() / refined.volume() > 10: raise RuntimeError, 'catastrophic change in unit cell volume' if refined.volume() / original.volume() > 10: raise RuntimeError, 'catastrophic change in unit cell volume' # record reindex operation used for future reference... this # is to trap trac #419 if 'reindex_op' in self._results: format = 'XDS applied reindex:' + 12 * ' %d' Debug.write(format % tuple(self._results['reindex_op'])) self._reindex_used = self._results['reindex_op'] # get the reflections to remove... for line in open(os.path.join( self.get_working_directory(), 'CORRECT.LP'), 'r').readlines(): if '"alien"' in line: h, k, l = tuple(map(int, line.split()[:3])) z = float(line.split()[4]) if not (h, k, l, z) in self._remove: self._remove.append((h, k, l, z)) return return XDSCorrectWrapper() if __name__ == '__main__': correct = XDSCorrect() directory = os.path.join(os.environ['XIA2_ROOT'], 'Data', 'Test', 'Images') correct.setup_from_image(os.path.join(directory, '12287_1_E1_001.img')) correct.set_data_range(1, 1) correct.set_background_range(1, 1) correct.add_spot_range(1, 1) correct.run()