source: trunk/imports/G2phase_CIF.py @ 3192

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

the cif reader will now read magnetic cif files as two phases - one nuclear with all atoms & the other magnetic with just the magnetic atoms. NB: spin flips are not yet properly assigned; working on that.

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