source: trunk/imports/G2phase_CIF.py @ 3163

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

fix typo in Tb neutron scattering factors; extra space
fix lst output of neutron resonant form factors - computations were correct
implement import of mcif files.

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