source: trunk/imports/G2phase_CIF.py @ 3168

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

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