source: trunk/imports/G2phase_CIF.py @ 3210

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

some fixes for magnetic incommensurate cifs - dealing with 2 different dictionary defns for mcif!

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