source: trunk/imports/G2phase_ISO.py @ 1147

Last change on this file since 1147 was 1147, checked in by toby, 9 years ago

Complete initial ISODISPLACE implementation; mod. phase initialization; change atom pointer init.; rework parameter display window

  • Property svn:eol-style set to native
File size: 19.0 KB
Line 
1########### SVN repository information ###################
2# $Date: 2012-02-13 11:33:35 -0600 (Mon, 13 Feb 2012) $
3# $Author: vondreele & toby $
4# $Revision: 482 $
5# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/G2importphase_CIF.py $
6# $Id: G2importphase_CIF.py 482 2012-02-13 17:33:35Z vondreele $
7########### SVN repository information ###################
8# Routines to import Phase information from CIF files
9import sys
10import random as ran
11import urllib
12import numpy as np
13import GSASIIIO as G2IO
14import GSASIIobj as G2obj
15import GSASIIspc as G2spc
16import GSASIIlattice as G2lat
17import GSASIIpy3 as G2p3
18import CifFile as cif # PyCifRW from James Hester
19
20class ISODISTORTPhaseReader(G2IO.ImportPhase):
21    varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz'}
22    def __init__(self):
23        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
24            extensionlist=('.CIF','.cif'),
25            strictExtension=False,
26            formatName = 'ISODISTORT CIF',
27            longFormatName = 'CIF from ISODISTORT import'
28            )
29    def ContentsValidator(self, filepointer):
30        filepointer.seek(0) # rewind the file pointer
31        lines = filepointer.read(10000) # scan the first 10K bytes
32        # quick tests for an ISODISTORT output file
33        if lines.find('ISODISTORT') != -1:
34            return True
35        if lines.find('data_isodistort-') != -1:
36            return True
37        if lines.find('_iso_') != -1:
38            return True
39        return False
40    def Reader(self,filename,filepointer, ParentFrame=None, usedRanIdList=[], **unused):
41        import re
42        self.Phase = G2IO.SetNewPhase(Name='new phase',SGData=G2IO.P1SGData) # create a new empty phase dict
43        # make sure the ranId is really unique!
44        while self.Phase['ranId'] in usedRanIdList:
45            self.Phase['ranId'] = ran.randint(0,sys.maxint)
46        returnstat = False
47        cellitems = (
48            '_cell_length_a','_cell_length_b','_cell_length_c',
49            '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',)
50        reqitems = (
51             '_atom_site_type_symbol',
52             '_atom_site_fract_x',
53             '_atom_site_fract_y',
54             '_atom_site_fract_z',
55            )
56        phasenamefields = (
57            '_chemical_name_common',
58            '_pd_phase_name',
59            '_chemical_formula_sum'
60            )
61        try:
62            self.ShowBusy() # this can take a while
63            ciffile = 'file:'+urllib.pathname2url(filename)
64            cf = cif.ReadCif(ciffile)
65            # scan blocks for structural info
66            str_blklist = []
67            for blk in cf.keys():
68                for r in reqitems+cellitems:
69                    if r not in cf[blk].keys():
70                        break
71                else:
72                    str_blklist.append(blk)
73            if not str_blklist:
74                selblk = None # no block to choose
75            elif len(str_blklist) == 1: # only one choice
76                selblk = 0
77            else:                       # choose from options
78                choice = []
79                for blknm in str_blklist:
80                    choice.append('')
81                    # accumumlate some info about this phase
82                    choice[-1] += blknm + ': '
83                    for i in phasenamefields: # get a name for the phase
84                        name = cf[blknm].get(i).strip()
85                        if name is None or name == '?' or name == '.':
86                            continue
87                        else:
88                            choice[-1] += name.strip()[:20] + ', '
89                            break
90                    na = len(cf[blknm].get("_atom_site_fract_x"))
91                    if na == 1:
92                        choice[-1] += '1 atom'
93                    else:
94                        choice[-1] += ('%d' % nd) + ' atoms'
95                    choice[-1] += ', cell: '
96                    fmt = "%.2f,"
97                    for i,key in enumerate(cellitems):
98                        if i == 3: fmt = "%.f,"
99                        if i == 5: fmt = "%.f"
100                        choice[-1] += fmt % cif.get_number_with_esd(
101                            cf[blknm].get(key))[0]
102                    sg = cf[blknm].get("_symmetry_space_group_name_H-M")
103                    if sg: choice[-1] += ', (' + sg.strip() + ')'
104                selblk = self.PhaseSelector(
105                    choice,
106                    ParentFrame=ParentFrame,
107                    title= 'Select a phase from one the CIF data_ blocks below',
108                    size=(600,100)
109                    )
110            if selblk is None:
111                returnstat = False # no block selected or available
112            else:
113                blknm = str_blklist[selblk]
114                blk = cf[str_blklist[selblk]]
115                SpGrp = blk.get("_symmetry_space_group_name_H-M")
116                if SpGrp:
117                     E,SGData = G2spc.SpcGroup(SpGrp)
118                if E:
119                    self.warnings += ' ERROR in space group symbol '+SpGrp
120                    self.warnings += ' N.B.: make sure spaces separate axial fields in symbol' 
121                    self.warnings += G2spc.SGErrors(E)
122                else:
123                    self.Phase['General']['SGData'] = SGData
124                # cell parameters
125                cell = []
126                for lbl in cellitems:
127                    cell.append(cif.get_number_with_esd(blk[lbl])[0])
128                Volume = G2lat.calc_V(G2lat.cell2A(cell))
129                self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
130                # read in atoms
131                atomlbllist = [] # table to look up atom IDs
132                atomloop = blk.GetLoop('_atom_site_label')
133                atomkeys = [i.lower() for i in atomloop.keys()]
134                if blk.get('_atom_site_aniso_label'):
135                    anisoloop = blk.GetLoop('_atom_site_aniso_label')
136                    anisokeys = [i.lower() for i in anisoloop.keys()]
137                else:
138                    anisoloop = None
139                    anisokeys = []
140                self.Phase['Atoms'] = []
141                G2AtomDict = {  '_atom_site_type_symbol' : 1,
142                                '_atom_site_label' : 0,
143                                '_atom_site_fract_x' : 3,
144                                '_atom_site_fract_y' : 4,
145                                '_atom_site_fract_z' : 5,
146                                '_atom_site_occupancy' : 6,
147                                '_atom_site_aniso_u_11' : 11,
148                                '_atom_site_aniso_u_22' : 12,
149                                '_atom_site_aniso_u_33' : 13,
150                                '_atom_site_aniso_u_12' : 14,
151                                '_atom_site_aniso_u_13' : 15,
152                                '_atom_site_aniso_u_23' : 16, }
153                ranIdlookup = {}
154                for aitem in atomloop:
155                    atomlist = ['','','',0,0,0,1.0,'',0,'I',0.01,0,0,0,0,0,0,0]
156                    atomlist[-1] = ran.randint(0,sys.maxint) # add a random Id
157                    while atomlist[-1] in ranIdlookup:
158                        atomlist[-1] = ran.randint(0,sys.maxint) # make it unique                       
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                    ranIdlookup[atomlist[0]] = atomlist[-1]
181                    if atomlist[0] in atomlbllist:
182                        self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
183                        returnstat = False # cannot use this
184                    else:
185                        atomlbllist.append(atomlist[0])
186                if len(atomlbllist) != len(self.Phase['Atoms']):
187                    raise Exception,'Repeated atom labels prevents ISODISTORT decode'
188                for lbl in phasenamefields: # get a name for the phase
189                    name = blk.get(lbl)
190                    if name is None:
191                        continue
192                    name = name.strip()
193                    if name == '?' or name == '.':
194                        continue
195                    else:
196                        break
197                else: # no name found, use block name for lack of a better choice
198                    name = blknm
199                self.Phase['General']['Name'] = name.strip()[:20]
200                #----------------------------------------------------------------------
201                # now read in the ISODISPLACE modes
202                #----------------------------------------------------------------------
203                modelist = []
204                shortmodelist = []
205                for lbl in blk.get('_iso_displacivemode_label'):
206                    modelist.append(lbl)
207                    # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
208                    # where SSSSS is the parent spacegroup, [x,y,z] is a location
209                    regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
210                    # this extracts the AAAAA and BBBBB parts of the string
211                    if regexp:
212                        lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
213                    G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
214                   
215                # read in the coordinate offset variables names and map them to G2 names/objects
216                coordVarLbl = []
217                G2varLbl = []
218                G2varObj = []
219                error = False
220                for lbl in blk.get('_iso_deltacoordinate_label'):
221                    coordVarLbl.append(lbl)
222                    if '_' in lbl:
223                        albl = lbl[:lbl.rfind('_')]
224                        vlbl = lbl[lbl.rfind('_')+1:]
225                    else:
226                        self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
227                        error = True
228                        continue
229                    if albl not in atomlbllist:
230                        self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
231                        error = True
232                        continue
233                    else:
234                        anum = atomlbllist.index(albl)
235                    var = self.varLookup.get(vlbl)
236                    if not var:
237                        self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
238                        error = True
239                        continue
240                    G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
241                    G2varObj.append(G2obj.G2VarObj(
242                        (self.Phase['ranId'],None,var,ranIdlookup[albl])
243                        ))
244                if error:
245                    raise Exception,"Error decoding variable labels"
246                   
247                if len(G2varObj) != len(modelist):
248                    print "non-square input"
249                    raise Exception,"Rank of _iso_displacivemode != _iso_deltacoordinate"
250               
251                error = False
252                ParentCoordinates = {}
253                for lbl,exp in zip(
254                    blk.get('_iso_coordinate_label'),
255                    blk.get('_iso_coordinate_formula'),
256                    ):
257                    if '_' in lbl:
258                        albl = lbl[:lbl.rfind('_')]
259                        vlbl = lbl[lbl.rfind('_')+1:]
260                    else:
261                        self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
262                        error = True
263                        continue
264                    if vlbl not in 'xyz' or len(vlbl) != 1:
265                        self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
266                        error = True
267                        continue
268                    i = 'xyz'.index(vlbl)
269                    if not ParentCoordinates.get(albl):
270                        ParentCoordinates[albl] = [None,None,None]
271                    if '+' in exp:
272                        val = exp.split('+')[0].strip()
273                    else:
274                        self.warnings += ' ERROR: _iso_coordinate_formula not parsed: '+lbl
275                        error = True
276                        continue
277                    val = G2p3.FormulaEval(val)
278                    if val is None:
279                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
280                        error = True
281                        continue
282                    ParentCoordinates[albl][i] = val
283                if error:
284                    raise Exception,"Error decoding variable labels"
285                # get mapping of modes to atomic coordinate displacements
286                displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
287                for row,col,val in zip(
288                    blk.get('_iso_displacivemodematrix_row'),
289                    blk.get('_iso_displacivemodematrix_col'),
290                    blk.get('_iso_displacivemodematrix_value'),):
291                    displacivemodematrix[int(row)-1,int(col)-1] = float(val)
292                # Invert to get mapping of atom displacements to modes
293                displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
294                # create the constraints
295                self.Constraints = []
296                for i,row in enumerate(displacivemodeInvmatrix):
297                    constraint = []
298                    for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
299                        if k == 0: continue
300                        constraint.append((k,G2varObj[j]))
301                    constraint += [shortmodelist[i],False,'f']
302                    self.Constraints.append(constraint)
303                # save the ISODISTORT info for "mode analysis"
304                self.Phase['ISODISTORT'] = {
305                    'IsoModeList' : modelist,
306                    'G2ModeList' : shortmodelist,
307                    'IsoVarList' : coordVarLbl,
308                    'G2VarList' : G2varObj,
309                    'ParentStructure' : ParentCoordinates,
310                    'Var2ModeMatrix' : displacivemodeInvmatrix,
311                    'Mode2VarMatrix' : displacivemodematrix,
312                    }
313                # make explaination dictionary
314                explaination = {}
315                for mode,shortmode in zip(modelist,shortmodelist):
316                    explaination[shortmode] = "IsoDisplace full name "+str(mode)
317                self.Constraints.append(explaination)
318
319                # # debug: show the mode var to mode relations
320                # for i,row in enumerate(displacivemodeInvmatrix):
321                #     l = ''
322                #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
323                #         if k == 0: continue
324                #         if l: l += ' + '
325                #         #l += lbl+' * '+str(k)
326                #         l += G2varLbl[j]+' * '+str(k)
327                #     print str(i) + ': '+shortmodelist[i]+' = '+l
328                # print 70*'='
329
330                # # debug: Get the ISODISPLACE offset values
331                # coordVarDelta = {}
332                # for lbl,val in zip(
333                #     blk.get('_iso_deltacoordinate_label'),
334                #     blk.get('_iso_deltacoordinate_value'),):
335                #     coordVarDelta[lbl] = float(val)
336                # modeVarDelta = {}
337                # for lbl,val in zip(
338                #     blk.get('_iso_displacivemode_label'),
339                #     blk.get('_iso_displacivemode_value'),):
340                #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
341
342                # print 70*'='
343                # # compute the mode values from the reported coordinate deltas
344                # for i,row in enumerate(displacivemodeInvmatrix):
345                #     l = ''
346                #     sl = ''
347                #     s = 0.
348                #     for lbl,k in zip(coordVarLbl,row):
349                #         if k == 0: continue
350                #         if l: l += ' + '
351                #         l += lbl+' * '+str(k)
352                #         if sl: sl += ' + '
353                #         sl += str(coordVarDelta[lbl])+' * '+str(k)
354                #         s += coordVarDelta[lbl] * k
355                #     print 'a'+str(i)+' = '+l
356                #     print '\t= '+sl
357                #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
358                #     print
359
360                # print 70*'='
361                # # compute the coordinate displacements from the reported mode values
362                # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
363                #     l = ''
364                #     sl = ''
365                #     s = 0.0
366                #     for j,k in enumerate(row):
367                #         if k == 0: continue
368                #         if l: l += ' + '
369                #         l += 'a'+str(j+1)+' * '+str(k)
370                #         if sl: sl += ' + '
371                #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
372                #         s += modeVarDelta[modelist[j]] * k
373                #     print lbl+' = '+l
374                #     print '\t= '+sl
375                #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
376                #     print
377
378                # determine the coordinate delta values from deviations from the parent structure
379                # for atmline in self.Phase['Atoms']:
380                #     lbl = atmline[0]
381                #     x,y,z = atmline[3:6]
382                #     if lbl not in ParentCoordinates:
383                #         print lbl,x,y,z
384                #         continue
385                #     px,py,pz = ParentCoordinates[lbl]
386                #     print lbl,x,y,z,x-px,y-py,z-pz
387
388                returnstat = True
389        except Exception as detail:
390            print 'CIF error:',detail # for testing
391            print sys.exc_info()[0] # for testing
392            import traceback
393            print traceback.format_exc()
394            returnstat = False
395        finally:
396            self.DoneBusy()
397        return returnstat
Note: See TracBrowser for help on using the repository browser.