source: trunk/imports/G2phase_CIF.py @ 3191

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

fix integer divides in responding to unit cell changes on Unit Cells window
show StarFile? errors from cif importer
add missing nonstandard space group symbols for face centered orthorhombics
begin implementation of cif import magnetic phases into 2 phase result

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