source: trunk/imports/G2phase_CIF.py @ 3137

Last change on this file since 3137 was 3137, checked in by vondreele, 5 years ago

replace old CifFile? with new py 2/7/3.6 compliant code
fix cif file import phase & powder file
fix CemComp? restraint editing

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