source: trunk/imports/G2phase_CIF.py @ 3710

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

fix issues about site sym & multilicity for mcif import of gray (incommensurate magnetic) groups

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