source: trunk/imports/G2phase_CIF.py @ 3172

Last change on this file since 3172 was 3172, checked in by vondreele, 6 years ago

change from wx.DefaultSize? to our defaults if window size is 'None'
A TODO in Integrate for saving the x,y --> 2th,azm map between images
complete fix to allowed super symmetries by lattice + pt. grp lookup
also fix operator check for complete super symmetry
finish cif import of super symmetry cases
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 33.6 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-12-04 22:22:31 +0000 (Mon, 04 Dec 2017) $
4# $Author: vondreele $
5# $Revision: 3172 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3172 2017-12-04 22:22:31Z vondreele $
8########### SVN repository information ###################
9'''
10*Module G2phase_CIF: Coordinates from CIF*
11------------------------------------------
12
13Parses a CIF using  PyCifRW from James Hester and pulls out the
14structural information.
15
16If a CIF generated by ISODISTORT is encountered, extra information is
17added to the phase entry and constraints are generated.
18
19'''
20# Routines to import Phase information from CIF files
21from __future__ import division, print_function
22import sys
23import random as ran
24import numpy as np
25import re
26import GSASIIIO as G2IO
27import GSASIIobj as G2obj
28import GSASIIspc as G2spc
29import GSASIIElem as G2elem
30import GSASIIlattice as G2lat
31import GSASIIpy3 as G2p3
32import GSASIIpath
33GSASIIpath.SetVersionNumber("$Revision: 3172 $")
34import CifFile as cif # PyCifRW from James Hester
35
36class CIFPhaseReader(G2obj.ImportPhase):
37    'Implements a phase importer from a possibly multi-block CIF file'
38    def __init__(self):
39        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
40            extensionlist=('.CIF','.cif','.txt','.mcif'),
41            strictExtension=False,
42            formatName = 'CIF',
43            longFormatName = 'Crystallographic Information File import'
44            )
45       
46    def ContentsValidator(self, filename):
47        fp = open(filename,'r')
48        return self.CIFValidator(fp)
49        fp.close()
50
51    def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused):
52        self.isodistort_warnings = ''
53        self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
54        # make sure the ranId is really unique!
55        while self.Phase['ranId'] in usedRanIdList:
56            self.Phase['ranId'] = ran.randint(0,sys.maxsize)
57        returnstat = False
58        cellitems = (
59            '_cell_length_a','_cell_length_b','_cell_length_c',
60            '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',)
61#        cellwaveitems = (
62#            '_cell_wave_vector_seq_id',
63#            '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z')
64        reqitems = (
65             '_atom_site_fract_x',
66             '_atom_site_fract_y',
67             '_atom_site_fract_z',
68            )
69        phasenamefields = (
70            '_chemical_name_common',
71            '_pd_phase_name',
72            '_chemical_formula_sum'
73            )
74        try:
75            cf = G2obj.ReadCIF(filename)
76        except:
77            msg = 'Unreadable cif file'
78            self.errors = msg
79            self.warnings += msg
80            return False
81        # scan blocks for structural info
82        self.errors = 'Error during scan of blocks for datasets'
83        str_blklist = []
84        for blk in cf.keys():
85            for r in reqitems+cellitems:
86                if r not in cf[blk].keys():
87                    break
88            else:
89                str_blklist.append(blk)
90        if not str_blklist:
91            selblk = None # no block to choose
92        elif len(str_blklist) == 1: # only one choice
93            selblk = 0
94        else:                       # choose from options
95            choice = []
96            for blknm in str_blklist:
97                choice.append('')
98                # accumumlate some info about this phase
99                choice[-1] += blknm + ': '
100                for i in phasenamefields: # get a name for the phase
101                    name = cf[blknm].get(i,'phase name').strip()
102                    if name is None or name == '?' or name == '.':
103                        continue
104                    else:
105                        choice[-1] += name.strip()[:20] + ', '
106                        break
107                na = len(cf[blknm].get("_atom_site_fract_x"))
108                if na == 1:
109                    choice[-1] += '1 atom'
110                else:
111                    choice[-1] += ('%d' % nd) + ' atoms'
112                choice[-1] += ', cell: '
113                fmt = "%.2f,"
114                for i,key in enumerate(cellitems):
115                    if i == 3: fmt = "%.f,"
116                    if i == 5: fmt = "%.f"
117                    choice[-1] += fmt % cif.get_number_with_esd(
118                        cf[blknm].get(key))[0]
119                sg = cf[blknm].get("_symmetry_space_group_name_H-M",'')
120                if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'')
121                sg = sg.replace('_','')
122                if sg: choice[-1] += ', (' + sg.strip() + ')'
123            selblk = G2IO.PhaseSelector(choice,ParentFrame=ParentFrame,
124                title= 'Select a phase from one the CIF data_ blocks below',size=(600,100))
125        self.errors = 'Error during reading of selected block'
126        if selblk is None:
127            returnstat = False # no block selected or available
128        else:
129            blknm = str_blklist[selblk]
130            blk = cf[str_blklist[selblk]]
131            E = True
132            Super = False
133            if blk.get("_space_group_ssg_name",''):
134                Super = True
135            magnetic = False
136            self.Phase['General']['Type'] = 'nuclear'
137            SpGrp = blk.get("_symmetry_space_group_name_H-M",'')
138            if not SpGrp:
139                SpGrp = blk.get("_space_group_name_H-M_alt",'')
140            if not SpGrp:
141                SpGrp = blk.get("_parent_space_group.name_H-M",'')
142                if SpGrp:   #TODO need to decide if read nuclear phase or magnetic phase
143                    magnetic = True
144                    self.Phase['General']['Type'] = 'magnetic'
145                    self.Phase['General']['AtomPtrs'] = [3,1,10,12]
146            if Super:
147                sspgrp = blk.get("_space_group_ssg_name",'').split('(')
148                SpGrp = sspgrp[0]
149                SuperSg = '('+sspgrp[1].replace('\\','')
150                Super = True
151                SuperVec = [[0,0,.1],False,4]
152            SpGrp = SpGrp.replace('_','')
153            # try normalizing the space group, to see if we can pick the space group out of a table
154            SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
155            if SpGrpNorm:
156                E,SGData = G2spc.SpcGroup(SpGrpNorm)
157            # nope, try the space group "out of the Box"
158            if E and SpGrp:
159                E,SGData = G2spc.SpcGroup(SpGrp)
160            if E:
161                if not SpGrp:
162                    self.warnings += 'No space group name was found in the CIF.'
163                    self.warnings += '\nThe space group has been set to "P 1". '
164                    self.warnings += "Change this in phase's General tab."
165                else:
166                    self.warnings += 'ERROR in space group symbol '+SpGrp
167                    if 'X' in SpGrp:
168                        self.warnings += '\nAd hoc incommensurate space groups not allowed in GSAS-II'
169                    self.warnings += '\nThe space group has been set to "P 1". '
170                    self.warnings += "Change this in phase's General tab."
171                    self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
172                    self.warnings += G2spc.SGErrors(E)
173                SGData = G2obj.P1SGData # P 1
174            self.Phase['General']['SGData'] = SGData
175            if Super:
176                E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
177                if E:
178                    self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
179                    self.warnings += '\n'+E
180                    SuperSg = SuperSg[:SuperSg.index(')')+1]
181                    self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
182                    E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
183                self.Phase['General']['SSGData'] = SSGData
184            # cell parameters
185            cell = []
186            for lbl in cellitems:
187                cell.append(cif.get_number_with_esd(blk[lbl])[0])
188            Volume = G2lat.calc_V(G2lat.cell2A(cell))
189            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
190            if Super:
191                if int(blk.get('_cell_modulation_dimension','')) > 1:
192                    msg = 'more than 3+1 super symmetry is not allowed in GSAS-II'
193                    self.errors = msg
194                    self.warnings += msg
195                    return False
196                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
197                waveDict = dict(waveloop.items())
198                SuperVec = [[float(waveDict['_cell_wave_vector_x'][0].replace('?','0')),
199                    float(waveDict['_cell_wave_vector_y'][0].replace('?','0')),
200                    float(waveDict['_cell_wave_vector_z'][0].replace('?','0'))],False,4]
201            # read in atoms
202            self.errors = 'Error during reading of atoms'
203            atomlbllist = [] # table to look up atom IDs
204            atomloop = blk.GetLoop('_atom_site_label')
205            atomkeys = [i.lower() for i in atomloop.keys()]
206            if not blk.get('_atom_site_type_symbol'):
207                self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
208            if magnetic:
209                magatomloop = blk.GetLoop('_atom_site_moment_label')
210                magatomkeys = [i.lower() for i in magatomloop.keys()]
211            if blk.get('_atom_site_aniso_label'):
212                anisoloop = blk.GetLoop('_atom_site_aniso_label')
213                anisokeys = [i.lower() for i in anisoloop.keys()]
214                anisolabels = blk.get('_atom_site_aniso_label')
215            else:
216                anisoloop = None
217                anisokeys = []
218                anisolabels = []
219            if Super:
220                occFloop = None
221                occCloop = None
222                occFdict = {}
223                occCdict = {}
224                displSloop = None
225                displFloop = None
226                dispSdict = {}
227                dispFdict = {}
228                UijFloop = None
229                UijFdict = {}
230                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
231                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
232                    occFdict = dict(occFloop.items())
233                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
234                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
235                    occCdict = dict(occCloop.items())
236                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
237                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
238                    displFdict = dict(displFloop.items())                           
239                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
240                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
241                    displSdict = dict(displSloop.items())
242                if blk.get('_atom_site_U_Fourier_atom_site_label'):
243                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
244                    UijFdict = dict(UijFloop.items())
245            self.Phase['Atoms'] = []
246            G2AtomDict = {  '_atom_site_type_symbol' : 1,
247                            '_atom_site_label' : 0,
248                            '_atom_site_fract_x' : 3,
249                            '_atom_site_fract_y' : 4,
250                            '_atom_site_fract_z' : 5,
251                            '_atom_site_occupancy' : 6,
252                            '_atom_site_aniso_u_11' : 11,
253                            '_atom_site_aniso_u_22' : 12,
254                            '_atom_site_aniso_u_33' : 13,
255                            '_atom_site_aniso_u_12' : 14,
256                            '_atom_site_aniso_u_13' : 15,
257                            '_atom_site_aniso_u_23' : 16, }
258            G2MagDict = {'_atom_site_moment_label': 0,
259                         '_atom_site_moment_crystalaxis_x':7,
260                         '_atom_site_moment_crystalaxis_y':8,
261                         '_atom_site_moment_crystalaxis_z':9}
262
263            ranIdlookup = {}
264            for aitem in atomloop:
265                if magnetic:
266                    atomlist = ['','','',0.,0.,0.,1.0,0.,0.,0.,'',0.,'I',0.01,0.,0.,0.,0.,0.,0.]
267                else:
268                    atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01,0.,0.,0.,0.,0.,0.]
269                for val,key in zip(aitem,atomkeys):
270                    col = G2AtomDict.get(key,-1)
271                    if col >= 3:
272                        atomlist[col] = cif.get_number_with_esd(val)[0]
273                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
274                    elif col is not None:
275                        atomlist[col] = val
276                    elif key in ('_atom_site_thermal_displace_type',
277                               '_atom_site_adp_type'):   #Iso or Aniso?
278                        if val.lower() == 'uani':
279                            if magnetic:
280                                atomlist[12] = 'A'
281                            else:
282                                atomlist[9] = 'A'
283                    elif key == '_atom_site_u_iso_or_equiv':
284                        uisoval = cif.get_number_with_esd(val)[0]
285                        if uisoval is not None: 
286                            if magnetic:
287                                atomlist[13] = uisoval                               
288                            else:   
289                                atomlist[10] = uisoval
290                if not atomlist[1] and atomlist[0]:
291                    typ = atomlist[0].rstrip('0123456789-+')
292                    if G2elem.CheckElement(typ):
293                        atomlist[1] = typ
294                    if not atomlist[1]:
295                        atomlist[1] = 'Xe'
296                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
297                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
298                    if magnetic:
299                        atomlist[12] = 'A'
300                    else:
301                        atomlist[9] = 'A'
302                    for val,key in zip( # load the values
303                            anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]), 
304                            anisokeys):
305                        col = G2AtomDict.get(key)
306                        if magnetic:
307                            col += 3
308                        if col:
309                            atomlist[col] = cif.get_number_with_esd(val)[0]
310                if magnetic:
311                    for mitem in magatomloop:
312                        matom = mitem[G2MagDict.get('_atom_site_moment_label',-1)]
313                        if atomlist[0] == matom:
314                            for mval,mkey in zip(mitem,magatomkeys):
315                                mcol = G2MagDict.get(mkey,-1)
316                                if mcol:
317                                    atomlist[mcol] = cif.get_number_with_esd(mval)[0]
318                            break                           
319                    atomlist[10],atomlist[11] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
320                else:
321                    atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
322                atomlist[1] = G2elem.FixValence(atomlist[1])
323                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
324                self.Phase['Atoms'].append(atomlist)
325                ranIdlookup[atomlist[0]] = atomlist[-1]
326                if atomlist[0] in atomlbllist:
327                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
328                else:
329                    atomlbllist.append(atomlist[0])
330                if Super:
331                    Sfrac = []
332                    Sadp = []                     
333                    Spos = np.zeros((4,6))
334                    nim = -1
335                    for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
336                        if item == atomlist[0]:
337                            waveType = 'Fourier'                               
338                            ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
339                            im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
340                            if im != nim:
341                                nim = im
342                            val = displFdict['_atom_site_displace_fourier_param_sin'][i]
343                            Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
344                            val = displFdict['_atom_site_displace_fourier_param_cos'][i]
345                            Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
346                    if nim >= 0:
347                        Spos = [[spos,False] for spos in Spos[:nim]]
348                    else:
349                        Spos = []
350                    if UijFdict:
351                        nim = -1
352                        Sadp = np.zeros((4,12))
353                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
354                            if item == atomlist[0]:
355                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
356                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
357                                if im != nim:
358                                    nim = im
359                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
360                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
361                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
362                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
363                        if nim >= 0:
364                            Sadp = [[sadp,False] for sadp in Sadp[:nim]]
365                        else:
366                            Sadp = []
367                   
368                    SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':[]}}
369                    atomlist.append(SSdict)
370            if len(atomlbllist) != len(self.Phase['Atoms']):
371                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
372            for lbl in phasenamefields: # get a name for the phase
373                name = blk.get(lbl)
374                if name is None:
375                    continue
376                name = name.strip()
377                if name == '?' or name == '.':
378                    continue
379                else:
380                    break
381            else: # no name found, use block name for lack of a better choice
382                name = blknm
383            self.Phase['General']['Name'] = name.strip()[:20]
384            self.Phase['General']['Super'] = Super
385            if magnetic:
386                self.Phase['General']['Type'] = 'magnetic'               
387            if Super:
388                self.Phase['General']['Type'] = 'modulated'
389                self.Phase['General']['SuperVec'] = SuperVec
390                self.Phase['General']['SuperSg'] = SuperSg
391                self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
392            if not self.isodistort_warnings:
393                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
394                    self.errors = "Error while processing ISODISTORT constraints"
395                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
396            else:
397                self.warnings += self.isodistort_warnings
398            returnstat = True
399        return returnstat
400
401    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
402        'Process ISODISTORT items to create constraints etc.'
403        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
404        'Maps ISODISTORT parm names to GSAS-II names'
405        #----------------------------------------------------------------------
406        # read in the ISODISTORT displacement modes
407        #----------------------------------------------------------------------
408        self.Constraints = []
409        explaination = {}
410        if blk.get('_iso_displacivemode_label'):
411            modelist = []
412            shortmodelist = []
413            for lbl in blk.get('_iso_displacivemode_label'):
414                modelist.append(lbl)
415                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
416                # where SSSSS is the parent spacegroup, [x,y,z] is a location
417                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
418                # this extracts the AAAAA and BBBBB parts of the string
419                if regexp:
420                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
421                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
422            # read in the coordinate offset variables names and map them to G2 names/objects
423            coordVarLbl = []
424            G2varLbl = []
425            G2varObj = []
426            error = False
427            for lbl in blk.get('_iso_deltacoordinate_label'):
428                coordVarLbl.append(lbl)
429                if '_' in lbl:
430                    albl = lbl[:lbl.rfind('_')]
431                    vlbl = lbl[lbl.rfind('_')+1:]
432                else:
433                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
434                    error = True
435                    continue
436                if albl not in atomlbllist:
437                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
438                    error = True
439                    continue
440                else:
441                    anum = atomlbllist.index(albl)
442                var = varLookup.get(vlbl)
443                if not var:
444                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
445                    error = True
446                    continue
447                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
448                G2varObj.append(G2obj.G2VarObj(
449                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
450                    ))
451            if error:
452                raise Exception("Error decoding variable labels")
453
454            if len(G2varObj) != len(modelist):
455                print ("non-square input")
456                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
457
458            error = False
459            ParentCoordinates = {}
460            for lbl,exp in zip(
461                blk.get('_iso_coordinate_label'),
462                blk.get('_iso_coordinate_formula'),
463                ):
464                if '_' in lbl:
465                    albl = lbl[:lbl.rfind('_')]
466                    vlbl = lbl[lbl.rfind('_')+1:]
467                else:
468                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
469                    error = True
470                    continue
471                if vlbl not in 'xyz' or len(vlbl) != 1:
472                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
473                    error = True
474                    continue
475                i = 'xyz'.index(vlbl)
476                if not ParentCoordinates.get(albl):
477                    ParentCoordinates[albl] = [None,None,None]
478                if '+' in exp:
479                    val = exp.split('+')[0].strip()
480                    val = G2p3.FormulaEval(val)
481                    if val is None:
482                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
483                        error = True
484                        continue
485                    ParentCoordinates[albl][i] = val
486                else:
487                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
488            if error:
489                print (self.warnings)
490                raise Exception("Error decoding variable labels")
491            # get mapping of modes to atomic coordinate displacements
492            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
493            for row,col,val in zip(
494                blk.get('_iso_displacivemodematrix_row'),
495                blk.get('_iso_displacivemodematrix_col'),
496                blk.get('_iso_displacivemodematrix_value'),):
497                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
498            # Invert to get mapping of atom displacements to modes
499            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
500            # create the constraints
501            for i,row in enumerate(displacivemodeInvmatrix):
502                constraint = []
503                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
504                    if k == 0: continue
505                    constraint.append([k,G2varObj[j]])
506                constraint += [shortmodelist[i],False,'f']
507                self.Constraints.append(constraint)
508            #----------------------------------------------------------------------
509            # save the ISODISTORT info for "mode analysis"
510            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
511            self.Phase['ISODISTORT'].update({
512                'IsoModeList' : modelist,
513                'G2ModeList' : shortmodelist,
514                'IsoVarList' : coordVarLbl,
515                'G2VarList' : G2varObj,
516                'ParentStructure' : ParentCoordinates,
517                'Var2ModeMatrix' : displacivemodeInvmatrix,
518                'Mode2VarMatrix' : displacivemodematrix,
519                })
520            # make explaination dictionary
521            for mode,shortmode in zip(modelist,shortmodelist):
522                explaination[shortmode] = "ISODISTORT full name "+str(mode)
523        #----------------------------------------------------------------------
524        # now read in the ISODISTORT occupancy modes
525        #----------------------------------------------------------------------
526        if blk.get('_iso_occupancymode_label'):
527            modelist = []
528            shortmodelist = []
529            for lbl in blk.get('_iso_occupancymode_label'):
530                modelist.append(lbl)
531                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
532                # where SSSSS is the parent spacegroup, [x,y,z] is a location
533                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
534                # this extracts the AAAAA and BBBBB parts of the string
535                if regexp:
536                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
537                lbl = lbl.replace('order','o')
538                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
539            # read in the coordinate offset variables names and map them to G2 names/objects
540            occVarLbl = []
541            G2varLbl = []
542            G2varObj = []
543            error = False
544            for lbl in blk.get('_iso_deltaoccupancy_label'):
545                occVarLbl.append(lbl)
546                if '_' in lbl:
547                    albl = lbl[:lbl.rfind('_')]
548                    vlbl = lbl[lbl.rfind('_')+1:]
549                else:
550                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
551                    error = True
552                    continue
553                if albl not in atomlbllist:
554                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
555                    error = True
556                    continue
557                else:
558                    anum = atomlbllist.index(albl)
559                var = varLookup.get(vlbl)
560                if not var:
561                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
562                    error = True
563                    continue
564                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
565                G2varObj.append(G2obj.G2VarObj(
566                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
567                    ))
568            if error:
569                raise Exception("Error decoding variable labels")
570
571            if len(G2varObj) != len(modelist):
572                print ("non-square input")
573                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
574
575            error = False
576            ParentCoordinates = {}
577            for lbl,exp in zip(
578                blk.get('_iso_occupancy_label'),
579                blk.get('_iso_occupancy_formula'),
580                ):
581                if '_' in lbl:
582                    albl = lbl[:lbl.rfind('_')]
583                    vlbl = lbl[lbl.rfind('_')+1:]
584                else:
585                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
586                    error = True
587                    continue
588                if vlbl != 'occ':
589                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
590                    error = True
591                    continue
592                if '+' in exp:
593                    val = exp.split('+')[0].strip()
594                    val = G2p3.FormulaEval(val)
595                    if val is None:
596                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
597                        error = True
598                        continue
599                    ParentCoordinates[albl] = val
600            if error:
601                raise Exception("Error decoding occupancy labels")
602            # get mapping of modes to atomic coordinate displacements
603            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
604            for row,col,val in zip(
605                blk.get('_iso_occupancymodematrix_row'),
606                blk.get('_iso_occupancymodematrix_col'),
607                blk.get('_iso_occupancymodematrix_value'),):
608                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
609            # Invert to get mapping of atom displacements to modes
610            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
611            # create the constraints
612            for i,row in enumerate(occupancymodeInvmatrix):
613                constraint = []
614                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
615                    if k == 0: continue
616                    constraint.append([k,G2varObj[j]])
617                constraint += [shortmodelist[i],False,'f']
618                self.Constraints.append(constraint)
619            #----------------------------------------------------------------------
620            # save the ISODISTORT info for "mode analysis"
621            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
622            self.Phase['ISODISTORT'].update({
623                'OccModeList' : modelist,
624                'G2OccModeList' : shortmodelist,
625                'OccVarList' : occVarLbl,
626                'G2OccVarList' : G2varObj,
627                'BaseOcc' : ParentCoordinates,
628                'Var2OccMatrix' : occupancymodeInvmatrix,
629                'Occ2VarMatrix' : occupancymodematrix,
630                })
631            # make explaination dictionary
632            for mode,shortmode in zip(modelist,shortmodelist):
633                explaination[shortmode] = "ISODISTORT full name "+str(mode)
634        #----------------------------------------------------------------------
635        # done with read
636        #----------------------------------------------------------------------
637        if explaination: self.Constraints.append(explaination)
638
639        # # debug: show the mode var to mode relations
640        # for i,row in enumerate(displacivemodeInvmatrix):
641        #     l = ''
642        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
643        #         if k == 0: continue
644        #         if l: l += ' + '
645        #         #l += lbl+' * '+str(k)
646        #         l += G2varLbl[j]+' * '+str(k)
647        #     print str(i) + ': '+shortmodelist[i]+' = '+l
648        # print 70*'='
649
650        # # debug: Get the ISODISTORT offset values
651        # coordVarDelta = {}
652        # for lbl,val in zip(
653        #     blk.get('_iso_deltacoordinate_label'),
654        #     blk.get('_iso_deltacoordinate_value'),):
655        #     coordVarDelta[lbl] = float(val)
656        # modeVarDelta = {}
657        # for lbl,val in zip(
658        #     blk.get('_iso_displacivemode_label'),
659        #     blk.get('_iso_displacivemode_value'),):
660        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
661
662        # print 70*'='
663        # # compute the mode values from the reported coordinate deltas
664        # for i,row in enumerate(displacivemodeInvmatrix):
665        #     l = ''
666        #     sl = ''
667        #     s = 0.
668        #     for lbl,k in zip(coordVarLbl,row):
669        #         if k == 0: continue
670        #         if l: l += ' + '
671        #         l += lbl+' * '+str(k)
672        #         if sl: sl += ' + '
673        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
674        #         s += coordVarDelta[lbl] * k
675        #     print 'a'+str(i)+' = '+l
676        #     print '\t= '+sl
677        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
678        #     print
679
680        # print 70*'='
681        # # compute the coordinate displacements from the reported mode values
682        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
683        #     l = ''
684        #     sl = ''
685        #     s = 0.0
686        #     for j,k in enumerate(row):
687        #         if k == 0: continue
688        #         if l: l += ' + '
689        #         l += 'a'+str(j+1)+' * '+str(k)
690        #         if sl: sl += ' + '
691        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
692        #         s += modeVarDelta[modelist[j]] * k
693        #     print lbl+' = '+l
694        #     print '\t= '+sl
695        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
696        #     print
697
698        # determine the coordinate delta values from deviations from the parent structure
699        # for atmline in self.Phase['Atoms']:
700        #     lbl = atmline[0]
701        #     x,y,z = atmline[3:6]
702        #     if lbl not in ParentCoordinates:
703        #         print lbl,x,y,z
704        #         continue
705        #     px,py,pz = ParentCoordinates[lbl]
706        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.