source: trunk/imports/G2phase_CIF.py @ 3211

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

fix integer divide in wave display
some fixes to mcif import

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 42.3 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2018-01-01 20:12:54 +0000 (Mon, 01 Jan 2018) $
4# $Author: vondreele $
5# $Revision: 3211 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3211 2018-01-01 20:12:54Z 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: 3211 $")
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                if not sg: sg = cf[blknm].get("_space_group_ssg_name",'')
126                if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'')
127                if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'')
128                #how about checking for super/magnetic ones as well? - reject 'X'?
129                sg = sg.replace('_','')
130                if sg: choice[-1] += ', (' + sg.strip() + ')'
131            selblk = G2IO.PhaseSelector(choice,ParentFrame=ParentFrame,
132                title= 'Select a phase from one the CIF data_ blocks below',size=(600,100))
133        self.errors = 'Error during reading of selected block'
134#process selected phase
135        if selblk is None:
136            returnstat = False # no block selected or available
137        else:   #do space group symbol & phase type first
138            blknm = str_blklist[selblk]
139            blk = cf[str_blklist[selblk]]
140            E = True
141            Super = False
142            magnetic = False
143            moddim = int(blk.get("_cell_modulation_dimension",'0'))
144            if moddim:      #incommensurate
145                if moddim > 1:
146                    msg = 'more than 3+1 super symmetry is not allowed in GSAS-II'
147                    self.errors = msg
148                    self.warnings += '\n'+msg
149                    return False
150                if blk.get('_cell_subsystems_number'):
151                    msg = 'Composite super structures not allowed in GSAS-II'
152                    self.errors = msg
153                    self.warnings += '\n'+msg
154                    return False
155                sspgrp = blk.get("_space_group_ssg_name",'')
156                if not sspgrp:          #might be incommensurate magnetic
157                    MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'')
158                    if not MSSpGrp:
159                        MSSpGrp = blk.get("_space_group.magn_ssg_name",'') 
160                    if not MSSpGrp:
161                        msg = 'No incommensurate space group name was found in the CIF.'
162                        self.errors = msg
163                        self.warnings += '\n'+msg
164                        return False                                                           
165                    if 'X' in MSSpGrp:
166                        msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II'
167                        self.warnings += '\n'+msg
168                        self.errors = msg
169                        return False
170                    magnetic = True
171                if 'X' in sspgrp:
172                    msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II'
173                    self.warnings += '\n'+msg
174                    self.errors = msg
175                    return False
176                Super = True
177                if magnetic:
178                    sspgrp = MSSpGrp.split('(')
179                    sspgrp[1] = "("+sspgrp[1]
180                    SpGrp = G2spc.StandardizeSpcName(sspgrp[0])
181                    MSpGrp = sspgrp[0]
182                    self.MPhase['General']['Type'] = 'magnetic'
183                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
184                else:
185                    sspgrp = sspgrp.split('(')
186                    sspgrp[1] = "("+sspgrp[1]
187                    SpGrp = sspgrp[0]
188                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
189                    self.Phase['General']['Type'] = 'nuclear'
190                if not SpGrp:
191                    print (sspgrp)
192                    self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II'
193                    return False
194                SuperSg = sspgrp[1].replace('\\','')
195                SuperVec = [[0,0,.1],False,4]
196            else:   #not incommensurate
197                SpGrp = blk.get("_symmetry_space_group_name_H-M",'')
198                if not SpGrp:
199                    SpGrp = blk.get("_space_group_name_H-M_alt",'')
200                if not SpGrp:   #try magnetic           
201                    MSpGrp = blk.get("_space_group.magn_name_BNS",'')
202                    if not MSpGrp:
203                        MSpGrp = blk.get("_space_group_magn.name_BNS",'')
204                        if not MSpGrp:
205                            msg = 'No recognizable space group name was found in the CIF.'
206                            self.errors = msg
207                            self.warnings += '\n'+msg
208                            return False                   
209                    SpGrp = MSpGrp.replace("'",'')
210                    SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','')   #get rid of screw '_'
211                    if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:]
212                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
213                    magnetic = True
214                    self.MPhase['General']['Type'] = 'magnetic'
215                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
216                    if not SpGrp:
217                        print (MSpGrp)
218                        self.warnings += 'No space group name was found in the CIF.'
219                        return False
220                else:
221                    SpGrp = SpGrp.replace('_','')
222                    self.Phase['General']['Type'] = 'nuclear'
223#process space group symbol
224            E,SGData = G2spc.SpcGroup(SpGrp)
225            if E and SpGrp:
226                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
227                if SpGrpNorm:
228                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
229            # nope, try the space group "out of the Box"
230            if E:
231                self.warnings += 'ERROR in space group symbol '+SpGrp
232                self.warnings += '\nThe space group has been set to "P 1". '
233                self.warnings += "Change this in phase's General tab."
234                self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
235                self.warnings += G2spc.SGErrors(E)
236                SGData = G2obj.P1SGData # P 1
237            self.Phase['General']['SGData'] = SGData
238
239            if magnetic and not Super:
240                SGData['SGFixed'] = True
241                try:
242                    sgoploop = blk.GetLoop('_space_group_symop_magn.id')
243                    sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
244                    opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
245                    centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]                   
246                except KeyError:        #old mag cif names
247                    sgoploop = blk.GetLoop('_space_group_symop.magn_id')
248                    sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id')
249                    opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1]
250                    centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1]
251                SGData['SGOps'] = []
252                SGData['SGCen'] = []
253                spnflp = []
254                for op in sgoploop:
255                    M,T,S = G2spc.MagText2MTS(op[opid])
256                    SGData['SGOps'].append([np.array(M,dtype=float),T])
257                    spnflp.append(S)
258                censpn = []
259                for cent in sgcenloop:
260                    M,C,S = G2spc.MagText2MTS(cent[centid])
261                    SGData['SGCen'].append(C)
262                    censpn += list(np.array(spnflp)*S)
263                self.MPhase['General']['SGData'] = SGData
264                self.MPhase['General']['SGData']['SpnFlp'] = censpn
265                self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp
266                MagPtGp = blk.get('_space_group.magn_point_group')
267                if not MagPtGp:
268                    MagPtGp = blk.get('_space_group_magn.point_group_name')
269                self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp
270#                GenSym,GenFlg = G2spc.GetGenSym(SGData)
271#                self.MPhase['General']['SGData']['GenSym'] = GenSym
272#                self.MPhase['General']['SGData']['GenFlg'] = GenFlg
273
274            if Super:
275                E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
276                if E:
277                    self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
278                    self.warnings += '\n'+E
279                    SuperSg = SuperSg[:SuperSg.index(')')+1]
280                    self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
281                    E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
282                self.Phase['General']['SSGData'] = SSGData
283                if magnetic:
284                    self.MPhase['General']['SGData'] = SGData
285#                    self.MPhase['General']['SGData']['SpnFlp'] = censpn
286                    self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','')
287                    self.MPhase['General']['SGData']['MagPtGp'] = blk.get('_space_group.magn_point_group')
288                    self.MPhase['General']['SSGData'] = SSGData
289    #                GenSym,GenFlg = G2spc.GetGenSym(SGData)
290    #                self.MPhase['General']['SGData']['GenSym'] = GenSym
291    #                self.MPhase['General']['SGData']['GenFlg'] = GenFlg
292                   
293
294            # cell parameters
295            cell = []
296            for lbl in cellitems:
297                cell.append(cif.get_number_with_esd(blk[lbl])[0])
298            Volume = G2lat.calc_V(G2lat.cell2A(cell))
299            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
300            if magnetic:
301                self.MPhase['General']['Cell'] = [False,]+cell+[Volume,]               
302            if Super:
303                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
304                waveDict = dict(waveloop.items())
305                SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
306                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
307                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
308
309            # read in atoms
310            self.errors = 'Error during reading of atoms'
311            atomlbllist = [] # table to look up atom IDs
312            atomloop = blk.GetLoop('_atom_site_label')
313            atomkeys = [i.lower() for i in atomloop.keys()]
314            if not blk.get('_atom_site_type_symbol'):
315                self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
316            if magnetic:
317                try:
318                    magmoment = '_atom_site_moment.label'
319                    magatomloop = blk.GetLoop(magmoment)
320                    magatomkeys = [i.lower() for i in magatomloop.keys()]
321                    magatomlabels = blk.get(magmoment)
322                    G2MagDict = {'_atom_site_moment.label': 0,
323                                 '_atom_site_moment.crystalaxis_x':7,
324                                 '_atom_site_moment.crystalaxis_y':8,
325                                 '_atom_site_moment.crystalaxis_z':9}
326                except KeyError:
327                    magmoment = '_atom_site_moment_label'
328                    magatomloop = blk.GetLoop(magmoment)
329                    magatomkeys = [i.lower() for i in magatomloop.keys()]
330                    magatomlabels = blk.get(magmoment)
331                    G2MagDict = {'_atom_site_moment_label': 0,
332                                 '_atom_site_moment_crystalaxis_x':7,
333                                 '_atom_site_moment_crystalaxis_y':8,
334                                 '_atom_site_moment_crystalaxis_z':9}
335                   
336            if blk.get('_atom_site_aniso_label'):
337                anisoloop = blk.GetLoop('_atom_site_aniso_label')
338                anisokeys = [i.lower() for i in anisoloop.keys()]
339                anisolabels = blk.get('_atom_site_aniso_label')
340            else:
341                anisoloop = None
342                anisokeys = []
343                anisolabels = []
344            if Super:
345                occFloop = None
346                occCloop = None
347                occFdict = {}
348                occCdict = {}
349                displSloop = None
350                displFloop = None
351                MagFloop = None
352                displSdict = {}
353                displFdict = {}
354                MagFdict = {}
355                UijFloop = None
356                UijFdict = {}
357                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
358                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
359                    occFdict = dict(occFloop.items())
360                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
361                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
362                    occCdict = dict(occCloop.items())
363                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
364                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
365                    displFdict = dict(displFloop.items())                           
366                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
367                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
368                    displSdict = dict(displSloop.items())
369                if blk.get('_atom_site_U_Fourier_atom_site_label'):
370                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
371                    UijFdict = dict(UijFloop.items())
372                if blk.get('_atom_site_moment_Fourier_atom_site_label'):
373                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label')
374                    MagFdict = dict(MagFloop.items())
375                    Mnames =  ['_atom_site_moment_fourier_atom_site_label',
376                               '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id',
377                               '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos']
378                elif blk.get('_atom_site_moment_Fourier.atom_site_label'):
379                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label')
380                    MagFdict = dict(MagFloop.items())                           
381                    Mnames =  ['_atom_site_moment_fourier.atom_site_label',
382                               '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id',
383                               '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos']
384            self.Phase['Atoms'] = []
385            if magnetic:
386                self.MPhase['Atoms'] = []
387            G2AtomDict = {  '_atom_site_type_symbol' : 1,
388                            '_atom_site_label' : 0,
389                            '_atom_site_fract_x' : 3,
390                            '_atom_site_fract_y' : 4,
391                            '_atom_site_fract_z' : 5,
392                            '_atom_site_occupancy' : 6,
393                            '_atom_site_aniso_u_11' : 11,
394                            '_atom_site_aniso_u_22' : 12,
395                            '_atom_site_aniso_u_33' : 13,
396                            '_atom_site_aniso_u_12' : 14,
397                            '_atom_site_aniso_u_13' : 15,
398                            '_atom_site_aniso_u_23' : 16, }
399
400            ranIdlookup = {}
401            for aitem in atomloop:
402                atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
403                for val,key in zip(aitem,atomkeys):
404                    col = G2AtomDict.get(key,-1)
405                    if col >= 3:
406                        atomlist[col] = cif.get_number_with_esd(val)[0]
407                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
408                    elif col is not None and col != -1:
409                        atomlist[col] = val
410                    elif key in ('_atom_site_thermal_displace_type',
411                               '_atom_site_adp_type'):   #Iso or Aniso?
412                        if val.lower() == 'uani':
413                            atomlist[9] = 'A'
414                    elif key == '_atom_site_u_iso_or_equiv':
415                        uisoval = cif.get_number_with_esd(val)[0]
416                        if uisoval is not None: 
417                            atomlist[10] = uisoval
418                if not atomlist[1] and atomlist[0]:
419                    typ = atomlist[0].rstrip('0123456789-+')
420                    if G2elem.CheckElement(typ):
421                        atomlist[1] = typ
422                    if not atomlist[1]:
423                        atomlist[1] = 'Xe'
424                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
425                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
426                    atomlist[9] = 'A'
427                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
428                        col = G2AtomDict.get(key)
429                        if col:
430                            atomlist[col] = cif.get_number_with_esd(val)[0]
431                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
432                atomlist[1] = G2elem.FixValence(atomlist[1])
433                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
434                self.Phase['Atoms'].append(atomlist)
435                ranIdlookup[atomlist[0]] = atomlist[-1]
436                if atomlist[0] in atomlbllist:
437                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
438                else:
439                    atomlbllist.append(atomlist[0])
440
441                if magnetic and atomlist[0] in magatomlabels:
442                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
443                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
444                        mcol = G2MagDict.get(mkey,-1)
445                        if mcol:
446                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
447                    self.MPhase['Atoms'].append(matomlist)
448                if Super:
449                    Sfrac = []
450                    Sadp = np.zeros((4,12))
451                    Spos = np.zeros((4,6))
452                    Smag = np.zeros((4,6))
453                    nim = -1
454                    waveType = 'Fourier'                               
455                    if displFdict:
456                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
457                            if item == atomlist[0]:
458                                waveType = 'Fourier'                               
459                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
460                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
461                                if im != nim:
462                                    nim = im
463                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
464                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
465                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
466                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
467                    if nim >= 0:
468                        Spos = [[spos,False] for spos in Spos[:nim]]
469                    else:
470                        Spos = []
471                    if UijFdict:
472                        nim = -1
473                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
474                            if item == atomlist[0]:
475                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
476                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
477                                if im != nim:
478                                    nim = im
479                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
480                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
481                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
482                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
483                    if nim >= 0:
484                        Sadp = [[sadp,False] for sadp in Sadp[:nim]]
485                    else:
486                        Sadp = []
487                    if MagFdict:
488                        nim = -1
489                        for i,item in enumerate(MagFdict[Mnames[0]]):
490                            if item == atomlist[0]:
491                                waveType = 'Fourier'                               
492                                ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
493                                im = int(MagFdict[Mnames[2]][i])
494                                if im != nim:
495                                    nim = im
496                                val = MagFdict[Mnames[3]][i]
497                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
498                                val = MagFdict[Mnames[4]][i]
499                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
500                        if nim >= 0:
501                            Smag = [[smag,False] for smag in Smag[:nim]]
502                        else:
503                            Smag = []
504                    SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
505                    if magnetic and atomlist[0] in magatomlabels:
506                        matomlist.append(SSdict)
507                    atomlist.append(SSdict)
508            if len(atomlbllist) != len(self.Phase['Atoms']):
509                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
510            for lbl in phasenamefields: # get a name for the phase
511                name = blk.get(lbl)
512                if name is None:
513                    continue
514                name = name.strip()
515                if name == '?' or name == '.':
516                    continue
517                else:
518                    break
519            else: # no name found, use block name for lack of a better choice
520                name = blknm
521            self.Phase['General']['Name'] = name.strip()[:20]
522            self.Phase['General']['Super'] = Super
523            if magnetic:
524                self.MPhase['General']['Type'] = 'magnetic'               
525                self.MPhase['General']['Name'] = name.strip()[:20]+' mag'
526                self.MPhase['General']['Super'] = Super
527                if Super:
528                    self.MPhase['General']['Modulated'] = True
529                    self.MPhase['General']['SuperVec'] = SuperVec
530                    self.MPhase['General']['SuperSg'] = SuperSg
531                    self.MPhase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
532            else:
533                self.MPhase = None
534            if Super:
535                self.Phase['General']['Modulated'] = True
536                self.Phase['General']['SuperVec'] = SuperVec
537                self.Phase['General']['SuperSg'] = SuperSg
538                self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
539            if not self.isodistort_warnings:
540                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
541                    self.errors = "Error while processing ISODISTORT constraints"
542                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
543            else:
544                self.warnings += self.isodistort_warnings
545            returnstat = True
546        return returnstat
547
548    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
549        'Process ISODISTORT items to create constraints etc.'
550        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
551        'Maps ISODISTORT parm names to GSAS-II names'
552        #----------------------------------------------------------------------
553        # read in the ISODISTORT displacement modes
554        #----------------------------------------------------------------------
555        self.Constraints = []
556        explaination = {}
557        if blk.get('_iso_displacivemode_label'):
558            modelist = []
559            shortmodelist = []
560            for lbl in blk.get('_iso_displacivemode_label'):
561                modelist.append(lbl)
562                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
563                # where SSSSS is the parent spacegroup, [x,y,z] is a location
564                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
565                # this extracts the AAAAA and BBBBB parts of the string
566                if regexp:
567                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
568                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
569            # read in the coordinate offset variables names and map them to G2 names/objects
570            coordVarLbl = []
571            G2varLbl = []
572            G2varObj = []
573            error = False
574            for lbl in blk.get('_iso_deltacoordinate_label'):
575                coordVarLbl.append(lbl)
576                if '_' in lbl:
577                    albl = lbl[:lbl.rfind('_')]
578                    vlbl = lbl[lbl.rfind('_')+1:]
579                else:
580                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
581                    error = True
582                    continue
583                if albl not in atomlbllist:
584                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
585                    error = True
586                    continue
587                else:
588                    anum = atomlbllist.index(albl)
589                var = varLookup.get(vlbl)
590                if not var:
591                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
592                    error = True
593                    continue
594                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
595                G2varObj.append(G2obj.G2VarObj(
596                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
597                    ))
598            if error:
599                raise Exception("Error decoding variable labels")
600
601            if len(G2varObj) != len(modelist):
602                print ("non-square input")
603                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
604
605            error = False
606            ParentCoordinates = {}
607            for lbl,exp in zip(
608                blk.get('_iso_coordinate_label'),
609                blk.get('_iso_coordinate_formula'),
610                ):
611                if '_' in lbl:
612                    albl = lbl[:lbl.rfind('_')]
613                    vlbl = lbl[lbl.rfind('_')+1:]
614                else:
615                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
616                    error = True
617                    continue
618                if vlbl not in 'xyz' or len(vlbl) != 1:
619                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
620                    error = True
621                    continue
622                i = 'xyz'.index(vlbl)
623                if not ParentCoordinates.get(albl):
624                    ParentCoordinates[albl] = [None,None,None]
625                if '+' in exp:
626                    val = exp.split('+')[0].strip()
627                    val = G2p3.FormulaEval(val)
628                    if val is None:
629                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
630                        error = True
631                        continue
632                    ParentCoordinates[albl][i] = val
633                else:
634                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
635            if error:
636                print (self.warnings)
637                raise Exception("Error decoding variable labels")
638            # get mapping of modes to atomic coordinate displacements
639            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
640            for row,col,val in zip(
641                blk.get('_iso_displacivemodematrix_row'),
642                blk.get('_iso_displacivemodematrix_col'),
643                blk.get('_iso_displacivemodematrix_value'),):
644                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
645            # Invert to get mapping of atom displacements to modes
646            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
647            # create the constraints
648            for i,row in enumerate(displacivemodeInvmatrix):
649                constraint = []
650                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
651                    if k == 0: continue
652                    constraint.append([k,G2varObj[j]])
653                constraint += [shortmodelist[i],False,'f']
654                self.Constraints.append(constraint)
655            #----------------------------------------------------------------------
656            # save the ISODISTORT info for "mode analysis"
657            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
658            self.Phase['ISODISTORT'].update({
659                'IsoModeList' : modelist,
660                'G2ModeList' : shortmodelist,
661                'IsoVarList' : coordVarLbl,
662                'G2VarList' : G2varObj,
663                'ParentStructure' : ParentCoordinates,
664                'Var2ModeMatrix' : displacivemodeInvmatrix,
665                'Mode2VarMatrix' : displacivemodematrix,
666                })
667            # make explaination dictionary
668            for mode,shortmode in zip(modelist,shortmodelist):
669                explaination[shortmode] = "ISODISTORT full name "+str(mode)
670        #----------------------------------------------------------------------
671        # now read in the ISODISTORT occupancy modes
672        #----------------------------------------------------------------------
673        if blk.get('_iso_occupancymode_label'):
674            modelist = []
675            shortmodelist = []
676            for lbl in blk.get('_iso_occupancymode_label'):
677                modelist.append(lbl)
678                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
679                # where SSSSS is the parent spacegroup, [x,y,z] is a location
680                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
681                # this extracts the AAAAA and BBBBB parts of the string
682                if regexp:
683                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
684                lbl = lbl.replace('order','o')
685                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
686            # read in the coordinate offset variables names and map them to G2 names/objects
687            occVarLbl = []
688            G2varLbl = []
689            G2varObj = []
690            error = False
691            for lbl in blk.get('_iso_deltaoccupancy_label'):
692                occVarLbl.append(lbl)
693                if '_' in lbl:
694                    albl = lbl[:lbl.rfind('_')]
695                    vlbl = lbl[lbl.rfind('_')+1:]
696                else:
697                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
698                    error = True
699                    continue
700                if albl not in atomlbllist:
701                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
702                    error = True
703                    continue
704                else:
705                    anum = atomlbllist.index(albl)
706                var = varLookup.get(vlbl)
707                if not var:
708                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
709                    error = True
710                    continue
711                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
712                G2varObj.append(G2obj.G2VarObj(
713                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
714                    ))
715            if error:
716                raise Exception("Error decoding variable labels")
717
718            if len(G2varObj) != len(modelist):
719                print ("non-square input")
720                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
721
722            error = False
723            ParentCoordinates = {}
724            for lbl,exp in zip(
725                blk.get('_iso_occupancy_label'),
726                blk.get('_iso_occupancy_formula'),
727                ):
728                if '_' in lbl:
729                    albl = lbl[:lbl.rfind('_')]
730                    vlbl = lbl[lbl.rfind('_')+1:]
731                else:
732                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
733                    error = True
734                    continue
735                if vlbl != 'occ':
736                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
737                    error = True
738                    continue
739                if '+' in exp:
740                    val = exp.split('+')[0].strip()
741                    val = G2p3.FormulaEval(val)
742                    if val is None:
743                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
744                        error = True
745                        continue
746                    ParentCoordinates[albl] = val
747            if error:
748                raise Exception("Error decoding occupancy labels")
749            # get mapping of modes to atomic coordinate displacements
750            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
751            for row,col,val in zip(
752                blk.get('_iso_occupancymodematrix_row'),
753                blk.get('_iso_occupancymodematrix_col'),
754                blk.get('_iso_occupancymodematrix_value'),):
755                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
756            # Invert to get mapping of atom displacements to modes
757            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
758            # create the constraints
759            for i,row in enumerate(occupancymodeInvmatrix):
760                constraint = []
761                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
762                    if k == 0: continue
763                    constraint.append([k,G2varObj[j]])
764                constraint += [shortmodelist[i],False,'f']
765                self.Constraints.append(constraint)
766            #----------------------------------------------------------------------
767            # save the ISODISTORT info for "mode analysis"
768            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
769            self.Phase['ISODISTORT'].update({
770                'OccModeList' : modelist,
771                'G2OccModeList' : shortmodelist,
772                'OccVarList' : occVarLbl,
773                'G2OccVarList' : G2varObj,
774                'BaseOcc' : ParentCoordinates,
775                'Var2OccMatrix' : occupancymodeInvmatrix,
776                'Occ2VarMatrix' : occupancymodematrix,
777                })
778            # make explaination dictionary
779            for mode,shortmode in zip(modelist,shortmodelist):
780                explaination[shortmode] = "ISODISTORT full name "+str(mode)
781        #----------------------------------------------------------------------
782        # done with read
783        #----------------------------------------------------------------------
784        if explaination: self.Constraints.append(explaination)
785
786        # # debug: show the mode var to mode relations
787        # for i,row in enumerate(displacivemodeInvmatrix):
788        #     l = ''
789        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
790        #         if k == 0: continue
791        #         if l: l += ' + '
792        #         #l += lbl+' * '+str(k)
793        #         l += G2varLbl[j]+' * '+str(k)
794        #     print str(i) + ': '+shortmodelist[i]+' = '+l
795        # print 70*'='
796
797        # # debug: Get the ISODISTORT offset values
798        # coordVarDelta = {}
799        # for lbl,val in zip(
800        #     blk.get('_iso_deltacoordinate_label'),
801        #     blk.get('_iso_deltacoordinate_value'),):
802        #     coordVarDelta[lbl] = float(val)
803        # modeVarDelta = {}
804        # for lbl,val in zip(
805        #     blk.get('_iso_displacivemode_label'),
806        #     blk.get('_iso_displacivemode_value'),):
807        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
808
809        # print 70*'='
810        # # compute the mode values from the reported coordinate deltas
811        # for i,row in enumerate(displacivemodeInvmatrix):
812        #     l = ''
813        #     sl = ''
814        #     s = 0.
815        #     for lbl,k in zip(coordVarLbl,row):
816        #         if k == 0: continue
817        #         if l: l += ' + '
818        #         l += lbl+' * '+str(k)
819        #         if sl: sl += ' + '
820        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
821        #         s += coordVarDelta[lbl] * k
822        #     print 'a'+str(i)+' = '+l
823        #     print '\t= '+sl
824        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
825        #     print
826
827        # print 70*'='
828        # # compute the coordinate displacements from the reported mode values
829        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
830        #     l = ''
831        #     sl = ''
832        #     s = 0.0
833        #     for j,k in enumerate(row):
834        #         if k == 0: continue
835        #         if l: l += ' + '
836        #         l += 'a'+str(j+1)+' * '+str(k)
837        #         if sl: sl += ' + '
838        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
839        #         s += modeVarDelta[modelist[j]] * k
840        #     print lbl+' = '+l
841        #     print '\t= '+sl
842        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
843        #     print
844
845        # determine the coordinate delta values from deviations from the parent structure
846        # for atmline in self.Phase['Atoms']:
847        #     lbl = atmline[0]
848        #     x,y,z = atmline[3:6]
849        #     if lbl not in ParentCoordinates:
850        #         print lbl,x,y,z
851        #         continue
852        #     px,py,pz = ParentCoordinates[lbl]
853        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.