source: trunk/imports/G2phase_CIF.py @ 3186

Last change on this file since 3186 was 3186, checked in by vondreele, 5 years ago

revisions to integrate images so that masks are reused if not changed in integrate all & auto integrate
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 33.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-12-09 14:54:59 +0000 (Sat, 09 Dec 2017) $
4# $Author: vondreele $
5# $Revision: 3186 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3186 2017-12-09 14:54:59Z 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: 3186 $")
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            E,SGData = G2spc.SpcGroup(SpGrp)
155            if E and SpGrp:
156                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
157                if SpGrpNorm:
158                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
159            # nope, try the space group "out of the Box"
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 = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
199                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
200                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[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                mc = 0
266                if magnetic:
267                    atomlist = ['','','',0.,0.,0.,1.0, 0.,0.,0.,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
268                    mc = 3
269                else:
270                    atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
271                for val,key in zip(aitem,atomkeys):
272                    col = G2AtomDict.get(key,-1)
273                    if col >= 3:
274                        atomlist[col] = cif.get_number_with_esd(val)[0]
275                        if col >= 11: atomlist[9+mc] = 'A' # if any Aniso term is defined, set flag
276                    elif col is not None and col != -1:
277                        atomlist[col] = val
278                    elif key in ('_atom_site_thermal_displace_type',
279                               '_atom_site_adp_type'):   #Iso or Aniso?
280                        if val.lower() == 'uani':
281                            atomlist[9+mc] = 'A'
282                    elif key == '_atom_site_u_iso_or_equiv':
283                        uisoval = cif.get_number_with_esd(val)[0]
284                        if uisoval is not None: 
285                            atomlist[10+mc] = uisoval
286                if not atomlist[1] and atomlist[0]:
287                    typ = atomlist[0].rstrip('0123456789-+')
288                    if G2elem.CheckElement(typ):
289                        atomlist[1] = typ
290                    if not atomlist[1]:
291                        atomlist[1] = 'Xe'
292                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
293                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
294                    atomlist[9+mc] = 'A'
295                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
296                        col = G2AtomDict.get(key)
297                        if col:
298                            atomlist[col+mc] = cif.get_number_with_esd(val)[0]
299                if magnetic:
300                    for mitem in magatomloop:
301                        matom = mitem[G2MagDict.get('_atom_site_moment_label',-1)]
302                        if atomlist[0] == matom:
303                            for mval,mkey in zip(mitem,magatomkeys):
304                                mcol = G2MagDict.get(mkey,-1)
305                                if mcol:
306                                    atomlist[mcol] = cif.get_number_with_esd(mval)[0]
307                            break                           
308                atomlist[7+mc],atomlist[8+mc] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
309                atomlist[1] = G2elem.FixValence(atomlist[1])
310                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
311                self.Phase['Atoms'].append(atomlist)
312                ranIdlookup[atomlist[0]] = atomlist[-1]
313                if atomlist[0] in atomlbllist:
314                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
315                else:
316                    atomlbllist.append(atomlist[0])
317                if Super:
318                    Sfrac = []
319                    Sadp = []                     
320                    Spos = np.zeros((4,6))
321                    nim = -1
322                    for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
323                        if item == atomlist[0]:
324                            waveType = 'Fourier'                               
325                            ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
326                            im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
327                            if im != nim:
328                                nim = im
329                            val = displFdict['_atom_site_displace_fourier_param_sin'][i]
330                            Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
331                            val = displFdict['_atom_site_displace_fourier_param_cos'][i]
332                            Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
333                    if nim >= 0:
334                        Spos = [[spos,False] for spos in Spos[:nim]]
335                    else:
336                        Spos = []
337                    if UijFdict:
338                        nim = -1
339                        Sadp = np.zeros((4,12))
340                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
341                            if item == atomlist[0]:
342                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
343                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
344                                if im != nim:
345                                    nim = im
346                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
347                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
348                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
349                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
350                        if nim >= 0:
351                            Sadp = [[sadp,False] for sadp in Sadp[:nim]]
352                        else:
353                            Sadp = []
354                   
355                    SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':[]}}
356                    atomlist.append(SSdict)
357            if len(atomlbllist) != len(self.Phase['Atoms']):
358                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
359            for lbl in phasenamefields: # get a name for the phase
360                name = blk.get(lbl)
361                if name is None:
362                    continue
363                name = name.strip()
364                if name == '?' or name == '.':
365                    continue
366                else:
367                    break
368            else: # no name found, use block name for lack of a better choice
369                name = blknm
370            self.Phase['General']['Name'] = name.strip()[:20]
371            self.Phase['General']['Super'] = Super
372            if magnetic:
373                self.Phase['General']['Type'] = 'magnetic'               
374            if Super:
375                self.Phase['General']['Type'] = 'modulated'
376                self.Phase['General']['SuperVec'] = SuperVec
377                self.Phase['General']['SuperSg'] = SuperSg
378                self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
379            if not self.isodistort_warnings:
380                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
381                    self.errors = "Error while processing ISODISTORT constraints"
382                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
383            else:
384                self.warnings += self.isodistort_warnings
385            returnstat = True
386        return returnstat
387
388    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
389        'Process ISODISTORT items to create constraints etc.'
390        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
391        'Maps ISODISTORT parm names to GSAS-II names'
392        #----------------------------------------------------------------------
393        # read in the ISODISTORT displacement modes
394        #----------------------------------------------------------------------
395        self.Constraints = []
396        explaination = {}
397        if blk.get('_iso_displacivemode_label'):
398            modelist = []
399            shortmodelist = []
400            for lbl in blk.get('_iso_displacivemode_label'):
401                modelist.append(lbl)
402                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
403                # where SSSSS is the parent spacegroup, [x,y,z] is a location
404                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
405                # this extracts the AAAAA and BBBBB parts of the string
406                if regexp:
407                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
408                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
409            # read in the coordinate offset variables names and map them to G2 names/objects
410            coordVarLbl = []
411            G2varLbl = []
412            G2varObj = []
413            error = False
414            for lbl in blk.get('_iso_deltacoordinate_label'):
415                coordVarLbl.append(lbl)
416                if '_' in lbl:
417                    albl = lbl[:lbl.rfind('_')]
418                    vlbl = lbl[lbl.rfind('_')+1:]
419                else:
420                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
421                    error = True
422                    continue
423                if albl not in atomlbllist:
424                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
425                    error = True
426                    continue
427                else:
428                    anum = atomlbllist.index(albl)
429                var = varLookup.get(vlbl)
430                if not var:
431                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
432                    error = True
433                    continue
434                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
435                G2varObj.append(G2obj.G2VarObj(
436                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
437                    ))
438            if error:
439                raise Exception("Error decoding variable labels")
440
441            if len(G2varObj) != len(modelist):
442                print ("non-square input")
443                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
444
445            error = False
446            ParentCoordinates = {}
447            for lbl,exp in zip(
448                blk.get('_iso_coordinate_label'),
449                blk.get('_iso_coordinate_formula'),
450                ):
451                if '_' in lbl:
452                    albl = lbl[:lbl.rfind('_')]
453                    vlbl = lbl[lbl.rfind('_')+1:]
454                else:
455                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
456                    error = True
457                    continue
458                if vlbl not in 'xyz' or len(vlbl) != 1:
459                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
460                    error = True
461                    continue
462                i = 'xyz'.index(vlbl)
463                if not ParentCoordinates.get(albl):
464                    ParentCoordinates[albl] = [None,None,None]
465                if '+' in exp:
466                    val = exp.split('+')[0].strip()
467                    val = G2p3.FormulaEval(val)
468                    if val is None:
469                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
470                        error = True
471                        continue
472                    ParentCoordinates[albl][i] = val
473                else:
474                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
475            if error:
476                print (self.warnings)
477                raise Exception("Error decoding variable labels")
478            # get mapping of modes to atomic coordinate displacements
479            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
480            for row,col,val in zip(
481                blk.get('_iso_displacivemodematrix_row'),
482                blk.get('_iso_displacivemodematrix_col'),
483                blk.get('_iso_displacivemodematrix_value'),):
484                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
485            # Invert to get mapping of atom displacements to modes
486            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
487            # create the constraints
488            for i,row in enumerate(displacivemodeInvmatrix):
489                constraint = []
490                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
491                    if k == 0: continue
492                    constraint.append([k,G2varObj[j]])
493                constraint += [shortmodelist[i],False,'f']
494                self.Constraints.append(constraint)
495            #----------------------------------------------------------------------
496            # save the ISODISTORT info for "mode analysis"
497            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
498            self.Phase['ISODISTORT'].update({
499                'IsoModeList' : modelist,
500                'G2ModeList' : shortmodelist,
501                'IsoVarList' : coordVarLbl,
502                'G2VarList' : G2varObj,
503                'ParentStructure' : ParentCoordinates,
504                'Var2ModeMatrix' : displacivemodeInvmatrix,
505                'Mode2VarMatrix' : displacivemodematrix,
506                })
507            # make explaination dictionary
508            for mode,shortmode in zip(modelist,shortmodelist):
509                explaination[shortmode] = "ISODISTORT full name "+str(mode)
510        #----------------------------------------------------------------------
511        # now read in the ISODISTORT occupancy modes
512        #----------------------------------------------------------------------
513        if blk.get('_iso_occupancymode_label'):
514            modelist = []
515            shortmodelist = []
516            for lbl in blk.get('_iso_occupancymode_label'):
517                modelist.append(lbl)
518                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
519                # where SSSSS is the parent spacegroup, [x,y,z] is a location
520                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
521                # this extracts the AAAAA and BBBBB parts of the string
522                if regexp:
523                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
524                lbl = lbl.replace('order','o')
525                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
526            # read in the coordinate offset variables names and map them to G2 names/objects
527            occVarLbl = []
528            G2varLbl = []
529            G2varObj = []
530            error = False
531            for lbl in blk.get('_iso_deltaoccupancy_label'):
532                occVarLbl.append(lbl)
533                if '_' in lbl:
534                    albl = lbl[:lbl.rfind('_')]
535                    vlbl = lbl[lbl.rfind('_')+1:]
536                else:
537                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
538                    error = True
539                    continue
540                if albl not in atomlbllist:
541                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
542                    error = True
543                    continue
544                else:
545                    anum = atomlbllist.index(albl)
546                var = varLookup.get(vlbl)
547                if not var:
548                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
549                    error = True
550                    continue
551                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
552                G2varObj.append(G2obj.G2VarObj(
553                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
554                    ))
555            if error:
556                raise Exception("Error decoding variable labels")
557
558            if len(G2varObj) != len(modelist):
559                print ("non-square input")
560                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
561
562            error = False
563            ParentCoordinates = {}
564            for lbl,exp in zip(
565                blk.get('_iso_occupancy_label'),
566                blk.get('_iso_occupancy_formula'),
567                ):
568                if '_' in lbl:
569                    albl = lbl[:lbl.rfind('_')]
570                    vlbl = lbl[lbl.rfind('_')+1:]
571                else:
572                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
573                    error = True
574                    continue
575                if vlbl != 'occ':
576                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
577                    error = True
578                    continue
579                if '+' in exp:
580                    val = exp.split('+')[0].strip()
581                    val = G2p3.FormulaEval(val)
582                    if val is None:
583                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
584                        error = True
585                        continue
586                    ParentCoordinates[albl] = val
587            if error:
588                raise Exception("Error decoding occupancy labels")
589            # get mapping of modes to atomic coordinate displacements
590            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
591            for row,col,val in zip(
592                blk.get('_iso_occupancymodematrix_row'),
593                blk.get('_iso_occupancymodematrix_col'),
594                blk.get('_iso_occupancymodematrix_value'),):
595                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
596            # Invert to get mapping of atom displacements to modes
597            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
598            # create the constraints
599            for i,row in enumerate(occupancymodeInvmatrix):
600                constraint = []
601                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
602                    if k == 0: continue
603                    constraint.append([k,G2varObj[j]])
604                constraint += [shortmodelist[i],False,'f']
605                self.Constraints.append(constraint)
606            #----------------------------------------------------------------------
607            # save the ISODISTORT info for "mode analysis"
608            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
609            self.Phase['ISODISTORT'].update({
610                'OccModeList' : modelist,
611                'G2OccModeList' : shortmodelist,
612                'OccVarList' : occVarLbl,
613                'G2OccVarList' : G2varObj,
614                'BaseOcc' : ParentCoordinates,
615                'Var2OccMatrix' : occupancymodeInvmatrix,
616                'Occ2VarMatrix' : occupancymodematrix,
617                })
618            # make explaination dictionary
619            for mode,shortmode in zip(modelist,shortmodelist):
620                explaination[shortmode] = "ISODISTORT full name "+str(mode)
621        #----------------------------------------------------------------------
622        # done with read
623        #----------------------------------------------------------------------
624        if explaination: self.Constraints.append(explaination)
625
626        # # debug: show the mode var to mode relations
627        # for i,row in enumerate(displacivemodeInvmatrix):
628        #     l = ''
629        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
630        #         if k == 0: continue
631        #         if l: l += ' + '
632        #         #l += lbl+' * '+str(k)
633        #         l += G2varLbl[j]+' * '+str(k)
634        #     print str(i) + ': '+shortmodelist[i]+' = '+l
635        # print 70*'='
636
637        # # debug: Get the ISODISTORT offset values
638        # coordVarDelta = {}
639        # for lbl,val in zip(
640        #     blk.get('_iso_deltacoordinate_label'),
641        #     blk.get('_iso_deltacoordinate_value'),):
642        #     coordVarDelta[lbl] = float(val)
643        # modeVarDelta = {}
644        # for lbl,val in zip(
645        #     blk.get('_iso_displacivemode_label'),
646        #     blk.get('_iso_displacivemode_value'),):
647        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
648
649        # print 70*'='
650        # # compute the mode values from the reported coordinate deltas
651        # for i,row in enumerate(displacivemodeInvmatrix):
652        #     l = ''
653        #     sl = ''
654        #     s = 0.
655        #     for lbl,k in zip(coordVarLbl,row):
656        #         if k == 0: continue
657        #         if l: l += ' + '
658        #         l += lbl+' * '+str(k)
659        #         if sl: sl += ' + '
660        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
661        #         s += coordVarDelta[lbl] * k
662        #     print 'a'+str(i)+' = '+l
663        #     print '\t= '+sl
664        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
665        #     print
666
667        # print 70*'='
668        # # compute the coordinate displacements from the reported mode values
669        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
670        #     l = ''
671        #     sl = ''
672        #     s = 0.0
673        #     for j,k in enumerate(row):
674        #         if k == 0: continue
675        #         if l: l += ' + '
676        #         l += 'a'+str(j+1)+' * '+str(k)
677        #         if sl: sl += ' + '
678        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
679        #         s += modeVarDelta[modelist[j]] * k
680        #     print lbl+' = '+l
681        #     print '\t= '+sl
682        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
683        #     print
684
685        # determine the coordinate delta values from deviations from the parent structure
686        # for atmline in self.Phase['Atoms']:
687        #     lbl = atmline[0]
688        #     x,y,z = atmline[3:6]
689        #     if lbl not in ParentCoordinates:
690        #         print lbl,x,y,z
691        #         continue
692        #     px,py,pz = ParentCoordinates[lbl]
693        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.