source: trunk/imports/G2phase_CIF.py @ 3197

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

Rearrange SingleStringDialog? window - put '?' at end of text box
modify space group stuff & cif importer to read any magnetic cif file.

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