source: trunk/imports/G2phase_CIF.py @ 3194

Last change on this file since 3194 was 3194, checked in by vondreele, 4 years ago

fix to CifFile?
fixes to cbf & tif importers for new Pilatus 2M detector
fixes to mag structure display

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