source: trunk/imports/G2phase_CIF.py @ 3219

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

fix import mag commensurate structures - all drawn correctly now
modulated mags not correct yet

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