source: trunk/imports/G2sfact_CIF.py @ 3245

Last change on this file since 3245 was 3245, checked in by vondreele, 4 years ago

fix problem with incommensurate magnetic structures from CIFhklReader
make "Edit Body" as "Edit Residue Body" in RB menu
change protein validation bar plots to be vertical
change 'twin' to 'flag' in Reflection List table & add '-2' as 'free' (for proteins - not yet implemented)
fix PDB phase import to get atom type from 76:78 of ATOM record
change importer names for single crystal TOF data to be more explicit (SNS vs ISIS)
change cif structure factor importer - F2, sig(F2) & Fcalc now allowed

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 12.4 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2018-01-27 14:49:12 +0000 (Sat, 27 Jan 2018) $
4# $Author: vondreele $
5# $Revision: 3245 $
6# $URL: trunk/imports/G2sfact_CIF.py $
7# $Id: G2sfact_CIF.py 3245 2018-01-27 14:49:12Z vondreele $
8########### SVN repository information ###################
9'''
10*Module G2sfact_CIF: CIF import*
11-----------------------------------
12Read structure factors from a CIF reflection table.
13
14'''
15# routines to read in structure factors from a CIF
16#
17from __future__ import division, print_function
18import numpy as np
19import os.path
20import GSASIIobj as G2obj
21import GSASIIIO as G2IO
22import GSASIIpath
23GSASIIpath.SetVersionNumber("$Revision: 3245 $")
24import CifFile as cif # PyCifRW from James Hester
25
26class CIFhklReader(G2obj.ImportStructFactor):
27    'Routines to import Phase information from a CIF file'
28    def __init__(self):
29        super(self.__class__,self).__init__( # fancy way to self-reference
30            extensionlist = ('.CIF','.cif','.FCF','.fcf','.HKL','.hkl'),
31            strictExtension = False,
32            formatName = 'CIF',
33            longFormatName = 'CIF format structure factor file (.cif or .hkl)'
34            )
35    # Validate the contents
36    def ContentsValidator(self, filename):
37        'Use standard CIF validator'
38        fp = open(filename,'r')
39        return self.CIFValidator(fp)
40        fp.close()
41
42    def Reader(self, filename, ParentFrame=None, **kwarg):
43        '''Read single crystal data from a CIF.
44        If multiple datasets are requested, use self.repeat and buffer caching.
45        '''
46        hklitems = [('_refln_index_h','_refln_index_k','_refln_index_l','_refln_index_m_1'),
47                    ('_refln.index_h','_refln.index_k','_refln.index_l','_refln.index_m_1')]
48        cellitems = [
49            ('_cell_length_a','_cell_length_b','_cell_length_c',
50             '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',),
51            ('_cell.length_a','_cell.length_b','_cell.length_c',
52             '_cell.angle_alpha','_cell.angle_beta','_cell.angle_gamma',),]
53
54        Fdatanames = ('_refln_f_meas','_refln.f_meas','_refln.f_meas_au',
55                      )
56       
57        F2datanames = ('_refln_f_squared_meas','_refln.f_squared_meas',
58            '_refln_intensity_meas','_refln.intensity_meas',
59                      )
60
61#        Idatanames = ('_refln_intensity_meas','_refln.intensity_meas',
62#                      ) # not used yet
63#
64#        Isignames = ('_refln_intensity_meas_sigma','_refln.intensity_meas_sigma',
65#                      ) # not used yet
66
67        Fcalcnames = ('_refln_f_calc','_refln.f_calc','_refln.f_calc_au',
68                      )
69       
70        F2calcnames = ('_refln_f_squared_calc','_refln.f_squared_calc',
71                      )
72
73        Fsignames = ('_refln_f_meas_sigma','_refln.f_meas_sigma','_refln.f_meas_sigma_au',
74                    '_refln_f_sigma',
75                      )
76       
77        F2signames = ('_refln_f_squared_meas_sigma','_refln.f_squared_meas_sigma',
78                      '_refln_f_squared_sigma','_refln.f_squared_sigma',
79                      '_refln_intensity_meas_sigma','_refln.intensity_meas_sigma',
80                      '_refln.intensity_sigma',)
81
82        phasenames = ('_refln_phase_calc','_refln.phase_calc',
83                      )
84
85
86        SGdataname = ('_symmetry_space_group_name_H-M', '_symmetry.space_group_name_H-M')
87                     
88        phasenamefields = (
89            '_chemical_name_common',
90            '_pd_phase_name',
91            '_chemical_formula_sum'
92            )
93        rdbuffer = kwarg.get('buffer')
94        cf = None
95        if self.repeat and rdbuffer is not None:
96            cf = rdbuffer.get('lastcif')
97            print ('Reusing previously parsed CIF')
98        if cf is None:
99            cf = G2obj.ReadCIF(filename)
100        # scan blocks for reflections
101        self.errors = 'Error during scan of blocks for datasets'
102        blklist = []
103        for blk in cf.keys(): # scan for reflections, F or F2 values and cell lengths.
104            # Ignore blocks that do not have structure factors and a cell
105            blkkeys = [k.lower() for k in cf[blk].keys()]
106            gotFo = False
107            im = 0
108            for i in range(2):
109                if hklitems[i][3] in blkkeys:   #Super lattice reflections h,k,l,m
110                    im = 1
111                if hklitems[i][0] in blkkeys and hklitems[i][1] in blkkeys and hklitems[i][2] in blkkeys:
112                    dnIndex = i
113                    break
114            else:
115                break # no reflections
116            for dn in Fdatanames: 
117                if dn in blkkeys:
118                    blklist.append(blk)
119                    gotFo = True
120                    break
121            if gotFo: break
122            for dn in F2datanames: 
123                if dn in blkkeys:
124                    blklist.append(blk)
125                    break
126            else:
127                break
128        if not blklist:
129            selblk = None # no block to choose
130        elif len(blklist) == 1: # only one choice
131            selblk = 0
132        else:                       # choose from options
133            choice = []
134            for blknm in blklist:
135                choice.append('')
136                # accumulate some info about this phase
137                choice[-1] += blknm + ': '
138                for i in phasenamefields: # get a name for the phase
139                    name = cf[blknm].get(i)
140                    if name is None or name == '?' or name == '.':
141                        continue
142                    else:
143                        choice[-1] += name.strip()[:20] + ', '
144                        break
145                s = ''
146                fmt = "%.2f,"
147                for i,key in enumerate(cellitems[dnIndex]):
148                    if i == 3: fmt = "%.f,"
149                    if i == 5: fmt = "%.f"
150                    val = cf[blknm].get(key)
151                    if val is None: break
152                    s += fmt % cif.get_number_with_esd(val)[0]
153                if s: choice[-1] += ', cell: ' + s
154                for dn in SGdataname:
155                    sg = cf[blknm].get(dn)
156                    if sg: 
157                        choice[-1] += ', (' + sg.strip() + ')'
158                        break
159            choice.append('Import all of the above')
160            if self.repeat: # we were called to repeat the read
161                selblk = self.repeatcount
162                self.repeatcount += 1
163                if self.repeatcount >= len(blklist): self.repeat = False
164            else:
165                selblk = G2IO.BlockSelector(
166                    choice,
167                    ParentFrame=ParentFrame,
168                    title='Select a dataset from one the CIF data_ blocks below',
169                    size=(600,100),
170                    header='Dataset Selector')
171        self.errors = 'Error during reading of selected block'
172        if selblk is None:
173            return False # no block selected or available
174        if selblk >= len(blklist): # all blocks selected
175            selblk = 0
176            self.repeat = True
177            if rdbuffer is not None:
178                rdbuffer['lastcif'] = cf # save the parsed cif for the next loop
179            self.repeatcount = 1
180        blknm = blklist[selblk]
181        blk = cf[blklist[selblk]]
182        self.objname = os.path.basename(filename)+':'+str(blknm)
183        self.errors = 'Error during reading of reflections'
184        # read in reflections
185        try:
186            refloop = blk.GetLoop(hklitems[0][0])
187            dnIndex = 0
188        except KeyError:
189            try:
190                refloop = blk.GetLoop(hklitems[1][0])
191                dnIndex = 1
192            except KeyError:
193                self.errors += "\nUnexpected: '_refln[-.]index_h not found!"
194                return False
195        itemkeys = {}
196        # prepare an index to the CIF reflection loop
197        for i,key in enumerate(refloop.keys()):
198            itemkeys[key.lower()] = i
199           
200        # scan for obs & sig data names:
201        F2dn = None
202        Fdn = None
203        F2cdn = None
204        Fcdn = None
205        F2sdn = None
206        Fsdn = None
207        Phdn = None
208        for dn in F2datanames:
209            if dn in itemkeys:
210                F2dn = dn
211                for dm in F2signames:
212                    if dm in itemkeys:
213                        F2sdn = dm
214                        break
215                break
216        else:
217            for dn in Fdatanames:
218                if dn in itemkeys:
219                    Fdn = dn
220                    for dm in Fsignames:
221                        if dm in itemkeys:
222                            Fsdn = dm
223                            break
224                    break
225            else:
226                msg = "\nno F or F2 loop value found in file\n"
227                msg += "A CIF reflection file needs to have at least one of\n"
228                for dn in F2datanames+Fdatanames:
229                    msg += dn + ', '
230                self.errors += msg                       
231                return False
232       
233        # scan for calc data names - might be different!
234        for dm in F2calcnames:
235            if dm in itemkeys:
236                F2cdn = dm
237                break
238        for dm in Fcalcnames:
239            if dm in itemkeys:
240                Fcdn = dm
241                break
242
243        for dn in phasenames:
244            if dn in itemkeys:
245                Phdn = dn
246                break
247           
248        # loop over all reflections
249        for item in refloop:
250            F2c = 0.0
251            sigF2 = 0.0
252            HKL = []
253            try:
254                for i in hklitems[dnIndex][:3+im]: # '_refln[._]index_[hkl]'
255                    num = itemkeys.get(i)
256                    try:
257                        HKL.append(int(item[num]))
258                    except:
259                        HKL.append('.')
260                #h,k,l,tw,dsp,Fo2,sig,Fc2,Fot2,Fct2,phase,Ext
261                if im:
262                    ref = HKL+[1,0,0,0,1, 0,0,0,0,0, 0,0] 
263                else:
264                    ref = HKL+[1,0,0,1,0, 0,0,0,0,0, 0] 
265                if F2dn:
266                    F2 = item[itemkeys[F2dn]]
267                    if '(' in F2:
268                        F2, sigF2 = cif.get_number_with_esd(F2)
269                        F2 = float(F2)
270                        sigF2 = float(sigF2)
271                    elif F2sdn:
272                        F2 = float(F2)
273                        sigF2 = float(item[itemkeys[F2sdn]])
274                    else:
275                        F2 = float(F2)
276                else:
277                    F = item[itemkeys[Fdn]]
278                    if '(' in F:
279                        F, sig = cif.get_number_with_esd(F)
280                    elif Fsdn:
281                        F = float(F)
282                        sig = float(item[itemkeys[Fsdn]])
283                    else:
284                        F = float(F)
285                        sig = 0.0
286                    F2 = F**2
287                    sigF2 = 2.0*F*sig
288                   
289                try:
290                    if F2cdn:
291                        F2c = float(item[itemkeys[F2cdn]])
292                except:
293                    pass
294                try:
295                    if Fcdn:
296                        Fc = float(item[itemkeys[Fcdn]])
297                        F2c = Fc*Fc
298                except:
299                    pass
300                           
301                ref[8+im] = F2
302                ref[5+im] = F2
303                ref[6+im] = sigF2
304                ref[9+im] = F2c
305                ref[7+im] = F2c
306                try:
307                    if Phdn:
308                        ref[10+im] = float(item[itemkeys[Phdn]])
309                except:
310                    pass
311            except:
312                continue # skip over incompletely parsed reflections
313            self.RefDict['RefList'].append(ref)
314#                self.RefDict['FF'].append({})
315        self.RefDict['RefList'] = np.array(self.RefDict['RefList'])
316        self.errors = 'Error during reading of dataset parameters'
317        Type = 'SXC'
318        if blk.get('_diffrn_radiation_probe'):
319            if blk['_diffrn_radiation_probe'] == 'neutron':
320                Type = 'SNC'
321        elif blk.get('_diffrn_radiation.probe'):
322            if blk['_diffrn_radiation.probe'] == 'neutron':
323                Type = 'SNC'
324        self.RefDict['Type'] = Type
325        self.RefDict['Super'] = im
326        if blk.get('_diffrn_radiation_wavelength'):
327            wave = float(blk['_diffrn_radiation_wavelength'])
328        elif blk.get('_diffrn_radiation.wavelength'):
329            wave = float(blk['_diffrn_radiation.wavelength'])
330        else:
331            wave = 0.70926
332        self.UpdateParameters(Type=Type,Wave=wave) # histogram type
333        return True
Note: See TracBrowser for help on using the repository browser.