source: trunk/imports/G2phase_CIF.py @ 3169

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

more complete fixes to G2phase_CIF.py importer for incommensurate structures
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 33.0 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-11-30 22:01:26 +0000 (Thu, 30 Nov 2017) $
4# $Author: vondreele $
5# $Revision: 3169 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3169 2017-11-30 22:01:26Z 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: 3169 $")
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:
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                    self.warnings += '\nThe space group has been set to "P 1". '
168                    self.warnings += "Change this in phase's General tab."
169                    self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
170                    self.warnings += G2spc.SGErrors(E)
171                SGData = G2obj.P1SGData # P 1
172            self.Phase['General']['SGData'] = SGData
173            # cell parameters
174            cell = []
175            for lbl in cellitems:
176                cell.append(cif.get_number_with_esd(blk[lbl])[0])
177            Volume = G2lat.calc_V(G2lat.cell2A(cell))
178            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
179            if Super:
180                print(blk.get('_cell_modulation_dimension',''))
181                if int(blk.get('_cell_modulation_dimension','')) > 1:
182                    msg = 'more than 3+1 super symmetry is not allowed in GSAS-II'
183                    self.errors = msg
184                    self.warnings += msg
185                    return False
186                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
187                waveDict = dict(waveloop.items())
188                SuperVec = [[float(waveDict['_cell_wave_vector_x'][0].replace('?','0')),
189                    float(waveDict['_cell_wave_vector_y'][0].replace('?','0')),
190                    float(waveDict['_cell_wave_vector_z'][0].replace('?','0'))],False,4]
191            # read in atoms
192            self.errors = 'Error during reading of atoms'
193            atomlbllist = [] # table to look up atom IDs
194            atomloop = blk.GetLoop('_atom_site_label')
195            atomkeys = [i.lower() for i in atomloop.keys()]
196            if not blk.get('_atom_site_type_symbol'):
197                self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
198            if magnetic:
199                magatomloop = blk.GetLoop('_atom_site_moment_label')
200                magatomkeys = [i.lower() for i in magatomloop.keys()]
201            if blk.get('_atom_site_aniso_label'):
202                anisoloop = blk.GetLoop('_atom_site_aniso_label')
203                anisokeys = [i.lower() for i in anisoloop.keys()]
204                anisolabels = blk.get('_atom_site_aniso_label')
205            else:
206                anisoloop = None
207                anisokeys = []
208                anisolabels = []
209            if Super:
210                occFloop = None
211                occCloop = None
212                occFdict = {}
213                occCdict = {}
214                displSloop = None
215                displFloop = None
216                dispSdict = {}
217                dispFdict = {}
218                UijFloop = None
219                UijFdict = {}
220                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
221                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
222                    occFdict = dict(occFloop.items())
223                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
224                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
225                    occCdict = dict(occCloop.items())
226                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
227                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
228                    displFdict = dict(displFloop.items())                           
229                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
230                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
231                    displSdict = dict(displSloop.items())
232                if blk.get('_atom_site_U_Fourier_atom_site_label'):
233                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
234                    UijFdict = dict(UijFloop.items())
235            self.Phase['Atoms'] = []
236            G2AtomDict = {  '_atom_site_type_symbol' : 1,
237                            '_atom_site_label' : 0,
238                            '_atom_site_fract_x' : 3,
239                            '_atom_site_fract_y' : 4,
240                            '_atom_site_fract_z' : 5,
241                            '_atom_site_occupancy' : 6,
242                            '_atom_site_aniso_u_11' : 11,
243                            '_atom_site_aniso_u_22' : 12,
244                            '_atom_site_aniso_u_33' : 13,
245                            '_atom_site_aniso_u_12' : 14,
246                            '_atom_site_aniso_u_13' : 15,
247                            '_atom_site_aniso_u_23' : 16, }
248            G2MagDict = {'_atom_site_moment_label': 0,
249                         '_atom_site_moment_crystalaxis_x':7,
250                         '_atom_site_moment_crystalaxis_y':8,
251                         '_atom_site_moment_crystalaxis_z':9}
252
253            ranIdlookup = {}
254            for aitem in atomloop:
255                if magnetic:
256                    atomlist = ['','','',0.,0.,0.,1.0,0.,0.,0.,'',0.,'I',0.01,0.,0.,0.,0.,0.,0.]
257                else:
258                    atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01,0.,0.,0.,0.,0.,0.]
259                for val,key in zip(aitem,atomkeys):
260                    col = G2AtomDict.get(key,-1)
261                    if col >= 3:
262                        atomlist[col] = cif.get_number_with_esd(val)[0]
263                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
264                    elif col is not None:
265                        atomlist[col] = val
266                    elif key in ('_atom_site_thermal_displace_type',
267                               '_atom_site_adp_type'):   #Iso or Aniso?
268                        if val.lower() == 'uani':
269                            if magnetic:
270                                atomlist[12] = 'A'
271                            else:
272                                atomlist[9] = 'A'
273                    elif key == '_atom_site_u_iso_or_equiv':
274                        uisoval = cif.get_number_with_esd(val)[0]
275                        if uisoval is not None: 
276                            if magnetic:
277                                atomlist[13] = uisoval                               
278                            else:   
279                                atomlist[10] = uisoval
280                if not atomlist[1] and atomlist[0]:
281                    typ = atomlist[0].rstrip('0123456789-+')
282                    if G2elem.CheckElement(typ):
283                        atomlist[1] = typ
284                    if not atomlist[1]:
285                        atomlist[1] = 'Xe'
286                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
287                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
288                    if magnetic:
289                        atomlist[12] = 'A'
290                    else:
291                        atomlist[9] = 'A'
292                    for val,key in zip( # load the values
293                            anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]), 
294                            anisokeys):
295                        col = G2AtomDict.get(key)
296                        if magnetic:
297                            col += 3
298                        if col:
299                            atomlist[col] = cif.get_number_with_esd(val)[0]
300                if magnetic:
301                    for mitem in magatomloop:
302                        matom = mitem[G2MagDict.get('_atom_site_moment_label',-1)]
303                        if atomlist[0] == matom:
304                            for mval,mkey in zip(mitem,magatomkeys):
305                                mcol = G2MagDict.get(mkey,-1)
306                                if mcol:
307                                    atomlist[mcol] = cif.get_number_with_esd(mval)[0]
308                            break                           
309                    atomlist[10],atomlist[11] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
310                else:
311                    atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
312                atomlist[1] = G2elem.FixValence(atomlist[1])
313                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
314                self.Phase['Atoms'].append(atomlist)
315                ranIdlookup[atomlist[0]] = atomlist[-1]
316                if atomlist[0] in atomlbllist:
317                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
318                else:
319                    atomlbllist.append(atomlist[0])
320                if Super:
321                    Sfrac = []
322                    Sadp = []                     
323                    Spos = np.zeros((4,6))
324                    nim = -1
325                    for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
326                        if item == atomlist[0]:
327                            waveType = 'Fourier'                               
328                            ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
329                            im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
330                            if im != nim:
331                                nim = im
332                            val = displFdict['_atom_site_displace_fourier_param_sin'][i]
333                            Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
334                            val = displFdict['_atom_site_displace_fourier_param_cos'][i]
335                            Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
336                    if nim >= 0:
337                        Spos = [[spos,False] for spos in Spos[:nim]]
338                    else:
339                        Spos = []
340                    if UijFdict:
341                        nim = -1
342                        Sadp = np.zeros((4,12))
343                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
344                            if item == atomlist[0]:
345                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
346                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
347                                if im != nim:
348                                    nim = im
349                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
350                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
351                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
352                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
353                        if nim >= 0:
354                            Sadp = [[sadp,False] for sadp in Sadp[:nim]]
355                        else:
356                            Sadp = []
357                   
358                    SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':[]}}
359                    atomlist.append(SSdict)
360            if len(atomlbllist) != len(self.Phase['Atoms']):
361                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
362            for lbl in phasenamefields: # get a name for the phase
363                name = blk.get(lbl)
364                if name is None:
365                    continue
366                name = name.strip()
367                if name == '?' or name == '.':
368                    continue
369                else:
370                    break
371            else: # no name found, use block name for lack of a better choice
372                name = blknm
373            self.Phase['General']['Name'] = name.strip()[:20]
374            self.Phase['General']['Super'] = Super
375            if magnetic:
376                self.Phase['General']['Type'] = 'magnetic'               
377            if Super:
378                self.Phase['General']['Type'] = 'modulated'
379                self.Phase['General']['SuperVec'] = SuperVec
380                self.Phase['General']['SuperSg'] = SuperSg
381                self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
382            if not self.isodistort_warnings:
383                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
384                    self.errors = "Error while processing ISODISTORT constraints"
385                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
386            else:
387                self.warnings += self.isodistort_warnings
388            returnstat = True
389        return returnstat
390
391    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
392        'Process ISODISTORT items to create constraints etc.'
393        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
394        'Maps ISODISTORT parm names to GSAS-II names'
395        #----------------------------------------------------------------------
396        # read in the ISODISTORT displacement modes
397        #----------------------------------------------------------------------
398        self.Constraints = []
399        explaination = {}
400        if blk.get('_iso_displacivemode_label'):
401            modelist = []
402            shortmodelist = []
403            for lbl in blk.get('_iso_displacivemode_label'):
404                modelist.append(lbl)
405                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
406                # where SSSSS is the parent spacegroup, [x,y,z] is a location
407                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
408                # this extracts the AAAAA and BBBBB parts of the string
409                if regexp:
410                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
411                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
412            # read in the coordinate offset variables names and map them to G2 names/objects
413            coordVarLbl = []
414            G2varLbl = []
415            G2varObj = []
416            error = False
417            for lbl in blk.get('_iso_deltacoordinate_label'):
418                coordVarLbl.append(lbl)
419                if '_' in lbl:
420                    albl = lbl[:lbl.rfind('_')]
421                    vlbl = lbl[lbl.rfind('_')+1:]
422                else:
423                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
424                    error = True
425                    continue
426                if albl not in atomlbllist:
427                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
428                    error = True
429                    continue
430                else:
431                    anum = atomlbllist.index(albl)
432                var = varLookup.get(vlbl)
433                if not var:
434                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
435                    error = True
436                    continue
437                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
438                G2varObj.append(G2obj.G2VarObj(
439                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
440                    ))
441            if error:
442                raise Exception("Error decoding variable labels")
443
444            if len(G2varObj) != len(modelist):
445                print ("non-square input")
446                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
447
448            error = False
449            ParentCoordinates = {}
450            for lbl,exp in zip(
451                blk.get('_iso_coordinate_label'),
452                blk.get('_iso_coordinate_formula'),
453                ):
454                if '_' in lbl:
455                    albl = lbl[:lbl.rfind('_')]
456                    vlbl = lbl[lbl.rfind('_')+1:]
457                else:
458                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
459                    error = True
460                    continue
461                if vlbl not in 'xyz' or len(vlbl) != 1:
462                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
463                    error = True
464                    continue
465                i = 'xyz'.index(vlbl)
466                if not ParentCoordinates.get(albl):
467                    ParentCoordinates[albl] = [None,None,None]
468                if '+' in exp:
469                    val = exp.split('+')[0].strip()
470                    val = G2p3.FormulaEval(val)
471                    if val is None:
472                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
473                        error = True
474                        continue
475                    ParentCoordinates[albl][i] = val
476                else:
477                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
478            if error:
479                print (self.warnings)
480                raise Exception("Error decoding variable labels")
481            # get mapping of modes to atomic coordinate displacements
482            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
483            for row,col,val in zip(
484                blk.get('_iso_displacivemodematrix_row'),
485                blk.get('_iso_displacivemodematrix_col'),
486                blk.get('_iso_displacivemodematrix_value'),):
487                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
488            # Invert to get mapping of atom displacements to modes
489            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
490            # create the constraints
491            for i,row in enumerate(displacivemodeInvmatrix):
492                constraint = []
493                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
494                    if k == 0: continue
495                    constraint.append([k,G2varObj[j]])
496                constraint += [shortmodelist[i],False,'f']
497                self.Constraints.append(constraint)
498            #----------------------------------------------------------------------
499            # save the ISODISTORT info for "mode analysis"
500            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
501            self.Phase['ISODISTORT'].update({
502                'IsoModeList' : modelist,
503                'G2ModeList' : shortmodelist,
504                'IsoVarList' : coordVarLbl,
505                'G2VarList' : G2varObj,
506                'ParentStructure' : ParentCoordinates,
507                'Var2ModeMatrix' : displacivemodeInvmatrix,
508                'Mode2VarMatrix' : displacivemodematrix,
509                })
510            # make explaination dictionary
511            for mode,shortmode in zip(modelist,shortmodelist):
512                explaination[shortmode] = "ISODISTORT full name "+str(mode)
513        #----------------------------------------------------------------------
514        # now read in the ISODISTORT occupancy modes
515        #----------------------------------------------------------------------
516        if blk.get('_iso_occupancymode_label'):
517            modelist = []
518            shortmodelist = []
519            for lbl in blk.get('_iso_occupancymode_label'):
520                modelist.append(lbl)
521                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
522                # where SSSSS is the parent spacegroup, [x,y,z] is a location
523                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
524                # this extracts the AAAAA and BBBBB parts of the string
525                if regexp:
526                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
527                lbl = lbl.replace('order','o')
528                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
529            # read in the coordinate offset variables names and map them to G2 names/objects
530            occVarLbl = []
531            G2varLbl = []
532            G2varObj = []
533            error = False
534            for lbl in blk.get('_iso_deltaoccupancy_label'):
535                occVarLbl.append(lbl)
536                if '_' in lbl:
537                    albl = lbl[:lbl.rfind('_')]
538                    vlbl = lbl[lbl.rfind('_')+1:]
539                else:
540                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
541                    error = True
542                    continue
543                if albl not in atomlbllist:
544                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
545                    error = True
546                    continue
547                else:
548                    anum = atomlbllist.index(albl)
549                var = varLookup.get(vlbl)
550                if not var:
551                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
552                    error = True
553                    continue
554                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
555                G2varObj.append(G2obj.G2VarObj(
556                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
557                    ))
558            if error:
559                raise Exception("Error decoding variable labels")
560
561            if len(G2varObj) != len(modelist):
562                print ("non-square input")
563                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
564
565            error = False
566            ParentCoordinates = {}
567            for lbl,exp in zip(
568                blk.get('_iso_occupancy_label'),
569                blk.get('_iso_occupancy_formula'),
570                ):
571                if '_' in lbl:
572                    albl = lbl[:lbl.rfind('_')]
573                    vlbl = lbl[lbl.rfind('_')+1:]
574                else:
575                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
576                    error = True
577                    continue
578                if vlbl != 'occ':
579                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
580                    error = True
581                    continue
582                if '+' in exp:
583                    val = exp.split('+')[0].strip()
584                    val = G2p3.FormulaEval(val)
585                    if val is None:
586                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
587                        error = True
588                        continue
589                    ParentCoordinates[albl] = val
590            if error:
591                raise Exception("Error decoding occupancy labels")
592            # get mapping of modes to atomic coordinate displacements
593            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
594            for row,col,val in zip(
595                blk.get('_iso_occupancymodematrix_row'),
596                blk.get('_iso_occupancymodematrix_col'),
597                blk.get('_iso_occupancymodematrix_value'),):
598                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
599            # Invert to get mapping of atom displacements to modes
600            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
601            # create the constraints
602            for i,row in enumerate(occupancymodeInvmatrix):
603                constraint = []
604                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
605                    if k == 0: continue
606                    constraint.append([k,G2varObj[j]])
607                constraint += [shortmodelist[i],False,'f']
608                self.Constraints.append(constraint)
609            #----------------------------------------------------------------------
610            # save the ISODISTORT info for "mode analysis"
611            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
612            self.Phase['ISODISTORT'].update({
613                'OccModeList' : modelist,
614                'G2OccModeList' : shortmodelist,
615                'OccVarList' : occVarLbl,
616                'G2OccVarList' : G2varObj,
617                'BaseOcc' : ParentCoordinates,
618                'Var2OccMatrix' : occupancymodeInvmatrix,
619                'Occ2VarMatrix' : occupancymodematrix,
620                })
621            # make explaination dictionary
622            for mode,shortmode in zip(modelist,shortmodelist):
623                explaination[shortmode] = "ISODISTORT full name "+str(mode)
624        #----------------------------------------------------------------------
625        # done with read
626        #----------------------------------------------------------------------
627        if explaination: self.Constraints.append(explaination)
628
629        # # debug: show the mode var to mode relations
630        # for i,row in enumerate(displacivemodeInvmatrix):
631        #     l = ''
632        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
633        #         if k == 0: continue
634        #         if l: l += ' + '
635        #         #l += lbl+' * '+str(k)
636        #         l += G2varLbl[j]+' * '+str(k)
637        #     print str(i) + ': '+shortmodelist[i]+' = '+l
638        # print 70*'='
639
640        # # debug: Get the ISODISTORT offset values
641        # coordVarDelta = {}
642        # for lbl,val in zip(
643        #     blk.get('_iso_deltacoordinate_label'),
644        #     blk.get('_iso_deltacoordinate_value'),):
645        #     coordVarDelta[lbl] = float(val)
646        # modeVarDelta = {}
647        # for lbl,val in zip(
648        #     blk.get('_iso_displacivemode_label'),
649        #     blk.get('_iso_displacivemode_value'),):
650        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
651
652        # print 70*'='
653        # # compute the mode values from the reported coordinate deltas
654        # for i,row in enumerate(displacivemodeInvmatrix):
655        #     l = ''
656        #     sl = ''
657        #     s = 0.
658        #     for lbl,k in zip(coordVarLbl,row):
659        #         if k == 0: continue
660        #         if l: l += ' + '
661        #         l += lbl+' * '+str(k)
662        #         if sl: sl += ' + '
663        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
664        #         s += coordVarDelta[lbl] * k
665        #     print 'a'+str(i)+' = '+l
666        #     print '\t= '+sl
667        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
668        #     print
669
670        # print 70*'='
671        # # compute the coordinate displacements from the reported mode values
672        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
673        #     l = ''
674        #     sl = ''
675        #     s = 0.0
676        #     for j,k in enumerate(row):
677        #         if k == 0: continue
678        #         if l: l += ' + '
679        #         l += 'a'+str(j+1)+' * '+str(k)
680        #         if sl: sl += ' + '
681        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
682        #         s += modeVarDelta[modelist[j]] * k
683        #     print lbl+' = '+l
684        #     print '\t= '+sl
685        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
686        #     print
687
688        # determine the coordinate delta values from deviations from the parent structure
689        # for atmline in self.Phase['Atoms']:
690        #     lbl = atmline[0]
691        #     x,y,z = atmline[3:6]
692        #     if lbl not in ParentCoordinates:
693        #         print lbl,x,y,z
694        #         continue
695        #     px,py,pz = ParentCoordinates[lbl]
696        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.