source: trunk/imports/G2phase_CIF.py @ 3465

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

deal with export/import issues for mag cif files
fix bug for new phase - append atom (space group change problem)

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