source: trunk/G2importphase_CIF.py @ 469

Last change on this file since 469 was 469, checked in by toby, 11 years ago

rework phase import
File size: 8.6 KB
Line 
1# Routines to import Phase information from CIF files
2import sys
3import random as ran
4import GSASIIIO as G2IO
5import GSASIIspc as G2spc
6import GSASIIlattice as G2lat
7import CifFile as cif # PyCifRW from James Hester
8
9class CIFPhaseReader(G2IO.ImportPhase):
10    def __init__(self):
11        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
12            #extensionlist=('.CIF','.cif'),
13            extensionlist=('.CIF','.cif','.pdb'), # just for test!
14            strictExtension=False,
15            formatName = 'CIF',
16            longFormatName = 'Crystallographic Information File import'
17            )
18    def ContentsValidator(self, filepointer):
19        filepointer.seek(0) # rewind the file pointer
20        for i,line in enumerate(filepointer):
21            if i >= 1000: break
22            ''' Encountered only blank lines or comments in first 1000
23            lines. This is unlikely, but assume it is CIF since we are
24            even less likely to find a file with nothing but hashes and
25            blank lines'''
26            line = line.strip()
27            if len(line) == 0:
28                continue # ignore blank lines
29            elif line.startswith('#'):
30                continue # ignore comments
31            elif line.startswith('data_'):
32                return True
33            else:
34                return False # found something else
35        return True
36    def Reader(self,filename,filepointer, ParentFrame=None):
37        cellitems = (
38            '_cell_length_a','_cell_length_b','_cell_length_c',
39            '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',)
40        reqitems = (
41             '_atom_site_type_symbol',
42             '_atom_site_fract_x',
43             '_atom_site_fract_y',
44             '_atom_site_fract_z',
45            )
46        phasenamefields = (
47            '_chemical_name_common',
48            '_pd_phase_name',
49            '_chemical_formula_sum'
50            )
51        try:
52#### development code (to speed testing)
53#            try:
54#                    fp = open(filename+'cP',"rb")
55#                    print("reading from "+filename+'cP')
56#                    cf = cPickle.load(fp)
57#                    fp.close()
58#            except:
59#                    cf = cif.ReadCif(filename)
60#                    fp = open(filename+'cP',"wb")
61#                    cPickle.dump(cf,fp)
62#                    fp.close()
63####
64#### end development code
65            self.ShowBusy() # this can take a while
66            cf = cif.ReadCif(filename)
67            # scan blocks for structural info
68            str_blklist = []
69            for blk in cf.keys():
70                for r in reqitems+cellitems:
71                    if r not in cf[blk].keys():
72                        break
73                else:
74                    str_blklist.append(blk)
75            self.DoneBusy()
76            if not str_blklist:
77                return False            # no blocks with coordinates
78            elif len(str_blklist) == 1: # no choices
79                selblk = 0
80            else:                       # choose from options
81                choice = []
82                for blknm in str_blklist:
83                    choice.append('')
84                    # accumumlate some info about this phase
85                    choice[-1] += blknm + ': '
86                    for i in phasenamefields: # get a name for the phase
87                        name = cf[blknm].get(i).strip()
88                        if name is None or name == '?' or name == '.':
89                            continue
90                        else:
91                            choice[-1] += name.strip()[:20] + ', '
92                            break
93                    na = len(cf[blknm].get("_atom_site_fract_x"))
94                    if na == 1:
95                        choice[-1] += '1 atom'
96                    else:
97                        choice[-1] += ('%d' % nd) + ' atoms'
98                    choice[-1] += ', cell: '
99                    fmt = "%.2f,"
100                    for i,key in enumerate(cellitems):
101                        if i == 3: fmt = "%.f,"
102                        if i == 5: fmt = "%.f"
103                        choice[-1] += fmt % cif.get_number_with_esd(
104                            cf[blknm].get(key))[0]
105                    sg = cf[blknm].get("_symmetry_space_group_name_H-M")
106                    if sg: choice[-1] += ', (' + sg.strip() + ')'
107                selblk = self.PhaseSelector(
108                    choice,
109                    ParentFrame=ParentFrame,
110                    title= 'Select a phase from one the CIF data_ blocks below',
111                    size=(600,100)
112                    )
113                if selblk is None: return False # User pressed cancel
114            blkmm = str_blklist[selblk]
115            blk = cf[str_blklist[selblk]]
116            SpGrp = blk.get("_symmetry_space_group_name_H-M")
117            if SpGrp:
118                 E,SGData = G2spc.SpcGroup(SpGrp)
119            if E:
120                self.warnings += ' ERROR in space group symbol '+SpGrp
121                self.warnings += ' N.B.: make sure spaces separate axial fields in symbol' 
122                self.warnings += G2spc.SGErrors(E)
123            else:
124                self.Phase['General']['SGData'] = SGData
125            # cell parameters
126            cell = []
127            for lbl in (
128                '_cell_length_a','_cell_length_b','_cell_length_c',
129                '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',
130                ):
131                cell.append(cif.get_number_with_esd(blk[lbl])[0])
132            Volume = G2lat.calc_V(G2lat.cell2A(cell))
133            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
134            # read in atoms
135            atomloop = blk.GetLoop('_atom_site_label')
136            atomkeys = [i.lower() for i in atomloop.keys()]
137            if blk.get('_atom_site_aniso_label'):
138                anisoloop = blk.GetLoop('_atom_site_aniso_label')
139                anisokeys = [i.lower() for i in anisoloop.keys()]
140            else:
141                anisoloop = None
142                anisokeys = []
143            self.Phase['Atoms'] = []
144            G2AtomDict = {  '_atom_site_type_symbol' : 1,
145                            '_atom_site_label' : 0,
146                            '_atom_site_fract_x' : 3,
147                            '_atom_site_fract_y' : 4,
148                            '_atom_site_fract_z' : 5,
149                            '_atom_site_occupancy' : 6,
150                            '_atom_site_aniso_u_11' : 11,
151                            '_atom_site_aniso_u_22' : 12,
152                            '_atom_site_aniso_u_33' : 13,
153                            '_atom_site_aniso_u_12' : 14,
154                            '_atom_site_aniso_u_13' : 15,
155                            '_atom_site_aniso_u_23' : 16, }
156            for aitem in atomloop:
157                atomlist = ['','','',0,0,0,1.0,'',0,'I',0.01,0,0,0,0,0,0]
158                atomlist.append(ran.randint(0,sys.maxint)) # add a unique label
159                for val,key in zip(aitem,atomkeys):
160                    col = G2AtomDict.get(key)
161                    if col >= 3:
162                        atomlist[col] = cif.get_number_with_esd(val)[0]
163                    elif col is not None:
164                        atomlist[col] = val
165                    elif key in ('_atom_site_thermal_displace_type',
166                               '_atom_site_adp_type'):   #Iso or Aniso?
167                        if val.lower() == 'uani':
168                            atomlist[9] = 'A'
169                    elif key == '_atom_site_u_iso_or_equiv':
170                        atomlist[10] =cif.get_number_with_esd(val)[0]
171                ulbl = '_atom_site_aniso_label'
172                if  atomlist[9] == 'A' and atomlist[0] in blk.get(ulbl):
173                    for val,key in zip(anisoloop.GetKeyedPacket(ulbl,atomlist[0]),
174                                       anisokeys):
175                        col = G2AtomDict.get(key)
176                        if col:
177                            atomlist[col] = cif.get_number_with_esd(val)[0]
178                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)
179                self.Phase['Atoms'].append(atomlist)
180            for lbl in phasenamefields: # get a name for the phase
181                name = blk.get(lbl).strip()
182                if name is None or name == '?' or name == '.':
183                    continue
184                else:
185                    break
186            else: # no name found, use block name for lack of a better choice
187                name = blknm
188            self.Phase['General']['Name'] = name.strip()[:20]
189            return True
190        except Exception as detail:
191            print 'CIF error:',detail # for testing
192            print sys.exc_info()[0] # for testing
193            return False
194        finally:
195            self.DoneBusy()
196
Note: See TracBrowser for help on using the repository browser.