source: trunk/imports/G2phase_CIF.py @ 3334

Last change on this file since 3334 was 3334, checked in by vondreele, 5 years ago

correct issues with import & display of incommensurate magnetic spin operations

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