source: trunk/imports/G2phase_CIF.py @ 3208

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

fix transformation bug abc -> abc* needed transform of matrix - G2phsGUI line 373
work on modullated magnetic cif import; fix some magnetic & incommensurate bugs - more to do.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 41.5 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-12-27 16:41:30 +0000 (Wed, 27 Dec 2017) $
4# $Author: vondreele $
5# $Revision: 3208 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3208 2017-12-27 16:41:30Z 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: 3208 $")
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            self.Phase['Atoms'] = []
375            if magnetic:
376                self.MPhase['Atoms'] = []
377            G2AtomDict = {  '_atom_site_type_symbol' : 1,
378                            '_atom_site_label' : 0,
379                            '_atom_site_fract_x' : 3,
380                            '_atom_site_fract_y' : 4,
381                            '_atom_site_fract_z' : 5,
382                            '_atom_site_occupancy' : 6,
383                            '_atom_site_aniso_u_11' : 11,
384                            '_atom_site_aniso_u_22' : 12,
385                            '_atom_site_aniso_u_33' : 13,
386                            '_atom_site_aniso_u_12' : 14,
387                            '_atom_site_aniso_u_13' : 15,
388                            '_atom_site_aniso_u_23' : 16, }
389
390            ranIdlookup = {}
391            for aitem in atomloop:
392                atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
393                for val,key in zip(aitem,atomkeys):
394                    col = G2AtomDict.get(key,-1)
395                    if col >= 3:
396                        atomlist[col] = cif.get_number_with_esd(val)[0]
397                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
398                    elif col is not None and col != -1:
399                        atomlist[col] = val
400                    elif key in ('_atom_site_thermal_displace_type',
401                               '_atom_site_adp_type'):   #Iso or Aniso?
402                        if val.lower() == 'uani':
403                            atomlist[9] = 'A'
404                    elif key == '_atom_site_u_iso_or_equiv':
405                        uisoval = cif.get_number_with_esd(val)[0]
406                        if uisoval is not None: 
407                            atomlist[10] = uisoval
408                if not atomlist[1] and atomlist[0]:
409                    typ = atomlist[0].rstrip('0123456789-+')
410                    if G2elem.CheckElement(typ):
411                        atomlist[1] = typ
412                    if not atomlist[1]:
413                        atomlist[1] = 'Xe'
414                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
415                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
416                    atomlist[9] = 'A'
417                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
418                        col = G2AtomDict.get(key)
419                        if col:
420                            atomlist[col] = cif.get_number_with_esd(val)[0]
421                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
422                atomlist[1] = G2elem.FixValence(atomlist[1])
423                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
424                self.Phase['Atoms'].append(atomlist)
425                ranIdlookup[atomlist[0]] = atomlist[-1]
426                if atomlist[0] in atomlbllist:
427                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
428                else:
429                    atomlbllist.append(atomlist[0])
430
431                if magnetic and atomlist[0] in magatomlabels:
432                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
433                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
434                        mcol = G2MagDict.get(mkey,-1)
435                        if mcol:
436                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
437                    self.MPhase['Atoms'].append(matomlist)
438                if Super:
439                    Sfrac = []
440                    Sadp = np.zeros((4,12))
441                    Spos = np.zeros((4,6))
442                    Smag = np.zeros((4,6))
443                    nim = -1
444                    waveType = 'Fourier'                               
445                    if displFdict:
446                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
447                            if item == atomlist[0]:
448                                waveType = 'Fourier'                               
449                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
450                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
451                                if im != nim:
452                                    nim = im
453                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
454                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
455                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
456                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
457                    if nim >= 0:
458                        Spos = [[spos,False] for spos in Spos[:nim]]
459                    else:
460                        Spos = []
461                    if UijFdict:
462                        nim = -1
463                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
464                            if item == atomlist[0]:
465                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
466                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
467                                if im != nim:
468                                    nim = im
469                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
470                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
471                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
472                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
473                    if nim >= 0:
474                        Sadp = [[sadp,False] for sadp in Sadp[:nim]]
475                    else:
476                        Sadp = []
477                    if MagFdict:
478                        nim = -1
479                        for i,item in enumerate(MagFdict['_atom_site_moment_fourier_atom_site_label']):
480                            if item == atomlist[0]:
481                                waveType = 'Fourier'                               
482                                ix = ['x','y','z'].index(MagFdict['_atom_site_moment_fourier_axis'][i])
483                                im = int(MagFdict['_atom_site_moment_fourier_wave_vector_seq_id'][i])
484                                if im != nim:
485                                    nim = im
486                                val = MagFdict['_atom_site_moment_fourier_param_sin'][i]
487                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
488                                val = MagFdict['_atom_site_moment_fourier_param_cos'][i]
489                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
490                        if nim >= 0:
491                            Smag = [[smag,False] for smag in Smag[:nim]]
492                        else:
493                            Smag = []
494                    SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
495                    if magnetic and atomlist[0] in magatomlabels:
496                        matomlist.append(SSdict)
497                    atomlist.append(SSdict)
498            if len(atomlbllist) != len(self.Phase['Atoms']):
499                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
500            for lbl in phasenamefields: # get a name for the phase
501                name = blk.get(lbl)
502                if name is None:
503                    continue
504                name = name.strip()
505                if name == '?' or name == '.':
506                    continue
507                else:
508                    break
509            else: # no name found, use block name for lack of a better choice
510                name = blknm
511            self.Phase['General']['Name'] = name.strip()[:20]
512            self.Phase['General']['Super'] = Super
513            if magnetic:
514                self.MPhase['General']['Type'] = 'magnetic'               
515                self.MPhase['General']['Name'] = name.strip()[:20]+' mag'
516                self.MPhase['General']['Super'] = Super
517                if Super:
518                    self.MPhase['General']['Modulated'] = True
519                    self.MPhase['General']['SuperVec'] = SuperVec
520                    self.MPhase['General']['SuperSg'] = SuperSg
521                    self.MPhase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
522            else:
523                self.MPhase = None
524            if Super:
525                self.Phase['General']['Modulated'] = True
526                self.Phase['General']['SuperVec'] = SuperVec
527                self.Phase['General']['SuperSg'] = SuperSg
528                self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
529            if not self.isodistort_warnings:
530                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
531                    self.errors = "Error while processing ISODISTORT constraints"
532                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
533            else:
534                self.warnings += self.isodistort_warnings
535            returnstat = True
536        return returnstat
537
538    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
539        'Process ISODISTORT items to create constraints etc.'
540        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
541        'Maps ISODISTORT parm names to GSAS-II names'
542        #----------------------------------------------------------------------
543        # read in the ISODISTORT displacement modes
544        #----------------------------------------------------------------------
545        self.Constraints = []
546        explaination = {}
547        if blk.get('_iso_displacivemode_label'):
548            modelist = []
549            shortmodelist = []
550            for lbl in blk.get('_iso_displacivemode_label'):
551                modelist.append(lbl)
552                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
553                # where SSSSS is the parent spacegroup, [x,y,z] is a location
554                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
555                # this extracts the AAAAA and BBBBB parts of the string
556                if regexp:
557                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
558                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
559            # read in the coordinate offset variables names and map them to G2 names/objects
560            coordVarLbl = []
561            G2varLbl = []
562            G2varObj = []
563            error = False
564            for lbl in blk.get('_iso_deltacoordinate_label'):
565                coordVarLbl.append(lbl)
566                if '_' in lbl:
567                    albl = lbl[:lbl.rfind('_')]
568                    vlbl = lbl[lbl.rfind('_')+1:]
569                else:
570                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
571                    error = True
572                    continue
573                if albl not in atomlbllist:
574                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
575                    error = True
576                    continue
577                else:
578                    anum = atomlbllist.index(albl)
579                var = varLookup.get(vlbl)
580                if not var:
581                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
582                    error = True
583                    continue
584                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
585                G2varObj.append(G2obj.G2VarObj(
586                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
587                    ))
588            if error:
589                raise Exception("Error decoding variable labels")
590
591            if len(G2varObj) != len(modelist):
592                print ("non-square input")
593                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
594
595            error = False
596            ParentCoordinates = {}
597            for lbl,exp in zip(
598                blk.get('_iso_coordinate_label'),
599                blk.get('_iso_coordinate_formula'),
600                ):
601                if '_' in lbl:
602                    albl = lbl[:lbl.rfind('_')]
603                    vlbl = lbl[lbl.rfind('_')+1:]
604                else:
605                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
606                    error = True
607                    continue
608                if vlbl not in 'xyz' or len(vlbl) != 1:
609                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
610                    error = True
611                    continue
612                i = 'xyz'.index(vlbl)
613                if not ParentCoordinates.get(albl):
614                    ParentCoordinates[albl] = [None,None,None]
615                if '+' in exp:
616                    val = exp.split('+')[0].strip()
617                    val = G2p3.FormulaEval(val)
618                    if val is None:
619                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
620                        error = True
621                        continue
622                    ParentCoordinates[albl][i] = val
623                else:
624                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
625            if error:
626                print (self.warnings)
627                raise Exception("Error decoding variable labels")
628            # get mapping of modes to atomic coordinate displacements
629            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
630            for row,col,val in zip(
631                blk.get('_iso_displacivemodematrix_row'),
632                blk.get('_iso_displacivemodematrix_col'),
633                blk.get('_iso_displacivemodematrix_value'),):
634                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
635            # Invert to get mapping of atom displacements to modes
636            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
637            # create the constraints
638            for i,row in enumerate(displacivemodeInvmatrix):
639                constraint = []
640                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
641                    if k == 0: continue
642                    constraint.append([k,G2varObj[j]])
643                constraint += [shortmodelist[i],False,'f']
644                self.Constraints.append(constraint)
645            #----------------------------------------------------------------------
646            # save the ISODISTORT info for "mode analysis"
647            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
648            self.Phase['ISODISTORT'].update({
649                'IsoModeList' : modelist,
650                'G2ModeList' : shortmodelist,
651                'IsoVarList' : coordVarLbl,
652                'G2VarList' : G2varObj,
653                'ParentStructure' : ParentCoordinates,
654                'Var2ModeMatrix' : displacivemodeInvmatrix,
655                'Mode2VarMatrix' : displacivemodematrix,
656                })
657            # make explaination dictionary
658            for mode,shortmode in zip(modelist,shortmodelist):
659                explaination[shortmode] = "ISODISTORT full name "+str(mode)
660        #----------------------------------------------------------------------
661        # now read in the ISODISTORT occupancy modes
662        #----------------------------------------------------------------------
663        if blk.get('_iso_occupancymode_label'):
664            modelist = []
665            shortmodelist = []
666            for lbl in blk.get('_iso_occupancymode_label'):
667                modelist.append(lbl)
668                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
669                # where SSSSS is the parent spacegroup, [x,y,z] is a location
670                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
671                # this extracts the AAAAA and BBBBB parts of the string
672                if regexp:
673                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
674                lbl = lbl.replace('order','o')
675                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
676            # read in the coordinate offset variables names and map them to G2 names/objects
677            occVarLbl = []
678            G2varLbl = []
679            G2varObj = []
680            error = False
681            for lbl in blk.get('_iso_deltaoccupancy_label'):
682                occVarLbl.append(lbl)
683                if '_' in lbl:
684                    albl = lbl[:lbl.rfind('_')]
685                    vlbl = lbl[lbl.rfind('_')+1:]
686                else:
687                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
688                    error = True
689                    continue
690                if albl not in atomlbllist:
691                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
692                    error = True
693                    continue
694                else:
695                    anum = atomlbllist.index(albl)
696                var = varLookup.get(vlbl)
697                if not var:
698                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
699                    error = True
700                    continue
701                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
702                G2varObj.append(G2obj.G2VarObj(
703                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
704                    ))
705            if error:
706                raise Exception("Error decoding variable labels")
707
708            if len(G2varObj) != len(modelist):
709                print ("non-square input")
710                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
711
712            error = False
713            ParentCoordinates = {}
714            for lbl,exp in zip(
715                blk.get('_iso_occupancy_label'),
716                blk.get('_iso_occupancy_formula'),
717                ):
718                if '_' in lbl:
719                    albl = lbl[:lbl.rfind('_')]
720                    vlbl = lbl[lbl.rfind('_')+1:]
721                else:
722                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
723                    error = True
724                    continue
725                if vlbl != 'occ':
726                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
727                    error = True
728                    continue
729                if '+' in exp:
730                    val = exp.split('+')[0].strip()
731                    val = G2p3.FormulaEval(val)
732                    if val is None:
733                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
734                        error = True
735                        continue
736                    ParentCoordinates[albl] = val
737            if error:
738                raise Exception("Error decoding occupancy labels")
739            # get mapping of modes to atomic coordinate displacements
740            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
741            for row,col,val in zip(
742                blk.get('_iso_occupancymodematrix_row'),
743                blk.get('_iso_occupancymodematrix_col'),
744                blk.get('_iso_occupancymodematrix_value'),):
745                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
746            # Invert to get mapping of atom displacements to modes
747            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
748            # create the constraints
749            for i,row in enumerate(occupancymodeInvmatrix):
750                constraint = []
751                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
752                    if k == 0: continue
753                    constraint.append([k,G2varObj[j]])
754                constraint += [shortmodelist[i],False,'f']
755                self.Constraints.append(constraint)
756            #----------------------------------------------------------------------
757            # save the ISODISTORT info for "mode analysis"
758            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
759            self.Phase['ISODISTORT'].update({
760                'OccModeList' : modelist,
761                'G2OccModeList' : shortmodelist,
762                'OccVarList' : occVarLbl,
763                'G2OccVarList' : G2varObj,
764                'BaseOcc' : ParentCoordinates,
765                'Var2OccMatrix' : occupancymodeInvmatrix,
766                'Occ2VarMatrix' : occupancymodematrix,
767                })
768            # make explaination dictionary
769            for mode,shortmode in zip(modelist,shortmodelist):
770                explaination[shortmode] = "ISODISTORT full name "+str(mode)
771        #----------------------------------------------------------------------
772        # done with read
773        #----------------------------------------------------------------------
774        if explaination: self.Constraints.append(explaination)
775
776        # # debug: show the mode var to mode relations
777        # for i,row in enumerate(displacivemodeInvmatrix):
778        #     l = ''
779        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
780        #         if k == 0: continue
781        #         if l: l += ' + '
782        #         #l += lbl+' * '+str(k)
783        #         l += G2varLbl[j]+' * '+str(k)
784        #     print str(i) + ': '+shortmodelist[i]+' = '+l
785        # print 70*'='
786
787        # # debug: Get the ISODISTORT offset values
788        # coordVarDelta = {}
789        # for lbl,val in zip(
790        #     blk.get('_iso_deltacoordinate_label'),
791        #     blk.get('_iso_deltacoordinate_value'),):
792        #     coordVarDelta[lbl] = float(val)
793        # modeVarDelta = {}
794        # for lbl,val in zip(
795        #     blk.get('_iso_displacivemode_label'),
796        #     blk.get('_iso_displacivemode_value'),):
797        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
798
799        # print 70*'='
800        # # compute the mode values from the reported coordinate deltas
801        # for i,row in enumerate(displacivemodeInvmatrix):
802        #     l = ''
803        #     sl = ''
804        #     s = 0.
805        #     for lbl,k in zip(coordVarLbl,row):
806        #         if k == 0: continue
807        #         if l: l += ' + '
808        #         l += lbl+' * '+str(k)
809        #         if sl: sl += ' + '
810        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
811        #         s += coordVarDelta[lbl] * k
812        #     print 'a'+str(i)+' = '+l
813        #     print '\t= '+sl
814        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
815        #     print
816
817        # print 70*'='
818        # # compute the coordinate displacements from the reported mode values
819        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
820        #     l = ''
821        #     sl = ''
822        #     s = 0.0
823        #     for j,k in enumerate(row):
824        #         if k == 0: continue
825        #         if l: l += ' + '
826        #         l += 'a'+str(j+1)+' * '+str(k)
827        #         if sl: sl += ' + '
828        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
829        #         s += modeVarDelta[modelist[j]] * k
830        #     print lbl+' = '+l
831        #     print '\t= '+sl
832        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
833        #     print
834
835        # determine the coordinate delta values from deviations from the parent structure
836        # for atmline in self.Phase['Atoms']:
837        #     lbl = atmline[0]
838        #     x,y,z = atmline[3:6]
839        #     if lbl not in ParentCoordinates:
840        #         print lbl,x,y,z
841        #         continue
842        #     px,py,pz = ParentCoordinates[lbl]
843        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.