source: trunk/imports/G2phase_CIF.py @ 3491

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

fix export/import of magnetic cifs

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